Análise de Estruturas Axissimétricas - Autenticação · de estruturas mencionadas: uma utiliza soluções analíticas e a outra o método dos elementos finitos. A primeira destas

Análise de Estruturas Axissimétricas

Aplicação a Reservatórios Circulares

Cristiano Yudok Chung Rodrigues

Dissertação para obtenção do Grau de Mestre em

Engenharia Civil

Júri

Presidente: Prof. Pedro Guilherme Sampaio Viola Parreira

Orientador: Prof. Manuel da Cunha Ritto Corrêa

Vogais: Prof. António Manuel Figueiredo Pinto da Costa

Outubro 2009

I

Agradecimentos

Em primeiro lugar, agradeço ao meu orientador, o Professor Manuel da Cunha Ritto Corrêa, que

demonstrou grande disponibilidade e interesse em ajudar. Agradeço todos os seus conselhos, ajuda,

explicações e correcções. Agradeço-lhe pelos seus ensinamentos e esclarecimento de dúvidas, pois

permitiram desenvolver os meus conhecimentos na área desta tese.

Agradeço ao Professor Carlos Manuel Tiago Tavares Fernandes, pela sua ajuda, no que diz respeito

à aplicação do método dos elementos finitos a estruturas axissimétricas, que constituiu um trabalho

para cadeira de Análise de Estruturas II. Agradeço todo o seu esclarecimento de dúvidas neste tema

e sugestões fornecidas.

Agradeço ao Professor José Paulo Moitinho de Almeida, por ter disponibilizado aulas de introdução à

programação do método dos elementos finitos e por ter disponibilizado o código fonte de um

programa exemplo. Algumas das ideias que utilizei no meu programa provêm dessas mesmas aulas

e do código fonte disponibilizado, e constituíram algumas das rotinas fundamentais para o

funcionamento do programa, como por exemplo: a criação e utilização do sistema de índices para

espalhamento de matrizes e vectores.

Finalmente, agradeço aos meus pais por todo o seu apoio na minha formação, agradeço à minha

família e amigos que directamente ou indirectamente contribuíram para esta dissertação.

II

III

Resumo

A presente dissertação aborda a análise linear de cascas axissimétricas finas sujeitas a um

carregamento axissimétrico.

O elemento estrutural axissimétrico base analisado é a casca cónica, tendo sido identificadas e

deduzidas as equações que regem o comportamento deste elemento. A escolha do elemento de

casca cónica é justificada pelo facto de ser um elemento axissimétrico geral de curvatura única,

devido à sua simetria de revolução. Definindo-se como incógnita o ângulo com a horizontal é

possível obter alguns casos particulares, como por exemplo: cascas em anel, circulares e cilíndricas.

Desenvolveram-se duas metodologias, baseadas no método dos deslocamentos, para resolver o tipo

de estruturas mencionadas: uma utiliza soluções analíticas e a outra o método dos elementos finitos.

A primeira destas metodologias permite obter resultados exactos (para a teoria estrutural adoptada)

apenas para algumas situações (cascas em anel, circulares e cilíndricas). Em contrapartida, a

segunda metodologia envolve aproximações mas é de aplicação geral a cascas cónicas.

As duas formulações desenvolvidas são implementadas num programa de cálculo com interface

gráfica permitindo uma apresentação e saída de resultados gráfica e numérica. No âmbito do método

dos elementos finitos, o programa também considera cascas esféricas utilizando uma aproximação

«poligonal» por elementos de casca cónica, cuja qualidade aumenta com o número de elementos.

Demonstrou-se as capacidades do programa com alguns exemplos representativos. Validaram-se os

resultados obtidos através da comparação com algumas soluções de referência e, com base numa

análise de erro, verificou-se que a convergência de resultados utilizando o método dos elementos

finitos é relativamente rápida – em muitos casos, bastaram cerca de cinco elementos para obter bons

resultados. Verificou-se que existe compatibilidade entre as duas metodologias e que ambas se

complementam enriquecendo as possibilidades do programa.

Palavras-chave: estruturas axissimétricas, carregamento axissimétrico, cascas finas, soluções

analíticas, método dos elementos finitos, programa de cálculo automático.

IV

V

Abstract

The present paper addresses the linear analysis of axially symmetric structures of thin shell subjected

to axisymmetric loads.

The basic axisymmetric structural element analyzed is the conical shell, and its governing equations

were derived and identified. The choice of the conical shell element is justified by the fact that it

constitutes a general axisymmetric element of single curvature, due to its symmetry of revolution.

Defining the angle with the horizontal as variable it is possible to obtain some particular cases, such

as: circular annular shell, circular shell and cylindrical shell.

Two methodologies were developed based on the displacement method, to solve the type of

structures mentioned: one uses analytical solutions and the other the finite element method. The first

methodology allows obtaining exact solutions (for the adopted structural theory) only for some cases

(circular annular shell, circular shell and cylindrical shell). On the other hand, the second methodology

involves approximations but it is of general application to conical shells.

The two formulations are implemented in a computer program with graphical interface allowing the

presentation and output of graphical and numerical data. Within the finite element method, the

program also considers spherical shells using a “polygonal” approximation by conical shell elements,

whose quality increases with the number of elements.

The capabilities of the program were demonstrated with some representative examples. The obtained

results were validated by comparison with some reference solutions and, by means of an error

analysis, it was verified that the convergence of results using the finite element method is relatively

fast – in many cases, about five elements were enough to get good results. It was possible to verify

that the two methods were compatible and complementary, which enriches the possibilities of the

program.

Keywords: axisymmetric structures, axisymmetric loading, thin shells, analytical solutions, finite

element method, computer program.

VI

VII

Índice

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

2 A condição axissimétrica ................................................................................................................. 2

3 Fundamentos da teoria de cascas cónicas ..................................................................................... 4

3.1 Hipótese Cinemática ............................................................................................................... 5

3.2 Deslocamentos e deformações ............................................................................................... 6

3.3 Equações de equilíbrio .......................................................................................................... 10

3.3.1 Domínio ......................................................................................................................... 10

3.3.2 Fronteira estática........................................................................................................... 13

3.3.3 Fronteira cinemática ...................................................................................................... 13

3.4 Relações constitutivas .......................................................................................................... 14

3.5 Dedução das equações de equilíbrio através do PTV .......................................................... 15

4 Implementação de soluções analíticas ......................................................................................... 18

4.1 Elemento de anel .................................................................................................................. 18

4.1.1 Matriz de rigidez ............................................................................................................ 18

4.1.2 Carregamento ............................................................................................................... 19

4.1.3 Resultados finais ........................................................................................................... 19

4.2 Casca em anel ...................................................................................................................... 19

4.2.1 Deformações ................................................................................................................. 20

4.2.2 Equações de equilíbrio .................................................................................................. 20


4.2.4 Vector das forças de fixação ......................................................................................... 24


4.2.6 Caso particular: elemento de casca circular ................................................................. 26

4.3 Casca cilíndrica ..................................................................................................................... 28

4.3.1 Deformações ................................................................................................................. 28

4.3.2 Equações de equilíbrio .................................................................................................. 28


4.3.4 Vector das forças de fixação ......................................................................................... 35


VIII

5 Implementação do método dos elementos finitos a uma casca cónica ....................................... 39

5.1 Funções de forma ................................................................................................................. 39

5.1.1 Polinómios de Lagrange ............................................................................................... 40

5.1.2 Polinómios de Hermite .................................................................................................. 41

5.2 Elemento finito de casca axissimétrica ................................................................................. 43

5.2.1 Campo de deslocamentos ............................................................................................ 43

5.2.2 Campo de deformações ................................................................................................ 44

5.2.3 Campo de tensões generalizado .................................................................................. 45


5.2.5 Vector de forças nodais equivalentes ........................................................................... 46

5.2.6 Obtenção dos campos de deslocamentos e tensões ................................................... 47

5.3 Regra de quadratura de Gauss............................................................................................. 47

5.3.1 Análise de erro associado ao número de pontos de Gauss ......................................... 48

6 Sistema global ............................................................................................................................... 51

6.1 Transformação Linear de Coordenadas ............................................................................... 51

6.2 Definição do carregamento nodal ......................................................................................... 53

6.3 Reunião de elementos .......................................................................................................... 53

7 Singularidades .............................................................................................................................. 55

7.1 Nós no eixo de simetria ......................................................................................................... 55

7.2 Método dos elementos finitos ............................................................................................... 56

8 Programação ................................................................................................................................. 59

8.1 Objectivo do programa .......................................................................................................... 59

8.2 Programas utilizados ............................................................................................................. 59

8.3 Grau de precisão ................................................................................................................... 60

8.4 Base de dados ...................................................................................................................... 61

8.5 Fluxograma ........................................................................................................................... 64

8.6 Cuidados relacionados com a programação/formulação ..................................................... 65

9 Exemplos de aplicação e erro ....................................................................................................... 67

9.1 Lajes circulares ..................................................................................................................... 69

9.1.1 Exemplo nº 1 ................................................................................................................. 70

9.1.2 Exemplo nº 2 ................................................................................................................. 73

IX

9.1.3 Exemplo nº 3 ................................................................................................................. 76

9.2 Casca cilíndrica – Exemplo nº 4 ........................................................................................... 79

9.3 Casca cónica – Exemplo nº 5 ............................................................................................... 83

9.4 Reservatório – Exemplo nº 6 ................................................................................................ 86

9.4.1 Peso próprio da estrutura .............................................................................................. 87

9.4.2 Acção da água .............................................................................................................. 87

9.4.3 Reacção do solo............................................................................................................ 88

9.4.4 Obtenção dos deslocamentos e esforços ..................................................................... 88

9.5 Análise dos resultados e do erro .......................................................................................... 92

9.5.1 Descontinuidades no campo de esforços ..................................................................... 93

9.5.2 Refinamento adaptativo ................................................................................................ 95

10 Conclusão ................................................................................................................................. 97

11 Bibliografia ............................................................................................................................... 100

A Anexos ...................................................................................................................................... A - 1

A.1 Manual de utilização ......................................................................................................... A - 3

A.1.1 Breve introdução ....................................................................................................... A - 3

A.1.2 Menu ......................................................................................................................... A - 4

A.1.3 Separadores .............................................................................................................. A - 7

A.2 Código fonte .................................................................................................................... A - 17

A.2.1 Unidade “Unit1” (principal) ...................................................................................... A - 17

A.2.2 Unidade “CAxisSim” ................................................................................................ A - 42

A.2.3 Unidade “MEFAxisSim” ........................................................................................... A - 53

A.2.4 Unidade “Equacoes” ............................................................................................... A - 58

A.2.5 Unidade “OpMat” ..................................................................................................... A - 63

A.2.6 Unidade “NPGaussL” .............................................................................................. A - 65

X

XI

Lista de Figuras

Figura 2.1 - Estrutura axissimétrica com a secção transversal .............................................................. 2

Figura 2.2 - Sistema de coordenadas cilíndricas (r,z,θ) ......................................................................... 2

Figura 3.1 - Sistema de coordenadas local ............................................................................................ 4

Figura 3.2 - Sistema de coordenadas global .......................................................................................... 4

Figura 3.3 - Sentido positivo dos esforços com representação no PSA ................................................ 4

Figura 3.4 - Sentido positivo dos esforços representados no plano horizontal ...................................... 4

Figura 3.5 - Deformação generalizada de uma fibra .............................................................................. 5

Figura 3.6 - Campo de deslocamentos em coordenadas globais .......................................................... 6

Figura 3.7 - Campo de deslocamentos em coordenadas locais............................................................. 6

Figura 3.8 - Extensão tangencial representada no PSA ......................................................................... 7

Figura 3.9 - Extensão tangencial projectada no plano perpendicular ao eixo de axissimetria ............... 7

Figura 3.10 - Curvatura tangencial ......................................................................................................... 8

Figura 3.11 - Elemento infinitesimal de uma casca cónica ................................................................... 10

Figura 3.12 - Elemento infinitesimal no plano de simplificação axissimétrica ...................................... 11

Figura 3.13 - Corte do elemento infinitesimal no plano perpendicular ao eixo de revolução ............... 11

Figura 3.14 - Deslocamentos virtuais, forças externas e carregamento .............................................. 15

Figura 4.1 - Elemento geral de casca em anel ..................................................................................... 19

Figura 4.2 - Elemento geral de casca cilíndrica .................................................................................... 28

Figura 4.3 - Configuração de sistemas para a implementação da solução de membrana .................. 31

Figura 4.4 – Plano de corte horizontal de uma casca cilíndrica ........................................................... 36

Figura 5.1 - Elemento e graus de liberdade considerados para deslocamentos longitudinais ............ 40

Figura 5.2 - Função de forma associada a u1 ....................................................................................... 40

Figura 5.3 - Função de forma associada a u2 ....................................................................................... 40

Figura 5.4 - Elemento e graus de liberdade considerados para deslocamentos transversais ............. 41

Figura 5.5 - Função de forma associada a w1 ...................................................................................... 41

Figura 5.6 - Função de forma associada a φ1 ...................................................................................... 41

Figura 5.7 - Função de forma associada a w2 ...................................................................................... 42

Figura 5.8 - Função de forma associada a φ2 ...................................................................................... 42

Figura 5.9 - Elemento de finito de casca cónica representado no PSA ............................................... 43

Figura 5.10 - Erro relativo do nº de pontos de Gauss ........................................................................... 50

Figura 6.1 - Transformação linear entre coordenadas locais e coordenadas globais .......................... 51

Figura 7.1 - Graus de liberdade em sistema de coordenadas locais ................................................... 56

Figura 7.2 - Graus de liberdade em sistema de coordenadas globais ................................................. 56

Figura 9.1 - Configuração geral de uma laje circular ............................................................................ 69

Figura 9.2 - Configuração da estrutura e do carregamento do Exemplo nº 1 ...................................... 70

Figura 9.3 - Deslocamentos e diagramas de esforços do Exemplo nº 1 .............................................. 71

Figura 9.4 - Erro dos deslocamentos e esforços do Exemplo nº 1 ....................................................... 72


XII


Figura 9.7 - Erro dos deslocamentos e esforços do Exemplo nº 2 ....................................................... 75



Figura 9.10 - Erro dos deslocamentos e esforços do Exemplo nº 3 ..................................................... 78

Figura 9.11 - Configuração da estrutura e do carregamento do Exemplo nº 4 .................................... 79

Figura 9.12 - Deslocamentos e diagramas de esforços do Exemplo nº 4 ............................................ 81

Figura 9.13 - Erro dos deslocamentos e esforços do Exemplo nº 4 ..................................................... 82




Figura 9.17 - Modelo de cálculo auxiliar para o elemento de anel ....................................................... 87


Figura 9.19 - Diagramas de esforço na zona de ligação entre o Arco nº1 e o Elemento nº3 .............. 91

Figura 9.20 – Análise de erros provenientes da formulação ................................................................ 94

Figura 9.21 - Aplicação de um refinamento adaptativo ao Exemplo nº 3 ............................................. 96

Lista de Figuras em Anexo

Figura A1 - Programa de cálculo automático “AxisSim” ................................................................... A - 3

Figura A2 - Menu “Ficheiro” .............................................................................................................. A - 4

Figura A3 - Menu “Ver” ..................................................................................................................... A - 5

Figura A4 - Menu "Ferramentas” ...................................................................................................... A - 6

Figura A5 - Menu "Ajuda" .................................................................................................................. A - 6

Figura A6 - Separador "Nós" ............................................................................................................. A - 7

Figura A7 - Eixos do sistema global de coordenadas ...................................................................... A - 7

Figura A8 - Separador "Anéis" .......................................................................................................... A - 8

Figura A9 - Separador "Cascas" ....................................................................................................... A - 9

Figura A10 - Separador "MEF" ....................................................................................................... A - 10

Figura A11 - Separador "Opções"................................................................................................... A - 12

Figura A12 - Separador "Estrutura" ................................................................................................ A - 14

Figura A13 - Caixa de selecção "Exibir" ......................................................................................... A - 14

XIII

Lista de Tabelas

Tabela 5.1 - Erro relativo dos termos da matriz de rigidez em função do nº de pontos de Gauss ...... 49

Tabela 9.1 - Erro dos deslocamentos e esforços do Exemplo nº 1 ...................................................... 72




Tabela 9.5 - Erro do refinamento uniforme e adaptativo aplicado ao Exemplo nº 3 ............................ 96

XIV

Lista de abreviaturas

PSA Plano de simplificação axissimétrica: plano que contém o eixo e a secção transversal de

revolução que geram uma estrutura axissimétrica.

MEF Método dos elementos finitos

PTV Principio dos trabalhos virtuais

Notação

Letras latinas:

Operador diferencial

Área da secção de revolução

Matriz elementar de deformação

Matriz elementar das relações constitutivas

Vector de deslocamentos nodais global;

Vector de deslocamentos nodais elementar

Vector de deslocamentos associado a um nó

Vector de deslocamentos nodais de uma solução particular

Módulo de elasticidade

Vector de forças nodais global;

Vector de forças nodais elementar

Vector elementar de forças de fixação das cargas de vão

Vector elementar de forças aplicadas nos nós

Vector de forças nodais de uma solução particular

Vector de forças associado a um nó

Espessura do elemento

Inércia

XV

Matriz de rigidez global;

Matriz de rigidez elementar;

Momento flector tangencial

Momento flector longitudinal

Esforço de membrana tangencial

Esforço de membrana longitudinal

Vector do carregamento de vão

Carregamento de vão transversal e longitudinal

Matriz elementar de transformação de coordenadas

Matriz de transformação de coordenadas num ponto

Vector elementar de deslocamentos

Deslocamentos de um ponto na fibra de uma secção no referencial global

Deslocamento transversal e longitudinal de um ponto na fibra de uma secção

Deslocamentos de uma secção no referencial global

Deslocamento transversal e longitudinal de uma secção

Esforço transverso

Letras gregas:

Vector elementar do campo de deformações

Extensão segundo a direcção tangencial

Deformação segundo a direcção tangencial

Extensão segundo a direcção longitudinal

Deformação segundo a direcção longitudinal

Coeficiente de Poisson

Vector elementar do campo de tensões generalizado

XVI

Tensão segundo a direcção longitudinal

Tensão segundo a direcção tangencial

Trabalho realizado pelas forças externas

Trabalho realizado pelas forças internas

Curvatura segundo a direcção longitudinal

Curvatura segundo a direcção tangencial

Matriz elementar das funções de forma

1

1 Introdução

A presente dissertação aborda a análise linear de cascas axissimétricas finas sujeitas a um

carregamento axissimétrico. Tirando partido da condição de axissimetria e das hipóteses pertinentes

ao comportamento de cascas finas é possível obter soluções de boa qualidade com base numa

formulação unidimensional.

A estrutura axissimétrica base a ser estudada é a casca cónica com um ângulo de inclinação geral

em relação à horizontal. Pretende-se estudar as equações que regem a casca cónica, determinando

nomeadamente:

Equações de compatibilidade, relacionando deformações com deslocamentos;

Equações de equilíbrio;

Relações constitutivas, relacionando as tensões com as deformações.

Após o estudo da casca cónica, adopta-se dois métodos de resolução de estruturas, ambos

conduzindo a uma formulação em que as variáveis independentes são deslocamentos, sejam eles:

Implementação de soluções analíticas, determinando-se soluções para elementos de anel e

casos particulares das cascas cónicas como: cascas em anel, circulares e cilíndricas;

Implementação do método dos elementos finitos ao elemento de casca tronco-cónica.

Pretende-se definir a formulação completa das duas metodologias para a resolução geral de um

problema axissimétrico compatível com as possibilidades dos métodos.

Pretende-se desenvolver um programa de cálculo automático que utiliza a formulação definida nas

duas metodologias para a resolução de estruturas axissimétricas. O programa terá uma componente

de entrada de dados, manual ou com recurso a ficheiro de dados. Também terá uma componente de

saída de dados, permitindo a obtenção de resultados a nível numérico e gráfico, com a possibilidade

de várias opções visuais. Neste trabalho, descreve-se o fluxograma do programa, com os

procedimentos básicos efectuados e o sistema de base de dados utilizado.

Para validar e demonstrar as capacidades do programa, recorre-se a alguns exemplos que serão

sujeitos a análise de erro de modo a analisar a convergência das soluções do método dos elementos

finitos e a validade das soluções obtidas. Exemplifica-se também a utilização conjunta das duas

metodologias.

2

2 A condição axissimétrica

As estruturas axissimétricas ou de revolução, são estruturas que podem ser representadas por uma

secção transversal que contém um eixo de revolução. Estas estruturas podem assim ser geradas

rodando a secção transversal 360° segundo o eixo de revolução.

Figura 2.1 - Estrutura axissimétrica com a secção transversal que a gerou por rotação segundo o eixo de revolução [1]

Dada a existência de um eixo de revolução, os problemas axissimétricos podem ser definidos em

coordenadas cilíndricas. A seguinte imagem representa o sistema de coordenadas globais para um

problema axissimétrico:

Figura 2.2 - Sistema de coordenadas cilíndricas (r,z,θ) [1]

Considerando uma estrutura e carregamento axissimétricos, o campo de deslocamentos depende

apenas das coordenadas r e z, e tem igualmente apenas duas componentes uma vez que não existe

deslocamento circunferencial.

Eixo de revolução

Secção transversal

Eixo de revolução

Ponto (r,z,θ)

3

Embora a casca axissimétrica se encontre definida no espaço tridimensional, esta casca pode ser

definida no plano que contém o eixo de revolução, como uma secção transversal de revolução. Uma

vez que neste trabalho apenas será abordado cascas axissimétricas sujeitas a um carregamento

axissimétrico, o problema axissimétrico pode ser totalmente representado no espaço bidimensional,

daqui em diante designado por Plano de Simplificação Axissimétrica (PSA) . Embora existam tensões

e deformações tangenciais perpendiculares ao plano em análise, estas tensões e deformações são

constantes ao longo da coordenada angular de revolução, θ, devido à simplificação axissimétrica e

às condições expostas anteriormente.

Dado o problema definido no espaço bidimensional, é possível efectuar uma equivalência entre

elementos de casca axissimétrica com elementos de barra assentes em fundação elástica ao longo

do seu domínio, sendo esta rigidez definida segundo o eixo radial r. Esta fundação elástica é

explicada devido a um “efeito de arco” existente em qualquer ponto da secção transversal

axissimétrica com raio não nulo.

4

3 Fundamentos da teoria de cascas cónicas

No que se refere a estruturas axissimétricas procede-se ao estudo de uma casca cónica com um

ângulo de inclinação geral relativamente à horizontal: α. Para α=0º ou α=180º o elemento comporta-

se como uma casca circular, que engloba também as lajes e placas circulares. Para α=90º ou α=270º

o elemento comporta-se como uma casca cilíndrica. Se α for diferente dos ângulos referidos o

elemento comporta-se como uma casca cónica.

Considera-se dois sistemas de coordenadas: um sistema global e um sistema local que depende do

elemento considerado. No PSA, os elementos de casca tronco-cónica, daqui em diante

frequentemente designados simplesmente por elementos de casca cónica, podem ser representados

por troços rectos. Assim, um elemento de casca cónica pode ser definido por um ponto inicial (1) e

um ponto final (2). Seguidamente, apresenta-se os eixos correspondentes aos dois sistemas de

coordenas:

Figura 3.1 - Sistema de coordenadas local

Figura 3.2 - Sistema de coordenadas global

Os eixos dos sistemas de coordenadas na Figura 3.1 e Figura 3.2 indicam o sentido positivo dos

deslocamentos. As seguintes figuras representam o sentido positivo dos esforços:

Figura 3.3 - Sentido positivo dos esforços com representação no PSA

Figura 3.4 - Sentido positivo dos esforços representados no plano horizontal

5

3.1 Hipótese Cinemática

Relativamente à hipótese cinemática, adopta-se a teoria de Kirchhoff que é válida para elementos de

casca fina. Esta teoria considera que as secções do elemento permanecem planas na posição

deformada, desprezando a deformação por esforço transverso. A seguinte imagem representa a

deformação generalizada de uma fibra:

Figura 3.5 - Deformação generalizada de uma fibra

Estado deformado

Estado indeformado

A variável representa o ângulo de distorção da fibra representada a vermelho na Figura 3.5, e pode

ser obtido através da seguinte expressão:

Segundo a teoria de Kirchhoff o ângulo de distorção é nulo, ou seja: . Esta situação está

representada na Figura 3.5 pela fibra azul.

Nestas condições é possível obter a expressão do campo de rotações do elemento:

(3.1)

É possível constatar através desta expressão que a rotação do elemento depende apenas dos

deslocamentos transversais.

6

3.2 Deslocamentos e deformações

Seguidamente, apresentam-se o campo de deslocamentos de uma casca cónica no PSA, segundo o

eixo de coordenadas locais e globais:

Figura 3.6 - Campo de deslocamentos em coordenadas globais

Figura 3.7 - Campo de deslocamentos em coordenadas locais

Na Figura 3.6 e Figura 3.7, v representa o deslocamento de uma secção do elemento de casca e u

representa o deslocamento de um ponto na fibra da secção.

O campo de deslocamentos em coordenadas locais (Figura 3.7) de um ponto de uma secção pode

ser obtido através da seguinte expressão:

(3.2)

A deformação segundo a direcção longitudinal de um elemento no PSA é a seguinte:

(3.3)

(3.4)

Em que:

A componente de extensão axial é:

(3.5)

7

A componente de curvatura é:

(3.6)

As deformações na direcção tangencial podem ser expressas pela seguinte expressão:

(3.7)

A determinação geométrica da deformação segundo a direcção tangencial será efectuada obtendo-

se as suas componentes separadamente, nomeadamente, a extensão tangencial e a curvatura

tangencial . A extensão tangencial pode ser deduzida a partir das seguintes imagens:

Figura 3.8 - Extensão tangencial representada no PSA

Figura 3.9 - Extensão tangencial projectada no plano perpendicular ao eixo de axissimetria

A extensão de um elemento pode ser obtida através da relação entre a variação do comprimento do

elemento e o seu comprimento inicial:

A componente de extensão tangencial é:

(3.8)

8

A dedução da variação da curvatura tangencial é efectuada com o auxílio da seguinte figura:

Figura 3.10 - Curvatura tangencial

A parte superior da Figura 3.10 representa um deslocamento transversal numa determinada secção

do elemento representado no PSA, a parte inferior representa um rebatimento do elemento de casca

cónica não devendo ser confundido com uma projecção no plano perpendicular ao eixo de

axissimetria. Trata-se pois de um plano de revolução que contém o eixo x, ou seja, uma

“planificação” do elemento de casca cónica.

9

Tendo em conta a Figura 3.10 a curvatura pode ser obtida através da seguinte expressão:

em que t representa o eixo tangencial, pretendendo-se obter a curvatura para t=0 que é a secção

considerada. Deve-se ter em atenção que a obtenção da curvatura através deste método é apenas

possível se o vector t estiver contido no plano da casca cónica, razão pela qual se efectuou o

rebatimento da casca. Ou seja, o vector t deve estar contido no plano da casca para representar

correctamente o caminho da dupla derivada.

Primeiro é necessário obter o campo de deslocamentos transversais em função do eixo tangencial t.

Dada a condição axissimétrica o campo deslocamentos transversais segundo a direcção tangencial

pode ser facilmente obtido através de um raio adulterado rt, sendo o deslocamento transversal nessa

secção: .

Finalmente, é possível obter a curvatura tangencial:

para

Curvatura tangencial:

(3.9)

10

Seguidamente, apresenta-se todas as expressões das extensões e curvaturas deduzidas para a

casca cónica:

(3.10)

3.3 Equações de equilíbrio

3.3.1 Domínio

Seguidamente, determinam-se as equações diferenciais de equilíbrio para elementos de casca

cónica. Na Figura 3.11 está representado um elemento infinitesimal de uma casca cónica com todos

os esforços internos e carregamentos considerados.

Figura 3.11 - Elemento infinitesimal de uma casca cónica

11

Os vectores px e py representam os carregamentos distribuídos em área do elemento infinitesimal

segundo os eixos x e y respectivamente. Estes carregamentos são considerados constantes e

uniformes na secção infinitesimal. Os vectores localizados nas fronteiras do elemento representam

esforços uniformemente distribuídos ao longo das fronteiras onde se encontram representados.

Tratando-se de um elemento de casca compreende-se que este se encontre definido no espaço

tridimensional, pelo que deverão existir 6 equações de equilíbrio correspondentes a 6 graus de

liberdade. No entanto, dada a condição axissimétrica, e uma vez que todos os problemas podem ser

representados no PSA, o número de graus de liberdade e o número de equações de equilíbrio passa

a ser 3, sendo as restantes equações trivialmente satisfeitas.

Dada a condição de axissimetria da estrutura e carregamento é possível deduzir que as

componentes radiais e são constantes.

Todas as componentes necessárias para determinar as equações de equilíbrio encontram-se na

Figura 3.12, que representa o elemento infinitesimal no PSA. No entanto, existem algumas

componentes representadas na Figura 3.12 que são obtidas a partir da Figura 3.11 e Figura 3.13.

Figura 3.12 - Elemento infinitesimal no plano de simplificação axissimétrica

Figura 3.13 - Corte do elemento infinitesimal no plano perpendicular ao eixo de revolução

Tendo em conta o elemento infinitesimal, o carregamento e esforços considerados, adoptam-se os

seguintes pressupostos para a obtenção das equações diferenciais de equilíbrio:

Despreza-se infinitésimos de 2ª ordem;

O seno de ângulos infinitesimais é igual ao próprio ângulo.

12

A partir da Figura 3.11 é possível relacionar o ângulo com da seguinte forma:

Através desta relação de ângulos é possível determinar a parcela da componente que vai

contribuir no PSA:

Recorrendo à Figura 3.13, obtém-se a parcela da componente que vai contribuir no PSA:

Finalmente, obtém-se as equações de equilíbrio através do somatório de forças ou momentos

segundo os graus de liberdade do problema. As equações de equilíbrio serão escritas em função da

coordenada global r:

(3.11)

(3.12)

(3.13)

13

Ao adoptar-se a hipótese de Kirchhoff, despreza-se a deformação por esforço transverso, pelo que

as equações de equilíbrio podem ser reescritas substituindo a variável da expressão (3.13) na

(3.12), resultando na seguinte equação de equilíbrio:

(3.14)

Nestas condições as equações de equilíbrio são constituídas pelas expressões (3.11) e (3.14).

3.3.2 Fronteira estática

Seguidamente, apresenta-se as condições de fronteira estática:

As variáveis são valores prescritos na fronteira, é o factor multiplicativo que define o

sinal do esforço e está associado à faceta do elemento, no caso da faceta negativa (nó inicial)

, no caso da faceta positiva (nó final) .

3.3.3 Fronteira cinemática

Seguidamente, apresenta-se as condições de fronteira cinemática:

As condições de fronteira cinemáticas permitem a definição de deslocamentos impostos na fronteira

ou restrições dos deslocamentos e consequentemente dos graus de liberdade quando é assumido

um valor nulo.

14

3.4 Relações constitutivas

Face à pequena espessura da casca axissimétrica pode admitir-se que o estado de tensão em cada

ponto é duplo, correspondendo a tensão principal nula à direcção y. Admitindo-se que os elementos

em estudo têm um comportamento elástico-linear, é possível estabelecer as seguintes relações entre

tensões e deformações recorrendo à Lei de Hooke:

(3.15)

(3.16)

Tendo em conta a hipótese cinemática adoptada e em particular as expressões (3.3) e (3.7) obtém-

se as seguintes relações constitutivas:

Relações entre esforços de membrana e extensões:

(3.17)

(3.18)

Relações entre esforços de flexão e curvaturas:

(3.19)

(3.20)

15

3.5 Dedução das equações de equilíbrio através do PTV

Outra forma possível de deduzir as equações de equilíbrio, além da forma física efectuada na

subsecção 3.3.1, deriva do PTV. O princípio dos trabalhos virtuais para um corpo deformável afirma

que esse corpo está em equilíbrio se e só se o trabalho realizado pelas forças externas e pelas

forças internas ao longo de qualquer campo de deslocamentos virtuais for compatível com as

ligações:

(3.21)

A seguinte imagem identifica os deslocamentos virtuais nodais, as forças externas prescritas na

fronteira e o carregamento a que está sujeito um elemento:

Figura 3.14 - Deslocamentos virtuais, forças externas e carregamento

Na dedução das equações de equilíbrio através do PTV, em geral o equilíbrio no domínio aparece na

forma final como o integral no elemento segundo a variável x, os restantes elementos dizem respeito

ao equilíbrio na fronteira.

Através da Figura 3.14, é possível obter o trabalho das forças externas:

(3.22)

Os termos a superescrito, 0 e L, indicam a coordenada x correspondente à secção inicial (nó 1) e

final (nó 2) do elemento, respectivamente.

16

O trabalho das deformações internas é dado por:

(3.23)

Recorrendo às expressões (3.3), (3.7) e (3.10), obtém-se:

(3.24)

Desenvolvendo a parcela: , obtém-se:

Recorrendo à secção 3.4 identifica-se os termos correspondentes aos esforços:

Utilizando a primitivação por partes obtém-se a seguinte expressão:

(3.25)

Desenvolvendo a parcela: , obtém-se:

Recorrendo à secção 3.4 identifica-se os termos correspondentes aos esforços:

Utilizando a primitivação por partes obtém-se a seguinte expressão:

(3.26)

17

Agrupando as expressões (3.22), (3.25) e (3.26) é possível encontrar a forma final do PTV expressa

em termos de esforços:

(3.27)

Como o campo de deslocamentos virtuais é arbitário, conclui-se que os termos entre parênteses

rectos têm que ser identicamente nulos, recuperando-se assim as equações de equilíbrio no domínio

e na fronteira.

18

4 Implementação de soluções analíticas

Na definição da formulação de todos os elementos utiliza-se o sistema de coordenadas locais,

excepto para o elemento de anel que envolve apenas um nó, neste elemento utiliza-se o sistema de

coordenadas globais.

A metodologia das soluções analíticas é baseada no método dos deslocamentos, e para todos os

elementos excepto o anel, é utilizada uma solução particular. O elemento de anel pode ser definido

por apenas um nó, e os elementos de: casca cilíndrica, casca em anel, incluindo o caso particular da

casca circular devem ser definidos por dois nós um inicial (nó 1) e um final (nó 2). A cada nó estão

associados três graus de liberdade no PSA, um deslocamento vertical, horizontal e uma rotação.

Também podem estar associados forças por unidade de comprimento circunferencial com o mesmo

sentido e direcção dos deslocamentos. A nível elementar a equação do método dos deslocamentos

pode ser escrita da seguinte forma:

(4.1)

Em que:

– Matriz de rigidez elementar;

– Vector de deslocamentos nodais elementar;

– Vector de forças nodais elementar;

– Vector elementar de forças aplicadas nos nós;

– Vector elementar de forças de fixação das cargas de vão.

4.1 Elemento de anel

O elemento de anel pode ser representado por apenas um nó no PSA e como solução analítica

utiliza apenas um nó.

4.1.1 Matriz de rigidez

Trata-se de um elemento com apenas três graus de liberdade pelo que a matriz de rigidez elementar

terá essa dimensão [6]:

(4.2)

(4.3)

Em que As e I são a área e a inércia da secção de revolução do elemento de anel, respectivamente.

A variável r é constante e corresponde ao raio do nó onde está definido o elemento.

19

4.1.2 Carregamento

Como o elemento é definido por apenas um nó no PSA, não existem cargas de vão, existindo apenas

cargas concentradas definidas no nó, que no espaço tridimensional são efectivamente cargas

linearmente distribuídas pelo perímetro do elemento de anel. Estas cargas são definidas a nível

global como especificado no capítulo 6.

4.1.3 Resultados finais

Os deslocamentos do elemento de anel correspondem ao vector de deslocamentos associado ao nó

em que está definido, estes são obtidos a partir da resolução do sistema de equações lineares do

método dos deslocamentos.

Os esforços do elemento de anel podem ser obtidos através da multiplicação do vector de forças

nodais do elemento (4.3) pelo raio a que se encontra. É possível observar que a solução dos

esforços compreende apenas o esforço de membrana tangencial e o momento flector tangencial:

(4.4)

4.2 Casca em anel

Nesta secção serão apresentadas as equações de deformação e de equilíbrio da casca em anel, que

podem ser obtidas através das expressões do capítulo 3, correspondentes ao elemento de casca

cónica, adoptando-se o caso particular de α=0°. Nesta secção também serão apresentadas soluções

analíticas para o elemento de casca em anel e para o caso particular do raio inicial (correspondente

ao nó 1) ser nulo constituindo uma casca circular.

Seguidamente, apresenta-se um elemento geral de casca em anel:

Figura 4.1 - Elemento geral de casca em anel

20

No elemento de casca em anel a variável r tem a mesma direcção e sentido que a variável x,

diferindo apenas na origem. A variável r tem origem fixa e a variável x tem origem no nó inicial do

elemento (nó 1). A dedução das soluções analíticas para o elemento de casca em anel é feita de

forma mais intuitiva se se obter as deformações e equações de equilíbrio em função da variável r.

4.2.1 Deformações

Seguidamente, apresenta-se o campo de deformações do elemento de casca em anel:

(4.5)

(4.6)

(4.7)

(4.8)

4.2.2 Equações de equilíbrio

Seguidamente, apresentam-se as equações de equilíbrio do elemento de casca em anel:

(4.9)

(4.10)

(4.11)


No que diz respeito à matriz de rigidez, efectua-se uma análise do elemento de casca em anel em

dois tópicos:

Comportamento de laje [5];

Comportamento de placa [7].

4.2.3.1 Comportamento de laje

A equação diferencial governativa correspondente à teoria de flexão de lajes circulares, com

espessura constante e adoptando a hipótese de Kirchhoff é [5]:

(4.12)

21

Para um carregamento transversal constante, , a solução da equação diferencial é [5]:

(4.13)

Sendo os Bs constantes a determinar.

Definindo-se da seguinte forma o vector {B}:

(4.14)

É possível estabelecer uma relação entre os deslocamentos nodais e {B}:

(4.15)

(4.16)

Com:

É também possível estabelecer uma relação entre as forças nodais e {B}:

(4.17)

(4.18)

Finalmente, é possível obter a matriz de rigidez da casca em anel, tendo em conta apenas o

comportamento de laje:

(4.19)

22

4.2.3.2 Comportamento de placa

O comportamento de placa considerando um carregamento nulo segundo a direcção longitudinal,

, é regido pela equação de equilíbrio (4.9) [7]. Recorrendo às relações constitutivas (3.17) e

(3.18) e às deformações (4.5) e (4.7), é possível obter a equação diferencial do deslocamento

longitudinal [7]:

Adoptando-se: e utilizando a seguinte regra de derivação em função da variável r:

, obtém-se:

Cuja solução é:

Para determinar as variáveis A e B, impõe-se as condições de fronteira dos deslocamentos nodais

segundo a direcção longitudinal:

Obtendo-se:

Finalmente, o deslocamento longitudinal da casca em anel é expresso por:

(4.20)

23

Recorrendo à relação constitutiva de (3.17), à extensão longitudinal (4.5) e tangencial (4.7) e à

expressão do deslocamento longitudinal (4.20), é possível obter o esforço de membrana longitudinal:

(4.21)

Uma vez deduzida a expressão do esforço de membrana longitudinal é possível obter as forças

nodais segundo essa direcção:

(4.22)

(4.23)

Identificando os termos que multiplicam pelos deslocamentos nodais longitudinais é possível

construir a matriz de rigidez, que contabiliza o efeito de placa da casca em anel:

(4.24)

4.2.3.3 Comportamento de casca

No que diz respeito à matriz de rigidez, para contabilizar o comportamento de casca deste elemento,

recorre-se aos termos das matrizes de rigidez definidos para o comportamento de laje (4.19) e de

placa (4.24). Uma vez que comportamento laje e placa são independentes, a matriz de rigidez

elementar do elemento de casca em anel é:

(4.25)

24

4.2.4 Vector das forças de fixação

Para determinar as forças de fixação recorre-se a uma solução particular.

Considerando um carregamento de vão transversal constante, , um carregamento de vão

longitudinal nulo, e adoptando o vector {B} nulo, obtém-se o seguinte vector de deslocamentos

nodais da solução particular [5]:

(4.26)

As forças nodais para a mesma solução particular são [5]:

(4.27)

Com:

Finalmente, o vector das forças de fixação para a solução particular adoptada é dado por:

(4.28)


Com base nas matrizes (4.16) e (4.18) para , nas equações de deformação, de equilíbrio e

relações constitutivas, define-se o vector da solução parcial de deslocamentos e esforços:

(4.29)

25

(4.30)

O vector pode ser obtido através da expressão (4.15), no entanto, a utilização de uma solução

particular para a obtenção de resultados da casca em anel, implica a subtracção dos deslocamentos

nodais obtidos através da resolução do método deslocamentos pelos deslocamentos da solução

particular:

(4.31)

No que diz respeito ao campo de deslocamentos, o deslocamento longitudinal é obtido através da

expressão (4.20), os restantes deslocamentos, transversal e rotação são obtidos através do vector

onde se adiciona os deslocamentos da solução particular com base no vector (4.26) para

.

O campo de deslocamentos da casca em anel é dado por:

(4.32)

Os esforços de membrana são obtidos com base nas equações de deformação e relações

constitutivas:

(4.33)

26

Os esforços de flexão e o esforço transverso são obtidos através do vector ao qual se subtrai

as parcelas correspondentes à solução particular com base no vector (4.27) para :

(4.34)

Com:

4.2.6 Caso particular: elemento de casca circular

É possível considerar o elemento de casca circular como um caso particular do elemento de casca

em anel atribuindo a um dos nós a coordenada de raio nula.

Nesta secção, analisa-se as implicações que ocorrem ao considerar a coordenada inicial com raio

nulo: . Examinando as soluções analíticas, é possível verificar que a variável constitui por

vezes um operador de divisão, ou é responsável pela obtenção de uma função logarítmica do valor

zero ( ). Nestes casos é necessário levantar as indeterminações, pelo que se deve efectuar uma

análise limite da variável a tender para zero por valores positivos.

Foi necessário efectuar uma análise limite às seguintes matrizes:

da expressão (4.29) para a obtenção do vector ;

da expressão (4.31) para a obtenção do vector ;

de forma a obter a matriz de rigidez elementar da casca circular

Constatou-se que só foi possível a obtenção de valores limite para a matriz de rigidez elementar

após a multiplicação dos termos da matriz de rigidez da casca em anel pelo raio do respectivo nó

correspondente ao grau de liberdade associado, como explicado na secção 6.3.

4.2.6.1 Resultados finais

Os resultados finais da casca circular podem ser obtidos como descrito para a casca em anel,

excepto quando se pretende obter os resultados para a secção de raio nulo, que neste caso

corresponde a . A metodologia utilizada na obtenção de resultados para esta secção pode

variar dependendo de determinados casos.

27

O campo de deslocamentos, na secção , pode ser sempre fornecido pelos deslocamentos

nodais correspondentes, que são obtidos a partir da resolução do sistema de equações lineares do

método dos deslocamentos:

(4.35)

Os esforços de membrana podem ser obtidos como descrito para o caso da casca em anel, mas

considerando :

(4.36)

No que diz respeito aos esforços de flexão e esforço transverso não existe valor limite no caso de

existir uma carga concentrada aplicada no nó de raio nulo, pois origina uma singularidade como é

explicado no capítulo 7. Também não existe limite no caso de existir um apoio no nó que restrinja o

deslocamento vertical, o que resulta numa reacção concentrada originando uma singularidade à

semelhança da carga concentrada.

Uma análise limite à matriz permite concluir que não existe valor limite na segunda coluna para

os esforços de flexão e na quarta coluna para os mesmos esforços e também para o esforço

transverso. Torna-se necessário analisar o vector , nomeadamente, o segundo e o quarto termo

de forma a verificar se é possível obter valores limite. Efectuando uma análise limite à matriz ,

verifica-se que o termo é obrigatoriamente nulo, no entanto, o mesmo já não acontece para o

termo . Com base na solução do esforço transverso para a casca em anel (4.34) e observando a

matriz (4.30), é possível verificar que para x=0 o esforço transverso depende apenas da

multiplicação do termo da quarta coluna de pelo termo . Sabendo que, por simetria e

equilíbrio, o esforço transverso é nulo para x=0 na ausência de carga concentrada ou apoio

concentrado no nó de raio nulo, deduz-se que nestas condições o termo é nulo, obtendo-se os

seguintes esforços:

(4.37)

28

4.3 Casca cilíndrica

Nesta secção serão apresentadas as equações de deformação e de equilíbrio da casca cilíndrica,

que podem ser obtidas através das expressões do capítulo 3, correspondentes ao elemento de casca

cónica, adoptando-se o caso particular de α=90° e raio constante. Nesta secção também serão

apresentadas as soluções analíticas para este elemento.

Seguidamente, apresenta-se um elemento geral de casca cilíndrica:

Figura 4.2 - Elemento geral de casca cilíndrica

4.3.1 Deformações

Seguidamente, apresenta-se o campo de deformações do elemento de casca cilíndrica:

(4.38)

(4.39)

(4.40)

(4.41)

4.3.2 Equações de equilíbrio

Seguidamente, apresentam-se as equações de equilíbrio do elemento de casca cilíndrica:

(4.42)

(4.43)

(4.44)

29


A análise do elemento de casca cilíndrica, no que diz respeito à matriz de rigidez, será abordada em

dois tópicos:

Fundamentos da teoria de flexão de cascas cilíndricas [5];

Implementação da solução de membrana e correcção da parcela de flexão [3].

4.3.3.1 Teoria de flexão de cascas cilíndricas

Nesta secção apresenta-se uma solução analítica considerando a casca cilíndrica apenas sujeita a

carregamentos transversais. Considera-se que a casca não está restringida a nível de

deslocamentos longitudinais e que consequentemente, segundo esta direcção, não existem esforços

de membrana.

Adoptando-se um esforço de membrana longitudinal nulo, , verifica-se o seguinte: a um

deslocamento transversal/radial de grandeza, , está associado um esforço de membrana

tangencial por unidade de comprimento circunferencial, dado pela seguinte expressão [5]:

(4.45)

Com base nas equações de equilíbrio (4.43) e (4.44), na relação constitutiva de (3.19) e na

expressão (4.45), obtém-se a seguinte equação diferencial do comportamento de flexão da casca

cilíndrica:

(4.46)

Analisando a expressão (4.46), é possível verificar que é análoga àquela que rege o comportamento

de uma viga sobre fundação elástica, a qual é dada por

(4.47)

sendo EI a rigidez de flexão da viga, k a rigidez da fundação elástica.

Uma solução geral para a equação diferencial (4.46) é [5]:

(4.48)

30

O parâmetro representa o inverso do comprimento característico:

(4.49)

As funções são adimensionais e têm as seguintes expressões:

(4.50)

Estas funções estão relacionadas com a sua primeira derivada da seguinte forma:

; ; ; ;

Definindo-se o seguinte vector :

(4.51)

É possível estabelecer uma relação entre os deslocamentos nodais e {B}:

(4.52)

(4.53)

31

É também possível estabelecer uma relação entre as forças nodais e {B}:

(4.54)

(4.55)

Finalmente, é possível obter a matriz de rigidez da casca cilíndrica, tendo em conta apenas os

deslocamentos transversais:

(4.56)

4.3.3.2 Solução de membrana e correcção da parcela de flexão

Para determinar os termos de rigidez que envolvem os graus de liberdade associados aos

deslocamentos longitudinais da casca cilíndrica, pretende-se obter um sistema com deslocamentos

nodais unitários segundo a direcção longitudinal e deslocamentos nodais nulos nos restantes graus

de liberdade. Uma solução possível consiste na sobreposição de dois sistemas de aplicação de

cargas de forma a obter apenas deslocamentos longitudinais não nulos. O sistema de sobreposição

de acções utilizado é o seguinte:

Figura 4.3 - Configuração de sistemas para a implementação da solução de membrana

Ambos os sistemas têm todos graus de liberdade nodais não restringidos. No sistema , considera-

se a casca cilíndrica sujeita apenas a esforço de membrana , isto origina não só deslocamentos

32

longitudinais, mas também transversais devido ao efeito de Poisson, no entanto, verifica-se que não

existem rotações no elemento. É necessário com o sistema determinar as forças a aplicar nos

graus de liberdade que vão anular o deslocamento transversal do sistema e manter as rotações

nodais nulas.

Seguidamente, procede-se à análise do sistema :

Verifica-se que o esforço de membrana tangencial é nulo:

Recorrendo às relações constitutivas (3.17) e (3.18), obtém-se a seguinte extensão tangencial

e longitudinal:

Uma vez obtidas as extensões é possível determinar os deslocamentos nodais longitudinais e

transversais através das seguintes expressões:

Seguidamente, procede-se à análise do sistema :

Neste sistema, pretende-se anular os deslocamentos transversais nodais do sistema , impondo a

seguinte condição:

Verifica-se que no sistema não existe esforço de membrana longitudinal:

Recorrendo às relações constitutivas (3.17) e (3.18), obtém-se a seguinte extensão

longitudinal:

33

Após a obtenção da expressão da extensão longitudinal do sistema , é possível determinar o

deslocamento longitudinal através da seguinte expressão:

Recorrendo à expressão (4.45) é possível obter o deslocamento longitudinal correspondente ao

sistema :

(4.57)

Verifica-se que as forças de fixação correspondentes ao sistema podem ser obtidas através de:

É possível observar que na expressão (4.57) é apenas necessário determinar as forças de fixação do

nó inicial (nó 1):

Tendo em conta as condições impostas ao sistema , em que: e , a expressão

pode ser reescrita da seguinte forma:

(4.58)

34

Finalmente, com a sobreposição de acções dos sistemas e , é possível obter a rigidez de

membrana segundo a direcção longitudinal:

(4.59)

Com:

(4.60)

A matriz de rigidez de um elemento de casca cilíndrica associada unicamente aos graus de liberdade

longitudinais apresenta-se de seguida:

(4.61)

Recorrendo às expressões (4.59) e (4.61) é possível determinar os termos da matriz de rigidez da

casca cilíndrica que relacionam os graus de liberdade transversais com os longitudinais:

(4.62)

A consideração dos deslocamentos longitudinais na casca cilíndrica e a introdução do seu efeito de

acoplamento nos deslocamentos transversais requer uma correcção da matriz de rigidez transversal

inicial (4.56) que despreza ou considera nulo o esforço de membrana longitudinal. Torna-se

necessário determinar a deformação longitudinal que ocorre quando um dos deslocamentos

transversais ou rotações é unitário sendo os restantes nulos.

Determina-se a deformação longitudinal recorrendo à expressão (4.57), sendo apenas necessário

determinar os deslocamentos e esforços no nó inicial:

35

A matriz de rigidez transversal corrigida, da casca cilíndrica, é:

(4.63)

Seguidamente, indica-se o deslocamento nodal com valor unitário que origina as funções , quando

os restantes deslocamentos são nulos:

(4.64)

Finalmente, apresenta-se a configuração da matriz de rigidez da casca cilíndrica, no referencial local,

com base nas matrizes (4.61), (4.62) e (4.63):

(4.65)

4.3.4 Vector das forças de fixação

É possível determinar o vector de forças de fixação da casca cilíndrica utilizando uma solução

particular para um determinado carregamento.

No que diz respeito aos esforços de membrana, uma solução particular possível é:

Esforço de membrana segundo a direcção longitudinal:

(4.66)

36

Para determinar uma solução particular do esforço de membrana tangencial pode-se efectuar dois

cortes: um perpendicular ao eixo de revolução numa determinada secção da casca cilíndrica, e outro

que contém o eixo de axissimetria. Efectuando-se somatório de forças no plano de corte horizontal é

possível obter o esforço de membrana tangencial em função do carregamento transversal, conforme

se ilustra de seguida:

Figura 4.4 – Plano de corte horizontal de uma casca cilíndrica

(4.67)

Seguidamente, determina-se uma solução particular dos deslocamentos através das extensões

longitudinais e tangenciais, recorrendo às equações constitutivas e utilizando as soluções

particulares dos esforços de membrana:

(4.68)

(4.69)

(4.70)

37

O vector de deslocamentos nodais da solução particular é:

(4.71)

As forças nodais da solução particular são:

(4.72)

Finalmente, o vector de forças fixação é:

(4.73)


Com base na equação (4.48), (4.57), nas matrizes (4.53) e (4.55), nas equações de deformação, de

equilíbrio e relações constitutivas, define-se o seguinte vector da solução parcial de deslocamentos e

esforços:

(4.74)

(4.75)

38

De forma a ser possível a utilização das soluções particulares na obtenção do campo de

deslocamentos e esforços, é necessário determinar as forças nodais existentes no elemento que

resultam da diferença dos deslocamentos nodais obtidos a partir da resolução do sistema de

equações lineares do método dos deslocamentos pelos deslocamentos da solução particular:

(4.76)

Às parcelas correspondentes aos deslocamentos no vector {B} subtrai-se igualmente os respectivos

deslocamentos da solução particular obtendo-se:

(4.77)

Finalmente, adiciona-se ao vector a solução particular dos deslocamentos, e utiliza-se a solução

particular dos esforços de membrana.

Os resultados finais para a casca cilíndrica apresentam-se de seguida:

Deslocamentos:

(4.78)

Esforços de membrana:

(4.79)

Esforços de flexão e esforço transverso:

(4.80)

39

5 Implementação do método dos elementos finitos a uma casca

cónica

O MEF é um método de resolução numérico de equações diferenciais parciais.

A implementação do MEF será efectuada no PSA, uma vez que este plano é suficiente para a

definição completa do problema. Na base da definição de um elemento finito estão os nós: inicial e

final.

Para a implementação do MEF a cascas cónicas procede-se à alteração de determinadas variáveis

para a simplificação do problema. A variável x tem a direcção do eixo do elemento, sentido do nó

inicial para o nó final, e tem início coincidente com o nó inicial. As seguintes variáveis mudam a

notação:

O MEF é baseado no método dos deslocamentos, no entanto, ao contrário das soluções analíticas

este método não utiliza soluções particulares, pois baseia-se em funções de aproximação do campo

de deslocamentos e a nível elementar a equação do método dos deslocamentos pode ser escrita da

seguinte forma:

(5.1)

Em que:

– Matriz de rigidez elementar;

– Vector de deslocamentos nodais elementar;

– Vector de forças nodais equivalentes elementar;

5.1 Funções de forma

As funções de forma estão associadas aos deslocamentos do elemento. Neste trabalho, para a

implementação do MEF, considera-se funções de forma polinomiais devido à sua simplicidade. A

solução homogénea dos deslocamentos de elementos de casca axissimétricos não segue funções

polinomiais, pelo que as funções de forma fornecem apenas uma aproximação do campo de

deslocamentos, incluindo os deslocamentos nodais.

A função de forma associada a um determinado grau de liberdade é obtida estabelecendo um valor

unitário do deslocamento pretendido e igualando os restantes a zero.

40

5.1.1 Polinómios de Lagrange

Considerando elementos unidimensionais, sujeitos apenas a esforço axial, originam-se tensões,

deformações e deslocamentos também unidimensionais que podem ser expressos segundo o próprio

eixo. Considera-se um elemento constituído por dois nós, inicial (1) e final (2), com os deslocamentos

nodais considerados:

Figura 5.1 - Elemento e graus de liberdade considerados para deslocamentos longitudinais

O campo de deslocamentos neste elemento pode ser descrito por funções de forma polinomiais do 1º

grau, tendo em conta o número de nós e condições de fronteira:

Para um comprimento geral do elemento, L, apresenta-se as funções de forma obtidas:

Figura 5.2 - Função de forma associada a u1

Figura 5.3 - Função de forma associada a u2

00,25

0,50,75

1ψ

xL

00,25

0,50,75

1ψ

xL

41

5.1.2 Polinómios de Hermite

Para determinar o campo de deslocamentos transversal e rotações apresenta-se um elemento geral

constituído por dois nós, inicial (1) e final (2), e os deslocamentos nodais considerados:

Figura 5.4 - Elemento e graus de liberdade considerados para deslocamentos transversais

O campo de deslocamentos deve garantir continuidade das funções de forma e as suas derivadas,

tendo como base funções polinomiais do terceiro grau, nomeadamente, os polinómios de Hermite. As

funções de forma têm por conseguinte a seguinte expressão geral:

w(x)= ψw(x)×d

Para um comprimento geral do elemento, L, as condições que permitem a definição das funções de

forma são as seguintes:

:

Figura 5.5 - Função de forma associada a w1

:

Figura 5.6 - Função de forma associada a φ1

0

0,25

0,5

0,75

1

ψ

xL

ψ

xL

42

:

Figura 5.7 - Função de forma associada a w2

:

Figura 5.8 - Função de forma associada a φ2

O campo de rotações é determinado a partir dos deslocamentos nodais do elemento através de:

φ(x)= ψφ(x)×d

As funções de forma do campo de rotações podem ser obtidas como é descrito na secção 3.1:

0

0,25

0,5

0,75

1

ψ

xL

ψxL

43

5.2 Elemento finito de casca axissimétrica

Seguidamente, apresenta-se um elemento finito de casca axissimétrica representado no PSA, com

todos os deslocamentos considerados. Trata-se de um elemento recto unido por um nó inicial (1) e

um nó final (2) com um ângulo de inclinação geral com a horizontal α:

Figura 5.9 - Elemento de finito de casca cónica representado no PSA

A posição de uma secção do elemento relativamente à coordenada r em função da coordenada x é a

seguinte:

Onde:

5.2.1 Campo de deslocamentos

Recorrendo à secção 5.1 é possível obter o campo de deslocamentos para um elemento finito de

casca axissimétrico constituído por dois nós apresentando cada um deles 3 graus de liberdade: dois

deslocamentos e uma rotação.

44

O campo de deslocamentos para o elemento finito de casca axissimétrica é o seguinte:

(5.2)

Apresentando as expressões de todas as funções de forma obtém-se:

(5.3)

5.2.2 Campo de deformações

O campo de deformações é definido pelo seguinte vector:

(5.4)

Este vector pode ser obtido através da seguinte expressão:

(5.5)

A matriz , por sua vez, tem a seguinte expressão:

(5.6)

A matriz [A] é um operador diferencial, sendo obtido através das expressões (3.10) e configuração do

vector de deformações, resultando na seguinte matriz:

(5.7)

45

A matriz toma a seguinte forma:

(5.8)

Com:

5.2.3 Campo de tensões generalizado

O campo de tensões generalizado é definido pelo seguinte vector agrupando os esforços de

membrana e de flexão:

(5.9)

Este vector pode ser obtido através da seguinte expressão:

(5.10)

A matriz diz respeito às relações constitutivas que relacionam as tensões com as deformações.

Esta matriz pode ser obtida recorrendo à secção 3.4, e tem a seguinte forma:

(5.11)

O esforço transverso pode ser obtido por equilíbrio através da expressão (3.13):

(5.12)

46


Nesta dissertação, não serão abordados elementos de casca assentes em fundação elástica, pelo

que a matriz de rigidez elementar é obtida através da seguinte expressão [2]:

(5.13)

Dada a condição de axissimetria a integração no domínio do elemento de casca, , pode ser

expresso em função da secção de revolução multiplicada pelo perímetro correspondente à

coordenada r de um ponto da secção. Logo, a matriz de rigidez elementar pode ser rescrita da

seguinte forma:

(5.14)

5.2.5 Vector de forças nodais equivalentes

No âmbito desta dissertação adoptaram-se os seguintes pressupostos:

Não se considera deformações iniciais;

Considera-se que uma secção de revolução representada por um elemento finito é

constituída por um material homogéneo, de peso volúmico constante e uniforme;

A espessura da casca de um elemento finito é constante no seu domínio.

O vector de forças nodais equivalentes de um elemento tem a seguinte expressão:

(5.15)

A integração no domínio pode ser rescrita em função da secção de revolução:

(5.16)

O vector descreve o carregamento no domínio, considerando um carregamento longitudinal

constante e um carregamento transversal linearmente variável, que no espaço tridimensional

representam carregamentos de área:

(5.17)

47

Finalmente, o vector de forças nodais equivalentes de um elemento é:

(5.18)

5.2.6 Obtenção dos campos de deslocamentos e tensões

A obtenção dos campos de deslocamentos e tensões é feita para cada elemento, onde é constituído

um vector de deslocamentos dos respectivos nós que definem o elemento. Uma vez que esse vector

de deslocamentos se encontra no sistema de coordenadas global é necessário fazer uma

transformação de coordenadas para o sistema local pois os campos de deslocamentos e tensões

foram definidos nesse sistema.

Finalmente, os deslocamentos num determinado ponto do elemento são obtidos como descrito na

subsecção 5.2.1. Como estes deslocamentos vêm no referencial local, é necessário aplicar uma

transformação linear de coordenadas para o sistema global, no entanto, deve-se ter em atenção que

a matriz de transformação linear de coordenadas diz respeito apenas a um ponto. O campo de

tensões generalizado é obtido como descrito na subsecção 5.2.3, e neste caso não é necessário

efectuar uma transformação de coordenadas uma vez que se pretende os esforços no referencial

local do elemento.

5.3 Regra de quadratura de Gauss

Para o cálculo da matriz de rigidez, é necessário efectuar a integração no domínio do elemento.

Como se pode observar na matriz , as funções a integrar não são polinomiais pois o raio é um

operador de divisão em algumas parcelas. As primitivas se possíveis não são fáceis de determinar,

pelo que se optou por aplicar um método numérico de integração: a quadratura de Gauss.

A regra de quadratura de Gauss fornece uma aproximação do valor do integral de uma função,

obtido pelo somatório de pesos aplicados a valores da função em pontos específicos, dentro do

intervalo de integração.

Uma regra de quadratura de Gauss de n pontos produz o resultado exacto de uma função polinomial

de grau 2n-1 ou menor, utilizando as abcissas xi e pesos wi adequados. Os pontos utilizados na

regra correspondem a raízes de um polinómio pertencente a uma classe de polinómios ortogonais.

48

Por convenção, a regra de quadratura de Gauss é aplicada a um intervalo de integração [-1,1] e tem

a seguinte expressão:

(5.19)

Uma vez que o intervalo de integração [-1,1] não é apropriado para a formulação definida, é

necessário generalizar o intervalo de integração de uma função e alterar a função de integração para

ser compatível com o intervalo convencional da quadratura de Gauss:

(5.20)

Aplicando a regra de quadratura de Gauss, obtém-se:

(5.21)

A integração definida no intervalo de um elemento [0,L] aplicando a regra de quadratura de Gauss é

a seguinte:

(5.22)

5.3.1 Análise de erro associado ao número de pontos de Gauss

Uma das opções a ser definida pelo utilizador é o número de pontos de Gauss a serem utilizados na

integração da matriz de rigidez elementar do elemento de casca cónica do MEF. Nesta secção,

procede-se a uma análise do erro da integração numérica da matriz de rigidez em função do número

de pontos de Gauss, com o objectivo de determinar uma ordem de grandeza do número de pontos

de Gauss aconselhável na utilização do programa.

Para avaliar o erro da integração, recorre-se a um exemplo de apenas um elemento de casca cónica,

com ângulo em relação à horizontal diferente de 0 ( ) de modo a minimizar a eliminação de

parcelas de cálculo da integração da matriz de rigidez. Para evitar problemas de não convergência

de termos da matriz de rigidez, nenhum dos nós que definem o elemento se encontram no eixo

axissimétrico. O exemplo utilizado será constituído por apenas um elemento finito, sendo necessário

calcular apenas uma matriz de rigidez elementar no sistema de coordenadas locais.

49

O elemento apresenta a configuração definida na Figura 5.9 da secção 5.2, onde se atribuiu os

seguintes valores:

Coordenadas: ; ; ;

Espessura:

Módulo de elasticidade:

Coeficiente de Poisson:

Com estes dados é possível calcular:

Comprimento do elemento:

;

Esbelteza do elemento:

A solução de referência será obtida através do software Mathematica, utilizando os dados definidos.

Este software possibilita a obtenção da solução exacta do integral tendo em conta o grau de precisão

com o qual trabalha desde que haja convergência dos integrais.

Todos os termos da matriz de rigidez serão analisados, tendo em consideração que os termos

simétricos resultarão na mesma informação.

Seguidamente, apresenta-se um quadro com a análise do erro relativo das várias componentes da

matriz de rigidez, em função do número de pontos de Gauss utilizados na integração:

Nº de pontos de Gauss

Erro relativo (%)

Máximo Médio Mínimo

1 174,961 61,776 1,327

2 34,257 11,586 0,213

3 5,516 1,283 3,24E-02

4 0,824 0,192 4,84E-03

5 0,122 2,84E-02 7,17E-04

6 1,80E-02 4,18E-03 1,06E-04

7 2,64E-03 6,14E-04 1,55E-05

8 3,87E-04 9,01E-05 2,28E-06

9 5,67E-05 1,32E-05 3,33E-07

10 8,31E-06 1,93E-06 4,88E-08

15 5,55E-10 1,29E-10 3,52E-12

20 1,32E-12 3,78E-13 0

30 3,65E-12 6,75E-13 0

Tabela 5.1 - Erro relativo dos termos da matriz de rigidez em função do nº de pontos de Gauss

50

Apresenta-se também o gráfico respectivo, com a análise do erro relativo, em função do número de

pontos de Gauss, em escala logarítmica:

Figura 5.10 - Erro relativo do nº de pontos de Gauss

É possível observar que a partir de 4 pontos de Gauss o erro relativo já é bastante baixo, sendo este

inferior a 1%. A partir de 8 pontos de Gauss o erro é relativamente insignificante e a partir de 20

pontos de Gauss a diferença é quase inexistente e verifica-se que o erro mínimo é 0% o que significa

que alguns termos da matriz de rigidez são iguais no grau de precisão definido. Considera-se que a

partir dos 20 pontos de Gauss a convergência é atingida pois com 30 pontos a ordem de erro é a

mesma, apresentando mesmo um erro mesmo ligeiramente superior, que está relacionado com o

grau de precisão utilizado.

Deve-se ter em consideração que é possível fazer deduções a partir deste problema adquirindo uma

ideia da grandeza do número de pontos de Gauss a utilizar na integração, mas a informação obtida

não pode ser rigorosamente generalizada para diferentes tipos de problemas. Em geral, cada

problema diferente apresentará erros relativos diferentes para um número não exagerado de pontos

de Gauss, embora possa existir uma correlação.

0

20

40

60

80

100

120

140

160

180

200

1 10 100

Er (%)

Número de pontos de Gauss

Erro relativo

Máximo

Médio

Mínimo

51

6 Sistema global

Este capítulo aborda a incorporação das soluções analíticas e do elemento finito apresentado num

modelo estrutural.

6.1 Transformação Linear de Coordenadas

A transformação linear de coordenadas diz respeito a todos os elementos excepto o elemento de

anel das soluções analíticas, uma vez que este contém apenas um nó e se encontra definido sempre

no sistema global de coordenadas.

Dadas as possibilidades de direcção e sentido dos elementos de casca axissimétricos, e uma vez

que ao nível do elemento são utilizadas as coordenadas locais, torna-se necessário efectuar uma

transformação de coordenadas para um sistema global. Como o problema axissimétrico pode ser

expresso no PSA a transformação de coordenadas é feita nesse plano.

Figura 6.1 - Transformação linear entre coordenadas locais e coordenadas globais

A transformação de coordenadas globais para coordenadas locais num ponto é a seguinte:

(6.1)

(6.2)

52

Como todos os elementos, excepto o de anel, são definidos por dois nós com três graus de liberdade

cada, a matriz de transformação de coordenadas é a seguinte:

(6.3)

Com a matriz de transformação é possível transformar o sistema de coordenadas das matrizes de

rigidez e dos vectores de forças nodais para um sistema global, tornando possível a interacção de

elementos com direcções e sentidos diferentes.

Transformação da matriz de rigidez:

(6.4)

Transformação do vector de forças nodais:

(6.5)

A resolução do sistema de equações em coordenadas globais resultará na obtenção de

deslocamentos no mesmo sistema de coordenadas. No entanto, como o campo de deslocamentos e

esforços são definidos no sistema de coordenadas locais, é necessário proceder a uma

transformação de coordenadas dos deslocamentos.

Vector de deslocamentos em sistema de coordenadas locais:

(6.6)

A obtenção do campo de deslocamentos, numa determinada secção, é feita no sistema de

coordenadas locais e é útil transformar os deslocamentos obtidos para um sistema de coordenadas

globais ao nível do ponto:

(6.7)

53

6.2 Definição do carregamento nodal

No que diz respeito à formulação concebida, no PSA, as cargas concentradas devem ser definidas

nos nós, pelo que a estrutura deverá ser discretizada em função da localização das cargas

concentradas, para além de outras condicionantes do problema, como por exemplo: diferentes

materiais, diferentes secções, diferentes carregamentos de vão, restrições, etc.

As cargas concentradas são definidas segundo os graus de liberdade dos nós no sistema global de

coordenadas. Dada a condição axissimétrica qualquer carga definida num nó com raio não nulo é no

espaço tridimensional uma carga uniformemente distribuída ao longo do raio que contém o nó,

apenas as cargas definidas para um nó com raio nulo são cargas concentradas. O carregamento

nodal não é definido para um elemento, mas é definido directamente para um determinado nó ni:

Para r=0:

Para r≠0:

6.3 Reunião de elementos

Os elementos definidos pela metodologia das soluções analíticas são compatíveis com o elemento

finito de casca cónica do MEF, pois ambas as metodologias se baseiam no método dos

deslocamentos. No entanto, como se pode observar na formulação do MEF, a matriz de rigidez e o

vector de forças nodais equivalentes são obtidos por integração do elemento no domínio de

revolução, o mesmo já não acontece para as soluções analíticas. Para a metodologia das soluções

analíticas ser compatível com a formulação do MEF é necessário multiplicar o grau de liberdade

correspondente ao respectivo termo da matriz de rigidez e do vector de forças de nodais pelo factor

, em que corresponde ao raio do nó associado ao grau de liberdade. O procedimento descrito

corresponde à integração do elemento ao longo da coordenada de revolução, .

Para resolver uma estrutura sujeita a um carregamento composta por vários elementos de casca

axissimétrica, é necessário proceder a uma assemblagem para construir a matriz de rigidez global e

do vector de forças nodais global a partir das matrizes e vectores elementares.

Utilizou-se um sistema de assemblagem com base em índices numéricos. Dada a definição da

estrutura, e constituição da base de dados dos nós, considera-se que cada nó tem três graus de

liberdade subtraído pelo número de restrições que pode ser atribuído a cada grau de liberdade. O

sistema de índices é construído segundo uma ordem definida pelos nós, e é apenas atribuído índices

a graus de liberdade, sendo que os deslocamentos restringidos não fazem parte do sistema de

equações. A cada grau de liberdade corresponde um índice próprio e com este sistema de índices é

feito o espalhamento da matriz de rigidez e do vector de forças nodais. As cargas concentradas

54

dizem respeito a um determinado grau de liberdade identificado por um índice que caso não esteja

restringido é directamente adicionado ao vector de forças nodais global.

Seguidamente, apresenta-se o processo de assemblagem adoptado:

1. Construção da base de dados dos nós e sistema de índices a partir da respectiva

identificação dos graus de liberdade não restringidos associados a um nó;

2. Espalhamento directo das cargas concentradas no vector de forças nodais global através da

identificação do grau de liberdade;

3. Espalhamento directo da matriz de rigidez dos elementos de anel na matriz de rigidez global,

através do sistema de índices associado ao nó;

4. Construção da matriz de rigidez e vector de forças nodais elementar;

5. Construção da matriz de transformação linear de coordenadas;

6. Transformação do sistema de coordenadas da matriz de rigidez e vector de forças nodais

para um sistema global;

7. Definição de um vector de índices elementar com base na identificação dos nós que

constituem o elemento;

8. Espalhamento da matriz de rigidez e do vector de forças nodais;

9. Repetição do passo 4 ao 7 até percorrer todos os elementos, excepto o elemento de anel,

que constituem a estrutura e obter a matriz de rigidez e vector de forças nodais global.

Após a finalização do processo de assemblagem é possível resolver o seguinte sistema de equações

lineares:

Com:

– Matriz de rigidez global;

– Vector de forças nodais global;

– Deslocamentos nodais da estrutura dos graus de liberdade não restringidos no sistema de

coordenadas globais;

O sistema de equações lineares é resolvido recorrendo ao método de eliminação de Gauss de forma

a obter o vector de deslocamentos d.

55

7 Singularidades

7.1 Nós no eixo de simetria

Avaliando a formulação e as metodologias adoptadas para a resolução de estruturas axissimétricas

tirando partido da simplificação axissimétrica, é possível verificar que existirão alguns problemas de

singularidades. A maioria dos problemas de singularidades diz respeito a nós definidos com raio

nulo. Este problema é evidente ao contemplar certos passos da formulação onde o raio é um

denominador ou é responsável pela obtenção de uma função logarítmica do valor zero ( ). Estas

situações são comuns na obtenção de determinadas parcelas da matriz de rigidez e de resultados.

Como condição axissimétrica de estrutura e carregamento, os nós definidos no eixo de axissimetria,

ou seja, com raio nulo, podem ter deslocamentos apenas segundo o eixo de axissimetria, z, pelo que

não podem ter deslocamentos radiais, segundo r, pois o resultado seria a criação de uma abertura

inexistente na estrutura. Também não podem ter rotação, , pois é necessário respeitar a condição

de axissimetria e continuidade da estrutura.

Em ambas as metodologias foi necessário efectuar análises limite para o caso de elementos com nós

definidos no eixo de axissimetria, para tal utilizou-se a ferramenta Mathematica. No entanto, existem

problemas de singularidades, que não têm valor limite, como é o caso da aplicação de cargas

concentradas na estrutura, ou apoios concentrados que de forma semelhante causam uma reacção

concentrada. As cargas concentradas são apenas possíveis se aplicadas segundo o eixo de

revolução devido à condição axissimétrica, isto porque o único grau de liberdade permitido no eixo de

revolução é segundo o próprio eixo. Analisando, por equilíbrio, a situação de uma carga concentrada,

caso se efectue um corte de raio constante numa secção do elemento, à medida que o raio do corte

diminui será necessário mais esforço transverso distribuído lateralmente na face do corte, de modo

equilibrar a carga concentrada. Na situação limite, quando o raio do corte efectuado ao elemento

tende para 0, já não existe área para equilibrar a carga pelo que o esforço resultante tende para

infinito.

56

7.2 Método dos elementos finitos

Nesta secção, descreve-se a análise limite efectuada para o caso do método dos elementos finitos.

Tendo em conta a condição axissimétrica de um elemento com um nó de raio nulo, no qual, como já

foi referido, só pode existir deslocamentos segundo o eixo de axissimetria, apresenta-se duas figuras

com os graus de liberdade desse nó:

Figura 7.1 - Graus de liberdade em sistema de coordenadas locais

Figura 7.2 - Graus de liberdade em sistema de coordenadas globais

No referencial global, tendo em consideração os deslocamentos impedidos, devido à condição

axissimétrica, é possível efectuar a seguinte dedução:

e

No referencial local, para , é possível perceber que os deslocamentos e têm projecção

segundo ambos os eixos z e r, pelo que os dois graus de liberdade serão afectados por uma parcela

que tende para infinito na matriz de rigidez. O grau de liberdade é equivalente ao grau de

liberdade . Dada a utilização de métodos numéricos e uma vez que a integração pela quadratura

de Gauss não utiliza o valor da função para a aplicação da transformação de coordenadas

para um sistema global da matriz de rigidez vai cancelar algumas parcelas que supostamente

tendem para infinito, permanecendo apenas os dois termos que tendem para infinito indicados para o

referencial global.

Uma vez que pela condição axissimétrica, num ponto com , e devem estar impedidos, ou

seja, com deslocamento nulo, as parcelas de rigidez que envolvem qualquer um dos dois graus de

liberdade mencionados não vão comparecer na matriz de rigidez global, e não farão parte do sistema

57

de equações a ser resolvido. Por conseguinte, os termos de rigidez que tendem para infinito nunca

serão utilizados e a matriz de rigidez global só terá parcelas com valores que têm limite.

Todas estas conclusões foram testadas e comprovadas a nível prático recorrendo a um software de

álgebra computacional, o programa Mathematica.

Para a obtenção do campo de tensões generalizado efectuou-se também uma análise das

singularidades. O campo de tensões generalizado é obtido através da expressão (5.9). Como já foi

referido anteriormente, a matriz contém como operador de divisão, o raio, pelo que se deverá

analisar o caso em que este é nulo.

Os deslocamentos obtidos após a resolução do sistema de equações vêm em coordenadas globais.

No entanto, os esforços estão definidos em coordenadas locais, mas para a simplificação do

problema utilizou-se a propriedade associativa do produto de matrizes de forma a poder utilizar os

deslocamentos no sistema de coordenadas globais, uma vez que se poderá tirar partido de ter

ambos os deslocamentos e nulos, resultando na seguinte expressão:

Ao produto resultante da matriz pela matriz de transformação de coordenadas , surge uma

nova matriz onde se deve definir o valor limite para por valores positivos de . Nesta matriz é

possível eliminar as colunas correspondentes aos deslocamentos nulos, uma vez que

inevitavelmente serão anuladas pelo produto de parcelas do vector de deslocamentos que têm valor

0. Mais uma vez utilizando o software Mathematica, obteve-se o valor limite da matriz resultante de

, onde se verificou que de facto as parcelas que tendem para infinito correspondem às

colunas eliminadas.

Seguidamente, apresenta-se as matrizes utilizadas, obtidas através do valor limite, para o cálculo do

campo de tensões generalizado num ponto com raio nulo:

Quando é o nó inicial que tem raio nulo, :

58

Quando é o nó final que tem raio nulo, :

É possível observar que o campo de tensões generalizado, à excepção do esforço transverso, não

apresentará problemas de singularidades, a nível de cálculo numérico.

O esforço transverso, como se pode observar pela expressão (3.13), contém o raio também como

operador de divisão, pelo que irá apresentar uma singularidade para . Após uma análise limite

relativamente à variável r tendendo para 0 por valores positivos, verificou-se que não foi possível

encontrar um valor limite para o esforço transverso. Como consequência, para r=0 adoptou-se o

procedimento de calcular o valor de esforço transverso muito próximo de r=0, mas com r diferente de

0.

No que diz respeito às cargas concentradas, exceptuando o caso do esforço transverso, é possível

obter os valores limites para os esforços quando r=0. Ou seja, à parte do esforço transverso, dada a

discretização da estrutura e refinamento de malha adoptado, os esforços apresentam um valor limite

para r=0 na aplicação do método dos elementos finitos. No entanto, consoante o refinamento de

malha, em particular o refinamento na vizinhança da carga concentrada, obtém-se esforços que não

convergem na secção, r=0. Os esforços tendem para infinito quanto mais refinada for a malha na

vizinhança da carga concentrada.

59

8 Programação

8.1 Objectivo do programa

Pretende-se conceber um programa de cálculo automático que tem por base duas metodologias:

soluções analíticas e o método dos elementos finitos. O programa tem por objectivo a resolução de

estruturas axissimétricas de casca fina sujeitas a um carregamento axissimétrico. O resultado final do

programa deverá consistir no seguinte:

Obtenção do desenho automatizado da deformada da estrutura e os esforços de membrana

e de flexão;

Obtenção de ficheiros de texto contendo os resultados obtidos.

A informação mais detalhada sobre as possibilidades e uso do programa encontra-se em anexo A.1 -

Manual de utilização.

8.2 Programas utilizados

No que diz respeito à programação foram utilizados dois programas:

Wolfram Mathematica v7

Lazarus v0.9.24

Breves descrições e uso dos programas:

Mathematica é um programa de computador que implementa um sistema de álgebra computacional.

Para além de uma linguagem de programação, contém diversas bibliotecas de programação prontas

a serem usadas para diversos fins e em várias áreas das ciências exactas.

O programa insere-se em diversas áreas da engenharia e matemática, além de servir como um

ambiente para desenvolvimento rápido de programas.

Nesta dissertação, o programa Mathematica foi usado com o objectivo de:

Efectuar a formulação das soluções analíticas e do MEF para cascas axissimétricas;

Análise de limites em pontos de singularidade;

Teste e comprovação dos resultados obtidos;

Exportação para linguagem de programação das matrizes utilizadas nas metodologias

adoptadas.

60

Lazarus é um ambiente de desenvolvimento integrado desenvolvido para o compilador Free Pascal.

O software tem por objectivo ser compatível com o Delphi e, ao mesmo tempo, suportar diversas

arquitecturas de processador e sistemas operacionais.

Free Pascal é um compilador de Object Pascal, e foi desenhado para compilar código com a sintaxe

Delphi ou dialectos Pascal e gerar executáveis para diferentes plataformas a partir de um mesmo

código-fonte.

Neste trabalho, o Lazarus foi utilizado para desenhar e conceber o programa de cálculo automático

de estruturas axissimétricas sujeitas a um carregamento axissimétrico através de duas metodologias:

soluções analíticas e o MEF.

Deve-se ter em atenção que o software Lazarus é uma ferramenta de programação em fase

experimental (beta), pelo que teoricamente o programa desenvolvido neste trabalho poderá

eventualmente apresentar defeitos (bugs) em determinadas situações e condições. O Lazarus

apresenta maiores limitações a nível de componentes gráficas, pelo que é neste campo que poderá

hipoteticamente ocorrer determinados problemas. No entanto, a experiência de utilização do

programa neste trabalho mostrou-se bastante estável e sem problemas significativos a apresentar.

8.3 Grau de precisão

Na programação todas as variáveis utilizadas para os cálculos intermédios e armazenamento de

resultados finais foram definidas como double. O double suporta aproximadamente 15 algarismos de

precisão e está definido no intervalo numérico de 2.23 x 10-308

a 1.79 x 10308

.

As soluções analíticas fornecem os resultados correspondentes à teoria elástico-linear de casca fina

axissimétrica, apresentando apenas erros que incorrem do grau de precisão definido para as

variáveis.

O MEF apenas fornece resultados exactos dos deslocamentos nodais se as expressões utilizadas no

campo de deslocamentos coincidirem com a solução homogénea dos deslocamentos da estrutura.

Na formulação adoptada, o campo de deslocamentos é construído por aproximação através de

funções polinomiais que não coincidem com a equação governativa do campo de deslocamentos das

cascas axissimétricas. Logo, os deslocamentos nodais obtidos serão aproximações, que serão mais

precisas quantos mais elementos forem utilizados numa configuração apropriada.

Embora a regra de quadratura de Gauss de n pontos produza o resultado exacto de uma função

polinomial de grau 2n-1 ou menor, algumas das funções a integrar no método dos elementos finitos

não são polinomiais devido à existência de um operador de divisão polinomial na matriz , o raio.

Consequentemente, a regra de quadratura de Gauss, apenas fornecerá uma aproximação do integral

das funções não polinomiais que será tanto mais preciso quantos mais pontos forem utilizados.

61

A nível de tratamento de dados, efectuou-se arredondamentos para a amostragem de resultados

relevantes na componente gráfica dos diagramas de esforços. Na saída de dados numéricos os

resultados finais vêm em formato de ficheiro “csv”, onde se optou por não alterar os dados de modo a

preservar a precisão. O “csv” (Comma-separated values) é um formato de ficheiro que armazena

dados em forma de tabela, onde os valores ou texto são separados por vírgulas. Por serem bastante

simples, os ficheiros de “csv” são comuns na maioria das plataformas de computador. Existe

bastante escolha de software para abrir ou editar estes ficheiros, como por exemplo: o Microsoft

Excel. Parte-se do princípio que o tratamento de dados, como arredondamentos ou formatos de

apresentação, é feito pelo utilizador através de software compatível com o formato de ficheiro “csv”.

8.4 Base de dados

O sistema de programação funciona com bases de dados, que podem ser distinguidas nas seguintes

partes fundamentais:

Identificação dos elementos, com a identificação dos nós que constituem cada elemento;

Identificação dos nós, com restrições dos graus de liberdade/identificação dos graus de

liberdade;

Diversos dados essenciais e complementares para a formulação das soluções analíticas e do

MEF;

Resultados finais obtidos;

Em geral, as bases de dados são organizadas num sistema de matrizes.

Seguidamente, será apresentada a estrutura das bases de dados.

Base de dados dos elementos:

Identificação do tipo de elemento

Identificação do elemento Identificação do nó inicial Identificação do nó final

No que diz respeito às soluções analíticas, a base de dados dos elementos é construída

directamente com os dados fornecidos ao programa pelo utilizador. No entanto, os elementos de anel

contêm apenas um nó, e na obtenção da base de dados destes elementos, estes identificam apenas

um nó. No que se refere ao MEF, os elementos são obtidos em função do número de elementos

definidos para cada troço ou arco. No caso do troço, os elementos são obtidos em função de uma

divisão igual do troço na sua extensão. No caso do arco, este é dividido igualmente em função do

número de elementos pretendidos, no entanto, os elementos apresentam-se rectos no PSA. Estes

elementos são definidos pelos pontos que dividem o arco em segmentos iguais. A casca esférica é

portanto uma aproximação, composta por cascas cónicas, pelo que é possível concluir que quanto

mais elementos forem utilizados melhor será a aproximação.

62

Base de dados inicial dos nós:

Coordenadas

Restrições dos graus de

liberdade

Identificação do nó r z

Os nós são obtidos em função dos nós definidos pelo utilizador e pelos elementos gerados para cada

troço e arco. As restrições são identificadas quando são colocados números negativos no respectivo

grau de liberdade, por exemplo: “-1”.

Esta base de dados irá ser alterada de forma a identificar os graus de liberdade não restringidos dos

respectivos nós.

Base de dados final dos nós:

Coordenadas

Identificação dos graus de

liberdade não restringidos


Os graus de liberdade restringidos contêm, na parcela respectiva, o valor “0”, os não restringidos são

identificados numericamente por um número natural maior ou igual a 1.

O número total de graus de liberdade não restringidos vai definir as dimensões da matriz de rigidez

global e do vector de forças nodais global. Consequentemente, define também o número de

equações e incógnitas a serem resolvidos. A identificação dos graus de liberdade vai possibilitar o

espalhamento da matriz de rigidez e do vector de forças nodais e a obtenção dos deslocamentos

nodais considerados no problema.

Após a resolução do sistema de equações lineares é gerada uma base de dados com os resultados

finais que têm a seguinte estrutura:

Deslocamentos nodais:

Coordenadas

Deslocamentos no sistema

global de coordenadas


63

Resultados gerais finais:

Identificação do tipo de elemento

Identificação do elemento r z

Na base de dados dos resultados gerais, estes são obtidos em função do número de pontos por

elemento, sendo este parâmetro definido pelo utilizador.

As restantes bases de dados guardam informações como:

Espessura do elemento;

Coeficiente de Poisson;

Módulo de elasticidade;

Carregamento de vão e cargas nos nós;

Coordenada do centro do raio apenas para arcos.

O carregamento de vão foi definido de forma a contemplar as situações mais comuns, como o peso

próprio da estrutura ou um líquido contido pela estrutura. Resultando nos seguintes tipos de

carregamento:

Carregamento trapezoidal ou uniforme com direcção normal ao elemento;

Carregamento uniforme com direcção tangencial ao elemento;

Carregamento uniforme com direcção gravítica.

Os primeiros dois tipos de carregamentos assumem o sentido do eixo local do elemento que

depende do nó inicial e do nó final.

Para a definição da localização dos nós, deve-se ter em consideração a utilização de coordenadas

polares: . No entanto, para definir a coordenada do centro do raio de um arco, r pode tomar

valores negativos.

64

8.5 Fluxograma

Seguidamente, apresenta-se um fluxograma com os passos fundamentais dos procedimentos

efectuados pelo programa:

Preparação preliminar de

dados

Abrir ficheiro de dados

Alterar dados

Guardar ficheiro de dados

Função de

desenho da

estrutura

Função de resolução da estrutura

Construção da matriz

de nós inicial e da

matriz de elementos

Desenho visual da estrutura

Processo de assemblagem dos

elementos

Resolução do sistema de equações

Obtenção do campo de deslocamentos e

tensões

Desenho visual do campo de

deslocamentos e tensões

Saída de ficheiros com dados numéricos dos resultados finais

1º

2º

65

Como é possível observar as funções fundamentais para a resolução de uma estrutura pelas

soluções analíticas e MEF encontram-se no menu Ferramentas do programa, e são as funções:

“Desenhar Estrutura” e “Resolver Estrutura”. As restantes funções são na maioria complementares e

utilizam a base de dados para aceder a informação dos resultados finais.

Os passos que dizem respeito ao processo de assemblagem dos elementos e resolução do sistema

de equações, estão descritos na secção 6.3. Os passos que dizem respeito à obtenção do campo de

deslocamentos e tensões estão descritos na formulação dos diversos elementos das soluções

analíticas e do MEF.

8.6 Cuidados relacionados com a programação/formulação

No que diz respeito às soluções analíticas, em particular o elemento de casca cilíndrica, é

aconselhável respeitar a condição: , de modo a evitar erros numéricos significativos que

são agravados pelo carácter exponencial das funções Y (4.50). Caso seja necessário é sempre

possível cumprir esta condição subdividindo o elemento de casca cilíndrica. O programa verifica esta

condição e alerta o utilizador caso não seja cumprida, apresentando a opção de continuar ou

cancelar a operação.

No que se refere às soluções analíticas, excluindo o elemento de anel, deve-se ter em consideração

que embora a matriz de rigidez e vector de forças nodais estejam definidos no referencial local e seja

aplicada a matriz de transformação de coordenadas para o referencial global, os elementos devem

ser definidos sempre com o mesmo sentido. Na casca cilíndrica o nó inicial deve se situar a uma cota

inferior à do nó final. Na casca em anel o nó inicial deve situar-se a uma coordenada de raio inferior à

do nó final. Embora se efectue uma transformação de coordenadas, os elementos foram definidos

com o sentido mencionado, e devido à formulação adoptada, a construção da matriz de rigidez local

para um sentido contrário ao definido resulta na obtenção de termos errados.

Na definição de um elemento com um nó de raio nulo é necessário restringir os graus de liberdade

da rotação e deslocamento horizontal nesse nó, pois devido à condição axissimétrica estes

deslocamentos encontram-se impedidos, como é explicado no capítulo 7. Este procedimento é

importante quando se utiliza o MEF, no entanto, no caso das soluções analíticas, nomeadamente, na

casca em anel para o caso particular em que o raio inicial é nulo, a formulação dispensa a restrição

dos graus de liberdade por parte do utilizador no nó inicial.

É necessário ter algum cuidado na análise e escrita da programação em Delphi/Pascal pois o

sistema de índices de matrizes e vectores utiliza por defeito “0” como índice inicial e não “1”, no

entanto, algumas rotinas utilizadas, como por exemplo o método de eliminação de Gauss utiliza o

valor 1 como índice inicial, pelo que se deve ter cuidado para não confundir os dois casos.

66

Embora a formulação do MEF utilize o ângulo do elemento com a horizontal, , este ângulo na

realidade nunca é calculado, sendo calculados directamente o e através das

expressões fornecidas na secção 5.2. A utilização de um único ângulo para as duas funções

trigonométricas resulta na obtenção do sinal contrário em determinadas situações.

Houve o cuidado com a saída gráfica dos diagramas de flexão pois é habitual seguir a convenção de

representar o diagrama de momentos do lado das fibras traccionadas.

No desenho dos diagramas de esforços, os valores calculados estão representados por um troço

recto e perpendicular ao elemento, que une o andamento do diagrama com o elemento.

67

9 Exemplos de aplicação e erro

Seguidamente, serão apresentados alguns exemplos de aplicação, onde se pretende efectuar uma

análise simples de estruturas sujeitas a um determinado carregamento. Neste capítulo não se

pretende efectuar o dimensionamento das estruturas ou verificar critérios de segurança estrutural,

embora o programa concebido tenha bastante utilidade nesse sentido. O objectivo deste capítulo

consiste em demonstrar as capacidades do programa, verificar a validade da formulação utilizada, e

efectuar uma análise de erro. No entanto, as estruturas apresentadas serão uma aproximação

simplificada da realidade.

No que diz respeito ao MEF, pretende-se efectuar uma análise ao refinamento h uniforme, ou seja,

um refinamento da malha da estrutura. O refinamento para os troços e arcos é uniforme pois trata-se

de uma função automática do programa.

Como meio para analisar o erro, utilizar-se-á a seguinte norma de erro, expressa em percentagem:

Com:

- Solução analítica

- Solução decorrente do MEF

Esta norma avalia o erro da solução aproximada obtida pelo MEF relativamente à solução analítica, a

qual sempre que exista é tomada como solução de referência.

Para a análise de erro efectuada não se contabiliza os valores em pontos de singularidades onde

não existe valor limite, ou seja, onde o valor tende para infinito.

Todos os gráficos apresentados nos exemplos, que dizem respeito à análise de erro, apresentam

dupla escala logarítmica no erro e no número de elementos utilizados para uma melhor percepção da

evolução do erro.

O sistema de unidades é independente do programa e cabe ao utilizador definir as unidades em que

pretende trabalhar. O programa não indica as unidades dos parâmetros a serem introduzidos nem as

unidades dos resultados. No entanto, deve haver coerência de valores a introduzir de forma a não

misturar sistemas de unidades diferentes, assim, os resultados obtidos serão compatíveis com o

sistema de unidades pretendido pelo utilizador. Deve-se ter em atenção que nos resultados os

esforços são obtidos por unidade de comprimento circunferencial.

68

Seguidamente, apresenta-se um exemplo da utilização de um sistema de unidades que pode ser

tomado como referência nos exemplos:

Dados a introduzir:

Comprimentos: m

Cargas concentradas: kN

Carregamentos com distribuição linear: kN/m

Carregamentos distribuídos em área: kN/m2 (kPa)

Resultados obtidos:

Deslocamentos: m

Esforços de membrana (Nx, Nθ): kN/m

Esforço transverso (Vx): kN/m

Esforços de flexão (Mx,Mθ): kN.m/m

Atendendo à análise efectuada na subsecção 5.3.1, adoptou-se 6 pontos de Gauss em todos os

exemplos que serão apresentados. Deste modo garante-se um pequeno erro de integração ao

mesmo tempo que não se exagera no número de pontos de Gauss, tornando assim o programa mais

eficiente.

Seguem-se os resultados relativos aos vários exemplos adiando-se a sua discussão crítica para a

secção 9.5.

69

9.1 Lajes circulares

Nesta secção, utilizar-se-á a mesma laje circular sujeita a carregamentos e condições de fronteira

diferentes. As características gerais da laje circular são:

Figura 9.1 - Configuração geral de uma laje circular

Dados:

h = 0,5 m (espessura);

L = 5 m (raio da laje);

ν = 0,2 (coeficiente de Poisson);

E = 30e6 kPa (módulo de elasticidade);

Esbelteza da laje:

70

9.1.1 Exemplo nº 1

Laje circular, simplesmente apoiada, sujeita a um carregamento uniforme em área.

Figura 9.2 - Configuração da estrutura e do carregamento do Exemplo nº 1

Para q = 50kN/m2 obteve-se os seguintes deslocamentos e diagramas de esforços:

a)

-0,007

-0,006

-0,005

-0,004

-0,003

-0,002

-0,001

0

0,001

0 1 2 3 4 5 6

w

r

1

2

5

10

20

50

100

200

500

Sol. Ref.

Número de elementos

71

b)

c)

d)

Figura 9.3 - Deslocamentos e diagramas de esforços do Exemplo nº 1 a) Deslocamentos transversais: w; b) Esforço de membrana: Nx; c) Esforço de membrana: Nθ; d) Momento flector: Mx; e) Momento flector: Mθ; f) Esforço transverso: Vx;

0

50

100

150

200

250

300

350

0 1 2 3 4 5 6

Mx

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

50

100

150

200

250

300

350

0 1 2 3 4 5 6

Mθ

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

20

40

60

80

100

120

140

0 1 2 3 4 5 6

Vx

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


72

Através dos resultados foi possível obter o erro relativo à solução de referência que se apresenta de

seguida:

Erro (%)

Nº de elementos

w Mx Mθ Vx

1 0,958 9,953 6,310 52,132

2 5,62E-02 2,472 1,011 22,463

5 1,52E-03 0,524 0,124 8,805

10 1,21E-04 0,238 4,80E-02 7,670

20 7,56E-06 5,90E-02 1,07E-02 3,791

50 1,95E-07 9,40E-03 1,58E-03 1,507

100 7,48E-08 2,35E-03 3,85E-04 0,752

200 1,95E-06 5,88E-04 9,70E-05 0,375

500 2,16E-05 7,25E-05 7,03E-06 0,150

Tabela 9.1 - Erro dos deslocamentos e esforços do Exemplo nº 1

Figura 9.4 - Erro dos deslocamentos e esforços do Exemplo nº 1

1,0E-08

1,0E-07

1,0E-06

1,0E-05

1,0E-04

1,0E-03

1,0E-02

1,0E-01

1,0E+00

1,0E+01

1,0E+02

1 10 100 1000

Erro(%)


Análise de erro

w

Mx

Mθ

Vx

73

9.1.2 Exemplo nº 2

Laje circular, encastrada no bordo, sujeita a uma carga concentrada aplicada no centro.


Para P = 500 kN obteve-se os seguintes deslocamentos e diagramas de esforços:

a)

-0,0009

-0,0008

-0,0007

-0,0006

-0,0005

-0,0004

-0,0003

-0,0002

-0,0001

0

0,0001

0 1 2 3 4 5 6

w

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


74

b)

c)

d)


-100

-50

0

50

100

150

200

250

300

350

400

0 1 2 3 4 5 6

Mx

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


-50

0

50

100

150

200

250

300

350

400

0 1 2 3 4 5 6

Mθ

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

2000

4000

6000

8000

10000

12000

14000

16000

18000

0 1 2 3 4 5 6

Vx

r

1

2

5

10

20

50

100

200

500

Sol. Ref.


75


seguida:

Erro (%)

Nº de elementos

w Mx Mθ Vx

1 7,706 34,783 30,530 62,166

2 1,009 15,781 11,642 45,084

5 6,77E-02 7,043 2,036 23,502

10 1,36E-02 3,702 0,593 24,812

20 1,70E-03 1,684 0,269 20,343

50 1,09E-04 0,628 0,101 16,357

100 1,37E-05 0,305 4,90E-02 14,218

200 1,78E-06 0,150 2,41E-02 12,562

500 3,59E-06 5,94E-02 9,55E-03 10,878



1,0E-06

1,0E-05

1,0E-04

1,0E-03

1,0E-02

1,0E-01

1,0E+00

1,0E+01

1,0E+02

1 10 100 1000

Erro(%)


Análise de erro

w

Mx

Mθ

Vx

76

9.1.3 Exemplo nº 3

Laje circular, encastrada no bordo, sujeita a um carregamento linearmente distribuído ao longo de um

perímetro circular com raio igual a b.


Seguidamente, apresentam-se as soluções analíticas de:

Para p = 500 kN/m e b=2,5 m obteve-se os seguintes deslocamentos e diagramas de esforços:

a)

-0,006

-0,005

-0,004

-0,003

-0,002

-0,001

0

0,001

0 1 2 3 4 5 6

w

r

2

6

10

20

50

100

200

500

Sol. Ref.


77

b)

c)

d)


-600

-500

-400

-300

-200

-100

0

100

200

300

0 1 2 3 4 5 6Mx

r

2

6

10

20

50

100

200

500

Sol. Ref.


-150

-100

-50

0

50

100

150

200

250

300

0 1 2 3 4 5 6

Mθ

r

2

6

10

20

50

100

200

500

Sol. Ref.


-100

0

100

200

300

400

500

600

0 1 2 3 4 5 6

Vx

r

2

6

10

20

50

100

200

500

Sol. Ref.


78


seguida:

Erro (%)

Nº de elementos

w Mx Mθ Vx

2 0,443 5,299 1,941 21,598

6 8,05E-03 1,244 0,347 8,183

10 6,81E-04 0,607 0,167 9,841

20 4,29E-05 0,153 4,21E-02 4,949

50 1,10E-06 2,46E-02 6,74E-03 1,983

100 8,66E-08 6,15E-03 1,69E-03 0,992

200 1,42E-07 1,54E-03 4,22E-04 0,496

500 4,61E-06 2,46E-04 7,11E-05 0,198



1,0E-08

1,0E-07

1,0E-06

1,0E-05

1,0E-04

1,0E-03

1,0E-02

1,0E-01

1,0E+00

1,0E+01

1,0E+02

1 10 100 1000

Erro(%)


Análise de erro

w

Mx

Mθ

Vx

79

9.2 Casca cilíndrica – Exemplo nº 4

Neste exemplo apresenta-se uma casca cilíndrica com rotação livre e deslocamentos impedidos na

base e deslocamentos livres no topo. Esta casca encontra-se sujeita a um carregamento transversal

triangular. A configuração geométrica da estrutura e carregamentos apresenta-se na seguinte

imagem:


Dados:

r = 18 m;

L = 6 m;




Esbelteza da casca:

O carregamento definido por Qyi e Qyf é representativo da pressão hidrostática da água. Considera-se

que a secção final se encontra à cota de superfície livre de água, obtendo-se os seguintes

carregamentos:

Qyi = 0 kN/m2

Qyf =60 kN/m2

80

Com a geometria, propriedades do material e os carregamentos definidos obtém-se os seguintes

deslocamentos e diagramas de esforços:

a)

b)

c)

0

1

2

3

4

5

6

7

0 0,0005 0,001

z

w

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

1

2

3

4

5

6

7

-100 -50 0 50

z

Nx

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

1

2

3

4

5

6

7

-200 0 200 400 600

z

Nθ

1

2

5

10

20

50

100

200

500

Sol. Ref.


81

d)

e)

f)


0

1

2

3

4

5

6

7

-200204060

z

Mx

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

1

2

3

4

5

6

7

-50510

z

Mθ

1

2

5

10

20

50

100

200

500

Sol. Ref.


0

1

2

3

4

5

6

7

-100 -50 0 50

z

Vx

1

2

5

10

20

50

100

200

500

Sol. Ref.

Número de

elementos

82


seguida:

Erro (%)

Nº de elementos

w Nθ Mx Mθ Vx

1 5,445 5,066 35,903 35,903 87,600

2 1,582 2,292 20,694 20,694 80,872

5 0,112 0,559 4,769 4,769 36,447

10 2,13E-02 0,521 2,183 2,183 32,527

20 5,46E-03 0,261 0,540 0,540 16,624

50 8,81E-04 0,104 8,63E-02 8,63E-02 6,684

100 2,21E-04 5,22E-02 2,16E-02 2,16E-02 3,343

200 5,50E-05 2,61E-02 5,39E-03 5,39E-03 1,672

500 1,08E-04 1,04E-02 1,01E-03 1,01E-03 0,669



1,0E-05

1,0E-04

1,0E-03

1,0E-02

1,0E-01

1,0E+00

1,0E+01

1,0E+02

1 10 100 1000

Erro(%)


Análise de erro

w

Nθ

Mx

Mθ

Vx

83

9.3 Casca cónica – Exemplo nº 5

Seguidamente, apresenta-se uma casca cónica simplesmente apoiada sujeita um carregamento

trapezoidal distribuído em área, com direcção perpendicular ao eixo da casca no PSA e um

carregamento uniformemente distribuído em área segundo a direcção gravítica. A configuração

geométrica da estrutura e os carregamentos apresentam-se na seguinte imagem:


Dados:




Esbelteza da casca:

O carregamento definido por Qyi e Qyf é representativo da pressão hidrostática da água, considera-se

que na secção final existe uma coluna de água de 5 metros, obtendo-se os seguintes carregamentos:

Qyi = 80 kN/m2

Qyf =50 kN/m2

O carregamento definido por Qvg é representativo do peso próprio da estrutura e adopta-se o

seguinte valor: Qvg = 15 kN/m2.

84

Com a geometria, propriedades do material e os carregamentos definidos obtém-se os seguintes

deslocamentos e diagramas de esforços:

a)

b)

c)

-2

-1

0

1

2

3

4

-1 0 1 2 3 4 5

z

r

Deformada

1

2

5

10

20

50

100

200

500

Indef.


-50

0

50

100

150

200

250

300

350

0 1 2 3 4 5 6

Nx

x

1

2

5

10

20

50

100

200

500


-800

-600

-400

-200

0

200

400

600

0 1 2 3 4 5 6Nθ

x

1

2

5

10

20

50

100

200

500


85

d)

e)

f)

Figura 9.15 - Deslocamentos e diagramas de esforços do Exemplo nº 5 a) Deformada da estrutura (com deslocamentos afectados por um factor de ampliação de 5000); b) Esforço de membrana: Nx; c) Esforço de membrana: Nθ; d) Momento flector: Mx; e) Momento flector: Mθ; f) Esforço transverso: Vx;

-30

-20

-10

0

10

20

30

40

50

60

70

80

0 1 2 3 4 5 6

Mx

x

1

2

5

10

20

50

100

200

500


-20

-10

0

10

20

30

40

0 1 2 3 4 5 6

Mθ

x

1

2

5

10

20

50

100

200

500


-60

-40

-20

0

20

40

60

80

100

120

140

0 1 2 3 4 5 6

Vx

x

1

2

5

10

20

50

100

200

500


86

9.4 Reservatório – Exemplo nº 6

Pretende-se demonstrar a utilidade do programa concebido num exemplo mais prático e corrente de

estruturas axissimétricas. Escolheu-se, portanto, um reservatório circular constituído por: um anel,

uma laje circular, uma casca cilíndrica, e uma abobada esférica. No que diz respeito a acções,

considera-se o peso próprio do reservatório, o peso da água que este contém e a reacção do solo.

Seguidamente, apresenta-se a configuração da estrutura e carregamentos considerados:


87

Dados:



γb = 25 kN/m3 (peso volúmico do betão);

γw = 10 kN/m3 (peso volúmico da água);

9.4.1 Peso próprio da estrutura

No caso do elemento de anel (elemento nº1) adoptou-se um modelo de cálculo em consola:

Peso próprio:

Figura 9.17 - Modelo de cálculo auxiliar para o elemento de anel

Carregamento nodal:

Elemento de casca circular (elemento nº 2):

Elemento de casca cilíndrica (elemento nº 3):

Arco nº 1:

9.4.2 Acção da água

Elemento de casca circular:

Elemento de casca cilíndrica:

88

9.4.3 Reacção do solo

Considera-se o solo suficientemente deformável para apresentar uma reacção uniformemente

distribuída do somatório das cargas verticais. Pretende-se obter a reacção por unidade de área da

laje circular:

Elemento de anel:

Elemento de casca circular:

Elemento de casca cilíndrica:

Arco nº 1:

Carga vertical a actuar no elemento de casca circular:

9.4.4 Obtenção dos deslocamentos e esforços

Para a obtenção dos deslocamentos e diagramas de esforços utilizou-se 100 elementos finitos de

casca cónica no Arco nº1. A parede lateral cilíndrica e a laje circular de fundo foram modeladas com

recurso às soluções analíticas. Com a geometria, propriedades do material, carregamentos e

refinamento definidos obtém-se os seguintes deslocamentos e diagramas de esforços:

a)

b) deformada

89

c) Nx

d) Nθ

e) Mx

f) Mθ

90

g) Vx

Figura 9.18 - Deslocamentos e diagramas de esforços do Exemplo nº 6 a) Orientação e configuração dos elementos e arco; b) Deformada da estrutura; c) Esforço de membrana: Nx; d) Esforço de membrana: Nθ; e) Momento flector: Mx; f) Momento flector: Mθ ; g) Esforço transverso: Vx;

As imagens da Figura 9.18 são baseadas na saída gráfica de dados do programa desenvolvido. No

entanto, por uma questão de apresentação foram alteradas para permitir uma melhor percepção dos

resultados obtidos. Na Figura 9.18 as escalas dos diagramas de esforços entre os elementos e o

arco são diferentes para uma melhor percepção do andamento dos diagramas de esforços, uma vez

que a amplitude de esforços pode ser muito diferente dificultando a visualização de alguns diagramas

quando é utilizada a mesma escala.

Na Figura 9.18 d) e f) os valores dos esforços a negrito (bold) simbolizam o esforço actuante no

elemento de anel.

91

Verificou-se que alguns diagramas de esforço apresentam oscilações na zona de ligação entre o

Arco nº1 e o Elemento nº3. Estas oscilações não são visíveis na Figura 9.18 por apresentarem

valores muito baixos, pelo que se apresenta de seguida alguns pormenores dos diagramas de

esforços na zona de ligação entre o Arco nº1 e o Elemento nº3, a uma escala mais elevada:

a) Mx

b) Mθ

c) Vx

Figura 9.19 - Diagramas de esforço na zona de ligação entre o Arco nº1 e o Elemento nº3 a) Momento flector: Mx; b) Momento flector: Mθ; c)Esforço Transverso: Vx

92

9.5 Análise dos resultados e do erro

Nesta secção, as referências à solução exacta dizem respeito à solução fornecida pela teoria

elástico-linear de casca fina axissimétrica.

A escolha de vários exemplos de lajes circulares deve-se ao facto de existirem expressões analíticas

simples obtidas da referência bibliográfica [10], pretendendo-se com estas expressões verificar a

validade dos resultados obtidos pelo programa tanto das soluções analíticas como do MEF.

Verificou-se a validade dos resultados com todas expressões analíticas utilizadas. Além disso, a

validade dos resultados obtidos pelo programa é apoiada pela convergência das soluções do MEF

nos resultados das soluções analíticas, verificando que os dois métodos diferentes atingem o mesmo

objectivo.

É possível verificar que a convergência de resultados é variável e pode depender não só da

estrutura, mas também das condições de fronteira, do tipo de carregamento e do tipo de resultado a

obter, nomeadamente, deslocamentos e esforços. No entanto, em geral, considera-se que os

resultados convergem relativamente depressa.

Observando os gráficos dos resultados, verifica-se que os deslocamentos convergem mais

rapidamente para a solução de referência, seguidos dos esforços de membrana e flexão obtidos por

deformação dos elementos e finalmente o esforço transverso obtido por equilíbrio. Esta situação

também se observa nos gráficos de análise de erro, onde o mesmo decresce, em geral, de forma

linear na dupla escala logarítmica e que, em geral, os deslocamentos decrescem a uma taxa superior

seguidos dos esforços de membrana e de flexão e por ultimo o esforço transverso. Isto acontece

porque o método dos elementos finitos baseia-se na aproximação de um campo de deslocamentos,

obtendo-se os esforços de membrana e flexão a partir do campo de deformações obtidos a partir do

campo de deslocamentos, e o esforço transverso é obtido, por sua vez, através de equilíbrio a partir

dos esforços de flexão. Quanto mais indirecto for o grau de obtenção de resultados pior será a

aproximação.

É possível observar nos gráficos de análise de erro que os deslocamentos tendem a apresentar um

erro crescente, invertendo a tendência, a partir de 100 ou 200 elementos, deduz-se que este

acontecimento se deve a erros de precisão. À medida que o erro se torna inferior, algo que acontece

rapidamente para o caso dos deslocamentos, a norma de erro utilizada fica cada vez mais

susceptível a erros de precisão, podendo gerar erros superiores mesmo utilizando mais elementos

finitos.

Verifica-se na Figura 9.13, correspondente ao gráfico de análise de erro do exemplo nº4, que o erro

dos esforços de flexão longitudinal e tangencial são coincidentes, isto deve-se ao facto de que no

elemento de casca cilíndrica a relação entre os dois esforços é a seguinte: . Ou seja, os

diagramas de momento flector são apenas afectados por uma constante, o coeficiente de Poisson,

resultando na obtenção do mesmo erro.

93

A nível geral, observando os gráficos dos resultados, e ignorando os casos de descontinuidade

significativa entre elementos nos esforços, verifica-se que a partir dos 5 elementos finitos já existe

uma boa convergência de resultados. Com base nos resultados obtidos, para esbeltezas da ordem

dos exemplos utilizados, verifica-se que estabelecendo como ordem de grandeza 20 elementos

finitos, obtém-se uma boa aproximação dos deslocamentos e esforços da solução de referência.

Verificou-se que as duas metodologias, as soluções analíticas e o MEF, são compatíveis. É possível

observar que os resultados obtidos no exemplo nº6 apresentam bastante conformidade. Também se

efectuou uma verificação não documentada com o exemplo nº1 e nº4 onde se subdividiu a respectiva

estrutura a metade no PSA, utilizando elementos das duas metodologias, constatou-se que os

resultados obtidos estão de acordo com as soluções de referência.

9.5.1 Descontinuidades no campo de esforços

Ao longo de todos os exemplos é possível observar que o campo de deslocamentos é sempre

contínuo, no entanto, o mesmo não se pode afirmar para os diagramas de esforços obtidos,

observando-se uma descontinuidade entre elementos que vai diminuindo à medida que se efectua

um maior refinamento. O campo de deslocamentos é contínuo, porque a formulação do MEF utilizada

baseia-se na aproximação de um campo de deslocamentos contínuo. O mesmo já não acontece para

o campo de tensões generalizado, porque é função de derivadas de ordem superior para as quais

não é estabelecida continuidade a priori. A satisfação do equilíbrio baseia-se no PTV que estabelece

que o trabalho interno de deformação é igual ao trabalho externo das forças aplicadas. Como o

campo de deslocamentos não é exacto, os diagramas de esforços obtidos apresentam

descontinuidades as quais são mais aparentes quando se usa um baixo número de elementos finitos.

À medida que o campo de deslocamentos se vai aproximando da solução exacta o campo de

tensões generalizado apresentará menos descontinuidades e ficará mais próximo também da

solução exacta.

As descontinuidades entre elementos nos esforços não devem ser confundidas com algumas

situações onde ocorrem erros provenientes da formulação utilizada, como é o caso do esforço de

membrana do exemplo nº 4. O exemplo nº 4 diz respeito a uma casca cilíndrica sujeita a uma

carga triangular transversal, não apresenta nenhum carregamento longitudinal e encontra-se livre no

topo, mas com deslocamentos impedidos e rotação livre na base. Nestas condições, é possível

deduzir que o esforço membrana é nulo, no entanto, esse não foi o resultado obtido como é

possível observar na Figura 9.12. Analisando a expressão do esforço de membrana (3.17) é

possível concluir que o erro provém da contribuição da extensão tangencial afectada pelo coeficiente

de Poisson e que se deve às funções de forma adoptadas na presente formulação do MEF. Nesta

formulação verifica-se que a extensão longitudinal pode ser constante ou nula, enquanto que a

extensão tangencial pode ter função cúbica. Acontece que, no exemplo demonstrado, a extensão

tangencial não era constante, mas de grau superior, logo, a função da extensão longitudinal não

permite a anulação da componente tangencial, resultando nas descontinuidades observadas no

94

diagrama de esforço de membrana . Conhecendo a origem do erro é possível efectuar algumas

demonstrações, que se apresentam de seguida:

1a)

1b)

1c)

2a)

2b)

3a)

3b)

Figura 9.20 – Análise de erros provenientes da formulação 1: Exemplo nº 4 – a) Deformada da estrutura; b) Esforço axial Nx; c) Esforço axial Nx com coeficiente Poisson nulo 2: Exemplo nº 4 com deslocamento horizontal livre na base – a) Deformada da estrutura; b) Esforço axial Nx 3: Exemplo nº 4 com deslocamento horizontal livre na base e carregamento transversal constante e não triangular – a) Deformada da estrutura; b) Esforço axial Nx

95

Para a obtenção dos resultados apresentados na Figura 9.20 utilizou-se um refinamento de malha

uniforme com 20 elementos finitos.

Na Figura 9.20, não foi apresentado o diagrama do esforço de membrana tangencial pois o

andamento é semelhante à deformada da estrutura. Verifica-se que na Figura 9.20 1c) o problema

fica resolvido ao eliminar a componente de extensão tangencial utilizando um coeficiente de Poisson

nulo. Na Figura 9.20 2a) verifica-se que uma vez que o andamento da deformada é linear o erro

obtido no esforço membrana é mais uniforme conforme se pode observar na Figura 9.20 2b). Na

Figura 9.20 3a) verifica-se que o andamento da deformada é constante pelo que o erro devido à

formulação não deveria existir, uma vez que a componente de extensão longitudinal deveria anular a

componente de extensão tangencial. No entanto, o que parece ser descontinuidades do esforço de

membrana , por defeito da formulação, apresentadas na Figura 9.20 3b) trata-se na realidade de

erros numéricos, pois os valores que o diagrama apresenta são de ordem muito reduzida, tratando-

se de um indicativo de que os valores deveriam ser nulos. Uma possível solução deste problema

será utilizar funções de forma que permitam uma função de extensão longitudinal com um grau

polinomial suficientemente elevado para anular a componente de extensão tangencial.

No que diz respeito ao MEF, no exemplo nº2 é possível verificar que a laje circular sujeita a uma

carga concentrada no eixo de revolução, resulta na obtenção de esforços que tendem para infinito na

secção de aplicação da carga à medida que se efectua um maior refinamento que abrange essa

mesma secção. Para cada malha de refinamento, a solução é determinada e tem limite, no entanto,

verifica-se que com um refinamento de malha sucessivamente maior a solução dos esforços na

secção da carga concentrada não converge, pelo que tende para infinito. Esta situação deve-se a

uma singularidade de esforços na secção onde está aplicada a carga concentrada, situação que se

encontra explicada no capítulo 7.

9.5.2 Refinamento adaptativo

Como foi referido o programa concebido permite um refinamento automático e uniforme da malha no

PSA. No entanto, é essencial conhecer o conceito de refinamento adaptativo. O refinamento

adaptativo é mais eficaz em determinadas situações, nomeadamente, em secções onde ocorrem

perturbações relativas à uniformidade da estrutura e carregamento. Geralmente estas situações

ocorrem nas ligações entre elementos onde as características se alteram, zonas onde varia:

espessura, material ou configuração geométrica, nas zonas sujeitas a carregamentos localizados ou

alteração da configuração do carregamento e nas zonas sujeitas a determinadas restrições. Nestas

situações é aconselhável a utilização de um refinamento adaptativo, pelo que as zonas que sofrem

as perturbações devem apresentar maior refinamento que as restantes para uma convergência mais

rápida e eficaz da solução.

96

Seguidamente, apresenta-se uma solução de refinamento adaptativo aplicada ao exemplo nº 3, com

o propósito de demonstrar o seu potencial. Neste exemplo é possível verificar que os esforços

obtidos são constantes para secções com raio inferior ao raio de aplicação da carga, constituindo um

bom exemplo para aplicação de um refinamento adaptativo que tira partido desta particularidade

tornando mais eficiente a obtenção de resultados. Tendo em conta este aspecto, optou-se por utilizar

a seguinte malha de refinamento adaptativo, utilizando 6 elementos finitos:

Figura 9.21 - Aplicação de um refinamento adaptativo ao Exemplo nº 3

É possível observar na Figura 9.21 que apenas se utiliza um elemento finito antes do ponto de

aplicação da carga, pois os diagramas de esforços são constantes até esse ponto, o elemento

seguinte (2) é o mais pequeno para captar melhor as perturbações da aplicação da carga na

vizinhança e os elementos seguintes (3 a 6) são uniformes.

Efectuando uma análise de erro com uma malha uniforme de 6 elementos finitos e a malha

apresentada na Figura 9.21 obteve-se os seguintes erros utilizando a solução de referência:

Erro (%)

w Mx Mθ Vx

Ref. uniforme 8,05E-03 1,244 0,347 8,183

Ref. adaptativo 8,66E-04 0,250 0,069 4,734

9,295 4,977 5,059 1,729

Tabela 9.5 - Erro do refinamento uniforme e adaptativo aplicado ao Exemplo nº 3

Observando a Tabela 9.5 verifica-se que a utilização do refinamento adaptativo escolhido resulta

num menor erro comparativamente ao refinamento uniforme. É possível observar que, relativamente

ao refinamento uniforme, no refinamento adaptativo o erro é cerca de: dez vezes inferior no caso dos

deslocamentos, cinco vezes inferior nos esforços dos momentos flectores e duas vezes inferior no

esforço transverso. De destacar que existem várias soluções de refinamento adaptativo possíveis e

que a utilizada não é, necessariamente, a que gera menor erro, podendo existir configurações de

refinamento mais eficazes. Deve-se ter em atenção que o exemplo utilizado é bastante simples e

apresenta bastante uniformidade, apenas se procurou demonstrar o potencial do refinamento

adaptativo. O refinamento adaptativo torna-se mais útil e apresenta melhorias significativas dos

resultados para configurações mais complexas da estrutura e carregamento.

97

10 Conclusão

Em problemas de estruturas axissimétricas sujeitas a um carregamento axissimétrico é essencial

conhecer a condição axissimétrica que permite a adopção de determinadas simplificações. Todos os

problemas podem ser representados no plano de simplificação axissimétrica (PSA) e o problema

torna-se muito mais simples, semelhante à resolução de um problema de barras num plano

bidimensional.

Como elemento de base utilizou-se a casca cónica para a resolução de estruturas axissimétricas

devido à sua generalidade e simplicidade. O elemento é definido com um ângulo geral, α, relativo à

horizontal. Para α=0º ou α=180º o elemento comporta-se como uma casca circular, que engloba

também as lajes e placas circulares. Para α=90º ou α=270º o elemento comporta-se como uma

casca cilíndrica. Se α for diferente dos ângulos referidos o elemento comporta-se como uma casca

cónica.

Verificou-se que as duas metodologias, soluções analíticas e o MEF, são compatíveis e

complementam-se, enriquecendo a funcionalidade do programa. As soluções analíticas têm algumas

vantagens em relação ao MEF, na medida em que necessitam de apenas um elemento que fornece

os resultados pretendidos, tornando mais eficiente a utilização de recursos do programa

comparativamente à obtenção de uma solução aceitável proveniente do MEF. No entanto, o MEF,

também tem algumas vantagens em relação às soluções analíticas, pois permite uma configuração

estrutural mais diversificada e possibilita sempre a utilização de um carregamento transversal

trapezoidal e um carregamento longitudinal constante.

O tipo de elementos estruturais adoptado no programa foi seleccionado com base nas estruturas

axissimétricas mais gerais. Também se adoptou os carregamentos para os casos mais gerais de

acções em estruturas axissimétricas, nomeadamente, carregamento trapezoidal transversal para

simular, por exemplo, a pressão hidrostática, carregamento uniforme gravítico para simular, por

exemplo, o peso próprio da estrutura e carregamento uniforme longitudinal. A ter em atenção que o

elemento de casca em anel permite apenas um carregamento de vão uniforme e transversal, por ser

o caso mais comum de carregamento de cascas em anel ou circulares.

No que diz respeito ao MEF, no PSA, os elementos de casca cónica são representados como rectas

unidas por um ponto inicial e um ponto final. A representação por elementos rectos é muito

importante pois qualquer estrutura axissimétrica pode ser poligonizada no PSA, sendo este um

processo simples de aproximação à configuração da estrutura real. Quanto maior for o refinamento,

mais próximo será o comportamento relativo à estrutura real. Tendo isto em mente, criou-se a

possibilidade de utilizar arcos no programa com representação no PSA, sendo esses arcos

poligonizados por elementos rectos. A utilização de arcos possibilita a utilização de cascas esféricas

ou abóbadas.

98

Nas metodologias adoptadas para a resolução de estruturas axissimétricas sujeitas a um

carregamento axissimétrico, não só é importante conhecer as possibilidades da formulação como

também é importante conhecer os defeitos. Neste sentido, é importante conhecer as possíveis

singularidades, descritas no capítulo 7, e efectuar uma análise de resultados e erro, descritos na

secção 9.5.

Na secção 9.5, Análise dos resultados e do erro, verificou-se que os resultados obtidos pelo método

dos elementos finitos convergem relativamente depressa e, nas situações verificadas, apresentam

resultados muito próximos das soluções de referência. Verificou-se também que a nível de resultados

obtidos, os deslocamentos convergem mais rapidamente, seguidos dos esforços de membrana e

flexão e por ultimo o esforço transverso.

Verifica-se que o programa concebido constitui uma óptima ferramenta para a resolução de

estruturas axissimétricas sujeitas a um carregamento axissimétrico, podendo ser muito útil para o

auxílio do dimensionamento deste tipo de estruturas.

No que diz respeito à programação, houve uma “disputa” entre repetição de código e consumo de

memória. Em geral, optou-se por consumir menos memória possível, guardando-se apenas variáveis

essenciais, como as bases de dados fundamentais à resolução do problema e à

obtenção/exportação de resultados. Considerou-se que para um refinamento de malha muito

elevado, os requisitos de consumo de memória são mais importantes que a repetição de

procedimentos, uma vez que os processadores da actualidade têm ciclos de processamento por

segundo muito elevados.

O programa desenvolvido para este trabalho é especialmente concebido para a resolução de

estruturas axissimétricas sujeitas a um carregamento axissimétrico. Neste âmbito, o programa tem

vantagens de simplificar a definição da estrutura e carregamento, e obtenção de resultados. A

concepção do programa teve por objectivo maximizar a simplificação axissimétrica com base no eixo

e secção de revolução. Após a definição da configuração da estrutura por elementos, troços e arcos,

do carregamento e definição das propriedades dos materiais, o programa permite um refinamento

automatizado e uniforme ao longo dos troços e arcos definidos. O processo de refinamento fica então

bastante simplificado e torna-se bastante prático, sendo fácil mudar o refinamento de malha

conforme o desejo do utilizador. No entanto, como se verificou na subsecção 9.5.2, o refinamento

uniforme em determinados casos é menos eficaz que um refinamento adaptativo, sendo necessário

uma configuração de troços e arcos específica para atingir este objectivo.

Em alguns diagramas de esforços observaram-se determinadas oscilações, sendo o exemplo nº 6 o

caso mais representativo. Neste exemplo verificou-se que as oscilações ocorriam no diagrama de

momentos flectores e de esforço transverso. Também se observou que as oscilações diminuem

substancialmente ao longo do elemento, sendo imperceptíveis a determinadas escalas compatíveis

com os valores máximos. Este fenómeno ocorre devido a uma propriedade das cascas axissimétricas

relacionada com a rigidez radial. Uma casca axissimétrica pode ser equiparada a uma viga assente

99

em fundação elástica apresentando rigidez ao longo do domínio segundo a direcção radial.

Analisando a casca cilíndrica na Figura 9.18 e), f) e g) e na Figura 9.19, verifica-se que o

carregamento triangular vai originar maiores perturbações na base. Equiparando a casca cilíndrica a

uma viga assente em fundação elástica, verifica-se que o carregamento vai impor curvaturas na

estrutura que são contrariadas por contra-curvaturas devido à fundação elástica, que vão diminuindo

à medida que diminui o carregamento, originando as oscilações.

Principais aspectos a melhorar no programa:

Um dos aspectos a melhorar no programa passa pela importação de um ficheiro com a configuração

da estrutura, permitindo a importação de configurações mais complexas. O ficheiro seria criado pelo

utilizador através de um programa de folha de cálculo como o Excel ou de desenho como o

AutoCAD.

Outro aspecto importante será a possibilidade de efectuar uma ampliação de uma zona da estrutura.

O programa permite aumentar a escala da deformada e de esforços e maximizar a janela obtendo-se

maior espaço de desenho, mas para estruturas complexas será útil a possibilidade de ampliação

específica e controlada de uma parte da estrutura.

A interface gráfica relativa à introdução de dados foi desenhada de forma a facilitar a implementação

de dados no código. Poderá ser útil o investimento numa interface gráfica mais acessível ao

utilizador, como a utilizada em programas comerciais, sendo um caso de exemplo: o SAP.

Para um número elevado de graus de liberdade, a matriz de rigidez global gerada pelo programa é

muito esparsa contendo muitas parcelas com valor nulo, além de ser uma matriz simétrica. A

obtenção da matriz de rigidez e método de resolução do sistema de equações podem ser

melhorados tendo em conta as características da matriz de rigidez global. Com a implementação de

rotinas especializadas, o código torna-se mais eficiente, diminuindo o tempo de cálculo e o consumo

de memória.

O programa não foi concebido para detectar todos os erros decorrentes do utilizador, pelo que o

utilizador deve ter cuidado na introdução de dados. A introdução de dados errados, em algumas

situações, pode causar uma mensagem de erro, a terminação do programa ou o sobre carregamento

do programa causando a sua paralisação e ausência de resposta perante o utilizador. Neste sentido,

poderá ser útil a implementação de rotinas de verificação da integridade de dados e a construção de

um sistema de interface menos propício a erros decorrentes do utilizador.

100

11 Bibliografia

Referências bibliográficas:

[1] Advanced Finite Element Methods (ASEN 6367), Aerospace Engineering Sciences,

University of Colorado at Boulder, Spring 2009

[2] Azevedo, A. F. M., "Método dos Elementos Finitos", Faculdade de Engenharia da

Universidade do Porto - Portugal, Abril de 2003

[3] Corrêa, Manuel C. R., “Cálculo de cascas axissimétricas pelo método dos deslocamentos –

Implementação de um programa de cálculo automático usando soluções analíticas”,

Mestrado em Engenharia de Estruturas, Disciplina de Estruturas Laminares, Dezembro 1993

[4] Flügge, W., "Stresses in shells", Berlin : Springer, 1973

[5] Ghali, A., “Circular storage tanks and silos”, Wiley & Sons, 1979

[6] Kelkar, Vasant S.; Sewell, Robert T., "Fundamentals of the analysis and design of shell

structures", Englewood Cliffs: Prentice-Hall, 1987

[7] Martins, J. A. C., “Teoria elástica linear de cascas finas”, Abril 1989

[8] Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T., “Numerical Recipes in

FORTRAN: The Art of Scientific Computing”, Cambridge University Press, second edition,

1992, Cambridge, England.

[9] Ross, C. T. F., "Finite element programs for axisymmetric problems in engineering",

Chichester : Ellis Horwood, 1984

[10] Timoshenko, S.; Woinowsky-Krieger, S., “Theory of plates and shells”, Auckland : McGraw-

Hill, 1981

Referências de páginas na internet (programação):

[11] http://delphi.about.com/

[12] http://www.delphibasics.co.uk/

[13] http://www.efg2.com/Lab/Library/Delphi/MathFunctions/index.html

http://delphi.about.com/

http://www.delphibasics.co.uk/

http://www.efg2.com/Lab/Library/Delphi/MathFunctions/index.html

A - 1

A Anexos

A - 2

A - 3

A.1 Manual de utilização

A.1.1 Breve introdução

Figura A1 - Programa de cálculo automático “AxisSim”

O programa denomina-se “AxisSim” e trata-se de um programa com interface gráfica que resolve

estruturas axissimétricas de casca fina sujeitas a um carregamento axissimétrico. É um programa de

cálculo numérico automático que tem por base duas metodologias: soluções analíticas e o método

dos elementos finitos. O programa permite um “output” de resultados gráfico e numérico.

Define-se plano de simplificação axissimétrica (PSA) como o plano que contém o eixo e a secção

transversal de revolução que geram uma estrutura axissimétrica.

O programa foi concebido de forma a utilizar a simplificação axissimétrica, pelo que a definição do

problema feita pelo utilizador deve ser efectuada no PSA, obtendo-se resultados representados

nesse mesmo plano.

A - 4

A.1.2 Menu

A.1.2.1 Ficheiro

Figura A2 - Menu “Ficheiro”

Nova Estrutura: Prepara a interface gráfica para uma nova estrutura. Elimina todos os elementos,

troços e arcos, mas cria 2 nós sem definições, pois trata-se do número mínimo de nós na definição

de uma estrutura.

Abrir Ficheiro de Dados: Abre um ficheiro de dados ".dat" que preenche a interface gráfica dos

seguintes separadores: “Nós”, ”Anéis”, “Cascas”, “MEF” e dentro de “Opções” as frames “Solução

Analítica” e “Método dos Elementos Finitos”.

Guardar Ficheiro de Dados: Guarda as definições dos separadores: “Nós”, ”Anéis”, “Cascas”,

“MEF” e dentro de “Opções” as frames “Solução Analítica” e “Método dos Elementos Finitos”, num

ficheiro de dados ".dat".

Exportar Imagem em Display: Exporta a imagem em display no separador "Estrutura", em formato

".png".

Exportar Resultados: Todos os resultados são exportados em formato de ficheiro ".csv".

Nós e Deslocamentos Nodais: Exporta a base de dados dos nós com as coordenadas e os

respectivos deslocamentos nodais no sistema de coordenadas globais.

Elementos: Exporta a base de dados dos elementos utilizados, identificando o tipo de

elemento e os respectivos: nó inicial e nó final, ou apenas um nó no caso do anel.

Resultados Gerais: Exporta os resultados finais em função do número de pontos por

elemento definido. O ficheiro contém as coordenadas e os deslocamentos nos pontos em

sistema de coordenadas global e os esforços nos respectivos pontos.

Sair: Termina o programa.

A - 5

A.1.2.2 Ver

Figura A3 - Menu “Ver”

Nós: Selecciona o separador "Nós".

Anéis: Selecciona o separador "Anéis".

Cascas: Selecciona o separador "Cascas".

MEF: Selecciona o separador "MEF".

Estrutura Indeformada: Selecciona o separador "Estrutura" e mostra a estrutura indeformada na

imagem em display.

Estrutura Deformada: Selecciona o separador "Estrutura" e mostra a estrutura deformada na

imagem em display.

Esforços na Estrutura: Selecciona o separador "Estrutura" e mostra na imagem em display os

seguintes diagramas de esforços na estrutura:

Membrana: Nx

Membrana: Nt ( )

Flector: Mx

Flector: Mt ( )

Transverso: Vx

A - 6

A.1.2.3 Ferramentas

Figura A4 - Menu "Ferramentas”

Opções: Selecciona o separador "Opções".

Desenhar Estrutura: Selecciona o separador "Estrutura" e efectua o desenho da estrutura

apresentando-o na imagem em display.

Resolver Estrutura: Selecciona o separador "Estrutura" e efectua a resolução da estrutura

apresentando o desenho da estrutura indeformada na imagem em display.

A.1.2.4 Ajuda

Figura A5 - Menu "Ajuda"

Manual de Utilização: Apresenta o manual de utilização que descreve as possibilidades e uso do

programa.

Sobre: Apresenta uma janela com uma breve descrição do programa e os respectivos créditos.

A selecção do item "Desenhar Estrutura" activa os seguintes itens:

Estrutura Indeformada

Exportar Imagem em Display

A selecção do item "Resolver Estrutura" activa os mesmos itens e ainda:

Estrutura Deformada

Esforços na Estrutura

A - 7

A.1.3 Separadores

A.1.3.1 Nós

Figura A6 - Separador "Nós"

Neste separador é possível definir o número de nós tendo em conta o número de elementos e

configuração a utilizar.

A 1ª tabela diz respeito às coordenadas de cada nó (r,z) com .

A 2ª tabela diz respeito a dois parâmetros distintos:

Cargas nodais;

Restrição dos graus de liberdade de um determinado nó.

Ambos os parâmetros dizem respeito ao sistema global de coordenadas:

Figura A7 - Eixos do sistema global de coordenadas

Para definir a restrição de um determinado grau de liberdade de um nó é necessário colocar um

número negativo no respectivo grau de liberdade, por exemplo: “-1”.

A - 8

A.1.3.2 Anéis

Figura A8 - Separador "Anéis"

Este separador diz respeito ao elemento de anel das soluções analíticas.

Permite definir o número de elementos de anel a utilizar e atribuir as seguintes características:

Nó (onde está definido o anel)

Área (do anel)

Inércia (no anel)

E – módulo de elasticidade

A - 9

A.1.3.3 Cascas

Figura A9 - Separador "Cascas"

Este separador diz respeito aos seguintes elementos definidos nas soluções analíticas:

Cascas em anel e o caso particular de cascas circulares quando o nó inicial (Nó i) tem raio

nulo;

Cascas cilíndricas.

Neste separador é possível definir o número de elementos a utilizar e os seguintes parâmetros:

Nó i – Nó inicial;

Nó f – Nó final;

h – espessura da casca;

niu ( ) – Coeficiente de Poisson;

E – Módulo de Elasticidade;

Qy i – Valor inicial do carregamento linear de vão com direcção transversal ao troço;

Qy f – Valor final do carregamento linear de vão com direcção transversal ao troço;

Qx – Valor do carregamento uniforme de vão com direcção longitudinal ao troço;

Qvg – Carregamento vertical uniforme de vão com direcção gravítica;

A - 10

A.1.3.4 MEF

Figura A10 - Separador "MEF"

Este separador diz respeito ao método dos elementos finitos é composto por duas tabelas, a primeira

tabela permite representar elementos que podem ser definidos por troços rectos no PSA, é portanto,

possível representar os seguintes elementos axissimétricos:

Cascas cónicas;

Cascas cilíndricas;

Cascas/lajes/placas circulares;

Cascas/lajes/placas em anel.

Esta tabela apresenta as seguintes variáveis a definir para cada troço:

Nó i – Nó inicial;

Nó f – Nó final;

h – espessura da casca;

niu ( ) – Coeficiente de Poisson;

E – Módulo de Elasticidade;

Qy i – Valor inicial do carregamento linear de vão com direcção transversal ao troço;

Qy f – Valor final do carregamento linear de vão com direcção transversal ao troço;

Qx – Valor do carregamento uniforme de vão com direcção longitudinal ao troço;

A - 11

Qvg – Carregamento vertical uniforme de vão com direcção gravítica;

Nº Elem – Número de elementos finitos a utilizar num determinado troço.

Os carregamentos de vão, transversal e longitudinal, têm os sentidos dos eixos locais do troço que

são definidos pelo nó inicial e final. O carregamento gravítico vertical é independente dos eixos

locais, apresentando sempre a direcção do eixo z e sentido contrário.

A segunda tabela diz respeito a elementos que podem ser representados por arcos no

PSA, sendo portanto, possível representar:

Cascas esféricas;

Abobadas.

Os arcos são obtidos por aproximação de troços rectos que resultam da união de pontos no seu

domínio definidos pelo número de elementos utilizados, pelo que quantos mais elementos forem

utilizados melhor será a aproximação a um arco.

Para além das variáveis da tabela anterior, é também necessário definir para esta tabela:

rC – Coordenada r do centro do arco definida no sistema cartesiano;

rZ – Coordenada z do centro do arco definida no sistema cartesiano.

A - 12

A.1.3.5 Opções

Figura A11 - Separador "Opções"

Definições:

– Selecciona as definições padrão presentes na Figura A11

– Guarda as definições

Exibir

Nós: Assinala todos os nós considerados com círculos azuis representados no desenho da

estrutura.

Identificação dos Nós: Identifica numericamente todos os nós.

Identificação dos Elementos: Identifica numericamente todos os elementos.

Valores Relevantes: No desenho dos diagramas de esforços são exibidos os valores mais

relevantes: valor inicial e final, mínimo e máximo de cada elemento/troço/arco.

Orientação dos Elementos: São exibidos os eixos locais de cada elemento das soluções

analíticas, excepto o elemento de anel, e para cada troço e arco do método dos elementos

finitos onde é possível deduzir a orientação dos elementos que os constituem.

Origem: Permite exibir sempre a origem do referencial correspondente ao ponto (r,z)=(0,0),

este ponto é assinalado com um círculo azul e a letra “O”. Esta função pode ter utilidade em

A - 13

determinados casos em que os elementos definidos estão muito afastados do eixo de

revolução e este deixa de ser visível.

Representação de Elementos Anel: Representa os elementos de anel com um círculo a cor

de laranja no desenho da estrutura.

Visualização

Manter proporções: As proporções geométricas da estrutura são mantidas no desenho.

Esticar desenho: O desenho da estrutura será esticado (distorcido) para maximizar o

espaço disponível.

Escala da Deformada/Escala dos esforços

Automática: A escala é determinada automaticamente pelo programa tendo em conta a

amplitude máxima de valores.

Factor de Escala: O utilizador pode escolher o factor de escala a ser utilizado.

Aconselha-se o utilizador a escolher inicialmente a escala “Automática”. Consoante a visualização da

deformada estrutura ou dos esforços, o factor de escala automático é indicado na caixa de texto

respectiva. Com o valor indicado é possível ter uma noção da grandeza da escala possibilitando uma

alteração de escala mais acertada.

Solução Analítica

Nº de pontos por elemento: Esta opção permite definir o número de pontos a obter do

campo de deslocamentos e esforços dos elementos que utilizam as soluções analíticas,

excepto o elemento de anel. O número mínimo de pontos a definir é 2, tratando-se do ponto

inicial e final.

Método dos Elementos Finitos

Nº de pontos por elemento: Esta opção permite definir o número de pontos para obtenção

do campo de deslocamentos e esforços num determinado elemento finito. O número mínimo

de pontos a definir é 2, tratando-se do ponto inicial e final.

Integração numérica, Nº de pontos de Gauss: Permite a definição do número de pontos de

Gauss a utilizar na integração numérica da matriz de rigidez elementar.

A - 14

A.1.3.6 Estrutura

Figura A12 - Separador "Estrutura"

Permite a visualização da estrutura, a deformada e os diagramas de esforços.

Figura A13 - Caixa de selecção "Exibir"

A - 15

Exibir: Permite seleccionar a visualização de:

Estrutura Indeformada

Estrutura deformada

Esforços: Nx

Esforços: Nt ( )

Esforços: Mx

Esforços: Mt ( )

Esforços: Vx

A selecção do item "Desenhar Estrutura" activa apenas o primeiro item da caixa de selecção “Exibir”,

a selecção do item "Resolver Estrutura" activa todos os itens.

– Actualiza a imagem

Esta função é particularmente útil caso se efectue alguma alteração de opções visuais, no separador

“Opções”.

Caso se pretenda alterar as seguintes opções:

Nº de pontos por elemento;

Integração numérica, Nº de pontos de Gauss.

Qualquer opção dos separadores: “Nós”, “Anéis”,“Cascas” e “MEF”

É necessária a selecção do item "Resolver Estrutura" para que tomem efeito.

Factor de escala (Deformada/Esforços): Permite a fácil e prática alteração do factor de escala da

visualização seleccionada, excepto quando se visualiza a estrutura indeformada.

– Incrementa o factor de escala

– Diminui o factor de escala

– Selecciona a opção de escala automática

A alteração do factor de escala é uma função de potência com base: 1.25 no incremento e 0.8 na

diminuição, o que resulta num aumento ou diminuição progressivo da escala em 25%.

A - 16

A - 17

A.2 Código fonte

A.2.1 Unidade “Unit1” (principal)

Algum texto da unidade foi removido devido à ocupação excessiva de espaço, e por conter apenas informação

trivial não necessária à compreensão do código fonte.

unit Unit1;

{$mode objfpc}{$H+}

interface

uses

Classes, SysUtils, LResources, Forms, Controls, Windows, Graphics, Dialogs,

ExtCtrls, StdCtrls, Menus, ComCtrls, Grids, Math, IniFiles, Equacoes, CAxisSim,

MEFAxisSim, NPGaussL, OpMat, Unit2;

type

{ TForm1 }

TMatriz = array of array of double;

TVector = array of double;

MatInt = array of array of integer;

TForm1 = class(TForm)

private

{ private declarations }

Nel: array[1..9] of TEdit;

Mdados: array[1..7] of TStringGrid;

xi, yi, dx, dy, zoomres: double;

ires: array[0..6] of integer;

resf: TMatriz;

public

{ public declarations }

end;

var

Form1: TForm1;

L, ca, sa, r1, r2, Qy1, Qy2, Qx, h, E, v, B: double;

nos, diag, Data: TMatriz;

El, ires: array of array of integer;

idat: integer;

implementation

{ TForm1 }

//Função do Menu: "Ferramentas" -> "Resolver Estrutura"

procedure TForm1.MenuItem9Click(Sender: TObject);

var

x: double;

ai, af, i, j, c, k, ik, ie, de, idat, te: integer;

ind: array of integer;

Kg, Ke, Fe, gauss, temp, F, d, B0: TMatriz;

Eq: CEquacao;

item: TMenuItem;

begin

MenuItem8Click(nil);

//Verificação se BL>20 nas cascas cilíndricas

Data := StrGrdToArr(5);

for idat := 0 to High(Data) do

begin

k := StrToInt(StringGrid5.Cells[0, idat+1]);

PrepDGeom(El, k-1);

PrepData(idat, 0, 0);

A - 18

B := (Power(3, 0.25)*Power(1-Power(v, 2), 0.25))/Sqrt(h*r1);

if B*L>20 then

begin

j := MessageDlg('O elemento nº'+IntToStr(k)+

' (casca cílindrica) apresenta um factor B*L>20.'+#13#10+

'Esta situação poderá resultar em erros numéricos significativos pelo que se aconselha

a subdividir o elemento mencionado.'

+#13#10+'Deseja cancelar a operação?', mtWarning, mbYesNo, 0);

if j = mrYes then

Exit;

end;

end;

temp := StrGrdToArr(2);

SetLength(nos, Length(nos), 5);

for i := 0 to High(temp) do

for j := 3 to High(temp[0]) do

nos[i, j-1] := temp[i, j];

//Identificação dos graus de liberdade na base de dados dos nós

k := 0;

for i := 0 to High(nos) do

for j := 2 to High(nos[0]) do

if nos[i, j]<0 then

nos[i, j] := 0

else

begin

k := k+1; //k define o número de graus de liberdade

nos[i, j] := k;

end;

SetLength(F, k, 1);

SetLength(d, k, 1);

SetLength(Kg, k, k);

//Espalhamento das cargas concentradas ou linearmente distribuídas


for j := 0 to 2 do

if nos[i, j+2]>0 then

if nos[i, 0] = 0 then

F[Round(nos[i, j+2])-1, 0] := temp[i, j]

else

F[Round(nos[i, j+2])-1, 0] := temp[i, j]*2*PI*nos[i, 0];

SetLength(Ke, 6, 6);

SetLength(ind, 6);

//Matriz de abcissas e pesos a utilizar na integração pela regra de quadratura de Gauss

gauss := gauleg(StrToInt(Edit9.Text));

idat := 0;

k := -1;

//Espalhamento directo da matriz de rigidez do elemento de anel



begin

k := k+1;

E := Data[idat, 3];

r1 := nos[Round(Data[idat, 0])-1, 0];

for i := 0 to 1 do

begin

ie := Round(nos[Round(Data[idat, 0])-1, 2+2*i])-1;

if ie>-1 then

Kg[ie, ie] := Kg[ie, ie]+E*Data[idat, 1+i]*2*Pi/r1;

end;

end;

A - 19

{

Construção da matriz de rigidez global e vector de forças nodais global

Os elementos podem ser identificados pela variável te:

te=4: Elementos de cascas em anel;

te=5: Elementos de cascas cilíndricas';

te=6: MEF: Troços;

te=7: MEF: Arcos;

}

for te := 4 to High(Mdados) do

begin

Data := StrGrdToArr(te);


begin

if te<6 then

af := 0

else

af := Round(Data[idat, 9])-1;

for ai := 0 to af do

begin

k := k+1;

PrepDGeom(El, k);

PrepData(idat, ai, af);

//Cálculo do vector de forças nodais elementar

if te = 4 then

Fe := AnelFe();

if te = 5 then

begin


VecY(L);

Fe := CilFe();

end;

if te>5 then

Fe := MatFe();

//Transformação de coordenadas locais para globais

Fe := MultMatTr(MatTr(2), Fe);

Ke := nil;

SetLength(Ke, 6, 6);

//Cálculo da matriz de rigidez elementar

if te = 4 then

Ke := AnelKe();

if te = 5 then

Ke := CilKe();

if te>5 then

for i := 0 to High(gauss) do

begin

x := gauss[i, 0]*L/2+L/2;

MatBe(x);

Ke := AddSubMat(Ke, MultEsc(MultMat(MultMatTr(Be, Dg), Be),

gauss[i, 1]*L*PI*(r1+((-r1+r2)*x)/L)), 1);

end;

//Transformação de coordenadas locais para globais

Ke := MultMatTr(MatTr(2), Ke);

Ke := MultMat(Ke, MatTr(2));

//Construção de um vector de índices para espalhamento de matrizes ou vectores

for i := 0 to 1 do

for j := 2 to 4 do

ind[i*3+j-2] := Round(nos[El[k, i], j])-1;

//Espalhamento da matriz de rigidez

for i := 0 to 5 do

for j := 0 to 5 do

A - 20

if (ind[i]>-1) and (ind[j]>-1) then

Kg[ind[i], ind[j]] := Kg[ind[i], ind[j]]+Ke[i, j];

//Espalhamento do vector de forças nodais

for i := 0 to 5 do

begin

if ind[i]>-1 then

F[ind[i], 0] := F[ind[i], 0]+Fe[i, 0];

Fe[i, 0] := 0;

end;

end;

end;

end;

//Resolução do sistema de equações lineares pelo método de eliminação de Gauss

Eq := CEquacao.Cria(Length(Kg), 1, False);

Eq.ZeraMatriz;

for i := 0 to High(Kg) do

for j := 0 to High(Kg) do

Eq.AtribuiAij(i+1, j+1, Kg[i, j]);

for i := 0 to High(F) do

Eq.AtribuiBj(i+1, F[i, 0]);

Eq.Solve;

for i := 0 to High(d) do

d[i, 0] := Eq.Solucao(1)[i+1];

//Adição dos deslocamentos nodais à base de dados dos Nós


for j := 2 to 4 do

if nos[i, j]>0 then

nos[i, j] := d[Round(nos[i, j])-1, 0];

d := nil;

SetLength(d, 6, 1);

j := ires[0];

for i := 3 to 4 do

j := j+ires[i];

resf := nil;

SetLength(resf, j, 11);

idat := 0;

k := -1;

//Construção da base de dados dos Resultados Gerais

//Elemento de anel



begin

k := k+1;

E := Data[idat, 3];

r1 := nos[Round(Data[idat, 0])-1, 0];

resf[k, 0] := k+1;

resf[k, 1] := r1;

resf[k, 2] := nos[Round(Data[idat, 0])-1, 1];

//Obtenção dos deslocamentos

for i := 2 to 4 do

resf[k, i+1] := nos[Round(Data[idat, 0])-1, i];

//Obtenção dos esforços

for i := 0 to 1 do

A - 21

resf[k, 7+2*i] := nos[El[k, 0], 2+2*i]*E*Data[idat, 1+i]/r1;

end;

ie := ires[0];

ik := -1;

de := ires[5];

//Obtenção dos resultados dos restantes elementos


begin


if te = 6 then

begin

ie := ie+ires[3];

de := ires[6];

ik := ik-ires[1];

end;


begin

if te<6 then

af := 0

else



begin

k := k+1;

ik := ik+1;

PrepDGeom(El, k);

PrepData(idat, ai, af);

//Construção do vector de deslocamentos de um elemento

for i := 0 to 1 do

for j := 2 to 4 do

d[i*3+j-2, 0] := nos[El[k, i], j];

//Preparação para obtenção de resultados da casca em anel e cilíndrica

if te = 4 then

begin

SetLength(temp, 4, 1);

temp[0, 0] := d[1, 0]-(3*Qy1*Power(r1, 4)*(1-Power(v, 2)))/(16.*E*Power(h, 3));




B0 := MultMat(MatInvC(), temp);

end;

if te = 5 then

begin

SetLength(B0, 4, 1);


VecY(L);

Ke := CilKe();

temp := MultMat(MatTr(2), d);

Fe := MultMat(Ke, AddSubMat(temp, Cilqp(), -1));

Fe := MultEsc(Fe, 1/(2*Pi*r1));

B0[0, 0] := temp[1, 0]-(r1*(2*Qy1*r1-2*Fe[0, 0]*v+L*Qx*v))/(2.*E*h);

B0[1, 0] := temp[2, 0]+((r1*(Qy1*r1-Qy2*r1+L*Qx*v))/(E*h*L));

B0[2, 0] := -Fe[2, 0];

B0[3, 0] := Fe[1, 0];

end;

A - 22

for c := 0 to de-1 do

begin

resf[ie+de*ik+c, 0] := k+1;

resf[ie+de*ik+c, 1] := r1+(r2-r1)*c/(de-1);

resf[ie+de*ik+c, 2] :=

nos[El[k, 0], 1]+(nos[El[k, 1], 1]-nos[El[k, 0], 1])*c/(de-1);

x := c*L/(de-1);

if te = 4 then

if (r1 = 0) and (x = 0) then

begin

//Obtenção dos deslocamentos de elementos de casca em anel

for i := 0 to 2 do

resf[ie+de*ik+c, 3+i] := d[i, 0];

//Obtenção dos esforços de elementos de casca em anel

resf[ie+de*ik+c, 6] := (d[3, 0]*E*h)/(r2-r2*v);

resf[ie+de*ik+c, 7] := (d[3, 0]*E*h)/(r2-r2*v);

if (Mdados[2].Cells[2, El[k, 0]+1] = '0') and

(Mdados[2].Cells[5, El[k, 0]+1] = '0') then

begin


-(E*Power(h, 3))/(6.*Power(r2, 2)*(-1+v))*B0[2, 0];


-(E*Power(h, 3))/(6.*Power(r2, 2)*(-1+v))*B0[2, 0];

resf[ie+de*ik+c, 10] := 0;

end;

end

else

begin

temp := MultMat(ResAnel(x), B0);

//Obtenção dos deslocamentos de elementos de casca em anel


(-(d[3, 0]*Power(r1, 2)*r2)+d[0, 0]*r1*Power(r2, 2))/

((-Power(r1, 2)+Power(r2, 2))*(r1+((-r1+r2)*x)/L))+

((-(d[0, 0]*r1)+d[3, 0]*r2)*(r1+((-r1+r2)*x)/L))/(-Power(r1, 2)+Power(r2, 2));


temp[0, 0]-(-3*Qy1*(1-Power(v, 2))*Power(r1+((-r1+r2)*x)/L,

4))/(16.*E*Power(h, 3));


temp[1, 0]-(-3*Qy1*(1-Power(v, 2))*Power(r1+((-r1+r2)*x)/L, 3))/(4.*E*Power(h,

3));

//Obtenção dos esforços de elementos de casca em anel

resf[ie+de*ik+c, 6] := -((E*h*(-(d[3, 0]*r2*(2*Power(L, 2)*Power(r1, 2)-

2*L*r1*(r1-r2)*

(1+v)*x+Power(r1-r2, 2)*(1+v)*Power(x, 2)))+d[0, 0]*r1*

(Power(L, 2)*(-(Power(r2, 2)*(-1+v))+Power(r1, 2)*(1+v))-

2*L*r1*(r1-r2)*(1+v)*x+Power(r1-r2, 2)*(1+v)*Power(x, 2))))/

((r1-r2)*(r1+r2)*(-1+Power(v, 2))*Power(L*r1+(-r1+r2)*x, 2)));

resf[ie+de*ik+c, 7] := -((E*h*(-(d[3, 0]*r2*(2*Power(L, 2)*Power(r1, 2)*v-

2*L*r1*(r1-r2)*

(1+v)*x+Power(r1-r2, 2)*(1+v)*Power(x, 2)))+d[0, 0]*r1*

(Power(L, 2)*(Power(r2, 2)*(-1+v)+Power(r1, 2)*(1+v))-2*L*r1*

(r1-r2)*(1+v)*x+Power(r1-r2, 2)*(1+v)*Power(x, 2))))/

((r1-r2)*(r1+r2)*(-1+Power(v, 2))*Power(L*r1+(-r1+r2)*x, 2)));


temp[2, 0]+(Qy1*(3+v)*Power(r1+((-r1+r2)*x)/L, 2))/16;


temp[3, 0]+(Qy1*(1+3*v)*Power(r1+((-r1+r2)*x)/L, 2))/16;

resf[ie+de*ik+c, 10] := temp[4, 0]-Qy1*(r1+((-r1+r2)*x)/L)/2;

end;

A - 23

if te = 5 then

begin

VecY(x);

temp := MultMat(ResCil(), B0);

//Obtenção dos deslocamentos de elementos de casca cilíndrica

resf[ie+de*ik+c, 3] := -(temp[1, 0]+(r1*(Power(L, 2)*Qx*v+2*(-Qy1+Qy2)*r1*x+2*L*

(Qy1*r1-v*(Fe[0, 0]+Qx*x))))/(2.*E*h*L));


(v/r1)*temp[0, 0]+(x*(-2*Fe[0, 0]*L+Power(L, 2)*Qx+(-Qy1+Qy2)*r1*v*x+L*

(2*Qy1*r1*v-Qx*x)))/(2.*E*h*L);


temp[2, 0]-((r1*(Qy1*r1-Qy2*r1+L*Qx*v))/(E*h*L));

//Obtenção dos esforços de elementos de casca cilíndrica

resf[ie+de*ik+c, 6] := Qx*(L/2-x)-Fe[0, 0];

resf[ie+de*ik+c, 7] := -(E*h/r1)*temp[1, 0]+r1*(-Qy1+((Qy1-Qy2)*x)/L);

resf[ie+de*ik+c, 8] := temp[4, 0];

resf[ie+de*ik+c, 9] := v*temp[4, 0];


end;

if te>5 then

begin

//Obtenção dos deslocamentos de elementos do MEF

temp := MultMat(MatPsi(x), MultMat(MatTr(2), d));

temp := MultMatTr(MatTr(1), temp);

for i := 0 to 2 do

resf[ie+de*ik+c, 3+i] := temp[i, 0];

//Obtenção dos esforços de elementos do MEF

if (r1+((-r1+r2)*x)/L = 0) then

MatBeSing() //Função que lida com as singularidades na matriz B

else

begin

MatBe(x);

Be := MultMat(Be, MatTr(2));

end;

temp := MultMat(Dg, MultMat(Be, d));

for i := 0 to 3 do

resf[ie+de*ik+c, 6+i] := temp[i, 0];

temp := MultMat(MatVe(x), MultMat(MatTr(2), d));


end;

end;

end;

end;

end;

if ComboBox1.Items.Count = 1 then

begin

ComboBox1.Items.Add('Estrutura deformada');

ComboBox1.Items.Add('Esforços: Nx');

ComboBox1.Items.Add('Esforços: Nt');

ComboBox1.Items.Add('Esforços: Mx');

ComboBox1.Items.Add('Esforços: Mt');

ComboBox1.Items.Add('Esforços: Vx');

end;

A - 24

for i := 18 to 21 do

begin

item := FindComponent(Format('MenuItem%d', [i])) as TMenuItem;

item.Enabled := True;

end;

end;

//Função do Menu: "Ferramentas" -> "Desenhar Estrutura"


var

i, j, a, c, ie, k, de, inos: integer;

r, phi, phi0, rC, zC: double;

item: TMenuItem;

temp: TMatriz;

begin

zoomres := 1;

//Construção da base de dados dos nós e elementos

nos := StrGrdToArr(1);

El := nil;

ires[0] := StrToInt(Edit2.Text);

//Elemento de anel

SetLength(El, ires[0], 2);

for i := 0 to StringGrid3.RowCount-2 do

begin

El[i, 0] := StrToInt(StringGrid3.Cells[1, i+1])-1;

El[i, 1] := StrToInt(StringGrid3.Cells[1, i+1])-1;

end;

//Casca cilíndrica e Casca em anel

for j := 4 to 5 do

begin

ie := Length(El);

SetLength(El, ie+StrToInt(Nel[j-1].Text), 2);

for i := 0 to Mdados[j].RowCount-2 do

for a := 0 to 1 do

El[ie+i, a] := StrToInt(Mdados[j].Cells[1+a, i+1])-1;

end;

ires[1] := Length(El)-ires[0];

//Troços e Arcos do MEF

for a := 6 to 7 do

begin

for i := 0 to Mdados[a].RowCount-2 do

begin

ie := Length(El);

c := StrToInt(Mdados[a].Cells[1, i+1])-1;

k := StrToInt(Mdados[a].Cells[2, i+1])-1;

de := StrToInt(Mdados[a].Cells[10, i+1]);

SetLength(El, ie+de, 2);

inos := Length(nos);

SetLength(nos, inos+de-1, 2);

for j := 0 to de-1 do

begin

El[ie+j, 0] := inos+j-1;

El[ie+j, 1] := inos+j;

end;

El[ie, 0] := c;

El[ie+de-1, 1] := k;

if a = 7 then

begin

rC := StrToFloat(Mdados[a].Cells[11, i+1]);

A - 25

zC := StrToFloat(Mdados[a].Cells[12, i+1]);

r := sqrt(sqr(nos[c, 0]-rc)+sqr(nos[c, 1]-zc));

phi0 := ArcCos((nos[c, 0]-rc)/r);

phi := ArcCos((nos[k, 0]-rc)/r)-phi0;

end;

for j := 1 to de-1 do

begin

if a = 6 then

begin

nos[inos+j-1, 0] := nos[c, 0]+(nos[k, 0]-nos[c, 0])*j/de;

nos[inos+j-1, 1] := nos[c, 1]+(nos[k, 1]-nos[c, 1])*j/de;

end else

begin

nos[inos+j-1, 0] := Cos(phi0+j*phi/de)*r+rC;

nos[inos+j-1, 1] := Sin(phi0+j*phi/de)*r+zC;

end;

end;

end;

end;

ires[2] := Length(El)-ires[1]-ires[0];



ires[3] := ires[1]*ires[5];

ires[4] := ires[2]*ires[6];

DesenharEstrutura();

ComboBox1.Items.Clear;

ComboBox1.Items.Add('Estrutura indeformada');

ComboBox1.ItemIndex := 0;


begin


item.Enabled := False;

end;

MenuItem17.Enabled := True;

MenuItem20.Enabled := True;

end;

//Desenho da estrutura

procedure TForm1.DesenharEstrutura();

var

max, min: TVector;

tmpint: MatInt;

sw: double;

i, j, c: integer;

Bitmap: TBitmap;

begin

if CheckGroup1.Checked[5] then

SetLength(nos, Length(nos)+1, Length(nos[0]));

max := MatMax(nos);

min := MatMin(nos);

dx := max[0]-min[0];

dy := max[1]-min[1];

//Função que permite manter as proporções ou esticar o desenho

if RadioButton1.Checked then

begin

if dx>dy then

A - 26

dy := dx*Image1.Height/Image1.Width

else

dx := dy*Image1.Width/Image1.Height;

end;

//Definição da quantidade de área a utilizar para o desenho da estrutura

if dx<>0 then

dx := 0.75*Image1.Width/dx;

if dy<>0 then

dy := 0.75*Image1.Height/dy;

//Coordenadas de referência do desenho

xi := 0.5*(Image1.Width-(max[0]+min[0])*dx);

yi := 0.5*(Image1.Height+(max[1]+min[1])*dy);

Bitmap := TBitmap.Create;

Bitmap.Width := Image1.Width;

Bitmap.Height := Image1.Height;

Bitmap.Canvas.Rectangle(0, 0, Image1.Width, Image1.Height);

Image1.Picture.Bitmap := Bitmap;

//Desenho dos eixos das coordenadas globais

ca := 1;

sa := 0;

DSeta(10, Image1.Height-11);

ca := 0;

sa := 1;

DSeta(11, Image1.Height-10);

Image1.Canvas.TextOut(42, Image1.Height-20, 'r');

Image1.Canvas.TextOut(10, Image1.Height-52, 'z');

//Desenho do símbolo de axissimetria

Bitmap.Handle := LoadBitmap(hInstance, 'SimSymb');

Image1.Canvas.Draw(Round(xi-Bitmap.Width div 2), 10, Bitmap);

Bitmap.Free;

with Image1.Canvas do

begin

i := 5;

while i<Image1.Height-25 do

begin

Line(Round(xi), i, Round(xi), i+10);

i := i+20;

end;

//Desenho dos elementos


Pixels[round(xi+nos[i, 0]*dx), round(yi-nos[i, 1]*dy)] := clBlack;

for i := ires[0] to High(El) do

Line(round(xi+nos[El[i, 0], 0]*dx), round(yi-nos[El[i, 0], 1]*dy),

round(xi+nos[El[i, 1], 0]*dx), round(yi-nos[El[i, 1], 1]*dy));

//Representação do elemento de anel

Brush.Color := $0080FF;


for i := 0 to ires[0]-1 do

Ellipse(round(xi+nos[El[i, 0], 0]*dx)-5,

round(yi-nos[El[i, 0], 1]*dy)-5, round(xi+nos[El[i, 0], 0]*dx)+

6, round(yi-nos[El[i, 0], 1]*dy)+6);

Brush.Color := clWhite;

Brush.Style := bsClear;

Font.Color := clBlue;

//Desenho da origem


A - 27

begin

Pen.Color := clBlue;

Pen.Width := 5;

TextOut(Round(xi-14), Round(yi+4), 'O');

Line(Round(xi), Round(yi), Round(xi), Round(yi));

SetLength(nos, Length(nos)-1, Length(nos[0]));

Pen.Width := 1;

end;

//Desenho dos nós


begin

Pen.Color := clBlue;

Pen.Width := 5;


Line(round(xi+nos[i, 0]*dx), round(yi-nos[i, 1]*dy),

round(xi+nos[i, 0]*dx),

round(yi-nos[i, 1]*dy));

Pen.Color := clBlack;

Pen.Width := 1;

end;

//Desenho da identificação dos nós



TextOut(round(xi+nos[i, 0]*dx)+5, round(yi-nos[i, 1]*dy)+

5, IntToStr(i+1));

Font.Color := clBlack;

//Desenho da identificação dos elementos


begin

for i := 0 to ires[0]-1 do

TextOut(round(xi+nos[El[i, 0], 0]*dx)-10,

round(yi-nos[El[i, 0], 1]*dy)-20, IntToStr(i+1));

for i := ires[0] to High(El) do

begin

PrepDGeom(El, i);

TextOut(round(xi+(nos[El[i, 0], 0]+nos[El[i, 1], 0])/2*dx+5*Abs(sa)),

round(yi-(nos[El[i, 0], 1]+nos[El[i, 1], 1])/2*dy+5*Abs(ca)),

IntToStr(i+1));

end;

end;

end;

Image1.Canvas.Pen.Color := $00C000;

//Desenho da orientação dos elementos


for c := 4 to High(Mdados) do

with Mdados[c] do

begin

SetLength(tmpint, RowCount-1, 2);

for j := 0 to 1 do

for i := 0 to RowCount-2 do

tmpint[i, j] := StrToInt(Cells[j+1, i+1])-1;

for i := 0 to high(tmpint) do

begin

PrepDGeom(tmpint, i);

DSeta(xi+dx*(nos[tmpint[i, 0], 0]+nos[tmpint[i, 1], 0])/

2, yi-dy*(nos[tmpint[i, 0], 1]+nos[tmpint[i, 1], 1])/2);

sw := ca;

ca := -sa;

sa := sw;

DSeta(xi+dx*(nos[tmpint[i, 0], 0]+nos[tmpint[i, 1], 0])/

2, yi-dy*(nos[tmpint[i, 0], 1]+nos[tmpint[i, 1], 1])/2);

A - 28

end;

end;

TabSheet8.TabVisible := True;

TabSheet8.Show;

end;

//Função do Menu: "Ficheiro" -> "Nova Estrutura"


var

i: integer;

item: TMenuItem;

begin

for i := 1 to 6 do

Nel[i].Text := '0';

Nel[1].Text := '2';

TabSheet1.Show;

TabSheet8.TabVisible := False;


begin


item.Enabled := False;

end;

end;

//Função do Menu: "Ficheiro" -> "Abrir Ficheiro de Dados"


var

d, i, a, b: integer;

ns, txt, loc: string;

myFile: TextFile;

begin

i := 0;

loc := Abrir('DAT (Ficheiro de Dados) (*.dat)|*.dat', self);

if loc<>'' then

begin

AssignFile(myFile, loc);

Reset(myFile);

ReadLn(myFile, txt);

while Length(txt)>0 do

begin

i := i+1;

d := Pos(';', txt);

ns := Copy(txt, 0, d-1);

Nel[i].Text := ns;

txt := Copy(txt, d+1, MaxInt);

end;

i := 0;

while not EOF(myFile) do

begin

i := i+1;

for b := 1 to Mdados[i].RowCount-1 do

begin

a := 0;

ReadLn(myFile, txt);

while Length(txt)>0 do

begin

a := a+1;

d := Pos(';', txt);

ns := Copy(txt, 0, d-1);

Mdados[i].Cells[a, b] := ns;

txt := Copy(txt, d+1, MaxInt);

end;

end;

end;

CloseFile(myFile);

A - 29

end;

end;

//Função do Menu: "Ficheiro" -> "Guardar Ficheiro de Dados"


var

i, a, b: integer;

texto: string;

begin

texto := '';

for i := 1 to High(Nel) do

texto := texto+Nel[i].Text+';';

texto := texto+#13#10;

for i := 1 to High(Mdados) do

with Mdados[i] do

for b := 1 to RowCount-1 do

begin

for a := 1 to ColCount-1 do

texto := texto+Cells[a, b]+';';


end;

Guardar('DAT (Ficheiro de Dados) (*.dat)|*.dat', 'dat', texto, self);

end;

//Função do Menu: "Ficheiro" -> "Exportar Imagem em Display"


var

saveDialog: TSaveDialog;

img: TPortableNetworkGraphic;

begin

saveDialog := TSaveDialog.Create(self);

saveDialog.Title := 'Guardar como';

saveDialog.InitialDir := GetCurrentDir;

saveDialog.Filter := 'Ficheiro de PNG (*.png)|*.png';

saveDialog.DefaultExt := 'png';

saveDialog.FilterIndex := 1;

if saveDialog.Execute then

begin

if FileExists(SaveDialog.FileName) then

if MessageDlg(Format('O ficheiro "%s" já existe. Pretende substituir?',

[ExtractFileName(SaveDialog.FileName)]), mtConfirmation, mbYesNo, 0) = idNo then

Exit;

img := TPortableNetworkGraphic.Create;

with Image1.Picture do

img.CreateFromBitmapHandles(Bitmap.Handle, Bitmap.MaskHandle,

Bounds(0, 0, Width, Height));

img.SaveToFile(SaveDialog.FileName);

img.Free;

end;

saveDialog.Free;

end;

//Função do Menu: "Ficheiro" -> "Exportar Resultados" -> "Nós e Deslocamentos Nodais"


begin

Guardar('CSV (Ficheiro de CSV) (*.csv)|*.csv', 'csv',

MatOutTC(nos, 'Nó;r;z;dr;dz;dt'), self);

end;

//Função do Menu: "Ficheiro" -> "Exportar Resultados" -> "Elementos"


var

temp: TMatriz;

i, j: integer;

A - 30

begin

SetLength(temp, Length(El), 3);

for i := 0 to High(El) do

for j := 0 to High(El[0]) do

temp[i, j+1] := El[i, j]+1;


temp[i, 0] := i+1;

Guardar('CSV (Ficheiro de CSV) (*.csv)|*.csv', 'csv', PrepOut(0, temp), self);

end;

//Função do Menu: "Ficheiro" -> "Exportar Resultados" -> "Resultados Gerais"


begin

Guardar('CSV (Ficheiro de CSV) (*.csv)|*.csv', 'csv', PrepOut(1, resf), self);

end;

//Função do Menu: "Ficheiro" -> "Sair"


begin

Application.Terminate;

end;

//Função do Menu: "Ver" -> "Nós" e "Cascas", "Ferramentas" -> "Opções"

procedure TForm1.MenuDados(Sender: TObject);

begin

PageControl1.ActivePageIndex := StrToInt(Copy((Sender as TMenuItem).Name, 10, 1));

end;

//Função do Menu: "Ver" -> "Estrutura Indeformada", "Estrutura Deformada", "Esforços na

Estrutura"

procedure TForm1.MenuResultados(Sender: TObject);

begin

TabSheet8.Show;

ComboBox1.ItemIndex := StrToInt(Copy((Sender as TMenuItem).Name, 10, 1));

ComboBox1Change(nil);

end;

//Função do Menu: "Ajuda" -> "Manual de Utilização"


begin

if FileExists('Manual.pdf') = False then

ShowMessage('Não foi possível localizar o manual de utilização: "Manual.pdf".')

else

ShellExecute(Form1.Handle, PChar('open'), PChar('Manual.pdf'),

PChar(''), PChar(''), 1);

end;

//Função do Menu: "Ajuda" -> "Sobre"


begin

Form2.Show;

end;

//Preparação de dados: propriedades dos materiais e carregamento

procedure TForm1.PrepData(idat, ai, af: integer);

begin

h := Data[idat, 2];

v := Data[idat, 3];

E := Data[idat, 4];

MatD();

if Length(Data[0]) = 6 then

begin

Qx := 0;

Qy1 := -Data[idat, 5]*ca;

Qy2 := Qy1;

end

else

A - 31

begin

Qx := -Data[idat, 8]*sa+Data[idat, 7];

Qy1 := -Data[idat, 8]*ca;

Qy2 := Qy1;

if af = 0 then

begin

Qy1 := Qy1+Data[idat, 5];

Qy2 := Qy2+Data[idat, 6];

end

else

begin

Qy1 := Qy1+Data[idat, 5]+(Data[idat, 6]-Data[idat, 5])*ai/(af+1);

Qy2 := Qy2+Data[idat, 5]+(Data[idat, 6]-Data[idat, 5])*(ai+1)/(af+1);

end;

end;

end;

//Preparação de dados geométricos

procedure TForm1.PrepDGeom(m: MatInt; k: integer);

begin

r1 := nos[m[k, 0], 0];

r2 := nos[m[k, 1], 0];

L := sqrt(sqr(r1-r2)+sqr(nos[m[k, 0], 1]-nos[m[k, 1], 1]));

ca := (r2-r1)/L;

sa := (nos[m[k, 1], 1]-nos[m[k, 0], 1])/L;

end;

//Incremento da escala de deslocamentos e diagramas de esforços

procedure TForm1.Button1Click(Sender: TObject);

begin

if ComboBox1.ItemIndex>0 then

begin

zoomres := zoomres*1.25;

Edit10.Text := NumOut(1.25*StrToFloat(Edit10.Text), 6);



end;

end;

//Diminuição da escala de deslocamentos e diagramas de esforços


begin


begin

zoomres := zoomres*0.8;




end;

end;

//Activação de escala automática de deslocamentos e diagramas de esforços


begin


begin

zoomres := 1;

RadioButton5.Checked := True;



end;

end;

//Actualização da imagem em display


begin


end;

//Função da caixa de selecção "Exibir" que permite a visualização da estrutura,

//a deformada e os esforços da estrutura.

A - 32

procedure TForm1.ComboBox1Change(Sender: TObject);

var

i, j, id, k, np, ie, te, af: integer;

valor: double;

max, min, esc: TVector;

temp: TMatriz;

ind: array of integer;

add: boolean;

begin

DesenharEstrutura();

id := ComboBox1.ItemIndex;

ie := 0;

//Determinação da escala da deformada e diagramas de esforços

if id>0 then

begin

np := ires[6];

max := MatMax(resf);

min := MatMin(resf);

for i := 0 to High(max) do

if max[i]<0 then

max[i] := 0;

for i := 0 to High(min) do

if min[i]>0 then

min[i] := 0;

SetLength(esc, Length(max));

for i := 0 to High(max) do

if (max[i]-min[i])<>0 then

esc[i] := 40/(max[i]-min[i]);

if (esc[3]<esc[4]) and (esc[3]<>0) then

esc[4] := esc[3];

if (esc[4]<esc[3]) and (esc[4]<>0) then

esc[3] := esc[4];

if esc[3] = 0 then

esc[3] := esc[4];

if esc[4] = 0 then

esc[4] := esc[3];

for i := 0 to High(esc) do

esc[i] := esc[i]*zoomres;

if RadioButton6.Checked = True then

for i := 3 to 4 do

esc[i] := StrToFloat(Edit10.Text);

if RadioButton8.Checked = True then

for i := 6 to 10 do

esc[i] := StrToFloat(Edit11.Text);

for i := 8 to 9 do

esc[i] := -esc[i]; //Convenção de sinais dos momentos

Image1.Canvas.Pen.Color := clRed;

end;

//Desenho da deformada da estrutura

if id = 1 then

begin

Edit10.Text := NumOut(abs(esc[3]), 6);

for j := 0 to 1 do

begin

np := ires[j+5];

ie := ie+ires[j*3];

for k := 0 to ires[j+1]-1 do

begin

A - 33

Image1.Canvas.MoveTo(round(xi+resf[ie+k*np, 1]*dx+

resf[ie+k*np, 3]*esc[3]),

round(yi-resf[ie+k*np, 2]*dy-resf[ie+k*np, 4]*esc[4]));

for i := 1 to np-1 do

Image1.Canvas.LineTo(round(xi+resf[ie+k*np+i, 1]*dx+resf[ie+k*np+i, 3]*esc[3]),

round(yi-resf[ie+k*np+i, 2]*dy-resf[ie+k*np+i, 4]*esc[4]));

end;

end;

end;

//Desenho dos diagramas de esforços da estrutura

if id>1 then

begin

id := id+4;

Edit11.Text := NumOut(abs(esc[id]), 6);

for j := 0 to 1 do

begin

np := ires[j+5];

ie := ie+ires[j*3];

for k := 0 to ires[j+1]-1 do

begin

PrepDGeom(El, Round(resf[ie+k*np, 0])-1);

Image1.Canvas.Line(round(xi+resf[ie+k*np, 1]*dx),

round(yi-resf[ie+k*np, 2]*dy),

round(xi+resf[ie+k*np, 1]*dx-resf[ie+k*np, id]*esc[id]*sa),

round(yi-resf[ie+k*np, 2]*dy-resf[ie+k*np, id]*esc[id]*ca));

for i := 1 to np-1 do

begin

Image1.Canvas.LineTo(

round(xi+resf[ie+k*np+i, 1]*dx-resf[ie+k*np+i, id]*esc[id]*sa),

round(yi-resf[ie+k*np+i, 2]*dy-resf[ie+k*np+i, id]*esc[id]*ca));

Image1.Canvas.Line(round(xi+resf[ie+k*np+i, 1]*dx-resf[ie+k*np+i,

id]*esc[id]*sa), round(yi-resf[ie+k*np+i, 2]*dy-resf[ie+k*np+i, id]*esc[id]*ca),

round(xi+resf[ie+k*np+i, 1]*dx), round(yi-resf[ie+k*np+i, 2]*dy));

Image1.Canvas.MoveTo(

round(xi+resf[ie+k*np+i, 1]*dx-resf[ie+k*np+i, id]*esc[id]*sa),

round(yi-resf[ie+k*np+i, 2]*dy-resf[ie+k*np+i, id]*esc[id]*ca));

end;

end;

end;

k := 0;

SetLength(temp, 2, 0);

SetLength(ind, 2);

ie := ires[0];

np := ires[5];

//Função de desenho dos "Valores relevantes"


begin

Image1.Canvas.Font.Bold := True;

for i := 0 to ie-1 do

if (id = 7) or (id = 9) then

Image1.Canvas.TextOut(round(xi+resf[i, 1]*dx)+10,

round(yi-resf[i, 2]*dy)-25, NumOut(resf[i, id], 6));

Image1.Canvas.Font.Bold := False;


begin


if te = 6 then

begin

ie := ie+ires[3];

np := ires[6];

k := k-ires[1];

end;


begin

if te<6 then

A - 34

af := 1

else

af := Round(Data[idat, 9]);

begin

ind[0] := ie+k*np;

k := k+af;

ind[1] := ie+k*np-1;

temp[0] := MatMinPos(resf, ind[0], ind[1]);

temp[1] := MatMaxPos(resf, ind[0], ind[1]);

for i := 0 to 1 do

begin

add := True;

for j := 0 to High(ind) do

if temp[i, id] = ind[j] then

add := False;

if add then

begin

SetLength(ind, Length(ind)+1);

ind[High(ind)] := Round(temp[i, id]);

end;

end;

for i := 0 to High(ind) do

begin

SetTextAlign(Image1.Canvas.Handle, TA_CENTER);

PrepDGeom(El, Round(resf[ind[i], 0])-1);

if i = 0 then

if ca>0 then

SetTextAlign(Image1.Canvas.Handle, TA_LEFT)

else

SetTextAlign(Image1.Canvas.Handle, TA_RIGHT);

if i = 1 then

if ca>0 then

SetTextAlign(Image1.Canvas.Handle, TA_RIGHT)

else

SetTextAlign(Image1.Canvas.Handle, TA_LEFT);

valor := resf[ind[i], id]*esc[id];

if valor<0 then

valor := valor-15

else

valor := valor+15;

Image1.Canvas.TextOut(round(xi+(resf[ind[i], 1])*dx-valor*sa),

round(yi-resf[ind[i], 2]*dy-valor*ca-7),

NumOut(resf[ind[i], id], 6));

end;

end;

end;

end;

end;

end;

SetTextAlign(Image1.Canvas.Handle, TA_LEFT);

end;

//Identificação do tipo de elemento

function TForm1.TEl(i: integer): string;

begin

case i of

3: Result := 'Elementos de anel';

4: Result := 'Cascas em anel';

5: Result := 'Cascas cilíndricas';

6: Result := 'MEF: Troço nº ';

7: Result := 'MEF: Arco nº ';

end;

end;

//Preparação de saída de dados dos resultados gerais e elementos

A - 35

function TForm1.PrepOut(tr: integer; M: TMatriz): string;

var

txt: string;

i, j, te, ie, de, ik, ai, af, c: integer;

temp: array[0..1, 0..1] of string;

begin

temp[0, 0] := 'Elemento;Nó';

temp[1, 0] := 'Elemento;Nó inicial;Nó final';

temp[0, 1] := 'Elemento;r;z;dr;dz;dt;Nt;Mt';

temp[1, 1] := 'Elemento;r;z;dr;dz;dt;Nx;Nt;Mx;Mt;Vx';

txt := '';

ie := ires[0];

de := 1;

if ie>0 then

begin

txt := txt+TEl(3)+#13#10+temp[0, tr]+#13#10;

for i := 0 to ie-1 do

begin

for j := 0 to 1+tr*4 do

txt := txt+FloatToStr(M[i, j])+';';

if tr = 1 then

txt := txt+FloatToStr(M[i, 7])+';'+FloatToStr(M[i, 9])+';';

end;

txt := txt+#13#10;

end;

ik := -1;

if tr = 1 then

de := ires[5];


begin


if te = 6 then

begin

if tr = 1 then

begin

ie := ie+ires[3];

de := ires[6];

end

else

ie := ie+ires[1];

ik := ik-ires[1];

end;

if (High(Data)>=0) and (te<6) then

txt := txt+#13#10+TEl(te)+#13#10+temp[1, tr]+#13#10;


begin

if te>5 then

txt := txt+#13#10+TEl(te)+IntToStr(idat+1)+#13#10+temp[1, tr]+#13#10;

if te<6 then

af := 0

else



begin

ik := ik+1;

for c := 0 to de-1 do

begin

for j := 0 to High(M[0]) do

txt := txt+FloatToStr(M[ie+de*ik+c, j])+';';

txt := txt+#13#10;

end;

end;

A - 36

end;

end;

Result := txt;

end;

//Preparação de saída de dados de "Nós"

function TForm1.MatOutTC(M: TMatriz; title: string): string;

var

texto: string;

i, j: integer;

begin

texto := '';

texto := title+#13#10;

for i := 0 to High(M) do

begin

texto := texto+IntToStr(i+1)+';';


texto := texto+FloatToStr(M[i, j])+';';


end;

Result := texto;

end;

//Determinação do valor máximo de cada coluna de uma matriz

function TForm1.MatMax(M: TMatriz): TVector;

var

vector: array of double;

i, j: integer;

begin

SetLength(vector, Length(M[0]));


begin

vector[j] := M[0, j];


if M[i, j]>vector[j] then

vector[j] := M[i, j];

end;

Result := vector;

end;

//Determinação do valor mínimo de cada coluna de uma matriz

function TForm1.MatMin(M: TMatriz): TVector;

var

vector: array of double;

i, j: integer;

begin



begin

vector[j] := M[0, j];


if M[i, j]<vector[j] then

vector[j] := M[i, j];

end;

Result := vector;

end;

//Determinação da posição do valor máximo de cada coluna de uma matriz

function TForm1.MatMaxPos(M: TMatriz; a, b: integer): TVector;

var

vector: TVector;

i, j: integer;

max: double;

begin



begin

A - 37

max := M[a, j];

vector[j] := a;

for i := a to b do

if (M[i, j]>max) and (abs(M[i, j]-max)>1e-4) then

begin

max := M[i, j];

vector[j] := i;

end;

end;

Result := vector;

end;

//Determinação da posição do valor mínimo de cada coluna de uma matriz

function TForm1.MatMinPos(M: TMatriz; a, b: integer): TVector;

var

vector: TVector;

i, j: integer;

min: double;

begin



begin

min := M[a, j];

vector[j] := a;

for i := a to b do

if (M[i, j]<min) and (abs(M[i, j]-min)>1e-4) then

begin

min := M[i, j];

vector[j] := i;

end;

end;

Result := vector;

end;

//Cópia dos dados de uma componente "StringGrid" (tabela gráfica) para uma matriz (Array)

function TForm1.StrGrdToArr(num: integer): TMatriz;

var

M: TMatriz;

i, j: integer;

begin

with Mdados[num] do

begin

SetLength(M, RowCount-1, ColCount-1);

for j := 0 to ColCount-2 do

for i := 0 to RowCount-2 do

M[i, j] := StrToFloat(Cells[j+1, i+1]);

end;

Result := M;

end;

//Tratamento do tipo de apresentação numérica

function TForm1.NumOut(n: double; e: integer): string;

var

txt: string;

begin

txt := FloatToStrF(n, ffFixed, 16, 3);

if abs(n)>=intpower(10, e) then

txt := FloatToStrF(n, ffExponent, 4, 1);

if abs(n)<intpower(10, -3) then

txt := FloatToStrF(n, ffExponent, 4, 1);

if n = 0 then

txt := '0';

Result := txt;

end;

//Desenho de uma seta

procedure TForm1.DSeta(ix, iy: double);

var

coord: TMatriz;

i: integer;

A - 38

begin

SetLength(coord, 4, 3);

coord[0, 0] := 25;

coord[2, 0] := coord[0, 0]-7;

coord[2, 1] := 4;

coord[3, 0] := coord[2, 0];

coord[3, 1] := -coord[2, 1];

coord := MultMat(coord, MatTr(1));

for i := 0 to 3 do

Image1.Canvas.Line(Round(ix+coord[0, 0]), Round(iy-coord[0, 1]),

Round(ix+coord[i, 0]), Round(iy-coord[i, 1]));

end;

//Função de importação de um ficheiro de dados

function TForm1.Abrir(Filter: string; Form: TForm): string;

var

openDialog: TOpenDialog;

begin

openDialog := TOpenDialog.Create(Form);

openDialog.InitialDir := GetCurrentDir;

openDialog.Options := [ofFileMustExist];

openDialog.Filter := Filter;

openDialog.FilterIndex := 1;

if openDialog.Execute then

Result := openDialog.FileName

else

Result := '';

openDialog.Free;

end;

//Função de exportação de ficheiros de dados

procedure TForm1.Guardar(Filter, Ext, texto: string; Form: TForm);

var

myFile: TextFile;

saveDialog: TSaveDialog;

begin

saveDialog := TSaveDialog.Create(Form);

saveDialog.Title := 'Guardar como';

saveDialog.InitialDir := GetCurrentDir;

saveDialog.Filter := Filter;

saveDialog.DefaultExt := Ext;

saveDialog.FilterIndex := 1;

if saveDialog.Execute then

begin

if FileExists(SaveDialog.FileName) then

if MessageDlg(Format('O ficheiro "%s" já existe. Pretende substituir?',

[ExtractFileName(SaveDialog.FileName)]), mtConfirmation, mbYesNo, 0) = idNo then

Exit;

AssignFile(myFile, saveDialog.FileName);

ReWrite(myFile);

Write(myFile, texto);

CloseFile(myFile);

end;

saveDialog.Free;

end;

//Selecciona as definições padrão do separador "Opções"


var

i: integer;

begin

for i := 0 to 6 do

CheckGroup1.Checked[i] := False;

A - 39


end;

//Selecciona as do separador "Opções"


var

appINI: TIniFile;

i: integer;

begin

appINI := TIniFile.Create(ChangeFileExt(Application.ExeName, '.ini'));

for i := 0 to 6 do

appINI.WriteBool('Exibir', CheckGroup1.Items[i], CheckGroup1.Checked[i]);

appINI.WriteBool('Visualização', RadioButton1.Caption, RadioButton1.Checked);

appINI.Free;

end;

//Funções complementares do funcionamento do programa

procedure TForm1.FormResize(Sender: TObject);

begin

if Form1.Width>527 then

begin

PageControl1.Width := Form1.Width+2;

Image1.Width := Form1.Width-28;

end;

if Form1.Height>595 then

begin

PageControl1.Height := Form1.Height-18;

Image1.Height := Form1.Height-96;

end;

end;

procedure TForm1.TabSheet0Show(Sender: TObject);

begin

if WindowState<>wsMaximized then

Form1.Height := 436;

end;

procedure TForm1.TabSheet8Show(Sender: TObject);

begin

if Form1.Height<596 then

Form1.Height := 596;

end;

procedure TForm1.RadioButton5Change(Sender: TObject);

begin

Edit10.Enabled := False;

end;


begin

Edit10.Enabled := True;

end;


begin

Edit11.Enabled := False;

end;


begin

Edit11.Enabled := True;

end;

procedure TForm1.StringGrid1ColRowInserted(Sender: TObject; IsColumn: boolean;

sIndex, tIndex: integer);

var

i, j: integer;

begin

with (Sender as TStringGrid) do

A - 40

for i := 1 to ColCount-1 do

for j := 1 to RowCount-1 do

if Cells[i, j] = '' then

if i = 10 then

Cells[i, j] := '1'

else

Cells[i, j] := '0';

end;

procedure TForm1.Edit0KeyPress(Sender: TObject; var Key: char);

begin

if not (Key in ['0'..'9', #8]) then

Key := #0;

end;

procedure TForm1.Edit0Change(Sender: TObject);

var

i, b: integer;

begin

if Nel[1].Text<>'' then

for i := 1 to 2 do

begin

Mdados[i].RowCount := StrToInt(Nel[1].Text)+1;

for b := 1 to StrToInt(Nel[1].Text) do

Mdados[i].Cells[0, b] := IntToStr(b);

end;

end;


var

i, b, k: integer;

begin

k := 0;

for i := 2 to 4 do

if Nel[i].Text<>'' then

begin

Mdados[i+1].RowCount := StrToInt(Nel[i].Text)+1;

for b := 1 to StrToInt(Nel[i].Text) do

begin

k := k+1;

Mdados[i+1].Cells[0, b] := IntToStr(k);

end;

end;

end;


var

i, b: integer;

begin

i := StrToInt(Copy((Sender as TEdit).Name, 5, 1));

if Nel[i].Text<>'' then

begin

Mdados[i+1].RowCount := StrToInt(Nel[i].Text)+1;

for b := 1 to StrToInt(Nel[i].Text) do

Mdados[i+1].Cells[0, b] := IntToStr(b);

end;

end;

procedure TForm1.FormCreate(Sender: TObject);

var

i: integer;

co: TComponent;

appINI: TIniFile;

begin

appINI := TIniFile.Create(ChangeFileExt(Application.ExeName, '.ini'));

for i := 0 to 6 do

CheckGroup1.Checked[i] := appINI.ReadBool('Exibir', CheckGroup1.Items[i], False);

RadioButton1.Checked := appINI.ReadBool('Visualização', RadioButton1.Caption, True);

appINI.Free;

if RadioButton1.Checked = False then


A - 41

for i := 0 to ComponentCount-1 do

begin

co := Components[i];

if co is TControl then

begin

TControl(co).Font.Name := 'Tahoma';

TControl(co).Font.Size := 8;

end;

end;

for i := 1 to High(Nel) do

Nel[i] := FindComponent(Format('Edit%d', [i])) as TEdit;

for i := 1 to High(Mdados) do

Mdados[i] := FindComponent(Format('StringGrid%d', [i])) as TStringGrid;

if DirectoryExists('Output') = False then

CreateDir('Output');

end;

initialization

{$I Unit1.lrs}

end.

A - 42

A.2.2 Unidade “CAxisSim”

//Matrizes resultantes da formulaçăo das soluçőes analíticas aplicada a

//estruturas axissimétricas

unit CAxisSim;

{$mode objfpc}{$H+}

interface

uses

Classes, SysUtils, Math, OpMat;

type

Vector=array of double;

Matriz=array of Vector;

procedure VecY(x: double);

function CilFe(): Matriz;

function CilKe(): Matriz;

function Cilqp(): Matriz;

function ResCil(): Matriz;

function AnelFe(): Matriz;

function AnelKe(): Matriz;

function MatInvC(): Matriz;

function ResAnel(x: double): Matriz;

var

Y: Array[1..4] of double;

implementation

uses

Unit1;

//Definiçăo das funçőes Y

procedure VecY(x: double);

begin

Y[1]:=Cos(B*x)*Cosh(B*x);

Y[2]:=(Cosh(B*x)*Sin(B*x) + Cos(B*x)*Sinh(B*x))/2.;

Y[3]:=(Sin(B*x)*Sinh(B*x))/2.;

Y[4]:=(Cosh(B*x)*Sin(B*x) - Cos(B*x)*Sinh(B*x))/4.;

end;

//Matriz de rigidez elementar da casca cilíndrica

function CilKe(): Matriz;

var

Kecil: Matriz;

begin

SetLength(KeCil,6,6);

Kecil[0,0]:=(8*B*E*h*Pi*r1*(Power (Y[3],2) - Y[2]*Y[4]))/(4*Power (v,2)*Power (Y[2],2)*Y[4] +

Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] + Power (v,2)*Power (Y[1],2) - 4*Power (v,2)*Power

(Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power (Y[3],2) + 2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power

(v,2)*Power (Y[4],3)));

Kecil[0,1]:=(-2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[2] +

4*Y[3]*Y[4]))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) -

2*Power (v,2)*Y[1] + Power (v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) +

4*(B*L*Power (Y[3],2) + 2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

A - 43

Kecil[0,2]:=(-2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[3] + 4*Power

(Y[4],2)))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) -



Kecil[0,3]:=(-8*B*E*h*Pi*r1*(Power (Y[3],2) - Y[2]*Y[4]))/(4*Power (v,2)*Power (Y[2],2)*Y[4] +



(v,2)*Power (Y[4],3)));

Kecil[0,4]:=(2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*(Power (Y[1],2)*Y[2] - Y[1]*(Y[2] -

8*Y[3]*Y[4]) + 4*(-(Y[2]*Power (Y[3],2)) + Power (Y[2],2)*Y[4] - Y[3]*Y[4] + 4*Power

(Y[4],3))))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) -



Kecil[0,5]:=(-2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Power (Y[2],2) -

Power (Y[1],2)*Y[3] - 4*Power (Y[3],3) + 8*Y[2]*Y[3]*Y[4] + Y[1]*(Y[3] - 4*Power




Kecil[1,0]:=(-2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[2] +




Kecil[1,1]:=(2*Power (B,3)*E*Power (h,3)*Pi*r1*(-(B*L*Y[1]*Y[2]) + Power (v,2)*Power (Y[2],2)

+ 4*Y[4]*(-(B*L*Y[3]) + Power (v,2)*Y[4])))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power

(Y[2],2)*Y[4] + Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] + Power (v,2)*Power (Y[1],2) - 4*Power

(v,2)*Power (Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power (Y[3],2) + 2*Power (v,2)*(-1 +

Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

Kecil[1,2]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*(-(B*L*Power (Y[2],2)) + Power (v,2)*Power

(Y[2],3) + Power (v,2)*Power (Y[1],2)*Y[4] - 4*Power (v,2)*Power (Y[3],2)*Y[4] +

Y[1]*(B*L*Y[3] - Power (v,2)*Y[4]) + Power (v,2)*Y[2]*((1 - 2*Y[1])*Y[3] + 4*Power




Kecil[1,3]:=(2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[2] +




Kecil[1,4]:=(2*Power (B,3)*E*Power (h,3)*Pi*r1*(B*L*Y[2] - Power (v,2)*Power (Y[2],2) -

4*Power (v,2)*Power (Y[4],2)))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] +



(v,2)*Power (Y[4],3))));

Kecil[1,5]:=(-2*Power (B,2)*E*Power (h,3)*Pi*r1*((B*L - Power (v,2)*Y[2])*Y[3] + Power

(v,2)*(-1 + Y[1])*Y[4]))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] +



(v,2)*Power (Y[4],3))));

Kecil[2,0]:=(-2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[3] + 4*Power




Kecil[2,1]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*(Power (v,2)*Y[2]*Y[3] + Y[4]*(Power (v,2) -

4*B*L*Y[4]) - Y[1]*(B*L*Y[3] + Power (v,2)*Y[4])))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Kecil[2,2]:=(B*E*Power (h,3)*Pi*r1*(2*Power (v,2)*Power (Y[1],2) - Power (v,2)*Power (Y[1],3)

- Y[1]*(Power (v,2) + 4*Power (v,2)*Power (Y[3],2) + (-4*B*L + 8*Power (v,2)*Y[2])*Y[4]) +

4*(Power (v,2)*Power (Y[2],2)*Y[3] + Y[2]*(-(B*L*Y[3]) + Power (v,2)*Y[4]) + Power

(v,2)*Y[3]*(Y[3] - 4*Power (Y[4],2)))))/(6.*(-1 + Power (v,2))*(4*Power (v,2)*Power

A - 44




Kecil[2,3]:=(2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[3] + 4*Power




Kecil[2,4]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*((B*L - Power (v,2)*Y[2])*Y[3] + Power (v,2)*(-

1 + Y[1])*Y[4]))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2)

- 2*Power (v,2)*Y[1] + Power (v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4])

+ 4*(B*L*Power (Y[3],2) + 2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power

(Y[4],3))));

Kecil[2,5]:=(B*E*Power (h,3)*Pi*r1*(Power (v,2) - 2*Power (v,2)*Y[1] + Power (v,2)*Power

(Y[1],2) + (-4*B*L + 4*Power (v,2)*Y[2])*Y[4]))/(6.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Kecil[3,0]:=(-8*B*E*h*Pi*r1*(Power (Y[3],2) - Y[2]*Y[4]))/(4*Power (v,2)*Power (Y[2],2)*Y[4] +



(v,2)*Power (Y[4],3)));

Kecil[3,1]:=(2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[2] +




Kecil[3,2]:=(2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Y[3] + 4*Power




Kecil[3,3]:=(8*B*E*h*Pi*r1*(Power (Y[3],2) - Y[2]*Y[4]))/(4*Power (v,2)*Power (Y[2],2)*Y[4] +



(v,2)*Power (Y[4],3)));

Kecil[3,4]:=(-2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*(Power (Y[1],2)*Y[2] - Y[1]*(Y[2]

- 8*Y[3]*Y[4]) + 4*(-(Y[2]*Power (Y[3],2)) + Power (Y[2],2)*Y[4] - Y[3]*Y[4] + 4*Power




Kecil[3,5]:=(2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Power (Y[2],2) - Power

(Y[1],2)*Y[3] - 4*Power (Y[3],3) + 8*Y[2]*Y[3]*Y[4] + Y[1]*(Y[3] - 4*Power (Y[4],2))))/(3.*(-1

+ Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] +

Power (v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power

(Y[3],2) + 2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

Kecil[4,0]:=(2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*(Power (Y[1],2)*Y[2] - Y[1]*(Y[2] -

8*Y[3]*Y[4]) + 4*(-(Y[2]*Power (Y[3],2)) + Power (Y[2],2)*Y[4] - Y[3]*Y[4] + 4*Power




Kecil[4,1]:=(2*Power (B,3)*E*Power (h,3)*Pi*r1*(B*L*Power (Y[1],2)*Y[2] - 4*B*L*Y[2]*Power

(Y[3],2) + 8*B*L*Y[1]*Y[3]*Y[4] - Power (Y[2],2)*(Power (v,2) - 4*B*L*Y[4]) + 4*Power

(Y[4],2)*(-Power(v,2) + 4*B*L*Y[4])))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Kecil[4,2]:=(-2*Power (B,2)*E*Power (h,3)*Pi*r1*(Power (v,2)*Power (Y[2],3) + 4*Power

(Y[3],2)*(B*L*Y[3] - Power (v,2)*Y[4]) + Power (Y[1],2)*(B*L*Y[3] + Power (v,2)*Y[4]) +

Y[2]*(4*Power (v,2)*Power (Y[4],2) + Y[3]*(Power (v,2) - 8*B*L*Y[4])) - Y[1]*(B*L*Power

(Y[2],2) + 2*Power (v,2)*Y[2]*Y[3] + Y[4]*(Power (v,2) - 4*B*L*Y[4]))))/(3.*(-1 + Power

(v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] + Power

A - 45

(v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power (Y[3],2) +

2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

Kecil[4,3]:=(-2*Power (B,4)*E*Power (h,3)*Pi*Power (r1,2)*v*(Power (Y[1],2)*Y[2] - Y[1]*(Y[2]

- 8*Y[3]*Y[4]) + 4*(-(Y[2]*Power (Y[3],2)) + Power (Y[2],2)*Y[4] - Y[3]*Y[4] + 4*Power




Kecil[4,4]:=(2*Power (B,3)*E*Power (h,3)*Pi*r1*(-(B*L*Y[1]*Y[2]) + Power (v,2)*Power (Y[2],2)

+ 4*Y[4]*(-(B*L*Y[3]) + Power (v,2)*Y[4])))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Kecil[4,5]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*(-(Power (v,2)*Y[2]*Y[3]) + Y[4]*(-Power(v,2) +

4*B*L*Y[4]) + Y[1]*(B*L*Y[3] + Power (v,2)*Y[4])))/(3.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Kecil[5,0]:=(-2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Power (Y[2],2) -

Power (Y[1],2)*Y[3] - 4*Power (Y[3],3) + 8*Y[2]*Y[3]*Y[4] + Y[1]*(Y[3] - 4*Power




Kecil[5,1]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*(Power (v,2)*Power (Y[2],3) + 4*Power

(Y[3],2)*(B*L*Y[3] - Power (v,2)*Y[4]) + Power (Y[1],2)*(B*L*Y[3] + Power (v,2)*Y[4]) +

Y[2]*(4*Power (v,2)*Power (Y[4],2) + Y[3]*(Power (v,2) - 8*B*L*Y[4])) - Y[1]*(B*L*Power

(Y[2],2) + 2*Power (v,2)*Y[2]*Y[3] + Y[4]*(Power (v,2) - 4*B*L*Y[4]))))/(3.*(-1 + Power

(v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] + Power

(v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power (Y[3],2) +

2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

Kecil[5,2]:=(B*E*Power (h,3)*Pi*r1*(-2*Power (v,2)*Power (Y[1],3) + Power (v,2)*Power (Y[1],4)

+ Power (Y[1],2)*(Power (v,2) + 8*Power (v,2)*Power (Y[3],2) - 4*(B*L - 4*Power

(v,2)*Y[2])*Y[4]) + 8*Y[1]*(-2*Power (v,2)*Power (Y[2],2)*Y[3] + Y[2]*(B*L*Y[3] - 2*Power

(v,2)*Y[4]) - Power (v,2)*Y[3]*(Y[3] - 8*Power (Y[4],2))) + 4*(-(B*L*Power (Y[2],3)) + Power

(v,2)*Power (Y[2],4) + Y[2]*Y[4]*(Power (v,2) - 16*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) +

2*Power (v,2)*Power (Y[2],2)*(Y[3] + 4*Power (Y[4],2)) + 4*(Power (v,2)*Power (Y[3],4) +

B*L*Power (Y[3],2)*Y[4] - 2*Power (v,2)*Y[3]*Power (Y[4],2) + 4*Power (v,2)*Power

(Y[4],4)))))/(6.*(-1 + Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) -



Kecil[5,3]:=(2*Power (B,3)*E*Power (h,3)*Pi*Power (r1,2)*v*((-1 + Y[1])*Power (Y[2],2) - Power

(Y[1],2)*Y[3] - 4*Power (Y[3],3) + 8*Y[2]*Y[3]*Y[4] + Y[1]*(Y[3] - 4*Power (Y[4],2))))/(3.*(-1

+ Power (v,2))*(4*Power (v,2)*Power (Y[2],2)*Y[4] + Y[2]*(Power (v,2) - 2*Power (v,2)*Y[1] +

Power (v,2)*Power (Y[1],2) - 4*Power (v,2)*Power (Y[3],2) - 4*B*L*Y[4]) + 4*(B*L*Power

(Y[3],2) + 2*Power (v,2)*(-1 + Y[1])*Y[3]*Y[4] + 4*Power (v,2)*Power (Y[4],3))));

Kecil[5,4]:=(2*Power (B,2)*E*Power (h,3)*Pi*r1*(B*L*Power (Y[2],2) - Power (v,2)*Power

(Y[2],3) - Power (v,2)*Power (Y[1],2)*Y[4] + 4*Power (v,2)*Power (Y[3],2)*Y[4] + Y[1]*(-

(B*L*Y[3]) + Power (v,2)*Y[4]) + Power (v,2)*Y[2]*((-1 + 2*Y[1])*Y[3] - 4*Power




Kecil[5,5]:=(B*E*Power (h,3)*Pi*r1*(2*Power (v,2)*Power (Y[1],2) - Power (v,2)*Power (Y[1],3)

- Y[1]*(Power (v,2) + 4*Power (v,2)*Power (Y[3],2) + (-4*B*L + 8*Power (v,2)*Y[2])*Y[4]) +

4*(Power (v,2)*Power (Y[2],2)*Y[3] + Y[2]*(-(B*L*Y[3]) + Power (v,2)*Y[4]) + Power

(v,2)*Y[3]*(Y[3] - 4*Power (Y[4],2)))))/(6.*(-1 + Power (v,2))*(4*Power (v,2)*Power




Result := Kecil;

end;

//Vector de forças nodais elementar da casca cilíndrica

A - 46

function CilFe(): Matriz;

var

Fecil: Matriz;

begin

SetLength(Fecil,6,1);

Fecil[0,0]:=L*Pi*Qx*r1;

Fecil[1,0]:=0;

Fecil[2,0]:=0;

Fecil[3,0]:=L*Pi*Qx*r1;

Fecil[4,0]:=0;

Fecil[5,0]:=0;

Fecil:=AddSubMat(Fecil,MultMat(CilKe(), Cilqp()),1);

Result := Fecil;

end;

//Vector de deslocamentos da soluçăo particular da casca cilíndrica

function Cilqp(): Matriz;

var

qp: Matriz;

begin

SetLength(qp,6,1);

qp[0,0]:=0;

qp[1,0]:=(r1*(2*Qy1*r1 + L*Qx*v))/(2.*E*h);

qp[2,0]:=-((r1*(Qy1*r1 - Qy2*r1 + L*Qx*v))/(E*h*L));

qp[3,0]:=(L*(Qy1 + Qy2)*r1*v)/(2.*E*h);

qp[4,0]:=(r1*(2*Qy2*r1 - L*Qx*v))/(2.*E*h);

qp[5,0]:=-((r1*(Qy1*r1 - Qy2*r1 + L*Qx*v))/(E*h*L));

Result := qp;

end;

//Vector de resultados parciais da casca cilíndrica

function ResCil(): Matriz;

var

RCil: Matriz;

begin

SetLength(RCil,5,4);

RCil[0,0]:=Y[2]/B;

RCil[0,1]:=Y[3]/Power(B,2);

RCil[0,2]:=(12*(1 - Power (v,2))*Y[4])/(Power (B,3)*E*Power (h,3));

RCil[0,3]:=(-3*(1 - Power (v,2))*(-1 + Y[1]))/(Power (B,4)*E*Power (h,3));

RCil[1,0]:=Y[1];

RCil[1,1]:=Y[2]/B;



RCil[2,0]:=-4*B*Y[4];

A - 47

RCil[2,1]:=Y[1];

RCil[2,2]:=(12*(1 - Power (v,2))*Y[2])/(B*E*Power (h,3));


RCil[3,0]:=(Power (B,3)*E*Power (h,3)*Y[2])/(3.*(1 - Power (v,2)));

RCil[3,1]:=(Power (B,2)*E*Power (h,3)*Y[3])/(3.*(1 - Power (v,2)));

RCil[3,2]:=4*B*Y[4];

RCil[3,3]:=-Y[1];

RCil[4,0]:=-(Power (B,2)*E*Power (h,3)*Y[3])/(3.*(1 - Power (v,2)));

RCil[4,1]:=-(B*E*Power (h,3)*Y[4])/(3.*(1 - Power (v,2)));

RCil[4,2]:=Y[1];

RCil[4,3]:=Y[2]/B;

Result := RCil;

end;

//Matriz de rigidez elementar da casca em anel

function AnelKe(): Matriz;

var

Keanel: Matriz;

begin

SetLength(Keanel,6,6);

if r1=0 then //Caso particular onde o nó inicial tem raio nulo constituindo uma casca

circular

begin

Keanel[0,0]:=(2*E*h*Pi)/(1 + v);

Keanel[0,1]:=0;

Keanel[0,2]:=0;

Keanel[0,3]:=0;

Keanel[0,4]:=0;

Keanel[0,5]:=0;

Keanel[1,0]:=0;

Keanel[1,1]:=(-4*E*Power (h,3)*Pi)/(3.*Power (r2,2)*(-1 + Power (v,2)));

Keanel[1,2]:=0;

Keanel[1,3]:=0;

Keanel[1,4]:=(4*E*Power (h,3)*Pi)/(3.*Power (r2,2)*(-1 + Power (v,2)));

Keanel[1,5]:=(2*E*Power (h,3)*Pi)/(3*r2 - 3*r2*Power (v,2));

Keanel[2,0]:=0;

Keanel[2,1]:=0;

Keanel[2,2]:=(E*Power (h,3)*Pi)/(6 + 6*v);

Keanel[2,3]:=0;

Keanel[2,4]:=0;

Keanel[2,5]:=0;

A - 48

Keanel[3,0]:=0;

Keanel[3,1]:=0;

Keanel[3,2]:=0;

Keanel[3,3]:=(-2*E*h*Pi)/(-1 + v);

Keanel[3,4]:=0;

Keanel[3,5]:=0;

Keanel[4,0]:=0;

Keanel[4,1]:=(4*E*Power (h,3)*Pi)/(3.*Power (r2,2)*(-1 + Power (v,2)));

Keanel[4,2]:=0;

Keanel[4,3]:=0;

Keanel[4,4]:=(-4*E*Power (h,3)*Pi)/(3.*Power (r2,2)*(-1 + Power (v,2)));

Keanel[4,5]:=(2*E*Power (h,3)*Pi)/(3.*r2*(-1 + Power (v,2)));

Keanel[5,0]:=0;

Keanel[5,1]:=(2*E*Power (h,3)*Pi)/(3*r2 - 3*r2*Power (v,2));

Keanel[5,2]:=0;

Keanel[5,3]:=0;

Keanel[5,4]:=(2*E*Power (h,3)*Pi)/(3.*r2*(-1 + Power (v,2)));

Keanel[5,5]:=-(E*Power (h,3)*Pi*(3 + v))/(6.*(-1 + Power (v,2)));

end

else

begin

Keanel[0,0]:=(2*E*h*Pi*((Power (r1,2) + Power (r2,2))/(-Power(r1,2) + Power (r2,2)) - v))/(1 -

Power (v,2));

Keanel[0,1]:=0;

Keanel[0,2]:=0;

Keanel[0,3]:=(-4*E*h*Pi*r1*r2)/((-Power(r1,2) + Power (r2,2))*(1 - Power (v,2)));

Keanel[0,4]:=0;

Keanel[0,5]:=0;

Keanel[1,0]:=0;

Keanel[1,1]:=(4*E*Power (h,3)*Pi*(Power (r1,2) - Power (r2,2)))/(3.*(-1 + Power (v,2))*(Power

(Power (r1,2) - Power (r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2)));

Keanel[1,2]:=(-2*E*Power (h,3)*Pi*r1*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln

(r1/r2)))/(3.*(-1 + Power (v,2))*(Power (Power (r1,2) - Power (r2,2),2) - 4*Power (r1,2)*Power

(r2,2)*Power (ln (r1/r2),2)));

Keanel[1,3]:=0;

Keanel[1,4]:=(-4*E*Power (h,3)*Pi*(Power (r1,2) - Power (r2,2)))/(3.*(-1 + Power (v,2))*(Power


Keanel[1,5]:=(-2*E*Power (h,3)*Pi*r2*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[2,0]:=0;

A - 49

Keanel[2,1]:=(-2*E*Power (h,3)*Pi*r1*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[2,2]:=(E*Power (h,3)*Pi*((Power (r1,2) - Power (r2,2))*(-(Power (r2,2)*(-1 + v)) +

Power (r1,2)*(3 + v)) - 8*Power (r1,2)*Power (r2,2)*ln (r1/r2) - 4*Power (r1,2)*Power

(r2,2)*(1 + v)*Power (ln (r1/r2),2)))/(6.*(-1 + Power (v,2))*(Power (Power (r1,2) - Power

(r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2)));

Keanel[2,3]:=0;

Keanel[2,4]:=(2*E*Power (h,3)*Pi*r1*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[2,5]:=(2*E*Power (h,3)*Pi*r1*r2*(-Power(r1,2) + Power (r2,2) + (Power (r1,2) + Power

(r2,2))*ln (r1/r2)))/(3.*(-1 + Power (v,2))*(Power (Power (r1,2) - Power (r2,2),2) - 4*Power

(r1,2)*Power (r2,2)*Power (ln (r1/r2),2)));

Keanel[3,0]:=(-4*E*h*Pi*r1*r2)/((-Power(r1,2) + Power (r2,2))*(1 - Power (v,2)));

Keanel[3,1]:=0;

Keanel[3,2]:=0;

Keanel[3,3]:=(2*E*h*Pi*((Power (r1,2) + Power (r2,2))/(-Power(r1,2) + Power (r2,2)) + v))/(1 -

Power (v,2));

Keanel[3,4]:=0;

Keanel[3,5]:=0;

Keanel[4,0]:=0;

Keanel[4,1]:=(4*E*Power (h,3)*Pi*(-Power(r1,2) + Power (r2,2)))/(3.*(-1 + Power (v,2))*(Power


Keanel[4,2]:=(2*E*Power (h,3)*Pi*r1*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[4,3]:=0;

Keanel[4,4]:=(4*E*Power (h,3)*Pi*(Power (r1,2) - Power (r2,2)))/(3.*(-1 + Power (v,2))*(Power


Keanel[4,5]:=(2*E*Power (h,3)*Pi*r2*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[5,0]:=0;

Keanel[5,1]:=(-2*E*Power (h,3)*Pi*r2*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[5,2]:=(2*E*Power (h,3)*Pi*r1*r2*(-Power(r1,2) + Power (r2,2) + (Power (r1,2) + Power

(r2,2))*ln (r1/r2)))/(3.*(-1 + Power (v,2))*(Power (Power (r1,2) - Power (r2,2),2) - 4*Power

(r1,2)*Power (r2,2)*Power (ln (r1/r2),2)));

Keanel[5,3]:=0;

Keanel[5,4]:=(2*E*Power (h,3)*Pi*r2*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln


(r2,2)*Power (ln (r1/r2),2)));

Keanel[5,5]:=-(E*Power (h,3)*Pi*((Power (r1,2) - Power (r2,2))*(Power (r1,2)*(-1 + v) - Power

(r2,2)*(3 + v)) + 8*Power (r1,2)*Power (r2,2)*ln (r1/r2) - 4*Power (r1,2)*Power (r2,2)*(1 +

v)*Power (ln (r1/r2),2)))/(6.*(-1 + Power (v,2))*(Power (Power (r1,2) - Power (r2,2),2) -

4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2)));

A - 50

end;

Result := Keanel;

end;

//Vector de forças nodais elementar da casca em anel

function AnelFe(): Matriz;

var

Feanel, qp, temp: Matriz;

begin

SetLength(Feanel,6,1);

SetLength(qp,6,1);

Feanel[0,0]:=0;

Feanel[1,0]:=-(Pi*Qy1*Power (r1,2));

Feanel[2,0]:=(Pi*Qy1*Power (r1,3)*(3 + v))/8.;

Feanel[3,0]:=0;

Feanel[4,0]:=Pi*Qy1*Power (r2,2);

Feanel[5,0]:=(Pi*Qy1*Power (r2,3)*(-3 - v))/8.;

qp[0,0]:=0;

qp[1,0]:=(3*Qy1*Power (r1,4)*(1 - Power (v,2)))/(16.*E*Power (h,3));


qp[3,0]:=0;



temp:=MultMat(AnelKe(), qp);

Feanel:=AddSubMat(Feanel,temp,1);

Result := Feanel;

end;

//Inverso da Matriz C da casca em anel

function MatInvC(): Matriz;

var

InvC: Matriz;

begin

SetLength(InvC,4,4);

if r1=0 then //Caso particular da casca circular

begin

InvC[0,0]:=1;

InvC[0,1]:=0;

InvC[0,2]:=0;

InvC[0,3]:=0;

InvC[1,0]:=0;

InvC[1,1]:=0;

InvC[1,2]:=0;

InvC[1,3]:=0;

InvC[2,0]:=-1;

A - 51

InvC[2,1]:=0;

InvC[2,2]:=1;

InvC[2,3]:=0;

InvC[3,0]:=2;

InvC[3,1]:=0;

InvC[3,2]:=-2;

InvC[3,3]:=r2;

end

else

begin

InvC[0,0]:=(Power (r2,2)*(-Power(r1,2) + Power (r2,2) - 2*Power (r1,2)*ln (r1/r2)))/(Power

(Power (r1,2) - Power (r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[0,1]:=(r1*Power (r2,2)*(Power (r1,2) - Power (r2,2))*ln (r1/r2))/(Power (Power (r1,2) -

Power (r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[0,2]:=(Power (r1,2)*(Power (r1,2) - Power (r2,2) + 2*Power (r2,2)*ln (r1/r2) - 4*Power

(r2,2)*Power (ln (r1/r2),2)))/(Power (Power (r1,2) - Power (r2,2),2) - 4*Power (r1,2)*Power

(r2,2)*Power (ln (r1/r2),2));

InvC[0,3]:=(2*Power (r1,2)*Power (r2,3)*Power (ln (r1/r2),2))/(Power (Power (r1,2) - Power

(r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[1,0]:=(-4*Power (r1,2)*Power (r2,2)*ln (r1/r2))/(Power (Power (r1,2) - Power (r2,2),2) -

4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[1,1]:=(r1*Power (r2,2)*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln (r1/r2)))/(Power


InvC[1,2]:=(4*Power (r1,2)*Power (r2,2)*ln (r1/r2))/(Power (Power (r1,2) - Power (r2,2),2) -

4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[1,3]:=(Power (r1,2)*r2*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln (r1/r2)))/(Power


InvC[2,0]:=(Power (r2,2)*(Power (r1,2) - Power (r2,2) + 2*Power (r1,2)*ln (r1/r2)))/(Power


InvC[2,1]:=(r1*Power (r2,2)*(-Power(r1,2) + Power (r2,2))*ln (r1/r2))/(Power (Power (r1,2) -

Power (r2,2),2) - 4*Power (r1,2)*Power (r2,2)*Power (ln (r1/r2),2));

InvC[2,2]:=(Power (r2,2)*(-Power(r1,2) + Power (r2,2) - 2*Power (r1,2)*ln (r1/r2)))/(Power


InvC[2,3]:=(-2*Power (r1,2)*Power (r2,3)*Power (ln (r1/r2),2))/(Power (Power (r1,2) - Power


InvC[3,0]:=(2*Power (r2,2)*(-Power(r1,2) + Power (r2,2)))/(Power (Power (r1,2) - Power


InvC[3,1]:=(r1*Power (r2,2)*(Power (r1,2) - Power (r2,2) - 2*Power (r2,2)*ln (r1/r2)))/(Power


InvC[3,2]:=(2*Power (r2,2)*(Power (r1,2) - Power (r2,2)))/(Power (Power (r1,2) - Power


InvC[3,3]:=(Power (r2,3)*(-Power(r1,2) + Power (r2,2) + 2*Power (r1,2)*ln (r1/r2)))/(Power


end;

Result := InvC;

end;

//Vector de resultados parciais da casca em anel

A - 52

function ResAnel(x: double): Matriz;

var

RAnel: Matriz;

begin

SetLength(RAnel,5,4);

RAnel[0,0]:=1;

RAnel[0,1]:=Ln ((r1 + ((-r1 + r2)*x)/L)/r2);

RAnel[0,2]:=Power (r1 + ((-r1 + r2)*x)/L,2)/Power(r2,2);

RAnel[0,3]:=(Power (r1 + ((-r1 + r2)*x)/L,2)*Ln ((r1 + ((-r1 + r2)*x)/L)/r2))/Power(r2,2);

RAnel[1,0]:=0;

RAnel[1,1]:=L/(L*r1 - r1*x + r2*x);

RAnel[1,2]:=(2*(r1 + ((-r1 + r2)*x)/L))/Power(r2,2);

RAnel[1,3]:=((r1 + ((-r1 + r2)*x)/L)*(1 + 2*Ln ((r1 + ((-r1 + r2)*x)/L)/r2)))/Power(r2,2);

RAnel[2,0]:=0;

RAnel[2,1]:=-(E*Power (h,3)*Power (L,2))/(12.*(1 + v)*Power (L*r1 + (-r1 + r2)*x,2));

RAnel[2,2]:=-(E*Power (h,3))/(6.*Power (r2,2)*(-1 + v));

RAnel[2,3]:=-(E*Power (h,3)*(3 + v + 2*(1 + v)*Ln ((r1 + ((-r1 + r2)*x)/L)/r2)))/(12.*Power

(r2,2)*(-1 + Power (v,2)));

RAnel[3,0]:=0;

RAnel[3,1]:=(E*Power (h,3)*Power (L,2))/(12.*(1 + v)*Power (L*r1 + (-r1 + r2)*x,2));

RAnel[3,2]:=-(E*Power (h,3))/(6.*Power (r2,2)*(-1 + v));

RAnel[3,3]:=-(E*Power (h,3)*(1 + 3*v + 2*(1 + v)*Ln ((r1 + ((-r1 + r2)*x)/L)/r2)))/(12.*Power

(r2,2)*(-1 + Power (v,2)));

RAnel[4,0]:=0;

RAnel[4,1]:=0;

RAnel[4,2]:=0;

RAnel[4,3]:=(E*Power (h,3)*L)/(3.*Power (r2,2)*(-1 + Power (v,2))*(L*r1 + (-r1 + r2)*x));

Result := RAnel;

end;

end.

A - 53

A.2.3 Unidade “MEFAxisSim”

//Matrizes resultantes da formulação do método dos elementos finitos aplicada a

//estruturas axissimétricas

unit MEFAxisSim;

{$mode objfpc}{$H+}

interface

uses

Classes, SysUtils, Math;

type

Matriz=array of array of double;

function MatFe(): Matriz;

function MatPsi(x: double): Matriz;

function MatTr(s: integer): Matriz;

function MatVe(x: double): Matriz;

procedure MatBe(x: double);

procedure MatBeSing();

procedure MatD();

var

Be,Dg: Matriz;

implementation

uses

Unit1;

//Matriz B elementar

procedure MatBe(x: double);

begin

SetLength(Be,4,6);

Be[0,0]:=-(1/L);

Be[0,1]:=0;

Be[0,2]:=0;

Be[0,3]:=1/L;

Be[0,4]:=0;

Be[0,5]:=0;

Be[1,0]:=((1 - x/L)*ca)/(r1 + ((-r1 + r2)*x)/L);

Be[1,1]:=-(((1 - (3*Power (x,2))/Power(L,2) + (2*Power (x,3))/Power(L,3))*sa)/(r1 + ((-r1 +

r2)*x)/L));

Be[1,2]:=-(((x - (2*Power (x,2))/L + Power (x,3)/Power(L,2))*sa)/(r1 + ((-r1 + r2)*x)/L));

Be[1,3]:=(x*ca)/(L*(r1 + ((-r1 + r2)*x)/L));

Be[1,4]:=-((((3*Power (x,2))/Power(L,2) - (2*Power (x,3))/Power(L,3))*sa)/(r1 + ((-r1 +

r2)*x)/L));

Be[1,5]:=-(((-(Power (x,2)/L) + Power (x,3)/Power(L,2))*sa)/(r1 + ((-r1 + r2)*x)/L));

Be[2,0]:=0;

Be[2,1]:=-6/Power(L,2) + (12*x)/Power(L,3);

Be[2,2]:=-4/L + (6*x)/Power(L,2);

Be[2,3]:=0;

Be[2,4]:=6/Power(L,2) - (12*x)/Power(L,3);

Be[2,5]:=-2/L + (6*x)/Power(L,2);

Be[3,0]:=0;

Be[3,1]:=(((-6*x)/Power(L,2) + (6*Power (x,2))/Power(L,3))*ca)/(r1 + ((-r1 + r2)*x)/L);

Be[3,2]:=((1 - (4*x)/L + (3*Power (x,2))/Power(L,2))*ca)/(r1 + ((-r1 + r2)*x)/L);

Be[3,3]:=0;

Be[3,4]:=(((6*x)/Power(L,2) - (6*Power (x,2))/Power(L,3))*ca)/(r1 + ((-r1 + r2)*x)/L);

Be[3,5]:=(((-2*x)/L + (3*Power (x,2))/Power(L,2))*ca)/(r1 + ((-r1 + r2)*x)/L);

end;

//Valores limite da Matriz B elementar para singularidades

procedure MatBeSing();

begin

A - 54

SetLength(Be,4,6);

if r1=0 then

begin

Be[0,0]:=0;

Be[1,0]:=0;

Be[2,0]:=0;

Be[3,0]:=0;

Be[0,1]:=-(sa/L);

Be[1,1]:=-((ca*sa)/r2);

Be[2,1]:=(-6*ca)/Power(L,2);

Be[3,1]:=(-6*Power (ca,2))/(L*r2);

Be[0,2]:=0;

Be[1,2]:=0;

Be[2,2]:=0;

Be[3,2]:=0;

Be[0,3]:=ca/L;

Be[1,3]:=Power (ca,2)/r2;

Be[2,3]:=(-6*sa)/Power(L,2);

Be[3,3]:=(-6*ca*sa)/(L*r2);

Be[0,4]:=sa/L;

Be[1,4]:=(ca*sa)/r2;

Be[2,4]:=(6*ca)/Power(L,2);

Be[3,4]:=(6*Power (ca,2))/(L*r2);

Be[0,5]:=0;

Be[1,5]:=0;

Be[2,5]:=-2/L;

Be[3,5]:=(-2*ca)/r2;

end;

if r2=0 then

begin

Be[0,0]:=-(ca/L);

Be[1,0]:=Power (ca,2)/r1;

Be[2,0]:=(-6*sa)/Power(L,2);

Be[3,0]:=(6*ca*sa)/(L*r1);

Be[0,1]:=-(sa/L);

Be[1,1]:=(ca*sa)/r1;

Be[2,1]:=(6*ca)/Power(L,2);

Be[3,1]:=(-6*Power (ca,2))/(L*r1);

Be[0,2]:=0;

Be[1,2]:=0;

Be[2,2]:=2/L;

Be[3,2]:=(-2*ca)/r1;

Be[0,3]:=0;

Be[1,3]:=0;

Be[2,3]:=0;

Be[3,3]:=0;

Be[0,4]:=sa/L;

Be[1,4]:=-((ca*sa)/r1);

Be[2,4]:=(-6*ca)/Power(L,2);

Be[3,4]:=(6*Power (ca,2))/(L*r1);

Be[0,5]:=0;

Be[1,5]:=0;

Be[2,5]:=0;

Be[3,5]:=0;

end;

end;

//Matriz das funções de forma elementar

function MatPsi(x: double): Matriz;

var

Psi: Matriz;

begin

SetLength(Psi,3,6);

Psi[0,0]:=1 - x/L;

Psi[0,1]:=0;

Psi[0,2]:=0;

A - 55

Psi[0,3]:=x/L;

Psi[0,4]:=0;

Psi[0,5]:=0;

Psi[1,0]:=0;

Psi[1,1]:=1 - (3*Power(x,2))/Power(L,2) + (2*Power(x,3))/Power(L,3);

Psi[1,2]:=x - (2*Power(x,2))/L + Power(x,3)/Power(L,2);

Psi[1,3]:=0;

Psi[1,4]:=(3*Power(x,2))/Power(L,2) - (2*Power(x,3))/Power(L,3);

Psi[1,5]:=-(Power(x,2)/L) + Power(x,3)/Power(L,2);

Psi[2,0]:=0;

Psi[2,1]:=(-6*x)/Power(L,2) + (6*Power(x,2))/Power(L,3);

Psi[2,2]:=1 - (4*x)/L + (3*Power(x,2))/Power(L,2);

Psi[2,3]:=0;

Psi[2,4]:=(6*x)/Power(L,2) - (6*Power(x,2))/Power(L,3);

Psi[2,5]:=(-2*x)/L + (3*Power(x,2))/Power(L,2);

Result := Psi

end;

//Vector de forças nodais equivalentes elementar

function MatFe(): Matriz;

var

Fe: Matriz;

begin

SetLength(Fe,6,1);

Fe[0,0]:=(L*Pi*Qx*(2*r1 + r2))/3.;

Fe[1,0]:=(L*Pi*(16*Qy1*r1 + 5*Qy2*r1 + 5*Qy1*r2 + 4*Qy2*r2))/30.;

Fe[2,0]:=(Power(L,2)*Pi*(Qy2*(r1 + r2) + Qy1*(2*r1 + r2)))/30.;

Fe[3,0]:=(L*Pi*Qx*(r1 + 2*r2))/3.;

Fe[4,0]:=(L*Pi*(4*Qy1*r1 + 5*Qy2*r1 + 5*Qy1*r2 + 16*Qy2*r2))/30.;

Fe[5,0]:=-(Power(L,2)*Pi*(Qy1*(r1 + r2) + Qy2*(r1 + 2*r2)))/30.;

Result := Fe;

end;

//Matriz de transformação de coordenadas

function MatTr(s: integer): Matriz;

var

Tr: Matriz;

i: integer;

begin

SetLength(Tr,3*s,3*s);

for i:=0 to s-1 do

begin

Tr[0+3*i,0+3*i]:=ca;

Tr[0+3*i,1+3*i]:=sa;

Tr[1+3*i,0+3*i]:=-sa;

Tr[1+3*i,1+3*i]:=ca;

Tr[2+3*i,2+3*i]:=1;

end;

A - 56

Result := Tr

end;

//Matriz das relações constitutivas

procedure MatD();

begin

SetLength(Dg,4,4);

Dg[0,0]:=(E*h)/(1 - Power (v,2));

Dg[0,1]:=(E*h*v)/(1 - Power (v,2));

Dg[0,2]:=0;

Dg[0,3]:=0;

Dg[1,0]:=(E*h*v)/(1 - Power (v,2));

Dg[1,1]:=(E*h)/(1 - Power (v,2));

Dg[1,2]:=0;

Dg[1,3]:=0;

Dg[2,0]:=0;

Dg[2,1]:=0;

Dg[2,2]:=(E*Power (h,3))/(12.*(1 - Power (v,2)));

Dg[2,3]:=(E*Power (h,3)*v)/(12.*(1 - Power (v,2)));

Dg[3,0]:=0;

Dg[3,1]:=0;

Dg[3,2]:=(E*Power (h,3)*v)/(12.*(1 - Power (v,2)));

Dg[3,3]:=(E*Power (h,3))/(12.*(1 - Power (v,2)));

end;

//Vector de esforço transverso, obtido por equilíbrio

function MatVe(x: double): Matriz;

var

Ve: Matriz;

begin

SetLength(Ve,1,6);

if (r1=0) and (x=0) then x:=1e-15;

if (r2=0) and (x=L) then x:=L-1e-15;

Ve[0,0]:=0;

Ve[0,1]:=(-((E*h*(r1 + ((-r1 + r2)*x)/L)*(Power (h,2)/Power(L,3) + (ca*Power (h,2)*v*(-

6/Power(L,2) + (12*x)/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L)) - (ca*Power (h,2)*(-r1 +

r2)*v*((-6*x)/Power(L,2) + (6*Power (x,2))/Power(L,3)))/(12.*L*Power (r1 + ((-r1 +

r2)*x)/L,2))))/(1 - Power (v,2))) + (ca*E*h*((Power (h,2)*v*(-6/Power(L,2) +

(12*x)/Power(L,3)))/12. + (ca*Power (h,2)*((-6*x)/Power(L,2) + (6*Power

(x,2))/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(1 - Power (v,2)) - (E*h*(-r1 +

r2)*((Power (h,2)*(-6/Power(L,2) + (12*x)/Power(L,3)))/12. + (ca*Power (h,2)*v*((-

6*x)/Power(L,2) + (6*Power (x,2))/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(L*(1 - Power

(v,2))))/(r1 + ((-r1 + r2)*x)/L);

Ve[0,2]:=(-((E*h*(r1 + ((-r1 + r2)*x)/L)*(Power (h,2)/(2.*Power (L,2)) + (ca*Power (h,2)*v*(-

4/L + (6*x)/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L)) - (ca*Power (h,2)*(-r1 + r2)*v*(1 -

(4*x)/L + (3*Power (x,2))/Power(L,2)))/(12.*L*Power (r1 + ((-r1 + r2)*x)/L,2))))/(1 - Power

(v,2))) + (ca*E*h*((Power (h,2)*v*(-4/L + (6*x)/Power(L,2)))/12. + (ca*Power (h,2)*(1 -

(4*x)/L + (3*Power (x,2))/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(1 - Power (v,2)) -

(E*h*(-r1 + r2)*((Power (h,2)*(-4/L + (6*x)/Power(L,2)))/12. + (ca*Power (h,2)*v*(1 - (4*x)/L

A - 57

+ (3*Power (x,2))/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(L*(1 - Power (v,2))))/(r1 +

((-r1 + r2)*x)/L);

Ve[0,3]:=0;

Ve[0,4]:=(-((E*h*(r1 + ((-r1 + r2)*x)/L)*(-(Power (h,2)/Power(L,3)) + (ca*Power

(h,2)*v*(6/Power(L,2) - (12*x)/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L)) - (ca*Power (h,2)*(-

r1 + r2)*v*((6*x)/Power(L,2) - (6*Power (x,2))/Power(L,3)))/(12.*L*Power (r1 + ((-r1 +

r2)*x)/L,2))))/(1 - Power (v,2))) + (ca*E*h*((Power (h,2)*v*(6/Power(L,2) -

(12*x)/Power(L,3)))/12. + (ca*Power (h,2)*((6*x)/Power(L,2) - (6*Power

(x,2))/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(1 - Power (v,2)) - (E*h*(-r1 +

r2)*((Power (h,2)*(6/Power(L,2) - (12*x)/Power(L,3)))/12. + (ca*Power

(h,2)*v*((6*x)/Power(L,2) - (6*Power (x,2))/Power(L,3)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(L*(1

- Power (v,2))))/(r1 + ((-r1 + r2)*x)/L);

Ve[0,5]:=(-((E*h*(r1 + ((-r1 + r2)*x)/L)*(Power (h,2)/(2.*Power (L,2)) + (ca*Power (h,2)*v*(-

2/L + (6*x)/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L)) - (ca*Power (h,2)*(-r1 + r2)*v*((-

2*x)/L + (3*Power (x,2))/Power(L,2)))/(12.*L*Power (r1 + ((-r1 + r2)*x)/L,2))))/(1 - Power

(v,2))) + (ca*E*h*((Power (h,2)*v*(-2/L + (6*x)/Power(L,2)))/12. + (ca*Power (h,2)*((-2*x)/L +

(3*Power (x,2))/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(1 - Power (v,2)) - (E*h*(-r1 +

r2)*((Power (h,2)*(-2/L + (6*x)/Power(L,2)))/12. + (ca*Power (h,2)*v*((-2*x)/L + (3*Power

(x,2))/Power(L,2)))/(12.*(r1 + ((-r1 + r2)*x)/L))))/(L*(1 - Power (v,2))))/(r1 + ((-r1 +

r2)*x)/L);

Result := Ve

end;

end.

A - 58

A.2.4 Unidade “Equacoes”

//Unidade de rotinas concebida por: Professor Manuel Ritto Correa

//Resolução de sistemas de equações lineares pelo método de eliminação de Gauss.

unit Equacoes;

interface

Type

Linha=array of double;

Coluna=array of double;

Matriz=array of Linha;

{ Equacao }

{ CEquacao }

CEquacao=Class

private

N,NR:Integer;

bt:Matriz;

xt:Matriz;

Preserve:Boolean;

procedure CriaMatriz(Var m:Matriz;nl,nc:Integer);

public

A:Matriz; //TODO: passar a private

constructor Cria(Neq,Nrhs:Integer;QualPreserve:Boolean);

procedure ZeraMatriz;

procedure AtribuiAij(i,j:Integer;valor:double);

procedure AtribuiBij(i,j: Integer; valor: double);

procedure AtribuiBj(j: Integer; valor: double);

procedure AtribuiColunaMatriz(j:Integer;bj:Coluna);

procedure AtribuiVectorIndependente(j:Integer;bj:Coluna);

function Solucao(j:Integer):Coluna;

destructor done;

procedure ocalhas;

procedure Simetriza;

procedure Solve;

procedure Verifica;

procedure mostra(var Flog:text);

end;

Procedure produtor(r:Double;m:Matriz;n1,n2:Integer;Var mr:Matriz);

Procedure produtot(m1,m2t:Matriz;n1,nc,n2:Integer;Var mr:Matriz);

Procedure produto(m1,m2:Matriz;n1,nc,n2:Integer;Var mr:Matriz);

Procedure soma(m1,m2:Matriz;n1,n2:Integer;Var mr:Matriz);

implementation

Constructor CEquacao.Cria(Neq,Nrhs:Integer;QualPreserve:Boolean);

begin

self.N:=Neq;

self.NR:=Nrhs;

self.Preserve:=QualPreserve;

self.CriaMatriz(self.A,self.N,self.N);

self.CriaMatriz(self.bt,self.NR,self.N);

self.CriaMatriz(self.xt,self.NR,self.N);

//self.ocalhas;

end;

procedure CEquacao.ZeraMatriz;

var j1,j2: integer;

begin

for j1:=1 to N do

for j2:=1 to N do

A[j1,j2]:=0;

end;

procedure CEquacao.AtribuiAij(i, j: Integer; valor: double);

begin

A - 59

A[i,j]:=valor;

end;

procedure CEquacao.AtribuiBij(i, j: Integer; valor: double);

begin

bt[i,j]:=valor;

end;

procedure CEquacao.AtribuiBj(j: Integer; valor: double);

begin

bt[1,j]:=valor;

end;

procedure CEquacao.AtribuiColunaMatriz(j: Integer; bj: Coluna);

var i: integer;

begin

for i:=1 to N do A[i,j]:=bj[i];

end;

procedure CEquacao.AtribuiVectorIndependente(j: Integer; bj: Coluna);

var i: integer;

begin

for i:=1 to N do bt[j,i]:=bj[i];

end;

function CEquacao.Solucao(j: Integer): Coluna;

var i: integer;

begin

Setlength(solucao,N+1);

for i:=1 to N do solucao[i]:=xt[j,i];

end;

Destructor CEquacao.done;

begin

end;

Procedure CEquacao.CriaMatriz(Var m:Matriz;nl,nc:Integer);

var i:Integer;

begin

Setlength(m,nl+1);

For i:=1 to nl do Setlength(m[i],nc+1);

end;

Procedure CEquacao.ocalhas;

var i,j:Integer;

begin

For i:=1 to self.N do

begin

For j:=1 to self.N do self.A[i,j]:=random;

For j:=1 to self.NR do self.bt[j,i]:=random;

For j:=1 to self.NR do self.xt[j,i]:=0;

end;

end;

Procedure CEquacao.mostra(var Flog:text);

var i,j:Integer;

begin


begin

For j:=1 to self.N do Write(FLog,self.A[i,j]:10:3);

Write(FLog,' *');

For j:=1 to self.NR do Write(FLog,self.xt[j,i]:13:6);

Write(FLog,' =');

For j:=1 to self.NR do Write(FLog,self.bt[j,i]:10:3);

Writeln(FLog);

end;

Readln;

end;

Procedure CEquacao.Simetriza;

var

AA:Matriz;

A - 60

i,j:Integer;

begin

self.CriaMatriz(AA,self.N,self.N);


For j:=1 to self.N do AA[i,j]:=self.A[i,j];


For j:=1 to self.N do self.A[i,j]:=(AA[i,j]+AA[j,i])/2;

end;

Procedure CEquacao.Solve;

var

i,ic,j,imax,jmax:Integer;

jc:Array of Integer;

VMax,aux,soma:Double;

AA,bbt:Matriz;

begin

// inicializacao

// -------------

Setlength(jc,self.N+1);

If self.Preserve then

begin

self.CriaMatriz(AA,self.N,self.N);

self.CriaMatriz(bbt,self.NR,self.N);


begin

For j:=1 to self.N do AA[i,j]:=self.A[i,j];

For j:=1 to self.NR do bbt[j,i]:=self.bt[j,i];

end;

end

else begin

AA:=self.A;

bbt:=self.bt;

end;

For j:=1 to self.N do jc[j]:=j;

//Substituicao

//------------

For ic:=1 to self.N do

begin

VMax:=0.0;

imax:=0;

jmax:=0;

For i:=ic to self.N do

begin

soma:=0.0;

For j:=ic to self.N do soma:=soma+abs(AA[i,j]);

For j:=ic to self.N do

If (abs(AA[i,j])/soma>VMax) then

begin

imax:=i;

jmax:=j;

VMax:=abs(AA[i,j])/soma;

end;

End;

// Deteccao de Casos Patologicos

// -----------------------------

If (imax=0) then

begin

Write('Matriz Nao Positiva Definida!',self.N-ic+1);

Readln;

End;

// troca de linhas ic e imax

// -------------------------

For j:=ic to self.N do

begin

aux:=AA[ic,j];

AA[ic,j]:=AA[imax,j];

AA[imax,j]:=aux;

End;

For j:=1 to self.NR do

begin

A - 61

aux:=bbt[j,ic];

bbt[j,ic]:=bbt[j,imax];

bbt[j,imax]:=aux;

end;

// troca de colunas ic e jmax

// --------------------------


begin

aux:=AA[i,ic];

AA[i,ic]:=AA[i,jmax];

AA[i,jmax]:=aux;

End;

j:=jc[ic];

jc[ic]:=jc[jmax];

jc[jmax]:=j;

// normalizacao da linha

// ---------------------

For j:=ic+1 to self.N do AA[ic,j]:=AA[ic,j]/AA[ic,ic];

For j:=1 to self.NR do bbt[j,ic]:=bbt[j,ic]/AA[ic,ic];

AA[ic,ic]:=1.0;

// condensacao

// -----------

For i:=ic+1 to self.N do

begin

For j:=ic+1 to self.N do AA[i,j]:=AA[i,j]-AA[ic,j]*AA[i,ic];

For j:=1 to self.NR do bbt[j,i]:=bbt[j,i]-bbt[j,ic]*AA[i,ic];

AA[i,ic]:=0.0;

End;

end;

//retrosubstituicao

//-----------------

For ic:=N downto 1 do

For i:=ic-1 downto 1 do

begin

For j:=1 to self.NR do bbt[j,i]:=bbt[j,i]-bbt[j,ic]*AA[i,ic];

AA[i,ic]:=0.0;

End;

//reordenacao

//-----------

For j:=1 to self.NR do For i:=1 to N do self.xt[j,jc[i]]:=bbt[j,i];

//If self.Preserve then

//begin

//self.Matrizdone(AA,self.N,self.N);

//self.Matrizdone(bbt,self.NR,self.N);

//end;

End;

Procedure CEquacao.Verifica;

Var

i,j:Integer;

r1,r2,r3:Matriz;

begin

self.CriaMatriz(r1,self.N,self.NR);



produtot(self.A,self.xt,self.N,self.N,self.NR,r1);

produtor(-1.0,self.bt,self.NR,self.N,r2);

soma(r1,r2,self.N,self.NR,r3);


For j:=1 to self.NR do

Writeln(r3[i,j]);

Readln;

end;

// Operacoes com Matrizes.

A - 62

Procedure produtor(r:Double;m:Matriz;n1,n2:Integer;Var mr:Matriz);

var i,j:integer;

begin

For i:=1 to n1 do

For j:=1 to n2 do

mr[j,i]:=r*m[i,j];

end;

Procedure produtot(m1,m2t:Matriz;n1,nc,n2:Integer;Var mr:Matriz);

var i,j,k:Integer;

begin

For i:=1 to n1 do

For j:=1 to n2 do

begin

mr[i,j]:=0;

For k:=1 to nc do mr[i,j]:=mr[i,j]+m1[i,k]*m2t[j,k];

end;

end;

Procedure produto(m1,m2:Matriz;n1,nc,n2:Integer;Var mr:Matriz);

var i,j,k:Integer;

begin

For i:=1 to n1 do

For j:=1 to n2 do

begin

mr[i,j]:=0;

For k:=1 to nc do mr[i,j]:=mr[i,j]+m1[i,k]*m2[k,j];

end;

end;

Procedure soma(m1,m2:Matriz;n1,n2:Integer;Var mr:Matriz);

var i,j:Integer;

begin

For i:=1 to n1 do

For j:=1 to n2 do

mr[i,j]:=m1[i,j]+m2[i,j];

end;

end.

A - 63

A.2.5 Unidade “OpMat”

//Operações envolvendo matrizes

unit OpMat;

{$mode objfpc}{$H+}

interface

uses

Classes, SysUtils;

type


function AddSubMat(m1,m2:Matriz;s:integer): Matriz;

function MultMat(m1,m2:Matriz): Matriz;

function MultMatTr(m2t,m1:Matriz): Matriz;

function MultEsc(m1:Matriz;k:double): Matriz;

implementation

//Multiplicação de matrizes

function MultMat(m1,m2:Matriz): Matriz;

var i,j,k,n1,nc,n2:Integer;

mr:Matriz;

begin

n1:=High(m1);

n2:=High(m2[0]);

nc:=High(m1[0]);

SetLength(mr,n1+1,n2+1);

For i:=0 to n1 do

For j:=0 to n2 do

begin

mr[i,j]:=0;

For k:=0 to nc do mr[i,j]:=mr[i,j]+m1[i,k]*m2[k,j];

end;

Result := mr;

end;

//Multiplicação de matrizes utilizando a matriz transposta de m2t

function MultMatTr(m2t,m1:Matriz): Matriz;

var i,j,k,n1,nc,n2:Integer;

mr:Matriz;

begin

n1:=High(m2t[0]);

n2:=High(m1[0]);

nc:=High(m2t);

SetLength(mr,n1+1,n2+1);

For i:=0 to n1 do

For j:=0 to n2 do

begin

mr[i,j]:=0;

For k:=0 to nc do mr[i,j]:=mr[i,j]+m1[k,j]*m2t[k,i];

end;

Result := mr;

end;

//Soma de matrizes

function AddSubMat(m1,m2:Matriz;s:integer): Matriz;

var i,j:Integer;

mr:Matriz;

begin

SetLength(mr,Length(m1),Length(m1[0]));

For i:=0 to High(m1) do

For j:=0 to High(m1[0]) do

begin

A - 64

mr[i,j]:=0;

mr[i,j]:=m1[i,j]+s*m2[i,j];

end;

Result := mr;

end;

//Multiplicação de um escalar por uma matriz

function MultEsc(m1:Matriz;k:double): Matriz;

var i,j:Integer;

mr:Matriz;

begin

SetLength(mr,Length(m1),Length(m1[0]));

For i:=0 to High(m1) do

For j:=0 to High(m1[0]) do

begin

mr[i,j]:=0;

mr[i,j]:=m1[i,j]*k;

end;

Result := mr;

end;

end.

A - 65

A.2.6 Unidade “NPGaussL”

//Rotina adaptada da referência bibliográfica [8].

//Cálculo de abcissas e pesos a utilizar na integração pela regra de quadratura

//de Gauss, em função do número de pontos definido.

unit NPGaussL;

{$mode objfpc}{$H+}

interface

uses

Classes, SysUtils;

type


function gauleg(n: integer): Matriz;

implementation

function gauleg(n: integer): Matriz;

var

x1,x2,xm,xl,z,z1,p1,p2,p3,pp: double;

i,j,m: integer;

xw: Matriz;

begin

SetLength(xw,n,2);

x1:=-1;x2:=1;

m:=(n+1) div 2;

xm:=0.5*(x2+x1);

xl:=0.5*(x2-x1);

for i:=1 to m do

begin

z:=cos(pi*(i-0.25)/(n+0.5));

z1:=z+1;

while abs(z-z1)>1e-15 do

begin

p1:=1;

p2:=0;

for j:=1 to n do

begin

p3:=p2;

p2:=p1;

p1:=((2*j-1)*z*p2-(j-1)*p3)/j;

end;

pp:=n*(z*p1-p2)/(z*z-1);

z1:=z;

z:=z1-p1/pp;

end;

xw[i-1,0]:=xm-xl*z;

xw[n-i,0]:=xm+xl*z;

xw[i-1,1]:=2*xl/((1-z*z)*pp*pp);

xw[n-i,1]:=xw[i-1,1];

end;

Result := xw;

end;

end.

Documents

Análise de Estruturas Axissimétricas - Autenticação · de estruturas mencionadas: uma utiliza soluções analíticas e a outra o método dos elementos finitos. A primeira destas