View
2
Download
0
Category
Preview:
Citation preview
Universidade de Aveiro
2020 Departamento de Engenharia Mecânica
DZMITRY ZHYLTSOU
Dinâmica Estrutural: Implementação Numérica e Análise Computacional por Elementos Finitos
Universidade de Aveiro
2020 Departamento de Engenharia Mecânica
DZMITRY ZHYLTSOU
Dinâmica Estrutural: Implementação Numérica e Análise Computacional por Elementos Finitos
Dissertação apresentada à Universidade de Aveiro para cumprimento dos requisitos necessários à obtenção do grau de Mestre em Engenharia Mecânica, realizada sob a orientação científica do Doutor Joaquim Alexandre Mendes de Pinho da Cruz, Professor Auxiliar do Departamento de Engenharia Mecânica da Universidade de Aveiro
o júri
Presidente Professor Doutor Ricardo José Alves de Sousa Professor Auxiliar com Agregação do Departamento de Engenharia Mecânica da Universidade de Aveiro
Vogais Doutor Carlos André Soares Couto Investigador Auxiliar do Departamento de Engenharia Civil da Universidade de Aveiro (arguente)
Professor Doutor Joaquim Alexandre Mendes de Pinho da Cruz Professor Auxiliar do Departamento de Engenharia Mecânica da Universidade de Aveiro (orientador)
agradecimentos
Ao Professor Joaquim Alexandre Mendes de Pinho da Cruz, pelo seu apoio e orientação prestada no desenvolvimento deste trabalho.
palavras-chave
Vibrações, estruturas, vigas, barras, método de elementos finitos, resposta transitória, análise modal, MATLAB, GiD, SAP2000
resumo
Este trabalho visa aplicar o método de elementos finitos para a obtenção da resposta transitória ou permanente da vibração dos sistemas constituídos por vigas ou barras, com implementação numérica em MATLAB e com o apoio do GiD para a definição do problema e a posterior visualização de resultados. O modelo matemático da dinâmica das estruturas é obtido diretamente a partir da Segunda Lei de Newton, sendo assim um sistema de equações diferenciais que neste trabalho é desacoplado, para permitir a sobreposição modal nos sistemas lineares. Para se obter esse sistema de equações, com a aplicação do Princípio dos Trabalhos Virtuais, são instituídas as propriedades dinâmicas que definem um elemento finito: a rigidez, a inércia e o amortecimento, que depois são sobrepostas, de tal modo que satisfaçam o equilíbrio dinâmico e a compatibilidade de deslocamentos em qualquer nó que ligue esses elementos entre si. Embora sejam elementos unidimensionais, além de permitirem uma compreensão inicial mais fácil da dinâmica das estruturas, a sua utilização também pode ser vantajosa para diminuir o esforço computacional nos sistemas complexos, já com bastantes graus de liberdade, desde que seja empregado em domínios em que a formulação desses elementos seja representativa da mecânica do problema.
keywords
Vibrations, structures, beams, bars, finite element method, transient response, modal analysis, MATLAB, GiD, SAP2000
abstract
The present work aims to apply the finite element method to obtain a transient or steady-state vibration response of structures made up of beams or bars, with numerical implementation in MATLAB and with support of GiD in definition of the problem and the following display of results. The mathematical model of dynamic of structures is directly obtained from Newton’s Second Law of Motion, being a system of differential equations that in this work is decoupled to enable the application of modal superposition in linear systems. In order to obtain that system of equations, using the Principle of Virtual Work, the dynamic properties that define a finite element are established: stiffness, inertia and damping, which are subsequently superposed in a way that satisfy the dynamic equilibrium and compatibility of displacements at each nodal point. Despite being one-dimensional elements, they allow an easier initial understanding of structural dynamics phenomena and, beyond that, their use can also be advantageous to reduce computational cost in complex systems, already with a great number of degrees of freedom, as long as used in the domains where the formulation of those elements takes into account the mechanics of the problem.
ÍNDICE
1 INTRODUÇÃO ........................................................................................................... 1
2 VIBRAÇÃO DE SISTEMAS COM 1 GRAU DE LIBERDADE .............................. 2
2.1 Vibração livre não-amortecida .................................................................................. 3
2.2 Vibração livre amortecida ......................................................................................... 4
2.3 Vibração forçada amortecida .................................................................................... 6
2.3.1 Carregamento sinusoidal ........................................................................................ 6
2.3.2 Carregamento periódico ......................................................................................... 8
3 VIBRAÇÃO DE SISTEMAS COM 2 OU MAIS GRAUS DE LIBERDADE .......... 9
3.1 Edifício de dois andares ............................................................................................ 9
3.1.1 Ortogonalidade dos modos de vibração ............................................................... 11
3.2 Estruturas com amortecimento ............................................................................... 12
3.2.1 Transformação de coordenadas............................................................................ 13
3.2.2 Resposta do sistema através da sobreposição modal ........................................... 14
3.3 Condensação estática .............................................................................................. 15
4 ELEMENTOS DE VIGA E BARRA EM DINÂMICA ESTRUTURAL ................. 17
4.1 Barra sob cargas axiais ........................................................................................... 18
4.2 Viga sob cargas transversais e momentos fletores .................................................. 19
4.3 Viga sob momentos torsores ................................................................................... 23
4.4 Coeficientes de amortecimento ............................................................................... 23
4.5 Treliças 3D sob forças concentradas nos nós ......................................................... 24
4.5.1 Transformação de coordenadas............................................................................ 25
4.6. Estruturas 3D constituídas por elementos do tipo viga-barra ................................ 26
5 MODELAÇÃO DO PROGRAMA PRINCIPAL EM MATLAB ............................. 29
6 PROGRAMA APLICADO À RESOLUÇÃO DE PROBLEMAS ........................... 42
6.1 Cálculo das frequências/períodos naturais de sistemas .......................................... 42
6.1.1 Estrutura 2D (bidimensional) articulada .............................................................. 42
6.1.2 Estrutura 1D (unidimensional) reticulada ............................................................ 43
6.1.3 Estrutura 3D reticulada ........................................................................................ 45
6.2 Análise transitória de sistemas sujeitos a cargas dinâmicas ................................... 47
6.2.1 Estruturas articuladas ........................................................................................... 47
6.2.1.1 Estrutura 2D sem amortecimento ..................................................................... 47
6.2.1.2 Estrutura 2D com amortecimento ..................................................................... 52
6.2.1.3 Estrutura 3D sem amortecimento ..................................................................... 55
6.2.2 Estruturas reticuladas ........................................................................................... 58
6.2.2.1 Estrutura 2D sem amortecimento ..................................................................... 58
6.2.2.2 Estrutura 2D com amortecimento ..................................................................... 62
6.2.2.3 Estruturas 3D sem amortecimento .................................................................... 69
6.2.2.3.1 Exemplo 1 ...................................................................................................... 69
6.2.2.3.2 Exemplo 2 ...................................................................................................... 73
7 CONCLUSÕES ......................................................................................................... 80
APÊNDICE A – Fluxograma do programa principal ................................................... 81
APÊNDICE B – Algoritmos ......................................................................................... 84
REFERÊNCIAS .......................................................................................................... 102
LISTA DE FIGURAS
Figura 1: Sistema de massa m, mola de rigidez k, amortecedor de coeficiente de
amortecimento viscoso c e carga dinâmica F(t), para um deslocamento u(t) ....................... 2 Figura 2: (a) pórtico sob a ação de uma carga lateral concentrada, (b) reservatório de água
sob ação do vento e (c) viga encastrada-livre sob a ação do desbalanceamento de motor na
extremidade ........................................................................................................................... 2 Figura 3: Diagrama de corpo livre do modelo analítico ........................................................ 3 Figura 4: Soma e subtração das constantes de integração no Plano Complexo .................... 6 Figura 5: Sistema modelado com 2 GL (adaptado (2)) ......................................................... 9 Figura 6: Elemento do tipo viga-barra ................................................................................. 17
Figura 7: Treliça 2D no SAP2000 (adaptado (6)) ............................................................... 42
Figura 8: Viga encastrada-livre no SAP2000 (adaptado (8)) .............................................. 43
Figura 9: Estrutura 3D reticulada no SAP2000 (adaptado (8)) ........................................... 46 Figura 10: Estrutura 2D articulada no SAP2000 (adaptado (9)) ......................................... 48 Figura 11: Deslocamento do nó 4 (1 no GiD) em X, SAP2000 vs. MATLAB ................... 49 Figura 12: Deslocamento do nó 4 (1 no GiD) em Z, SAP2000 vs. MATLAB ................... 49
Figura 13: Deslocamento do nó 2 em X, SAP2000 vs. MATLAB ..................................... 49 Figura 14: Deslocamento do nó 2 em Z, SAP2000 vs. MATLAB ...................................... 50
Figura 15: Tensões axiais nas barras nos instantes 0,009 s e 0,282 s, SAP2000 ................ 50 Figura 16: Tensões axiais absolutas nas barras nos instantes 0,009 s e 0,282 s, MATLAB
............................................................................................................................................ .51
Figura 17: Deformada da estrutura no instante 0,282 s, SAP2000 vs. GiD ........................ 51 Figura 18: Estrutura 2D articulada no SAP2000 (adaptado (2)) ......................................... 52
Figura 19: Deslocamento do nó 1 (3 no GiD) em X, SAP2000 vs. MATLAB ................... 52
Figura 20: Deslocamento do nó 1 (3 no GiD) em Z, SAP2000 vs. MATLAB ................... 53
Figura 21: Deslocamento do nó 2 (4 no GiD) em Z, SAP2000 vs. MATLAB ................... 53 Figura 22: Tensões axiais nas barras no instante 1 s, SAP2000 .......................................... 54
Figura 23: Tensões axiais absolutas nas barras no instante 1 s, MATLAB ........................ 54 Figura 24: Deformada da estrutura no instante 1 s, SAP2000 vs. GiD ............................... 54 Figura 25: Estrutura 3D articulada no SAP2000 (adaptado (2)) ......................................... 55
Figura 26: Deslocamento do nó 3 (4 no GiD) em X, SAP2000 vs. MATLAB ................... 56 Figura 27: Deslocamento do nó 3 (4 no GiD) em Y, SAP2000 vs. MATLAB ................... 56 Figura 28: Deslocamento do nó 3 (4 no GiD) em Z, SAP2000 vs. MATLAB ................... 56
Figura 29: Tensões axiais nas barras no instante 0,1675 s, SAP2000 ................................. 57 Figura 30: Tensões axiais absolutas nas barras no instante 0,1675 s, MATLAB ............... 57
Figura 31: Deformada da estrutura no instante 0,1675 s, SAP2000 vs. GiD ...................... 58 Figura 32: Estrutura 2D reticulada no SAP2000 (adaptado (9)) ......................................... 58
Figura 33: Deslocamento do nó 2 (3 no GiD) em Z, SAP2000 vs. MATLAB ................... 59 Figura 34: Rotação do nó 2 (3 no GiD) em Y, SAP2000 vs. MATLAB ............................ 59 Figura 35: Deslocamento do nó 3 (5 no GiD) em Z, SAP2000 vs. MATLAB ................... 60 Figura 36: Rotação do nó 3 (5 no GiD) em Y, SAP2000 vs. MATLAB ............................ 60 Figura 37: Tensões axiais máximas absolutas no instante 1,916 s, MATLAB ................... 61
Figura 38: Tensões axiais nos pontos críticos, no instante 1,916 s, SAP2000 .................... 61 Figura 39: Deformada da estrutura no instante 1,916 s, GiD .............................................. 62 Figura 40: Deformada da estrutura no instante 1,916 s, SAP2000...................................... 62 Figura 41: Estrutura 2D reticulada no SAP2000 (adaptado (10)) ....................................... 63
Figura 42: Deslocamento do nó 3 (2 no GiD) em X, Análise 1, SAP2000 vs. MATLAB . 65 Figura 43: Deslocamento do nó 3 (2 no GiD) em X, Análise 2, SAP2000 vs. MATLAB . 66 Figura 44: Deslocamento do nó 3 (2 no GiD) em X, Análise 3, SAP2000 vs. MATLAB . 66 Figura 45: Tensões axiais máximas absolutas no instante 0,348 s, Análise 1, MATLAB .. 66 Figura 46: Tensões axiais nos pontos críticos no instante 0,348 s, Análise 1, SAP2000 .... 67
Figura 47: Deformada da estrutura no instante 0,348 s, Análise 1, SAP2000 vs. GiD ....... 67 Figura 48: Tensões axiais máximas absolutas no instante 5,993 s, Análise 2, MATLAB .. 67 Figura 49: Tensões axiais nos pontos críticos no instante 5,993 s, Análise 2, SAP2000 .... 68 Figura 50: Deformada da estrutura no instante 5,993 s, Análise 2, SAP2000 vs. GiD ....... 68 Figura 51: Tensões axiais máximas absolutas no instante 0,164 s, Análise 3, MATLAB .. 68
Figura 52: Tensões axiais nos pontos críticos no instante 0,164 s, Análise 3, SAP2000 .... 69 Figura 53: Deformada da estrutura no instante 0,164 s, Análise 3, SAP2000 vs. GiD ....... 69 Figura 54: Estrutura 3D reticulada (exemplo 1) no SAP2000 (adaptado (2)) ..................... 70
Figura 55: Deslocamento do nó 1 (2 no GiD) em X, SAP2000 vs. MATLAB ................... 71 Figura 56: Deslocamento do nó 1 (2 no GiD) em Z, SAP2000 vs. MATLAB ................... 71 Figura 57: Rotação do nó 1 (2 no GiD) em Y, SAP2000 vs. MATLAB ............................ 71 Figura 58: Tensões axiais máximas absolutas no instante 0,0444 s, MATLAB ................. 72
Figura 59: Tensões axiais nos pontos críticos no instante 0,0444 s, SAP2000 ................... 72 Figura 60: Deformada da estrutura no instante 0,0444 s, SAP2000 vs. GiD ...................... 72
Figura 61: Estrutura 3D reticulada (exemplo 2) no SAP2000 (adaptado (2)) ..................... 73 Figura 62: Deslocamento do nó 7 (9 no GiD) em Y, Análise 1, SAP2000 vs. MATLAB . 76
Figura 63: Deslocamento do nó 7 (24 no GiD) em Y, Análise 2, SAP2000 vs. MATLAB…
............................................................................................................................................. 76 Figura 64: Deslocamento do nó 7 (82 no GiD) em Y, Análise 3, SAP2000 vs. MATLAB
............................................................................................................................................ .77 Figura 65: Deslocamento do nó 5 (4 no GiD) em Y, Análise 1, SAP2000 vs. MATLAB . 77
Figura 66: Deslocamento do nó 5 (23 no GiD) em Y, Análise 2, SAP2000 vs. MATLAB
............................................................................................................................................ .77
Figura 67: Deslocamento do nó 5 (83 no GiD) em Y, Análise 3, SAP2000 vs. MATLAB
............................................................................................................................................ .78 Figura 68: Deformada da estrutura no instante 0,0538 s, Análise 3, SAP2000 vs. GiD ..... 78
Figura 69: Tensões axiais máximas absolutas no instante 0,0551 s, Análise 1, MATLAB
........................................................................................................................................... ..79 Figura 70: Diagrama da tensão axial máxima absoluta ao longo do elemento 7, no instante
0,0551 s, Análise 1, GiD ..................................................................................................... 79 Figura 71: Diagramas das tensões axiais nos 4 pontos (vertíces da secção quadrangular) ao
longo do elemento 7, no instante 0,0551s, Análise 1, SAP2000 ......................................... 79
LISTA DE TABELAS
Tabela 1: Frequências naturais do sistema da Figura 7 e 𝑑. 𝑝. (Referência (6) vs. Programa)
............................................................................................................................................. 43
Tabela 2: Períodos naturais do sistema da Figura 8 e 𝑑. 𝑝. (Referência (8) vs. Programa)..44 Tabela 3: Coordenadas e massas concentradas dos nós da estrutura da Figura 9 ............... 46 Tabela 4: Coordenadas e massas concentradas dos nós da estrutura da Figura 9 (cont.) .... 46
Tabela 5: Frequências naturais do sistema da Figura 9 e 𝑑. 𝑝. (SAP2000 vs. Programa) ... 47
Tabela 6: Frequências naturais do sistema da Figura 10 e 𝑑. 𝑝. (SAP2000 vs. Programa e
Referência (9) vs. Programa) ............................................................................................... 48 Tabela 7: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 10 e 𝑑. 𝑝. (SAP2000 vs. Programa) ......................................... 50
Tabela 8: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 18 e 𝑑. 𝑝. (SAP2000 vs. Programa) ......................................... 53 Tabela 9: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 25 e 𝑑. 𝑝. (SAP2000 vs. Programa) ......................................... 56
Tabela 10: Magnitude da força em vários instantes (histórico de carregamento) ............... 59 Tabela 11: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 32 e 𝑑. 𝑝. (SAP2000 vs. Programa) ......................................... 59
Tabela 12: Magnitude da força periódica não-sinusoidal em vários instantes nas 3 análises
diferentes ............................................................................................................................. 63 Tabela 13: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
1, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa) ........................ 64 Tabela 14: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
2, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa) ........................ 64
Tabela 15: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
3, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa) ........................ 65 Tabela 16: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 54 e 𝑑. 𝑝. (SAP2000 vs. Programa) ......................................... 70
Tabela 17: Magnitude da força triangular em vários instantes ............................................ 73 Tabela 18: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
1, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa) ............. 74
Tabela 19: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
2, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa) ............. 74 Tabela 20: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
3, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa) ............. 75
LISTA DE SÍMBOLOS
𝑚 massa equivalente do modelo
𝑐 coeficiente de amortecimento equivalente do modelo
𝑘 rigidez equivalente do modelo
𝜉 fator de amortecimento
𝑔 magnitude da força gravítica
N força normal
[Me] matriz de massa elementar
[Ce] matriz de amortecimento elementar
[Ke] matriz de rigidez elementar
[M] matriz de massa global
[C] matriz de amortecimento global
[K] matriz de rigidez global
{a𝑖} modo de vibração (vetor próprio) i
{φ𝑖} modo de vibração (vetor próprio) i normalizado
[A] matriz de modos de vibração do sistema
[Φ] matriz de modos de vibração do sistema normalizados
𝜆𝑖 valor próprio i
[Λ] matriz diagonal com valores próprios do sistema
𝜔n frequência natural não-amortecida
𝜔d frequência natural amortecida
𝜔𝑖 frequência natural do modo i
[Ω] matriz diagonal com frequências naturais do sistema
ω frequência de excitação
Χ fator de amplificação dinâmica
𝑦𝑖 fator de participação do modo i
𝜉𝑖 fator de amortecimento do modo i
𝑁𝑖 função de forma no grau de liberdade i
𝑃𝑖 esforço no grau de liberdade i do elemento no referencial local
𝛿𝑖 deslocamento no grau de liberdade i do elemento no referencial local
��𝑖 esforço no grau de liberdade i do elemento no referencial global
𝛿�� deslocamento no grau de liberdade i do elemento no referencial global
E módulo de Young
A área de secção
L comprimento longitudinal
ν coeficiente de Poisson
G módulo de corte
𝐼y segundo momento de área em relação ao eixo Y local
𝐼z segundo momento de área em relação ao eixo Z local
𝐼0 segundo momento polar de área
J constante de torção
�� densidade linear
LISTA DE ABREVIATURAS
MEF método de elementos finitos
GL grau de liberdade
MV modo de vibração
1D unidimensional
2D bidimensional
3D tridimensional
1
1 INTRODUÇÃO
Em muitas situações do dia a dia, as estruturas estão sujeitas a carregamentos que variam
com o tempo ou que são aplicados abruptamente, embora nestes últimos possam ser
constantes, pelo que surgem forças dissipativas e de inércia, induzindo as estruturas a
vibrar. Estes tipos de carregamento não podem ser desconsiderados na resposta transitória
do sistema, tal como na resposta permanente, se a frequência de excitação for superior a
cerca de um terço da frequência natural mais baixa da estrutura. Assim, a análise estática
das estruturas é válida para componentes de carregamento constantes ou com baixas
frequências na resposta permanente (1), sendo possível sobrepô-la nos sistemas lineares
com a resposta transitória associada aos efeitos dinâmicos, se estes existirem.
As estruturas em vibração apresentam tensões variáveis no tempo em virtude de
movimentos oscilatórios entre as partículas materiais que as constituem. Por conseguinte, o
projetista ou o engenheiro de estruturas deverá ter em conta, além dos critérios de falha
como o de von Mises, a resistência à fadiga dos elementos envolvidos e, assim,
dimensionar o sistema de modo condizente com as solicitações dinâmicas e estáticas
características ao fim a que se destina.
Para ser possível analisar o comportamento dinâmico das estruturas através do MEF
(método de elementos finitos), foram assumidas algumas idealizações, neste trabalho, para
as vigas e barras. Sendo elementos contínuos tridimensionais, estes podem ser descritos
apenas pelo seu eixo longitudinal e aproximados a sistemas discretos, contendo dois nós
por elemento, cada um com um ou vários GL (graus de liberdade) associados. O material
que os constitui foi considerado como homogéneo e isotrópico, respondendo linearmente a
cargas no regime elástico. Admitindo que estes elementos sofrem pequenas deformações,
aplicou-se a hipótese de sobreposição modal para se encontrar a resposta dinâmica das
estruturas. Geralmente este método numérico só se aplica em problemas lineares, em
estruturas com baixo amortecimento e é mais eficiente nos problemas com frequências de
excitação relativamente baixas. A integração direta é outra alternativa a este método e que
será imprescindível quando o primeiro não se adequa ou seja computacionalmente
impraticável.
No texto que se segue a esta introdução, em primeiro lugar, serão apresentados alguns
conceitos-base para os sistemas modelados com um único GL, mas que serão fundamentais
para a compreensão e a implementação do método numérico para o maior número de GL.
Posteriormente, será analisado um modelo simplificado de edifícios, com o número de GL
igual ao número de pisos, e, assim, apresentar-se-ão algumas propriedades dos sistemas
lineares que permitam aplicar a sobreposição modal. Prosseguir-se-á com a formulação de
matrizes elementares de rigidez, de amortecimento e de massa, dos elementos finitos de
viga e barra, com as devidas transformações de coordenadas, necessárias para se formar as
matrizes globais.
Enfim, o principal objetivo deste trabalho será o de obter a resposta dinâmica de algumas
estruturas, através do MEF com a sobreposição modal, comparando-a com a solução obtida
no programa SAP2000, e, assim, testar e validar o cálculo numérico do código
desenvolvido em MATLAB. A aplicação deste código com o pré-processador e pós-
-processador GiD facilitará a definição dos problemas e a posterior visualização dos
resultados.
2
2 VIBRAÇÃO DE SISTEMAS COM 1 GRAU DE LIBERDADE
Existem vários sistemas contínuos, e, portanto, com um número infinito de GL, cujo
comportamento possa ser definido apenas por um único GL (Figura 1). É o caso dos
modelos analíticos simplificados da Figura 2.
Figura 1: Sistema de massa m, mola de rigidez k, amortecedor de coeficiente de
amortecimento viscoso c e carga dinâmica F(t), para um deslocamento u(t)
Figura 2: (a) pórtico sob a ação de uma carga lateral concentrada, (b) reservatório de água
sob ação do vento e (c) viga encastrada-livre sob a ação do desbalanceamento de motor na
extremidade
Aplicando a 2.ª Lei de Newton ou o Princípio de D’Alembert ao modelo analítico, com
recurso ao diagrama de corpo livre (Figura 3), chega-se à equação diferencial completa:
𝑚�� + 𝑐�� + 𝑘𝑢 = 𝐹(𝑡), (2.1)
em que a única incógnita independente é o deslocamento em relação à posição de
equilíbrio 𝑢, multiplicada por termos constantes, pois é assumido que a força elástica
3
(−𝑘𝑢) resulta da ação de uma mola linear, diretamente proporcional ao seu alongamento.
De forma semelhante, a força dissipativa (−𝑐��) é diretamente proporcional à sua
velocidade de alongamento, visto que é considerado um amortecimento viscoso. Este tipo
de amortecimento não descreve o mecanismo real de dissipação ou de transformação de
energia mecânica das estruturas noutras formas de energia, tais como o calor, que ocorre
devido a fenómenos como o movimento nas juntas e a histerese no material. No entanto, é
em muitos casos uma aproximação aceitável para facilitar a abordagem matemática do
modelo.
Figura 3: Diagrama de corpo livre do modelo analítico
2.1 Vibração livre não-amortecida
Quando 𝐹(𝑡) e o coeficiente 𝑐 são iguais a zero, e o centro de massa do modelo é afastado
da sua posição de equilíbrio, o sistema vibra livremente, entre essa posição inicial extrema
e outro extremo do lado oposto ao inicial, se partir com uma velocidade inicial nula, ou,
entre os extremos mais afastados que no primeiro caso, se a velocidade inicial não for nula.
Nesta situação, a Equação (2.1) fica como:
𝑚�� + 𝑘𝑢 = 0. (2.2)
Duas soluções particulares independentes desta equação podem ser, por exemplo:
𝑢1 = 𝐴cos(𝜔n𝑡) e 𝑢2 = 𝐵sin(𝜔n𝑡).
Como a Equação (2.2) é diferencial homogénea linear de segunda ordem, a sua solução
geral homogénea é a sobreposição das duas soluções particulares:
𝑢GH(𝑡) = 𝑢1 + 𝑢2= 𝐴cos(𝜔n𝑡) + 𝐵sin(𝜔n𝑡).
As constantes de integração A e B dependem das condições iniciais do sistema, e como se
verá, de forma semelhante, na resposta da vibração amortecida. Substituindo 𝑢 na Equação
(2.2) pela solução particular 𝑢1, obtém-se:
𝐴cos(𝜔n𝑡)(−𝑚𝜔n
2 + 𝑘) = 0. (2.3)
𝑚��
𝑐��
4
Como a função 𝐴cos(𝜔n𝑡) não é sempre nula, então a equação
(−𝑚𝜔n
2 + 𝑘) = 0 (2.4)
tem de ser válida em qualquer instante. Logo, 𝜔n = √𝑘
𝑚 .
2.2 Vibração livre amortecida
Na prática, os sistemas reais têm sempre um amortecimento associado. O sistema vibra
com a diminuição da sua amplitude de oscilação até parar ou apenas volta até a sua posição
de equilíbrio após atingir a posição extrema. No primeiro caso o sistema está
subamortecido e noutro ou está sobreamortecido ou no limite criticamente amortecido.
Nesta situação, a Equação (2.1) fica como:
𝑚�� + 𝑐�� + 𝑘𝑢 = 0. (2.5)
Uma solução que satisfaz a Equação (2.5) poderá ser:
𝑢 = 𝐶𝑒𝑠𝑡 . (2.6)
Com a Equação (2.6) substitui-se na Equação (2.5) a incógnita 𝑢 e fica:
𝐶𝑒𝑠𝑡(𝑚𝑠2 + 𝑐𝑠 + 𝑘) = 0. (2.7)
Para a Equação (2.7) ser verdadeira, a equação caraterística tem de ser válida:
(𝑚𝑠2 + 𝑐𝑠 + 𝑘) = 0. (2.8)
É possível encontrar a solução da Equação (2.8), através da fórmula resolvente:
𝑠1,2 =−𝑐
2𝑚±√(
𝑐
2𝑚)2
−𝑘
𝑚. (2.9)
A solução geral homogénea da Equação (2.5) é:
𝑢GH(𝑡) = 𝐶1𝑒
𝑠1𝑡 + 𝐶2𝑒𝑠2𝑡, (2.10)
em que 𝐶1 e 𝐶2 são as constantes de integração, determinadas a partir das condições
inicias. Neste caso, as soluções, 𝑠1 e 𝑠2, são distintas. No caso do sistema criticamente
amortecido, as soluções coincidem:
5
𝑠1,2 =−𝑐
2𝑚= 𝑠0, (2.11)
pois
(𝑐
2𝑚)2
=𝑘
𝑚. (2.12)
Logo, é preciso encontrar outra solução geral homogénea para este caso, como:
𝑢GH(𝑡) = 𝐶1𝑒
𝑠0𝑡 + 𝐶2𝑡𝑒𝑠0𝑡. (2.13)
A partir da Equação (2.12) é obtido o coeficiente de amortecimento crítico:
𝑐cr = 2√𝑚𝑘 = 2𝑚𝜔n. (2.14)
Assim, o fator de amortecimento é definido através de:
𝜉 =
𝑐
𝑐cr. (2.15)
A maioria das estruturas tem 0,02 < 𝜉 < 0,1, pelo que será usada a Equação (2.10).
Assim, o termo dentro do radical da Equação (2.9) será negativo:
𝜔n
2(𝜉2 − 1) < 0
e então
𝑠1,2 = −𝜉𝜔n ± 𝑗𝜔n√1 − 𝜉2. (2.16)
Sendo
𝜔d = 𝜔n√1 − 𝜉2, (2.17)
substituindo 𝑠1 e 𝑠2 da Equação (2.9) na Equação (2.10), com a fórmula de Euler obtém-se:
𝑢GH(𝑡) = 𝑒
−𝜉𝜔n𝑡((𝐶1 + 𝐶2) cos(𝜔d𝑡) + 𝑗(𝐶1 − 𝐶2) sin(𝜔d𝑡)). (2.18)
Logo, 𝐶1 e 𝐶2 têm de ser números complexos conjugados, visto que 𝑢GH é uma grandeza
física real. Assim, a Equação (2.18) pode ser reescrita como:
𝑢GH(𝑡) = 𝑒
−𝜉𝜔n𝑡(𝐴 cos(𝜔d𝑡) + 𝐵 sin(𝜔d𝑡)). (2.19)
6
Através da Figura 4 sabe-se que:
𝐶1 + 𝐶2 = 𝐶cos(𝛼) e 𝑗(𝐶1 − 𝐶2) = 𝐶sin(𝛼).
Figura 4: Soma e subtração das constantes de integração no Plano Complexo
Assim, a Equação (2.18) pode ser reescrita de forma ainda mais compacta:
𝑢GH(𝑡) = 𝑒
−𝜉𝜔n𝑡𝐶cos(𝜔d𝑡 − 𝛼), (2.20)
em que as constantes 𝐶 e 𝛼 são obtidas a partir das condições iniciais, i.e., conhecendo a
posição inicial 𝑢0 e a velocidade inicial 𝑣0 do centro de massa do modelo:
𝛼 = arctan (𝑣0 + 𝑢0𝜉𝜔n𝜔d𝑢0
) 𝐶 = √𝑢02 +
(𝑣0 + 𝑢0𝜉𝜔n)2
𝜔d2.
2.3 Vibração forçada amortecida
A solução geral da resposta do modelo é a sobreposição linear da resposta transitória
𝑢GH(𝑡) e da resposta permanente up(𝑡):
𝑢(𝑡) = 𝑢GH(𝑡) + 𝑢p(𝑡). (2.21)
A solução da Equação (2.1) é a própria solução particular 𝑢p(𝑡), que coincide com a
solução geral 𝑢(𝑡) no regime permanente.
2.3.1 Carregamento sinusoidal
Se houver um carregamento periódico sinusoidal 𝐹(𝑡), sendo 𝐹0 a sua amplitude e 𝑤 a sua
frequência angular, a Equação (2.1) fica como:
𝑚�� + 𝑐�� + 𝑘𝑢 = 𝐹0 sin(𝑤𝑡). (2.22)
𝑗(𝐶1 − 𝐶2)
𝐶1 + 𝐶2
𝐶
𝛼
7
A solução particular da Equação (2.22) deverá ter a mesma natureza que 𝐹(𝑡) e é, por isso,
também sinusoidal e com a mesma frequência de excitação, tendo a amplitude 𝑈, mas
devido ao amortecimento tem um atraso (−𝜑) de resposta relativamente a esta:
𝑢p(𝑡) = 𝑈sin(𝑤𝑡 − (−𝜑)). (2.23)
A Equação (2.23) pode ainda ser escrita de outra forma equivalente:
𝑢p(𝑡) = 𝑈1sin(𝑤𝑡) + 𝑈2cos(𝑤𝑡), (2.24)
em que 𝑈1 = 𝑈cos(𝜑) e 𝑈2 = 𝑈sin(𝜑). Desta forma, as derivadas de primeira e de segunda ordem da solução particular são:
��p(𝑡) = 𝑈1𝑤cos(𝑤𝑡) − 𝑈2𝑤sin(𝑤𝑡) (2.25)
e
��p(𝑡) = −𝑈1𝑤2sin(𝑤𝑡) − 𝑈2𝑤
2cos(𝑤𝑡). (2.26)
Substituindo na Equação (2.22) o termo 𝑢 e as suas derivadas com as de solução particular,
encontram-se as constantes 𝑈1 e 𝑈2:
𝑈1 =(𝑘 −𝑚𝑤2)𝐹0
(𝑘 − 𝑚𝑤2)2 + (𝑐.𝑤)2 𝑈2 =
−(𝑐.𝑤)𝐹0(𝑘 − 𝑚𝑤2)2 + (𝑐. 𝑤)2
.
Como 𝑈 = √𝑈12 + 𝑈2
2, então:
𝑈 =𝐹0
√(𝑘 −𝑚𝑤2)2 + (𝑐.𝑤)2 ou 𝑈 =
𝐹0𝑘𝛸,
em que 𝛸 é o fator de amplificação dinâmica e 𝑟 a razão de frequências:
𝛸 =1
√(1 − 𝑟2)2 + (2𝜉𝑟)2 𝑟 =
𝑤
𝜔n.
Como se pode verificar, o deslocamento máximo do sistema não é só a razão estática de 𝐹0
e 𝑘, mas sim esta é “corrigida” pelo fator 𝛸, que a relaciona com o efeito das forças de
inércia. É de notar que quando 𝑟 tende para um e 𝜉 para zero este fator tende, no limite,
para o infinito. As excessivas amplitudes de vibração podem constituir um problema, mas
poderão ser atenuadas elevando-se a frequência natural do sistema.
Finalmente, o termo (−𝜑) da solução particular é obtido de seguinte forma:
8
tan(𝜑) =𝑈2𝑈1⇔ (−𝜑) = −arctan (
−𝑐.𝑤
𝑘 −𝑚𝑤2) = arctan (
2𝜉𝑟
1 − 𝑟2) =𝜃.
2.3.2 Carregamento periódico
Se o carregamento atuante na estrutura não for sinusoidal, embora seja periódico com
frequência angular 𝑤, a sua análise pode ser aproximada pela sobreposição de respostas
individuais de carregamentos sinusoidais que compõem uma série de Fourier.
Sendo o período deste carregamento periódico 𝑇 =2𝜋
𝑤, então
𝐹(𝑡) =
𝑎02+∑(
∞
𝑛=1
𝑎𝑛 cos(𝑛𝑤𝑡) + 𝑏𝑛 sin(𝑛𝑤𝑡)) (2.27)
e
𝑢p(𝑡) =
𝑎02𝑘+1
𝑘∑{𝑎𝑛 cos(𝑛𝑤𝑡 − 𝜃) + 𝑏𝑛 sin(𝑛𝑤𝑡 − 𝜃)}
∞
𝑛=1
𝛸, (2.28)
em que
𝑎0 =2
𝑇∫ 𝐹(𝑡)𝑑𝑡𝑇
0
, 𝑎𝑛 =2
𝑇∫ 𝐹(𝑡) cos(𝑛𝑤𝑡) 𝑑𝑡𝑇
0
, 𝑏𝑛 =2
𝑇∫ 𝐹(𝑡) sin(𝑛𝑤𝑡) 𝑑𝑡𝑇
0
𝜃 = arctan (2𝜉𝑛𝑟
1 − (𝑛𝑟)2) e 𝛸 =
1
√(1 − (𝑛𝑟)2)2 + (2𝜉𝑛𝑟)2.
9
3 VIBRAÇÃO DE SISTEMAS COM 2 OU MAIS GRAUS DE LIBERDADE
Se a resposta do sistema não apresentar apenas uma forma de vibrar, não se pode reduzir o
movimento desse sistema a um modelo com apenas um GL. De facto, uma estrutura
constituída por um material contínuo terá um número infinito de MV (modos de vibração),
enquanto o seu modelo discreto terá tantos MV quantos GL com inércia associada.
Neste capítulo, será analisada a vibração lateral de um edifício através do seu modelo
simplificado: o conjunto de colunas de cada andar comporta-se como se fosse uma única
viga, a massa em cada andar está concentrada nos pisos, estes só podem ter um movimento
horizontal e no início é desprezado o amortecimento. Ao mesmo tempo, serão introduzidos
alguns conceitos fundamentais para a análise de estruturas.
3.1 Edifício de dois andares
Como exemplo de ilustração, apresenta-se na Figura 5 um edifício de dois andares sujeito a
duas cargas que variam com o tempo (a), que foi reduzido a um sistema de dois GL (b).
Figura 5: Sistema modelado com 2 GL (adaptado (2))
Através do diagrama de corpo livre da Figura 5 (c), são obtidas as equações de movimento:
𝑚1��1 + 𝑘1𝑢1 − 𝑘2(𝑢2 − 𝑢1) − 𝐹1(𝑡) = 0 (3.1)
e
𝐹2(𝑡)
𝐹1(𝑡) 𝑢1(𝑡)
𝑢2(𝑡)
𝐹1(𝑡)
𝐹2(𝑡) 𝐹2(𝑡)
𝐹1(𝑡)
𝑢2(𝑡)
𝑢1(𝑡)
𝑢1(𝑡)
𝑢2(𝑡)
𝑘2
𝑘1
𝑘1
𝑘2
𝑘1𝑢1
𝑘2(𝑢2 − 𝑢1)
𝑚1��1
𝑚2��2
𝑘2(𝑢1 − 𝑢2)
𝑚1
𝑚2
𝑚1
𝑚2
10
𝑚2��2 − 𝑘2(𝑢1 − 𝑢2) − 𝐹2(𝑡) = 0. (3.2)
As Equações (3.1) e (3.2) podem ser representadas na forma matricial como:
[M]{u} + [K]{u} = {F}, (3.3)
em que
{u} = {𝑢1𝑢2} , {u} = {
��1��2} , {F} = {
𝐹1(𝑡)𝐹2(𝑡)
}.
Se o vetor das cargas {F(t)} for igual a zero e se as condições iniciais não forem nulas,
tem-se uma vibração livre. Quando uma estrutura é afastada da sua posição de equilíbrio
com a configuração deformada de um dos seus modos naturais de vibração ({a𝑖}), ela vibra
abandonada a si mesma com a configuração daquele modo e com uma frequência
característica (𝑤𝑖) daquele MV (3). Assim, para cada MV i:
{u𝑖} = {a𝑖}cos(𝑤𝑖𝑡 − 𝛼𝑖), (3.4)
em que 𝛼𝑖 = 0, pois no instante inicial: {u𝑖} = {a𝑖} e {u𝑖} = {0 0}T.
Com a hipótese de sobreposição modal, a deformada da estrutura pode ser obtida pela
combinação linear dos modos de vibração, cada um dos quais afetado pelo fator de
participação na resposta dinâmica:
{u} = ∑ 𝑦𝑖{a𝑖}
𝑁GL=2
𝑖=1
. (3.5)
Substituindo {u} e {u} na Equação (3.3) por {u𝑖} da Equação (3.4) e sua segunda derivada
{u𝑖} = −𝑤𝑖2{a𝑖}cos(𝑤𝑖𝑡),
obtém-se
[[K] − 𝑤𝑖2[M]]{a𝑖}cos(𝑤𝑖𝑡) = 0. (3.6)
A solução trivial é {a𝑖}={0 0}T. De outra forma, para a Equação (3.6) ser válida:
[[K] − 𝑤𝑖2[M]]{a𝑖} = 0 (3.7)
ou
[[K] − 𝜆𝑖[M]]{a𝑖} = 0. (3.8)
A solução da Equação (3.7) requer que o determinante da matriz fator do {a𝑖} seja zero:
|[K] − 𝑤𝑖2[M]| = 0. (3.9)
11
No problema dado, da Equação (3.9) resulta um polinómio de 2.º grau, cuja solução são as
duas frequências naturais de vibração, 𝑤1 e 𝑤2. Com estas pode-se então encontrar os
respetivos MV, {a1} e {a2}, arbitrando para cada vetor {a𝑖} um único valor de
deslocamento para o sistema da Equação (3.7) ser possível e determinado. Deste modo,
pode-se concluir que o vetor {a𝑖} possui valores que são relativos um ao outro, mas que na
resposta dinâmica será multiplicado pelo fator 𝑦𝑖, como já foi visto na Equação (3.5). Este
método direto de resolução dos problemas de valores e vetores próprios, como o da
Equação (3.8), revela-se adequado para um sistema simples, mas para um sistema com
muitos GL torna-se impraticável e é comum usar métodos numéricos mais eficientes. É o
caso da função eig da ferramenta computacional do MATLAB. Importa referir também
que na presença do amortecimento, sendo este relativamente baixo (𝜉𝑖 < 0,1) na maior
parte das estruturas, através da Equação (2.17) pode-se concluir que 𝑤d,𝑖(frequência
natural amortecida do modo i) ≈ 𝑤𝑖, e, assim, de forma semelhante, {ad,𝑖}(modo de
vibração amortecido i) ≈ {a𝑖}.
3.1.1 Ortogonalidade dos modos de vibração
A Equação (3.7) pode ser rescrita de seguinte forma:
[K]{a𝑖} = 𝑤𝑖2[M]{a𝑖}. (3.10)
Esta equação dá uma perceção estática do problema dinâmico, considerando que {a𝑖} é a
deformada estática causada pelas forças do lado direito desta equação. Esta interpretação
possibilita o uso do teorema geral de estática em sistemas lineares.
A Equação na forma matricial (3.10) é equivalente a duas equações para o problema dado:
(𝑘1 + 𝑘2)𝑎1,𝑖 − 𝑘2𝑎2,𝑖 = 𝑤𝑖2𝑚1𝑎1,𝑖
−𝑘2(𝑎1,𝑖 − 𝑎2,𝑖) = 𝑤𝑖2𝑚2𝑎2,𝑖,
em que os índices do deslocamento modal 𝑎𝑛,𝑖 indicam o seu GL e o MV, respetivamente.
De acordo com o teorema de Betti, numa estrutura sob dois sistemas de forças com os
deslocamentos correspondentes, o trabalho realizado pelo primeiro sistema de forças nos
deslocamentos provocados pelo segundo sistema de forças é igual ao trabalho realizado
pelo segundo sistema de forças nos deslocamentos provocados pelo primeiro sistema de
forças. Neste caso, cada MV representa um sistema. Aplicando o teorema de Betti:
𝑤12𝑚1𝑎1,1𝑎1,2 + 𝑤1
2𝑚2𝑎2,1𝑎2,2 = 𝑤22𝑚1𝑎1,2𝑎1,1 + 𝑤2
2𝑚2𝑎2,2𝑎2,1,
que é igual a
(𝑤12 − 𝑤2
2)(𝑚1𝑎1,1𝑎1,2 +𝑚2𝑎2,1𝑎2,2) = 0.
Se 𝑤1 for diferente de 𝑤2, então
(𝑚1𝑎1,1𝑎1,2 +𝑚2𝑎2,1𝑎2,2) = 0,
12
que é a relação ortogonal entre os MV para o sistema de 2 GL. Assim, se [M] for diagonal:
∑𝑚𝑘𝑎𝑘,𝑖𝑎𝑘,𝑗
𝑁GL
𝑘=1
= 0, 𝑖 ≠ 𝑗.
No caso geral:
{{a𝑗}
T[M]{a𝑖} = 0
{a𝑗}T[K]{a𝑖} = 0
, 𝑖 ≠ 𝑗
e
{{a𝑗}
T[M]{a𝑖} = 𝑚𝑖
{a𝑗}T[K]{a𝑖} = 𝑘𝑖
, 𝑖 = 𝑗,
em que 𝑚𝑖 e 𝑘𝑖 são, respetivamente, a massa e a rigidez generalizada do modo 𝑖 de
vibração. Os MV podem ser normalizados pela massa e se [M] for diagonal, então:
𝜑𝑛,𝑖 =𝑎𝑛,𝑖
√∑ 𝑚𝑘𝑎𝑘,𝑖2𝑁GL𝑘=1
.
No caso geral:
𝜑𝑛,𝑖 =𝑎𝑛,𝑖
√{a𝑖}T[M]{a𝑖}.
De forma semelhante a {a𝑖}, estes MV normalizados ({φ𝑖}) têm as seguintes propriedades:
{{φ𝑗}
T[M]{φ𝑖} = 0
{φ𝑗}T[K]{φ𝑖} = 0
, 𝑖 ≠ 𝑗
e
{{φ𝑗}
T[M]{φi} = 1
{φ𝑗}T[K]{φ𝑖} = 𝑤𝑖
2, 𝑖 = 𝑗
ou {[Φ]T[M][Φ] = [I]
[Φ]T[K][Φ] = [Ω].
3.2 Estruturas com amortecimento
Se o amortecimento da estrutura não for desprezável, a equação do movimento na forma
matricial fica como:
[M]{u} + [C]{u} + [K]{u} = {F}. (3.11)
Para tornar a Equação (3.11) num sistema de equações independentes, é preciso que a
matriz de amortecimento também satisfaça a condição de ortogonalidade:
13
{{φ𝑗}
T[C]{φ𝑖} = 0, 𝑖 ≠ 𝑗
{φ𝑗}T[C]{φ𝑖} = 𝑐𝑖 = 2𝜉𝑖𝑤𝑖 , 𝑖 = 𝑗 ou [Φ]T[C][Φ] = [c],
em que 𝑐𝑖 é o coeficiente de amortecimento generalizado do modo 𝑖 de vibração e [c] é
uma matriz diagonal com esses coeficientes.
Uma forma simples de obter [C] que satisfaça tal condição é pela combinação linear das
matrizes de rigidez e de massa, usando constantes de amortecimento 𝛼 e β:
[C] = 𝛼[K] + 𝛽[M], (3.12)
sendo
{𝛼𝛽} = 2 [
1 𝑤12
1 𝑤max2]
−1
{𝜉1𝜉max
}. (3.13)
O amortecimento 𝜉𝑖 pode ser determinado experimentalmente durante a vibração livre
amortecida para cada MV, provocando essa vibração com a frequência natural desse modo.
Neste tipo de vibração, o decaimento de oscilação entre os sucessivos picos é aproximado
a uma função exponencial (𝑒−𝜉𝑖ω𝑖𝑡), tal como se pode verificar através da Equação (2.20).
Assim, aplicando a operação logarítmica na relação entre os sucessivos picos de
deslocamento ou de aceleração obtidos experimentalmente, tem-se para cada MV:
𝛿𝑖 = ln (𝑒−𝜉𝑖ω𝑖𝑡1
𝑒−𝜉𝑖ω𝑖(𝑡1+𝑇d,𝑖)) = 𝜉𝑖ω𝑖𝑇d,𝑖 =
2𝜋𝜉𝑖
√1 − 𝜉𝑖2
,
que com 𝜉𝑖 < 0,1 pode ser aproximada por
𝛿𝑖 ≈ 2𝜋𝜉𝑖 , donde
𝜉𝑖 =𝛿𝑖2𝜋.
No entanto, o Amortecimento de Rayleigh, definido na Equação (3.12), só tem duas
constantes, pelo que só será preciso calcular o amortecimento de dois modos: sendo 𝜉1 e
𝜉max o amortecimento presente na estrutura durante a vibração com uma frequência natural
mais baixa (𝑤1) e com uma frequência natural superior a frequências de interesse para a
análise (𝑤max), respetivamente.
3.2.1 Transformação de coordenadas
A transformação de coordenadas já foi definida, para o caso particular de 2 GL, na
Equação (3.5), correspondendo no caso geral a:
{u} = ∑ 𝑦𝑖(𝑡){a𝑖}
𝑁modos
𝑖=1
= [A]{y(t)}, 𝑁modos ≤ 𝑁GL (3.14)
14
ou
{u} = ∑ 𝑞𝑖(𝑡){φ𝑖}
𝑁modos
𝑖=1
= [Φ]{q(t)}, 𝑁modos ≤ 𝑁GL. (3.15)
Usando a transformação de coordenadas da Equação (3.15) na Equação (3.11), tem-se:
[M][Φ]{q(t)} + [C][Φ]{q(t)} + [K][Φ]{q(t)} = {F}. (3.16)
Multiplicando à esquerda por [Φ]T, tem-se
[Φ]T [M][Φ]{q(t)} + [Φ]T [C][Φ]{q(t)} + [Φ]T[K][Φ]{q(t)} = [Φ]T{F} (3.17)
e da propriedade de ortogonalidade dos modos de vibração resulta
[I]{q(t)} + [c]{q(t)} + [Ω]{q(t)} = [Φ]T{F} (3.18)
ou
��𝑖(𝑡) + 2𝜉𝑖𝑤𝑖��𝑖(𝑡) + 𝑤𝑖
2𝑞𝑖(𝑡) = {φ𝑖}T{F} = 𝑃𝑖(𝑡), 𝑖 = 1, 2, . . . , 𝑁modos. (3.19)
A Equação (3.18) é um sistema de equações desacopladas de movimento em função das
coordenadas modais, tal como se pode constatar na Equação (3.19).
3.2.2 Resposta do sistema através da sobreposição modal
Para ser possível calcular a resposta global do sistema ({u}) pela Equação (3.15), é
necessário calcular o fator de participação de cada modo (𝑞𝑖(𝑡)), ou seja, de cada sistema
generalizado da Equação (3.19). Este fator é igual à própria solução da resposta geral de
uma vibração amortecida de um GL (ver a Equação (2.21)). Neste trabalho, foi admitido
que o carregamento atuante na estrutura é geral, mas que foi aproximado por um conjunto
de 𝑁s carregamentos sinusoidais da série de Fourier, com um período T igual ao intervalo
de tempo da análise. Deste modo, a solução particular foi baseada na Equação (2.28),
obtendo-se, assim, a solução geral para cada MV a partir da Equação (2.21):
𝑞𝑖(𝑡) = 𝑒−𝜉𝑖ω𝑖𝑡(𝐴𝑖 cos(ωd,𝑖𝑡) + 𝐵𝑖 sin(ωd,𝑖𝑡)) +𝑎0,𝑖2𝑘𝑖
+1
𝑘𝑖∑{𝑎𝑛,𝑖 cos(𝑛𝑤𝑡 − 𝜃𝑖) + 𝑏𝑛,𝑖 sin(𝑛𝑤𝑡 − 𝜃𝑖)}
𝑁s
𝑛=1
𝛸𝑛,𝑖, (3.20)
em que
𝜔d,𝑖 = 𝜔𝑖√1 − 𝜉𝑖2 , 𝑘𝑖 = 𝑤𝑖
2; , 𝑎0,𝑖 =2
𝑇∫ 𝑃𝑖(𝑡)𝑑𝑡𝑇
0
15
𝑎𝑛,𝑖 =2
𝑇∫ 𝑃𝑖(𝑡) cos(𝑛𝑤𝑡) 𝑑𝑡𝑇
0
, 𝑏𝑛,𝑖 =2
𝑇∫ 𝑃𝑖(𝑡) sin(𝑛𝑤𝑡) 𝑑𝑡𝑇
0
, 𝜃𝑖 = arctan (2𝜉𝑖𝑟𝑖
1 − (𝑛𝑟𝑖)2)
e onde as constantes 𝐴𝑖 e 𝐵𝑖 são obtidas a partir das condições iniciais:
𝐴𝑖 = {φ𝑖}T[M]{u0} −
𝑎0,𝑖2𝑘𝑖
−1
𝑘𝑖∑{𝑎𝑛,𝑖 cos(−𝜃𝑖) + 𝑏𝑛,𝑖 sin(−𝜃𝑖)}
𝑁s
𝑛=1
𝛸𝑛,𝑖
𝐵𝑖 = 𝜉𝑖𝜔𝑖𝐴𝑖 + {φ𝑖}T[M]{u0} −
𝑤
𝑘𝑖𝜔d,𝑖∑{𝑏𝑛,𝑖 cos(−𝜃𝑖) − 𝑎𝑛,𝑖 sin(−𝜃𝑖)}𝑛
𝑁s
𝑛=1
𝛸𝑛,𝑖.
O código usado para implementar esta solução de resposta no MATLAB é apresentado na
sub-rotina “fator_part_v_amort.m”. A solução geral durante a vibração livre amortecida
para cada modo pode-se obter de uma forma mais simples, visto que esta é igual à solução
geral homogénea. É sugerido usar um valor 𝑁s tão alto quanto maior forem:
• As frequências naturais do sistema;
• A frequência local da excitação;
• A divergência da excitação em relação às sinusoides;
• O intervalo de tempo da análise.
É também possível fazer várias análises consecutivas, se T for demasiado alto, diminuindo
desta forma o valor 𝑁s e, assim, o esforço computacional.
3.3 Condensação estática
Geralmente os modelos contêm um grande número de GL, contendo o mesmo número de
MV, mas não é praticável obter todos estes últimos para se alcançar uma razoável precisão
no cálculo da resposta dinâmica. Deste modo, ao se diminuir o número de MV, diminui-se
o cálculo numérico a efetuar, nomeadamente, a resolver-se o problema de valores e vetores
próprios, e a determinar-se os fatores de participação para cada modo.
Um dos métodos para atingir esse objetivo consiste em escolher os GL, que além de terem
a sua própria massa associada, irão também concentrar a massa pertencente aos GL livres
das forças de inércia e estes estarão restritos a movimentar-se apenas com uma relação
estática com os primeiros. Desta forma, após a assemblagem das matrizes de massa e de
rigidez globais, é preciso reorganizar os termos que as constituem, de forma a isolar na
matriz de massa uma submatriz com os termos não-nulos. Nesta condição, a Equação
(3.10) fica como:
[[K1,1] [K1,2]
[K2,1] [K2,2] ] {{a𝑖}1{a𝑖}2
} = 𝑤𝑖2 [[M1,1] [0]
[0] [0] ] {{a𝑖}1{a𝑖}2
}, (3.21)
𝑤 =2𝜋
𝑇 , 𝑟𝑖 =
𝑤
𝑤𝑖 e 𝛸𝑛,𝑖 =
1
√(1 − (𝑛𝑟𝑖)2)2 + (2𝜉𝑖𝑛𝑟𝑖)2
16
em que {a𝑖}1 e {a𝑖}2 são deslocamentos dos GL com massa e sem massa, respetivamente.
A Equação (3.21) é equivalente a
{[K1,1]{a𝑖}1 + [K1,2]{a𝑖}2 = 𝑤𝑖
2[M1,1]{a𝑖}1
[K2,1]{a𝑖}1 + [K2,2]{a𝑖}2 = 0⇔ {a𝑖}2 = −[K2,2]−1[K2,1]{a𝑖}1
. (3.22)
(3.23)
Substituindo os termos {a𝑖}2 da Equação (3.23) na Equação (3.22), obtém-se {a𝑖}1:
{a𝑖}1 = 𝑤𝑖2[Kaux]
−1[M1,1]{a𝑖}1,
em que
[Kaux] = [[K1,1] − [K1,2][K2,2]−1[K2,1]].
O problema de valores e vetores próprios ficará, desta maneira, reduzido a:
[[Kaux]−𝜆𝑖[M1,1]] {a𝑖}1 = 0.
17
4 ELEMENTOS DE VIGA E BARRA EM DINÂMICA ESTRUTURAL
Neste capítulo, são expostos os métodos para determinar as matrizes de rigidez e de massa
globais que descrevem o comportamento dinâmico de um sistema composto por elementos
unidimensionais de vigas e barras com dois nós. A dissipação de energia mecânica estará
relacionada com a matriz de amortecimento global definida no capítulo anterior, pois, seria
muito difícil de determinar as características de amortecimento de cada elemento
isoladamente.
Enquanto nas estruturas articuladas, como é o caso das treliças, existe apenas a propagação
de esforços normais, nas estruturas reticuladas existe também a transmissão dos momentos
fletores e torsores, e dos esforços transversos.
Na Figura 6, representa-se o caso geral de um elemento, no referencial local 𝑂𝑥1𝑦1𝑧1 e
coincidente com o referencial global 𝑂𝑋𝑌𝑍, com os esforços nodais internos e os
deslocamentos nos seus GL. Já o outro elemento, no referencial local 𝑂𝑥2𝑦2𝑧2, está com
uma rotação de θ no eixo 𝑍 em relação ao referencial global. Os esforços 𝑃𝑖 ou os deslocamentos 𝛿𝑖 de um elemento no seu referencial local estão organizados por
nó de tal forma a serem considerados primeiro as forças ou os deslocamentos lineares em
𝑥, 𝑦 e 𝑧, e depois os momentos ou as rotações em 𝑥, 𝑦 e 𝑧.
Figura 6: Elemento do tipo viga-barra
𝑃1, 𝛿1 𝑃4, 𝛿4 𝑃7, 𝛿7 𝑃10, 𝛿10
𝑃3, 𝛿3 𝑋 ≡ 𝑥1
𝑃6, 𝛿6
𝑃2, 𝛿2
𝑃5, 𝛿5
𝑃8, 𝛿8 𝑃9, 𝛿9
𝑃11, 𝛿11 𝑃12, 𝛿12
𝑍 ≡ 𝑧1
𝑌 ≡ 𝑦1
𝑧2
𝑥2
𝑦2
𝜃
18
4.1 Barra sob cargas axiais
Quando sobre a barra atuam as forças 𝑃1 e 𝑃7, e ocorrem os deslocamentos 𝛿1e 𝛿7, a
solução da equação diferencial que rege este problema para uma barra uniforme é:
𝑢 =
𝑃. 𝑥
𝐴. 𝐸+ 𝑐1 = 𝑐2𝑥 + 𝑐1, (4.1)
em que 𝑃 é o esforço axial na barra, 𝑢 é o deslocamento longitudinal em 𝑥, e 𝑐1 e 𝑐2 são as
constantes de integração. Assim, obtêm-se as seguintes funções de forma para cada GL:
{𝑢(0) = 𝛿1 = 1
𝑢(𝐿) = 𝛿7 = 0 ⇒ 𝑁1 = 𝑢 = 1 −
𝑥
𝐿 (4.2)
{𝑢(0) = 𝛿1 = 0
𝑢(𝐿) = 𝛿7 = 1 ⇒ 𝑁7 = 𝑢 =
𝑥
𝐿. (4.3)
Assumindo que a barra está em equilíbrio estático, com 𝛿1 = 1 e 𝛿7 = 0, havendo um
deslocamento virtual 𝛿7 no nó 2, o trabalho externo é igual a
𝑊E = 𝑃7 𝛿7 = (𝑘7,1𝛿1)𝛿7 = 𝑘7,1𝛿7
e o trabalho interno a
𝑊I = ∫ 𝑑𝑊I
𝐿
0
= ∫ 𝑃. 𝑑𝐿
0
(𝑁7𝛿7) = ∫ (𝐴. 𝐸. 𝑁1′). (𝑁7
′𝛿7𝑑𝑥)𝐿
0
= 𝛿7∫ 𝐴. 𝐸.𝑁1′𝑁7
′𝑑𝑥𝐿
0
.
De acordo com o Princípio dos Trabalhos Virtuais:
𝑊E = 𝑊I⇔𝑘7,1𝛿7 = 𝛿7∫ 𝐴. 𝐸.𝑁1′𝑁7
′𝑑𝑥𝐿
0
,
pelo que
𝑘7,1 = ∫ 𝐴. 𝐸.𝑁1′𝑁7
′𝑑𝑥𝐿
0
.
Da mesma forma, obtêm-se os coeficientes de rigidez para cada GL:
𝑘𝑖,𝑗 = ∫ 𝐴. 𝐸.𝑁𝑖′𝑁𝑗
′𝑑𝑥𝐿
0
, 𝑖, 𝑗 = {1,7}
e se a barra for uniforme
𝑘𝑖,𝑗 = 𝐴. 𝐸 ∫ 𝑁𝑖′𝑁𝑗
′𝑑𝑥𝐿
0
, 𝑖, 𝑗 = {1,7},
em que 𝑘𝑖,𝑗 representa o esforço que é exercido no GL 𝑖, para um deslocamento unitário do
GL 𝑗, enquanto nos restantes GL o deslocamento é nulo.
19
A matriz de rigidez para uma barra uniforme fica como:
[Ke] = [
𝑘1,1 𝑘1,7𝑘7,1 𝑘7,7
] =𝐴. 𝐸
𝐿|e
[1 −1−1 1
].
De forma semelhante, assumindo que a barra está em equilíbrio dinâmico, com ��1 = 1 e
��7 = 0, havendo um deslocamento virtual 𝛿7 no nó 2:
𝑊E = 𝑃7 𝛿7 = (𝑚7,1��1)𝛿7 = 𝑚7,1𝛿7.
Sendo 𝑓I a força de inércia por unidade de comprimento, sabe-se que
𝑑𝑊I = (𝑓I𝑑𝑥)(𝛿7𝑁7) = (��𝑁1��1𝑑𝑥)(𝛿7𝑁7) = 𝛿7��𝑁1𝑁7𝑑𝑥,
pelo que
𝑊I = ∫ 𝛿7��𝑁1𝑁7𝑑𝑥 = 𝛿7∫ ��𝑁1𝑁7𝑑𝑥
𝐿
0
.𝐿
0
Sendo
𝑊E = 𝑊I⇔𝑚7,1𝛿7 = 𝛿7∫ ��𝑁1𝑁7𝑑𝑥,𝐿
0
então
𝑚7,1 = ∫ ��𝑁1𝑁7𝑑𝑥.𝐿
0
Se a barra for uniforme, então
𝑚𝑖,𝑗 = ��∫ 𝑁𝑖𝑁𝑗𝑑𝑥, 𝑖, 𝑗 = {1,7} 𝐿
0
e
[Me] = [
𝑚1,1 𝑚1,7
𝑚7,1 𝑚7,7] =
(��. 𝐿)e
6[2 11 2
],
em que 𝑚𝑖,𝑗 representa o esforço que é exercido no GL 𝑖, tendo o GL 𝑗 uma aceleração
unitária e os restantes GL uma aceleração nula.
4.2 Viga sob cargas transversais e momentos fletores
Quando sobre a viga atuam os esforços 𝑃2, 𝑃6, 𝑃8 e 𝑃12, e ocorrem os deslocamentos 𝛿2,
𝛿6, 𝛿8 e 𝛿12, as seguintes equações diferenciais descrevem o problema:
𝐸. 𝐼z
𝑑2𝑢
𝑑𝑥2= 𝑀(𝑥) (4.4)
𝑑(𝑀(𝑥))
𝑑𝑥= 𝑉(𝑥) (4.5)
20
𝐸. 𝐼z
𝑑4𝑢
𝑑𝑥4= 𝑝(𝑥), (4.6)
em que 𝑢 é o deslocamento transversal no eixo 𝑦, 𝑀(𝑥) é o momento de flexão no eixo 𝑧,
𝑉(𝑥) é a força transversal de corte no eixo 𝑦 e 𝑝(𝑥) é a carga transversal distribuída no
eixo 𝑦. Se 𝑝(𝑥) = 0, então a Equação (4.6) fica como:
𝑑4𝑢
𝑑𝑥4= 0. (4.7)
A integração da Equação (4.7) resulta em:
𝑢 =
1
6𝑐1𝑥
3 +1
2𝑐2𝑥
2 + 𝑐3𝑥 + 𝑐4, (4.8)
em que 𝑐1, 𝑐2, 𝑐3 e 𝑐4 são as constantes de integração, que dependem das condições de
fronteira. Assim, as funções de forma obtêm-se como:
{
𝑑𝑢
𝑑𝑥|𝑥=0
𝑑𝑢
𝑑𝑥|𝑥=𝐿
𝑢(𝑥 = 0)
𝑢(𝑥 = 𝐿)}
= {
𝑎1𝑎2𝑎3𝑎4
}⇒ [
0 0 1 0𝐿2/2 𝐿 1 00
𝐿3/60
𝐿2/20𝐿
11
] {
𝑐1𝑐2𝑐3𝑐4
} = {
𝑎1𝑎2𝑎3𝑎4
} (4.9)
{
𝑐1𝑐2𝑐3𝑐4
} = [
0 0 1 0𝐿2/2 𝐿 1 00
𝐿3/60
𝐿2/20𝐿
11
]
−1
{
𝑎1𝑎2𝑎3𝑎4
}. (4.10)
Se 𝛿2 = 1, 𝛿6 = 0, 𝛿8 = 0 e 𝛿12 = 0, então 𝑢 = 𝑁2 e
{
𝑎1𝑎2𝑎3𝑎4
} = {
0010
}.
Através da Equação (4.10), encontram-se as constantes de integração e a função de forma
para o GL 2:
𝑁2 = 1 − 3(𝑥
𝐿)2
+ 2(𝑥
𝐿)3
.
As restantes funções de forma obtêm-se de forma semelhante:
𝑁6 = 𝑥 (1 −𝑥
𝐿)2
21
𝑁8 = 3(𝑥
𝐿)2
− 2(𝑥
𝐿)3
𝑁12 = (𝑥
𝐿)2
(𝑥
𝐿− 1).
Assumindo que a viga está em equilíbrio estático, com 𝛿6 = 1 e restantes GL
deslocamentos nulos, havendo um deslocamento virtual 𝛿2 no nó 1, o trabalho externo é
igual a:
𝑊E = 𝑃2 𝛿2 = (𝑘2,6𝛿6)𝛿2 = 𝑘2,6𝛿2
e o trabalho interno
𝑊I = ∫ 𝑀(𝑥)𝑑𝜃𝐿
0
.
Sabe-se que
𝑀(𝑥) = 𝐸. 𝐼z𝑑2(𝑁6)
𝑑𝑥2
e
𝑑𝜃
𝑑𝑥=𝑑2(𝑁2𝛿2)
𝑑𝑥2= 𝑁2
′′𝛿2,
pelo que
𝑊I = ∫ 𝐸. 𝐼z. 𝑁6′′𝑁2
′′𝛿2𝑑𝑥 = 𝛿2∫ 𝐸. 𝐼z. 𝑁6′′𝑁2
′′𝑑𝑥𝐿
0
𝐿
0
.
De acordo com o Princípio dos Trabalhos Virtuais:
𝑊E = 𝑊I⇔𝑘2,6𝛿2 = 𝛿2∫ 𝐸. 𝐼z. 𝑁6′′𝑁2
′′𝑑𝑥𝐿
0
,
pelo que
𝑘2,6 = ∫ 𝐸. 𝐼z. 𝑁6′′𝑁2
′′𝑑𝑥𝐿
0
.
Os coeficientes de rigidez para cada GL são:
𝑘𝑖,𝑗 = ∫ 𝐸. 𝐼z. 𝑁𝑖′′𝑁𝑗
′′𝑑𝑥𝐿
0
, 𝑖, 𝑗 = {2,6,8,12}.
Se a barra for uniforme, então
𝑘𝑖,𝑗 = 𝐸. 𝐼z∫ 𝑁𝑖′′𝑁𝑗
′′𝑑𝑥𝐿
0
, 𝑖, 𝑗 = {2,6,8,12}
e a matriz de rigidez fica como:
22
[Ke] =
[ 𝑘2,2 𝑘2,6 𝑘2,8 𝑘2,12
𝑘6,2 𝑘6,6 𝑘6,8 𝑘6,12𝑘8,2𝑘12,2
𝑘8,6𝑘12,6
𝑘8,8𝑘12,8
𝑘8,12𝑘12,12]
=𝐸. 𝐼z𝐿3|e
[
12 6𝐿 −12 6𝐿6𝐿 4𝐿2 −6𝐿 2𝐿2
−126𝐿
−6𝐿2𝐿2
12−6𝐿
−6𝐿4𝐿2
].
A matriz de massa consistente é determinada de forma semelhante. Assumindo que a viga
está em equilíbrio dinâmico, com ��6 = 1 e nos restantes GL acelerações nulas, havendo
um deslocamento virtual 𝛿2 no nó 1:
𝑊E = 𝑃2 𝛿2 = (𝑚2,6��6)𝛿2 = 𝑚2,6𝛿2
𝑑𝑊I = (𝑓I𝑑𝑥)(𝛿2𝑁2) = (��𝑁6��6𝑑𝑥)(𝛿2𝑁2) = 𝛿2��𝑁6𝑁2𝑑𝑥,
pelo que
𝑊I = ∫ 𝛿2��𝑁6𝑁2𝑑𝑥 = 𝛿2∫ ��𝑁6𝑁2𝑑𝑥.
𝐿
0
𝐿
0
Então
𝑊E = 𝑊I⇔𝑚2,6𝛿2 = 𝛿2∫ ��𝑁6𝑁2𝑑𝑥,𝐿
0
pelo que
𝑚2,6 = ∫ ��𝑁6𝑁2𝑑𝑥.𝐿
0
Se a barra for uniforme, então
𝑚𝑖,𝑗 = ��∫ 𝑁𝑖𝑁𝑗𝑑𝑥, 𝑖, 𝑗 = {2,6,8,12} 𝐿
0
e
[Me] = [
𝑚2,2 𝑚2,6𝑚2,8 𝑚2,12
𝑚6,2 𝑚6,6𝑚6,8 𝑚6,12
𝑚8,2
𝑚12,2
𝑚8,6
𝑚12,6
𝑚8,8
𝑚12,8
𝑚8,12
𝑚12,12
] =��. 𝐿
420|e
[
156 22𝐿 54 −13𝐿22𝐿 4𝐿2 13𝐿 −3𝐿2
54−13𝐿
13𝐿−3𝐿2
156−22𝐿
−22𝐿4𝐿2
].
Se a deformação da viga ocorrer no plano 𝑂𝑥𝑧, e não no plano 𝑂𝑥𝑦 como aconteceu até
agora, i.e., quando sobre a viga atuam os esforços 𝑃3, 𝑃5, 𝑃9 e 𝑃11, e ocorrem os
deslocamentos 𝛿3, 𝛿5, 𝛿9 e 𝛿11, a formulação dos coeficientes de rigidez e de massa é
semelhante, havendo a necessidade de efetuar a troca de sinais nestes, porque as funções de
forma do 3.º e 11.º GL mudam de sinal face às anteriores, ficando as matrizes definidas
como:
[Ke] =
[ 𝑘3,3 𝑘3,5 𝑘3,9 𝑘3,11
𝑘5,3 𝑘5,5 𝑘5,9 𝑘5,11𝑘9,3𝑘11,3
𝑘9,5𝑘11,5
𝑘9,9𝑘11,9
𝑘9,11𝑘11,11]
=𝐸. 𝐼y
𝐿3|e
[
12 −6𝐿 −12 −6𝐿−6𝐿 4𝐿2 6𝐿 2𝐿2
−12−6𝐿
6𝐿2𝐿2
126𝐿
6𝐿4𝐿2
]
e
23
[Me] = [
𝑚3,3 𝑚3,5𝑚3,9 𝑚3,11
𝑚5,3 𝑚5,5𝑚5,9 𝑚5,11
𝑚9,3
𝑚11,3
𝑚9,5
𝑚11,5
𝑚9,9
𝑚11,9
𝑚9,11
𝑚11,11
] =��. 𝐿
420|e
[
156 −22𝐿 54 13𝐿−22𝐿 4𝐿2 −13𝐿 −3𝐿2
5413𝐿
−13𝐿−3𝐿2
15622𝐿
22𝐿4𝐿2
].
4.3 Viga sob momentos torsores
Quando sobre a viga atuam os momentos torsores 𝑃4 e 𝑃10, e ocorrem os deslocamentos
angulares 𝛿4 e 𝛿10, assumindo que o centro de corte da secção transversal coincide com o
centróide desta, a solução da equação diferencial que governa este problema para uma viga
uniforme é:
𝜃 =
𝑇. 𝑥
𝐽. 𝐺+ 𝑐1 = 𝑐2𝑥 + 𝑐1, (4.11)
em que 𝑇 é o momento torsor na viga, 𝜃 é o deslocamento angular sobre o eixo 𝑥, e 𝑐1 e 𝑐2
são as constantes de integração. Como a Equação (4.11) é similar à Equação (4.1), neste
caso, as funções de forma serão iguais às que foram determinadas para os efeitos axiais:
𝑁4 = 1 −𝑥
𝐿
e
𝑁10 =𝑥
𝐿.
Por este motivo, os coeficientes de rigidez e de massa serão também semelhantes:
𝑘𝑖,𝑗 = ∫ 𝐽. 𝐺. 𝑁𝑖′𝑁𝑗
′𝑑𝑥𝐿
0
, 𝑖, 𝑗 = {4,10}
e
𝑚𝑖,𝑗 = ∫ 𝐼m𝑁𝑖𝑁𝑗𝑑𝑥, 𝑖, 𝑗 = {4,10},𝐿
0
em que 𝐼m = ��𝐼0/𝐴 é o segundo momento polar de massa. Assim, para as vigas uniformes
[Ke] = [𝑘4,4 𝑘4,10𝑘10,4 𝑘10,10
] =𝐽. 𝐺
𝐿|e
[1 −1−1 1
]
e
[Me] = [𝑚4,4 𝑚4,10
𝑚10,4 𝑚10,10] =
1
6
��𝐼0𝐿
𝐴|e
[2 11 2
].
4.4 Coeficientes de amortecimento
Os coeficientes de amortecimento também seriam obtidos de forma análoga aos
coeficientes de massa:
24
𝑐𝑖,𝑗 = ∫ 𝑐.𝑁𝑖𝑁𝑗𝑑𝑥, 𝑖, 𝑗 = {1,2,3, … ,12},𝐿
0
em que 𝑐 é o coeficiente de amortecimento por uma unidade de comprimento, que na
prática é difícil de obter isoladamente para cada elemento, como já foi referido
anteriormente.
4.5 Treliças 3D sob forças concentradas nos nós
As treliças são estruturas articuladas, sujeitas a forças aplicadas apenas nos nós (ou juntas)
dos elementos que as constituem, transmitindo, por conseguinte, apenas os esforços
normais. Por outro lado, as forças de inércia irão causar a flexão destes elementos, com
uma curvatura nula, sendo, por isso, caracterizados pelas mesmas funções de forma
determinadas para os efeitos axiais. No caso mais geral, numa treliça no espaço atuam as
forças 𝑃1, 𝑃2, 𝑃3, 𝑃7, 𝑃8 e 𝑃9, e ocorrem os deslocamentos 𝛿1, 𝛿2, 𝛿3, 𝛿7, 𝛿8 e 𝛿9. Neste
caso, os coeficientes da matriz de massa ficam como:
𝑚𝑖,𝑗 = ∫ ��𝑁𝑖𝑁𝑖𝑑𝑥, 𝑖, 𝑗 = {1,2,3,7,8,9}.𝐿
0
Para uma barra uniforme
𝑚𝑖,𝑗 = ��∫ 𝑁𝑖𝑁𝑗𝑑𝑥, 𝑖, 𝑗 = {1,2,3,7,8,9},𝐿
0
pelo que
[Me] =��. 𝐿
6|e
[ 2 0 00 2 00 0 2
1 0 00 1 00 0 1
1 0 00 1 00 0 1
2 0 00 2 00 0 2]
.
Como nestes elementos só existem os esforços normais, considerando que as deformações
são pequenas, para se ter os deslocamentos 𝛿2, 𝛿3, 𝛿8 ou 𝛿9, e como os nós dos elementos
deslocam-se livremente nessas direções, não são necessárias as forças 𝑃2, 𝑃3, 𝑃8 ou 𝑃9, sendo estas independentes das forças 𝑃1 e 𝑃7, logo, a matriz de rigidez fica como:
[Ke] =
[ 𝑘1,1 0 0
0 0 00 0 0
𝑘1,7 0 0
0 0 00 0 0
𝑘7,1 0 0
0 0 00 0 0
𝑘7,7 0 0
0 0 00 0 0]
=𝐴. 𝐸
𝐿|e
[ 1 0 0 0 0 0 0 0 0
−1 0 0 0 0 0 0 0 0
−1 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0]
.
25
4.5.1 Transformação de coordenadas
Os elementos que constituem as treliças podem ter o seu referencial local não-coincidente
com o referencial global, tal como se mostra na Figura 6. Por isso, as matrizes elementares
têm de ser transformadas em termos de coordenadas globais, para ser coerente assemblá-
-las nas matrizes globais. Sabe-se que
{𝑃1𝑃2𝑃3
} = [T1] {
��1��2��3
} {𝑃4𝑃5𝑃6
} = [T1] {
��4��5��6
}
e
{
𝛿1𝛿2𝛿3
} = [T1] {
𝛿1𝛿2𝛿3
} {
𝛿4𝛿5𝛿6
} = [T1] {
𝛿4𝛿5𝛿6
},
em que [T1] é a matriz de transformação das coordenadas:
[T1] = [
𝑐1 𝑐2 𝑐3𝑐4 𝑐5 𝑐6𝑐7 𝑐8 𝑐9
] = [
cos(∡(��, ��)) cos(∡(��, ��)) cos(∡(��, ��))
cos(∡(��, ��)) cos(∡(��, ��)) cos(∡(��, ��))
cos(∡(��, ��)) cos(∡(��, ��)) cos(∡(��, ��))
].
Os termos de [T1] são cossenos de ângulos entre os versores do referencial global e os
versores do referencial local. Sendo
[T] = ⌈[T1], [T1]⌋,
em que o operador ⌈ ⌋ indica que [T] é uma matriz diagonal constituída pelas submatrizes [T1], as forças e os deslocamentos podem ser representados de forma mais compacta:
{P} = [T]{P} (4.12)
e
{δ} = [T]{δ}. (4.13)
As forças elásticas que agem sobre um elemento nas coordenadas globais são:
{P} = [Ke]{δ} (4.14)
e nas coordenadas locais
{P} = [Ke]{δ}. (4.15)
Substituindo {P} e {δ} da Equação (4.12) e (4.13) na Equação (4.15), resulta:
26
[T]{P} = [Ke][T]{δ}. (4.16)
Como [T] é uma matriz ortogonal ([T]−1 = [T]T), multiplicando à esquerda a Equação
(4.16) pela transposta de [T]:
[T]T[T]{P} = [T]T[Ke][T]{δ}, (4.17)
tem-se
{P} = [T]T[Ke][T]{δ}, (4.18)
pelo que
[Ke] = [T]T[Ke][T]. (4.19)
Da mesma forma, encontra-se a matriz de massa de um elemento nas coordenadas globais:
[Me] = [T]T[Me][T]. (4.20)
A sub-rotina “Transf_Cord_3D.m” foi desenvolvida para permitir obter [T1], como
também o comprimento de cada elemento e a matriz de transformação do referencial
global para o local, tendo como os argumentos: as coordenadas de cada nó, os nós globais
de cada elemento, o tipo de estrutura e os nós auxiliares. No entanto, no caso das treliças
3D, só serão necessários os cossenos da primeira linha de [T1], ou seja, 𝑐1, 𝑐2 e 𝑐3, para
definir a matriz de rigidez elementar:
[Ke] =𝐴. 𝐸
𝐿|e
[ 𝑐1
2 𝑐1𝑐2 𝑐1𝑐3 𝑐1𝑐2 𝑐2
2 𝑐2𝑐3 𝑐1𝑐3 𝑐2𝑐3 𝑐3
2
−𝑐12 −𝑐1𝑐2 −𝑐1𝑐3
−𝑐1𝑐2 −𝑐22 −𝑐2𝑐3
−𝑐1𝑐3 −𝑐2𝑐3 −𝑐32
−𝑐12 −𝑐1𝑐2 −𝑐1𝑐3
−𝑐1𝑐2 −𝑐22 −𝑐2𝑐3
−𝑐1𝑐3 −𝑐2𝑐3 −𝑐32
𝑐12 𝑐1𝑐2 𝑐1𝑐3
𝑐1𝑐2 𝑐22 𝑐2𝑐3
𝑐1𝑐3 𝑐2𝑐3 𝑐32 ]
.
De forma ainda mais simples, a matriz de massa nas coordenas globais fica: [Me] = [Me]. As sub-rotinas “Rigidez_Glob_truss_3D.m” e “Mass_Consist_Glob_truss_3D.m”
permitem assemblar essas matrizes elementares nas matrizes de rigidez [K] e de massa [M] do sistema nas coordenas globais, com a possibilidade de nesta última se incluir as massas
concentradas nos nós, se estas existirem num dado problema.
4.6. Estruturas 3D constituídas por elementos do tipo viga-barra
Quando sobre os elementos do tipo viga-barra de uma estrutura reticulada 3D atuam todos
os esforços e ocorrem todos os deslocamentos representados na Figura 6, os coeficientes de
massa e de rigidez elementar, para este caso mais geral, são assemblados a partir dos
27
coeficientes que já foram sendo determinados isoladamente para os efeitos causados por
conjuntos específicos de esforços. Sabe-se que
{P1P2P3
} = [T1] {
P1P2P3
} {P4P5P6
} = [T1] {
P4P5P6
} {P7P8P9
} = [T1] {
P7P8P9
} {P10P11P12
} = [T1] {
P10P11P12
}
e
{
δ1δ2δ3
} = [T1] {
δ1δ2δ3
} {
δ4δ5δ6
} = [T1] {
δ4δ5δ6
} {
δ7δ8δ9
} = [T1] {
δ7δ8δ9
} {
δ10δ11δ12
} = [T1] {
δ10δ11δ12
}.
Como já foi visto na Equação (4.12) e (4.13), as forças e os deslocamentos podem ser
representados de forma mais compacta. Sendo a matriz
[T] = ⌈[T1], [T1], [T1], [T1]⌋,
em que [T1] é a submatriz que já foi definida anteriormente. As matrizes de rigidez e de
massa elementar são estabelecidas como:
[Ke] =
[ 𝑘1,1 0 0
0 𝑘2,2 0
0 0 𝑘3,3
0 0 00 0 𝑘2,60 𝑘3,5 0
0 0 00 0 𝑘5,30 𝑘6,2 0
𝑘4,4 0 0
0 𝑘5,5 0
0 0 𝑘6,6
𝑘1,7 0 0
0 𝑘2,8 0
0 0 𝑘3,9
0 0 00 0 𝑘2,120 𝑘3,11 0
0 0 00 0 𝑘5,90 𝑘6,8 0
𝑘4,10 0 0
0 𝑘5,11 0
0 0 𝑘6,12𝑘7,1 0 0
0 𝑘8,2 0
0 0 𝑘9,3
0 0 00 0 𝑘8,60 𝑘9,5 0
0 0 0 0 0 𝑘11,3 0 𝑘12,2 0
𝑘10,4 0 0
0 𝑘11,5 0
0 0 𝑘12,6
𝑘7,7 0 0
0 𝑘8,8 0
0 0 𝑘9,9
0 0 00 0 𝑘8,120 𝑘9,11 0
0 0 00 0 𝑘11,90 𝑘12,8 0
𝑘10,10 0 0
0 𝑘11,11 0
0 0 𝑘12,12]
e
[Me] =
[ 𝑚1,1 0 0
0 𝑚2,2 0
0 0 𝑚3,3
0 0 00 0 𝑚2,6
0 𝑚3,5 0
0 0 00 0 𝑚5,3
0 𝑚6,2 0
𝑚4,4 0 0
0 𝑚5,5 0
0 0 𝑚6,6
𝑚1,7 0 0
0 𝑚2,8 0
0 0 𝑚3,9
0 0 00 0 𝑚2,12
0 𝑚3,11 0
0 0 00 0 𝑚5,9
0 𝑚6,8 0
𝑚4,10 0 0
0 𝑚5,11 0
0 0 𝑚6,12
𝑚7,1 0 0
0 𝑚8,2 0
0 0 𝑚9,3
0 0 00 0 𝑚8,6
0 𝑚9,5 0
0 0 0 0 0 𝑚11,3
0 𝑚12,2 0
𝑚10,4 0 0
0 𝑚11,5 0
0 0 𝑚12,6
𝑚7,7 0 0
0 𝑚8,8 0
0 0 𝑚9,9
0 0 00 0 𝑚8,12
0 𝑚9,11 0
0 0 00 0 𝑚11,9
0 𝑚12,8 0
𝑚10,10 0 0
0 𝑚11,11 0
0 0 𝑚12,12]
.
As matrizes [Ke] e [Me] são obtidas diretamente através das Equações (4.19) e (4.20) nas
sub-rotinas “Rigidez_Glob_bbeam_3D.m” e “Mass_Consist_Glob_bbeam_3D.m”, como
também, estas permitem assemblá-las nas correspondentes matrizes de rigidez [K] e de
28
massa [M] do sistema nas coordenas globais, com a possibilidade de nesta última incluir as
massas concentradas nos nós, se estas existirem num dado problema.
29
5 MODELAÇÃO DO PROGRAMA PRINCIPAL EM MATLAB
O programa principal deste trabalho “estrutura3D_com_GiD.m” foi feito para ser possível
resolver problemas de vibração de estruturas reticuladas e articuladas, com implementação
numérica de elementos finitos em ambiente de programação MATLAB, funcionando não
só com as funções auxiliares que já foram introduzidas nos capítulos anteriores mas
também com outras funções necessárias para fazer a ponte com o GiD e que serão
apresentadas ao longo do texto que se segue. Estas últimas foram baseadas na aplicação já
existente para problemas estáticos com elementos finitos 3D (tridimensionais) (ver a
referência (4)).
Este programa foi elaborado de modo a ser sistemático durante a sua execução, i.e., ao
utilizador são apresentadas automática e sequencialmente janelas com perguntas por
responder e tabelas por preencher, relativamente à estrutura de um dado problema e à sua
análise. Desta forma, torna-se mais fácil de compreendê-lo, ao guiar o utilizador na
definição das variáveis estritamente necessárias em cada tipo de problema. Contudo, este
programa poderá ser futuramente alterado, otimizado ou adaptado, sendo, por isso, o
objetivo deste capítulo explicar a sua constituição. Além disso, o seu fluxograma foi
adicionado no Apêndice “A”.
O programa começa por adicionar o caminho relativo dos programas auxiliares (sub-
-rotinas), estando estes, neste caso, dentro da pasta funcoes_auxiliares, que está na mesma
pasta do programa principal. O número de nós por elemento é fixado pela constante
Nnos_el igual a 2. Na linha 5 é pedido para se inserir o nome do ficheiro que foi exportado
do GiD. A leitura deste ficheiro é realizada pela função “read_GiD_file.m”, extraindo as
coordenadas de cada nó, que são atribuídas à variável PP_nos e os nós de cada elemento à
matriz de conectividade nos_Stm. O número de elementos N_elem e os nós ordenados do
sistema nos_u obtêm-se a partir desta matriz de conectividade, enquanto o número total de
nós Nnos_t é obtido através da variável anterior PP_nos. Na linha 11 é escolhido o Sistema
de Unidades em função do qual nas linhas 15-36 definem-se as unidades das grandezas do
problema e com as quais se mostram posteriormente os resultados da análise. Na linha 13 é
pedido para se inserir os nós que tenham pelo menos um GL livre, ficando estes guardados
na variável no_l e na variável Nno_l o número desses nós. Na linha 37 é definida a variável
Estrutura, que determina se a estrutura é articulada ou reticulada. Em função desta última,
nas linhas 39-55 a constante Ng_no assume um valor predefinido igual ao número total dos
GL por nó dessa estrutura, os vetores constantes ind1 e ind2 indicam os GL locais
possíveis de cada nó do elemento. É também estabelecido um vetor com as funções de
forma paramétricas Nx_f_3D e a variável Nnos toma um valor igual ao número de nós
apenas do sistema, ou seja, sem os nós auxiliares. São também iniciados os vetores dos nós
auxiliares nos_aux e os vetores das propriedades mecânicas de cada elemento, que no caso
geral são: a área de secção transversal AA, o módulo de elasticidade EE, a densidade linear
mm, o segundo momento de área em relação ao eixo local x e y do elemento, II_x e II_y,
respetivamente, a constante de torção JJ e o módulo de corte GG. Para agilizar o passo
anterior, nas linhas 56-66 é atribuído o tipo de secção a cada elemento da estrutura. Tal
como se pode ver nas linhas 67-97, se a estrutura for articulada serão definidas apenas as
três primeiras propriedades mencionadas anteriormente, por outro lado serão definidas
todas aquelas, como também a variável nos_aux, no caso de ser uma estrutura reticulada,
embora em alguns problemas nem todas as propriedades mecânicas terão influência no
comportamento da estrutura, dependendo isto das condições de fronteira e das condições
30
iniciais do sistema. De seguida, é mostrado o trecho do código que foi referido até agora,
que são as linhas 1-97 do programa:
close all; clear; clc; addpath('funcoes_auxiliares/');
%% Dados do Sistema
Nnos_el=2;
dad='Introduza o nome do ficheiro com os dados do GiD : ';
str=inputdlg(dad,['Dados do Sistema: ' 'nome do ficheiro '],[1 80]);
nome_ficheiro=str{1};
[PP_nos, nos_Stm]=read_GiD_file(nome_ficheiro); nos_Stm=sort(nos_Stm,2);
N_elem=size(nos_Stm,1);
nos_u=unique(nos_Stm(:))'; % nós do sistema sem os auxiliares
Nnos_t=size(PP_nos,2); % número total de nós (com os auxiliares)
Unidades=questdlg('Que unidades pretende usar?','Unidades','Imperiais','Internacionais','Imperiais');
dad={'Nós do sistema com pelo menos um GL livre:'};
str=inputdlg(dad,'Dados do Sistema: nós com GL livres',[1 65]); no_l=str2num(str{1});
Nno_l=length(no_l);
switch Unidades
case 'Imperiais'
unid.comp='[in]'; % unidade de comprimento
unid.area='[in^2]'; % unidade de área
unid.forca='[lb]'; % unidade de força
unid.mom='[lb.in]'; % unidade de momento
unid.massa='[lb.s^2/in]'; % unidade de massa
unid.densid_lin='[lb.s^2/in^2]'; % unidade de densidade linear
unid.tempo='[s]'; % unidade de tempo
unid.mom_a='[in^4]'; % unidade de segundo momento de área
unid.young='[lb/in^2]'; % unidade de módulo de Young
case 'Internacionais'
unid.comp='[mm]'; % unidade de comprimento
unid.area='[mm^2]'; % unidade de área
unid.forca='[N]'; % unidade de força
unid.mom='[N.mm]'; % unidade de momento
unid.massa='[N.s^2/mm]'; % unidade de massa
unid.densid_lin='[N.s^2/mm^2]'; % unidade de densidade linear
unid.tempo='[s]'; % unidade de tempo
unid.mom_a='[mm^4]'; % unidade de segundo momento de área
unid.young='[N/mm^2]'; % unidade de módulo de Young
end
Estrutura=questdlg('A estrutura é articulada ou reticulada?','Tipo de estrutura', ...
'articulada','reticulada','articulada');
switch Estrutura
case 'articulada'
Ngb_no=3; AA=zeros(1,N_elem); EE=AA; mm=AA;
ind1=1:3; ind2=4:6;
Nx_f_3D=@(x,Le)[1-(x/Le); 1-(x/Le);1-(x/Le);(x/Le);(x/Le);(x/Le)];
Nnos=Nnos_t; nos_aux=[]; def=repmat({'0'},1,3);
case 'reticulada'
Ngb_no=6; def=repmat({'0'},1,10);
AA=zeros(1,N_elem); EE=AA; mm=AA; II_y=AA; II_z=AA; IIm=AA; JJ=AA; GG=AA;
max_y=AA; max_z=AA;
nos_aux=zeros(N_elem,1);
ind1=1:6; ind2=7:12;
Nx_f_3D=@(x,Le)[1-x/Le; (2*x.^3)/Le^3-(3*x.^2)/Le^2+1; (2*x.^3)/Le^3-(3*x.^2)/Le^2+1;...
1-x/Le; -x.*(x/Le-1).^2; x-(2*x.^2)/Le+x.^3/Le^2; x/Le; (3*x.^2)/Le^2-(2*x.^3)/Le^3;...
(3*x.^2)/Le^2-(2*x.^3)/Le^3; x/Le; x.^2/Le-x.^3/Le^2; x.^3/Le^2-x.^2/Le];
Nnos=length(nos_u); % número de nós sem os auxiliares
31
end
dad={'Quantas secções transversais com propriedades mecânicas diferentes existem?'};
str=inputdlg(dad,'Dados do Sistema: secções dos elementos ', [1 80] ); N_sec=str2double(str(1));
Elem_sec=cell(1,N_sec); dad=Elem_sec;
if N_sec~=1
for p=1:N_sec
dad(p)={['A secção n.º ' num2str(p) ' corresponde aos elementos: ']};
end
str=inputdlg(dad,'Dados do Sistema: secções dos elementos', [1 80] ); Elem_sec=str;
else
Elem_sec{1}=num2str(1:N_elem);
end
for p=1:N_sec
for elem=str2num(Elem_sec{p})
if Estrutura=="articulada"
dad={['Área de secção transversal ' unid.area ':'],...
['Módulo de Young ' unid.young ':'],...
['Densidade linear ' unid.densid_lin ':']};
str=inputdlg(dad,['Dados do Sistema: ' 'elemento ' num2str(elem) ...
' (com a secção n.º ' num2str(p) ')' ],[1 80],def);
AA(elem)=str2double(str(1)); EE(elem)=str2double(str(2)); mm(elem)=str2double(str(3));
def(:)=str(:);
elseif Estrutura=="reticulada"
dad={['Área de secção transversal ' unid.area ':'], ...
['Módulo de Young ' unid.young ':'], ...
['Densidade linear ' unid.densid_lin ':'], ...
['Segundo momento de área em relaçao ao eixo local y ' unid.mom_a ':'],...
['Segundo momento de área em relaçao ao eixo local z ' unid.mom_a ':'],...
['Constante de torção ' unid.mom_a ':'], ...
['Módulo de corte ' unid.young ':'], ...
['Posição máxima absoluta de um ponto da secção no eixo local y ' unid.comp ':'], ...
['Posição máxima absoluta de um ponto da secção no eixo local z ' unid.comp ':'], ...
'Um nó coincidente com o plano local Oxy:'};
str=inputdlg(dad,['Dados do Sistema: ' 'elemento ' num2str(elem) ...
' (com a secção n.º ' num2str(p) ')' ],[1 80],def);
AA(elem)=str2double(str(1)); EE(elem)=str2double(str(2)); mm(elem)=str2double(str(3));
II_y(elem)=str2double(str(4)); II_z(elem)=str2double(str(5)); JJ(elem)=str2double(str(6));
GG(elem)=str2double(str(7)); max_y(elem)=str2double(str(8)); max_z(elem)=str2double(str(9));
nos_aux(elem)=str2double(str(10)); def{10}=num2str(nos_aux(elem));
IIm(elem)=mm(elem)*(II_y(elem)+II_z(elem))/AA(elem); def(:)=str(:);
end
end
end
Na linha 99 é pedido ao utilizador para escolher se a matriz de massa de cada elemento é
Consistent, tal como já foi definida no capítulo anterior, ou Lumped, caso em que a inércia
rotacional é desprezada e a inércia translacional agregada em cada nó é desacoplada das
restantes. Posteriormente, nas linhas 101-117 introduzem-se os GL locais livres de cada
nó, que são guardados na variável gl_l_no. Estes GL podem ter um valor inteiro de 1 até ao
valor máximo da variável Ngb_no e, no caso geral, quando esta variável é igual a 6,
representam os deslocamentos lineares no eixo global X, Y e Z, e as 3 rotações nesses
mesmos eixos, respetivamente. Em função da variável Inercia definida anteriormente, são
deduzidos os GL livres com e sem a inércia nas variáveis gl_l_no_m e gl_l_no_sm,
respetivamente, como também o número de MV na variável Nmodos.
32
De seguida, na linha 118 são definidos os nós com uma massa concentrada na variável
nos_Mc e o seu valor é associado ao vetor Mc. Nas linhas 121-138, a partir dos dados que
foram especificados anteriormente relativamente ao sistema, são definidas outras variáveis,
tais como: o número total de GL livres Ngl_l, o número de GL por elemento Ngl_el e o
número de GL do sistema Ngl_Stm. São também criadas as matrizes auxiliares, sendo:
gl_no1 a matriz com os GL locais de cada nó, por linha, de um elemento apenas; gl_no3 a
matriz com os GL de cada nó, por linha, do sistema; gl_Stm a matriz com os GL de cada
elemento, por linha, do sistema. A variável gl_l é definida como um vetor com todos os
GL livres do sistema e gl_Mc como um vetor que contém os GL do sistema relativos aos 3
GL dos deslocamentos lineares de cada nó com uma massa concentrada. Assim, as linhas
98-138 correspondem à seguinte fração código:
nn=1; gl_l_no=cell(1,Nno_l); gl_l_no_m=gl_l_no; gl_l_no_sm=gl_l_no; def={''};
Inercia=questdlg('A matriz de massa de cada elemento da estrutura é:','Inércia','Lumped',...
'Consistent','Lumped');
switch Inercia
case 'Consistent'
for no=no_l
dad={['GL livres (de 1 a ' num2str(Ngb_no) ') do nó ' num2str(no) ':']};
str=inputdlg(dad,['Dados do Sistema: GL livres do nó ' num2str(no)],[1 65],def);
gl_l_no{nn}=str2num(str{1}); def=str; bb=gl_l_no{nn}; gl_l_no_m{nn}=gl_l_no{nn}; nn=nn+1;
end
Nmodos=length([gl_l_no_m{:}]);
case 'Lumped'
for no=no_l
dad={['GL livres (de 1 a ' num2str(Ngb_no) ') do nó ' num2str(no) ':']};
str=inputdlg(dad,['Dados do Sistema: GL livres do nó ' num2str(no)],[1 65],def);
gl_l_no{nn}=str2num(str{1}); def=str; bb=gl_l_no{nn}; gl_l_no_m{nn}=bb(gl_l_no{nn}<4);
gl_l_no_sm{nn}=bb(~(gl_l_no{nn}<4)); nn=nn+1;
end
Nmodos=length([gl_l_no_m{:}]);
end
dad={'Nós do sistema com massa concentrada:', ['Massa concentrada nesses nós ' unid.massa ':']};
str=inputdlg(dad,'Dados do Sistema: massa concentrada',[1 65]); nos_Mc=str2num(str{1});
Mc=str2num(str{2});
Ngl_l=length([gl_l_no{:}]); Ngl_l_m=Nmodos; Ngl_l_sm=Ngl_l-Ngl_l_m; Ngl_el=Nnos_el*Ngb_no;
Ngl_Stm=Nnos*Ngb_no; gl_no1=zeros(Nnos_el,Ngb_no)'; gl_no1(1:Ngl_el)=1:Ngl_el; gl_no1=gl_no1';
gl_Stm=zeros(N_elem,Ngl_el); gl_no2=zeros(Nnos,Ngb_no)'; gl_no2(1:Ngl_Stm)=1:Ngl_Stm;
gl_no2=gl_no2'; gl_no3=zeros(Nnos_t,Ngb_no); gl_no3(nos_u,:)=gl_no2; gl_l=zeros(1,Ngl_l);
gl_l_m=zeros(1,Ngl_l_m); gl_l_sm=zeros(1,Ngl_l_sm); n0=1; n00=1; n000=1;
for no=1:Nno_l
nf=n0+length(gl_l_no{no})-1; nff=n00+length(gl_l_no_m{no})-1;
nfff=n000+length(gl_l_no_sm{no})-1;
gl_l(1,n0:nf)=gl_no3(no_l(no),gl_l_no{no}); gl_l_m(1,n00:nff)=gl_no3(no_l(no),gl_l_no_m{no});
gl_l_sm(1,n000:nfff)=gl_no3(no_l(no),gl_l_no_sm{no});
n0=nf+1; n00=nff+1; n000=nfff+1;
end
for elem=1:N_elem
for no_e=1:Nnos_el
gl_Stm(elem,gl_no1(no_e,:))=gl_no3(nos_Stm(elem,no_e),:);
end
end
gl_Mc_aux=gl_no3(nos_Mc,:); gl_Mc=gl_Mc_aux(:,1:3);
33
Na próxima secção, é feito o cálculo numérico das frequências naturais w_n e dos MV
normalizados Phi do sistema, a partir das matrizes de rigidez K_gl e de massa M_gl da
estrutura, assembladas em função do tipo desta e da variável Inercia, da rotação dos
elementos, definidos na hipermatriz T1, do comprimento dos elementos LL e de outras
variáveis definidas anteriormente. Caso haja GL sem inércia associada recorre-se à
condensação estática para diminuir o esforço computacional nesse cálculo, tal como se
pode ver nas linhas 157-174. A hipermatriz de transformação dos elementos do referencial
global para o local TT1_inv será útil mais adiante, tal como se verá.
%% Cálculo das frequências naturais e dos MV normalizados do Sistema
[T1, LL, TT1_inv]=Transf_Cord_3D(PP_nos,nos_Stm,Estrutura,nos_aux);
if Estrutura=="articulada"
[K_gl,K_l]=Rigidez_Glob_truss_3D(LL,EE,T1,AA,gl_Stm,Ngl_Stm);
if Inercia=="Consistent"
[M_gl]=Mass_Consist_Glob_truss_3D(LL,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
elseif Inercia=="Lumped"
[M_gl]=Mass_Lump_Glob_truss_3D(LL,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
end
elseif Estrutura=="reticulada"
[K_gl,K_l]=Rigidez_Glob_bbeam_3D(LL,T1,EE,GG,JJ,II_y,II_z,AA,gl_Stm,Ngl_Stm);
if Inercia=="Consistent"
[M_gl]=Mass_Consist_Glob_bbeam_3D(LL,T1,mm,IIm,AA,gl_Stm,Ngl_Stm,gl_Mc,Mc);
elseif Inercia=="Lumped"
[M_gl]=Mass_Lump_Glob_bbeam_3D(LL,T1,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
end
end
K_aa=K_gl(gl_l_m,gl_l_m); Maa=M_gl(gl_l_m,gl_l_m);
if Ngl_l_sm~=0
K_ac=K_gl(gl_l_m,gl_l_sm); K_ca=K_gl(gl_l_sm,gl_l_m); K_cc=K_gl(gl_l_sm,gl_l_sm);
K_A=K_aa-K_ac/K_cc*K_ca;
[phi_a,D]=eig(K_A,Maa); phi_c=-K_cc\K_ca*phi_a; phi=zeros(Ngl_l,Ngl_l_m); aa=1; bb=1; nn=1;
for gl=gl_l
if gl==gl_l_m(aa)
phi(nn,:)=phi_a(aa,:);
if ~(gl==gl_l_m(end))
aa=aa+1;
end
elseif gl==gl_l_sm(bb)
phi(nn,:)=phi_c(bb,:);
if ~(gl==gl_l_sm(end))
bb=bb+1;
end
end
nn=nn+1;
end
elseif Ngl_l_sm==0
[phi,D]=eig(K_aa,Maa);
end
[w_n, I]=sort(diag(sqrt(D))); phi=phi(:,I); T=2*pi./w_n; f_n=1./T; Phi=zeros(Ngl_l,Nmodos);
for i=1:Ngl_l
for j=1:Nmodos
Phi(i,j)=phi(i,j)/sqrt(phi(:,j)'*M_gl(gl_l,gl_l)*phi(:,j));
end
end
34
Na seguinte secção, são definidos os dados da análise, que dependem essencialmente do
tipo de problema e da forma como se deseja obter e visualizar os resultados. Na primeira
análise do problema o programa pressupõe que no instante inicial o sistema está em
repouso, i.e., os deslocamentos U0 e as velocidades iniciais DU0 de cada GL são nulos. A
partir da linha 186 o programa passa a ser iterativo, ou seja, no final do programa é dada a
possibilidade ao utilizador de corrigir ou alterar os dados da análise e depois voltar a ver ou
exportar os resultados com essa alteração. O amortecimento de cada MV é definido no
vetor eps. Nas linhas 195-202 são definidas as variáveis em relação ao tempo total da
análise tf, ao intervalo de tempo dt, ao vetor de tempo discretizado t e ao número de
sinusoides da série de Fourier Ns para aproximar as forças modais equivalentes a um
conjunto de sinusoides sobrepostos. A variável dx é o intervalo de comprimento que é
aplicado aos elementos na sua discretização, calculada automaticamente para tomar um
valor mínimo e ser ainda possível importar os resultados na versão estudante do GiD, com
um número máximo de nós limitado a 1010. Na linha 204 é escolhido ambiente em que se
deseja ver os resultados, ficando esta opção armazenada na variável pos_proc. É possível
observar nas linhas 205-212 que se for escolhido o pós-processador do GiD para ver os
resultados, o programa acaba logo que dois ficheiros de saída forem criados, pois quer a
variável Ampl_ok quer a variável Analise_ok estarão verdadeiras e determinarão o fim dos
dois ciclos repetitivos. Por outro lado, se for escolhido o ambiente do MATLAB, além de
ser necessário definir a variável de amplificação Ampl_n, essas variáveis estarão falsas
nesse momento e terão, por isso, que ser validadas mais tarde, para ser possível mudar os
dados relativamente à análise, se for assim decidido.
Nas linhas 213-260 é estabelecida a matriz F, que representa o esforço que é aplicado nos
GL livres, e, portanto, sob ou sem esforço em cada instante do tempo da análise. Esta
matriz é obtida de forma individual para cada GL livre sob esforço, sendo possível definir
para cada um destes um tipo de esforço diferente: constante com uma certa magnitude;
sinusoidal com uma magnitude, uma fase e uma frequência angular específica; ou não-
-uniforme, obtido de um ficheiro externo, onde em cada linha o primeiro termo é o instante
de tempo e o segundo termo, separado do primeiro pela vírgula ou espaço, é a magnitude
do esforço nesse instante, que são armazenados na variável t_esf e interpolados
linearmente para os instantes definidos pela variável do tempo t. Esta secção do algoritmo
encontra-se representada na seguinte lista:
%% Dados da Análise
U0=zeros(Ngl_l,1); DU0=U0; Analise_ok=false;
while Analise_ok==false
dad={'O amortecimento presente (ou não) é igual para a maior parte dos MV a: ',...
['Os modos de vibração (de 1 a ' num2str(Nmodos) ') com um amortecimento diferente são: ']};
str=inputdlg(dad,'Dados da Análise: Amortecimento', [1 80] );
eps=str2double(str(1))*ones(Nmodos,1); Nmv_am2=str2num(str{2}); n0=1;
for n_m=Nmv_am2
dad={['O amortecimento presente (ou não) no MV ' num2str(Nmv_am2(n0)) ' é igual a:' ]}; n0=n0+1;
str=inputdlg(dad,'Dados da Análise: Amortecimento', [1 80] ); eps(n_m)=str2double(str(1));
end
dad={['Analisar a resposta do sistema desde 0 até ' unid.tempo ':'], ...
['Com um intervalo de tempo discretizado (dt) de ' unid.tempo ':'], ....
'Com um número de sinusoides (Ns) da série de Fourier igual a:'};
str=inputdlg(dad,'Dados da Análise', [1 65] );
t0=0; tf=str2double(str(1));
dt=str2double(str(2));
t=t0:dt:tf; lgt_t=length(t);
35
Ns=str2double(str(3));
dx=sum(LL)/(1010-1);
pos_proc=questdlg('Onde deseja ver os resultados?','Pós-proces.','No MATLAB','No GiD','No MATLAB');
switch pos_proc
case 'No MATLAB'
dad={'Usar uma razão de amplificação da deformada do sistema de:'};
str=inputdlg(dad,'Dados da Análise: deformada do sistema', [1 70] );
Ampl_n=[1 str2double(str(1))]; n_c=2; Ampl_ok=false;
case 'No GiD'
Ampl_n=1; n_c=1; Ampl_ok=true; Analise_ok=true;
end
dad='Quais são os nós livres do sistema que estão sob algum esforço externo?';
str=inputdlg(dad,'Dados da Análise: nós livres sob esforço',[1 75]); no_F=str2num(str{1});
no_F=sort(no_F); nn=1; nn2=1; gl_F=cell(1,Nno_l); F=zeros(Ngl_l,lgt_t); m00=0;
for no=no_l
if nn2<=length(no_F)
if no_F(nn2)==no
dad=['Entre os GL locais livres (' num2str(gl_l_no{nn}) ') do nó ' num2str(no) ...
', os que estão sob um esforço são:' ];
str=inputdlg(dad,'Dados da Análise: GL livres sob esforço',[1 70]);
gl_F{nn}=str2num(str{1});
for gl=gl_F{nn}
Esforco=questdlg(['O esforço sobre o GL local ' num2str(gl) ' do nó ' num2str(no) ' é:'],...
'Dados da Análise: tipo de esforço','Constante','Sinusoidal',...
'Não-uniforme (interpolado linearmente)','Constante');
dad2=['Dados da Análise: Esforço no GL ' num2str(gl_no3(no,gl)) ' do Sistema'];
m0=m00+find(gl==gl_l_no{nn});
switch Esforco
case 'Constante'
dad=['A magnitude do esforço constante ' unid.forca ' ou ' unid.mom ' no GL local '...
num2str(gl) ' do nó ' num2str(no) ' é:'];
str=inputdlg(dad,dad2,[1 80]); magn=str2double(str(1)); F(m0,:)=magn*ones(1,lgt_t);
case 'Sinusoidal'
dad={['A magnitude do esforço sinusoidal ' unid.forca ' ou ' unid.mom ...
' no GL local ' num2str(gl) ' do nó ' num2str(no) ' é:'], ...
'Com uma fase de:', 'E com uma frequência angular de:'};
str=inputdlg(dad,dad2,[1 80]);
amplt=str2double(str(1)); fase=str2double(str(2)); w_a=str2double(str(3));
F(m0,:)=amplt*sin(w_a*t+fase);
case 'Não-uniforme (interpolado linearmente)'
dad=['Introduza o nome do ficheiro com o histórico de carregamento no GL local ' ...
num2str(gl) ' do nó ' num2str(no) ' (instante ' unid.tempo ', esforço [ ' ...
unid.forca ' ou ' unid.mom ' ] (por linha) ) :' ];
str=inputdlg(dad,dad2,[1 145]); t_esf=load(str{1}); N_esf_l=size(t_esf,1); n0=1;
for n=1:(N_esf_l-1)
tt=t_esf(n,1):dt:(t_esf(n+1,1)-dt); lgt_tt=length(tt);
if lgt_tt>1
nf=n0+length(tt)-1;
m_b=polyfit([tt(1), tt(end)], [t_esf(n,2), t_esf(n+1,2)],1); m=m_b(1); b=m_b(2);
F(m0,n0:nf)=m*tt+b; n0=nf+1;
end
end
end
end
nn2=nn2+1;
end
end
m00=m00+length(gl_l_no{nn}); nn=nn+1;
36
end
Desta forma, o pré-processamento termina na linha 260. Já na linha 262 é calculado o fator
de participação yi de cada MV do sistema, em função das variáveis definidas
anteriormente, com o qual é determinada a resposta nos GL livres do sistema, definida na
variável U, através da sobreposição modal nas linhas 264-266. Caso seja escolhido ver os
resultados no MATLAB, neste momento, são representadas as figuras para cada nó, como
se pode constatar nas linhas 269-293, em cada uma das quais aparecem os gráficos de
resposta para cada GL livre desse nó, tal como se pode ver nas linhas 261-294 do código,
na seguinte lista:
%% Cálculo e visualização da resposta do sistema nos seus GL
[yi] =fator_part_v_amort(M_gl(gl_l,gl_l),w_n,Phi,eps,F,Ns,t,dt,U0,DU0);
U=zeros(Ngl_l,lgt_t);
for k=1:Nmodos
U=yi(k,:).*Phi(:,k)+U; % resposta vista dos GL do sistema
end
if n_c==2
no1=1; n0=0; m0=0;
for no=no_l
f_no=figure; set(f_no,'Name',['Deslocamentos sem amplificação do nó ' num2str(no) ], ...
'NumberTitle','off'); nn=1;
for gl=gl_l_no{no1}
n0=m0+nn;
subplot(2,3,nn);
plot(t,U(n0,:));
grid on; xlabel(['t ' unid.tempo]);
if gl==1
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em X']);
elseif gl==2
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em Y']);
elseif gl==3
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em Z']);
elseif gl==4
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo X']);
elseif gl==5
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo Y']);
elseif gl==6
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo Z']);
end
nn=nn+1;
end
m0=m0+length(gl_l_no{no1}); no1=no1+1;
end
end
Na próxima parte do programa, o principal objetivo será o de obter a deformada da
estrutura em função da resposta nos GL do sistema, armazenada anteriormente na matriz
U. Esta é amplificada pela variável Ampl_n na linha 298, se o ciclo for da linha 296 ocorrer
pela segunda vez, para depois ser mais fácil detetar a deformada da estrutura no MATLAB.
Interessa, por isso, saber a posição de cada ponto de cada elemento discretizado entre os
seus nós. Este objetivo é executado em várias etapas. É criada a hipermatriz U_e com os
deslocamentos de cada elemento, vistos nos seus referenciais locais, através da
multiplicação da matriz de transformação T1_no pelos deslocamentos dos GL dos nós
37
desse elemento, vistos no referencial global e definidos na hipermatriz Ustm para esse
efeito. São também calculados esforços internos P1 que atuam em cada elemento por causa
da deformação elástica. As linhas 295-333 ilustram esta parte do código:
cic1=true;
for n=1:n_c
while Ampl_ok==false || cic1==true
U_ampl=Ampl_n(n)*U;
%% Alguns dados "rescritos" noutra forma
Ustm=zeros(Ngb_no,lgt_t,Nnos_t);
U_e=zeros(Ngl_el,lgt_t,N_elem);
i1=1;
for no=1:length(no_l)
i2=i1+length(gl_l_no{no})-1;
Ustm(gl_l_no{no},:,no_l(no))=U_ampl(i1:i2,:);
i1=i2+1;
end
for elm=1:N_elem
if Estrutura=="articulada"
T1no=T1(:,:,elm);
elseif Estrutura=="reticulada"
T1no=zeros(Ngb_no); T1no(1:3,1:3)=T1(:,:,elm); T1no(4:6,4:6)=T1(:,:,elm);
end
for no=no_l
no_e=1;
for nos_elm=nos_Stm(elm,:)
if no==nos_elm
if no_e==1
U_e(ind1,:,elm)=T1no*Ustm(:,:,nos_elm);
elseif no_e==2
U_e(ind2,:,elm)=T1no*Ustm(:,:,nos_elm);
end
end
no_e=no_e+1;
end
end
end
%% Esforços elementares elásticos nos seus GL
if n==1
P1=U_e;
for elm=1:N_elem
P1(:,:,elm)=K_l(:,:,elm)*U_e(:,:,elm);
end
Na seguinte secção, são inicializadas algumas variáveis não só para um melhor
desempenho do código mas também para definir algumas variáveis, tais como a posição
inicial (Xe0, Ye0, Ze0), a rotação inicial (Rot_ye0, Rot_ze0) de cada nó de cada elemento,
vistas nos seus referenciais locais. As funções de forma são também discretizadas para
cada GL e definidas para cada elemento na variável Nx. É também inicializada a variável
sigm_x_max, que é a tensão axial máxima em absoluto em cada secção transversal de cada
elemento. Esta parte do código encontra-se nas linhas 334-365:
%% Inicialização de alguns dados e resultados
Np=zeros(N_elem,1);
Nx=repmat({zeros(Ngl_el,1)},N_elem,1);
38
X_Y_Z_g=cell(N_elem,1); Xe=X_Y_Z_g; Ye=Xe; ums=Xe;
Xe0=zeros(N_elem,Nnos_el); Ye0=Xe0; Ze0=Xe0;
Rot_ye0=Xe0; Rot_ze0=Xe0; Xe0(:,2)=LL;
my=X_Y_Z_g; mz=X_Y_Z_g;
sigm_x_1=zeros(N_elem,lgt_t); sigm_x_max=X_Y_Z_g; sum1=0; nf=0;
for elm=1:N_elem
Le=LL(elm);
if Estrutura=="articulada"
n_xx=2;
elseif Estrutura=="reticulada"
n_xx=floor(Le/dx);
end
xx=linspace(0,Le,n_xx);
Np(elm)=length(xx);
ums{elm}=ones(1,Np(elm));
Xe{elm}=xx;
Ye{elm}=zeros(1,Np(elm));
X_Y_Z_g{elm}=zeros(4,Np(elm));
Nx{elm}=Nx_f_3D(xx,Le);
my{elm}=zeros(Np(elm)-1,1); mz{elm}=my{elm};
sigm_x_max{elm}=zeros(Np(elm)-1,lgt_t);
nf=Np(elm)+nf-1; n0=nf-Np(elm)+2;
Faces(n0:nf,1:2)=[(1+sum1):((Np(elm)-1)+sum1); (2+sum1):(Np(elm)+sum1)]';
sum1=Np(elm)+sum1;
end
Nt_p=sum(Np);
sigm_x_max_T=zeros(Nt_p+2*N_elem,lgt_t); cor=sigm_x_max_T;
Ze=Ye; ind=Ye; P_abs=zeros(3,Nt_p,lgt_t); N_cor=9; C=jet(N_cor); Cor=zeros(Nt_p,3,lgt_t);
min_sigm=zeros(1,lgt_t); max_sigm=min_sigm;
end
Desta forma, procede-se ao cálculo da tensão axial e da deformada da estrutura. Esta
última é calculada para cada elemento através da sobreposição da deformada causada por
deslocamentos individuais de cada GL desse elemento (U_e). As coordenadas de cada
ponto do elemento discretizado, definidas nas variáveis Xe, Ye e Ze, obtêm-se, assim, em
primeiro lugar, vistas no seu referencial local, multiplicando, em cada instante, as funções
de forma pelas correspondentes coordenadas e rotações dos nós, tal como se pode ver nas
linhas 371-373 para as estruturas articuladas e nas linhas 376-380 para as estruturas
reticuladas. Em seguida, é feita a multiplicação dessas coordenadas de cada ponto do
elemento, dispostas numa matriz homogénea, pela hipermatriz de transformação
geométrica TT1_inv, ficando, assim, expressas em relação ao referencial global, nas linhas
374 e 381. Estas últimas coordenadas dos pontos foram armazenadas na hipermatriz P_abs
e sob outra forma na Vert. Nas linhas 386-403 é calculada a tensão axial referida
anteriormente, em função do esforço axial interno equivalente e da área de cada elemento,
somado com a tensão axial resultante da flexão, em função da curvatura em cada ponto e
do módulo de Young de cada elemento.
Em seguida, nas linhas 409-430 são atribuídas as cores a cada elemento subdividido em
cada instante na variável Cor, em função da tensão calculada anteriormente, apresentando
uma cor azul (mais fria) onde esta tensão é mínima até uma cor vermelha (mais quente)
onde esta é máxima. Esta parte do programa está nas linhas 366-431:
%% Cálculo da deformada do sistema e da tensão axial
for p=1:lgt_t
39
sum1=0; sum2=0;
for elm=1:N_elem
if Estrutura=="articulada"
Xe{elm}=Nx{elm}(1,:)*(U_e(1,p,elm)+Xe0(elm,1))+Nx{elm}(4,:)*(U_e(4,p,elm)+Xe0(elm,2));
Ye{elm}=Nx{elm}(2,:)*(U_e(2,p,elm)+Ye0(elm,1))+Nx{elm}(5,:)*(U_e(5,p,elm)+Ye0(elm,2));
Ze{elm}=Nx{elm}(3,:)*(U_e(3,p,elm)+Ze0(elm,1))+Nx{elm}(6,:)*(U_e(6,p,elm)+Ze0(elm,2));
X_Y_Z_g{elm}=TT1_inv(:,:,elm)*[Xe{elm};Ye{elm};Ze{elm};ums{elm}];
elseif Estrutura=="reticulada"
Xe{elm}=Nx{elm}(1,:)*(U_e(1,p,elm)+Xe0(elm,1))+Nx{elm}(7,:)*((U_e(7,p,elm))+Xe0(elm,2));
Ye{elm}=Nx{elm}(2,:)*(U_e(2,p,elm)+Ye0(elm,1))+Nx{elm}(8,:)*(U_e(8,p,elm)+Ye0(elm,2))+...
Nx{elm}(6,:)*(U_e(6,p,elm)+Rot_ze0(elm,1))+Nx{elm}(12,:)*(U_e(12,p,elm)+Rot_ze0(elm,2));
Ze{elm}=Nx{elm}(3,:)*(U_e(3,p,elm)+Ze0(elm,1))+Nx{elm}(9,:)*(U_e(9,p,elm)+Ze0(elm,2))+...
Nx{elm}(5,:)*(U_e(5,p,elm)+Rot_ye0(elm,1))+Nx{elm}(11,:)*(U_e(11,p,elm)+Rot_ye0(elm,2));
X_Y_Z_g{elm}=TT1_inv(:,:,elm)*[Xe{elm};Ye{elm};Ze{elm};ums{elm}];
end
sum1=Np(elm)+sum1; nf=sum1;
n0=nf-Np(elm)+1;
P_abs(:,n0:nf,p)=[X_Y_Z_g{elm}(1,:); X_Y_Z_g{elm}(2,:); X_Y_Z_g{elm}(3,:)];
if n==1
sigm_x_1(elm,p)=abs(-P1(1,p,elm)/AA(elm));
if Estrutura=="reticulada"
for nn=2:(Np(elm)-1)
my{elm}(nn)=-EE(elm)*(Ze{elm}(nn-1)+Ze{elm}(nn+1)-2*Ze{elm}(nn))/dx^2;
mz{elm}(nn)=EE(elm)*(Ye{elm}(nn-1)+Ye{elm}(nn+1)-2*Ye{elm}(nn))/dx^2;
sigm_x_max{elm}(nn,p)=+abs(my{elm}(nn)*max_z(elm))+abs(mz{elm}(nn)*max_y(elm))...
+sigm_x_1(elm,p);
end
sigm_x_max{elm}(1,p)=sigm_x_max{elm}(2,p);
elseif Estrutura=="articulada"
sigm_x_max{elm}(:,p)=sigm_x_1(elm,p);
end
elseif n==2
sum2=Np(elm)+sum2;
nf=sum2+elm*2; n0=nf-Np(elm)-1;
sigm_x_max_T(n0:nf,p)=[NaN sigm_x_max{elm}(:,p)' NaN NaN];
end
end
end
Vert=permute(P_abs,[2 1 3]);
%% Atribuição de cor a deformada do sistema
if n==2
for p=1:lgt_t
min_sigm(p)=min(sigm_x_max_T(:,p)); max_sigm(p)=max(sigm_x_max_T(:,p));
v_sigm=max_sigm(p)-min_sigm(p);
if v_sigm~=0
cor(:,p)=round((N_cor-1)*((sigm_x_max_T(:,p)-min_sigm(p))/v_sigm))+1;
ccc=cor(:,p);
ind=ccc(~isnan(ccc)); CC=C(ind,:);
n0=1; nf=0; n00=0; nff=0;
for elm=1:N_elem
nf=Np(elm)+n0-2;
n00=1+nff; nff=Np(elm)+n00-2;
Cor(n0:nf,:,p)=CC(n00:nff,:);
Cor(nf+1,:,p)=CC(nff,:);
n0=nf+2;
end
else
for n_p=1:Nt_p
40
Cor(n_p,:,p)=C(1,:);
end
max_sigm(p)=(1.0e-6)*(min_sigm(p)+1);
end
end
end
Se for escolhido ver os resultados no MATLAB, é apresentada uma figura com a
deformada da estrutura a cores ao longo dos instantes de tempo discretizados, tal como se
pode confirmar nas linhas 434-447. No final da sua visualização, na linha 448, é
questionado se o utilizador deseja alterar os dados da análise no geral. Caso seja escolhida
a opção afirmativa, é ainda questionado se na nova análise as condições iniciais serão as
condições finais da análise anterior e o programa volta a executar o algoritmo a partir da
linha 188 e, por isso, são novamente inseridos os dados relativamente à análise e efetuados
os cálculos, até se chegar novamente à questão anterior ou são apenas exportados os
ficheiros com os resultados, se anteriormente for decidido visualizá-los no GiD. Por outro
lado, se a opção escolhida for negativa face a essa questão anterior, é questionado na linha
464 se o utilizador pretende ver a deformada com um novo valor de amplificação. De
forma semelhante, se esta última questão for respondida afirmativamente o programa volta
a executar o algoritmo a partir da linha 297, de forma contrária o programa acaba. Esta
secção corresponde às linhas 432-475:
%% Visualização da deformada e da tensão axial no MATLAB
if n==2
fig_d=figure; p_max=max(PP_nos(:)); p_min=min(PP_nos(:)); v_p=p_max-p_min;
p2_max=(p_max+0.1*v_p); p2_min=(p_min-0.1*v_p);
axis([p2_min p2_max p2_min p2_max p2_min p2_max]); hold on; grid on;
xlabel(['Xaxis ' unid.comp]); ylabel(['Yaxis ' unid.comp]); zlabel(['Zaxis ' unid.comp]);
p1=patch('Faces',Faces,'Vertices',Vert(:,:,1),'FaceVertexCData',Cor(:,:,1),'EdgeColor',...
'interp','FaceColor','none','LineWidth',2);
cbar=colorbar; colormap(C); cbar.Label.String = ['Tensão axial máxima absoluta ' unid.young];
txt=text(10,10,0,['Instante: ' '0' ' s']); cbar.CDataMapping='scaled'; Div=N_cor+1; pause(2);
for p=1:lgt_t
set(p1,'Vertices',Vert(:,:,p),'FaceVertexCData',Cor(:,:,p));
txt.String=([ 'Instante: ' num2str(t(p)) ' s']); caxis([min_sigm(p) max_sigm(p)]);
cbar.Limits=[min_sigm(p) max_sigm(p)]; cbar.Ticks=linspace(min_sigm(p),max_sigm(p),Div);
pause(dt);
end
Alt_result=questdlg('Deseja alterar os dados da análise no geral?', ...
'Dados da Análise','Sim','Não','Não');
switch Alt_result
case 'Sim'
Analise_ok=false; Ampl_ok=true;
Dad_ant=questdlg(['Deseja usar os deslocamentos e as velocidades iniciais nos GL '...
'obtidos do resultado da análise anterior (Opc1) ou considerá-los como nulos (Opc2)?'],...
'Dados da Análise: condições iniciais','Opc1','Opc2','Opc2');
switch Dad_ant
case 'Opc1'
U0=U(:,end); DU0=(U(:,end)-U(:,end-1))/dt;
case 'Opc2'
U0=zeros(Ngl_l,1); DU0=U0;
end
case 'Não'
Analise_ok=true;
41
Aum_ampl=questdlg(['Deseja visualizar novamente a deformada do sistema com uma '...
'amplificação diferente?'],'Dados da Análise: deformada do sistema','Sim','Não','Não');
switch Aum_ampl
case 'Sim'
dad={'Usar uma razão de amplificação da deformada do sistema de:'};
str=inputdlg(dad,'Dados da Análise: deformada do sistema', [1 80] );
Ampl_n=[1 str2double(str(1))]; close(fig_d);
case 'Não'
Ampl_ok=true;
end
end
end
Na última secção do código, as tensões e os deslocamentos em cada instante da análise são
exportados pela função “ToGiD_Trusses3D_resposta.m" para as estruturas articuladas, na
linha 482, ou pela função “ToGiD_BBeams3D_resposta.m” para as estruturas reticuladas,
na linha 488. Estas duas funções foram criadas de forma a gerarem um ficheiro com esses
dados e um ficheiro com as coordenadas dos nós e a conetividade, ambos compatíveis com
o pós-processador do GiD. Esta última parte do código encontra-se nas linhas 476-498:
%% Exportação dos resultados compatíveis com o GiD
if n==1
if Estrutura=="articulada"
if n_c==1
Ut=zeros(Ngl_Stm,lgt_t);
Ut(gl_l,:)=U_ampl;
ToGiD_Trusses3D_resposta(nome_ficheiro,PP_nos,nos_Stm,Ut,sigm_x_1,N_elem,Nnos_el,Nnos,...
lgt_t);
end
break;
elseif Estrutura=="reticulada"
if n_c==1
ToGiD_BBeams3D_resposta(nome_ficheiro,P_abs,Np,sigm_x_max,N_elem,lgt_t);
elseif n_c==2
cic1=false;
break;
end
end
end
cic1=false;
end
end
end
Após vários testes deste programa, foi observado que há uma evidente vantagem do pós-
-processador do GiD em relação ao ambiente do MATLAB para as estruturas articuladas,
pois no GiD é mais intuitiva e cómoda a observação da deformada do sistema, sendo
possível facilmente alterar a velocidade com que esta e as tensões normais são
representadas, como também voltar a vê-las a partir de qualquer instante selecionado. No
entanto, para estruturas reticuladas notou-se que o tempo de geração dos ficheiros de
entrada para o GiD é bastante elevado para um intervalo de tempo (dt) pequeno e que estes
ficheiros ocupam um considerável espaço na memória do disco, já que também existe uma
maior quantidade de nós para ser possível representar a curvatura dos elementos destas
estruturas, que nas estruturas articuladas era naturalmente nula.
42
6 PROGRAMA APLICADO À RESOLUÇÃO DE PROBLEMAS
Para testar e validar o programa “estrutura3D_com_GiD.m”, e outras funções associadas a
este, foram resolvidos vários casos particulares de problemas das diferentes fontes de
referência. Enquanto em alguns destes é pedido para se calcular apenas as
frequências/períodos naturais, noutros são pedidas as respostas do sistema a solicitações
dinâmicas em alguns dos seus GL. Para além do programa desenvolvido neste trabalho
produzir a solução para esses problemas, nos problemas com cargas dinâmicas haverá a
possibilidade de visualizar a deformada e a tensão axial a cores, representativas da sua
magnitude, no MATLAB ou no GiD, conforme a escolha do utilizador. Por fim, estes
resultados serão posteriormente comparados com os resultados da simulação no SAP2000
e/ou com a solução das fontes de referência. É ainda importante referir que quer no
programa quer no SAP2000 as deformações de corte foram desprezadas, já que o modelo
matemático escolhido para o comportamento da viga foi obtido a partir da teoria de Euler-
-Bernoulli, com resultados razoáveis para as estruturas delgadas com um material
isotrópico (5).
6.1 Cálculo das frequências/períodos naturais de sistemas
6.1.1 Estrutura 2D (bidimensional) articulada
Na Figura 7 está representada uma treliça, na qual os seus elementos e nós estão
numerados. Cada quadrado da grelha auxiliar tem 4 pés de largura. As propriedades dos
elementos são: 𝐸 = 29000 ksi, 𝛾(peso específico) = 490 lb/ft3, 𝐴 = 2 in2 (elemento 4 e
5) e 𝐴 = 3 in2 (elemento 6). Existe também uma massa concentrada no nó 6 de
1,0402 [lb. s2/in].
Figura 7: Treliça 2D no SAP2000 (adaptado (6))
43
Para se obter as frequências ou os períodos naturais de um sistema só é preciso executar o
programa até a linha 178. No problema dado é pedido para se utilizar o método em que a
matriz de massa de cada elemento é Consistent, enquanto no SAP2000 esta é sempre
Lumped (7) no programa é possível escolher uma destas. A diferença percentual entre os
resultados foi obtida de seguinte forma:
𝑑. 𝑝. = 100 ∗ (𝑟𝑒𝑠𝑢𝑙𝑡𝑎𝑑𝑜 𝑜𝑏𝑡𝑖𝑑𝑜
𝑟𝑒𝑠𝑢𝑙𝑡𝑎𝑑𝑜 𝑑𝑒 𝑟𝑒𝑓𝑒𝑟ê𝑛𝑐𝑖𝑎− 1)
Tabela 1: Frequências naturais do sistema da Figura 7 e 𝑑. 𝑝. (Referência (6) vs. Programa)
MV Parâmetro SAP2000 Ref. (6) Massa Programa Ref. d.p. (%)
1 Frequência
angular
natural
[rad/s]
487,7613 512,9 Cons. 517,3338 (6) 0,864
Lump. 487,7613 SAP2000 0
2 659,7340 694,6 Cons. 699,7330 (6) 0,739
Lump. 659,7340 SAP2000 0
Através da Tabela 1 é possível verificar que a solução obtida pelo programa está bastante
próxima da solução da referência (6) e é idêntica a solução do SAP2000, e, por isso, existe
grande probabilidade de a formulação do modelo matemático da estrutura ser igual neste
último ao do programa.
6.1.2 Estrutura 1D (unidimensional) reticulada
Na Figura 8 está representada uma viga encastrada-livre com 96 polegadas de
comprimento. Esta viga é discretizada no programa e no SAP2000 com 1, 2, 4, 6, 8, 10 e
96 elementos finitos para averiguar se para uma malha cada vez mais refinada a solução
obtida tende a se aproximar do resultado teórico que também foi usado na referência (8).
Deste modo, pretende-se determinar os primeiros períodos naturais dos MV de flexão desta
apenas, ao se bloquear o 1.º e o 4.º GL de cada nó (e todos os GL do nó fixo), excluindo,
assim, os MV torsores e axiais. As propriedades desta viga são: 𝐴 = 216 in2, 𝐸 = 3600
k/in2, 𝜌 (massa específica) = 2,3 × 10−7 k. s2/in4, 𝐼y = 5832 in4 e 𝐼z = 2592 in4.
Figura 8: Viga encastrada-livre no SAP2000 (adaptado (8))
44
Visto que os resultados do SAP2000 e do programa deste trabalho são iguais, quando a
matriz de massa de cada elemento é Lumped, tal como se pode notar na Tabela 2, a
diferença percentual só é relativa aos resultados do programa e à solução da referência (8),
para cada tipo de matriz de massa considerado e número de elementos finitos.
Tabela 2: Períodos naturais do sistema da Figura 8 e 𝑑. 𝑝. (Referência (8) vs. Programa)
MV Parâmetro N. El. SAP2000 Ref. (8) Massa Programa Ref. d.p. (%)
1
Flexão no
Plano XoY
Período
natural [s]
1 0,054547
0,038005
Cons. 0,0378210
(8)
−0,484
Lump. 0,054547 43,526
2 0,042333 Cons. 0,037983 −0,058
Lump. 0,042333 11,388
4 0,039090 Cons. 0,038000 −0,013
Lump. 0,039090 2,855
6 0,038485 Cons. 0,038001 −0,011
Lump. 0,038485 1,263
8 0,038273 Cons. 0,038001 −0,011
Lump. 0,038273 0,705
10 0,038175 Cons. 0,038001 −0,011
Lump. 0,038175 0,447
96 0,038003 Cons. 0,038001 −0,011
Lump. 0,038003 −0,005
2
Flexão no
Plano XoZ
Período
natural [s]
1 0,036364
0,025337
Cons. 0,025214
(8)
−0,485
Lump. 0,036364 43,521
2 0,028222 Cons. 0,025322 −0,059
Lump. 0,028222 11,387
4 0,026060 Cons. 0,025333 −0,016
Lump. 0,026060 2,854
6 0,025657 Cons. 0,025334 −0,012
Lump. 0,025657 1,263
8 0,025516 Cons. 0,025334 −0,012
Lump. 0,025516 0,706
10 0,025450 Cons. 0,025334 −0,012
Lump. 0,025450 0,446
96 0,025335 Cons. 0,025334 −0,012
Lump. 0,025335 −0,008
3
Flexão no
Plano XoY
Período
natural [s]
1 N/A
0,006064
Cons. 0,003839
(8)
−36,692
Lump. N/A N/A
2 0,008218 Cons. 0,006013 −0,841
Lump. 0,008218 35,521
4 0,006651 Cons. 0,006057 −0,115
Lump. 0,006651 9,680
6 0,006330 Cons. 0,006062 −0,033
Lump. 0,006330 4,387
8 0,006214 Cons. 0,006063 −0,016
Lump. 0,006214 2,474
10 0,006160 Cons. 0,006064 0
Lump. 0,006160 1,583
96 0,006065 Cons. 0,006064 0
Lump. 0,006065 0,016
45
4
Flexão no
Plano XoZ
Período
natural [s]
1 N/A
0,004043
Cons. 0,002559
(8)
−36,705
Lump. N/A N/A
2 0,005479 Cons. 0,004008 −0,866
Lump. 0,005479 35,518
4 0,004434 Cons. 0,004038 −0,124
Lump. 0,004434 9,671
6 0,004220 Cons. 0,004042 −0,025
Lump. 0,004220 4,378
8 0,004143 Cons. 0,004042 −0,025
Lump. 0,004143 2,473
10 0,004107 Cons. 0,004042 −0,025
Lump. 0,004107 1,583
96 0,004043 Cons. 0,004043 0
Lump. 0,004043 0
5
Flexão no
Plano XoY
Período
natural [s]
1 N/A
0,002165
Cons. N/A
(8)
N/A
Lump. N/A N/A
2 N/A Cons. 0,001778 −17,875
Lump. N/A N/A
4 0,002511 Cons. 0,002149 −0,739
Lump. 0,002511 15,982
6 0,002321 Cons. 0,002162 −0,139
Lump. 0,002321 7,206
8 0,002254 Cons. 0,002164 −0,046
Lump. 0,002254 4,111
10 0,002222 Cons. 0,002165 0
Lump. 0,002222 2,633
96 0,002166 Cons. 0,002166 0,046
Lump. 0,002166 0,046
Com o aumento do número de elementos finitos, através da Tabela 2 pode-se verificar uma
convergência significativamente mais rápida para uma massa Consistent do que para uma
massa Lumped, como era de esperar, à custa de um maior esforço computacional.
Enquanto para uma massa Consistent só são precisos (MV – 1) elementos nos MV
superiores a 1 para se ser ter uma d.p. inferior a 1%, quando se considera uma massa
Lumped é preciso um número de elementos igual ao triplo dos MV para se ter uma d.p.
sempre inferior a 5%.
6.1.3 Estrutura 3D reticulada
Na Figura 9 está representada uma estrutura reticulada encastrada nos 4 apoios da base, na
qual se considera que a massa está apenas concentrada em alguns dos nós indicados pelos
números. As propriedades dos elementos desta estrutura são: 𝐴 = 1,0745 in2, 𝐸 = 27900
k/in2, 𝐺 = 10730,769 k/in2, 𝐼y = 0,6657 in4, 𝐼z = 0,6657 in4 e 𝐽 = 1,3315 in4. As
coordenadas e as massas concentradas dos nós estão definidas na Tabela 3 e na Tabela 4,
onde 𝑀1 = 25,35533 × 10−3 lb. s2/in e 𝑀2 = 8,942228 × 10−3 lb. s2/in.
46
Tabela 3: Coordenadas e massas concentradas dos nós da estrutura da Figura 9
Nó 6 7 8 9 10 11 12 13 14
X [in] 0 0 0 0 0 8,625 18,625 27,25 27,25
Y [in] 0 0 0 8,625 17,25 17,25 17,25 17,25 8,625
Z [in] 0 10 18,625 18,625 18,625 18,625 18,625 18,625 18,625
Massa 0 𝑀2 𝑀1 𝑀2 𝑀1 𝑀2 𝑀2 𝑀1 𝑀2
Tabela 4: Coordenadas e massas concentradas dos nós da estrutura da Figura 9 (cont.)
Nó 15 16 17 18 19 20 21 22 23
X [in] 27,25 18,625 8,625 27,25 27,25 27,25 27,25 0 0
Y [in] 0 0 0 0 0 17,25 17,25 17,25 17,25
Z [in] 18,625 18,625 18,625 0 10 0 10 0 10
Massa 𝑀1 𝑀2 𝑀2 0 𝑀2 0 𝑀2 0 𝑀2
Figura 9: Estrutura 3D reticulada no SAP2000 (adaptado (8))
Após o cálculo das frequências naturais deste sistema no SAP2000 e no programa, definiu-
-se a Tabela 5 com esses valores, em que se pode constatar uma d.p. abaixo de 0,004%, o
que indica uma precisão elevada do programa relativamente ao Software comercial. A
razão mais provável de haver divergências dos resultados, embora pouco significativas,
deve-se a diferenças de arredondamento e de truncatura na definição do problema e nos
métodos numéricos adoptados pelo MATLAB ou pelo SAP2000.
47
Tabela 5: Frequências naturais do sistema da Figura 9 e 𝑑. 𝑝. (SAP2000 vs. Programa)
MV Parâmetro SAP2000 Programa d.p. (%)
1
Frequência
natural
[Hz]
115,9221 115,9180 −0,0035 2 122,0175 122,0132 −0,0035 3 143,3283 143,3239 −0,0031 4 224,7677 224,7602 −0,0033 5 421,7819 421,7690 −0,0031 6 449,8606 449,8453 −0,0034 7 481,7900 481,7756 −0,0030 8 596,9226 596,9025 −0,0034 9 801,1479 801,1215 −0,0033 10 827,6355 827,6116 −0,0029 11 939,3811 939,3498 −0,0033 12 988,4880 988,4549 −0,0033 13 1000,4459 1000,4181 −0,0028 14 1015,13048 1015,1000 −0,0030 15 1060,1745 1060,1438 −0,0029 16 1091,2775 1091,2493 −0,0026 17 1095,4768 1095,4451 −0,0029 18 1103,5177 1103,4855 −0,0029 19 1122,1950 1122,1569 −0,0034 20 1159,9266 1159,8907 −0,0031 21 1184,0426 1184,0099 −0,0028 22 1288,9716 1288,9468 −0,0019 23 1324,5057 1324,4634 −0,0032 24 1335,7947 1335,7484 −0,0035 25 1443,9345 1443,8857 −0,0034 26 1467,6785 1467,6395 −0,0027 27 1521,4047 1521,3580 −0,0031 28 1552,3287 1552,2987 −0,0019 29 1570,5755 1570,5297 −0,0029 30 1635,1709 1635,1360 −0,0021 31 1965,7911 1965,7555 −0,0018 32 2005,8865 2005,8472 −0,0020 33 3656,7285 3656,6746 −0,0015 34 3663,5704 3663,5160 −0,0015 35 4514,3438 4514,2776 −0,0015 36 4515,0634 4514,9971 −0,0015 37 4516,4054 4516,3390 −0,0015 38 4517,1792 4517,11280 −0,0015 39 4816,5024 4816,4316 −0,0015 40 4817,4167 4817,3459 −0,0015 41 5312,4586 5312,3807 −0,0015 42 5313,2956 5313,2176 −0,0015
6.2 Análise transitória de sistemas sujeitos a cargas dinâmicas
6.2.1 Estruturas articuladas
6.2.1.1 Estrutura 2D sem amortecimento
48
Na Figura 10 está representada uma treliça, na qual os seus elementos e nós estão
numerados. Cada quadrado da grelha auxiliar tem 5 pés de largura. Os nós 4 e 2 (1 e 2 no
GiD) estão livres, sendo o principal objetivo o de determinar os deslocamentos horizontal e
vertical que podem ocorrer nesses nós ao longo do tempo, como também a tensão axial nos
elementos. Ao nó 4 (1 no GiD) está acoplado um motor que pesa 2000 libras e que com
um desbalanceamento rotativo origina uma força vertical sinusoidal nesse mesmo nó, com
uma intensidade igual a 100 cos(7𝜋𝑡) libras. Os elementos deste sistema têm as mesmas
propriedades: 𝐴 = 1,25 in2, 𝐸 = 30 × 106 lb/in2 e �� = 0,9173 × 10−3 lb. s2/in2.
Através da Tabela 6 é notório que os valores obtidos para as frequências do sistema são
sensíveis ao tipo de matriz de massa considerado, como também já foi visto anteriormente,
mas como o SAP2000 possibilita dar um maior conjunto de resultados do que as próprias
referências dos problemas com soluções, os resultados obtidos no programa serão apenas
comparados com este Software adiante.
Tabela 6: Frequências naturais do sistema da Figura 10 e 𝑑. 𝑝. (SAP2000 vs. Programa e
Referência (9) vs. Programa)
MV Parâmetro SAP2000 Ref. (9) Massa Programa Ref. d.p. (%)
1
Frequência
angular
natural
[rad/s]
55,5337 55,835 Cons. 55,835 (9) 0
Lump. 55,5337 SAP2000 0
2 234,2341 235,46 Cons. 235,45 (9) −0,004
Lump. 234,2341 SAP2000 0
3 1304,8599 1598,4 Cons. 1598,4 (9) 0
Lump. 1304,8599 SAP2000 0
4 1609,1526 1971,6 Cons. 1971,6 (9) 0
Lump. 1609,1526 SAP2000 0
Figura 10: Estrutura 2D articulada no SAP2000 (adaptado (9))
Dependendo da escolha do utilizador, os resultados aos problemas podem ser tanto vistos
no GiD ou no MATLAB. Neste trabalho, optou-se por obter os gráficos dos deslocamentos
49
nos GL e a tensão axial absoluta, juntamente com a deformada do sistema, do próprio
MATLAB, mas também se irá apresentar estes resultados vistos do GiD, nomeadamente
nos casos mais convenientes, de modo a demonstrar que as duas escolhas são possíveis e
válidas em termos de pós-processamento.
Figura 11: Deslocamento do nó 4 (1 no GiD) em X, SAP2000 vs. MATLAB
Como se pode ver pelas Figura 11, 12, 13, 14 e pela Tabela 7, os deslocamentos nos GL
livres estão de acordo com os resultados do SAP2000, pois os gráficos obtidos são muito
semelhantes em ambos e a d.p. é muito pequena. Visto que não é possível obter os gráficos
das tensões no SAP2000, mas é possível ver os diagramas destas em instantes bem
definidos, estas serão comparadas desta forma com as do programa e para alguns instantes
apenas.
Figura 12: Deslocamento do nó 4 (1 no GiD) em Z, SAP2000 vs. MATLAB
Figura 13: Deslocamento do nó 2 em X, SAP2000 vs. MATLAB
50
Figura 14: Deslocamento do nó 2 em Z, SAP2000 vs. MATLAB
Tabela 7: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 10 e 𝑑. 𝑝. (SAP2000 vs. Programa)
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa d.p. (%)
2 Desl.
[in],
Inst.
[s]
X Min −1,310 × 10−3 4,510 × 10−1 −1,311 × 10−3 4,510 × 10−1 0,076 0
Max 1,408 × 10−3 2,820 × 10−1 1,408 × 10−3 2,820 × 10−1 0 0
Z Min −4,601 × 10−4 4,520 × 10−1 −4,601 × 10−4 4,520 × 10−1 0 0
Max 4,949 × 10−4 2,820 × 10−1 4,947 × 10−4 2,820 × 10−1 -0,040 0
4
X Min −3,102 × 10−3 8,530 × 10−1 −3,102 × 10−3 8,530 × 10−1 0 0
Max 2,981 × 10−3 4,470 × 10−1 2,982 × 10−3 4,470 × 10−1 0,034 0
Z Min −1,285 × 10−2 4,500 × 10−1 −1,285 × 10−2 4,500 × 10−1 0 0
Max 1,363 × 10−2 2,830 × 10−1 1,363 × 10−2 2,830 × 10−1 0 0
Figura 15: Tensões axiais nas barras nos instantes 0,009 s e 0,282 s, SAP2000
Através da Figura 15 e da Figura 16 é possível confirmar que os valores das tensões axiais
calculadas no programa correspondem em valor absoluto a solução do SAP2000. A
comparação da deformada da estrutura, sobreposta ao seu estado sem deformação, é feita
com os resultados importados do MATLAB no GiD, tal como se pode ver na Figura 17,
sendo possível confirmar visualmente a coincidência no instante considerado.
51
Figura 16: Tensões axiais absolutas nas barras nos instantes 0,009 s e 0,282 s, MATLAB
Figura 17: Deformada da estrutura no instante 0,282 s, SAP2000 vs. GiD
52
6.2.1.2 Estrutura 2D com amortecimento
Na Figura 18 está representada uma treliça, onde os seus elementos e nós estão numerados.
Cada quadrado da grelha auxiliar tem 60 polegadas de largura. Os nós 1 e 2 (3 e 4 no GiD)
estão livres. No nó 1 são exercidas forças horizontal e vertical, com intensidades iguais a
10. 𝑡 e 5. 𝑡2 libras, respetivamente, durante 1 segundo.
Figura 18: Estrutura 2D articulada no SAP2000 (adaptado (2))
Os elementos deste sistema têm as mesmas propriedades: 𝐴 = 10 in2, 𝐸 = 1 × 104 ksi e
�� = 0,5 k. s2/in2. O amortecimento é igual a 0,1 para todos os MV. De seguida, mostram-
-se os resultados obtidos, tal como no exemplo anterior.
Figura 19: Deslocamento do nó 1 (3 no GiD) em X, SAP2000 vs. MATLAB
53
Figura 20: Deslocamento do nó 1 (3 no GiD) em Z, SAP2000 vs. MATLAB
Figura 21: Deslocamento do nó 2 (4 no GiD) em Z, SAP2000 vs. MATLAB
O deslocamento do nó 2 em X é muito pequeno em relação aos outros deslocamentos, pelo
que não foi representado. Como os esforços aplicados na estrutura são aplicados de uma
forma gradual, principalmente no eixo Z, o efeito da inércia e do amortecimento no
comportamento dinâmico da estrutura é relativamente pequeno, sendo, por isso, a curva de
resposta nesse eixo semelhante a uma parábola e, portanto, aproximadamente igual a
própria forma da curva do esforço.
Tabela 8: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 18 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Nó Parâmetro Eixo SAP2000 Programa d.p. (%)
1 Deslocamento
máximo [in],
Instante [s]
X 2,104 × 10−6 1 2,105 × 10−6 1 0,048 0
Z 1,802 × 10−6 1 1,804 × 10−6 1 0,111 0
2 Z 7,812 × 10−7 1 7,821 × 10−7 1 0,115 0
Pela Figura 19, 20, 21 e pela Tabela 8, os deslocamentos nos GL livres estão de acordo
com os resultados do SAP2000.
54
Figura 22: Tensões axiais nas barras no instante 1 s, SAP2000
Figura 23: Tensões axiais absolutas nas barras no instante 1 s, MATLAB
É possível verificar que os valores das tensões da Figura 22 calculados no SAP2000 estão
dentro ou muito próximos dos valores da Figura 23 calculados pelo programa deste
trabalho.
Figura 24: Deformada da estrutura no instante 1 s, SAP2000 vs. GiD
55
A comparação da deformada da estrutura, sobreposta ao seu estado sem deformação, é
feita através da Figura 24, sendo possível averiguar visualmente a coincidência dos
resultados no instante considerado.
6.2.1.3 Estrutura 3D sem amortecimento
Na Figura 25 está representada uma treliça, na qual os seus elementos e nós estão
numerados. Cada cubo da grelha auxiliar tem 50 polegadas de largura. O nó 3 (4 no GiD)
está livre e sobre o qual foi exercida repentinamente uma força no eixo Z, com uma
intensidade constante igual a 5000 libras. Os elementos deste sistema têm as mesmas
propriedades: 𝐴 = 10 in2, 𝐸 = 3 × 106 psi e �� = 0,1 lb. s2/in2.
Figura 25: Estrutura 3D articulada no SAP2000 (adaptado (2))
De seguida, mostram-se os resultados obtidos, tal como nos casos anteriores, desde 0 s até
0,2 s com a respetiva carga aplicada no sistema.
56
Figura 26: Deslocamento do nó 3 (4 no GiD) em X, SAP2000 vs. MATLAB
Figura 27: Deslocamento do nó 3 (4 no GiD) em Y, SAP2000 vs. MATLAB
Figura 28: Deslocamento do nó 3 (4 no GiD) em Z, SAP2000 vs. MATLAB
Tabela 9: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 25 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Nó Parâmetro Eixo Desl. SAP2000 Programa d.p. (%)
3
Desl.
[in],
Instante
[s]
X Min −1,545 × 10−3 1,858 × 10−1 −1,545 × 10−3 1,858 × 10−1 0 0
Max 8,389 × 10−4 1,810 × 10−2 8,389 × 10−4 1,810 × 10−2 0 0
Y Min −2,457 × 10−3 1,830 × 10−2 −2,457 × 10−3 1,830 × 10−2 0 0
Max 1,783 × 10−4 1,504 × 10−1 1,783 × 10−4 1,504 × 10−1 0 0
Z Min 0 0 0 0 0 0
Max 2,611 × 10−3 1,675 × 10−1 2,611 × 10−3 1,675 × 10−1 0 0
57
Através dos gráficos da Figura 26, 27, 28 e pela Tabela 9, pode-se concluir que os
deslocamentos nos GL livres estão de acordo com os resultados do SAP2000.
Figura 29: Tensões axiais nas barras no instante 0,1675 s, SAP2000
Figura 30: Tensões axiais absolutas nas barras no instante 0,1675 s, MATLAB
Pelas duas figuras anteriores pode-se verificar que o valor máximo da tensão na estrutura é
igual nestas e as tensões noutras barras da Figura 29 estão dentro dos valores da Figura 30.
A examinação da deformada da estrutura, sobreposta ao seu estado sem deformação e vista
de uma perspetiva isométrica, é feita através da Figura 31, sendo possível confirmar a
coincidência desta no instante considerado.
58
Figura 31: Deformada da estrutura no instante 0,1675 s, SAP2000 vs. GiD
6.2.2 Estruturas reticuladas
6.2.2.1 Estrutura 2D sem amortecimento
Na Figura 32 está representada uma estrutura reticulada plana, na qual os seus elementos e
nós estão numerados. Cada quadrado da grelha auxiliar tem 1000 milímetros de largura. Os
nós 2 e 3 (3 e 5 no GiD) estão livres. As propriedades dos elementos deste sistema são:
𝐴 = 240 mm2, 𝐼y = 2000 mm4 e 𝐸 = 2 × 105 N/mm2 para todos os elementos;
�� = 2,2269 × 10−5 N. s2/mm2 para os elementos 1 e 3; �� = 1,8816 × 10−6 N. s2/
/mm2 para o elemento 2.
Figura 32: Estrutura 2D reticulada no SAP2000 (adaptado (9))
59
Sobre o nó 2 é exercida uma força vertical, cuja magnitude é definida ao longo do tempo
pela Tabela 10, sendo depois interpolada linearmente para os restantes instantes da análise.
Tabela 10: Magnitude da força em vários instantes (histórico de carregamento)
Instante da análise [s] 0 0,6 0,601 0,85 1,2 1,4 10
Magnitude da Força [N] 20 50 20 20 −50 0 0
Como se trata de uma estrutura plana, existem 3 GL por cada nó livre, no entanto, o
deslocamento no eixo X é relativamente muito pequeno e, por isso, não será analisado.
Tal como foi feito para as estruturas articuladas, mostram-se agora os resultados obtidos,
neste caso desde 0 s até 2 s da análise, pois, a partir deste instante final, o sistema já se
encontra numa vibração livre não-amortecida e a resposta é praticamente uniforme.
Figura 33: Deslocamento do nó 2 (3 no GiD) em Z, SAP2000 vs. MATLAB
Figura 34: Rotação do nó 2 (3 no GiD) em Y, SAP2000 vs. MATLAB
Tabela 11: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 32 e 𝑑. 𝑝. (SAP2000 vs. Programa)
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa d.p. (%)
2
Desl. Z
[mm]/
Rot. Y
[rad],
Z Min −2,154 × 101 1,277 × 100 −2,153 × 101 1,277 × 100 -0,046 0
Max 1,934 × 101 1,961 × 100 1, 934 × 101 1,961 × 100 0 0
Y Min −1,934 × 10−2 1,961 × 100 −1,934 × 10−2 1,961 × 100 0 0
Max 2,154 × 10−2 1,277 × 100 2,153 × 10−2 1,277 × 100 -0,046 0
60
3
Instante
[s]
Z Min −2,154 × 101 1,277 × 100 −2,153 × 101 1,277 × 100 -0,046 0
Max 1,934 × 101 1,961 × 100 1,934 × 101 1,961 × 100 0 0
Y Min −2,154 × 10−2 1,277 × 100 −2,153 × 10−2 1,277 × 100 -0,046 0
Max 1,934 × 10−2 1,961 × 100 1,934 × 10−2 1,961 × 100 0 0
Figura 35: Deslocamento do nó 3 (5 no GiD) em Z, SAP2000 vs. MATLAB
Figura 36: Rotação do nó 3 (5 no GiD) em Y, SAP2000 vs. MATLAB
Pela tabela anterior e pela Figura 33, 34, 35, 36, pode-se verificar que os deslocamentos
nos GL livres estão de acordo com os resultados do SAP2000. Apesar de na referência
deste problema não se especificar o tipo de secção transversal de cada elemento, para ser
possível calcular a tensão axial resultante da flexão e da tração ou da compressão, assume-
-se que esta está contida num quadrado com 30 milímetros de largura, estando, por isso, a
tensão máxima absoluta situada numa das suas arestas ou num dos seus cantos.
61
Figura 37: Tensões axiais máximas absolutas no instante 1,916 s, MATLAB
Figura 38: Tensões axiais nos pontos críticos, no instante 1,916 s, SAP2000
As tensões representadas sob a forma de diagramas na Figura 38 correspondem a qualquer
ponto da face inferior das vigas horizontais ou da face esquerda da viga vertical. Assim, os
valores absolutos destas estão dentro dos valores da Figura 37 e com um valor máximo nas
zonas onde a estrutura é encastrada.
A comparação da deformada da estrutura, sobreposta ao seu estado sem deformação, é
feita com a Figura 39 e a Figura 40, podendo-se verificar a compatibilidade entre estas
duas em relação ao perfil da deformada.
62
Figura 39: Deformada da estrutura no instante 1,916 s, GiD
Figura 40: Deformada da estrutura no instante 1,916 s, SAP2000
6.2.2.2 Estrutura 2D com amortecimento
Na Figura 41 está representada uma estrutura reticulada plana, com 180 polegadas de
largura e 144 polegadas de altura, na qual os seus elementos e nós estão numerados. Os nós
2 e 3 (4 e 2 no GiD) estão livres, e os nós 1 e 4 (3 e 1 no GiD) têm apenas o GL de rotação
em Y livre.
63
Figura 41: Estrutura 2D reticulada no SAP2000 (adaptado (10))
As propriedades do material que constitui os elementos da estrutura são: 𝐸 = 29 × 106 lb/
/in2 e 𝜈 = 0,3. Os elementos 1 e 3 têm o mesmo perfil de secção transversal, definido pelo
código W8 × 48, enquanto que o do elemento 2 é W12 × 26. Assim, assumindo que o
eixo menor da secção de cada elemento coincide com o eixo Y, de forma a haver uma
maior rigidez no plano da estrutura, as propriedades dos elementos deste sistema são: 𝐼y =
= 182,161 mm4, �� = 1,0360 × 10−2 lb. s2/in2 e 𝐴 = 13,9627 in2 para os elementos 1
e 3; 𝐴 = 7,5682 in2, 𝐼y = 201,7692 mm4 e �� = 5,6118 × 10−3 lb. s2/in2 para o
elemento 2. Nos nós 2 e 3 existe também uma massa agregada de 20 [lb. s2/in]. O
amortecimento é de 0,05 para todos os MV.
Sobre o nó 2 é exercida uma força horizontal periódica não-sinusoidal, cuja magnitude é
definida ao longo do tempo na Tabela 12, durante 1 período, para as 3 frequências
diferentes, sendo depois interpolada linearmente para os restantes instantes da análise,
desde 0 s até 6 s.
Tabela 12: Magnitude da força periódica não-sinusoidal em vários instantes nas 3 análises
diferentes
Instante da Análise 1 [s] 0 0,25 0,75 1
Instante da Análise 2 [s] 0 0,125 0,375 0,5
Instante da Análise 3 [s] 0 0,0625 0,185 0,25
Magnitude da Força [Kip] 0 10 −10 0
Nas seguintes 3 tabelas, apresentam-se os deslocamentos mínimos e máximos, tal como os
instantes em que ocorrem, para cada um dos 8 GL livres, obtidos para os 3 tipos de força.
64
Tabela 13: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
1, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 1 – Força com 1 Hz
Nó Parâ-
metro Eixo Desl. SAP2000 Programa
d.p.
(%)
1
Desl. X
[in]/
Desl. Z
[in]/
Rot. Y
[rad],
Instante
[s]
Y Min −1,717 × 10−2 6,650 × 10−1 −1,717 × 10−2 6,650 × 10−1 0 0
Max 1,917 × 10−2 3,480 × 10−1 1,917 × 10−2 3,480 × 10−1 0 0
2
X Min −1,875 × 100 6,650 × 10−1 −1,875 × 100 6,650 × 10−1 0 0
Max 2,094 × 100 3,480 × 10−1 2,094 × 100 3,480 × 10−1 0 0
Z Min −3,608 × 10−3 6,640 × 10−1 −3,608 × 10−3 6,640 × 10−1 0 0
Max 4,030 × 10−3 3,480 × 10−1 4,030 × 10−3 3,480 × 10−1 0 0
Y Min −4,724 × 10−3 6,650 × 10−1 −4,724 × 10−3 6,650 × 10−1 0 0
Max 5,276 × 10−3 3,480 × 10−1 5,276 × 10−3 3,480 × 10−1 0 0
3
X Min −1,872 × 100 6,640 × 10−1 −1,872 × 100 6,640 × 10−1 0 0
Max 2,0910 × 100 3,480 × 10−1 2,0910 × 100 3,480 × 10−1 0 0
Z Min −4,030 × 10−3 3,480 × 10−1 −4,030 × 10−3 3,480 × 10−1 0 0
Max 3,608 × 10−3 6,640 × 10−1 3,608 × 10−3 6,640 × 10−1 0 0
Y Min −4,712 × 10−3 6,640 × 10−1 −4,712 × 10−3 6,640 × 10−1 0 0
Max 5,265 × 10−3 3,480 × 10−1 5,265 × 10−3 3,480 × 10−1 0 0
4 Y Min −1,714 × 10−2 6,640 × 10−1 −1,714 × 10−2 6,640 × 10−1 0 0
Max 1,915 × 10−2 3,480 × 10−1 1,915 × 10−2 3,480 × 10−1 0 0
Tabela 14: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
2, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 2 – Força com 2 Hz
Nó Parâmetro Eixo Desl. SAP2000 Programa d.p.
(%)
1
Desl. X
[in]/
Desl. Z
[in]/
Rot. Y
[rad],
Instante
[s]
Y Min −1,076 × 10−1 5,993 × 100 −1,076 × 10−1 5,993 × 100 0 0
Max 1,071 × 10−1 5,743 × 100 1,071 × 10−1 5,743 × 100 0 0
2
X Min −1,175 × 101 5,993 × 100 −1,175 × 101 5,993 × 100 0 0
Max 1,170 × 101 5,743 × 100 1,170 × 101 5,743 × 100 0 0
Z Min −2,263 × 10−2 5,993 × 100 −2,263 × 10−2 5,993 × 100 0 0
Max 2,254 × 10−2 5,743 × 100 2,254 × 10−2 5,743 × 100 0 0
Y Min −2,958 × 10−2 5,993 × 100 −2,958 × 10−2 5,993 × 100 0 0
Max 2,947 × 10−2 5,743 × 100 2,947 × 10−2 5,743 × 100 0 0
3
X Min −1,175 × 101 5,993 × 100 −1,175 × 101 5,993 × 100 0 0
Max 1,170 × 101 5,743 × 100 1,170 × 101 5,743 × 100 0 0
Z Min −2,254 × 10−2 5,743 × 100 −2,254 × 10−2 5,743 × 100 0 0
Max 2,263 × 10−2 5,993 × 100 2,263 × 10−2 5,993 × 100 0 0
Y Min −2,958 × 10−2 5,993 × 100 −2,958 × 10−2 5,993 × 100 0 0
Max 2,947 × 10−2 5,743 × 100 2,947 × 10−2 5,743 × 100 0 0
4 Y Min −1,076 × 10−1 5,993 × 100 −1,076 × 10−1 5,993 × 100 0 0
Max 1,071 × 10−1 5,743 × 100 1,071 × 10−1 5,743 × 100 0 0
65
Tabela 15: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
3, nos GL livres do sistema da Figura 41 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 3 – Força com 4 Hz
Nó Parâ-
metro Eixo Desl. SAP2000 Programa
d.p.
(%)
1
Desl. X
[in]/
Desl. Z
[in]/
Rot. Y
[rad],
Instante
[s]
Y Min −8,092 × 10−3 3,230 × 10−1 −8,092 × 10−3 3,230 × 10−1 0 0
Max 8,694 × 10−3 1,640 × 10−1 8,694 × 10−3 1,640 × 10−1 0 0
2
X Min −8,835 × 10−1 3,230 × 10−1 −8,835 × 10−1 3,230 × 10−1 0 0
Max 9,492 × 10−1 1,640 × 10−1 9,492 × 10−1 1,640 × 10−1 0 0
Z Min −1,708 × 10−3 3,230 × 10−1 −1,708 × 10−3 3,230 × 10−1 0 0
Max 1,834 × 10−3 1,640 × 10−1 1,834 × 10−3 1,640 × 10−1 0 0
Y Min −2,223 × 10−3 3,230 × 10−1 −2,223 × 10−3 3,230 × 10−1 0 0
Max 2,388 × 10−3 1,640 × 10−1 2,388 × 10−3 1,640 × 10−1 0 0
3
X Min −8,865 × 10−1 3,220 × 10−1 −8,865 × 10−1 3,220 × 10−1 0 0
Max 9,519 × 10−1 1,640 × 10−1 9,519 × 10−1 1,640 × 10−1 0 0
Z Min −1,834 × 10−3 1,640 × 10−1 −1,834 × 10−3 1,640 × 10−1 0 0
Max 1,708 × 10−3 3,230 × 10−1 1,708 × 10−3 3,230 × 10−1 0 0
Y Min −2,236 × 10−3 3,220 × 10−1 −2,236 × 10−3 3,220 × 10−1 0 0
Max 2,400 × 10−3 1,640 × 10−1 2,400 × 10−3 1,640 × 10−1 0 0
4 Y Min −8,116 × 10−3 3,220 × 10−1 −8,116 × 10−3 3,220 × 10−1 0 0
Max 8,716 × 10−3 1,640 × 10−1 8,716 × 10−3 1,640 × 10−1 0 0
Através da Figura 13, 14 e 15, pode-se verificar novamente que os deslocamentos obtidos
pelo programa estão dentro dos limites do SAP2000, ou seja, entre os seus valores
máximos e mínimos, como também que os instantes em que estes ocorrem coincidem com
os deste Software. Dado que a resposta nos GL do sistema tem a mesma natureza nestes
todos, nas três seguintes figuras será apenas comparado o deslocamento em X do nó 3 para
cada tipo de força.
Figura 42: Deslocamento do nó 3 (2 no GiD) em X, Análise 1, SAP2000 vs. MATLAB
66
Figura 43: Deslocamento do nó 3 (2 no GiD) em X, Análise 2, SAP2000 vs. MATLAB
Figura 44: Deslocamento do nó 3 (2 no GiD) em X, Análise 3, SAP2000 vs. MATLAB
Pela Figura 42, 43, 44 pode-se visualizar que os deslocamentos em X do nó 3 estão de
acordo com os resultados do SAP2000 para as 3 frequências diferentes de força. É possível
também observar o efeito de ressonância, limitada pelo amortecimento presente, na Figura
43, pois a frequência natural mais baixa da estrutura (2,0079 Hz) é quase igual à da força.
Figura 45: Tensões axiais máximas absolutas no instante 0,348 s, Análise 1, MATLAB
67
Figura 46: Tensões axiais nos pontos críticos no instante 0,348 s, Análise 1, SAP2000
Figura 47: Deformada da estrutura no instante 0,348 s, Análise 1, SAP2000 vs. GiD
Figura 48: Tensões axiais máximas absolutas no instante 5,993 s, Análise 2, MATLAB
68
Figura 49: Tensões axiais nos pontos críticos no instante 5,993 s, Análise 2, SAP2000
Figura 50: Deformada da estrutura no instante 5,993 s, Análise 2, SAP2000 vs. GiD
Figura 51: Tensões axiais máximas absolutas no instante 0,164 s, Análise 3, MATLAB
69
Figura 52: Tensões axiais nos pontos críticos no instante 0,164 s, Análise 3, SAP2000
Figura 53: Deformada da estrutura no instante 0,164 s, Análise 3, SAP2000 vs. GiD
A comparação das Figura 45 e 46, 48 e 49, e 51 e 52 permite confirmar que os valores das
tensões axiais máximas absolutas estão de acordo com a solução do SAP2000, pois os
valores obtidos pelo programa estão próximos desta e distribuídos de forma semelhante ao
longo da estrutura, com a tensão máxima global nas extremidades da viga horizontal. No
caso da Análise 2, os valores máximos das tensões da Figura 48 excedem largamente a
tensão de cedência do material, e, de facto, isto acontece, ainda que em menor escala, logo
desde o primeiro pico de oscilação da resposta, o que revela que o efeito de ressonância
neste caso é de todo prejudicial.
Na Figura 47, 50 e 53 é feita a comparação da deformada da estrutura para cada tipo de
força, sobreposta ao seu estado sem deformação, pelos quais se pode observar a
correspondência.
6.2.2.3 Estruturas 3D sem amortecimento
6.2.2.3.1 Exemplo 1
Na Figura 54 está representada uma estrutura reticulada a 3D, onde os seus elementos e
nós estão numerados. Cada cubo da grelha auxiliar tem 200 polegadas de largura. O nó 1
70
(2 no GiD) está livre e sobre o qual foi exercida repentinamente uma força no eixo Z, com
uma intensidade constante igual a 5000 libras durante 0,1 s. As propriedades do material
que constitui os elementos da estrutura são: 𝐸 = 30 × 106 lb/in2 e 𝐺 = 12 × 106 lb/in2.
Outras propriedades são: 𝐼y = 𝐼z = 200 in4, �� = 0,2 lb. s2/in2, 𝐴 = 50 in2, 𝐼0 = 205
in4 e 𝐽 = 40 in4 para os elementos 1 e 2; 𝐼y = 𝐼z = 64 in4, �� = 0,1 lb. s2/in2, 𝐴 = 28
in2, 𝐼0 = 68 in4 e 𝐽 = 12,8 in4 para os elementos 3 e 4. Como a referência deste problema
não especifica as dimensões da secção transversal de cada elemento, foi admitido que esta
está contida num quadrado com 12 polegadas de largura.
Figura 54: Estrutura 3D reticulada (exemplo 1) no SAP2000 (adaptado (2))
Existem 6 GL livres no total, no entanto, o deslocamento no eixo Y, a rotação no eixo X e
Z do nó livre são muito pequenos, e, desta forma, o SAP2000 não possibilita a
representação dos gráficos destes, pelo que também não serão analisados neste trabalho.
Tal como foi feito anteriormente, a partir daqui são apresentados os resultados obtidos,
neste caso desde 0 s até 0,1 s da análise.
Tabela 16: Deslocamentos mínimos e máximos, e os instantes em que ocorrem, nos GL
livres do sistema da Figura 54 e 𝑑. 𝑝. (SAP2000 vs. Programa)
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa
d.p.
(%)
1
Desl. X
[in]/
Desl. Z
[in]/
Rot. Y
[rad],
Instante
[s]
X Min −4,734 × 10−6 8,340 × 10−2 −4,734 × 10−6 8,340 × 10−2 0 0
Max 4,606 × 10−6 9,310 × 10−2 4,606 × 10−6 9,310 × 10−2 0 0
Z Min 0 0 0 0 0 0
Max 1,331 × 10−3 4,440 × 10−2 1,331 × 10−3 4,440 × 10−2 0 0
Y Min −4,963 × 10−6 4,440 × 10−2 −4,963 × 10−6 4,440 × 10−2 0 0
Max 0 0 0 0 0 0
71
Figura 55: Deslocamento do nó 1 (2 no GiD) em X, SAP2000 vs. MATLAB
Figura 56: Deslocamento do nó 1 (2 no GiD) em Z, SAP2000 vs. MATLAB
Figura 57: Rotação do nó 1 (2 no GiD) em Y, SAP2000 vs. MATLAB
Pela tabela anterior e pela Figura 55, 56 e 57, pode-se verificar que os deslocamentos nos
GL livres coincidem com os do SAP2000.
Como se pode ver através da Figura 58 e da Figura 59, a tensão máxima absoluta ocorre
nos extremos de cada elemento, mas é o elemento 1, coincidente com a direção da força,
que apresenta os valores máximos no geral, embora, ainda assim, estes estejam bastante
abaixo da tensão de cedência do material.
Pela Figura 60 é também possível discernir que a deformada obtida pelo programa está
correta no instante considerado, através da comparação visual com a solução do SAP2000.
72
Figura 58: Tensões axiais máximas absolutas no instante 0,0444 s, MATLAB
Figura 59: Tensões axiais nos pontos críticos no instante 0,0444 s, SAP2000
Figura 60: Deformada da estrutura no instante 0,0444 s, SAP2000 vs. GiD
73
6.2.2.3.2 Exemplo 2
Na Figura 61 está representada uma estrutura reticulada 3D, na qual os seus elementos e
nós estão numerados. O seu modelo de análise será igual ao de esta ou terá um número
maior de elementos finitos por cada elemento da estrutura. O cubo da grelha auxiliar que a
contém tem 120 polegadas de largura. Os nós 2, 3, 4, 5 e 8 estão livres e o nó 7 pode
deslocar-se apenas na direção do eixo Y, sobre o qual foi exercida uma força triangular
nesse mesmo eixo, cuja variação da magnitude está respresentada na Tabela 17.
As propriedades do material que constitui os elementos da estrutura são: 𝐸 = 30 × 106 lb/
/in2 e 𝜈 = 0,25. Outras propriedades são: 𝐼z = 400 mm4, 𝐼y = 300 mm4, �� = 0,2 lb.
. s2/in2, 𝐴 = 20 in2 e 𝐽 = 500 in4. Nos nós 2, 3, 4, 5 e 7 existe também uma massa
agregada de 10 [lb. s2/in]. O amortecimento é considerado nulo.
Como a referência deste problema não especifica as dimensões da secção transversal de
cada elemento nem as respetivas orientações, foi admitido que a anterior está contida num
quadrado com 12 polegadas de largura e as orientações dos referenciais locais de cada
elemento são os mesmos que o SAP2000 assume por defeito (ver a referência (7)).
Figura 61: Estrutura 3D reticulada (exemplo 2) no SAP2000 (adaptado (2))
Tabela 17: Magnitude da força triangular em vários instantes
Instante [s] 0 0,02 0,06 0,08
Magnitude da Força [Kip] 0 1 −1 0
74
Nas seguintes 3 tabelas, apresentam-se os deslocamentos mínimos e máximos, obtidos para
uma malha com o número de elementos finitos variável, tal como os instantes em que
ocorrem, para alguns dos nós da Figura 61, sendo, assim, uma amostra de resultados não
muita extensa, mas que, ainda assim, é considerada suficiente para comprovar precisão do
programa relativamente ao SAP2000.
Tabela 18: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
1, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 1 – 1 elemento finito
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa d.p. (%)
7
Desl. Y
[in],
Instante
[s]
Y
Min 0 0 -1,627×10-8 3,00×10-4 N/A N/A
Max 2,334×10-3 5,51×10-2 2,334×10-3 5,51×10-2 0 0
4 Desl. X
[in]/
Desl. Y
[in]/
Desl. Z
[in]/
Rot. X
[rad]/
Rot. Y
[rad]/
Rot. Z
[rad],
Instante
[s]
X Min -1,502×10-6 2,32×10-2 -1,506×10-6 2,33×10-2 N/A N/A
Max 2,520×10-4 6,98×10-2 2,520×10-4 6,98×10-2 0 0
Y Min 0 0 -1,881×10-10 7,00×10-4 N/A N/A
Max 2,414×10-3 5,54×10-2 2,414×10-3 5,54×10-2 0 0
Z Min 0 0 -1,575×10-12 9,00×10-4 N/A N/A
Max 2,465×10-3 5,90×10-2 2,465×10-3 5,90×10-2 0 0
X Min -5,345×10-8 7,60×10-3 -5,417×10-8 7,60×10-3 N/A N/A
Max 2,222×10-5 5,93×10-2 2,222×10-5 5,93×10-2 0 0
Y Min -1,025×10-5 6,01×10-2 -1,025×10-5 6,01×10-2 0 0
Max 0 0 1,114×10-13 7,00×10-4 N/A N/A
Z Min 0 0 -6,685×10-13 7,00×10-4 N/A N/A
Max 9,218×10-6 5,31×10-2 9,217×10-6 5,31×10-2 -0,011 0
5
X Min -1,034×10-3 8,00×10-2 -1,034×10-3 8,00×10-2 0 0
Max 0 0 3,192×10-15 1,20×10-3 N/A N/A
Y Min -1,215×10-4 6,68×10-2 -1,215×10-4 6,68×10-2 0 0
Max 7,475×10-5 2,40×10-2 7,480×10-5 2,40×10-2 0,067 0
Z Min 0 0 -3,837×10-12 1,00×10-3 N/A N/A
Max 2,453×10-3 5,91×10-2 2,453×10-3 5,91×10-2 0 0
X Min -5,460×10-9 4,40×10-3 -5,637×10-9 4,40×10-3 N/A N/A
Max 2,199×10-5 5,88×10-2 2,199×10-5 5,88×10-2 0 0
Y Min -9,898×10-6 6,18×10-2 -9,897×10-6 6,18×10-2 -0,010 0
Max 0 0 1,202×10-11 3,10×10-3 N/A N/A
Z Min -2,134×10-8 3,46×10-2 -2,176×10-8 3,46×10-2 N/A N/A
Max 7,098×10-6 8,00×10-2 7,098×10-6 8,00×10-2 0 0
Tabela 19: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
2, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 2 – 3 elementos finitos
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa d.p. (%)
7
Desl. Y
[in],
Instante
[s]
Y
Min 0 0 0 0 0 0
Max 2,506×10-3 5,40×10-2 2,506×10-3 5,40×10-2 0 0
75
4 Desl. X
[in]/
Desl. Y
[in]/
Desl. Z
[in]/
Rot. X
[rad]/
Rot. Y
[rad]/
Rot. Z
[rad],
Instante
[s]
X Min -5, 147×10-6 2,27×10-2 -5, 147×10-6 2,27×10-2 0 0
Max 3,499×10-4 6,64×10-2 3,499×10-4 6,64×10-2 0 0
Y Min 0 0 0 0 0 0
Max 2,594×10-3 5,48×10-2 2,594×10-3 5,48×10-2 0 0
Z Min 0 0 0 0 0 0
Max 2,632×10-3 5,65×10-2 2,632×10-3 5,65×10-2 0 0
X Min -1,089×10-12 8,00×10-4 -1,074×10-12 8,00×10-4 N/A N/A
Max 2,129×10-5 5,08×10-2 2,129×10-5 5,08×10-2 0 0
Y Min -1,146×10-5 5,44×10-2 -1,146×10-5 5,44×10-2 0 0
Max 0 0 0 0 0 0
Z Min 0 0 0 0 0 0
Max 8,837×10-6 5,07×10-2 8,837×10-6 5,07×10-2 0 0
5
X Min -1,125×10-3 7,51×10-2 -1,125×10-3 7,51×10-2 0 0
Max 9,377×10-14 1,80×10-3 7,926×10-20 2,00×10-4 N/A N/A
Y Min -7,513×10-5 6,68×10-2 -7,513×10-5 6,68×10-2 0 0
Max 5,871×10-5 2,49×10-2 5,871×10-5 2,49×10-2 0 0
Z Min -1,218×10-11 2,10×10-3 -1,229×10-11 2,10×10-3 N/A N/A
Max 2,634×10-3 5,61×10-2 2,634×10-3 5,61×10-2 0 0
X Min -1,403×10-8 8,60×10-3 -1,403×10-8 8,60×10-3 0 0
Max 2,478×10-5 5,64×10-2 2,478×10-5 5,64×10-2 0 0
Y Min -9,898×10-6 6,18×10-2 -9,897×10-6 6,18×10-2 -0,010 0
Max 9,888×10-12 3,00×10-3 2,911×10-12 2,60×10-3 N/A N/A
Z Min -3,177×10-8 3,08×10-2 -3,177×10-8 3,08×10-2 0 0
Max 7,884×10-6 7,85×10-2 7,884×10-6 7,85×10-2 0 0
Tabela 20: Deslocamentos mínimos e máximos, e os instantes em que ocorrem na Análise
3, em alguns GL livres do sistema da Figura 61 e 𝑑. 𝑝. (SAP2000 vs. Programa)
Análise 3 – 10 elementos finitos
N
ó
Parâ-
metro Eixo Desl. SAP2000 Programa d.p. (%)
7
Desl. Y
[in],
Instante
[s]
Y
Min 0 0 0 0 0 0
Max 2,526×10-3 5,38×10-2 2,526×10-3 5,38×10-2 0 0
4
Desl. X
[in]/
Desl. Y
[in]/
Desl. Z
[in]/
Rot. X
[rad]/
Rot. Y
[rad]/
Rot. Z
[rad],
Instante
[s]
X Min -5, 148×10-6 2,34×10-2 -5, 148×10-6 2,34×10-2 0 0
Max 3,501×10-4 6,58×10-2 3,501×10-4 6,58×10-2 0 0
Y Min -1,143×10-12 1,50×10-3 -8,365×10-13 1,50×10-3 N/A N/A
Max 2,614×10-3 5,47×10-2 2,614×10-3 5,47×10-2 0 0
Z Min -3,908×10-11 2,20×10-3 -3,550×10-11 2,20×10-3 N/A N/A
Max 2,650×10-3 5,64×10-2 2,650×10-3 5,64×10-2 0 0
X Min -1,149×10-11 1,60×10-3 -9,018×10-12 1,60×10-3 N/A N/A
Max 2,107×10-5 4,99×10-2 2,107×10-5 4,99×10-2 0 0
Y Min -1,144×10-5 5,35×10-2 -1,144×10-5 5,35×10-2 0 0
Max 9,339×10-14 1,60×10-3 8,984×10-14 1,50×10-3 N/A N/A
Z Min -8,322×10-14 1,50×10-3 -5,495×10-14 1,50×10-3 N/A N/A
Max 8,741×10-6 5,07×10-2 8,741×10-6 5,07×10-2 0 0
5
X Min -1,121×10-3 7,46×10-2 -1,121×10-3 7,46×10-2 0 0
Max 3,314×10-13 2,20×10-3 1,284×10-13 2,20×10-3 N/A N/A
Y Min -7,211×10-5 6,62×10-2 -7,212×10-5 6,62×10-2 0,014 0
Max 5,824×10-5 2,50×10-2 5,824×10-5 2,50×10-2 0 0
76
Z Min -4,227×10-10 3,20×10-3 -4,060×10-10 3,20×10-3 N/A N/A
Max 2,650×10-3 5,60×10-2 2,650×10-3 5,60×10-2 0 0
X Min -1,295×10-9 3,90×10-3 -1,291×10-9 3,90×10-3 N/A N/A
Max 2,484×10-5 5,57×10-2 2,484×10-5 5,57×10-2 0 0
Y Min -1,072×10-5 5,41×10-2 -1,072×10-5 5,41×10-2 0 0
Max 3,544×10-11 3,20×10-3 3,466×10-11 3,20×10-3 N/A N/A
Z Min -2,682×10-8 2,98×10-2 -2,682×10-8 2,98×10-2 0 0
Max 7,915×10-6 7,76×10-2 7,915×10-6 7,76×10-2 0 0
Através da Tabela 18, 19 e 20, pode-se corroborar efetivamente que os deslocamentos
máximos e mínimos obtidos pelo programa estão de acordo com os do SAP2000, tal como
se pode observar que com o aumento do número de elementos finitos os resultados tendem
a convergir. Devido a erros numéricos associados ao cálculo, verifica-se alguma diferença
percentual ou discrepância em alguns resultados que, apesar de serem relativamente
pequenos, se destacam quando os respetivos valores se aproximam muito de zero. Os
instantes em que estes ocorrem também coincidem com os deste Software, exceto,
novamente, para os valores de deslocamento muito próximos de zero. Nas seis seguintes
figuras será apenas comparado o histórico do deslocamento em Y do nó 7 e do nó 5, para
cada uma das análises.
Figura 62: Deslocamento do nó 7 (9 no GiD) em Y, Análise 1, SAP2000 vs. MATLAB
Figura 63: Deslocamento do nó 7 (24 no GiD) em Y, Análise 2, SAP2000 vs. MATLAB
77
Figura 64: Deslocamento do nó 7 (82 no GiD) em Y, Análise 3, SAP2000 vs. MATLAB
Figura 65: Deslocamento do nó 5 (4 no GiD) em Y, Análise 1, SAP2000 vs. MATLAB
Figura 66: Deslocamento do nó 5 (23 no GiD) em Y, Análise 2, SAP2000 vs. MATLAB
Desde a Figura 62 até a Figura 64 e desde a Figura 65 até a Figura 67, com o número de
elementos finitos crescente, verifica-se que os gráficos também têm tendência a estabilizar,
o que acontece de forma mais notória em relação ao deslocamento do nó 7 em Y, como
também existe concordância relativamente aos gráficos presentes em cada uma destas
figuras. É de notar que enquanto o mesmo nó do sistema permanece com a mesma
numeração no SAP2000 após a geração da malha, esta muda no GiD, daí a variação desta
de análise para análise apenas no programa deste trabalho.
78
Figura 67: Deslocamento do nó 5 (83 no GiD) em Y, Análise 3, SAP2000 vs. MATLAB
Figura 68: Deformada da estrutura no instante 0,0538 s, Análise 3, SAP2000 vs. GiD
Na Figura 68 pode-se observar que a deformada obtida pelo programa está certa no
instante considerado.
Como se pode ver através da Figura 69, a tensão máxima absoluta ocorre no elemento 7,
junto do nó 7. Tal como no caso da estrutura anterior, no exemplo 1, esta está bastante
abaixo da tensão de cedência do material. Também se verifica que quando se acha o valor
máximo absoluto de um dos quatro valores da Figura 71 em cada posição axial bem
definida do elemento 7, pode-se notar que este está muito próximo da tensão máxima
absoluta que se regista no diagrama da Figura 70, o que comprova que o cálculo desta pelo
programa foi feito com sucesso.
79
Figura 69: Tensões axiais máximas absolutas no instante 0,0551 s, Análise 1, MATLAB
Figura 70: Diagrama da tensão axial máxima absoluta ao longo do elemento 7, no instante
0,0551 s, Análise 1, GiD
Figura 71: Diagramas das tensões axiais nos 4 pontos (vertíces da secção quadrangular) no
elemento 7, no instante 0,0551s, Análise 1, SAP2000
80
7 CONCLUSÕES
Este trabalho permitiu ganhar alicerces da teoria subjacente ao fenómeno das vibrações de
estruturas, das técnicas numéricas envolvidas com o MEF e as maiores dificuldades que se
encontram na sua análise.
Os códigos desenvolvidos foram testados com sucesso para alguns problemas das fontes
bibliográficas indicadas, com ou sem a presença do amortecimento nos sistemas lineares,
através da comparação com as soluções presentes nestas ou pelos resultados obtidos no
SAP2000. A implementação do programa deste trabalho permitiu também a definição dos
problemas neste Software comercial ser mais compreensível e se ter um melhor
entendimento da solução modal e histórica aos esforços dinâmicos.
O amortecimento nas estruturas, quando presente, foi considerado como viscoso, por ser
vantajoso do ponto de vista matemático, apesar da sua natureza não estar totalmente bem
compreendida atualmente e haver incertezas quanto a outras abordagens. No entanto, o seu
valor pode ser estimado ou correlacionado a partir dos testes laboratoriais feitos em
diferentes materiais ou em estruturas já existentes. Sabe-se que este cresce nos modos mais
altos, depende do tipo de estrutura, dos elementos estruturais, do período e da magnitude
da vibração (2). Porém, na ausência de informações específicas acerca deste, é comum
adaptar um amortecimento entre 0,01 e 0,10 para avaliar a resposta dinâmica dos sistemas
excitados por carregamentos não-uniformes (2). Estes últimos podem ser interpretados
como vários carregamentos harmónicos sobrepostos, podendo algumas destas ter
frequências que coincidam com as frequências naturais do sistema, o que pode causar
elevadas amplitudes de vibração, devido a amplificações dinâmicas que surgem deste
modo. As frequências de excitação mais baixas são as que causam maiores amplitudes de
vibração, sendo, por isso, frequente elevar as frequências naturais mais baixas de estruturas
como medida de prevenção, alterando as características de inércia e de rigidez destas,
ainda na fase de projeto.
A continuidade deste trabalho poderia incluir vários tópicos, nomeadamente, para ser
possível:
• O programa distinguir automaticamente se a análise do problema é 2D ou 3D, pois,
embora o último caso fosse o adotado neste trabalho por ser mais universal, o primeiro
requer uma menor quantidade de parâmetros e um menor volume de cálculos a efetuar;
• Todo o pré-processamento ser feito na interface do GiD;
• Para as estruturas com um número elavado de elementos finitos e com os componentes
harmónicos de baixa frequência nos esforços aplicados, reduzir o número de GL com a
inércia associada, através da condensação estática ou dinâmica;
• Analisar uma estrutura com diferentes tipos de junta, i.e., ser possível definir os GL
rotacionais de cada junta que estejam articulados ou reticulados, e também com diferentes
orientações no espaço, ou seja, que não sejam apenas as coincidentes com o referencial
global;
• Importar ou descrever graficamente os históricos de carregamento distribuídos ao longo
de cada elemento para alguns instantes da análise, que posteriormente seriam interpolados
para os restantes instantes desta;
• O programa calcular as propriedades dinâmicas de elementos não-uniformes a partir da
geometria da secção transversal e de outras propriedades mecânicas que estes possam
apresentar variável ao longo do seu eixo longitudinal. Neste caso, os coeficientes que
definem as matrizes de massa e de rigidez elementares teriam de ser obtidos
individualmente, com a integração numérica no caso mais geral.
81
sim
sim
APÊNDICE A – Fluxograma do programa principal
A.1 – Código “estrutura3D_com_GiD.m” – 1/3
Importação do ficheiro com a mesh do GiD, contendo as coordenadas dos nós da
estrutura (e dos nós auxiliares se esta for reticulada) e da conectividade dos elementos PP_nos, nos_Stm
Escolha do Sistema de Unidades (Imperiais ou Internacionais), definição dos nós com GL livres Unidades, Nno_l
Tipo de estrutura Estrutura
Definição das propriedades mecânicas de cada elemento da estrutura
Definição dos GL livres de cada nó, nós com massa concentrada e o valor dessas massas gl_l_no, nos_Mc, Mc
Articulada? AA, EE, mm
Reticulada? AA, EE, mm, nos_aux, II_y, II_z, JJ, GG, max_y, max_z
max_y,max_y,
Obtenção das matrizes de rotação, de comprimento e de transformação geométrica de cada elemento T1, LL, TT1_inv
Montagem das matrizes globais de rigidez e de massa da estrutura K_gl, M_gl
Cálculo das frequências naturais e dos modos de vibração normalizados do sistema w_n, Phi
1
Definição do amortecimento presente em cada MV, do tempo total da
análise, do intervalo de tempo discretizado e do número de sinusoides tf, dt, Ns
4
82
sim
sim
sim
A.1 – Código “estrutura3D_com_GiD.m” – 2/3
1
Escolha do pós-processador e o valor da amplificação da resposta se for escolhido o MATLAB Pos_proc, Ampl_n
Definição dos nós e dos seus GL livres sob esforços externos no_F, gl_F
Tipo de esforço em cada GL Esforco
Constante?
Sinusoidal?
Não-uniforme?
Magnitude magn, F
Frequência angular w_a, F
Fase fase
Magnitude amplt
Importação do ficheiro com o histórico do esforço t_esf
Interpolação linear do esforço nos instantes do tempo da análise F
Cálculo do fator de participação de cada modo e da resposta do sistema nos seus GL livres yi, U_amp
Pos_proc=MATLAB
Representação dos gráficos da resposta do sistema nos seus GL
Cálculo da deformada do sistema sem
amplificação e da tensão axial sigm_x_max, P_abs
Pos_proc=GiD
2
Cálculo da deformada do
sistema com amplificação P_abs
3
5
83
A.1 – Código “estrutura3D_com_GiD.m” – 3/3
Visualização da figura com a deformada do sistema
2
3
Exportação dos resultados compatíveis com o GiD
Fazer uma nova análise do sistema desde o início? Analise_ok, U0, DU0
Sim Não
4 Fazer uma nova análise do sistema, mas apenas
com uma nova amplificação da deformada? Ampl_n, Ampl_ok
Sim Não
5
Fim do programa
84
APÊNDICE B – Algoritmos
B.1 – Código “estrutura3D_com_GiD.m”
close all; clear; clc; addpath('funcoes_auxiliares/');
%% Dados do Sistema
Nnos_el=2;
dad='Introduza o nome do ficheiro com os dados do GiD : ';
str=inputdlg(dad,['Dados do Sistema: ' 'nome do ficheiro '],[1 80]);
nome_ficheiro=str{1};
[PP_nos, nos_Stm]=read_GiD_file(nome_ficheiro); nos_Stm=sort(nos_Stm,2);
N_elem=size(nos_Stm,1);
nos_u=unique(nos_Stm(:))'; % nós do sistema sem os auxiliares
Nnos_t=size(PP_nos,2); % número total de nós (com os auxiliares)
Unidades=questdlg('Que unidades pretende usar?','Unidades','Imperiais','Internacionais','Imperiais');
dad={'Nós do sistema com pelo menos um GL livre:'};
str=inputdlg(dad,'Dados do Sistema: nós com GL livres',[1 65]); no_l=str2num(str{1});
Nno_l=length(no_l);
switch Unidades
case 'Imperiais'
unid.comp='[in]'; % unidade de comprimento
unid.area='[in^2]'; % unidade de área
unid.forca='[lb]'; % unidade de força
unid.mom='[lb.in]'; % unidade de momento
unid.massa='[lb.s^2/in]'; % unidade de massa
unid.densid_lin='[lb.s^2/in^2]'; % unidade de densidade linear
unid.tempo='[s]'; % unidade de tempo
unid.mom_a='[in^4]'; % unidade de segundo momento de área
unid.young='[lb/in^2]'; % unidade de módulo de Young
case 'Internacionais'
unid.comp='[mm]'; % unidade de comprimento
unid.area='[mm^2]'; % unidade de área
unid.forca='[N]'; % unidade de força
unid.mom='[N.mm]'; % unidade de momento
unid.massa='[N.s^2/mm]'; % unidade de massa
unid.densid_lin='[N.s^2/mm^2]'; % unidade de densidade linear
unid.tempo='[s]'; % unidade de tempo
unid.mom_a='[mm^4]'; % unidade de segundo momento de área
unid.young='[N/mm^2]'; % unidade de módulo de Young
end
Estrutura=questdlg('A estrutura é articulada ou reticulada?','Tipo de estrutura', ...
'articulada','reticulada','articulada');
switch Estrutura
case 'articulada'
Ngb_no=3; AA=zeros(1,N_elem); EE=AA; mm=AA;
ind1=1:3; ind2=4:6;
Nx_f_3D=@(x,Le)[1-(x/Le); 1-(x/Le);1-(x/Le);(x/Le);(x/Le);(x/Le)];
Nnos=Nnos_t; nos_aux=[]; def=repmat({'0'},1,3);
case 'reticulada'
Ngb_no=6; def=repmat({'0'},1,10);
AA=zeros(1,N_elem); EE=AA; mm=AA; II_y=AA; II_z=AA; IIm=AA; JJ=AA; GG=AA;
max_y=AA; max_z=AA;
nos_aux=zeros(N_elem,1);
ind1=1:6; ind2=7:12;
Nx_f_3D=@(x,Le)[1-x/Le; (2*x.^3)/Le^3-(3*x.^2)/Le^2+1; (2*x.^3)/Le^3-(3*x.^2)/Le^2+1;...
1-x/Le; -x.*(x/Le-1).^2; x-(2*x.^2)/Le+x.^3/Le^2; x/Le; (3*x.^2)/Le^2-(2*x.^3)/Le^3;...
(3*x.^2)/Le^2-(2*x.^3)/Le^3; x/Le; x.^2/Le-x.^3/Le^2; x.^3/Le^2-x.^2/Le];
85
Nnos=length(nos_u); % número de nós sem os auxiliares
end
dad={'Quantas secções transversais com propriedades mecânicas diferentes existem?'};
str=inputdlg(dad,'Dados do Sistema: secções dos elementos ', [1 80] ); N_sec=str2double(str(1));
Elem_sec=cell(1,N_sec); dad=Elem_sec;
if N_sec~=1
for p=1:N_sec
dad(p)={['A secção n.º ' num2str(p) ' corresponde aos elementos: ']};
end
str=inputdlg(dad,'Dados do Sistema: secções dos elementos', [1 80] ); Elem_sec=str;
else
Elem_sec{1}=num2str(1:N_elem);
end
for p=1:N_sec
for elem=str2num(Elem_sec{p})
if Estrutura=="articulada"
dad={['Área de secção transversal ' unid.area ':'],...
['Módulo de Young ' unid.young ':'],...
['Densidade linear ' unid.densid_lin ':']};
str=inputdlg(dad,['Dados do Sistema: ' 'elemento ' num2str(elem) ...
' (com a secção n.º ' num2str(p) ')' ],[1 80],def);
AA(elem)=str2double(str(1)); EE(elem)=str2double(str(2)); mm(elem)=str2double(str(3));
def(:)=str(:);
elseif Estrutura=="reticulada"
dad={['Área de secção transversal ' unid.area ':'], ...
['Módulo de Young ' unid.young ':'], ...
['Densidade linear ' unid.densid_lin ':'], ...
['Segundo momento de área em relaçao ao eixo local y ' unid.mom_a ':'],...
['Segundo momento de área em relaçao ao eixo local z ' unid.mom_a ':'],...
['Constante de torção ' unid.mom_a ':'], ...
['Módulo de corte ' unid.young ':'], ...
['Posição máxima absoluta de um ponto da secção no eixo local y ' unid.comp ':'], ...
['Posição máxima absoluta de um ponto da secção no eixo local z ' unid.comp ':'], ...
'Um nó coincidente com o plano local Oxy:'};
str=inputdlg(dad,['Dados do Sistema: ' 'elemento ' num2str(elem) ...
' (com a secção n.º ' num2str(p) ')' ],[1 80],def);
AA(elem)=str2double(str(1)); EE(elem)=str2double(str(2)); mm(elem)=str2double(str(3));
II_y(elem)=str2double(str(4)); II_z(elem)=str2double(str(5)); JJ(elem)=str2double(str(6));
GG(elem)=str2double(str(7)); max_y(elem)=str2double(str(8)); max_z(elem)=str2double(str(9));
nos_aux(elem)=str2double(str(10)); def{10}=num2str(nos_aux(elem));
IIm(elem)=mm(elem)*(II_y(elem)+II_z(elem))/AA(elem); def(:)=str(:);
end
end
end
nn=1; gl_l_no=cell(1,Nno_l); gl_l_no_m=gl_l_no; gl_l_no_sm=gl_l_no; def={''};
Inercia=questdlg('A matriz de massa de cada elemento da estrutura é:','Inércia','Lumped',...
'Consistent','Lumped');
switch Inercia
case 'Consistent'
for no=no_l
dad={['GL livres (de 1 a ' num2str(Ngb_no) ') do nó ' num2str(no) ':']};
str=inputdlg(dad,['Dados do Sistema: GL livres do nó ' num2str(no)],[1 65],def);
gl_l_no{nn}=str2num(str{1}); def=str; bb=gl_l_no{nn}; gl_l_no_m{nn}=gl_l_no{nn}; nn=nn+1;
end
Nmodos=length([gl_l_no_m{:}]);
case 'Lumped'
for no=no_l
dad={['GL livres (de 1 a ' num2str(Ngb_no) ') do nó ' num2str(no) ':']};
86
str=inputdlg(dad,['Dados do Sistema: GL livres do nó ' num2str(no)],[1 65],def);
gl_l_no{nn}=str2num(str{1}); def=str; bb=gl_l_no{nn}; gl_l_no_m{nn}=bb(gl_l_no{nn}<4);
gl_l_no_sm{nn}=bb(~(gl_l_no{nn}<4)); nn=nn+1;
end
Nmodos=length([gl_l_no_m{:}]);
end
dad={'Nós do sistema com massa concentrada:', ['Massa concentrada nesses nós ' unid.massa ':']};
str=inputdlg(dad,'Dados do Sistema: massa concentrada',[1 65]); nos_Mc=str2num(str{1});
Mc=str2num(str{2});
Ngl_l=length([gl_l_no{:}]); Ngl_l_m=Nmodos; Ngl_l_sm=Ngl_l-Ngl_l_m; Ngl_el=Nnos_el*Ngb_no;
Ngl_Stm=Nnos*Ngb_no; gl_no1=zeros(Nnos_el,Ngb_no)'; gl_no1(1:Ngl_el)=1:Ngl_el; gl_no1=gl_no1';
gl_Stm=zeros(N_elem,Ngl_el); gl_no2=zeros(Nnos,Ngb_no)'; gl_no2(1:Ngl_Stm)=1:Ngl_Stm;
gl_no2=gl_no2'; gl_no3=zeros(Nnos_t,Ngb_no); gl_no3(nos_u,:)=gl_no2; gl_l=zeros(1,Ngl_l);
gl_l_m=zeros(1,Ngl_l_m); gl_l_sm=zeros(1,Ngl_l_sm); n0=1; n00=1; n000=1;
for no=1:Nno_l
nf=n0+length(gl_l_no{no})-1; nff=n00+length(gl_l_no_m{no})-1;
nfff=n000+length(gl_l_no_sm{no})-1;
gl_l(1,n0:nf)=gl_no3(no_l(no),gl_l_no{no}); gl_l_m(1,n00:nff)=gl_no3(no_l(no),gl_l_no_m{no});
gl_l_sm(1,n000:nfff)=gl_no3(no_l(no),gl_l_no_sm{no});
n0=nf+1; n00=nff+1; n000=nfff+1;
end
for elem=1:N_elem
for no_e=1:Nnos_el
gl_Stm(elem,gl_no1(no_e,:))=gl_no3(nos_Stm(elem,no_e),:);
end
end
gl_Mc_aux=gl_no3(nos_Mc,:); gl_Mc=gl_Mc_aux(:,1:3);
%% Cálculo das frequências naturais e dos MV normalizados do Sistema
[T1, LL, TT1_inv]=Transf_Cord_3D(PP_nos,nos_Stm,Estrutura,nos_aux);
if Estrutura=="articulada"
[K_gl,K_l]=Rigidez_Glob_truss_3D(LL,EE,T1,AA,gl_Stm,Ngl_Stm);
if Inercia=="Consistent"
[M_gl]=Mass_Consist_Glob_truss_3D(LL,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
elseif Inercia=="Lumped"
[M_gl]=Mass_Lump_Glob_truss_3D(LL,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
end
elseif Estrutura=="reticulada"
[K_gl,K_l]=Rigidez_Glob_bbeam_3D(LL,T1,EE,GG,JJ,II_y,II_z,AA,gl_Stm,Ngl_Stm);
if Inercia=="Consistent"
[M_gl]=Mass_Consist_Glob_bbeam_3D(LL,T1,mm,IIm,AA,gl_Stm,Ngl_Stm,gl_Mc,Mc);
elseif Inercia=="Lumped"
[M_gl]=Mass_Lump_Glob_bbeam_3D(LL,T1,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc);
end
end
K_aa=K_gl(gl_l_m,gl_l_m); Maa=M_gl(gl_l_m,gl_l_m);
if Ngl_l_sm~=0
K_ac=K_gl(gl_l_m,gl_l_sm); K_ca=K_gl(gl_l_sm,gl_l_m); K_cc=K_gl(gl_l_sm,gl_l_sm);
K_A=K_aa-K_ac/K_cc*K_ca;
[phi_a,D]=eig(K_A,Maa); phi_c=-K_cc\K_ca*phi_a; phi=zeros(Ngl_l,Ngl_l_m); aa=1; bb=1; nn=1;
for gl=gl_l
if gl==gl_l_m(aa)
phi(nn,:)=phi_a(aa,:);
if ~(gl==gl_l_m(end))
aa=aa+1;
end
elseif gl==gl_l_sm(bb)
phi(nn,:)=phi_c(bb,:);
if ~(gl==gl_l_sm(end))
87
bb=bb+1;
end
end
nn=nn+1;
end
elseif Ngl_l_sm==0
[phi,D]=eig(K_aa,Maa);
end
[w_n, I]=sort(diag(sqrt(D))); phi=phi(:,I); T=2*pi./w_n; f_n=1./T; Phi=zeros(Ngl_l,Nmodos);
for i=1:Ngl_l
for j=1:Nmodos
Phi(i,j)=phi(i,j)/sqrt(phi(:,j)'*M_gl(gl_l,gl_l)*phi(:,j));
end
end
%% Dados da Análise
U0=zeros(Ngl_l,1); DU0=U0; Analise_ok=false;
while Analise_ok==false
dad={'O amortecimento presente (ou não) é igual para a maior parte dos MV a: ',...
['Os modos de vibração (de 1 a ' num2str(Nmodos) ') com um amortecimento diferente são: ']};
str=inputdlg(dad,'Dados da Análise: Amortecimento', [1 80] );
eps=str2double(str(1))*ones(Nmodos,1); Nmv_am2=str2num(str{2}); n0=1;
for n_m=Nmv_am2
dad={['O amortecimento presente (ou não) no MV ' num2str(Nmv_am2(n0)) ' é igual a:' ]}; n0=n0+1;
str=inputdlg(dad,'Dados da Análise: Amortecimento', [1 80] ); eps(n_m)=str2double(str(1));
end
dad={['Analisar a resposta do sistema desde 0 até ' unid.tempo ':'], ...
['Com um intervalo de tempo discretizado (dt) de ' unid.tempo ':'], ....
'Com um número de sinusoides (Ns) da série de Fourier igual a:'};
str=inputdlg(dad,'Dados da Análise', [1 65] );
t0=0; tf=str2double(str(1));
dt=str2double(str(2));
t=t0:dt:tf; lgt_t=length(t);
Ns=str2double(str(3));
dx=sum(LL)/(1010-1);
pos_proc=questdlg('Onde deseja ver os resultados?','Pós-proces.','No MATLAB','No GiD','No MATLAB');
switch pos_proc
case 'No MATLAB'
dad={'Usar uma razão de amplificação da deformada do sistema de:'};
str=inputdlg(dad,'Dados da Análise: deformada do sistema', [1 70] );
Ampl_n=[1 str2double(str(1))]; n_c=2; Ampl_ok=false;
case 'No GiD'
Ampl_n=1; n_c=1; Ampl_ok=true; Analise_ok=true;
end
dad='Quais são os nós livres do sistema que estão sob algum esforço externo?';
str=inputdlg(dad,'Dados da Análise: nós livres sob esforço',[1 75]); no_F=str2num(str{1});
no_F=sort(no_F); nn=1; nn2=1; gl_F=cell(1,Nno_l); F=zeros(Ngl_l,lgt_t); m00=0;
for no=no_l
if nn2<=length(no_F)
if no_F(nn2)==no
dad=['Entre os GL locais livres (' num2str(gl_l_no{nn}) ') do nó ' num2str(no) ...
', os que estão sob um esforço são:' ];
str=inputdlg(dad,'Dados da Análise: GL livres sob esforço',[1 70]);
gl_F{nn}=str2num(str{1});
for gl=gl_F{nn}
Esforco=questdlg(['O esforço sobre o GL local ' num2str(gl) ' do nó ' num2str(no) ' é:'],...
'Dados da Análise: tipo de esforço','Constante','Sinusoidal',...
'Não-uniforme (interpolado linearmente)','Constante');
dad2=['Dados da Análise: Esforço no GL ' num2str(gl_no3(no,gl)) ' do Sistema'];
88
m0=m00+find(gl==gl_l_no{nn});
switch Esforco
case 'Constante'
dad=['A magnitude do esforço constante ' unid.forca ' ou ' unid.mom ' no GL local '...
num2str(gl) ' do nó ' num2str(no) ' é:'];
str=inputdlg(dad,dad2,[1 80]); magn=str2double(str(1)); F(m0,:)=magn*ones(1,lgt_t);
case 'Sinusoidal'
dad={['A magnitude do esforço sinusoidal ' unid.forca ' ou ' unid.mom ...
' no GL local ' num2str(gl) ' do nó ' num2str(no) ' é:'], ...
'Com uma fase de:', 'E com uma frequência angular de:'};
str=inputdlg(dad,dad2,[1 80]);
amplt=str2double(str(1)); fase=str2double(str(2)); w_a=str2double(str(3));
F(m0,:)=amplt*sin(w_a*t+fase);
case 'Não-uniforme (interpolado linearmente)'
dad=['Introduza o nome do ficheiro com o histórico de carregamento no GL local ' ... num2str(gl) ' do nó ' num2str(no) ' (instante ' unid.tempo ', esforço [ ' ... unid.forca ' ou ' unid.mom ' ] (por linha) ) :' ];
str=inputdlg(dad,dad2,[1 145]); t_esf=load(str{1}); N_esf_l=size(t_esf,1); n0=1;
for n=1:(N_esf_l-1)
tt=t_esf(n,1):dt:(t_esf(n+1,1)-dt); lgt_tt=length(tt);
if lgt_tt>1
nf=n0+length(tt)-1;
m_b=polyfit([tt(1), tt(end)], [t_esf(n,2), t_esf(n+1,2)],1); m=m_b(1); b=m_b(2);
F(m0,n0:nf)=m*tt+b; n0=nf+1;
end
end
end
end
nn2=nn2+1;
end
end
m00=m00+length(gl_l_no{nn}); nn=nn+1;
end
%% Cálculo e visualização da resposta do sistema nos seus GL
[yi] =fator_part_v_amort(M_gl(gl_l,gl_l),w_n,Phi,eps,F,Ns,t,dt,U0,DU0);
U=zeros(Ngl_l,lgt_t);
for k=1:Nmodos
U=yi(k,:).*Phi(:,k)+U; % resposta vista dos GL do sistema
end
if n_c==2
no1=1; n0=0; m0=0;
for no=no_l
f_no=figure; set(f_no,'Name',['Deslocamentos sem amplificação do nó ' num2str(no) ], ...
'NumberTitle','off'); nn=1;
for gl=gl_l_no{no1}
n0=m0+nn;
subplot(2,3,nn);
plot(t,U(n0,:));
grid on; xlabel(['t ' unid.tempo]);
if gl==1
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em X']);
elseif gl==2
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em Y']);
elseif gl==3
ylabel(['deslocamento ' unid.comp]); title(['Deslocamento do nó ' num2str(no) ' em Z']);
elseif gl==4
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo X']);
elseif gl==5
89
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo Y']);
elseif gl==6
ylabel('rotação [rad]'); title(['Rotação do nó ' num2str(no) ' no eixo Z']);
end
nn=nn+1;
end
m0=m0+length(gl_l_no{no1}); no1=no1+1;
end
end
cic1=true;
for n=1:n_c
while Ampl_ok==false || cic1==true
U_ampl=Ampl_n(n)*U;
%% Alguns dados "rescritos" noutra forma
Ustm=zeros(Ngb_no,lgt_t,Nnos_t);
U_e=zeros(Ngl_el,lgt_t,N_elem);
i1=1;
for no=1:length(no_l)
i2=i1+length(gl_l_no{no})-1;
Ustm(gl_l_no{no},:,no_l(no))=U_ampl(i1:i2,:);
i1=i2+1;
end
for elm=1:N_elem
if Estrutura=="articulada"
T1no=T1(:,:,elm);
elseif Estrutura=="reticulada"
T1no=zeros(Ngb_no); T1no(1:3,1:3)=T1(:,:,elm); T1no(4:6,4:6)=T1(:,:,elm);
end
for no=no_l
no_e=1;
for nos_elm=nos_Stm(elm,:)
if no==nos_elm
if no_e==1
U_e(ind1,:,elm)=T1no*Ustm(:,:,nos_elm);
elseif no_e==2
U_e(ind2,:,elm)=T1no*Ustm(:,:,nos_elm);
end
end
no_e=no_e+1;
end
end
end
%% Esforços elementares elásticos nos seus GL
if n==1
P1=U_e;
for elm=1:N_elem
P1(:,:,elm)=K_l(:,:,elm)*U_e(:,:,elm);
end
%% Inicialização de alguns dados e resultados
Np=zeros(N_elem,1);
Nx=repmat({zeros(Ngl_el,1)},N_elem,1);
X_Y_Z_g=cell(N_elem,1); Xe=X_Y_Z_g; Ye=Xe; ums=Xe;
Xe0=zeros(N_elem,Nnos_el); Ye0=Xe0; Ze0=Xe0;
Rot_ye0=Xe0; Rot_ze0=Xe0; Xe0(:,2)=LL;
my=X_Y_Z_g; mz=X_Y_Z_g;
sigm_x_1=zeros(N_elem,lgt_t); sigm_x_max=X_Y_Z_g; sum1=0; nf=0;
for elm=1:N_elem
Le=LL(elm);
90
if Estrutura=="articulada"
n_xx=2;
elseif Estrutura=="reticulada"
n_xx=floor(Le/dx);
end
xx=linspace(0,Le,n_xx);
Np(elm)=length(xx);
ums{elm}=ones(1,Np(elm));
Xe{elm}=xx;
Ye{elm}=zeros(1,Np(elm));
X_Y_Z_g{elm}=zeros(4,Np(elm));
Nx{elm}=Nx_f_3D(xx,Le);
my{elm}=zeros(Np(elm)-1,1); mz{elm}=my{elm};
sigm_x_max{elm}=zeros(Np(elm)-1,lgt_t);
nf=Np(elm)+nf-1; n0=nf-Np(elm)+2;
Faces(n0:nf,1:2)=[(1+sum1):((Np(elm)-1)+sum1); (2+sum1):(Np(elm)+sum1)]';
sum1=Np(elm)+sum1;
end
Nt_p=sum(Np);
sigm_x_max_T=zeros(Nt_p+2*N_elem,lgt_t); cor=sigm_x_max_T;
Ze=Ye; ind=Ye; P_abs=zeros(3,Nt_p,lgt_t); N_cor=9; C=jet(N_cor); Cor=zeros(Nt_p,3,lgt_t);
min_sigm=zeros(1,lgt_t); max_sigm=min_sigm;
end
%% Cálculo da deformada do sistema e da tensão axial
for p=1:lgt_t
sum1=0; sum2=0;
for elm=1:N_elem
if Estrutura=="articulada"
Xe{elm}=Nx{elm}(1,:)*(U_e(1,p,elm)+Xe0(elm,1))+Nx{elm}(4,:)*(U_e(4,p,elm)+Xe0(elm,2));
Ye{elm}=Nx{elm}(2,:)*(U_e(2,p,elm)+Ye0(elm,1))+Nx{elm}(5,:)*(U_e(5,p,elm)+Ye0(elm,2));
Ze{elm}=Nx{elm}(3,:)*(U_e(3,p,elm)+Ze0(elm,1))+Nx{elm}(6,:)*(U_e(6,p,elm)+Ze0(elm,2));
X_Y_Z_g{elm}=TT1_inv(:,:,elm)*[Xe{elm};Ye{elm};Ze{elm};ums{elm}];
elseif Estrutura=="reticulada"
Xe{elm}=Nx{elm}(1,:)*(U_e(1,p,elm)+Xe0(elm,1))+Nx{elm}(7,:)*((U_e(7,p,elm))+Xe0(elm,2));
Ye{elm}=Nx{elm}(2,:)*(U_e(2,p,elm)+Ye0(elm,1))+Nx{elm}(8,:)*(U_e(8,p,elm)+Ye0(elm,2))+...
Nx{elm}(6,:)*(U_e(6,p,elm)+Rot_ze0(elm,1))+Nx{elm}(12,:)*(U_e(12,p,elm)+Rot_ze0(elm,2));
Ze{elm}=Nx{elm}(3,:)*(U_e(3,p,elm)+Ze0(elm,1))+Nx{elm}(9,:)*(U_e(9,p,elm)+Ze0(elm,2))+...
Nx{elm}(5,:)*(U_e(5,p,elm)+Rot_ye0(elm,1))+Nx{elm}(11,:)*(U_e(11,p,elm)+Rot_ye0(elm,2));
X_Y_Z_g{elm}=TT1_inv(:,:,elm)*[Xe{elm};Ye{elm};Ze{elm};ums{elm}];
end
sum1=Np(elm)+sum1; nf=sum1;
n0=nf-Np(elm)+1;
P_abs(:,n0:nf,p)=[X_Y_Z_g{elm}(1,:); X_Y_Z_g{elm}(2,:); X_Y_Z_g{elm}(3,:)];
if n==1
sigm_x_1(elm,p)=abs(-P1(1,p,elm)/AA(elm));
if Estrutura=="reticulada"
for nn=2:(Np(elm)-1)
my{elm}(nn)=-EE(elm)*(Ze{elm}(nn-1)+Ze{elm}(nn+1)-2*Ze{elm}(nn))/dx^2;
mz{elm}(nn)=EE(elm)*(Ye{elm}(nn-1)+Ye{elm}(nn+1)-2*Ye{elm}(nn))/dx^2;
sigm_x_max{elm}(nn,p)=+abs(my{elm}(nn)*max_z(elm))+abs(mz{elm}(nn)*max_y(elm))...
+sigm_x_1(elm,p);
end
sigm_x_max{elm}(1,p)=sigm_x_max{elm}(2,p);
elseif Estrutura=="articulada"
sigm_x_max{elm}(:,p)=sigm_x_1(elm,p);
end
elseif n==2
sum2=Np(elm)+sum2;
91
nf=sum2+elm*2; n0=nf-Np(elm)-1;
sigm_x_max_T(n0:nf,p)=[NaN sigm_x_max{elm}(:,p)' NaN NaN];
end
end
end
Vert=permute(P_abs,[2 1 3]);
%% Atribuição de cor a deformada do sistema
if n==2
for p=1:lgt_t
min_sigm(p)=min(sigm_x_max_T(:,p)); max_sigm(p)=max(sigm_x_max_T(:,p));
v_sigm=max_sigm(p)-min_sigm(p);
if v_sigm~=0
cor(:,p)=round((N_cor-1)*((sigm_x_max_T(:,p)-min_sigm(p))/v_sigm))+1;
ccc=cor(:,p);
ind=ccc(~isnan(ccc)); CC=C(ind,:);
n0=1; nf=0; n00=0; nff=0;
for elm=1:N_elem
nf=Np(elm)+n0-2;
n00=1+nff; nff=Np(elm)+n00-2;
Cor(n0:nf,:,p)=CC(n00:nff,:);
Cor(nf+1,:,p)=CC(nff,:);
n0=nf+2;
end
else
for n_p=1:Nt_p
Cor(n_p,:,p)=C(1,:);
end
max_sigm(p)=(1.0e-6)*(min_sigm(p)+1);
end
end
end
%% Visualização da deformada e da tensão axial no MATLAB
if n==2
fig_d=figure; p_max=max(PP_nos(:)); p_min=min(PP_nos(:)); v_p=p_max-p_min;
p2_max=(p_max+0.1*v_p); p2_min=(p_min-0.1*v_p);
axis([p2_min p2_max p2_min p2_max p2_min p2_max]); hold on; grid on;
xlabel(['Xaxis ' unid.comp]); ylabel(['Yaxis ' unid.comp]); zlabel(['Zaxis ' unid.comp]);
p1=patch('Faces',Faces,'Vertices',Vert(:,:,1),'FaceVertexCData',Cor(:,:,1),'EdgeColor',...
'interp','FaceColor','none','LineWidth',2);
cbar=colorbar; colormap(C); cbar.Label.String = ['Tensão axial máxima absoluta ' unid.young];
txt=text(10,10,0,['Instante: ' '0' ' s']); cbar.CDataMapping='scaled'; Div=N_cor+1; pause(2);
for p=1:lgt_t
set(p1,'Vertices',Vert(:,:,p),'FaceVertexCData',Cor(:,:,p));
txt.String=([ 'Instante: ' num2str(t(p)) ' s']); caxis([min_sigm(p) max_sigm(p)]);
cbar.Limits=[min_sigm(p) max_sigm(p)]; cbar.Ticks=linspace(min_sigm(p),max_sigm(p),Div);
pause(dt);
end
Alt_result=questdlg('Deseja alterar os dados da análise no geral?', ...
'Dados da Análise','Sim','Não','Não');
switch Alt_result
case 'Sim'
Analise_ok=false; Ampl_ok=true;
Dad_ant=questdlg(['Deseja usar os deslocamentos e as velocidades iniciais nos GL '...
'obtidos do resultado da análise anterior (Opc1) ou considerá-los como nulos (Opc2)?'],...
'Dados da Análise: condições iniciais','Opc1','Opc2','Opc2');
switch Dad_ant
case 'Opc1'
U0=U(:,end); DU0=(U(:,end)-U(:,end-1))/dt;
92
case 'Opc2'
U0=zeros(Ngl_l,1); DU0=U0;
end
case 'Não'
Analise_ok=true;
Aum_ampl=questdlg(['Deseja visualizar novamente a deformada do sistema com uma '...
'amplificação diferente?'],'Dados da Análise: deformada do sistema','Sim','Não','Não');
switch Aum_ampl
case 'Sim'
dad={'Usar uma razão de amplificação da deformada do sistema de:'};
str=inputdlg(dad,'Dados da Análise: deformada do sistema', [1 80] );
Ampl_n=[1 str2double(str(1))]; close(fig_d);
case 'Não'
Ampl_ok=true;
end
end
end
%% Exportação dos resultados compatíveis com o GiD
if n==1
if Estrutura=="articulada"
if n_c==1
Ut=zeros(Ngl_Stm,lgt_t);
Ut(gl_l,:)=U_ampl;
ToGiD_Trusses3D_resposta(nome_ficheiro,PP_nos,nos_Stm,Ut,sigm_x_1,N_elem,Nnos_el,Nnos,...
lgt_t);
end
break;
elseif Estrutura=="reticulada"
if n_c==1
ToGiD_BBeams3D_resposta(nome_ficheiro,P_abs,Np,sigm_x_max,N_elem,lgt_t);
elseif n_c==2
cic1=false;
break;
end
end
end
cic1=false;
end
end
end
B.2 – Código “read_GiD_file.m”
function [Coordinates, Elements]=read_GiD_file(nome_ficheiro)
fid = fopen(nome_ficheiro);
str=fgetl(fid);
if str=="MESH dimension 3 ElemType Linear Nnode 2"
str=fgetl(fid);
if str=="Coordinates"
str=fgetl(fid);
while str~="End Coordinates"
num=str2num(str);
Coordinates(:,num(1))=num(2:4);
str=fgetl(fid);
end
93
end
fgetl(fid); str=fgetl(fid);
if str=="Elements"
str=fgetl(fid);
while str~="End Elements"
num=str2num(str);
Elements(num(1),:)=num(2:3);
str=fgetl(fid);
end
end
end
fclose(fid);
B.3 – Código “Rigidez_Glob_truss_3D.m”
function [K_gl,K_l]=Rigidez_Glob_truss_3D(LL,EE,T1,AA,gl_Stm,Ngl_Stm)
n_glb=size(gl_Stm,2);
K_gl=zeros(Ngl_Stm);
K_l=zeros(6,6,length(LL));
for elem=1:length(LL)
A=AA(elem);
E=EE(elem);
Le=LL(elem);
c1=T1(1,1,elem); c2=T1(1,2,elem); c3=T1(1,3,elem);
Ke=(A*E/Le)*[
[ c1^2, c1*c2, c1*c3, -c1^2, -c1*c2, -c1*c3]
[ c1*c2, c2^2, c2*c3, -c1*c2, -c2^2, -c2*c3]
[ c1*c3, c2*c3, c3^2, -c1*c3, -c2*c3, -c3^2]
[ -c1^2, -c1*c2, -c1*c3, c1^2, c1*c2, c1*c3]
[ -c1*c2, -c2^2, -c2*c3, c1*c2, c2^2, c2*c3]
[ -c1*c3, -c2*c3, -c3^2, c1*c3, c2*c3, c3^2]];
Ke_loc=(A*E/Le)*[1 0 0 -1 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0;
-1 0 0 1 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0];
K_l(:,:,elem)=Ke_loc;
for g_i=1:n_glb
i=gl_Stm(elem,g_i);
for g_j=1:n_glb
j=gl_Stm(elem,g_j);
K_gl(i,j)=Ke(g_i,g_j)+K_gl(i,j);
end
end
end
end
B.4 – Código “Rigidez_Glob_bbeam_3D.m”
function [K_gl,K_l]=Rigidez_Glob_bbeam_3D(LL,T1,EE,GG,JJ,II_y,II_z,AA,gl_Stm,Ngl_Stm)
n_glb=size(gl_Stm,2);
K_gl=zeros(Ngl_Stm);
K_l=zeros(12,12,length(LL));
for elem=1:length(LL)
L=LL(elem);
94
E=EE(elem);
Iy=II_y(elem);
Iz=II_z(elem);
A=AA(elem);
G=GG(elem);
J=JJ(elem);
E_Iz=E*Iz; E_Iy=E*Iy; Lp3=L^3; Lp2=L^2;
AA1=(A*E)/L;
BB1=(12*E_Iz)/Lp3;
CC1=(6*E_Iz)/Lp2;
DD1=(12*E_Iy)/Lp3;
EE1= -(6*E_Iy)/Lp2;
FF1=(G*J)/L;
GG1=(4*E_Iy)/L;
HH1=2*E_Iy/L;
II1=(4*E_Iz)/L;
JJ1=2*E_Iz/L;
T=zeros(12);
T11=T1(:,:,elem);
T(1:3,1:3)=T11; T(4:6,4:6)=T11; T(7:9,7:9)=T11; T(10:12,10:12)=T11;
Ke_loc=[[AA1, 0, 0, 0, 0, 0, -AA1, 0, 0, 0, 0, 0]
[0, BB1, 0, 0, 0, CC1, 0, -BB1, 0, 0, 0, CC1]
[0, 0, DD1, 0, EE1, 0, 0, 0, -DD1, 0, EE1, 0]
[0, 0, 0, FF1, 0, 0, 0, 0, 0, -FF1, 0, 0]
[0, 0, EE1, 0, GG1, 0, 0, 0, -EE1, 0, HH1, 0]
[0, CC1, 0, 0, 0, II1, 0, -CC1, 0, 0, 0, JJ1]
[-AA1, 0, 0, 0, 0, 0, AA1, 0, 0, 0, 0, 0]
[0, -BB1, 0, 0, 0, -CC1, 0, BB1, 0, 0, 0, -CC1]
[0, 0, -DD1, 0, -EE1, 0, 0, 0, DD1, 0, -EE1, 0]
[0, 0, 0, -FF1, 0, 0, 0, 0, 0, FF1, 0, 0]
[0, 0, EE1, 0, HH1, 0, 0, 0, -EE1, 0, GG1, 0]
[0, CC1, 0, 0, 0, JJ1, 0, -CC1, 0, 0, 0,II1]];
K_l(:,:,elem)=Ke_loc;
Ke =T'*Ke_loc*T;
for g_i=1:n_glb
i=gl_Stm(elem,g_i);
for g_j=1:n_glb
j=gl_Stm(elem,g_j);
K_gl(i,j)=Ke(g_i,g_j)+K_gl(i,j);
end
end
end
end
B.5 – Código “Mass_Lump_Glob_truss_3D.m”
function [M_gl,M_l]=Mass_Lump_Glob_truss_3D(LL,mm,glb_Stm,Nglb_Stm,glb_Mc,Mc)
n_glb=size(glb_Stm,2);
M_gl=zeros(Nglb_Stm);
I=eye(6);
for elem=1:length(LL)
m=mm(elem);
Le=LL(elem);
Me=0.5*Le*m*I;
for g_i=1:n_glb
95
i=glb_Stm(elem,g_i);
for g_j=1:n_glb
j=glb_Stm(elem,g_j);
M_gl(i,j)=Me(g_i,g_j)+M_gl(i,j);
end
end
end
n0=1;
for nn=Mc
ind=sub2ind([Nglb_Stm Nglb_Stm],glb_Mc(n0,:),glb_Mc(n0,:));
M_gl(ind)=M_gl(ind)+Mc(n0);
n0=n0+1;
end
B.6 – Código “Mass_Consist_Glob_truss_3D.m”
function [M_gl,M_l]=Mass_Consist_Glob_truss_3D(LL,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc)
n_glb=size(gl_Stm,2);
M_gl=zeros(Ngl_Stm);
for elem=1:length(LL)
m=mm(elem);
Le=LL(elem);
Me=(m*Le/6)*[
[ 2, 0, 0, 1, 0, 0]
[ 0, 2, 0, 0, 1, 0]
[ 0, 0, 2, 0, 0, 1]
[ 1, 0, 0, 2, 0, 0]
[ 0, 1, 0, 0, 2, 0]
[ 0, 0, 1, 0, 0, 2]];
for g_i=1:n_glb
i=gl_Stm(elem,g_i);
for g_j=1:n_glb
j=gl_Stm(elem,g_j);
M_gl(i,j)=Me(g_i,g_j)+M_gl(i,j);
end
end
end
n0=1;
for nn=Mc
ind=sub2ind([Ngl_Stm Ngl_Stm],gl_Mc(n0,:),gl_Mc(n0,:));
M_gl(ind)=M_gl(ind)+Mc(n0);
n0=n0+1;
end
end
B.7 – Código “Mass_Lump_Glob_bbeam_3D.m”
function [M_gl]=Mass_Lump_Glob_bbeam_3D(LL,T1,mm,gl_Stm,Ngl_Stm,gl_Mc,Mc)
n_glb=size(gl_Stm,2);
M_gl=zeros(Ngl_Stm);
for elem=1:length(LL)
m=mm(elem);
96
L=LL(elem);
T=zeros(12);
T11=T1(:,:,elem);
T(1:3,1:3)=T11; T(4:6,4:6)=T11; T(7:9,7:9)=T11; T(10:12,10:12)=T11;
Me_loc=0.5*L*m*[1 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0];
Me=T'*Me_loc*T;
for g_i=1:n_glb
i=gl_Stm(elem,g_i);
for g_j=1:n_glb
j=gl_Stm(elem,g_j);
M_gl(i,j)=Me(g_i,g_j)+M_gl(i,j);
end
end
end
n0=1;
for nn=Mc
ind=sub2ind([Ngl_Stm Ngl_Stm],gl_Mc(n0,:),gl_Mc(n0,:));
M_gl(ind)=M_gl(ind)+Mc(n0);
n0=n0+1;
end
end
B.8 – Código “Mass_Consist_Glob_bbeam_3D.m”
function [M_gl]=Mass_Consist_Glob_bbeam_3D(LL,T1,mm,II0,AA,gl_Stm,Ngl_Stm,gl_Mc,Mc)
n_glb=size(gl_Stm,2);
M_gl=zeros(Ngl_Stm);
for elem=1:length(LL) % para cada elemento
m=mm(elem);
L=LL(elem);
I0=II0(elem);
A=AA(elem);
AA1 =22*L;
BB1 =-13*L;
CC1 =(140*I0)/A;
DD1 =(70*I0)/A;
EE1 =4*L^2;
FF1 =-3*L^2;
T=zeros(12);
T11=T1(:,:,elem);
T(1:3,1:3)=T11; T(4:6,4:6)=T11; T(7:9,7:9)=T11; T(10:12,10:12)=T11;
Me_loc=((m*L)/420)*[
[140, 0, 0, 0, 0, 0, 70, 0, 0, 0, 0, 0]
[0, 156, 0, 0, 0, AA1, 0, 54, 0, 0, 0, BB1]
97
[0, 0, 156, 0, -AA1, 0, 0, 0, 54, 0, -BB1, 0]
[0, 0, 0, CC1, 0, 0, 0, 0, 0, DD1, 0, 0]
[0, 0, -AA1, 0, EE1, 0, 0, 0, BB1, 0, FF1, 0]
[0, AA1, 0, 0, 0, EE1, 0, -BB1, 0, 0, 0, FF1]
[70, 0, 0, 0, 0, 0, 140, 0, 0, 0, 0, 0]
[0, 54, 0, 0, 0, -BB1, 0, 156, 0, 0, 0, -AA1]
[0, 0, 54, 0, BB1, 0, 0, 0, 156, 0, AA1, 0]
[0, 0, 0, DD1, 0, 0, 0, 0, 0, CC1, 0, 0]
[0, 0, -BB1, 0, FF1, 0, 0, 0, AA1, 0, EE1, 0]
[0, BB1, 0, 0, 0, FF1, 0, -AA1, 0, 0, 0,EE1]];
Me=T'*Me_loc*T;
for g_i=1:n_glb % para cada glb do elemento
i=gl_Stm(elem,g_i);
for g_j=1:n_glb
j=gl_Stm(elem,g_j);
M_gl(i,j)=Me(g_i,g_j)+M_gl(i,j);
end
end
end
n0=1;
for nn=Mc
ind=sub2ind([Ngl_Stm Ngl_Stm],gl_Mc(n0,:),gl_Mc(n0,:));
M_gl(ind)=M_gl(ind)+Mc(n0);
n0=n0+1;
end
end
B.9 – Código “Transf_Cord_3D.m”
function [T1, LL, TT1_inv]=Transf_Cord_3D(PP_nos,nos_Stm,Estrutura,Nos_aux)
N_elem=size(nos_Stm,1);
T1=zeros(3,3,N_elem);
LL=zeros(1,N_elem);
TT1_inv=zeros(4,4,N_elem);
for elem=1:N_elem
Pi=PP_nos(:,nos_Stm(elem,1)); Pj=PP_nos(:,nos_Stm(elem,2));
L=Pj-Pi;
LL(elem)=sqrt(L'*L);
i_r_x=L(1)/LL(elem);
i_r_y=L(2)/LL(elem);
i_r_z=L(3)/LL(elem);
i_r=[i_r_x,i_r_y,i_r_z];
if Estrutura=="articulada"
P2=[i_r_y,i_r_z,-i_r_x];
elseif Estrutura=="reticulada"
P2=PP_nos(:,Nos_aux(elem))-Pi;
end
Z=cross(i_r,P2);
L_Z_r=sqrt(Z*Z');
k_r_x=Z(1)/L_Z_r;
k_r_y=Z(2)/L_Z_r;
k_r_z=Z(3)/L_Z_r;
k_r=[k_r_x,k_r_y,k_r_z];
j_r=cross(k_r,i_r);
T1(:,:,elem)=[i_r;j_r;k_r];
98
TT1_inv(:,:,elem)=[T1(:,:,elem)', Pi; [0 0 0 1]];
end
end
B.10 – Código “fator_part_v_amort.m”
function [yi] = fator_part_v_amort(M,w_n,phi,eps,F,N,t,h,U0,DU0)
nn=length(t);
t0=t(1); tf=t(nn);
T=tf-t0;
aa=[1 repmat([4 2], 1, floor(0.5*nn - 1)) 4 1];
DT=(2/T)*(h/3);
w_exc=2*pi/T;
n_modos=length(w_n);
yi=zeros(n_modos,length(t));
y0=zeros(n_modos,1);
dy0=zeros(n_modos,1);
for m=1:n_modos
P=phi(:,m)'*F;
y0(m)=phi(:,m)'*M*U0;
dy0(m)=phi(:,m)'*M*DU0;
ki=w_n(m).^2;
r=w_exc/w_n(m);
w_d=w_n(m)*sqrt(1-eps(m)^2);
a0=DT*(aa*P');
sum1=0;
sum2=0;
for n=1:N
P1= P.*cos(n*w_exc*(t-t0));
an=DT*(aa*P1');
P2= P.*sin(n*w_exc*(t-t0));
bn=DT*(aa*P2');
Xampl=1/sqrt((1-(n^2)*r^2 )^2 + ( 2*eps(m)*n*r )^2);
theta=atan2((2*eps(m)*n*r),(1-(n^2)*(r^2)));
sum1=((an*cos(n*w_exc*(t-t0)-theta)+...
+bn*sin(n*w_exc*(t-t0)-theta))*Xampl)+sum1;
sum2=(n*(bn*cos(n*w_exc*(t-t0)-theta)+...
-an*sin(n*w_exc*(t-t0)-theta))*Xampl)+sum2;
end
A=y0(m)-(a0/(2*ki) + (1/ki)*sum1(1));
B=(eps(m)*w_n(m)*A+dy0(m)-(w_exc/ki)*sum2(1))/w_d;
yi(m,:)=exp(-eps(m)*w_n(m)*(t-t0)).*(A*cos(w_d*(t-t0))+...
+B*sin(w_d*(t-t0)))+(a0/(2*ki))+(1/ki)*sum1;
end
end
B.11 – Código “ToGiD_Trusses3D_resposta.m”
function ToGiD_Trusses3D_resposta(nome_ficheiro,PP_nos,nos_Stm,Ut,sigm,N_elem,Nnos_el,Nnos,lgt_t)
msh_file = strcat(nome_ficheiro,'.flavia.msh');
res_file = strcat(nome_ficheiro,'.flavia.res');
eletyp='Linear';
99
fid = fopen(msh_file,'w');
fprintf(fid,'### \n');
fprintf(fid,'# MATfem Trusses3D.m \n');
fprintf(fid,'# \n');
fprintf(fid,'MESH dimension %3.0f Elemtype %s Nnode %2.0f \n \n',3,eletyp,Nnos_el);
fprintf(fid,'coordinates \n');
for i = 1 : Nnos
fprintf(fid,'%6.0f %12.5d %12.5d %12.5d \n',i,PP_nos(1,i),PP_nos(2,i),PP_nos(3,i));
end
fprintf(fid,'end coordinates \n \n');
fprintf(fid,'elements \n');
for i = 1 : N_elem
fprintf(fid,'%6.0f %6.0f %6.0f 1 \n',i,nos_Stm(i,1), nos_Stm(i,2));
end
fprintf(fid,'end elements \n \n');
status = fclose(fid);
fid = fopen(res_file,'w');
fprintf(fid,'Gid Post Results File 1.0 \n');
fprintf(fid,'### \n');
fprintf(fid,'# MAT-fem-Trusses3D \n');
fprintf(fid,'# \n');
fprintf(fid,'GaussPoints "GPTR" Elemtype Linear \n');
fprintf(fid,'Number of Gauss Points: 2 \n');
fprintf(fid,'Natural Coordinates: Internal \n');
fprintf(fid,'end gausspionts \n');
fprintf(fid,'# \n');
for n = 1 : lgt_t
fprintf(fid,['Result "Displacement" "Load Analysis" ' num2str(n) ' Vector OnNodes \n']);
fprintf(fid,'Values \n');
for i = 1 : Nnos
fprintf(fid,'%6.0i %13.5d %13.5d %13.5d \n',i,Ut(i*3-2,n),Ut(i*3-1,n),Ut(i*3,n));
end
fprintf(fid,'End Values \n');
end
fprintf(fid,'GaussPoints "internal" ElemType Linear \n');
fprintf(fid,'Number of Gauss Points: 1 \n');
fprintf(fid,'Natural Coordinates: Internal \n');
fprintf(fid,'end Gausspionts \n');
for n = 1 : lgt_t
fprintf(fid,['Result "Stresses on Elements Linear" "Load Analysis" ' num2str(n) ...
' Scalar OnGaussPoints "internal" \n']);
fprintf(fid,'ComponentNames "sigma_xx" \n');
fprintf(fid,'Values \n');
for i = 1 : N_elem
fprintf(fid,'%6.0i %13.5d \n',i,sigm(i,n));
end
fprintf(fid,'End Values \n');
fprintf(fid,'# \n');
end
status = fclose(fid);
B.12 – Código “ToGiD_BBeams3D_resposta.m”
function ToGiD_BBeams3D_resposta(nome_ficheiro,P_abs,Np,sigm_x_max,N_elem,lgt_t)
100
pnts=size(P_abs,2);
msh_file = strcat(nome_ficheiro,'.flavia.msh');
res_file = strcat(nome_ficheiro,'.flavia.res');
eletyp='Linear';
fid = fopen(msh_file,'w');
fprintf(fid,'### \n');
fprintf(fid,'# MATfem Trusses3D.m \n');
fprintf(fid,'# \n');
fprintf(fid,'MESH dimension %3.0f Elemtype %s Nnode %2.0f \n \n',3,eletyp,2);
fprintf(fid,'coordinates \n');
for i = 1 : pnts
fprintf(fid,'%6.0f %12.5d %12.5d %12.5d \n',i,P_abs(1,i,1),P_abs(2,i,1),P_abs(3,i,1));
end
fprintf(fid,'end coordinates \n \n');
fprintf(fid,'elements \n');
sum1=1; n=1;
for i = 1 : N_elem
for pnt=sum1:(Np(i)-2+sum1)
fprintf(fid,'%6.0f %6.0f %6.0f 1 \n',n,pnt,pnt+1);
n=n+1;
end
sum1=sum1+sum(Np(i));
end
fprintf(fid,'end elements \n \n');
status = fclose(fid);
fid = fopen(res_file,'w');
fprintf(fid,'Gid Post Results File 1.0 \n');
fprintf(fid,'### \n');
fprintf(fid,'# MAT-fem-Trusses3D \n');
fprintf(fid,'# \n');
fprintf(fid,'GaussPoints "GPTR" Elemtype Linear \n');
fprintf(fid,'Number of Gauss Points: 2 \n');
fprintf(fid,'Natural Coordinates: Internal \n');
fprintf(fid,'end gausspionts \n');
fprintf(fid,'# \n');
for n = 1 : lgt_t
fprintf(fid,['Result "Displacement" "Load Analysis" ' num2str(n) ' Vector OnNodes \n']);
fprintf(fid,'Values \n');
for i = 1 : pnts
fprintf(fid,'%6.0i %13.5d %13.5d %13.5d \n',i,P_abs(1,i,n)-P_abs(1,i,1),...
P_abs(2,i,n)-P_abs(2,i,1),P_abs(3,i,n)-P_abs(3,i,1));
end
fprintf(fid,'End Values \n');
end
fprintf(fid,'GaussPoints "internal" ElemType Linear \n');
fprintf(fid,'Number of Gauss Points: 1 \n');
fprintf(fid,'Natural Coordinates: Internal \n');
fprintf(fid,'end Gausspionts \n');
for p = 1 : lgt_t
i=1;
fprintf(fid,['Result "Stresses on Elements Linear" "Load Analysis" ' num2str(p) ...
' Scalar OnGaussPoints "internal" \n']);
fprintf(fid,'ComponentNames "sigma_xx_max" \n');
fprintf(fid,'Values \n');
for elm = 1 : N_elem
for nn=1:(Np(elm)-1)
101
fprintf(fid,'%6.0i %13.5d \n',i,sigm_x_max{elm}(nn,p));
i=i+1;
end
end
fprintf(fid,'End Values \n');
fprintf(fid,'# \n');
end
status = fclose(fid);
102
REFERÊNCIAS
1. F. Teixeira-Dias, J. Pinho-da-Cruz, R.A.F. Valente, R.J.A. Sousa. Método dos
Elementos Finitos: Técnicas de Simulação Numérica em Engenharia. 1.ª ed., Lisboa,
Portugal: ETEP – Edições Técnicas e Profissionais, 2010.
2. Mario Paz, Young Hoon Kim. Structural Dynamics: Theory and Computation. 6.ª ed.,
Cham, Suíça: Springer, 2019.
3. Avelino Alves Filho. Elementos Finitos: A Base da Tecnologia CAE/ Análise Dinâmica.
1.ª ed., São Paulo, Brasil: Érica, 2005.
4. CIMNE. MAT-fem 3D. [Em linha] Versão 1.2, Barcelona, Espanha : CIMNE. [Consult.
Maio 2020]. Disponível em WWW:<URL:http://www.cimne.com/projects/MAT-
fem/descarga.aspx?f=40>.
5. Maria Augusta Neto, Ana Amaro, Luis Roseiro, José Cirne, Rogério Leal.
Engeering Computation of Structures: The Finite Element Method. 1.ª ed., Cham, Suíça:
Springer, 2015.
6. Franklin Y. Cheng. Matrix Analysis of Structural Dynamics: Applications and
Earthquake Engineering. 1.ª ed., Nova Iorque, Estados Unidos da América: Marcel
Dekker, 2001.
7. Computers & Structures, Inc. CSI Analysis Reference Manual. [Livro eletrónico]
Versão 22, Berkeley, Califórnia, Estados Unidos da América : Computers & Structures,
Inc., 2017. [Consult. Julho 2020]. Disponível no Sofware SAP2000 22 Eval:
Help\Documentation\Manuals\Analysis Reference Manual.
8. Computers & Structures, Inc. Analysis Verification Examples. [Recurso eletrónico]
Berkeley, Califórnia, Estados Unidos da América : Computers & Structures, Inc. [Consult.
Agosto 2020]. Disponível no Sofware SAP2000 22 Eval:
Help\Documentation\Verification\Analysis\Frames\(1-014 Eigenvalue Problem | 1-023
ASME Eigenvalue).
9. M. Asghar Bhatti. Fundamental Finite Element Analysis and Applications: with
Mathematica and MATLAB Computations. 1.ª ed., Hoboken, Nova Jérsia, Estados Unidos
da América: Wiley, 2005.
10. Nguyen Lan. SAP2000 Features & A to Z Problems. [Em linha] 1.ª ed., Da Nang,
Vietname : Internal Circulation, 2007. [Consult. Agosto 2020]. Disponível em
WWW:<URL:http://www.erbakan.edu.tr/storage/files/department/insaatmuhendisligi/edito
r/DersSayfalari/Insaat%20Muh%20Bil%20Des%20Coz/sap2000_featuters_and_a_to_z_pr
oblems_book.pdf>.
Recommended