UM ALGORÍTMO DE ORDEM SUPERIOR PARA A INTEGRAÇÃO DIRETA DAS

EQUAÇÕES DA DINÂMICA ESTRUTURAL

Helio José Corrêa Barbosa

TESE SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS PROGRAMAS DE

PÔS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE FEDERAL DO RIO DE JA­

NEIRO COMO PARTE DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO

GRAU DE MESTRE EM CIÊNCIAS (M.Sc.).

Aprovada por:

e

~L----~a nancio Filho

Agus ~n Juan Ferrante

4 Lu~ ::vfa;~;=

RIO DE JANEIRO, RJ - BRASIL

MARÇO DE 1978

i

A meus pais e

meus irmaos

ii

AGRADECIMENTOS

Ao Professor Fernando Venâncio Filho pela suges­

tao do tema e orientação dada ao presente trabalho.

A Abimael F.D. Loula e Nelson F.F. Ebecken pela

amizade, estimulo e colaboração dispensados.

Aos Professores da COPPE/UFRJ, na pessoa de seu

Diretor Paulo Alcântara Gomes, pelos ensinamentos ministrados.

Aos funcionãrios da Biblioteca Central do Centro

de Tecnologia e do Núcleo de Computação Eletrônica pela atençao.

A CAPES e COPPE pelo apoio financeiro.

A Helena Santos de Oliveira pelo esmerado traba­

lho de datilografia.

A Maria Mareia Correa Barbosa, pela dedicação na

confecção de toda a parte gráfica, meu especial agradecimento.

V

ÍNDICE

I - INTRODUÇÃO ........................................... 1

II - FORMULAÇÃO DAS EQUAÇÕES DE MOVIMENTO . .. .......... .. .. 2

2.1 - Introdução ..................................... 2

2.2 - Formulação das Equações de Movimento ........... 2

2.2.1 - Emprego das Leis de Newton . .... ........ 3

2. 2. 2 - Princípio dos Trabalhos Virtuais . . . . . . . 3

2.2.3 - Princípio de Hamilton ... .... .. .. .... .. . 4

2. 3 - Discretização Espacial . . . • . . . . . . . . . . . . . . . . . . . . . 6

2.3.1 - Mitodos Residuais . ........ ... ... ... ... . 7

2.3.2 - Mitodos Variacionais - Mitodo de Ritz .. 10

2.3.3 - Funções Localizadas - Mitodo dos Elemen-tos Finitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

III - M~TODOS DE SOLUÇÃO DAS EQUAÇÕES DE MOVIMENTO DO MODELO

DISCRETO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1 - Introdução . . . . . . . .. .. . . . . . . . . . . . . . .. . . . . . . . .. . . 14

3.2 - Características do Modelo Discreto .. ..... ... ... 14

3.2.1 - Introdução do Amortecimento .......... .. 18

3.3 - Superposição Modal 19

3. 4 - Integração Direta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.5 - AlgorÍtmo.s Baseados na Interpolação das Forças de Inircia . . . . . . . . . . . . . . . . . . • . . . . . . . . . . . . . . . . . . . . . 31

IV - PROCEDIMENTOS IMPLEMENTADOS - EXEMPLOS - CONCLUSÕES 56

4 .1 - Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.2 - Elementos Implantados .................. ........ 56

4. 2. 1 - Matrizes de Massa . . . . . . . . . . . . . . . . . . . . . . 58

4.2.2 - Matrizes de Rigidez ................. ... 65

4.3 - Elementos Escalares 70

4.4 - Matriz de Amortecimento ........................ 71

vi

4.5 - Vetor de Excitações Externas ..••........••.••••

4.6 - Critirio de Convergincia •...........•....•••..•

4.7 - Fluxograma ...•.•••.•....•............•........•

4. 8 - Exemplos de Aplicação •..•..•.•••.•...•••...••.•

4.9 - Conclusões •...•••.•••..•......••...•..•••••....

BIBLIOGRAFIA •••••.•••..•....•.............•...•.•••..••••••

APÊ:NDICE •••.•••..••••..•.••.••..•...•••..•..••....•..•.....

74

76

78

82

109

111

115

1

I - INTRODUÇÃO

A obtenção da resposta de sistemas estruturais sub

metidos a solicitações de natureza dinâmica é um problema que tem

merecido continua atenção por parte dos pesquisadores. Isto se

deve, por um lado, ã prÕpria necessidade de se resolver problemas

físicos cada vez mais complexos e por outro as dificuldades rela­

cionadas com a integração das equações diferenciais de movimento.

A procura de algoritmos que combinem precisao com

rapidez (eficiência computacional) continua para que problemas n~

vos e maiores possam ser abordados e resolvidos a um custo razoa-

vel.

No presente trabalho estuda-se um algoritmo de in­

tegraçao passo a passo desenvolvido em Q., 2, TI para a obtenção

da resposta, a solicitações de natureza determinística, de estru-

turas em que se considera não-linearidade geométrica

elásticos lineares.

e mateiiai.s

No Capitulo II a formulação das equaçoes de movi­

mento e alguns processos de discretização espacial são comentados

chegando-se às equações de movimento na forma discreta.

No Capitulo III métodos de solução dessas equaçoes

de movimento do modelo discreto são revistos.

Apresenta-se e discute-se o algoritmo implementa-

do, deixando-se para o Capitulo IV a descrição mais detalhada dos

procedimentos implantados bem como alguns exemplos de aplicação

do algoritmo e as conclusões obtidas.

2

II- FORMULAÇÃO DAS EQUAÇÕES DE MOVIMENTO

2.1 - INTRODUÇÃO

Os sistemas mecânicos, para fins de análise, podem

ser classificados em dois grupos, segundo a natureza das

veis cinemáticas (deslocamentos, velocidades e acelerações)

variaçao ao longo do tempo deseja-se conhecer.

Um sistema é dito d{~e~eto quando seu estado

variá-

cuja

fica

perfeitamente definido a partir do conhecimento dos valores assu­

midos por aquelas variáveis cinemáticas em um conjunto finito de

pontos desse sistema.

tempo.

O domínio das variáveis independentes é o

Diz-se eontZnuo quando suas propriedades mecanicas

sao funções contínuas de variáveis espaciais e a caracterização

do seu estado é feita mediante a adoção de variáveis

que variam continuamente com a posiçao e o tempo.

cineffiâticas

As equaçoes diferenciais de movimento sao ordiná-

rias no primeiro caso e parciais no caso de sistemas contínuos.

Normalmente os sistemas contínuos sao estudados por

meio de técnicas numéricas que conduzem a um sistema de

diferenciais ordinárias no tempo: modelo discreto.

Neste capítulo discute-se a formulação das

equaçoes

equa-

çoes de movimento de sistemas contínuos e algumas técnicas de dis

cretização.

2.2 - FORMULAÇÃO DAS EQUAÇÕES DE MOVIMENTO

As equaçoes de movimento de sistemas mecanicos po­

dem ser obtidas através de:

3

2.2.1 - Emprego da• Leis de N~~ton

De acordo com as leis do movimento de Newton, em

um referencial inercial, a variaçao da quantidade de movimento li

near de uma partícula e igual a resultante das forças aplicadas

sobre ela; e a variaçao da quantidade de movimento angular, com

relação à origem dos eixos coordenados e igual a resultante do mo

menta de todas as forças em relação a essa mesma origem.

A aplicação geral desses princípios a meios contí­

nuos é feita de maneira a assegurar que o movimento de cada parti

cula desse meio obedeça às leis de Newton e assim sao

equações diferenciais parciais de movimento

obtidas as

Entretanto e no caso de sistemas discretos ou sis­

temas contínuos simples que a aplicação direta das leis de Newton

e usualmente feita.

Introduzindo-se o conceito de forças de inércia,

pelo princípio de d 1 Alembert, as equaçoes de movimento passam a

ser obtidas de maneira análoga ao problema estático: a partir de

condições de equilíbrio.

As forças atuantes no sistema podem incluir forças

de restituição de diversos tipos e forças externas dependentes do

tempo e/ou de deslocamentos que, juntamente com as forças de iné~

eia e de ligação, deverão ser consideradas vetorialmente para o

estabelecimento do equilíbrio.

Em sistemas simples esta e, frequentemente, a for­

ma mais direta e conveniente de se obter as equações de movimento.

2.2.2 - Princípio dos Trabalhos Virtuais

O princípio dos trabalhos virtuais pode ser expre~

so como se segue. Se um sôlido, em equilíbrio sob a ação de for-

4

ças de volume e de superfície, e submetido a um deslocamento vir-

tual, isto e, um campo de deslocamentos contínuo infinitesimal,

arbitrário, porem compatível com as restrições impostas ao sÕlid~

o trabalho virtual realizado pelas forças externas e igual ao tra

balho virtual das solicitações internas. Ou seja:

( 2. 1)

A equaçao acima e pois uma afirmação de equilíbrio

estático. Sua extensao ao problema dinimico e feita através do

princípio de d'Alembert com a inclusão das forças de inercia na

expressao do trabalho virtual das forças externas.

O princípio tem a vantagem de nao fazer uso das

forças de ligação, jã que não realizam trabalho em deslocamentos

virtuais, o que não pode ser feito quando da aplicação direta das

leis de Newton. Observa-se que se trata de um princípio escalar.

2.2.3 - Princípio de Hamilton

Dado um sistema dinâmico sujeito a forças externas

e restrições geométricas, chama-se de configuração admissível do

sistema qualquer configuração que atenda as restrições impostas.

Ao conjunto de todas as configurações reais assumi

das continuamente pelo sistema entre dois instantes cha

ma-se caminho dinâmico4 Uma variaçao virtual em relação ao cami-

nho dinâmico gera um caminho variado.

Considerando-se apenas os caminhos que coincidem

com o caminho real em o princípio de Hamilton afirma en

tao que o caminho real tomado pelo sistema é tal que satisfaz a :

onde:

5

]' ~,, - ,, ' aw,J dt = o ( 2. 2)

tl

T = energia cinética total do sistema

V energia potencial total do sistema, incluindo energia de

deformação e o potencial de todas as forças externas con

servativas

trabalho virtual das forças não-conservativas

no sistema

atuantes

Para sistemas conservativos trata-se de tornar estacionário o fun

ciona 1:

(T - V) dt ( 2. 3)

que, sob certas condições [ 4 J ' ~ . assume um m1n1mo.

Esta formulação prescinde do conhecimento explíci­

to das forças atuantes e usa grandezas escalares, sendo a de em­

prego mais generalizado em problemas complexos.

Deve ser lembrado que os três procedimentos de for

mulação apresentados levam às mesmas equaçoes diferenciais de mo­

vimento ou equações de equilíbrio dinâmico.

Fazendo-se uso das equaçoes de compatibilidade pa­

ra ~istemas continuas e das equações constitutivas para materiais

elásticos obtém-se, num referencial inercial equações de movi

6

menta na forma:

+ M(r) u(r, t) =F(r,t) ( 2. 4)

onde:

D e um operador diferencial constituído de derivadas em re-

lação às coordenadas espaciais r e representando as ca-

M

racteristicas de rigidez do sistema.

e uma matriz função das coordenadas espaciais

rizando a distribuição de massa do sistema.

F e o vetor de ações aplicadas ao sistema

u e o vetor de deslocamentos incógnitos

r caracte

Para sistemas naturalmente discretos as equaçoes de movimento se-

riam:

Q ( t)

onde os pontos indicam derivação em relação ao tempo e:

F e o vetor de forças elásticas de restituição

M e a matriz de massa diagonal do sistema

g e o vetor de açoes externas aplicadas

g e o vetor de deslocamentos generalizados incógnitos

2.3 - DISCRETIZAÇÃO ESPACIAL

( 2 • 5)

Embora para alguns sistemas contínuos simples com

condições de contorno também simples seja possível encontrar so-

luçÕes exatas, tais casos constituem exceçao. Nos problemas reais

7

quando é possível formular as equaçoes diferenciais parciais de

movimento nem sempre e possível resolvê-las; e, quando se canse-

gue resolvê-las, podem nao ser satisfeitas as condições de contor

no.

são usadas entao têcnicas numéricas que conduZem :a

am sistema de equações diferenciais ordinárias no tempo

as equações de um sistema discreto.

análogas

Entre os diversos processos de discretização[}, ~

será feita referência apenas àqueles que empregam funções de apr~

ximaçao.

Estes processos assumem que a resposta do sistema

contínuo pode ser aproximada por uma série finita na forma:

onde <j,. ( r) -1 -

u (r , t) -n -

= n E

i=l <j,. ( r) -1 -

q . ( t) 1

sao as funções de aproximaçao dependentes

( 2. 6)

somente

da variável espacial r e q. ( t) 1

são coordenadas generalizadas

incógnitas funções do tempo.

O objetivo é escolher funções que sejam conve-

nientes do ponto de vista computacional e proporcionem uma boa

aproximaçao para a solução u

De acordo com o critério de aproximação adotado,

tais processos podem ser chamados de residuais ou variacionais.

2.3.1 - Métodos Residuais

A característica principal dos métodos residuais é

fazer uso direto da equação diferencial de movimento (2.4) intro

<luzindo nela a solução aproximada (2.6) que, em princípio, deverá

satisfazer a todas as condições de contorno.

8

Como a solução

origem a um resíduo:

-nao e a exata essa substituição dá

s(r, t) ~~J + M(r) t) ( 2. 7)

A idéia básica comum aos métodos residuais ê deter

minar a aproximação u -n

de maneira a minimizar, de alguma forma,

o erro 1:;<E, t) cometido.

Assim o diferente tratamento dado ao resíduo

s(r, t) dá origem a diversos processos tais como Galerkin e colo

cação, vistos a seguir.

a) - Método de Galerkin -

Nesse método o resíduo s (r, t) ê ponderado pelas

n funções de aproximação e integrado sobre o domínio íl:

T . <jJ. (r) s(r, t) díl -1 -

i = 1, 2 ... n ( 2. 8)

Igualando-se essas n integrais a zero obtêm-se

n equações diferenciais ordinárias nas coordenadas

forma:

n i::

j=l nlq,~ - L'.' iJ

díl q.(t) + J

<j,~ F(r, t) díl -1

ou ainda:

n i::

j=l

i

~~ M ~. díl ij.(t) = -1 - J J

1, 2 ... n

q. ( t) J

na

(2 • 9)

9

g ( t) (2.10)

Os coeficientes das matrizes seriam:

m .. = l <P :' M <P • díl iJ -1 ., J

(2.lla)

k .. = 1 <P:' D~~ díl 1J -1 - -.L_

íl

(2.llb)

Q. = l <P :' ~e::: ' t) díl 1 -1

( 2. llc)

Observe-se que as matrizes obtidas serao

cas no caso de operadores auto-adjuntos [s]

simétri-

As funções <j). foram escolhidas de maneira a aten­-1

der as condições de contorno associadas ao problema, e, conseque~

temente, os resíduos foram construídos somente sobre o domínio íl.

Essas condições podem ser relaxadas [7] adotando-se funções <P • -1

que somente satisfaçam às condições de contorno geométricas, des­

de que se considerem resíduos ponderados tambêm no contorno e que

sejam adicionados ã expressão (2.8).

b) - Método_da_Coloca5ão -

O método da colocação admite inicialmente uma elas

sificação de acordo com a natureza das funções de aproximação <P • -1

i) Método do Contorno - quando as funções <P • -1

satisfazem somen

te ã equação diferencial

ii) Método do Interior - quando as funções <P • -1

satisfazem as con

iii) Método Misto

diçÕes de contorno mas não satisfazem à

equaçao diferencial

quando as funções <f:>. nao satisfazem nem a -1

equaçao diferencial nem às condições de con -

torno.

10

O método consiste em se tomar pontos no contorno e/

ou no domínio, de acordo com o caso i), ii) ou iii), e igualar a

zero, em cada ponto, uma ou mais componentes do erro E ( r. , t) g~ - ;_, 1

rando um total de n equações diferenciais o~di~ãrias nas coorde

nadas q. (t) que, considerando o e.aso ii), se escrevem: J

~ D.G<j,.(r.)~ . 1-1 -J -1 J=

q.(t) + J

n l: M.(r.)

. -1 -1 J=l

<j, . ( r. ) -J -1

êj.(t) J

F.(r.,t) 1 -1

i=l,2 n

(2.12)

D. e M. so contem a linha correspondente ã componente de E considerada e -1 -1.

o elemento de F(r. , t) - -1

correspondente e denotado por

matricial:

onde os coeficientes sao:

m •. = Mi ( r. ) <j, . ( r. ) 1) - -1 -J -1

k .. =Di ÍcP.(r.0 1J _ cJ -1'.J

Q. =F. (r. ,t) 1 1. -1.

F.(r.,t). 1 -1

Em forma

( 2. 13)

( 2. 14a)

(2.14b)

( 2.14c)

De imediato nota-se que os coeficientes acima sao

obtidos diretamente, sem integração como no método de Galerkin

Por outro lado nao se pode garantir a simetria das matrizes de

massa e rigidez o que, alem de complicar o problema de auto-valor,

prejudica a eficiência dos métodos numéricos de resolução das equa­

çoes diferenciais ordinárias resultantesª

2.3.2 - Métodos Variacionais - Método de Ritz

A característica principal dos métodos variacionais

11

e nao operar diretamente sobre as equaçoes diferenciais de movi­

mento (2.4), como fazem os mêtodos residuais, e sim sobre uma for

ma variacional equivalente (2.2).

O mais conhecido desses mêtodos e o mêtodo de Rit&

No mêtodo de Ritz a aproximação

u -n

n E

i=l

u -n

<P. ( r) -1 -

Dessa maneira a solução aproximada

adotada ê da forma:

q. ( t) 1

u -n

estarã, em cada

(2.15)

instante

t restrita ao sub-espaço s n

de dimensão n gerado pelas fun-

çoes <j).(i=l,2 ... n) -1

escolhidas e nao no espaço de dimensão

infinita em que o problema foi originalmente proposto. A aproxi-

maçao u -n

é levada ao princípio variacional, no caso o princípio

de Hamilton (2.2), e sao obtidas equaçoes diferenciais ordinárias

nas coordenadas generalizadas q. ( t) 1

As funções <P. -1

sao linearmente independentes e d~

vem atender à condição de admissibilidade e formarem um conjunto

completo para se garantir a convergência do processo a medida que

n cresce e

A condição de admissibilidade implica em que as

<P. -1

devam ter derivadas contínuas até ordem m - 1 on-. funções

de m e a maior ordem de derivação que surge no funcional (2.2),

e satisfaçam às condições de contorno geomêtricas.

O conjunto de funções <j).(i=l,2 -1

n)

completo [a] num espaço s quando para qualquer elemento

tencente a s existe uma combinação u -n

n E

i=l a. <P .

1 -1 que

e dito

u P eE_

conver-

ge para u , na norma de S , quando n tende para infinito.

12

Kitodo dos Elementos Finitos

Ati aqui todos os mitodos faziam uso da aproxima-

çao (2.6) onde as funções cj). ( r) -1 -

alim de atender pelo menos as

condições de contorno geométricas, deveriam se extender por

o domínio.

todo

A construçao das funções cj). -1

torna-se frequente-

mente impraticável em problemas práticos onde elas devem tambim

representar adequadamente contornos irregulares e propriedades ma

teriais complexas.

A idiia de se adotar funções de aproximaçao que se

anulam em todo o domínio exceto em sub-regiões leva ao mitodo dos

elementos finitos [7]

No mitodo dos elementos finitos o domínio i dividi

do em sub-domínios, ou elementos, de formato conveniente, defini­

dos por um conjunto de pontos desse sub-domínio, chamados nos.

Sobre cada elemento sao assumidas as funções de

aproximaçao, ·em geral na forma de polin~mios simples, que se anu-

1am fora do elemento considerado.

Alim disso os deslocamentos (e eventualmente suas

derivadas em relação a variâveis espaciais) são tomadas como coar

denadas generalizadas que passam assim a ter significado físico

claro.

Dessa maneira as funções de aproximação se tornam

funções de interpolação tendo valor unitário no nó correspondente

e nulo em todos os demais.

A maior vantagem [9] da discretização pelo mitodo

dos elementos finitos reside no fato de se poder aproximar domí-

nios arbitrários com propriedades físicas complexas usando um con

13

junto de elementos de formas e propriedades simples.

A introdução das condições de contorno geometricas,

feita apos a montagem do sistema de equaçoes, e facilitada e e

possível alto grau de automatização do metodo. Para se obter uma

melhor aproximaçao pode-se proceder agora de duas formas:

a) Aumenta-se o numero de elementos em que o domínio e dividido

mantendo-se as funções de interpolação em cada elemento.

b) Mantem-se o numero de elementos em que o domínio e dividido e

refinam-se as funções de interpolação em cada elemento.

Observe-se que fora do metodo dos elementos

tos a Única alternativa·e aumentar o numero de funções

nidas em todo o domínio, na serie (2.6).

fini-

defi

14

III - MÉTODOS DE SOLUÇÃO DAS EQUAÇÕES DE MOVIMENTO DO MODELO DIS­

CRETO

3.1 - INTRODUÇÃO

A solução das equaçoes de movimento do modelo dis-

ereto, formuladas no capitulo anterior, poderia ser obtida por

qualquer processo numérico clássico para a resolução de problemas

de valor inicial.

A prática mostrou entretanto que na maioria dos

problemas de dinâmica estrutural tais processos sao ineficientes.

Os processos mais usados em dinâmica estrutural sao

a) superposiçao modal

b) integração direta

As características dos dois processos serao comen-

tadas neste capitulo. Antes porem serao definidas as principais

propriedades do modelo discreto.

3.2 - CARACTERÍSTICAS DO MODELO DISCRETO

A aplicação dos metodos de discretização espacial

vistos até aqui, conduz sempre a um sistema de equações diferenciais

ordinárias que pode ser colocado na forma:

M u + R(u) F ( t) (3.la)

para sistemas nao lineares, e na forma:

M Ü + K u = F(t) (3.lb)

para sistemas lineares.

15

As matrizes M e ~, de massa e rigidez respecti

vamente, e o vetor F de ações externas, são obtidos de acordo

com o processo adotado. Note-se que o emprego do método de Ritz

ou Galerkin, em conjunto com o método dos elementos finitos, faz

com que essas matrizes sejam simétricas e em banda o que é muito

vantajoso do ponto de vista computacional.

As características principais do modelo

surgem ao se procurar soluções da forma:

u(t) X f(t)

para o problema de vibrações livres:

M u + K u o

Substituindo (3.2) em (3.3) obtém-se:

M x f(t) + K x f(t) = 0

-ou seja, n equaçoes:

n .. n

discreto

e 3. 2 l

e 3. 3 l

e 3. 4 l

l: m .. X. f(t) + l: k .. X. f(t) o i=l,2 ... n j=l lJ J j=l lJ J

e 3. s l

nas _quais podem ser separadas as variiveis:

n l: k .. x . ..

f(t) j=l lJ J - f(t)

i 1 , 2 ... n (3. 6) n l: m .. x.

j=l lJ J

16

Como a expressao ã esquerda da igualdade (3;6) in-

depende do Índice i e a expressao ã direita independe do tempo,

ambas devem ser constantes e assim pode-se escrever:

n I

j=l

f(t) + w2 f(t)

(k .. - w2 m •• ) iJ iJ

o ( 3. 7)

i = 1, 2 .... n ( 3. 8)

Ao se escolher uma constante positiva w2 obtém-se

soluções harmônicas na equação (3.7) da forma:

f ( t) A cos (wt - <j)) ( 3. 9)

onde A e sao constantes arbitrárias e os valores que w p~

de assumir sao dados pelas equaçoes (3.8) que constituem um pro­

blema de auto-valor que se escreve matricialmente:

(K - w2 M) x o (3.10)

Para matrizes M e K simétricas e positivas definidas sao en-

contrados n valores para w2 e as respectivas raízes quadradas

positivas sao as chamadas frequências naturais do sistema (3.lb).

A cada frequência natural w está associado um au

to-vetor ~ , obtido a menos de uma constante multiplicativa, cha

mado modo normal de vibração, que satisfaz a equação (3.10).

Chama-se de matriz modal X do sistema a matriz

cujas colunas sao os auto-vetores do sistema ordenados ,segundo as

respectivas frequências, de maneira crescente.

Os modos normais podem ser normalizados em relação

17

a matriz de massa pela relação:

T X M X 1 (3.11)

que determina de maneira Única o mÜdulo das componentes de x

A propriedade mais importante dos modos normais e

a sua ortogonalidade em relação às matrizes de massa e rigidez:

I (3.12a)

(3.12b)

onde I e a matriz identidade de ordem n e w e a matriz dia-

gonal contendo as frequências naturais do sistema.

Observe-se qti~ enquanto o sistem~ contínuo aprese~

ta um espectro de frequências discreto que parte da fundamental

(a mais baixa) e cresce indefinidamente, o modelo discreto (3.lb)

que o aproxima tem um espectro limitado.

·A frequência mais alta desse espectro e chamada

frequência de corte (cut-off frequency) e representa o limite su­

perior de excitação do modelo discreto; excitações em frequências

mais altas produzem ruído ou resposta espacial aleatória [1Q]

Outra característica comum aos modelos discretos

[li] e a tendência de deterioração da resposta à medida que as

frequências de excitação se aproximam da frequência de corte.

Isto se deve aos erros no espectro de frequências

e nos modos normais do sistema discreto, que são maiores nas fre­

quências e modos mais altos.

18

Quanto a solução do problema da resposta

e de interesse identificar quatro tipos de solução [l(U

dinâmica

1) A solução exata das equaçoes de movimento do sistema contínuo

2) A solução exata das equaçoes de movimento do modelo discreto

(3.1) conservando número ilimitado de dígitos

aritméticas.

nas operaçoes

3) A solução exata das equaçoes (3.1) porem com numero limitado de

digitas nas operações aritméticas.

4) A solução aproximada das equações (3.1) também com numero limi

tado de dígitos nas operações aritméticas.

Define-se entao erro de truncamento como aquele ca~

sado pela substituição de um sistema contínuo por um modelo

ereto e é a diferença entre as soluções 1 e 2

dis-

Erro de arredondamento (round-off errar) depende do

grau de precisão dos cálculos aritméticos e é a diferença entre

as soluções 2 e 3. A diferença das soluções 3 e 4 dependerá do

algorítmo de integração empregado; é o erro de integração.

Será vista a seguir a caracterização do amorteci -

menta.

3.2.1 - Introdução do Amortecimento

Todo sistema dinâmico real apresenta alguma forma

de dissipação de energia. A definição desses mecanismos de dissi

paçao de energia é muito difícil e normalmente são empregados mo-

delas matemáticos convenientes para sua representaçao .de forma

aproximada. são citados a seguir três desses modelos.

Chama-se de amortecimento de Coulomb ou atrito se­

co aquele que dá origem a uma força de módulo constante e que se

opõe sempre ao movimento.

19

Amortecimento estrutural [s1 e outro modelo também

usado e estã relacionado com a energia dissipada

material da estrutura.

internamente no

Entretanto, o modelo mais usado, devido a sua rela

tiva comodidade, é o amortecimento viscoso linear que dâ origem a

uma força proporcional a velocidade e de sentido oposto a ela.

Neste caso, introduzindo-se o amortecimento direta

mente no modelo discreto, obtém-se:

M u +Cu+ K u = F(t) (3.13)

No Capítulo IV formas usuais de construçao da ma-

triz de amortecimento e serao abordadas.

3.3 - SUPERPOSIÇÃO MODAL

O processo de superposiçao modal [s, 6] e adequado

ao caso em que as matrizes M , C e K sao constantes no tempo - -e busca o desacoplamento do sistema de equaçoes (3.13) através de

uma transformação linear conveniente.

Sua eficiência estã ligada à forma da .matriz e

que dá origem a dois procedimentos distintos conforme o amorteci­

mento seja do tipo proporcional ou não-proporcional.

a) Amortecimento_ProEorcional

Ocorre quando os modos normais sao ortogonais a ma

triz de amortecimento e a transformação linear procurada é:

li ( t) X 11(t) (3.14)

20

onde nCtl sao as novas coordenadas generalizadas e X e a ma-

triz modal do sistema.

Introduzindo-se a transformação (3.14) nas equa-

çoes de movimento (3.13) e levando em conta as relaç~es de ortog~

nalidade (3.12) e a nova relação:

onde

com:

C*

C*

e diagonal com seus elementos dados por:

e~. 11

2 i;. w. 1 1

w. = frequência natural do modo i 1

i; . 1

fração do amortecimento crítico do modo i

as equaçoes se desacoplam na forma:

n. + 21;. w. n. + w~ n. = N. 1 1 1 1 1 1 1

i 1, 2, ... n

(3.15)

(3.16)

(3.17)

onde os N. 1

sao elementos do vetor de forças generalizadas:

N ( t) (3.18)

A solução das equaçoes (3.17) pode ser obtida pela

transformada de Laplace. No caso mais comum de ~. < 1 1

temos:

21

t

11.(t) 1

J Ni (T) exp [- çiwi (t - T~ sen wid (t - T) dT +

o

com:

ç . 1

1 (1 - ç;:)2

1

+ t_!_ exp(- ç.w. t) sen w.d J n1.(0) ~id 1 1 1 ~

11 (O)

l w.(1-ç:)2

1 1

XTMu(O)

n(O) = XT M ~(O)

(3 .19)

(3.20)

onde a integral de convolução deverá, em geral, ser calculada nu-

mericamente.

Uma alternativa é aplicar sobre as equações (3.17)

um algoritmo de integração direta como os que

adiante.

-serao vistos

A solução final é dada pela superposição das

postas correspondentes a cada modo através de (3.14).

-mais

res-

Dessa maneira, o processo de superposição modal en

volve a resolução de um problema de autovalor para a determinação

de uma nova base para o espaço dos vetores deslocamentos onde as

equaçoes se apresentem desacopladas para, em seguida, ser feita a

integração no tempo.

22

A grande vantagem dessa mudança de base e que,para

virios tipos de problema,ê possivel se obter uma boa

çao para a resposta considerando-se apenas uma pequena

aproxima­

fração do

número total de equações desacopladas, normalmente aquelas corre~

pondentes aos primeiros modos. Isto permite também que apenas os

primeiros autovalores e autovetores sejam calculados.

A resposta ê truncada, só se considerando as con-

tribuiçÕes dos p primeiros modos onde p ê determinado pelas carac

terísticas da·estrutura, pela distribuição espacial do carregame~

to e pelo seu conteúdo de frequências.

b) Amortecimento_Não-ProEorcional

Neste caso a matriz C* da equaçao (3.15) nao e

diagonal e as equações de movimento na base formada pelos modos

normais sao ainda acopladas pela matriz de amortecimento C*

Entretanto, sempre que a resposta estiver essenci-

almente contida nos p primeiros modos, o processo mais eficien-

te [12) parece ser integrar diretamente o sistema truncado forma-

do pelas p primeiras equações na base dos modos normais:

N < t) (3.21)

onde a matriz C* e obtida a partir de C* desprezando-se o ac~

plamento entre os p primeiros modos considerados na resposta e

os restantes (n - p) modos superiores.

Em (3.21) as matrizes sao de ordem (p X p) e os

vetores de ordem p Uma maneira de desacoplar o sistema consis

te em inicialmente transformar o sistema de n equações diferen-

ciais de 2~ ordem num sistema de 2n equaçoes diferenciais de 1~

23

ordem como mostrado em [s]

A partir dai a solução e obtida de maneira seme-

lhante ao caso não-amortecido, envolvendo um problema de auto-va-

lorde ordem 2n e a determinação da matriz modal associada

desacopla o sistema. Resultam 2n equaçoes diferenciais de

que

1 !!

ordem desacopladas cuja solução pode ser obtida por meio da trans

formada de Laplace. Entretando o fato de se dobrar a ·.·ordem do

problema de autovalor e o uso de aritmética complexa acarretam um

aumento de esforço computacional que,normalmente,não

vantagem de se ter equaçoes desacopladas.

compensa a

Nesse caso e preferível a integração direta de to-

das as equaçoes (3.13) ou a integração das

nas coordenadas modais, equaçao (3.21).

p primeiras equaçoes

Qualquer que seja o tipo de amortecimento, o pro-

cesso de superposição modal exige que as matrizes ~ , ~ e K se

jam constantes no tempo o que limita o seu uso a problemas linea-

res, embora em [13] sua extensao ao caso não-linear

proposta.

tenha sido

Note-se que o Único erro inerente ao processo é de

truncamento, quando as contribuiçiies dos modos mais al·tos sao des

prezadas. Pelo visto a maior eficácia da superposição modal se

dâ em problemas lineares cuja resposta pode ser bem representada

pela superposição das respostas correspondentes a uma pequena fra

ção do número total de modos da estrutura, geralmente os mais bai

xos.

24

3.4 - INTEGRAÇÃO DIRETA

Os métodos de integração direta recebem esta deno-

minaçao por não efetuarem nenhuma transformação prévia nas

-çÕes de movimento.

O sistema de equaçoes em sua forma original,

equa-

aco-

plada, e integrado numericamente por um procedimento passo a pas­

so a partir do conhecimento das condições iniciais do problema.

A instituição dos algorÍtmos de integração direta

pode ser feita a partir de pelo menos tris ideias bisicas vistas

a seguir.

A primeira idéia consiste em se discretizar o ope­

rador diferencial com o uso de fórmulas de diferenças finitas co­

mo no caso do método da diferença central [14]

Uma outra maneira consiste em assumir formas deva

riaçao para as acelerações, velocidades e deslocamentos ao longo

do tempo em cada intervalo de integração e satisfazer as equaçoes

de movimento nos pontos extremos desses intervalos de tempo. É o

caso do método da aceleração linear e o método - 8 de Wilson [15].

Um procedimento alternativo seria assumir formas de variaçao para

acelerações, velocidades e deslocamentos ao longo do intervalo de

integração e aplicar um critério de aproximação como o de Galer-

kin [6, 16]

Embora a instituição dos métodos possa ser distin­

ta chega-se sempre a uma relação de recorrência que permite o ava~

ço da solução no tempo e pode geralmente ser escrita na forma:

D u(t -o - n t n-1

••• ) + p (3.22)

com:

~ ( t 1) - n+

onde D -o

e ~l

u ( t 1) - n+

e

25

~(t ,t 1, ... ) = - n n-

u ( t ) - n

sao matrizes que definem a transformação dos ve-

tores de aceleração, velocidade e deslocamento nos tempos

t n

t n-1

mento do tempo

çao externa.

nos vetores de aceleração, velocidade e desloca-

t n+l p

As matrizes

e um vetor que leva em conta a excita-

D -o

e ~l sao funções das proprieda-

des do sistema (matrizes de massa, amortecimento e rigidez) e do

intervalo de integração ti t

Supondo-se que seja inversível, a equação(3.22)

se escreve:

ou ainda:

onde A = D-l D -1 -o

D-l D u(t -1 -o - n

= A ;;(t n t l' n-

t n-1 -1

••• ) + ~l p

-1 ••• ) + ~l p

e o operador de integração.

(3.23a)

(3.23b)

Quando as matrizes de massa, amortecimento erigi­

dez nao sao funções dos deslocamentos ou de suas derivadas a trans

formação é dita linear.

Se o vetor u do 29 membro da relação (3.22) con-

26

têm informações apenas do instante anterior t o método e n

dito

de um passo (one-step); se contêm informações de tempos anterio-

res ê chamado de múltiplos passos (multi-step).

Sempre que a matriz ~l puder ser colocada em for

ma triangular o algoritmo e dito explícito; caso contrário implÍ-

cito.

Outra classificação importante surge a partir do

estudo da estabilidade do operador de integração [16, 17, 18, 19]

e a esse respeito serio ·feitas a seguir algumas observaç~es base~

das na referência [16] onde o problema ê abordado de maneira mais

rigorosa e completa.

Um algorítmo de integração é dito estável se todos

os operadores do conjunto An (n = 1, 2, ... ) forem uniformemen-

te limitados, ou seja, existem constantes positivas K n

li An V li < K li V li n

onde li v li e a norma de v

Definindo raio espectral de A [18] como:

p = max 1 ". 1 1

onde À. 1

sao os autovalores de A , a condição anterior

exige que p seja limitado pela unidade.

Pode-se entao resumir o ~omportamento

aproximada de um oscilador simples como se segue:

da

- se

- se

p > 1

p < 1

a solução e amplificada ao longo do tempo;

a solução e amortecida ao longo do tempo

tal que:

(3.24)

(3.25)

(3.24)

solução

27

- se p = 1 nao hã introdução de erro nas amplitudes da solução

numérica obtida.

Um operador é dito incondicionalmente estável se a

solução para quaisquer condições iniciais nao cresce indefinida

mente, qualquer que seja o valor do intervalo de integração usado

Quando isto so ocorre para valores de IH menores que um

determinado valor crítico, o operador é dito condicionalmente es-

tãvel e se não ocorre para qualquer At > O

condicionalmente instável.

o processo e dito in

Os operadores de integração introduzem, normalmen­

te, dois tipos de erro na solução do problema oscilatório: um er­

ro nas amplitudes e um erro no período.

Quando p < 1 hã introdução de um amortecimento

artificial na solução numérica e a amplitude do deslocamento de-

cai com o tempo. Com relação ao período hã, em geral, um aumen-

to do mesmo na solução numérica.

Para a representaçao da solução oscílatÕria os auto

valores de A

onde:

devem ser complexos conjugados da forma:

a ± bi

<j, are tg(b/a)

±i<j, p e (3.26)

(3.27à)

(3.27b)

Podem entao ser tomados como medidas dos erros introduzidos na so

lução numérica os parâmetros [20]

28

amortecimento artificial do algoritmo

T - T T

w 6t _ l -- l . ~ t erro reativo no periodo

onde T e T sio os períodos do oscilador e da soluçio numérica

aproximada, respectivamente.

São comentados a seguir dois operadores muito usa­

dos na prática e de características distintas.

a) Método_da_Diferensa_Central_ll41

Adota para acelerações e velocidades as fÕrmulas de

diferenças finitas:

Ü(t) 1

t(t - 6 t) 2 u(t) + u(t 6 t ;l = - + 6 t 2 J

(3.28a)

;, ( t) 1 t u(t 6 t) u(t + 6t~ = rn - + (3.28b)

que levadas a equaçao de movimento no tempo t

M Ü(t) + e {,(t) + K u(t) = F(t) (3.29)

fornecem a relaçio de recorrência do método:

e: 2 M - 2 lút ~ u ( t - 6 t) (3.30)

onde se usa para a inicialização do processo:

29

u ( - lu) = u (o) - ln ~ (o) + li~ 2

ü (o) (3.31)

-Da expressao (3.30) conclui-se que o uso de uma

matriz de massa diagonal e, no caso de problemas com amorteCimen

to, uma matriz de amortecimento também diagonal, tornam o método

explÍcito,não havendo ·necessidade de se resolver um sistema de

equações simultâneas e nem de se montar as matrizes globais da es

trutura, podendo-se trabalhar a nível de elemento [14}.

A desvantagem do método reside no fato de ser ape-

nas condicionalmente estâvel, exigindo intervalos

tit<tit .· - cr1.t.

No caso linear:

tit<i'it. cr1t.

= T .

m1n 'IT

onde T . e o período mínimo do modelo discreto. m1n

de integração

(2.32)

Tem sido aplicado para a obtenção de resposta tra~

siente em problemas lineàres [21, 22] e não lineares [23, 24].

b) Metodo_de_Newmark_l2s1

O método pode ser visto como uma extensao do méto­

do da aceleração linear e adota as expressões:

~<t + titl = i'.ict) + [o -yl · ~<tl + y ~ct + tit~ tit (3.33a)

u(t + i'it) = 2(t) + i'it ~(t) + [<-}- 13) Ü(t) + 13 Ü(t + tit~ tit 2

(3.33b)

onde y e 13 são parâmetros que definem a estabilidade e precisão

30

do algoritmo.

Com a equaçao. de movimento, tomada no tempo t + 1:it,

e as expressoes (3.33) chega-se ã relação de recorrência do méto-

do:

~+

+

+

onde:

a o

- y

- y

a o

M

e

1

1 2

1 2

M + -

~o

~l

e

e

ª1 ~ u(t

~ ( t)

u(t)

+ ª2

+ ª4

_Y_

s /:it

+ ti t)

;, e t) + ª3

;, ( t) + ªs

1

s 1:it

f(t + l:it) +

~(t~ +

~(t~ {3.34)

1 ª3 = - - 1

2S

Dois casos de interesse sao:

s

s

1 6

1 4

corresponde ao método da ac·eleração linear

e o método da aceleração média constante.

O primeiro e condicionalmente estável e o

incondicionalmente estável.

segundo

Na equaçao (3.34) observa-se que o método é implÍ-

cito e que em problemas onde K não é constante é necessário efe

tuar uma nova triangulação a cada intervalo de tempo.

Uma aplicação a problemas não-lineares e apresent~

da em [26].Além dos métodos de Newmark e de diferença central,são

também comuns aplicações do método - e de Wilson [27] e do méto-

31

do de Houboult [28]

De uma maneira geral pode-se dizer que os métodos

de integração direta não fazem restrições as matrizes de massa,

amortecimento e rigidez do sistema, sendo assim indicados para os

problemas não lineares.

Nos problemas lineares a integração direta e usada

naquelas situações em que a superposição modal perde suas vanta -

gens como no caso de amortecimento não-proporcional e resposta ri

caem frequências.

3.5 - ALGORÍTMOS BASEADOS NA INTERPOLAÇÃO DAS FORÇAS DE INÉRCIA

A equaçao de movimento do modelo discreto pode ser

escrita na forma:

R = M u !'. ( t) - C u - ~E (3. 35)

onde:

R e o vetor das forças de inércia

~E e o vetor das forças elásticas de restituição

O problema serã discretizado agora ao longo do te~

po adotando-se intervalos de tempo ,=t.1-t. i+ i

lo pode ser considerado um elemento finito de tempo

qual se interpolará o vetor de forças de inércia R

Esse interva

ao longo do

A sequencia de elementos de tempo forma o interva­

·lo de tempo total a ser estudado.

O vetor das forças de inércia R sera interpolado

32

por polinômios de Hermite a partir dos valores assumidos por R e

suas derivadas em relação ao tempo nos pontos extremos d? interva

lo T

Assim sao obtidas interpolações lineares, cúbicas,

etc., das quais as duas primeiras serão comentadas a seguir.

a) Intereolação_Linear

Neste caso R se determina a partir de seus valo-

res nos extremos do intervalo (elemento) de tempo

de interpolação:

com:

ou:

onde os Índices

1 - E; e

t T

A interpolação se escreve:

R

R

o e 1 correspondem ao início

pelas

(E; = O)

nal (E;= 1), respectivamente, do intervalo de tempo T

b) Intereolação_CÚbica

funções

(3.36)

(3.37a)

(3.37b)

e fi-

Agora as funções de interpolação sao os polinômios

de Hermite de 39 grau:

33

1 - 3 1; 2 + 2 1; 3 1;-21;2+1;' </>3 = 3 1;2 - 2 t;'

e a interpolação se escreve:

R (3.38a)

ou:

R (3.38b)

Já que 1

di; = d t . T

-Integrando a equaçao de movimento:

Mu=R(t) (3.39)

de t a um t genérico dentro do intervalo T obtém-se: o

i;

M u = M u + J R T di; - -o (3.40)

o

Integrando-se (3.40) chega-se a:

i; i;

M u ~ ~o + i; T ~ ~o + J f R T 2

di;di; (3. 41)

o o

-As equaçoes (3.39), (3.40) e (3.41) em conjunto

com as interpolações (3.37) ou (3.38) resultam em dois algoritmos

de integração:

34

a) Intereolasão_Linear

Introduzindo-se (3.37b) em (3.40) e (3.41), efetu

ando-se as integrações e fazendo-se ç = 1 obtém-se:

ou ainda:

M ~1 = M u - -o

M '.:1 M u + T M u - -O - -O

T '.:1 = u + 2 -o

'.:1 = u + T u + -o -o

T2 + -3 R -o

(Ü ~l) + -o

T2 .. T2 T u + 6 -o ~l

- -

(3.42a)

(3.42b)

(3.43a)

(3.43b)

Observe-se que as expressoes acima sao as mesmas

que se obtém fazendo

método de Newmark.

1 y = 2 e

1 8 = 6 nas expressoes (3.33) do

Assim a interpolação linear das forças de inér

eia conduz ao método das acelerações lineares, como era de se es-

perar.

Neste trabalho apenas o caso da interpolação cúbi­

ca sera estudado.

b) Intereola~ão_CÚbica

Seguindo-se o mesmo procedimento do caso anterior

obtém-se as -expressoes:

. M u - -o

T + TI <6 ~º + T R

-O (3.44a)

• Mu +TMu - -o - -o

.T 2 • • + 6Õ ( 21 ~o + 3 ·T !l:o + 9 !l:i - 2 T !l:i)

(3.44b)

ou ainda:

~l = u -o

'( -1 + TI ~

~l '( 2 -1

= u + T u + M -o -o 60 -

35

(6 R -o

+ '( . R -o

. (21 R + 3 T R -o -o

(3, 45a)

(3.45b)

Este algoritmo assume para as acelerações uma va­

riação cúbica, superior portanto, a dos métodos usuais.

Observe-se que as equaçoes (3.44) ou (3.45) nao p~

dem ser resolvidas diretamente por nao se conhecer o vetor das for

ças de inércia ;1

tempo.

. e sua derivada ; 1

ao final do intervalo de

Os processos de solução das equaçoes (3.44) depen-

derão da forma em que os vetores ;l puderem ser coloca -

dos.

É interessante portanto nesse ponto se separar os

problemas lineares e os não lineares.

i) Problemas Lineares

Nesse caso e possível se escrever:

;l = !1 - e ~l - K ~l (3.46a)

. ~l = !1 - e ~l - K ~l (3.46b)

onde ~l e o vetor das excitações externas ao final do intervalo

de tempo considerado.

A introdução das expressoes (3.46) nas equaçoes

36

(3.45), levando em conta que:

~l -1

= M ~l (3.47)

resulta num sistema de equaçoes que pode ser escrito na forma:

onde os vetores:

r = -o

u -o

u -o

D r + P -o ..... o

e :1

. ~ . l.Ill.Cl.O representam o estado do sistema no

~l

~l

e no final,

mente, do intervalo de tempo e as demais matrizes são

por:

(3.48)

(3.49)

respectiva-

definidas

, 9 2 2 3 ( C M-1 60 T C - 60 T K - º)

~l =

D = -o

.I K T2 -1 2-+12~·'.'.! K

T T2 1 - - K +~CM- K

2 - 12 - -

(3. 50a)

T T2 C M-1 M - 2 ~ - 12 (K - C)

(3.50b)

P=

~ I 2 -

T 2 1 - CM-12 - -

37

. 2 1 ~ I + !..._ C M-2 - 12 - -

F -o

+

(3.50c)

e I e a matriz identidade de ordem n.

-çao sucessiva

A solução ê obtida ao longo do tempo pela aplica-

-da equaçao (3.48) que,

r -o

colocada na forma:

-1 + !'.1 p (3. 51)

só envolve multiplicações matriciais apos a inversão de !'.i . Ob-

serve-se que as matrizes em (3.48) -sao de ordem 2n, onde n e o

número de grãus de liberdade do modelo discreto,e que alêm disso

nao sao simétricas.

Além dessas duas desvantagens existe ainda o pro­

blema da estabilidade do algoritmo.

A aplicação das equações (3.45) a um sistema com

um grau de liberdade

rador de integração é:

(n 1), sem amortecimento, mostra que o op~~

1 -13 2 2 1 wt+Tti. 30 W T + 8Õ

1 + J:._ W2T2 1 4 4

15 + 240 W T

A =

2 1 2 2 - W T(l - ·lO W T)

l + J:._ W2T2 1 4 4 15 + 240 W T

T(l _ J:._ W2T2 + 1 4 4 . 10 720 W T

1 1322 1 - 30 W T + - W4T4 80

(3.52)

38

onde w e a frequência natural do oscilador.

Os autovalores À da matriz A da equaçao (3.52)

sao complexos e de mÕdulo unitário podendo ser colocados na forma:

com:

1 WT (1 - W2T2

cj) 10 tg

(1 -

desde que se adote:

À e ±icj)

1 W4 T 4 )! + 720

13 W2T2 1 3Õ + 80

T<0.503T o

onde T o

e o período do oscilador.

(3.53)

1 .1

(1 - W2T2)2 10

(3. 54)

W 4 T 4 )

(3.55)

A condição (3.55) caracteriza o algoritmo definido

pela interpolação cúbica das forças de inércia como apenas condi­

cionalmente estável.

Para sistemas com vários graus de liberdade, T o

deve ser substituído na expressao (3.55) por T . , período mini min

mo do sistema, o que acarreta, normalmente, intervalos de integr~

çao desnecessariamente pequenos.

É possível, entretanto, modificando ligeiramente

os coeficientes das expresssÕes (3.45), tornar o algoritmo incon­

dicionalmente estável como mostrado em [2, 3].

As novas expressoes, correspondentes as expressoes

(3.45) sao:

~l u -o

~l

+ T U -o

u -o

39

+ T M-l (6 R + T R0

6 R R0

)

12 - -o -o+ -1 - T -1

(20 R -o

. + 2.5 T R

-o

(3.56a)

(3.56b)

Note-se que apenas a fÕrmula para os deslocamentos

teve seus coeficientes alterados e assim o algoritmo continua de­

terminando exatamente a velocidade de uma partícula submetida a

uma força resultante cúbica no tempo. Devido à alteração dos coe

ficientes de (3.45b) para (3.56b) os deslocamentos sao exatos ap~

nas para forças quadráticas no tempo.

Substituindo -as expressoes (3. 45) pelas -expressoes

(3.56) as matrizes na equação (3.51) passam a ser:

T 2 T 3 -1 6 e - 24 (K - e M C)

T T 2 1 K + C M- K

2 - 12 T 2 -1 12

(K - CM C)

(3.57a)

T 2 T 3 - C M-1 T M - } C - 24 (! C)

T T2 -1

M - z ~ 12 (K - CM C)

(3.57b)

40

-T2 T' 2 T' T 3 • -1 ~I -1 - F ) -I 24 .<: M + 24 _<; M F 24 (~º 3 - - 6 - - -o -1 p = +

2 T2 -1 T2 _!_ I - ~ C M-l _!_ I + TI_<; M ~l 12 <F - F ) 2 - 12 - - 2 - - -o -1

(3.57c)

O operador de integração para o caso de um grau de

liberdade se escreve agora:

A 1

(3.58)

Os autovalores da matriz A são complexos e de mi dulo unitário para qualquer intervalo de integração

T adotado e ~ e tal que

tg ~ (3.59) 5 2 2 W

4T

4

(1 - l2 w T + 144)

Dessa maneira o intervalo de integração passa a

ser escolhido de acordo com a distribuição espacial do carregame~

to e seu conteúdo de frequências, sem a preocupação com a instabi

lidade do algoritmo na avaliação da resposta correspondente aos

modos mais altos4

Note-se que o algoritmo obtido com a interpolação

linear das forças de inércia e também apenas condicionalmente es-

tâvel: T < 0.551 T o

(Newmark 1 y = 2 e

41

ff interessante notar que é possfvel modifici-lo de

maneira a torná-lo incondicionalmente estável e o algoritmo resu!

Cante [3] e o método da aceleração média constante da famflia de

Newmark, equaçoes (3.33), com 1

y = 2 e 1 s = 4 Como alternati-

va de solução das equações obtidas ·a partir da interpolação ciibi­

ca das forças de inércia (3.45) ou (3.56), que não dobre o niimero

de equações do problema e trabalhe com matrizes simétricas, esti

a solução iterativa que e usada nos problemas não lineares e sera

vista adiante.

A partir da definição de erro relativo no período

é possfvel traçar curvas (Fig. 3.1) que relacionem esse erro (ex-

pressa em percentagem) com o parâmetro T/T o

para os algoritmos

apresentados.

A fig. 3.1 foi retirada da referência [2]

Aumento no

Perto d o ('%)

7

6

5

4

3

2

o

Wilson e=1.4/ /

I I . '

.' / Newmark

/ / V:o.5 I , ·

1 l' '0.25 .' ,'

/1 • I /1 ·/ /1 •/ / 1 ., ,, ·/ t,

/ /

/

/

/ .. I

., ,, i'

__ ,,,, .. .. -./

0.1 0.2 0.3 0.4

FIG. 3.1

/

/ I

/ Alo, lncond. ·Estovel. Eq.( 3.58)

Al9. Cond. Esravel. Eq. < 3. 52 >

0.5

42

Devido à não-linearidade na relação entre as for-

ças elásticas ~E

se escrevera:

e os deslocamentos u a equação de movimento

R = M u = F(t) - Cu - ~E (3.60)

onde:

(3.61)

Derivando (3.60) em relação ao tempo obtém-se:

. R M u F(t) - eu - K(u) u (3.62)

j â que:

. ~E

u K(u) u (3.63)

Linearizando (3.61) no intervalo de tempo T obtém-se:

~El

K -o

(3.64a)

': 1 (3.64b)

Observe-se que em (3.64b) -e mais consistente usar-se a matriz de

de rigidez tangente avaliada em ':i : ~l = ~(~ 1 ) no lugar

K = K(u) -o - -o

Isto não será feito para evitar o armazenamento si-

multâneo de duas matrizes de rigidez.

43

A linearização definida pelas expressoes (3.64)

possibilita uma solução direta para ~l e ~l resolvendo um sis

tema de ordem 2n análogo a (3. 48).

Ocorre que para cada intervalo de tempo uma nova

triangularização deverã ser feita, jâ que K nao -o e mais constante

ao longo do tempo, o que torna o procedimento pouco atrativo.

O algoritmo iterativo que sera usado aqui [2] uti-

liza as equaçoes (3.64) com a matriz

triz tangente modificada:

K -o

substituída por uma ma-

onde

K* -o

e a matriz de rigidez tangente avaliada apos

iterações com as equaçoes em sua forma original (3.64).

(3.65)

algumas

Com o uso da matriz modificada K* -o

é possível, g~

ralmente, se calcular as forças elásticas ~E ao final de cada

intervalo de tempo, por acumulação do valor . ~ .

no 1n1c1.o do interva

lo, ~EO com o produto em vez de se recalcular

a partir dos deslocamentos finais ~l

O esquema iterativo e descrito a seguir.

Inicialmente o vetor das forças de inércia ~l e

. sua derivada !l:i ao final do intervalo de tempo s ao estimados

por uma fÕrmula 11 predictor"; no caso:

. . !l:1 = R

-o (3.66a)

. !l:1 R + T R

-o -o (3.66b)

44

Esses valores sao levados às expressoes (3.44) ob­

tendo-se uma primeira aproximaçao para ~l , ~l e '::1 M-l R - -1

Novos valores para ~l podem ser agora cal

culados por:

~l e '::1 (3.67a)

. . ~l !1 - ~º '::1 - e '::1 (3.67b)

e levados novamente a (3.44) para a obtenção de uma nova aproxim~

çao para '::i , '::i e '::1

O processo se repete ati ser satisfeito o critirio

de convergência adotado.

Ressalte-se que, apos algumas iterações, K -o

nas

expressoes (3.67) i substituída por K* , j ã definida anteriormen -o

te. As condições ao final do intervalo de tempo são tomadas como

condições iniciais para o próximo intervalo e o processo se repe­

te ati ser completada a resposta no tempo,

A convergencia do processo iterativo descrito sera

estudada mais adiante.

O processo iterativo pode ser aplicado tambem a

problemas lineares desde que se use sempre as expressoes (3.46)no

lugar de (3 .67).

Observando-se as operaçoes:

M(~ - u ) - -1 -o

M '::1

T = TI (6 ~o+ T

. R -o

(3.68a)

(3.68b)

45

M(u - u - T ~) - -1 -o -o

(3.68c)

que sao efetuadas a cada iteração, nota-se q~e, para matrizes de

massa constantes no tempo, apenas uma triangularização da matriz

de massa é feita.

Observe-se também que as matrizes, ao contrário da

equaçao (3.48), conservam a ordem n e são simétricas o que e

vantajoso tanto no ponto de vista de memória interna quanto de

tempo de processamento.

Duas alternativas existem entao para problemas li-

neares:

- usar o algoritmo modificado incondicionalmente estável através

da equação (3.51),de ordem 2n ,com um intervalo de integração

função do problema especifico

usar o algoritmo condicionalmente estável, iterativo, com um in

tervalo de tempo definido pelo critério de convergência do pro­

cesso iterativo.

Tudo indica que a primeira alternativa seria pref~

rível no caso em que os modos mais baixos dominam a resposta, já

que intervalos de integração maiores poderiam ser usados.

Para problemas de transientes curtos, como choques

e respostas ricas em frequências, que requerem intervalos de inte

gração da ordem de grandeza do período mínimo do modelo discreto,

a segunda alternativa é a indicada.

O algoritmo correspondente a primeira alternativa

ainda não foi aplicado a grandes sistemas estruturais embora pro­

porcione_, nos exemplos apresentados em [3], resultados mais preci:_

sos do que os obtidos com o método de Wilson.

46

Para analisar a convergencia do processo iterativo

substitui-se as expressões (3.45) e (3.47) em (3.67) obtendo-se:

i+l ~1

~1

onde:

B = --o

=

A -o

B -o

~EO

+

A = -o

- K -o

- e

-

G

T 2 1 -K M-12 -o

[~o T

K + TI -o

T2 -1 u + 60 M -o -

M -1

-

(21

r;, + L:. º

T M-l (6 R TI - -o

(6 R -o

R + -o

que representa o processo iterativo:

-i+l -i R = P + A R

onde A e a matriz de iteração.

i )] . + T + ~1 -o

3 i ~ T --o

A condição necessária e suficiente para se

• i ~1

~1

(3.69)

(3.70a)

(3.70b)

(3. 71)

obter

convergência nesse processo e que os módulos dos autovalores da

matriz A sejam limitados pela unidade ~6]

Procurando alcançar resultados mais Úteis do que

aqueles apresentados em [2], a convergência serâ estudada para o

caso de um sistema massa-mola com amortecimento.

47

Nesse caso a matriz A se escreve:

6 T - (a 12 + b)

A

9 T2 6 ( + b _T) - ª 60 12

onde, por simplicidade, se tomou:

a =

os autovalores de A são complexos conjugados:

Desde que se verifique:

o que sempre ocorre nos casos de amortecimento positivo

o modulo dos autovalores e então:.

11- 1 = P

-A expressao acima pode ainda se escrever:

240 p 2

(3. 72)

(3.73a)

(3. 73b)

(3.74)

(3.75)

(b > O).

(3.76a)

(3. 76b)

48

relacionando o mÔdulo dos autovalores da matriz de iteração com

as características do sistema, eq. (3. 73), e o intervalo de tempo

adotado T •

Nesse ponto e interessante introduzir as relações:

a =

T = 8 T

(21T)2 T

(3.77a)

(3. 77b)

e tomar o amortecimento como uma fração do amortecimento crítico:

e ~ e . crit. ~ 2 M w

O parâmetro b passa a se escrever:

b 2 ~ w 41T ~ T

(3.77c)

(3.77d)

Substituindo-se o conjunto de relações (3.77) em (3.76b) ·obtém-se:

2 1T282 ~2 + 5 1T8 ç + ----zo- 3 p2

4 1T282 o (3.78)

que relaciona a taxa de amortecimento Ç

mÕdulo dos autovalores de A

e o parâmetro e com o

As raízes da equação do 29 grau em ç

1T8 5

rr 2 8 2 - ---

25

(3. 7 8) s ao:

(3. 79)

A condição para que o modulo dos autovalores de A

seja limitado por p acarreta:

49

(3. 80)

o que ocorre para os valores de ~ no intervalo entre as

da equação (3.78).

Introduzindo-se agora o conjunto de relaçÕes(3.77)

na desigualdade (3. 75) obtém-se:

(3.81)

-A expressao acima se anula em:

= rre (- ~ (3.82)

Portanto a região de autovalores complexos e aquela onde a desi -

gualdade (3.81) se verifica:

da forma:

~ < ~ 1 = - 1.61939 e

~ > ~ 2 - o.89388 e

Fora dessa região os autovalores de A -sao

(3.83a)

(3.83b)

reais

(3.84)

A condição para que tais autovalores tenham módulo

igual a p e:

50

1200(2rr) 2 6 2 1; 2 +. ~600(2rr) p 8 + 240(2rr) 3 8~ /; +

+ ~600 p 2 + 240 (2rr) 2 P 8 2 + 15 (2rr) "8~ o

(3.85)

Essa equaçao nao foi resolvida explicitamente, po-

rem, para cada par de valores p e 8 na região de interesse,

e possível calcular raízes /; 1 e

Os resultados obtidos ati aqui podem ser visualiza

dos fazendo-se sua representaçao no semi-plano i; - 8 (8 > O) , na

figura 3. 2 onde apenas o caso limite de convergência, p = 1 , está

representado.

A partir das semi-retas -,. oc e

-,. OB cujas equaçoes

sao, respectivamente:

- 1. 61939 8 (3.86a)

l;=-0.893888 (3.86b)

Definem-se as regioes:

- autovalores reais: região interior as semi-retas -,. oc e

-,. OB

- autovalores complexos: região externa dada pelas desigualdades

(3. 83)

Na região de autovalores complexos as curvas Al e

A2 sao a representação das ~

ra1zes e respectivamente,

dadas pela expressao (3.79), no caso limite p = 1

Essas curvas tendem assintoticamente para o eixo

8 = o Cada ramo da bifurcação que se observa entre B e. e cor

responde à condição de um dos dois autovalores reais de A apr~

51

sentar módulo unitârio.

Dessa maneira a regiao de convergencia e aquela in

terior as curvas Al e A2 e ao ramo inferior da bifurcação en-

tre B e e .

Para valores de p < 1 obtém-se curvas semelhan-

tes e interiores a Al e A2 cada vez mais próximas da origem,

com a propriedade de se encontrarem sempre sobre a semi-reta

cuja equação:

_,_ OA

Tr s e (3.87)

e obtida a partir de (3.79)

Para o caso de amortecimento positivo a figura 3.3

mostra as curvas Al para alguns valores de p

No problema nao amortecido a situação se simplifi-

ca. Fazendo-se i; = o em (3.78) obtêm-se:

p 92 (3.88)

que dâ origem ao grâfico da figura 3.4.

O valor limite de e para que haja convergencia

corresponde ao caso p 1 e, de (3.88) obtêm-se:

e . = o.626 cr1t. ( 3. 89)

e o intervalo de integração deverá ser:

T < 0.626 T (3.90)

52

Observe-se que a convergencia em um intervalo de

tempo e mais rápida a medida que se diminui T (fig. 3.2, 3.3)

porem isto acarreta um numero maior de intervalos de

para se descrever o mesmo intervalo total.

integração

Para cada problema especifico deverá existir entao

um intervalo de integração Ótimo, com relação ao tempo de execu-

çao, que pode ser determinado por tentativas a partir de (3.90).

Normalmente, este intervalo de tempo atende também

a condição de estabilidade, eq. (3.55).

As conclus~es relativas ~ converg~ncia,obtidas aqui,

podem ser extendidas ao caso não-linear levando-se em conta que a

matriz de iteração ~ , devido à variação da rigidez da mola, se-

ra apenas seccionalmente constante ao longo da integração.

e= .r.. T

1.0

0.9

o.~

0.4

03

0.2

0.1

- 2.2 -1.8 -L6 -1.4 - 1.2 - 2. O -1.0

- 0.8 -0.8 -0.4 -0.2 o.o 0.2

FIG. 3.2

0.4 0.8 o.e 1.0

1.2 1.4 1.8 LI 2.0

2.2 2.4 2.8

v> w

J' '1.0

0.3

0.2

0.1

o.o'----------------------------------------------+------0.2 0.4 o.e o.a 1.0 1.2 1.4 1.8 ,.

e ~ ·-e crtt.

FIG. 3.3

'·" .,_,

1.0

0.9

o.a

0.7

0.6

O.OI

0.4

0.3

0.2

0.1

o.o

f-lÀI

8: Í '---===::....,..----,---~-----.----..----~-- T

0.1 0.2 0.3 0.4 o.a 0.8

FIG. 5.4

56

IV - PROCEDIMENTOS IMPLEMENTADOS~ EXEMPLOS - CONCLUSÕES

4.1 - INTRODUÇÃO

Neste capítulo ê feita a aplicação do algoritmo

condicionalmente estável, obtido pela interpolação ciibica das for

ças de inércia, visto no capítulo anterior, ao caso da resposta

dinâmica de estruturas reticuladas submetidas a solicitações de

natureza determinista.

Adotou-se o processo iterativo que usa a matriz de

rigidez tangente modificada

nearidade geométrica.

K* -o e é considerado o caso da não-li

As equações de movimento sao ,obtidas a partir do

princípio de Hamilton com a discretização espacial feita pelo me­

todo dos elementos finitos.

4.2 - ELEMENTOS IMPLANTADOS

Foram implantados cinco tipos de elementos finitos

derivados de dois tipos básicos: o pórtico plano e a treliça esp~

cial. são identificados, no programa, pela variável ITIPO

!TIPO = 1 - Pórtico Plano

ITIPO 2 - Treliça Plana

ITIPO 3 - Cabo Plano

ITIPO = 4 - Treliça Espacial

!TIPO = 5 - Cabo Espacial

Ós elementos de cabo sao usados apenas em análise

não-linear e diferem dos elementos de treliça somente no fato de

nao oferecerem resistência a compressao. Suas matrizes de rigi-

dez sao relacionadas por:

57

Kdàbd = H(E) Ktteliça X

onde H(E ) X

e a funçio degrau de Heaviside:

e E X

(t-t)/t o o

tante considerado.

H(E ) X

H(E ) X

o

1

E < Ü X

E > 0 X

e a deformação unitária do elemento no

( 4. l)

(4.2)

ins

Os elementos podem ser combinados em dois grupos

a) elementos dos tipos 1, 2 e 3

b) elementos dos tipos 4 e 5

permitindo a anilise de p5rticos planos enrigecidos com cabos e/6u

elementos de treliça, estruturas planas formadas por elementos de

cabo e/ou treliça (a) e estruturas espaciais formadas por elemen­

tos de cabo e/ou treliça (b)

As matrizes de massa sao idênticas para os elemen-

tos de cabo e treliça.

de ordem

6 X 6

Todos os elementos apresentam matrizes

correspondentes a 3 grius de liberdade por no. Os elemen-

tos planos de treli·ça e cabo tem as suas matrizes com as linhas e

colunas correspondentes ao terceiro grau de liberdade nulas para

poderem ser combinados diretamente com elementos de pÕrtico plano.

No que se segue, apenas serao mostradas as expres-

soes correspondentes aos elementos básicos: pÔrtico plano e tre-

liça espacial.

58

4.2.1 - Matrizes de Massa

As matrizes de massa dos elementos surgem, quando

da aplicação do princípio de Hamilton, a partir da expressão da

energia cinética:

T ( 4. 3)

onde:

µ - massa específica do material

v = velocidade do elemento dV

V= volume do elemento finito

a) Elemento de PÕrtico Plano

Considera-se o elemento de pÕrtico plano de eixo

reto de comprimento L e seção transversal constante de area A,

simétrica em relação ao plano de carregamento.

Supondo-se a conservaçao das seçoes planas durante

a flexão e desprezando-se a deformação por cisalhamento tem-se p~

ra o deslocamento longitudinal dos pontos de ordenada

ção ao centro de gravidade da seção:

u = u* - z av ãx

z

onde u* = ui z=O

transversal e v

e o deslocamento longitudinal do C.G. da

e o deslocamento transversal da seção.

A energia cinética (4.3) se escreve:

dAdx

em rela

( 4. 4)

seçao

( 4. 5)

59

Levando (4. 4) a (4. 5) e integrando-se sobre a are a

da seção transversal obtim-se:

T 1 2 µ A f +

L

1 2 µ A f

L

• 2 u* 1

dx + 2

onde J e o momento de inêrcia da seção, e a energia cinética to

tal estã separada em parcelas que correspondem, respectivamente,

as translações vertical e horizontal e à rotação da seção trans -

ver sal.

Introduzindo as funções de interpolação:

u* = N q -u -

V ~V 9

na expressão da energia cinitica (4.6) obtim-se:

T • T s

onde N - ;V' X

d N dx -v

N -v + N T

-u

Desenvolvendo-se a parcela:

do princípio de Hamilton tem-se:

N + N T -v -u

que, integrada por partes, resulta em:

g dx

s óg dt

( 4. 7 a)

(4.7b)

( 4. 8)

( 4. 9)

(4.10)

ou ainda:

N -v

+ N T -u

60

N l + J Íi, ! N J } dx g -~ ~v,x -~,x

Óq dt (4.lla)

(4.llb)

onde Me e a matriz de massa, ji que, pelo princípio de Hamilton:

= o (4.12)

Empregando-se funç~és de interpolaçio de Hermite,

lineares para o deslocamento u* :

o o X

L o

e cúbicas para o deslocamento transversal v

~V = [o ( X 2

X 1 - -) L

o

(4.13a)

-(- - 1) x2

X J L L

(4.13b)

e considerando os Uesiocamentos nodais g segundo a fig . (4.1):

. .'J .. ---i)>---------q-& __ J_._q-5

-,,~------s,:, X, u

ql 2 q4

FIG. 4.1

61

obtém-se a matriz de massa consistente para o elemento de pórtico

plano:

11 o o 1 o o l rc o o o o o l 6

1

13 . 11 9 13 6 .l_ 6 1 1 35 210 L o 70 420 L 5L 10

o. - 5L 10 1 1 2 13 1 2

1 2L 1 L 105 L o 420 L - 140 L TI o 10 - 30

Me \J AL 1 1 + )J J

3 o o o o o S I M • s I M •

13 11 1 6 1 35 - 210 L 1

l 5L 10

1

1 j 2L J 105 L 15

(4.14)

onde a primeira parcela corresponde a inércia de translação

segunda ã inércia de rotação.

e a

A matriz de massa é dita consistente quando sao

usadas na avaliação da energia cinética as mesmas funções de in­

terpolação que são utilizadas no cálculo da energia de deformação

no principio de Hamilton.

O conceito de matriz de massa consistente se deve

a Archer [29] e antes desse seu trabalho apenas matrizes de massa

diagonais, obtidas concentrando-se massas nos pontos nodais, eram

empregadas.

A formulação consistente mantém a propriedade do

método de Ritz de fornecer as frequências do modelo discreto sem-

premais altas do que os valores exatos.

Entre as matrizes de massa nao consistentes, sao

de particular interesse aquelas que se apresentam em forma diago­

nal jã que tornam explicito o algoritmo de integração aqui empre-

62

gado. Ressalte-se tambêm que as matrizes co-õ.sistentes sempre

acarretam maior esforço computacional mas nem sempre levam a me­

lhores resultados [9, 22, 24]

O procedimento descrito em [22] e aplicado t anib êm

em [24] para a obtenção de uma matriz de massa diagonal para ele-

mentas isoparamétricos quadráticos resume-se em distribuir a mas­

sa total do elemento proporcionalmente aos coeficientes da diago­

nal principal da matriz consistente do elemento.

A matriz de massa diagonal para o elemento de por-

tico plano empregada nesse trabalho ê obtida a partir da matriz

consistente segundo o procedimento, descrito em [21], que se es-

creve:

D + mll mll ml4

D + m25 m22 = m22

D + m36 m33 m33

( 4. 15) D

m44 m41 + m44

D + mss = m52 mss

D + m66 m36 m66

A matriz resultante e:

Me G JJ AL 1

AL 1 3 1 1 1 JJ AL~ - jJ 420 JJ AL - JJ AL - JJ AL 420 + 2 2 2

+ [ o o 1 10 JJ JL o o 110 ~ JL J (4.16),

onde a segunda parcela provem da inércia de rotação.

63

No programa em apindice duas variiveis, MASSA e

INROT, controlam a c,nstruçao da matriz de massa da estrutura de

acordo com os valores que assumem:

MASSA= 1

MASSA= 2

MASSA= 3

INROT o

INROT 1 O

e montada a matriz de massa consistente

e montada a matriz de massa diagonal

no caso da matriz de massa global ser diagonal po-

rem construida a partir de elementos escalares li-

dos pelo programa. Nesse caso a massa específica

dos elementos MESP(I) deve ser nula.

nao e levada em conta a parcela correspondente a

in~rcia de rotaçio na montagem das matrizes consis­

tentes ou diagonais dos elementos.

a parcela de inercia de rotaçao e considerada.

Observe-se que, em todos os casos, sao ainda leva­

das em consideração as massas discretas introduzidas como elemen-

tos escalares pela matriz ED(I, J) conforme o item 4.3.

b) Elemento_de_Trelisa_EsEacial

A matriz de massa para o elemento de treliça espa­

cial pode ser obtida utilizando-se a expressão da energia cineti-

ca:

onde U , V e w

direções usuais

T 1 2 ~µ(~2

V

dV (4.17)

sao as componentes do deslocamento segundo as

X ' y e z e as funções de interpolação:

64

u N g -u

V = N g (4.18) -v

w = N g -w

Levando-se (4.18) a (4.17) e por procedimentos análogos aos já

vistos no caso do elemento de p6rtico plano obtim-se a matriz de

massa do elemento de treliça espacial:

~ 11 J ~u T Nu + Nv T Nv + Nw T NJ dV

V

(4.19)

Empregando-se as funções de interpolação de Hermi-

te lineares para .os deslocamentos u, v e w e adotando-se os

deslocamentos nodais segundo a fig. 4.2:

FIG. 4.2

obtêm-se a matriz de massa consistente do elemento:

1 o .o 1 o o 3 6

1 o o 1 o 3 6

1 o o 1 3 6

Me µ AL 1 (4.20)

3 o o

s I M 1 3 o

1 _\ 3

65

A construçao de uma matriz diagonal i Õbvia, nesse

caso, bastando alocar metade da massa total do elemento a cada um

dos nÕs. A mesma matriz pode ser obtida substituindo-se as fun-

çÕes de interpolação lineares (4.13a) por funções

___ v_i_s_i_d_e ,__f_~g .• _4. 3_. __

degrau de Hea-

+- L/2 -t- L/2--+-

FIG. 4.3

N1 = H(x) - H(x - ~)

H(x - !'.) 2

(4.20a)

(4.20b)

Ainda a mesma matriz e obtida adotando-se o procedimento usado an

teriormente no caso do pÕrtico plano.

A matriz resultante e:

4.2.2 - Matrizes de Rigidez

1 2

1 2

1 2

1 2

A matriz de rigidez do elemento provem,

aplicação do princípio de Hamilton, da parcela

deformação associada ao deslocamento virtual ou

ou

(4.21)

quando da

energia de

Serão deduzidas a matriz de rigidez para o elemen­

to de pÕrtico plano e, a partir dela, a matriz para o elemento de

treliça espacial.

A energia de deformação e definida por:

u = i f V

T E: o dV (4.22)

66

onde:

E e o vetor de deformações específicas

cr e o vetor de tensoes associado

Para materiais elásticos lineares a relação entre tensoes e defor

mações ê dada pela lei de Hooke generalizada:

O = E E (4.23)

onde E e uma matriz de constantes elásticas e a expressão(4.22)

se escreve:

u 1 2 J ET E E dV

V

(4.24)

De acordo com as hipóteses assumidas no item 2.2.1

a respeito do elemento de pórtico, o vetor de deformações especí­

ficas se reduz a E e (4.23) e (4.24) se escrevem: X

tensoes:

au ax e mantendo

u 1 E 2

E E X

2 E dV

X

(4.25)

(4.26)

A partir da expressao da teoria do estado plano de

E XX

Desprezando-se o termo .!. (~) 2 2 ax em

(4.27)

presença de

já que e a contribuição de ordem mais

67

baixa de v(x) incluindo deformação flexão ô 2 v , e por - z e dX 2

deformação inicial E tem-se [ 3 O] o

dU 1 c"v) 2 d 2 v (4.28) E = E + dX

+ 2 - z X o ax

dX 2

Substituindo-se (4.28) em (4.26) e assumindo fun-

çÕes de interpolação para os campos de deslocamento:

N N1 i (4.29a) u = :! = q -u u

N Ni i (4.29b) V = :! = q -v V

onde :! sao os deslocamentos nodais da fig. 3.1 e vale a conven­

ção da soma sobre os Índices repetidos, tem-se a energia de defoE

mação expressa em função dos deslocamentos nodais:

locando-se a sua variação na forma:

ou= a q. a q. q. o q.

J l J l

ou, em notaçao matricial:

ou

-

U = U(q.) .Co l

(4.30a)

(4.30b)

chega-se a expressao dos elementos da matriz de rigidez:

-

k .. = lJ a q. a q.

l J (4.31)

Desprezando-se na expressao da energia de deforma-

çao (4.26) os termos de grau menor que 2 nos deslocamentos nodais

;! por não contribuírem para a matriz de rigidez (4.31) obtém-se:

68

u 1 E

~ ru2 + 2 2

v2 ~ dV + 2

.Z V + E L ,x ,xx o ,x

1 E ~ [t v4

v2 ~ dV + + - + u 2 ,x ,x ,x

1 + 2 E f . rz.v L ,xx

(2 u ,x

+ v 2 . jl , X '.J dV (4.32)

V

onde as vfrgulas denotam derivação.

A Última integral e nula. Desprezando-se a segun-

da integral em presença da primeira, introduzindo-se as funções

de interpolação (4.29) e usando a expressão (4.31) obtem-se:

k .. ]. J = EA J

L

+ EA E o

i N• N j u,x u,x dx + EJ J

L

i N v,x

N J dx v,x

i N v,xx

dx +

Usando-se as funções de interpolação (4.13),

se, em forma matricial explícita:

1 o o - 1 o o

o o o o o

EA o o o o

Ke = + L 1 o o

s I M o o

o

(4.33)

tem-

69

o o o o o o o o o o o o

l 12 6 o · 12 . 6 6 .. 1 6 1 L' i2 L' L2 5L 10

o 5L 10'.

4 o ·5 2 . 2. 1 L L - i2 L 15 L o 10 + EJ + EA E

30 o

o o o o o o 12 6 6 1 1' -i! 51 10

4 2 . L 15 L

As duas primeiras parcelas correspondem ã matriz de rigi­

dez linear do elemento de pórtico plano e a Última parcela , e a

chamada matriz de rigidez geométrica, onde:

N = EA (4.35)

i o esforço normal calculado no inrcio do intervalo de integraçio.

b) Elemento_de_Treli~a_Eseacial

Sua matriz de rigidez pode ser obtida pelo procedi

a2 v menta visto anteriormente retirando-se a parcela - z -- cor-

ax2 respondente a flexão, na expressão (4.28) e introduzindo o termo

Conservando-se apenas os termos quadráticos nos

deslocamentos nodais g na expressão da energia de

tem-se:

U = _l E 1 tu 2 + E 2 ,x o

V

(v2 + w2 ;l ,x ,x J dV

"deformação

(4.36)

Introduzindo-se as funções de interpolação linea-

res e aplicando a expressão (4.31) obtém-se, em forma matricial:

explícita:

1 o

o

o

o

o

- 1

o

o

1

4.3 - ELEMENTOS ESCALARES

o

o

o

o

o

o

o

o

o

o

oj

70

EA E: · . o +---

L

ro o

1

o

o

1

o

o

o

o

o

- 1

o

o

1

o

o

- 1

o

o

1

(4.37)

são chamados de elementos escalares, ou discretos,

as molas,amortecedores e massas discretas que podem ser

das à estrutura nos seus pontos nodais.

associa-

No programa implementado, atuam sempre segundo di­

reçoes globais estando associadas a eles forças elâsticas, dissi­

pativas ou de inércia, atuantes nessas mesmas direções.

As molas sao supostas lineares de constante k e 1

os amortecedores do tipo viscoso linear de constante ci•

As constantes dos elementos escalares (m , c , k)

sao somadas diretamente ao termo da diagonal principal da matriz

global correspondente (M , C , !) , na linha associada ã direção

em que o elemento atua.

A introdução dos elementos escalares ê feita atra-

ves do vetor NO(I) que guarda o número dos nos com elementos

escalares, e da matriz ED (I , J) cuja .linha I armazena as cons

tantes dos elementos do no de numero NO(I) na seguinte ordem:

71

J = 1, 2, 3 - massas segundo as direç;es globais 1, 2 e 3

J 4, 5, 6 - constantes dos amortecedores segundo as

ç;es globais 1, 2 e 3

J = 7, 8, 9 - constantes das molas segundo as direç;es

b ais · 1, 2 e 3

dire-

glo-

I varia de 1 atê NNED numero de nos com algum elemento discre

to.

4.4 - MATRIZ DE AMORTECIMENTO

Conforme visto no item 3.2.1 o modelo de amorteci-

menta mais usado é o amortecimento viscoso linear que, introduzi-

do diretamente no modelo discreto, dâ origem ao termo

equação de movimento (3.13)

c u na

matriz c

Para que o amortecimento seja dito proporcional, a

deverá ser expressa pela serie de Caughey [5]

c M l: a. ~-1 0 i i - :L l_:: :J (4. 38)

Em [3~] dois processos sao apresentados para o cálculo de c . o

primeiro determina os coeficientes a. l

de maneira a introduzir

uma fração estipulada do amortecimento crítico em cada modo. O

segundo, mais direto, expressa a matriz c como uma soma de ma-

trizes, cada uma delas introduzindo o amortecimento desejado em

cada modo.

Como a matriz obtida nao e em banda, e sugerida em

[32] a zeragem dos termos fora de uma certa banda.

Uma forma muito usada para a matriz e e a corres

pondente ao amorteciemnto de Rayleigh e pode ser obtida a partir

dos dois primeiros termos (i o i 1) da serie de ~Cauihey:

72

e= a M + a1

K o -

(4.39)

Quando o amortecimento varia sensivelmente ao lon-

go da estrutura é possível [12] se estabelecer diferentes valo-

re s de

frações

a o

e para diferentes .r·egiÕes da estrutura adotando

do amortecimento crítico convenientes. O amortecimen

to obtido nesse caso não ê proporcional.

Muitas vezes ê vantajoso, do ponto de vista compu­

tacional, o uso de uma matriz de amortecimento diagonal e em [33]

sao apresentados procedimentos para a sua obtenção.

Em todos os casos os valores numéricos para os fa-

tores de amortecimento em cada modo ( i; . ) são estimados ou basea 1

dos em dados experimentais da estrutura ou de o"utras estruturas de

geometria e materiais semelhantes.

O programa incluido em apêndice apresenta as se-

guintes alternativas de acordo com o valor assumido pela variável

de controle IAMOR:

IAMOR = 1 - Não se considera amortecimento, nao

as operações correspondentes.

sendo

IAMOR = 2 - O programa monta a matriz de amortecimento

C =ALFA* M +BETA* K onde ALFA e BETA

efetuadas

na f arma

sao pa-

râmetros lidos e as matrizes M e K são globais e

calculadas em t = o são tambêm considerados os amortecedores discretos in

traduzidos pela matriz ED

· IAMOR = 3 - A matriz de amortecimento, considerada simêtrica e de

mesma largura de banda de M e ; , e lida como um

vetor C(I) que contêm apenas os elementos incluídos

em sua banda superior.

73

Uma outra alternativa que pode ser implantada sem

maiores dificuldades e a formação de matrizes de amortecimento a

nível de elemento da forma onde ALFA

e BETA são fornecidos para cada elemento. A matriz global seria

montada por acumulação como no caso de M e K

No caso de se tomar, para um sistema de um grau de

liberdade, o amortecimento na forma:

c = a m + S k (4.40)

e exprimindo-o como uma fração do amortecimento crítico:

c = i; c • cr1t.

= 2 I; m w (4.41)

conclui-se, de (4.40) e (4.41), que:

"Í (~ + S w) (4.42)

Os resultados obtidos a partir do estudo da conver

gência do processo iterativo no caso linear e para um grau de li­

berdade podem ser extendidos ao caso de vários graus de liberdade

desde que o amortecimento seja do tipo proporcional. Nesse caso

o sistema se desacopla e, no modo i o amortecimento introduzi-

do e dado por (4.42) com w = w. 1

Adotando-se um intervalo de integração T e poss.!.

vel, com o uso da fig. 3.2, verificar se os pontos

vos da integração de cada modo, p. ( i;. 1 1

1 (~ 2 w.

1

+ B w.) 1

8.) 1

com:

representa ti-

(4.43a)

e. 1

7 !,

·, T

T. = 2rr wi 1

estao na regiao de convergencia.

(4.43b)

Normalmente o modo crítico ê o modo correspondente

a frequência de corte bastando verificar o ponto 8 ) n

i;n 1 (-ª- + B w ) 2 w max

(4.44a) max

8 T

2rr w n max (4-44b)

Muitas vezes e necessário adotar um intervalo de

integração menor do que o correspondente ao caso nio-amortecido o

que constitui uma desvantagem inerente ao processo iterativo de

solução e nao ao algoritmo de integração em si.

4.5 - VETOR DE EXCITAÇÕES EXTERNAS

O vetor F de excitações externas do modelo dis-

ereto (3.13) provêm, quando da aplicação do princípio de Hamilm~

da parcela trab3lho virtual das forças não conservativas.

sua forma geral, a nível de elemento, e:

J vol.

dV + J sup.

(4.45)

onde b e t sao forças de volume e superfície respectivamente.

O vetor de forças nodais equivalentes e dito

consistente se são usadas em (4.45) as mesmas funções de interpo­

lação usadas na expressão da energia de deformação.

Assim, uma formulação e dita consistente quando

75

sao usadas as mesmas funções de interpolaçio na discretizaçio de

todas as parcelas do princ{pio de Hamilton,

O vetor global F e montado somand~-se as contri-

buiçÕes de cada elemento.

Conforme jã visto no Cap{tulo III o algoritmo em-

pregado faz uso dos vetores F e • d F = - F - d t -

Ambos os vetores deverão ser montados diretamente

no referencial global através da subrotina EXCITA

são mostradas a seguir as expressoes de

para um elemento de pÕrtico plano horizontal percorrido

e • e F

por uma

carga mõvel concentrada vertical, de valor P , movendo-se da es­

querda para a direita.

ças se escreve:

onde s ( t) e a

figura ( 4. 4) .

(4.46) se reduz

Fazendo uso da funçio o de Dirac o vetor de for-

L

f N T P o -V G - s (t~ dx

o

abcissa da carga medida a

y

l p

• ( t) ,t 1 2

'L

FIG. 4.4

Da definiçio da funçio o

a:

Fe = p NT 1 -v

'

x=s(t)

partir do

) X

de Dirac

(4.46)

no esquerdo,

a int.egral

(4.47)

Derivando-se o

pe

ou ainda:

onde V ( t)

= p

ds dt

76

vetor Fe obtem-se:

r FS(,J

d ~VT

= dt

Fe = P v(t) d ds

~; SªS ,,J p d

ds

e a velocidade da carga móvel.

ds dt (4.48)

(4.49)

As expressoes (4.47) e (4.49) sao calculadas auto­

maticamente pela subrotina EXCITA listada em apêndice.

O trecho da estrutura percorrido pela carga móvel

deverá ser horizontal e ter seus elementos numerados sequencial-

mente, da esquerda para a direita, de 1 a NEP (número total de ele

mentos percorridos).

querda.

Em cada elemento, o nó 1 deverá ser o da es

4.6 - CRITÉRIO DE CONVERG~NCIA

Em todos os processos iterativos e necessário deci

dir se o resultado de uma iteração está suficientemente próximo

da solu-ção exata, desconhecida, ou s·e e necessário efetuar-se mais

uma iteração. Adota-se entao um critério de convergencia para

que essa decisão seja tomada.

Dados dois vetores e correspondentes

as i teraçÕes J e j + 1 respectivamente~ o critério de conver

gincia adotado se· escreve:

7 7

li vj+l - vj li

li vj li < y ( 4. 50)

onde y é a tolerância especificada e 11 • 11 indica a norma ve-

torial empregada, que no caso foi a norma euclideana:

li V IIE =

~

n li E V~ .

i=l i J (4.51)

Mais considerações sobre critérios de convergencia

na solução de problemas estruturais não-lineares

iterativos são encontradas na ref. [34]

processos

O critério adotado, equaçao (4.50), ê aplicado, em

análise não-linear, ao vetor incremento de deslocamentos llu , P.'.:

ra o cálculo da matriz K* -o

1 - (K + K_

1) (vide fluxograma).

2 -o

Em todas as análises, para se interromper o prece~

so iterativo, o mesmo critério e aplicado ao vetor R = M u , de-

rivada das forças de inércia, por estar ligado ã derivada mais al

ta, em relação ao tempo, do vetor de deslocamentos u

78

4.7 FLUXOGRAMA

1 Inicio 1 !

1 Leitura dos dados 1

1 DeflniqÔ o 41 conatonU•

! IIOt1ta1en1 dat matrlr.11 globais ! • !o

• do vetor !!,•o

l ; 1 1 ,3

1 IAMOR

•2 Leitura da Matrt1 g 1

11,n,a1•• de E zsol!+f!'!o

l .FALSE.

: NNED;i.t 1

.TRUE.

1 lotfodução do

oaortecedort• llscretos

-~ ,, o.

1•T= O

BOOL•.fALSE.

,, 1 CARIIEG 1

•2

•••clollzac;ão •• vortáwel•

r1lot1va1 o cario ' •• .,.1 ., .• ......... do f( ti o !ltl

! l llllrcNhlção ••• eoa•1i;:ÕH •• cNtono .. M 1

1 o

1

7 9

• !

1 llnpreuõo dOI Matr1.z .. !! ,, ~ e !§. 1 1

Cálculo da aceleração iatctol otrave, de

R = IIÂ.i =F •• -.. •• l

= 1 IAIIOR

I 1

1 CCÍlculo do ;Jo .1 .•

1 . . .!to='o·Siio

Q INT= INT-t-1 t : INT -:M TAU

BOOL=.FALSE.

ITERA • O

Fórmula PREDICTOR . , ..tt, =!o

1

11 = !o+ G ~ o

l Cálculo auxiliar

Ti-< 8]o+l; ~o>

2 • ~(21 R0+3G!l.,) 6 - -

l = 1 1

CARREGI '

• 2

Atu olizaçõo •• parâmetro,

retat1vo1 a corto móvel

' • Monto9e1n •• f( t) • flt)

b

80

•

ITERA , ITERA+ 1

!'ITERA

_ .. 1 .1 1 Obttni;ao de y .4!! • A! otravea de:

Míii = R1 --· ... , . i 1i , 1 • 1 llt.u , _ < e•-,1õ•-H 1-1iiR1 > "" - 12 .!.'O ~"O- "' ~

i • ..,_2 • 1 ·1 11(6u-1õu ),-'!.. (21 R+51õR•9R-21óR) "' "' -. 80 -o -o -, ....

,, ~----'-~IAIIOR

eâi1cu10 dt Cii1 • ci1 ...... , ...... ,

.TRUE, IOOL,OR,(TIP,.N,EQ,I) 1------ol

.f'ALSE,

.FALSE. ITERA<t

.TRUE,

IJJ191 JlyHIJ ---· T_R_U_E_ • ..., . 1 <;roL

!AJ( j

.FALIE,

ITERA. GE, ITIIAX ,f'ALSE.

.TRUE,

6

,i , l•I !, = I! 1

Ri_ R i+I ... 1 - ... 1

a_J-'= a .... i

2

81

A111at1acu;ão da Geo1n.tria

• 1 K 0 , K0=-( K0+K1) - - 2 - -

IOOL•.TRU!.

< TOL21--·f_A_LS_!_. ___ ~

ITERA.GE.ITIIAX .FALSE.

.TRUE.

lmprHt ÕO doa resulta dos H

INT é MÚitipio .. IIIPRES

lmpr111ão doa rHultados;

proce11a•ento interroapido

.TRUE. t ,> TF 1--------------il F 111

.FALSE.

R , R~I -o ... ,

!0=6:·1

JJ.=s: !o=s!

TIPAN

Re:Re+K6• 1

"' o .., o ........

'1~-"'---, 11-------i TIPAN

,z

Atualização da Geoaetrla

Monto1•• clt !o

,ao confl1uracão atual

,5 TI PA N 1-------,

Monta1e11 de !o• !•o

na conHguroção atual

lapre11ão de Hforooa ••

INT ~ m~ltiplo de UIPRES

,_ __ .,

82

4.8 - EXEMPLOS DE APLICAÇÃO

são apresentados a seguir alguns resultados obti­

dos utilizando-se o algorrtmo ji descrito e implementado no pro­

grama listado em Apêndice.

Foi desenvolvido também um programa auxiliar atra-

ves do qual as caracteristicas din~micas dos modelos discretos

(frequências e modos naturais) foram determinadas.

EXEMPLO I

Como primeiro exemplo tomou-se o problema de vibr~

çoes longitudinais de uma haste homogênea, de seçao transversal

constante, mostrada na figura 4.5 onde os dados estao num sistema

de unidades coerentes. Uma das extremidades da haste e fixa e a

outra foi submetida a uma carga cuja variação no tempo e dada tam

bem na figura 4.5.

8 E=0.2xl0

,.. o. 8

AI: 0. 1 F( t)

:!-------------- -

+------,~oo------+

FIG. 4.5

;;, -4

1 10 s

;, t

A discretização espacial foi feita por meio de 40

elementos de treliça .plana de mesmo comprimento.

O modelo discreto obtido com o uso da matriz de

massa consistente apresenta o petrodo fundamental

= 0.4 10- 2 s e o período 4 •

m1.n1.mo, correspondente a

T max T•

f =

frequência de

83

corte, T . m1.n

. -4 = 0,4537 10 s Os valores exatos (modelo contínuo)

[s] sao, respectivamente., -2

Tf = .0.4 x 10. s e -4 r

40 = 0.5063 10 s.

Foram analisados dois intervalos de· tempo T = 0.1 -4

10 s,

correspondente a 0,0025 Tf e 0.220 T . m1.n e

-4 T~= 0.2 10 s,corres

pondente a 0.005 Tf e 0.441 T . m1.n

Os diagramas do esforço normal ao longo da haste, para

alguns valores de t encontram-se nas figuras 6.. 7 e 4. 8 onde coe~

forço normal de compressio foi tomado como positivo. A soluçio exa

ta do problema contínuo ê tambêm apresentada na mesma figura.

Ambos os intervalos de integração levam a resulta­

dos satisfatórios nao se observando uma melhora significativa ao

se diminuir o intervalo para -4

T = 0.1 10 S

Adotando-se a matriz de massa diagonal, o modelo

discreto obtido apresenta o período fundamental

e o período mínimo T . m1.n

-4 = 0.7855 10 s

Tmax = Tf

Foram analisados os mesmos intervalos de

T = 0.1 10- 4 s = 0.127 T . m1.n

e T -4

0.2 10 s 0.255 T . m1.n

= 0.4

tempo:

e al-

guns resultados encontram-se nas figuras 4.9 e 4.10. Tambêm nes-

se caso os resultados para os dois intervalos de integração sao

satisfatórios, confirmando as conclusões da referência [35] segu~

do as quais a matriz de massa consistente não apresenta, nesse ca

so, melhor convergência, para as frequências e modos naturais, do

que a matriz de massa diagonal adotada.

Observe-se que o uso da matriz de massa diagonal,

além de tornar o algoritmo explícito, dá origem a um modelo dis-

ereto cuja frequência de corte ê inferior à do modelo consistente

o que torna mais rápida a convergência do processo iterativo (fi-

gura3.2). Isto acarretou, para -4

T = O. 2 10 s , um tempo de processamento·

84

da ordem de 40% do te~po do processamento gasto ao se utilizar a matriz consis

tente.

Em todos os casos adotou-se TOL2 = 0.1 10- 3

Em seguida considerou-se o mesmo problema para o

caso em que a haste ê composta de dois materiais com mÕdulos de

elasticidade diferentes como indicado na figura 4.6.

8 E

1 :0.8 X 10

c1

=10 000

FIG. 4.6

8 E=0.2xl0 2

e = 5 ao o 2

As características dinâmicas dos modelos sao agora:

massa consistente T = 0.2552 l0- 2 s T . min

massa diagonal

max

T = 0.2552 l0- 2 s max T • min

-4 0.3938 10 s

Para o caso de matriz de massa diagonal foram ana-

lisados os intervalos de integração -4 T = 0.1 10 s = 0.254 T . e min

T = 0.2 10- 4 s = 0.508 T . e alguns resultados estao, min

com a solução exata, nas figuras 4.12 e.4.13.

Já ao se fazer uso da matriz consistente

juntamente

apenas o

intervalo T = 0.1 10- 4 s = 0.437 T . min pode ser analisado (figura

4.11) uma vez que o intervalo T = 0.2 l0- 4 s = 0.875 T . está fo min

ra da região de convergência do processo iterativo (figura 3.2)

85

DIAGRAMAS DO ESFORÇO NORMAL AO LONGO DA HASTE

2.0

1.,

1.0

1.0

o.

FIG.4.7

Mo11a Conslltente ;,:0.00001 1

-3 t= 0.4 a 10 1

-3 t=o.e x 10 ,

-3 t: 1. 2 1 10 1

-3 t=2.1110 1

86

DIAGRAMAS DO ESFORÇO NORMAL AO LONGO DA HASTE

2.0

L~

1.0

1.0

o.

o.o,r1,~-=-=.-r:=----'=,ç--------=.,.....---

-0.

-1.0

FIG. 4.1

Con1l1tent, T,:Q.00002 •

-3 hO.I a 10 1

t:1.2 x 10-3

,

t::1.e • 10 -3,

87

DIAGRAMAS DO ESFORÇO NORMAL AO LONGO DA HASTE

Massa Dtscr.ta W = 0.00001 1

te 0.4. ,o-3,

t:O.Bal0-31

f:2.1 X 10 -:SI

LO

o.

FIG. 4.9

88

DIAGRAMAS DO ESFORÇO NORMAL AO LONGO DA HASTE

LO

O.&

FIG. 4.10

Ma•a Dttcreta "t,:Q,00002 s

-3 h-0.4 lll 'º 1

-3 t=o.ax10 1

-3 t = 1.2 x 10 1

f:2.1 X 10 -3

1

1: 2.6 a 10·3,

89

DIAGRAMAS 00 ESFORÇO NORMAL AO LONGO DA HASTE 1.

Mana Con1t1tente 'G: 0.00001 ,

o.o -3

f=0.4,: 10 1

o.ol-------------=lf-..-=-------~

1.8

1.0

o.o

o.o

1.8

1.0

o.e

o.o

2.0

2.

2:.oJ==~

LI!

1.0

o.

EI E2

FIG. 4. li a

t:0.1 a 10 -3

•

-3 t=o.exto,

2.

I.S

1.0

o.

o.o

-O.

-1.0

FIG. 4.11 b

90

-3 t:1.1 1t 10 1

-3 t:1.2 110 1

-3 hl.!xlO 1

-3 f:L4xl0 1

9 1

DIAGRAMAS 00 ESFORÇO NORMAL AO LONGO DA HASTE

,.e

FIG. 4.12 a

Mana Discreta W=0.00001 ,

-3 t:o.4xl0 1

f=0.71 IO -!.

-3 t:o.e x 10 ,

f:0.9 1 10 -! a

f:t,0 x 10 -! 1

I.&

'·º

1.0

-O.li

LO

O.&

- 0.!11----=--

-1.0

FIG. 4.12 b

92

-3 t:1.t a 10 •

-3 1::1.2 1 10 1

-3 t:1.3 110 1

-3 hl.4 & 10 1

93

DtAG~AMAS DO ESFORÇO NORMAL AO LONGO DA HASTE

1.e

1.e

1.0

o.•

FIG. 4. 13a

MaHa O..creta ,.,= 0.00002 •

.3 t:Q.4 ll 10 1

f:Q.7 x 10 -3

1

t=o.e x ,o-3,

f:Q.9 1C 10 -3

1

.3 f:1.0 x 10 1

I.&

1.0

1.0

-0.&

1.0

• 1.0

FIG. 4.13b

94

-3 t:1.1110 1

-3 t:1.2 1 10 1

-3 t:1.3 1 10 1

-3 t=IAalO 1

95

EXEMPLO II

Considerou-se o comportamento aio-linear de uma ·vi

ga bi-engastada, estudada na referência [36] e outras, submetida

a uma carga concentrada vertical aplicada subitamente no meio do

vao.

Dada a simetria do problema, apenas a metade da vi

ga e do carregamento foram utilizados.

A discretização espacial foi feita atravês de qua­

tro elementos de pórtico plano, de mesmo comprimento, conforme es

quematizado na figura 4.14.

((t) b=I.O ln ~.,_-------------~

e 2 E=O.I x 10 lty'in

)A::0.7971 lb ,r1n4 F( lb l

320~.

o h=o.21n 2 a e

->------ 10 ln ~-----+

FIG. 4.14

Os períodos máximo ~ . e m1n1mo do modelo discreto (fi

gura 4.14) em sua configuração inicial, para os casos de

consistente e massa discreta, levando-se em conta ou nao a

eia de rotação, estao indicados na Tabela 4.1.

massa

inér-

96

.Tf .Tmi'n

S/tn;Rot .. . . -1

0.4983. -4

MASSA 0.1737 .. 10 .. 10 ..

CONSISTENTE .C/In.Rot .. 0.1737 10- 1 0.4983. 10- 4

. S/In.Rot .. .0 .. 1738 . -1 . 0.7590 10- 4 MASSA 10 ..

DISCRETA C/In.Rot.. 0 ... 1738 -1 .0 ... 7 590 10- 4

. . .1.0 .. •,. ..

TABELA 4.1

Dada a pequena influência da inércia de rotaçao so

bre as frequências, nesse problema, ela nao foi

análise que se segue.

considerada na

Adotando-se o intervalo de integração T = 0.00001s

os valores de T/Tf e

estão na tabela 4.2.

T/T . min

MASSA CONSISTENTE

MASSA DISCRETA

Para os dois modelos

T/Tf T/T . min

0.576 10-3 0.201

0.575 10- 3 0.132

TABELA 4.2

discretos

Para se observar o -comportamento da matriz de mas-

sa discreta empregada, eq. (4.16), o deslocamento vertical do no

5 é apresentado na figura 4.15 para os casos de matrizes de massa

consistente e discreta. Empregou-se T = 0.00001s e TIPAN = 3 ,

ou seja, as forças elásticas ao início de um intervalo de tempo

são determinadas a partir da configuração deformada ao final do

97

intervalo anterior.

Usando-se a matriz de .consistente e

T = 0.00001s, a figura 4.16 mostra o deslocamento vertical do no

5 para os casos de TIPAN = 3 e TIPAN = 2 , quando as forças elás-

ticas ~E são calculadas por acumulação da parcela K* bu

-o com o

valor de no inicio do intervalo de integração (vide fluxogr~

ma) .

Os deslocamentos verticais do no 5 ' tomando-se

T = 0.000025s, são práticamente coincidentes com os obtidos com

T = O.OOOOls,impedindo uma representação gráfica clara.

da forma

e ct M onde

Finalmente introduziu-se amortecimento

M é a matriz de massa consistente e a= 200.

O valor de ct adotado corresponde, na análise li-

near, a frações do amortecimento crítico, equação (4.43a), de

!; = 0.276 1 e /;10 = 0.793 10-

3 no primeiro e Último (dêcimo) mo

dos, respectivamente. Os pontos

presentativos da integração dos modos 1 e 10 sao:

pl F; 1 0.276 81 T 0.576 10-3

(Tl Tf) = ... = Tl

plO /; 1 O 0.793 10- 3 8 10

T 0.201 (T 10 T . ) = = TlO m1n

Pela figura 3.3 verifica-se que P1 e estao

na região de con~ergencia e o intervalo de integração escolhido

pode ser usado.

Na análise não-linear a ~ . convergenc.1.a do processo

iterativo estará assegurada enquanto a matriz de iteração, equa -

çao (3.71), variável no tempo, tiver seu raio espectral

pela unidade.

limitado

98

Adotou-se matriz de massa consistente e TIPAN = 3

e o deslocamento vertical do nÕ 5 estâ mostrado na figura 4.17.

.Em todos os casos fez-se TOLl = TOL2 -4 0.1 10 .

DESLOCAMENTO

v9 (tnl

0.9

o.a

0.7

o.e

o.e

0.4

03

0.2

0.1

99

' TRANSVERSAL DO NO 5

Mlaua CoHllttntt

Ma11a Discreta

,.., __ _

' ' ' ' ' ' '

x 10-2$ o ""'---~-----r---~---~---~---~--~

0.1 0.2

V 9

( tn >

0.9

o.a

0.7

o.e

o.e

0'4

03

0.2

0.1

,' ', ' ' ' '

0.3

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

FIG, 4.le

' \ \

.,

0.4 0.15

TIPAN : 3

TIPAN = 2

,'

o.a

' /

0.7

-2 x 10 S o "'-----,----~----~----~----~----~---~

0.1 0.2 o.a 0.4 o.e 0.6 0.7

FIG. 4.16

100

DESLOCAMENTO TRANSVERSAL DO NÓ !5

s' 2 O O M ~

ve (ln)

Q9

o.e

0.7

0.6

o.e

0.4

0.3

0.2

0.1

-2

''º s o 0.2 0.4 0.6 0.8 1.0 1.2 1.4

FIG. 4.17

101

EXEMPLO III

Foi estudada a estrutura de cabos apresentada em

[z] e esquematizada na figura 4.18 submetida

transversal uniforme aplicado subitamente.

a um

Dada a

carregamento

simetria do

problema, analisou-se um octante da estrutura conforme indicado

na figura 4.18.

y

y

X

3 6 'º X

FIG. 4.18

A tabela 4.3 dá as coordenadas dos nos do modelo

discreto bem como a massa concentrada e o carregamento segundo a

direção z , associados a cada um deles.

tabela 4.4.

rentes.

As características dos elementos encontram-se na

Todos os parâmetros estão expressos em unidades coe-

102

N.Õ .X y z .. - _m F z

·1 O.DOO O.DOO o.o 0.6029 10- 2 312.5

.2 200.150 200.150 o.o 0.2408 10- 1 1250.0

3 200.556 0.000 o.o 0.2410 10- 1 1250.0

-. 4 401.050 401.0SO o.o 0.1929 10- 1 1250.0

5 400.299 200.522 o.o 0.4097 10- 1 2500.0

6 401.115 0.000 o.o 0.1989 10- 1 1250.0

7 600.000 600.000 o.o - -8 521.365 399.577 o.o 0.2594 1250.0

9 480.355 199.932 o.o 0.2483 1250.0

10 461.284 0.000 o.o 0.1226 312. 5

TABELA 4.3 - COORDENADAS, MASSAS NODAIS E CARREGAMENTO

ELEMENTO NÕ l NÕ 2 EA TRAÇÃO INICIAL

l 3 2 0.210 10 7 0.1900 10 4

2 6 5 0.210 10 7 0.3093 10 4

3 5 4 0.210 10 7 0.3110 10 4

4 10 9 0.210 10 8 0.1708 10 5

5 9 8 0.210 10 8 0.1734 10 5

6 8 7 0.210 10 8 O. 1821 10 5

7 4 8 0.210 10 7 0.3161 10 4

8 2 5 0.210 10 7 0.1892 10 4

9 5 9 0.210 10 7 0.1868 10 4

10 l 3 O. 105 10 7 0.1605 10 4

11 3 6 0.105 10 7 0.1609 10 4

12 6 10 0.105 10 7 0.1622 10 4

TABELA 4. 4 - CARACTERÍSTICAS DOS ELEMENTOS

O modelo discreto assim definido apresenta,em sua

configuração inicial, um período mâximo Tf = 0.7459s e um perío-

do mínimo T • m1n

-2 = 0.5479 10 s •

O intervalo de integração escolhido foi T

-2 que corresponde a 0.134 10 Tf e 0.1825 Tmin

O.OOls

103

A figura 4.19 apresenta o deslocamento vertical do

no 1 para os dois tipos de análise, T'ÍPAN = 2 e TIPAN = 3

do se utilizam elementos de treliça.

quan-

Na figura 4.20 a mesma comparaçao ê feita, utili -

- -zando-se elementos de cabo que nao atuam em compressao.

A diferença entre as duas soluções ê maior no se-

gundo caso jâ que o procedimento de cálculo das forças elásticas

~E correspondente a TIPAN = 2 introduz erros mais sensíveis sem-

pre que algum elemento deixar de ser tracionado. Quando isto ocor

re a matriz de rigidez do elemento nao contribui para a rigidez

da estrutura. Por outro lado a acumulação do produto K* /iu -o

ao

vetor !_lE no início do intervalo faz com que os elementos nao

tracionados continuem contribuindo para o vetor global !_lE

Finalmente na figura 4.21 o deslocamento vertical

do nÕ 1 ê apresentado para os modelos obtidos com elementos de tre

liça e cabo, utilizando-se o cálculo mais preciso (TIPAN = 3) das

forças elásticas ~E

O comportamento dos modelos ê idêntico atê que se

desenvolva compressão em algum elemento. A partir daí o elemento

de cabo fornece uma aproximação mais realista, desde que se use o

procedimento correspondente a TIPAN = 3

Adotou-se em todos os casos TOLl TOL2 -4

0.1 10 .

104

DESLOCAMENTO TRANSVERSAL DO NO

200 •,

Elem•nto1 de Treliça

180

160

140

120

100

80

60

40

20

o 0.1 02 0.3 0.4

FIG. 4.19

20 O •, Elementos de Cabo

180 ,, ' ' 1 ' 160 1 ' I ' ' 1

140 ' 1 I ' 1 ' I ' 120 1 ' ' ' ' 'J ' ' '

100

' 1

' 1 1 80 ' ' 1

' ' ' 1

1

60 1

' ' 1 1

40 1 1 1

' ' 20 I

' ~

o 0.1 o.2 o.a 0.4

_ 20 FIG. 4.20

0.5

' I

' ' ' ' 1

' ' ' t·.J 1

' ' ' ' '

1

' ' ' ' ' 1' 1' V

0.5

---- TIPAN = 2

-TlPAN=3

' . ' ,_,

0.6

,. '' ' '

' ' ' \,

---- TI PAN

--TIPAN

,''\

' ' '

= 2

= 3

' ' ' ' ' ',, I '·' 1

' 1

' ' 1

' ' 1 1

' ' ' ' ' 0.6

,,

T (s) 0.7

T (s)

0.7

105

' DESLOCAMENTO TRANSVERSAL DO NO

200 w,e

180

1&0

140

120

100

80

60

40

20

' ' -o

0.1 0.2 0.3 0.4

FIG. 4. 21

Treliça

Calo,

' ' '

' /

' ' ' ' ' ' ' ' '

0.15

,·, ' ' 1

,.,

' ' '

0.6

' '·

T (s)

07

106

EXEMPLO IV

A fim de se analisar um problema de .carga mÕvel

[40, 41, 42, 43, 44], não levando em conta a inércia do carrega -

mento [43, 45], mas considerando a não-linearidade geométrica, to

mou-se a viga bi-engastada estudada no ~xemplo II.

A discretização espacial foi feita por meio de 4

elementos de pÕrtico plano de igual comprimento, conforme indica-

do na figura 4.22, e usou-se a matriz de massa consistente.

fundamental

2 3

20 ln

FIG. 4. 22

4

O modelo discreto assim obtido apresenta o período

Tf = 0.01735s e o período ~ . m1.n1mo T. = 0.9965 m1n·-

-4 10 s,

não tendo sido considerada a inércia de rotaçao.

A carga mÕvel, de 500 ib, percorreu a viga com ve-

locidade constante definida pelo parâmetro a= T /T f trav. onde

T e o tempo gasto pela carga para atravessar toda a trav ext ens ao

da viga.

A tabela 4.5 indica os valores de a adotados, os

tempos de travessia e as velocidades da carga (in/s) corresponde~

tes.

107

.a. T trav .V

2.0 0.008675 2305.4735

1. 5 0.011567 1729.1066

1. o 0.01735 1152. 7377

·º· . .5. 0.034.70 .5.7.6. 36.89

TABELA 4.5 - TEMPOS DE TRAVESSIA E VELOCIDADES

O intervalo de integração usado foi -4 T = 0.2313 10 S

correspondente a 0.133 10- 2 Tf e 0.232 T . min

Os resultados da análise linear estao na figura

4.23 que fornece, em ordenadas, o deslocamento vertical do no cen

tral e, em abcissas, o tempo dividido pelo tempo de travessia cor

respondente ao valor de a. adotado.

Os resultados da anãlise não-linear, através do

procedimento correspondente a TIPAN = 2, encontram-se na figura

4.24.

Observa~se que o deslocamento miximo se deu para

a.= 2.0 que corresponde a uma velocidade maior do que aquela que

dão mãximo no caso linear (a.= 1.0)

Adotou-se neste exemplo TOLl = TOL2 -5 0.1 10 .

108

' DESLOCAMENTO VERTICAL 00 NO 3

v3 <1•1 c,l.:1,0

Resposta Ltneor

5.0-l

i 4.0j

' 3.0

2.

1.0

o.o T

T tra V.

-1.0

-ll.Oj

FIG. 4.23

V (ln) Resposta Não-Lrn•or 3

0.8 c,l.:1.5

0.6

0.4

0.2

T o.o

-0.2

-0.4

-o.e FIG. 4.24

109

4. 9 - CONCLUSÕES

O operador de integração discutido, em conjunto

com o método dos elementos finitos, constituem ferramenta

de tratar o problema da resposta dinâmica de estruturas.

capaz

No domínio linear o algorítmo implementado nao

parece ser, no caso geral, competitivo, do ponto de vista de tem­

po de processamento, por ser iterativo e ter o intervalo de inte~

gração limitado pela condição de convergência desse processo ite­

rativo, fig. 3.2.

O algoritmo modificado, incondicionalmente está­

vel, permite intervalos de integração arbitrários porem envolve a

solução direta de um sistema não-simétrico, de ordem 2n, em cada

. intervalo. Embora preciso, tal esquema ainda não foi aplicado a

problemas medias e grandes.

No domínio não-linear, o uso da matriz K* -o per-

mite, em geral, o calculo das forças elásticas por acumulação. Em

alguns casos, entretanto, como estruturas de cabos em que algum

elemento deixa de ser tracionado, o calculo das forças

devera ser feito a partir da configuração deformada.

elásticas

O uso da matriz tangente modificada K* -o e das

expressoes (3.64) fornecem, para os intervalos usualmente empreg~

dos [47], uma boa aproximaçao.

Em ambos os domínios, a maior eficiência dos al­

goritmos se dá em problemas de transientes curtos em que a contri

buição dos modos mais altos ê significativa, larguras de banda re

lativamente grandes e onde o uso de matriz de massa diagonal tor-

na trivial a solução do sistema de equações, independentemente da

forma da matriz de amortecimento adotada.

110

A convergencia do processo iterativo e definida,

para o caso de amortecimento proporcional, pelas figuras 3.2 e

3. 3. Esses gráficos possibilitam uma aproximação para o caso de

amortecimento não-proporcional sendo Úteis também em problemas

nao lineares, onde a matriz ·de iteraçio e seccionalmente constan-

te.

A introdução de amortecimento dificulta a conver

gencia do processo iterativo sendo, por vezes, necessário usar um

intervalo de integração menor do que aquele do problema sem amor-

te cimento. Trata-se de uma característica inerente ao processo

iterativo de resolução e não ao algorítmo de integraçao em si.

A discretização espacial devera ser cuidad6sa,

evitando um refinamento desnecessário que daria origem a uma fre­

quência de corte elevada e, consequentemente, um intervalo de in­

tegraçao muito pequeno, onerando a solução.

~ interessante se ter pelo menos uma estimativa

do período mínimo do modelo discreto para, com auxílio das figu-

ras 3.2 e 3.3, se fazer a escolha do intervalo de integraçao.

Na sua regiao de utilização, em termos de inter­

valos de integração, e para os problemas analisados aqui e em ~71, o algoritmo não se mostrou muito sensível a variaçoes no interva-

lo de integraçao.

O uso da matriz tangente K -o

em problemas nao-

lineares é uma alternativa que merece ser analisada no que diz· res

peito ã precisao e tempo de processamento.

A pesquisa de técnicas numéricas mais eficientes,

para fazer frente a novas necessidades, continua e é atestada pe­

lo grande volume de trabalhos na área.

111

BIBLIOGRAFIA

1 ARGYRIS, J.H. and eHAN, A.S.L. - Applications of Finite Ele-

ments in Space and Time - Ingenieur - Archiv 41, 235-

257, 1972.

2 ARGYRI S, J. H. ; DUNNE, P. e. and ANGELOPOULOS, T. - Non-Line-

3

ar Oscillations Using the Finite Element Technique - eo~

puter Methods in Apllied Mechanics and Engineering,Vol.

2, 203-250, 1973.

ARGYRIS, J.H.; DUNNE, P.e. and ANGELOPOULOS, T. - Dynamic

Response by Large Step Integration - Earthquake Engineer­

ing and Structural Dynamics, Vol. 2, 185-203, 1973.

4 GELFAND, I.M. and FOMIN, S.V. - Calculus df Váriations

Prentice-Hall, Inc., 1963.

5 MEIROVITCH, L. - Analytical Methods 1n Vibrations - The Mac-

Millan eompany, 1967.

6 LOULA, A.F.D. e GALEÃO, A.e. - Vibração de Sistemas Elásti-

cos Lineares - PDD 3/76, eOPPE/UFRJ, 1976.

7 BREBBIA, e.A. and FERRANTE, A.J. (Editors) - The Finite Ele-

ment Technique - Editora da UFRGS, 1975.

8 FEIJOO, R.A. - Aplicacion del Metodo de Ritz a Funcionales

Relajados en Mecanica de los Solides

eOPPE/UFRJ, 1973. Tese de M.Sc.,

9 eLOUGH; R.W. - Analysis of Structural Vibrations and Dynamic

Response - Recent Advances in Matrix Methocls of Structu­

ral Analysis and Design - U.S.-Japan Seminar, Univ. of Alabama Press, 441-486, 1969.

10 NieKEL, R.E. - Direct Integration Methods 1n Structural Dy­

namics - ASeE, Vol. 99, No. EM2, 303-317, April, 1973.

11 DUNHAM, R.S.; NieKELL, R.E. and STieKLER, D.e. - Integration

Operators for Transient Structural Response - eomputer

& Structures, Vol. 2, 1-15, 1972.

12 eLOUGH, R.W. and MOJTAHEDI, S. - Earthquake Response Analy­

sis Considering Non-Proportional Damping Earthquake

Engineering and Structural Dynamics, Vol. 4, 489-496, 1976.

112

13 NICKELL, R.E. - Nonlinear Dynamics by Mode Superposition

Computer Methods in Applied Mechanics and Engineering,

7, 107-129, 1976.

14 BATHE, K.J. and WILSON, E.L. - Numerical Methods in Finite

Element Analysis - Prentice-Hall, 1976.

15 CLOUGH, R.W: and PENZIEN, J. - Dynamics of Structures

Graw-Hill Inc., 1975. Me

16 GUIMARÃES, S.S. - Análise da Estabilidade de Métodos de Inte

graçio Passo a Passo - Tese de M.Sc., COPPE/UFRJ, 1977.

17 NICKELL, R.E. - On the .Stability of AppI'.oximation Operators in

Problems of Structural Dynamics - International Journal

of Solids and Structures, Vol. 7, 301-319, 1971.

18 BATHE, K.J. and WILSON, E.L. - Stability and Accuracy Analy­

sis of Direct Integration Methods - Earthquake Enginee~

ing and Structural Dynamics, Vol. l, No. 3, 283-291,1973.

19 GOUDREAU, G.L. and TAYLOR, R.L. - Evaluation of Numerical Me

thods in Elastodynamics - Computer Methods in Applied

Mechanics and Engineering, Vol. 3, 69-97, 1972.

20 HILBER, H.M.; HUGHES, T.J.R. and TAYLOR, R.L. - Improved Nu­

merical Dissipation for Time Integration Algorithms in

Structural Dynamics - Earthquake Engineering and Struc­

tural Dynamics, Vol. 5, 283-292, 1977.

21 KEY, S.W. and BEISINGER, Z.E. - The Transient Dynamic Analy­

sis of Thin Shells by the Finite Element Method - 479-

518, Proceedings of the Third Conference on Matrix Me­

thods in Structural Mechanics, Wright-Paterson Air For­ce Base, Ohio, 1971.

22 HINTON, E.; ROCK, T.A. and ZIENKIEWICZ, o.e. - A Note on Mass

Lumping and Related Processes in the Finite Element Me­

thod - Earthquake Engineering and Structural Dynamics, Vol. 4, 245-249, 1976.

23 ODEN, J.T.; KEY, J.E. and FOST, R.B. - A Note on the Analy­sis of Nonlinear Dynamics of Elastic Membranes by the

Finite Element Method - Computers & Structures Vol. 4, 445-452, 1974.

113

24 SHANTARAM, D.; OWEN, D.R.J. and ZIENKIEWICZ, o.e. - Dynamic Transient Behaviour of Two and Three-Dimensional Structures In­

cluding Plasticity, Large Deformation Effects and Fluid

Interaction - Earthquake Engineering and Structural Dy­

namics, Vol. 4, 561-578, 1976.

25 NEWMARK, N.M. - A Method of Computation for Structural Dyna­

mics - ASCE, Vol. 85, No. EM3, 67~94, 1959.

26 NAGARAJAN, S. and POPOV, E.P. - Elastic-Plastic Dynamic Ana­

lysis of Axisymmetric Solids - Computers & Structures,

Vol. 4, 1117-1134, 1974.

27 WILSON, E.L.; FARHOOMAND, I. and BATHE, K.J. - Nonlinear Dy­

namic Analysis of Complex Structures - Earthquake Engi­

neering and Structural Dynamics, Vol. l, 241-252, 1973.

28 McNAMARA, J.F. and MARCAL, P.V. - Incremental Stiffness Me­

thuds for Finite Element Analysis of the Nonlinear Dyna­

mic Problem - Numerical and Computer Methods in Structu­

ral Dynamics, 353-376, Symposium held at Urbana, Illi­

nois, Sept. 1971.

29 ARCHER, J.S. - Consistent Mass Matrix for Distributéd

Systems - ASCE, Vol. 85, No. ST4, 161-178, 1963.

Mass

30 EBNER, A.M. and UCCIFERRO, J.J. - A Theoretical and Numeri­

cal Comparison of Elastic Nonlinear Finite Element Me­

thods - Computers & Structures, Vol. 2, 1043-1061, 1972.

31 WILSON, E.L. and PENZIEN, J. - Evaluation of Orthogonal DamE

ing Matrices - International Journal for Numerical Me­

thods in Engineering, Vol. 4,, No. 1, 5-10, 1972.

32 JONES, A.T. and GRANT, J.E. - Generation and Properties of

Banded Viscous Damping Matrices - Computers & Structu­res, Vol. 3, 1037-1051, 1973.

33 THOMSON, W.T.; CALKINS, T. and CARAVANI, P. A Numerical Study of Damping - Earthquake Engineering and Structu­ral Dynamics, Vol. 3, 97-103, 1974.

34 BERGAN, P.G. and CLOUGH, R.W. - Convergence Criteria for !te

rative Processes - A.I.A.A. Journal, Vol. 10,

1107-1108, 1972.

No. 8,

35

114

TONG, P.; PIAN, T.H.H. and BUCCIARELLI, L.L. - Mede Shapes

and Frequencies by Finite Element Method Using Consistent

and Lumped Masses - Computers & Structures, Vol. l, 623

638, 1971.

36 WEEKS, G. - Temporal Operators for Nonlinear Structural Dyna

mies Problems - ASCE, EM~5, 1087-1104, October, 1972.

37 FERRANTE, A.J. - Notas de Aula de Mecânica das Estruturas II

- COPPE/UFRJ, 1975.

38 PRATES DE LIMA, E.C. - LORANE DINA - Uma

para Análise Dinâmica de Estruturas

COPPE/UFRJ, 1977.

Linguagem Orientada

Tese de D.Se.,

39 EBECKEN, N.F.F. - LORANE-NL - Uma Linguagem Orientada à Aná­

lise Estrutural Não-Linear - Tese de D.Se. - COPPE/UFRJ,

1977.

40 VENÂNCIO FILHO, F. - Dynamic Influence Lines of Beams and Frames - ASCE, No. 92, ST-2, 371-386, April, 1966.

41 BRUCH, Y.A. - Análise Dinâmica de Placas pelo Método dos Ele mentes Finitos - Tese de M.Sc. - COPPE/UFRJ, 1973.

42 JUN"G, M.P. - Dinâmica de Estruturas Reticuladas sob Cargas

Móveis pelo Método dos Elementos Finitos - Tese de M.Sc.

- COPPE/UFRJ, 1973.

43 FALABELLA, J.E. - Resposta Dinâmica de Estruturas Reticula­

das a Cargas Móveis pelo Método dos Elementos Finitos -Tese de M.Sc. - COPPE/UFRJ, 1975.

44 MIANA, P.R. - Análise Dinâmica de Vigas de Pontes pelo Méto­do dos Elementos Finitos - Aplicação a Vigas Gerber Tese de M.Sc. - COPPE/UFRJ, 1975.

45 GUERREIRO, J.N.C. - Dinâmica de Estruturas Prismáticas Lami­nares - Consideração de Cargas Móveis - Tese de M.Sc. -COPPE/UFRJ, 1977.

46

47

ALBRECHT, P. - Análise Numérica: Um Curso Moderno

Técnicos e Científicos Editora S.A. - Rio, 1973.

Livros

CORREA BARBOSA, H.J.; MAIA DA COSTA, A.; EBECKEN, N.F.F. e

VENÂNCIO FILHO, F. - Métodos de Integração Direta para a Análise Dinâmica Não-Linear - ·XIX Jornadas Sudamerica­

nas de Ingenieria Estructural, Santiago, Chile, 1978. A ser pu­

blicado.

115

APllNDICE

a) CARACTERÍSTICAS GERAIS DO PROGRAMA

O programa, desenvolvido em FORTRAN IV, segue o

fluxograma mostrado no item IV.7 e os procedimentos relativos a

montagem de matrizes globais e introdução de condições de contor­

no foram baseados no fluxograma apresentado em [37]

As condições de contorno incluem o caso de apoios

inclinados com relação ao referencial global. Todas as

são levadas em conta na integração ao longo do tempo,

aquelas correspondentes às direções impedidas.

As condições iniciais são supostas nulas:

equaçoes

inclusive

u (O) =

= ;,_(O) = O Entretanto, condições mais gerais poderão ser incor

paradas sem maiores dificuldades.

As características dos elementos implantados enco~

tram-se nos itens IV.2 e IV.3 salientando-se que foram despreza­

das a deformação por corte e a variação da seçao transversal e

que as deformações sao supostas pequenas.

O vetor de solicitações é montado automaticamente,

para o caso de carga mõvel, conforme indicado no item IV.5. Pro-

cedimentos mais gerais para esta montagem, como os empregados no

sistema LORANE [38, 39], podem ser implementados, de acordo com

as necessidades do usuârio.

b) SUBROTINAS EMPREGADAS

1 - Subrotina DADOS

Lê e imprime todos os dados de entrada do programa.

Calcula e imprime a largura de banda e o com~rimento

de cada elemento.

inicial

116

2 - Subrotina ROT

Forma a matriz de rotaçao de ordem (3x3) que e usada para a

introdução de apoios inclinados e mudança de referenciais (lo

cal e global).

3 - Subrotina ROTM

Transforma as matrizes do elemento (6x6) do referencial local

para o referencial global a partir da matriz de rotação obti­

da pela subrotina ROT.

4 - Subrotina MULTV

Efetua a multiplicação de urna matriz

um vetor.

A de ordem (NlxN2) por

Caso o parametro ITRANS seja diferente de zero e a matriz A

seja quadrada, serã obtido o produto de AT pelo mesmo vetor.

5 - Subrotina MULTM

Efetua o produto matricial

tangulares e conformes.

6 - Subrotina SACAP e BOTAP

c A B onde as matrizes sao re-

São subrotinas auxiliares para retirar e colocar, respectiva­

mente, vetores, correspondentes a um nó particular, em um ve­

tor global da estrutura [37]

7 - Subrotina SACA e BOTA

São subrotinas auxiliares para retirar e colocar, respectiva­

mente, uma submatriz de uma matriz global da estrutura [37]

8 - Subrotina MONTA

Forma as matrizes globais de massa e ·rigidez da estrutura com

a semi-banda armazenada em forma unidimensional [37].

117

9 - Subrotina EXCITA

E uma subrotina,construÍda para cada caso particular,que for

ma os vetores globais de e~citação da estrutura, "F", e sua

derivada em relação ao tempo "D.F "

10 - Subrotina APOIO

Introduz as condições de contorno, do tipo deslocamento pre~

crito nulo, segundo eixos que podem ser inclinados em rela­

ção ao referencial global [37]

11 - Subrotina TRIANG

Triangulariza uma matriz quadrada simétrica com sua semi-ban

da armazenada em forma de vetor.

12 - Subrotina RESOB

Resolve o sistema de equaçoes lineares, apos a triangulariz~

ção da matriz dos coeficientes pela subrotina TRIANG, para

vários vetores independentes armazenados sequencialmente.

13 - Subrotina RESOD

Resolve diretamente o sistema de equaçoes lineares

trivial em que a matriz dos coeficientes ê diagonal.

14 - Subrotina GLOBAL

no caso

Transforma os deslocamentos, referidos ao sistema local, de

um apoio inclinado para o referencial global.

15 - Subrotina MULTB

Efetua o produto de uma matriz quadrada simétrica, com sua

semi-banda armazenada em forma unidimensional, por um vetor.

118

16 - Subrotina CONV

Verifica se o critêrio de convergência adotado jâ foi satis­

feito.

Vide item IV.6.

17 - Subrotina MASSP

Forma a matriz de massa (consistente ou diagonal) do elemen-

to de pórtico plano no referencial local e transforma

subrotinas ROT e ROTM) para o referencial global.

18 - Subrotina MASSC

(via

O mesmo que a subrotina MASSP para os elementos de treliça e

cabo,planos ou espaciais.

19 ~ Subrotina RIGEP

Forma a matriz de rigidez usual, linear, e a matriz geomêtri

capara o elemento de pórtico plano, no referencial local, e

transforma para o referencial global.

Calcula o vetor de forças elásticas nodais no referencial lo

cal.

No caso de TIPAN = 3, o vetor global de forças elásticas e

montado por acumulação das contribuições dos elementos.

20 - Subrotina RIGEC

O mesmo que a subrotina RIGEP para os elementos de treliça e

cabo, planos ou espaciais.

21 - Subrotina SAIDA

Imprime os deslocamentos, velocidades e acelerações dos nos

selecionados e os esforços nos nÕs dos elementos

dos, nos intervalos de tempo especificados.

seleciona-

119

22 - Subrotina ESFLIN

Calcula e imprime, para cada el·emento selecionado, os esfor­

ços seccionais nos nôs no caso de análise linear~

c) MANUAL DE ENTRADA DO PROGRAMA

N9 DE N9 DE VARIÁVEIS FORMATO ORDEM CARTÕES

1 1 NNOS,NELEM,NAP,NNED,TIPAN,ESTAB,CARREG,MASSA 815

2 l INROT,IAMOR,ALFA,BETA,TOL1,TOL2,TF,TAU,IMPRES 2I5, 6Fl0.0, IS

3 1 PV, VI, A, NEP 3Fl0.0, IS

4 NNOS I, X(I), Y(I), Z (I), IMPDVA(I) IS, 3Fl0 . O, IS

5 NELEM I,ITIPO(I),CONEC(I,l),CONEC(I,2),E(I),MESP(I),

AX(I), IZ(I), FZERO(I), IMPFOR(I) 4IS, SFlO.O, IS

6 NNED NO(I), (ED(I,J), J = 1,9) IS, 9F8.0

7 NAP (MCC(I,J), J = 1,5), (RLOCAL(I, J), J = 1,9) SIS, 9F6.0

8 IT ~ 8 C (I) I = 1, IT IT = NNOS * NGLN * LBNOS * NGLN 8Fl0.0

DEFINIÇÕES E COMENTÁRIOS

1 - NNOS - N9 de nos da estrutura

NELEM - N9 de elementos

NAP - N9 de nos com alguma direção impedida

NNED - N9 de nos com algum elemento discreto

TIPAN - Indica o tipo de análise a ser efetuado:

TIPAN = 1 - Análise linear

TIPAN = 2 - Análise não linear onde o vetor de forças

elásticas ao final do intervalo de inte -

gração e calculado por acumulação de

K* 6u_ .. -o

120

TIPAN = 3 - Análise nao linear onde o vetor de forças

elásticas ao final do intervalo de inte -

gração e calculado a partir da configura­

ção deformada.

ESTAB - Indica os coeficientes da expressão dos deslocamentos

(3.45b) ou (3.56b).

ESTAB = O - São usados os coeficienies correspondeu-

tes ao algoritmo condicionalmente estável

ESTAB

(3.45b).

1 - são.usados os coeficientes correspondeu-

tes ao algoritmo incondicionalmente está­

vel (3.56b).

Como ê usado o processo

ESTAB = O [2]

iterativo deve-se adotar

CARREG - Igual a 1, indica carregamento qualquer definido pela

subrotina EXCITA.

Igual a 2 para o caso de carga móvel como definida em

IV. 5

MASSA - Controla a montagem da matriz de massa. Vide item

IV. 2. 1.

2 - INROT - Controla a inclusão do efeito de inercia de rotaçao

!AMOR,

ALFA

na matriz de massa. Vide Ítem IV.2.1.

BETA - Parâmetros que controlam a construçao da matriz de

TOLl,

TOL2

amortecimento. Vide Ítem IV.4.

Tolerâncias para os testes de convergencia respecti~~

mente do vetor de deslocamentos e do vetor de deriva-

das das forças de inercia.

121

TF - Tempo atê o qual se quer a resposta da estrutura.

TAU - Intervalo de integração

IMPRES - Intervalo de impressão dos resultados

3 - PV - Valor da carga móvel concentrada vertical

VI - Velocidade inicial da carga móvel

A - Aceleração constante da carga móvel

NEP N9 total de elementos percorridos pela carga

No caso de CARREG = 1 este cartao e deixado em branco.

4 - I - Número do no

X( I),

y ( I) ,

z ( I) - Coordenadas do no I. Para estruturas planas, Z = O.

IMPDVA(I)- Deve ser feito igual a 1 quando ·se deseja imprimir

o deslocamento, velocidade e aceleração do nó I.

5 - I - Número do elemento

ITIPO(I) - Tipo do elemento I. Vide Ítem IV.2.

CONEC(I, 1)

CONEC(I,2)-NÚmero dos nos 1 e 2 do elemento I .

E(I),

MESP(I) - Módulo de Young e massa específica do elemento I.

AX ( I) ,

IZ(l) Ãrea e momento de inércia da seçao transversal do elemento I.

FZERO(I) - Esforço normal inicial no elemento I, positivo se for de traçao.

H1PFOR(.I)- Deve ser feito igual a l quando se deseja imprimir

os esforços no elemento I

6 - Cartões fornecidos apenas no caso de NNED > 1

NO (1) Número do nÕ com elemento discreto

ED(I,J) - Matriz das constantes dos elementos discretos. Vide

IV. 3.

7 - MCC(I,J)

122

- Matriz de condições de contorno

MCC(I,l) - Numero do no com restrição

MCC(I,2) - Quando igual a 1 indica restrição na

direção 1

MCC(I,3) - Quando igual a 1 indica restrição na

direção 2

MCC(I,4) - Quando igual a 1 indica restrição na

direção 3

Quando a direção ê livre o valor for

necido deve ser zero.

MCC(I,5) - Quando diferente de zero indica que

as direções acima não são globais e

os deslocamentos locais

são relacionados por

do a matriz de rotação

em RLOCAL.

~L R

e globais

= R ~G s en

contida

RLOCAL(I,J) - Contêm na linha I os elementos da matriz de ro-

tação do apoio inclinado

do suas linhas.

MCC(I,l) dados segun-

Se MCC(I,5) =O, nao sao fornecidos.

8 - Cartões fornecidos apenas no caso de !AMOR 3

C(I) - Matriz de amortecimento lida conforme item IV.4

d) DIMENSIONAMENTO

O programa listado aqui estâ dimensionado da seguinte

maneira:

- n'? de nos - 100

- n'? de elementos - 100

n'? de ·-nos com apoio 20

- n'? de nos com elementos discretos - 50

123

Os vetores que armazenam os elementos da semi-ban­

da superior das matrizes gleba.is de massa, amortecimento e rigi­

dez, estão dimensionados para:

(NNOS * NGLN) * (LBNOS * NGLN) < 10.000

onde:

NNOS - N9 de nos

NGLN - N9 de graus de liberdade por no, no caso, igual a 3 .

LBNOS - Largura de banda em nôs

e) LISTAGEM

124

SUBROUTINE DADOSCNO,ED 1 MCC,RLOCALl INTEGER CONEC,TIPAN,ESTAS,CARREG REAL MESP,IZ,Ml,M2,KMOLA,LZERO COMMON /GERAL/lR,I",IT 1 NNE,NGLN,NNOS,NELEM,NAP1NNEO,LBNOS,

•TIPAN,ESTAB,INROT,MASSA,IAMOR,ALFA,8ETA,TOLl,TOL2,TF,TAU, •lMPRfS,CARREG

COMMON /PROP/AX(IOO),MESP(IOOJ,CONEC(I00,2),X(IOOJ,YCIOO), *Z(IOO),E(IOOl,IZCIOOJ,ITIPOCIOOl,LZER0(100)1XZERO(l00)1 -*YZEROC!OO),ZZERO(IOO),ANGZC!OOl,FZERO(!OO)

COMMON /IMPRE/IMPSl,lMPOVA(lOO),lMPS2,IMPFOR(IOO) COMMON /MOVEL/ PV,VI,A 1 NEP - . -DlMENSION NO(SO),EDCSO,q),MCCC20,S),RLOCAL(20,qJ ~RITECIW,qOO) .

900 FQRMAT(//,20X,t•UM ALGQRJ~MQ DE ORDEM 5UPERIOR PARA A INTE' •,'GRACAO DIRETA DAS EGUACOES DA DINAMICA ESTRUTUR•L•'l/,lOi •, 1 TESE DE MESTRA00',12X,'HELI0 JOSE CORREA BARB0SA',i2X, . •'COPPE/UFRJ' 1 1/) . .

RE•D(IR,IOO)NNOS,NE.LEM1NAP,NNEO,TIPAN1ESTAB,CiRREG1MASSA 100 FQRMAT(BIS)

WRITE(X~,101)NNOS,NELEM 1 NAP,NNED,TIPAN,ESTA6,CARREG,MASSA !OI FORM,q (/1IOX1 1NNQS! ,SX, INELEMI ,SX, INAPt ,SX, INNE.Q! ,SX, 'TIPAI

• , 1 N 1 , s X, ' E s TA a t , s X' ! e ARRE G 1 , 5 X, ! MA s 5 A ' , / /, 1 o X 1 4 ( Í LI, s X ) , l :s 1 . •9X,l1,1ox,r1,qx,11l

READ(IR,200JINROT1ÍAMOR,ALFA1BETA,TOL1 1 TOL2,TF1TAU1lMPRES 200 FORMAT(2I5,bf10,0,IS)

IFC INROT ,NE, O l INROT:\ xrc I•MoR ,LT, 1 ,oR, IAMQR ,GT, l l lAMoR=t

WRITE(IW,20!)INROT 1 1AMOR,ALFA 1 BETA 1 T0Lt,TOL2,TF 1 TAU,IMPRES 2 O 1 f DR MA T ( / /; I O)(, ' I NR O T ' , S X 1 1 I AMOR 1 , 5 )( 1

1 .A 1.. F A I t \ O X Í ' 6E TA ' , 1 O X, * 'T OL 1 1 , 1 O X 1 1 T OI,. 2 1 , 1 1 X, 1 T F 1 , 1 l X, ! TA Ú t I b X, t X MP RE S t , / ! 1 12 XI II 1 •9X1ll1b(4X,E10,Ll),2X1Ill

READ(IR,300)PV,Yl,A,NEP 300 FORMAT(3f10,0,I5)

IF( CARREG ,EQ, 2 ) WRITE(IW,JO!l PV,VI,A,NEP 30! FoRMÁ,TC/1,lOX,tC AR G A MO V E L',11,tBX,'PV',lBX,IVJ! *, I q X 1 ! A 1 , 1 2 X, ! NE P 1 , / /, 1 O X, 3 ( F 1 5, 4, S X ) ,IS) . . . .

READ(IR,QOO)(I,X(Il,Y(Il,Z(I),IMPDVACIJ,K=t,NNOS) aoO FORMAT(l5,3F10 1 0,I5J . . -

NRITE(lW,401lCI,X(ll,Y(l),z(I),IMPDVA(Il,I=l,NNOSl 401 FoRMAT(l/,10X, 1 C O ORDENA D AS NO D l I S',ll,11X

•,'I',BX, IX(I) 1 ,IOX, 'Y(IJ 1 ,IDX, 1 Z(l)',1tX,1IMPDiA(I)',II, . ~c1ox,13,ax,,10,a,ax,F10,a,ux 1 F10;a;12x,11)l ·

READCIR,500)CI,ITJPO(t),CONECCI,i),CONEC(I,2)1ECil,MESP(l), *AX(l),IZ(I),FZERO(IJ 1 IMPFOR(l),K=l,NELEMj -

soa FORMAT{•IS,SFIO,O,I5) -Ll>NOS:O JFIM=NNE•1 lMP52::o DO 14 I=l,NELEM

125

IMPS2=IMPS2+IMPFORCI) 00 13 J=l,JFtM . !<IN=J+l DO 12 K::KIN,.NNE L=IASS(CONECCI,Jl~CONEC(I,!<)) IF(L6NOS~Lll!,12,12

11 LBNOS=L 12 CONTINUE \3 cuuru,uE liJ CG•illNU~

t6r-J()S:;,.t;N05-+' 1 IMPs1::u DO JS I=!,NNOS lMPSl=IMPSl+IMPOVACI) XZEROCl)=XCI) YZERO(Il=YCil

15 ZZERO(I)~Z(l) DO 1ó I=!,NELEM DX:X(CONEC(I,2))~xccoNEC(J,1ll DV~Y(CONECCI,2))•Y\CONEC(I,ll) DZ=Z(CONEC(I 1 2)l•ZCCONEC(I,!)) LZERO(t):SQRTCON~Dk+DV•OVtDZ~OZ) IFCTIPAN ,GT, !l LZEROCtl=LZERO(Il/C!,+FZERO(I)IAX(t)IE(Il) ~NGZCI)=ATAN2COV,OX) . . . -

16 CONTINUE WRITE(IW,SO!)CI,lTlPO(ll,CONEC(J,ll,CONEC(I,2),E(I),MESPCIJ

"• A.X(IJ I tZ(J) ,FZERO( I J 1 1.ZERO(Í) 1 ANGZ ( I)-,IMPFQR(I) 1 I.:1,NELéMl 501 fORMAT(//,lOX,'C A.R AC TER l S TI C A 5 D OS Ei

*,'LEME N TO S!,//,aX,'I ITIPOCI) NO l NO 2',óXIE(i *,'I)',gX,•MESP(l) 1 18X1 'AX(I) 1 ,QX,'IZCI) 1 17X, 1fZEROCI)'ibX,. ~ 11. ;z E RO ( I l 1 , 7 X, 1 A NG Z (I ) ! , 3 X, ' l M PI' OR ' 1 / I 1 ( 3 X, I l, 3 X, I 2, ó X, I 3 1

•4X,l3,7C4X 1 EI0,4),qX,I1)) IF(NNEDl2,2, 1

1 READ(IR,óOO)(NO(I),(ED(J,JJ,J:1,Q),I.:t,NNED) 600 fORMAT(I5,qfa,Ol

WRITE(IW,b01lCNO(Il,(ED(I1Jl1J=t,Ql,I=1,NNEOl bOl FORMAT(l/,IOX,'E L.E ME N TOS . DISCRETO SI,//,

.. 4 X' 1 NO ' ' 8 )( 1 1 M 1 ! ' 1 2 X, ! M 2 1 , 12 X 1 1 M 3 1 , 12 X' ! e 1 1 , 12 X 1 1 e: 2 1 ' 12 X 1 'e ' f! 1 1 3 ' 1 12 X, ! K I ' , ! 2 X, li( 2 1 1 ! 2 ~, 1 K 3 1 , / / t C 3 X, I 3, 'H l! XI E I O~ ti ) l l . .

2 Rf.AP(IR,700)((MCCC!,J),J;!,S),(RLOCALCI,J),J.:t,q),I:t,NAP) 700 FORMAT(5!5 1 QFb,O) . . . -

WRITE(JW,70l)((MCC(I,J),J;\,Sl,(R~ocAL(I,J),J=1,q),I=l1NAP) 701 FORMAT(//,IOX,'C O N D I COES D E CONTO R N 01

,.,/l,4X 11 NO OIR, 1 OIR, 2 DlR ,l ROT,lr;A0',22X,'M A'

•,'TRIZ D E ROTA C A 0 1 ,//(]X,I3,5x,11,l(8X,Ii ol),4X,9ft0,5))

WRITECIW,800) LBNOS 800 FQRM-TC// 1 10X,ILARGUR4 DE B-NDA EM NOS 1',Il,//)

RETURN END

1 O

2

3

126

SUBROUTINE ROTCR,RLOCAL,CX,CY,CZ,K) DIMENSION R(3,3l,RLOCAL(20,q) . IF( K ,EO, O l GOTO 2 DO 10 1=1,:S RC1,Jl:RLOCAL(K,1l RC2,l)=RLOCAL(K,I+3) R(l,Il=RLOCAb(K,I•ól RETURN Q:SQRTCCX*CX+CZ•CZJ !F( Q ,GT, 0,000) ) GOTO l R C 1 , 1 ) :.o, . RC!,2Jt:CY R(l,3)::0, R(2, 1 l=~CY RC2,2l=Q, R(2,l)"O, R(3,l)=O, R(3 1 2)=0, R(3 1 3)iq, RETURN R(1,!l=CX R(l,2J=CY R(1,3)o:CZ R(2,tl::•C~•CY/Q R(2,2)=0 IH2,3):•CY•CZ/G R(3, 1 ):: .. CZ/0 R(J,2).:;o, R(3,3J:CX/Q RETURN ENO

/

127

SU6ROUTINE ROTM(RTSR 1 S,R) DIMENSION RTSR(b1bl,6Cb,b),R(3,3) DO 10 I=l,3 DO 10 J"l,3 E:: O• · f:O, G:: O, DO 20 N=1,3 RT=R(N, Il DO 20 ~=1,3 E:E+RT•S(N,K)tR(K,Jl F=F•RT•SCN,K+3)•R(K,Jl

20 G=G+RT~S{N+3 1 K+3)*RCK,J) RTSR(I,J):E RTSR(I,J+3l=F

10 RTSR(I+3,J+3)=G DO ,o I=l,S II=I+I DO 30 J::II,b

30 RTSR(J 1 T.):RTSR(I1Jl RETURN -END

SUBROUTINE MULTVCYT!,A,VT2,Nl,N21ITRANS) DIMENSION VTl(Nl),VT2(N2),A(Nl,N21

_JF( l!~ANS ,NE, O ) GD TO 2 DO 10 I:1,NI VTl(I):O, DO 20 J=t,N2

20 VT1(l):VT1tI)+A(I,J)tVT2(J) 10 CONTINUE - .

RETURN 2 00 30 I=l,NI

VTl(Il:o, Do ao J=!,N2

40 VTICI):VTÍ(I)+ACJ,I)~VT2CJ) 30 CONT INLJE

RETURN ENP

128

SUBROUTINE MULTM(C 1 A 1 B1 Nl,N2 1 N3l DIM~NSION C(N!,N3),A(N1,N2),B(N2,N3),V(3) DO 10 l=!,N! DO 20 J:1,Nl V(Jl=O• DO 20 K=l,N2

20 VCJl=V(J)+A(I,Kl•B(K,J) DO 30 J:t,Nl

30 C(I,Jl=V(Jl 10 CONTINUE

RETURN END

SUBROUTINF. SACAP(I 1 P 1 N1 VE·T,NV) COMMON /GERAL/lR,IN,lT1NNE,N;LN,NNOS DIMENSION P(NGLN),VETCNV) Kt:NGLN•(NNOS•(l•ll•N•ll DO 10 K:t,NGLN . KK:1(1+1<

10 P(K):VET(KKl

10

RETURN EtJD

SU8ROUTINE BOTAP(I,P,N,VET,NVl COMMON /GER~LIIR1IW,IT,NNF.,NGLN,NNOS .D I MENS I ON _P {_NGL.N t, V_E T (N_VJ_ l<t:NGLN•(NNOS•(IMl)+N-1) 00 10 K=l,NGLN . KK:Kt+K VF.T(KK)=P(I<) RETURN E:ND

129

SUBROUTINE SACACTEMP 1 N1,N2,VETl COMMON /GERALIIR,IW,IT,NNE,NGLN,NNOS,NELEM,NAP,NNED1LBNOS DlMENSION TEMP(NGLN,NGLN),VET(IT) Kt:::NGLN•(N!•I) K2=NGLN*(N2•Nl) IF(Nl•N2)!0111,12

10 DO 20 l=l,NGLN 00 20 J:1,NGLN lNO=NGLN•LBNOS•CK!+!•!)+K2tJ•J+!

20 TEMPCl,Jl:VfT(IND) RETURN

11 DO 30 I=1,NGLN DO 30 J=I,NGLN lNO:NGLN•LBNOS~(Kt+I•!)+J•I+I TEMP(I,J):VET(IND)

30 TEMP(J,Il=TEMP(I,J) RETURN .

12 WRITE(IW,100)N1,N2 100 FORMAT(/11QX!ERRO NA SUBRT, SACA NI= 'Il' N2= 1 Ill

CALL EXIT END

SUBROUTINE BOTA (TEMP,Nt,N2,VET) COMMON IGERALIIR,I~,IT1NNE,NGLN,NNOS,NELEM,NAP,NNED,LBNOS OIMENSlON TEMP(NGLN,NGLN),VETCIT) Kl=NGLN•(Nl•I) . K2:NGLN• (N2•NI) I F ( N 1 "N2) 1 O, 11 1 12

10 00 20 1=1,NGLN DO 20 J=!,NGLN INO=NGLN•LBN0St(Kl+I•l)tK2tJ•I•1

20 VETCIND)=TEMPCt,Jl RETURN

11 DO 30 1=1,NGLN 00 30 J=I,NGLN IND:NGLN•LBNOS•(Kt+I•l)+J•I+!

30 VET(lND):TEMPCl,J) RETURN

12 WRITE(IW,100)N1,N2 IDO FORMATC//IOX'ERRO NA SUBRT, BOTA NI= 'Il 1 N2= 'I3)

C4Ll. EXIT (NO

130

SUBROUTINE MONTA(MAT 1 NO,ED 1 MQUK) LOGlCAL BOOL1IMPRJM INJEGER TIPEL,CONEC,TlPAN,ESTAB REAL MAT,IZ,MESP,LZERO CQMMQN /GERALllR 1 IW,tT,NNE,NGLN,NNOS 1 NELEM,NAP1NNED1LBNOS,

*TIPAN,EST•B,INROT,MASSA COMMON /PROPIAX(lOO),MESPClOOl,CONECClo0,2),XCIOO),Y(IOOl,

tZ(!OOl,ECI00) 1 1Z(!OO),ITIPO(IOOl,LZERO(IOO),XZER0(100) 1 .

•YZERO(IOO),ZZERO(!OO),ANGZ(!OOJ,FZERO(!OO),RIGEL(IOO,b,b), •FORCA(IOO,b)1D2(3D0l 1 FEC30D) . . .

CoMMON /LOGICO/ BOOL,IMPRJM,JNT DIMENS{ON RIGE(b 1 b)1TEMP(1,3),N0(50)1EDCS01Q),MAT(lTl IFC MOU!< ,NE, 1 ) GOTO 111 . . IF( MASSA ,EQ, 3) GQ TO ~22

111 DO 20 NEL=l,NELEM TIPE~=lTlPOCNEL) GOTO (!,2) 1 MQUK

! GOTO (3 1 Q1 U,4 1 q),TIPEL 3 CALL MASSP(NEL,RlGE,lNROT,MASSA)

GOTO 7 q CALL MASSC(NEL,RIGE,MAS5Al

GOTO 7 2 GOTO CS,ó,& 1 &,b),TIPEL 5 CALL RIGEP(NEL,RXGE)

GOTO 7 b CALL RlSEC(NEL,RIGE) 7 DO JO I~l,NNE

N1=CONEC(NEL 1 I) DO 110 J=I,NNE N2:CONECCNEL 1 J) IF(N!•N2)11,12,12

11 Ntt::N2 N22=N! GO TO 13

12 Nl1=NI N?,2:N2

13 CALL SACA(TEMP 1 N22 1 N!! 1 MAT) IF!Nl•N2)15,lb,IQ

14 DO 50 KK:1,NGLN tK:NGLN~(I•l)+KK DO 50 L=t,NGLN. IL=NGLN~(J .. l)+L

50 TEMP(L,KKl=TEMPCb,KK)+RtGECIK,lLl GOTO 71 .

IS DO bQ KK=1,NGLN JKCNGLN*(I•l)+KK DO &O L=l,NGLN IL=NGLN, CJ•I) tL

bO TEMP(KK,L)=TEMPCKK,L)+RIGE(IK,IL) GOTO 71

lb DO 70 KK=t,NGLN

ll<=NGLN•(l•! l+KI< 00 7() L=KK,NGLN lL=NGLN•0.,1 l•L

131

70 TEMP(KK,L)=TEMP(KK,L)tRIGE(IK,ILJ 71 CALL 80T4(TEMP,N22,Nlt,MAT) -<lo CONTINUE .39 CONTINUE 20 CONTINUE 222 JF( NNEO ,LT, 1 ) RETURN

JJ:.2i>NGLN IF(MOUl<,EG,I) JJ=O DO 80 J;;\,NNED NOI=NO(l)"I 00 '10 Jr:t 1 NGLN INO=LBNOS•NGLN•(NGLN•NOitJ .. l)+l MAT(IND)=MAT(IND)+ED(I,J+JJ)

<IO CONTINUE 80 CONTINUE

IF( TIPAN ,NE, 3 ,OR, MOUK ,NE, 2 ,OR, ,NOT, BOOL l RETURN DO 82 t=l,NNED 11:NGLNi>NO(ll+I DO 82 J=!,NGLN

82 FECII•J):FECII•J)+D2(1J•Jl•ED(J,JtJJl RETURN . - - - -END

132

SUBROUTINI EXCITA(T,F,Df,NEQ) HHEGER CONEC REAL MESP DIMENSION F(NEQJ,OF(NEG) COMMON /PROP/AX(lOÕ),MESP(lOO),CONECCI0012)1X(IOO) COMMON /MDVEL/ PV,VI,A,NEP 1 NELP,xLOCAL,XEL- -XL;Xl;L NEL=NEL.P DO 10 I=l,NEQ F<JJ:o,

10 DF(I):o, !F(T ,NE, O,) GOTO 11 XL:X(CONECCl,2)) NEL':t

li E=XLOCALIXL E2=E•E E.3=E2*E V=VI+A*T I=CCONEC(NEL,ll•l)•3 J=CCONEC(NEL,2)•))*3 F(I+2);Py•Cl,•3,*E2+2,*E3) F(l+3):PV•~L•l~(!,•E)•Cl,•El F(J+2)=PV*(3,•E2•2,*E3l FCJ+3):PV•XL•E2•CE•!,) DF(lt2):PV•V•(•l,*E+b,*E2)/XL _DF (1+3)=pv~v·c 1,•'l. *E•3,*E2l DF(J+2l=•DFCI•2l DFCJ+3):PV•V•C3,•E2•2,•El RF.TURfll END

133

SUBROUTINE APOIO(MAT,VET,NEQ,NVET,RLOCAL,MCC,KOD/ REAL MAT COMMON /GERAL/JR,JW,IT,NNE,N.GLN.,NNOS,NELEM,NAP,NNED,LBNOS DIMENSION PC3l,P1(3) 1 C(3,3),CT(3,3),RC3,3l,RT(J,3),TEMP(3,3

•l,RL0CAL(20,q),MCC(20,S),REMP(3,3)~SEMP(3,3l,VET(NEQ), •MAT(ITl . .

DO 10 MM=t,NAP N.llMCÇ(MM,!) DO 20 I:;t,NGLN DO 30 J=t,NGLN C<I,J):io,

30 RCI,J):O, C(J,Il=t•MCC(MM,l+l)

20 R(l,I):1, lF(MCC(MM,NGLN.+2))1,211

1 CALL ROTCR,RLOCA~,CX,CY,CZ,MM) CALL MULTMCC,C,R,NGLN,NGLN,NGLNl

2 GO TO(U2,21,21),KOD 21 00 22 I=l,NGLN

DO 22 J=1,NGLN RTCJ,I):R(I,Jl

22 CT(J,I):C(I,J) IFIM=LBN0S•1 DO 40 1=1,r•IM J::N•I IFCJ)S,S,3

3 C1LL ~iCA~TEMP 1 J 1 N,MATl CALL MULTM($EMP,TEMP1CT,NGLN,NGLN,N.GLN) CALL BOTA(SEMP,J,N,MAT)

5 J::IHI lF(NIIQ5 .. Jl40,b,b

b CALL S~CACTEMP,N,J,MAT) CALL MULTM(SEMP,C,TEMP,NGLN,NGLN,NGLNl C~LL BOTA(SEMP,N,J,MATl

40 CONTINUE CALL 5ACA(TEMP 1 N1 N1 MATl GO T0(42,42,43),KOD

42 00 120 KK:1,NVET CALL SACAP(KK,P,N,VET,NEQl CALL MULTV(Pt,C,P,NGLN,NGLN,O) c•LL BOTAP(KK,Pt,N,VET,NEQ)

120 CONTINUE ~l GO TO(l0,121,121),KOD 121 C•LL MULTM(S~MP,TEMP 1 CT,NGLN 1 NGLN,NGLN)

CALL MULTM(TEMP,C,SEMP,NGLN,NGLN,NGLN) DO lóO LL~l,NGLN

lbO TEMPCLL,LLl=TEMP(LL,LLl•MCCCHM,LL•tl CALL BOTACTEMP,N.,N,MAT)

10 CONTINUE GO TO(lB1,11,!1),KOD

11 DO 180 I•l,NNOS

134

IND:LBNOS•NGI.N•CNGLN~(I•!)+2)+1 IF(MAT(IND)ll80;)70,l80 -

170 MAT(INOJ:1. 180 CONTINUE 181 RE1URN

END

SUBROUTINE TRIANGCA,NEA,NEQ) DIMHJSION A CNEA) LB:NEA/NEQ OQ 30 N=l,NEQ l=N DO 29 L:.2,LB I=I+I 1 NO:;L + Ci',• 1) *LB IF(A(IN0))24 1 2q,24

24 c:A(IND)/A(l+INO•L) J=o 00 27 !<=l.,LB J:J+1 H,O}:;~+CN•l l•LB IF(A(lN01))2b,27,2&

2& lND2=J+(I•1)~L6 A(IND2):A(IND2)~C•A(IND1)

27 CONTINUE __ _A_( I NO l =e

2q CONTINUE 30 CONTINUE

RETURN END

135

SUBROUTINE RES08(A 1 NEA,B,NEB 1 NVET) DIMENSION A(NEA),B(NEB) NEQ:NEB/NVET . L8:fJEA/NEG DO 30 N=l,NEQ I=N DO 29 L=2,L8 I=I+I IND:.L• (N"'l ),>l,.B IF(A(INDl)28,2Q,28

28 00 281 IV=l,NVET INC=(IV•l)i>NEQ

281 BCI•INC):B(I+INC)•A(JND)•B(NtINC) 2Q CONTINUE

00 291 IV=l,N\iET JNDV:N+(IV•!l*NEG

2QI B(lNDV):B(INDV)IA(l+(N•ll•LB) 30 CONTINUE -

N=NEQ 35 NCN'"l

1FíN)50,50,lé 3é L=N

Do 40 K=2,LS L=L+l IND=K+(N•ll•LB IF(A(XND))31,ao,37

37 DO 371 IV=t,NVET INC=(Iv•tl•NEQ

371 B(N+JNC):B(N+INC)•A(JNO)tB(L+INC) 40 CONTINUE .

GOTO 35 50 RETURN

ENO

SUBRDUTINE RES00(A 1 NEA 1 B1 NEB 1 NVET) DIMENSION A(NEA) 1 B(NE8) NEG=Nf.8/NVET . LB=NEA/NEG DO 20 IV=l,NVET IVET:(JV•ll•NEQ DO lO t=1,NEf.l

10 B(I+IVtT)~BCI+IViTl/A((I•ll~LB+ll 20 CONTINUE -

RETURN ENO

136

SUBROUTINE GLOBAL(P,NP,MCC,RLOCALl CQMMQN /GERALltR 1 IW,lT 1 NNE 1 NGLN,NN05,NELEM,NAP DIMENSION P(NP),RLOCALC20,q),RC3,3),U1Cl),U2(l),MCC(20,5l NVET=NP/(NGLN•NNOS) - . - . DO 10 !=1,NAP . If(MCCCI,5))1,10,1

1 CALL ROT(R1RLOCAL,CX,CY,CZ,I) DO 20 J=i,NVET c•LL SACAP(J,U),MCCCI,1)1P,NP) CALL MULTVCU2,R,U!,NGLN,NGLN,1)

2D CALL 80lAP(J 1 U2 1 MCC(I,ll,P,NP) -10 CONTINUE

9 1 O

1 1

12 19 2Q 30

RETURN ENO

SUBROUTJNE MULTB(A 1 V,P 1 NA,NV) OIMENSION A(NA),V(NV),PCNV) LB=N•.INV 00 30 I=l,NV P(I)cO, DO 20 J=l,NV J8=J•I+1 IF(JB•LB)q,q,30 IFCT•Jl 10, 10, 11 INOcJBt<l•l )•LB GOTO 19 IB=l•Jtl IFíJB•LBl12,12,20 IND=IB+(J•I )•l.a PCI):P(Il+A(lNOl~V(J) CONTINUE . CONTINUE RETURN END

137

SUBROUiINE CONVCVl,V2 1 N,TOL,ICONVl REAL NORMl 1 NORMD COMMON /GERALIIR1IW OIMENSION Vl(Nl,V2(N) lCONV:1 NORMl=O, NORMD=O, DO 10 I=l,N NORMt=NORM!tlil (I)•VI (I)

10 NORMD=NORMO+(V2(I)~V1(I)l**2 NORMl=SQRTCNORM!) NORMO:SQRT(NORMD) IF(NORM!•l,E•l5)101,101,102

101 WRITE(IW,100) 100 FORMAT(/ 1 10X1'NORM- DO VETOR VI MENOR QUE 1,E•IS'.l

NCRMl=l, 102 IF ( CNORMO/NORMI )•TO!.) 12, 12, 11 11 lCQNy:z 12 RETURfl

END

1

10

2

20

138

SUBROUTINE MASSP(NEL,RTSR,INROT,MASSA) tNTEGER CONEC REAL MESP,tZ COMMON /PROP/AX(IOO),HESP(IOO),CONECCI00,2),X(IOO),YCIOO),

•Z(!OO),E(!OOJ,I2(100l,ITIP0Cl00) . . . DIMENSJON RIGECb,b),RTSR(b,b),R(3,ll1RLOCAL(201q) DX:X(CONEC(NEL,2)l~X(CONEC(NEL,ll) -DY=Y(CONEC(NEL,2)) 8 YCCONECCNEL,1)) XL=SQRT(DX*DX+DY•DYl -CX::oDX/XL CY=D\'IXL f:MESP(NELJ•AX(NEL)•XL G:MESP(NELJ*IZ(NEL)*INROT GOTO (1,2ill,MAS5l RlGE ( 11 I ):F/3, RJGE(l,2)=0, RIGEC1,3):0, RIGE ( 1,,q:F/ó, RIGEC1 1 5)=0, RIGECl,&):O, RIGE(2 1 2)=F•l3,l35, R!GE(2 1 3l=F•!!,1210,*XL RIGE(2,4)::0, RIGEC2,Sl=F•9,/70, RIGE(Z,b)C•F*l3,/q20,•XL

+ +

.. +

R J GE O ,3) =F '1 os. otX_L~XL ___ + RIGEC3,ll):0, RIGE(3 1 5):•RIGE(2,ó) RIGE(3,b)=•F/IQO,*XL•XL RIGE(Q,4):RJGEC!,1) RIGE (!l 1 5)i=O, RJGE(ll,E>)=O, RIGE(5,5):RIGEC212l RIGECS,b)=nRIGE(2,3l RIGE(E> 1 b)=RJGEC3,3) 00 10 1=1,5 . II=ltt oo 10 J=II,b RIGE(J,Y)=RIGECI,J)

..

CALL RoT(R,RLOCAL,éX,CY,o,o,o) C~LL ROTM(RTSR,RlGE,R) RETURN DO 20 1=1,b DO 20 J=t,b RTSRCI,J):O, RTSR(l,t)=F•0,5 Rr5R(2,2)=RrsR(1,1) RTSR(3 1 3):F•XL•XL/4?.0, RTSR(ll 1 1l):RTSRC!,1) RTsRCS,Sl=RrsRC!,1) RTSRCb,b):RTSRC3,3)

t

G*l,2/~L G•o, 1

c;•t,2/XL G*O,I G•2,115,,.XL

139

3 RETURN END

140

SUSROUTJNE MASSC(NEL,RTSR,MASSAl INTF.GER CONEC,TIPEL REM. MESP COMMON /PROP/AX(I00) 1 MESP(100l,CONEC(I00,2)1X(100),YCIOO),

•ZC!OO),E(IOO),IZ(IOO),ITIPO(lOOl . . . DIMENSION RIGE(ó,ól,RTSR(&,&),R(J,3),RLOC•LC20,9l TIPEL=ITIPO(NEL) . DX::X(CONEC(NEL 1 2))•X(CONEC(NEL,1)) DY:Y(CONECCNEL 1 2J)•Y(CONEC(NEL,l)l DZ=Z(CONEC(NEL12)l•ZCCONEC(NEL,1ll IF( TIPEL ,LT, a) DZ=o. . XL::SQRT(DX*DX•DY•DY+DZ•DZ) CX::PX/)(L CY=DYl)(L CZ;oDZ/XL GOTO (111,222 1 333),MASSA

Ili f::MfSP(NEL)*AX(NEL)•XL/3, RJGE(1 1 1):F -RIGE(1 1 2l=O, RIGE(1 1 3l=o, RIGE (1,4l=F*0,5 RIGE0,5)=0, RIGE(1,b)::o, RJGEC2,2)::F RIGE(2 1 3)=0,

_RIGEC_2_ 1 11):o, RIGE(2,5)=F•o.s RIGEC2 1 b)::O, IU GE (3, 3) =O, RJGEC.J,u)=o, RIGE(3,S)i;O, RIGEt3,E>)=O, ?IGE(4,l.!)=F RJGECll,S)cO, RIGE(ü,b)=O, RIGE(S,S)=F RlGE (5,&);'•0, RIGECE>,b)=O, Go roC70,2,2,1,1l,rtPEL

1 RIGE(3,3)= f RIGE(3 1 6)= F•0,5 RIGE(b 1 t,):; F 00 lO I=\,5 II=I+I DO 10 J=II,E>

10 RIGE(J,ll=RIGECI,J) CAL~ ROT(R,RLOCtL,éX,CY,CZ,O) CALL ROTMCRTSR,RIGE,Rl RETURN

222 00 20 Ict,b 00 20 J=l,t>

20

30

333 70

141

RTSRCI,J):o, RTSR(l,1):MESPCNEL)*AX(NEL)*XL*0,5 RTSR(2 1 2)=RTSRCl1l) . RTSR(~,4):RTSR(t,!i RTSR(S,S)=RTSR(l,!) GOTO (70,333,333,iD,3Dl,TIPEL RTSR(313):RTSR(l,1) RTSR(&,6):RJSR(t,ll RETURN . C~LL EXIT END

142

SUBROUTINE RIGEP(NEL,RTSR) LOGtCAL 600L,1MPRIM INTEGER CONEC,TIPAN REAL IZ,MESP,LZERO COMMON /GERALIIR,I~,IT,NNE,NGLN1NNOS,NELEM,NAP1NNEO,LBNOS,

•TJPAN COMMON /?RQP/AX(!OO),MESPC)OO),CONECCI00,2),X(IOOl,Y(tOO),

•ZCIOOl,E(IOOl,IZCIOOl,ITIPO(!OOl1LZERO(IOO),XZERO(lDOl, -•YZERO(IDO),ZZERO(IOO),ANGZ(!OOl,FZEROC!OO),RiGEL(lOO,b,b), •FORCA(IOO,&),D2(300) 1 FE(300) . -

COMMON !LOGICO/BOOL,IMPRIM,!NT OlMENSION RIGE(b,&l,RTSR(b,b),R(3,3l,RLOCAL(201ql,

~ F!Cl),F213l,V1(3l,V2(3l,FEL(&l,VGCbl E Q UI V!< L EN CE C FEL ( \ ) , F ! Cl )') , ( F E(. ( 4 l Í F 2 ( 1 J )

* (\IGCi),V!Clll1CVG(tll,V2(!)) OX=~(CONEC(Nf,L,2))~xccoNEC(NEL,t)) OY:V(CONEC(N~L,2))•Y(CONEC(NF-L,1)) XL=SQRT(O~•DX+PY•DY) -CX=DX/XL CV:OY/XL CA~L ROT(R,RLOCAL,CX,CY,o,o,ol EA=E(NELl•••CNELl EI:F.(NEL)*IZ(NEL) NEGl=NGLN•NNOS GO TO(t,2,2),TtPAN

1 F:FZERO(NELl GOTO 3

2 f:(~L-LZERO(NEL)l/LZEROCNEL)tEA 3 RrGE(l11l=EA/XL

RIGE(1,2)=o, RI GF. ( l 1 3) "O• RIGE(l,ll):~RIGE(l,I) RIGE(l,S)=O, RíGE(l,bl=o, RIGE(2,2)=12,*EI/(XL**3) + F•l,2/XL RIGE(2,3)~ó,•El/(XL•XL) • FtO,I RH,Eç2,ai::o, RIGE(2 1 5l=~RIGE(2,2) RIGE(2,b)=RIGE(2,3l RIGEC3,3)=Q,*EI1XL + F•2,l15,*XL RrGE(l,ll)=o, RIGE(l 1 5):eRIGE(2,3) RIGEC3 1 ó)"2,*EI1XL - P/30,*XL RIGE(a,a)=RIGECl,1) RIGE(ll,S)=o, RtGE (11,ó)"O, RJGE(S,5)=RIGE(2,2) RIGECS1b):RJGF.(3iSi RlGECó,b)=RIGE(l,3l ºº 1 O l = 1, 5 -II=I+I

143

DO 10 J:II,& 10 RIGECJ,I)"RIGE(I,Jl

GOTO (23,21,22),TlPAN 2\ !F(.NOT, tMPRIM) GOTO 1.1

22 CALL SACAP(l,F!,CONEC(NEL,1),02,NEQ) CALL SACAPCl1F2 1 CONECCNEL,2l 1 D2,NEQ) CALL MULTV(V!,R,fl,3 1 31 0) CALL MULTVCV2,R,f2,3,3,0) F'EL(tl=•F FEL(2):(VIC2)•V2C2))*12,~EI/(XL*~3)+(V1(3J+V2(l))*&,•EI/

• . . . . (XL~XL) . FELC3l=(V!C2)•V2(21)•b,•EJ/CXL•XLl+CV2(3)+2,•VIC31)~2,*EI/

* XL . l'EL(ll)i: F FEL(S):"rEl.(2) FEL(6):(y1(2l•v2(2)l•b,•EI/(XL*XL)+(vl(ll+2,*V2(J)l*2,•EI/

* XL 00 bO I=l,6

=O FORCA(NEL,Il=FEL(Il GOTO 4 .

23 00 1.10 I=l,b 00 1.10 J:1,=

ao RIGEL(NEL,I,J)=RIGECI,J) 4 CALL RQTM(RTSR 1 RIGE,RJ

IF( lNT ,EO, O ) GOTO 5 IF( TlPAN ,NE, 3 ,OR, ,NOT, BOOL) RE_lURN

5 ·t~L~·MULT~(Vl,R,Fl,3,J,11 CALL MULTV(V2,R,F2,l 1 31 !) CALL SACAP(l1F!1CONEC(NEL,1),FE,NEG) CALL SACAPCl,F2,CONEC(NEL 1 2l,FE,NE~) oo 30 r=1,1i

30 FEL(I):FELCl)+YGCI) CALL BOTAPC!,Fl,CO~EC(NEL,1),FE,NEQ) CALL BoraPC!,f2,CONEC(NEL,a),FE,NEG) RETURN . . . END

50

1

ºº

bO li

144

SUBROUTINE RIGECCNEL 1 RTSR) LOGICAL BOOL lNTEGER CONEC,TIPAN,TIPEI. REAL MESP,IZ,LZERO COMMON /GERALIIR,IN,IT,NNE,NGLN,NNOS,NELEM,NAP,NNED,LBNOS,

•TIPAN COMMON /pROP/AX(t00) 1 MESP(!OO),CONEC(I00,2),X(IOO),v(IOO),

~Z(\00) 1E(!OO),IZ(!OO),ITIP0(!00)1LZER0(100),XZERO(!OO), . *YZERO(IOOl,ZZERO(IOO),ANGZ(lOOl,fZERO(!OOl1RiGEL(l001b,ól, •FQRCA(!00 1 ó),D2(300),FEC300) . .

COMMON /LOGitO/BOOL,IMPRIM,INT DlMENSIDN RIGE(b,ól,RTSR(b,b),RC3,3l,RLOCAL(201ql,

• Ft(3l,F213),Vl(3),v2(3),ffL(ól,vG(b) . EGUIVALENCE (FELC!l,Fl(ll'l,CFEL('ll,F2(t)l

~ (VG(i),V!C!)l 1 (VG(4),V2(!l) TIPEL=ITIPOCNEL.) DX:X(CONEC(NEL,2))-X(CONECCNEL,1l) DY:Y(CONEC(NEL,2l)•Y(CONEC(NEL,l)l oz=Z(CONEC(NEL,2))•Z(CONEC(NEL,l)l IF( TIPEL ,LT, u) OZ:O, XL:SQRT(DX•DX+OY•DY+Dz,oz) CX=OXIXL CY=OY/XL CZ:DZ/XL. DO 50 I=t•ó FEL(I):O, FORCtdNEL., I ):O, 00 50 J:1,ó R J Gf;'. ( I , .J l : O , GOTO (2,1,1),TIPAt; F:(XL•L.ZEROCNEL))/LZERO(NELl*E(NELJ•AX(NELl/XL !FC F ,GT, 01 ) GOTO 3 . . . GOTO C70,3,20,3,2ol,TIPEL F:FZEROCNELl/XI,. RIGE(l1ll=ECNEL)*AXCNELl/XL RXGE(l 1 4)=•RIGE(l,1l RIGE(4 1 !):RIGE(l,4) RtGE(a,a)=RIGECl,t) RlGE(2 1 2)=F . RJGE(2,5l=~F RIGE(S,2J: .. F RJGE(S,Sl=F GOTO (70,óO,ó0,40,40),TIPEL RIGEC3,3)= F RIGE (3 1 bl=•F RIGE(b,3):•F RIGECb,b)= F Go TO cs,o,~),T!PAN U2=CXL•LZEROCNELl)/LZEROCNELJ•XL f0RCA(NEL,ll=U2~RIGE(t,aJ FoRCA(NEL,a)=•FoRCA(NEL,ll

FEL(tl:fORCA(NEL,ll fEL(al=•FELCI) . GOTO b

5 DO 10 !=1,b DO 10 J=l,b

145

10 RlGELCNEL,I,JlcRIGE(I,Jl b CILL ROT(R,RLOCAL,CX,c~.cz,oJ

CALL ROTM(RT$R,RIGE,R) IF( INT 1 EQ, O l GOTO 7 IF( rIPAN ,NE, 3 ,OR, ,NOT, BOOLl RETURN

7 NEQ=NGLN*NNOS CALL MULTV(V! 1 R1 F!,3 1 J,11 CALL MULTVCV2 1 R,F2 1 3,J 1 1) CALL S•CAP(l,f2,CONEC(NEL,2l,fE,NEQ) CALL SACAP(!,F!,CONEC(NEL,l);FE,NEQ) DO 30 I=t,b . -

3Q FEL(I)=FELCI)+yGCil CILL BOTAP(t,F1,co~EC(NEL,ll,FE,NEG) CALL BOTAPCl,F21CONEC(NEL,2J,FE,NE~)

20 RETURN -70 CALL EXIT

ENO

SUBROUTINE SAIDA(T,lTERAD,ITERA1D2 1 V2 1 A2 1 NEQ) JNTEGER ESTAB,TIPAN,ClRREG COMMON /GERAL/IR,l~,IT,NNE 1 NGLN1NNOS,NELEM 1 NAP,NNE0 1 LBNOS 1

*TIPMJ,ESTAB, INROT ,MASSA; IA MOR, AL}~,BETA, TOL!, TOL.2, TF, TAU, •IMPRES,CARREG

COMMON /IMPRE/IMPSl,IMPDVA(IOO) OIMENSION D2(NEQJ,V2(NEG),A2(NEQ) GOTO (1,2),CARREG

l WRITE(I~,lOO)T,ITERAD,lTERA GOTO 3

2 TR::2.•T/TF WRITE(JW 1 l00)TR,ITERAD 1 ITERA

100 FOFIMI\T(/IOXITEMPO t,El5,7,5~,!ITERAC:0ES !I2,• , '121) 3 IF(IMPSI ,EQ, O) RETURN

WR!Tf(IW,200) 200 foRMAT(l/,JX, 'N0',1ox,1011,11x, 1v1•,11x,141 ',t3X, 102t,11x,

••v2•,11x, 'A2',13~, IOJ',tlX, •v31,11x, IA3•,,1 -DO 10 J=l,NNOS -1=31tJ~2 IF(IMPDVA(J),EQ,1) WRITE(IW,300) J,D2CI),V2(I)1A2(I),

* . D2CI+ll,V2(Itl),A2(It1l,D2(lt2),V2CI+2),A2<I+2) 300 FoRMÀTCI4,~X,3E!3,S,2X,3~13,S,2X;3E13,S) . 10 CONTINUE .

RE TlJRN END

êOO

30

20 1.10

300 1 O

SU8ROUTINE ESFLIN INTEGER CONEC REAL MESP,1Z,L2ERO

146

COMMQN /GERAL/lR,IW,IT,NNE,NGLN,NNOS,NELEM COMMDN /PROP/AX(IOOl,MESP(lDO),CONEÇ(I00,2),X(tOO),YCIOO),

•Z(IOOl,E(I00)1lZ(IOO),lTIPO(IOO),LZERD(IOO),XZERO(lOO), . •VZERO(lOOl,ZZEROC!OO),ANGZ(lOD),~ZERO(tOO),RIGEL(lOO,b1b)1 •FORCA(1DO,b),D2(300J .

COMMON IIMPRE/IMPS1,IMPDVA(IOO),JMPS2,IMPFOR(100) DIMENSION RIGE(b 1 b),Vt(3),V2(3),V3(3l,V4(3),V(ó)íR(3,3),

• F(bJ,UCó),RLOCAL(2Q,q) -EQUIVALENCE (V3(1),V(l)l,(V4(1),V(4l),(Vl(1),U(!ll,

•(V2(!) 1 UCU)) - -NEQ:NGLN~NNOS ' wRITE(lW,200) FoRMAT(/l,2X, 'EI..EM,

*'ESF, 3 1 ,8~,'NO 21 NO li ESF, 1',12X 1 'ESF, 2 1

1 12X, ESF, 1',!2X,'ESF, 2•;12x,,ESf, 31,n

DO 10 I=!,NELEM tF(IMPFOR(ll,NE,I) GQ TO 10

- DO 30 J= 1, b DO 30 K=l,ó RIGE(J,K)cRJGEL(I,J,K) C~LL S~CAPCl1Y!,CONEC:(I,1),D2,NEQ) 'ALL SACAP(l,Y2,CONEC(I,21,D2,NEQ) IF(aNGZ(J) ,EQ, O, ,AND, ITIPQ(I) ,LT, Ul GOTO 20 CX=JXZE:RO{CONEC U,2J )'"XZJBO{CONEC: ( I, U) )_11._lEBO_UJ. ________ _ cv=(YZERO(CONEC(I,2JJ~yzERO(CONEC(I,1)))/LZERO(l) CZ=(ZZERQ(CONEC(I,2l)•ZZERQ(CONEC(I,1)))/LZEROiI) IF( ITIPO(I) ,LT, li) CZ:O, , . CALL ROT(R,RLOCAL,CX,CY,CZ,O) CALL MULTV(V],R,v1,NGLN,NGLN,o) CALL MULTV(Vl.l1R,V21NGLN,NGLN10) CAL\.. MULTVCF,RIGE,V,b,6,0l GOTO 110 CALL MULTVCF1RIGE,U,b,b1D) F(!l=F(t)eFZERO(I) F(u):F(aJ+FZe:Ro(I) WRITECIW,300) I,(FCJ),J:1,b) FORMAT(III,' ;,,, ',3E!6,7,UX 1

1 ~ 1 ,3E16,7l CONTINUE RETURN ENO

147

LOGICAL BOOL,IMPRIM INTEGER CONEC1TIPAN,ESTAB,CARREG REAL M(IOOOO),K(IOOOO),MESP,IZ,Ml,M2,KMOLA,LZERO COMMON /GERALIIR 1 Iw,IT 1 NNE 1 NGLN,NNOS,NELEM,NAP1NNED,LBNOS,

•TIPAN,ESTAB,INROT,MASSA,IAMOR,ALFA,BETA,TOLl,TOL2,TF,TAU1 •IMPRES,CARREG

COMMON /PROP/AXCIOOJ,MESP(lDO),CONECCI0012),X(I00)1YCI00)1 •ZC!OOJ,E(IOO),IZ(\00),ITIPO(lOOl,LZERO(lOO),IZEROílObl, -*VZERO(t00)1ZZER0(!00)1ANGZC!OO),FZEROC100),RIGEL(!OO,ó,bl1 *FORC~(IOO,&l,D2C300l,FEC300l . . -

COMMON ILOGICO/ 600L 1 IMPR1M,INT COMMON IIMPRE/JMPS1,IMPOVA(IOOl,IMPS2,IMPFOR(IOO) COMMON /MOVEL/ PV,Vl 1 A,NEP,NELP,XLOCAL,MEL . DIMENSION ED(50,q),MCC(20,5J,RLOCALC20,q),N0(50),CADC8l,

• C(IOOOO),P(ÇQO),FC300),DFC300l,D1Cl00)1V1(300), * FI1(300),DFJ!(300),FI2(300),DFI2CJOO),DDCJOO), * DDZC300) 1 Y2(JOO),A2(300l1FÍ22(300),DPI22C3QOJ, * AA(l00),6B(300),CCC300l,AA1(300),BBl(3ool,CCi ~ (300) . . . .

DATA CAD/21,,3,,q,,2,,20.,2.s,1c,,2,s1 TEMPO=TIMEC2) tR:6 I \oi 115 1 TMA;r:=so NNE:;2 NGLti::3 CALL DAOOS(NO,ED,MCC,RLOCAL) Cl=H,U•O,S C2=Tt,U~TAU/12, C3:TAU~TAU/&O,•CAD(4~ESTAB+I) C4=TAU•TAUlóO,•TAU•CAo(a•ESTABt2) C5:TAU•TAU/&O,tCAD(Q~ESTAB+3) Có=TAU•TAU/óO,tTAU•CADCa•ESTAB+4) INCL1=0 -DO 5 1:1,NilP IF(MCC(I,S))ó,S,ó

5 CONTINUE GOTO 7

ó INCLI=I 7 NEQ:NGLN*NNOS

NC:NGLN•LBNOS IT:NEQ•NC DC 10 l=ldT MCI)=o, C(1):o,

1 O K C I l :o, 800Li=.FALSE, 1MPRP1=, TRUE, INT:o CALL MONTA(~,NO,ED,2) IMPRIM::,FALSE,

148

CALL MONTA(M,NO,ED,ll WRITE(IW,3qoz)CFECI),I=t,NEQ)

3902 FORMAT(/lOX,tFORCAS ELASTICAI EM T:0',1,C!OX,3El5,7)) GOTO (103,101,102),IAMOR

!01 DO 11 1=1,IT . 11 C(!)=ALFA•M(I)+SETA*K(Il

IF(NNED,LT,I) GOTO !03 DO 110 I=t,NNED NOI=NO(l)•1 DO 110 J:),NGI.N IND=LBNOS•NGLN*CNGLN•NOJ+J•ll+l

· 110 C(IND)=CCINDl+EDCJ,JtNGLN) GOTO 103

102 READ(IR,400)(C!Il,I=t,IT)' aoo FORMAT(BFIO,Ol 103 f:O,

GOTO (10ó,IQ4),CARREG

104 COMP:X(CONECCNEP 1 2)) NELP:o )(o::o, )(f:L=O,

c~~·~••••e•~~w·~~••~•~~~-~•-•~~~-·~~--ft~•~••••--• CARG~ MOVEL •~•• ! o t\, Q

30

100

óOO

'1 o o

31 32

33 31.1

"º 1 O O O

C~LL EXCITA(T,F,Df,NF.Q) DO 30 I=1,NEQ AA1(l)::O, BBI CI ):O, Dl(I):o. V!(I):O, FUC!l=•CIJ CALL APOIO(M,F,NEG,1,RLOCAL,MCC,2l WRITE(IW,IOO)(H(I);I:t,ITJ . FORMAT(tOX,IM .ATRIZ D E MAS S A',l/,(bElS,7,/ll {F( IAMOR ,NE, 1 l WRITE(IW 1 bQO)(C(IJ,I:1,IT) . . fORMAT(//,IOX,'M ATRIZ D f. AMOR T-E C IM EN Tt

*• ! O' ,li, (bEIS,7 1 I)) WRITE(Iw,qoo)(~Ct),I=t,IT) FORM/\T(//• IOX, IM A T R I Z D E Ri G IDE 2 1 ,//• , . .

•CbEIS,7,/)) IF(MASSA ,GT, () Go TO 31 CALL TRIANG(M,?T,NEQJ CALL RESOB(M,IT1F,NEQ,I) GOTO 32 CALL RESODCM1IT,F,NfQ,ll IF(INCLJ ,GT, O) CALL GLOBAL(F,NEQ,MCC,RLOCAL) GOTO (34,33,33),IAMOR CALL MULTBCC,F,AA,lT,NEQ) Do ao t=t,NEQ DFI1CI)=DF(l)•AA(I) INT=INTH T=HiT*TAU

ITERA:O BOOL:,FALSE, IMPRIM:.fALSE,

149

lf(INT/lMPRES•IMPRES ,EQ, INT) IMPRIM=,TRUE, C•••~--~P~••••••~•AGUI SE INSERE O ~CHEC~"~•w•~•••••••••~"·•~•~••

DO 13 I=1,NEQ DF I 2 (I l =DF I 1 (I) FJ2(J):FJ1(I)+TAU•DFI1(!) AA(I)=Ct•FI1CI)+C2•DFI!(1)

13 68(I)=C3•Fil(l)+Çü•DFil(I) GO TO (t14,131) 1 CM~REG

c~~-~•~~-~~~~~~--~~-~~ea••~~••••~~•8••••••·~-•-~ CARGA MOVEL ~-•~ 13! s=vt•T+O,S•A•T•T

XLOCAL=S•XO IF(S•COMP)43,43,41

41 DO U2 I=t,N(Q F<t>=o.

112 Of(l):O, GOTO 2000

113 IF(XLOCAL,LE,XEL) GOTO 4ü NEl.P=NELP+l IF( NELP ,GT, NEP l GOTO UI xo:xo+xEL XLOCAL=s .. xo XEL:X(CONEC(NELP,2ll•XCCONEC(NELP,1)) WRITECIW 1 60l)T 1 S1 XEL,NELP

bOl FORMAT(/,!8X,!Tl,1qx,1s1,1ex,,xEL•,12x,tNELP',/l1!0X, *3(EIS.7,5X),15 1 /) c•~P-~~-~-~~---~~~w,~e~-~~---~~~-·~e~---~----··· CARGA MOVEL ••••

aa CALL EXCITACT,F,OF,NEQ) 2000 ITERA=ITERA+1

DO 14 I=l,NEQ P(ll:FI2(I) P(ltNEQJ:AA(I)+Cl•Ft2(1)•C2•DFl2(Il

IU P(l+2*NEQ):BB(IltC5*FI2(I)•Cb*DFI2(1l CALL APOIO(M,P,l•NEQ,3,RLOCAL,MCC,l) IF(MASSA ,GT, 1) GOTO 15 C•LL RES08(M 1 IT,P,3•NEG,3) GO TO 1 b

15 CALL RESOD(M,IJ,P,3*NEQ,3l lb tF(INCLI ,GT, O) CALL GLoBAL(P,3•NtG,McC,RLOCAL)

DO 17 I=1,NEQ DD(I):P(l+2•NEQ)+TAU•~11I) D2CI)=D1(Il+DD(Il V2(I)=Vl(ll+P(I+NEQl

17 A2(I):P(Il 200 FORMATl/,bf-15,7)

GOTO (172,\71,171l,IAMOR 171 C•LL MULT8CC,A2,AAl,IT,NEG)

CALL MULTB(c,v2,BB1,IT,NEQ) 172 C~LL MULTB(K,V2,CC,IT,NEQ)

18

202

- 29 li

207 208

210

20

24

2.Q 1 300

25 2501

500

251 252

1S91 ;?b 27

28

150

CALL MULTBCK,DD,CCl,IT,NEGl DO 18 I:::1,NEQ DFI22(Il=DF(I)•CC(I)MAA1(I1 Ft22(Il=FCil-FE(I)•CC\(tl~B6!(Il IF(BOOL.OR,TIPAN,EG.1) GOTO 20 IPC ITERA .LE, 1 l GOTO 208 CALL CONV(DDZ,OD,NEQ,TOLl,ICONV) GOTO (202,2Q7),ICONV 00 201.1 t=!,NNOS X(I):::XCI)+00(3*I•2) Y(I):Y(.I)•OD(3~T•I) ZCl)=Z(I)+ODl3~I) -lTER~D=ITERA CALl MONTA(K1NO,E0 1 2J DO 201:i I=t,?T K(I)::o,5~1<:CI) BOOL.=, TRUE, GO TO 172

/

lf(ITERA,GT,ITMAX) GOTO 241 DO 210 I=l,NEQ DFI2(I):DfI22CI) FX2CI)=FI22(Il DDZ C t) :DO C I) GOTO 2000 CALl CONVCDF12,DFI22,NEQ,TOL2,1CONV) G o T O C 2 5 , 2 1 l , I C Q.N V IF( ITER~ ,GE, ITMAX ) GOTO 241 DO 24 l=l,NEQ 0FI2(I)=OFI22(I) FI2(I):Fl22Cl) GOTO 2000 WRITE(JW,300)!TMAX FORMAT(// 1 10X'NAO HOUVE CONVERGENCIA APOS 1 131 ITERACOES,t

_, EXECUCAO INTERROMPIDA!///,IOX'VALORES FINAIS!!) GOTO 2501 lFC,NOT, IMPRIM) GOTO 251 C~Ll SAI0A(T,ITERAD,ITERA,D2,V2,A2,NEQ) IFCIMPRIM ,ANO, TIPAN ,EG, 1 ,ANO, IMPS2 ,NE, O)CALL ESFLIN IF(,NOT,BOOL,AND,TIPAN,GT,I) WRITE(IW,500) FORMAT(//,!OXI NAO HOUVE CONV~RGENCIA NOS DESLOCAMENTOS,') GOTO (251,2é),lCONV . IF(T•TFl27 1 252 1 252 TEMPO=(TIMEC2l•TEMPO)/éO, WRITE OW, 15'11) TEMPO FORMATt//,IOX,ITEMPO DE PROCESSAMENTO 1',E12,5/ c~u. EXlT DO 28 1=1,NEQ DFII (I)=oFI22(Il Fl!Ol=FI22(!) D1(I):02(I) V! C!)=V2(Il

52

54

55 70 90

800

700 óO

GOTO (2q,zq,so),TIPAN DO 291 I1:t ,Nf;.Q FE(lJ=FECI)•CC1Ctl

151

GOTO (1000,50,SOJ,TIPAN DO 52 I=l,NNOS X(t)=XZERO(I)+P2(3*I•2l Y(l)=YZERO(Il•D2(3*I~t) Z(I)=ZZERO(Il+02(3•Il DO 54 I=!,IT l<(I):O, GOTO (1000,90 1 55),Tl?AN DO 70 1"'1,NEGI FECI)cO, C4LL MONTA(l<,NO,ED,2) IF(,NOT, JMPRJM ,OR, tMPS2 ,EG, Ol GQ TO 1000 WRITE(IW,8QO) FORMAT C/ I, 2X, 1 EL.EM,

•lf:;SF, 3',SX,'NO 21 DO óO I:l,NELEM

NO 11 ESF', \!,12x,•e:sir, 2'.t2X, ESF, 11,12x,'ESF, 2•;12x,1e:sF, 3 1 ,n

IF(JMPFOR(I) ,EQ, ll WRITECiw,700)1,CFQRCA(l,Jl,J:1,bl FORMAT(I4,! "" 1,3El8.7,IIX,'*' 1 3E18,7) CONTINUE GO TO 1000 END