Uma Introdução ao Método dos Elementos Finitos Marco

Uma Introdução ao Método dos Elementos Finitose suas aplicações

Marco Alexandre Claudino

Trabalho de Conclusão de Curso apresentadoao

Instituto de Matemática e Estatísticada

Universidade de São Paulopara

obtenção do títulode

Bacharel em Matemática Aplicada e Computacional

Orientador: Prof. Dr. Nelson Mugayar Kuhl

Durante o desenvolvimento deste trabalho o autor recebeu auxílio nanceiro da Reitoria da

Universidade de São Paulo.

São Paulo, 13 de dezembro de 2012

Uma Introdução ao Método dos Elementos Finitos e suas aplicações

Esta é a versão original do trabalho de conclusão de curso elaborado pelo

candidato Marco Alexandre Claudino, tal como

submetido à banca avaliadora.

Banca avaliadora:

• Prof. Dr. Nelson Mugayar Kuhl (orientador) - IME - USP

• Prof. Dr. Alexandre Megiorin Roma - IME - USP

• Prof. Dr. Saulo Rabello Maciel de Barros - IME - USP

Agradecimentos

Primeiramente, agradeço a Deus pela saúde, pela coragem e pela fé que sempre me guiou frente

aos momentos mais complicados da vida.

Aos meus pais, Marco e Oneide, e minha irmã, Cristiane pelo apoio incondicional dado durante

toda a graduação.

A minha noiva Stefany pelo carinho, atenção e paciência tidos comigo.

Aos companheiros de graduação Guilherme Silva Salomão, Willian Hans Corrêa (79), Lucas Reis

e Eli Enrico Carnette pela companhia durante a jornada acadêmica.

Ao professor Nelson Mugayar Kuhl, pelos valorosos conselhos dados durante as reuniões de

iniciação cientíca.

Aos professores Alexandre Megiorin Roma, Celso Massatoshi Furukawa (POLI), Daniel Victor

Tausk, Fabio Armando Tal, Fabio Gagliardi Cozman (POLI), Leonidas de Oliveira Brandão, Luis

Carlos de Castro Santos, Maria Izabel Ramalho Martins, Newton Maruyama (POLI), Orlando

Francisco Lopes, Sônia Regina Leite Garcia e Zara Issa Abud cujos conhecimentos transmitidos em

aula foram muito importantes para o meu crescimento acadêmico. Destaco também os alunos de

pós graduação Álvaro Junio Pereira Franco e Rodrigo Roque Dias com os quais tive contato nos

cursos de verão do IME.

i

"Toda essa multidão tem no peito amor e procura a paz

E apesar de tudo a esperança não se desfaz

Olhando a or que nasce no chão daquele que tem amor

Olho pro céu e sinto crescer a fé no meu salvador

Jesus Cristo"

Roberto Carlos

Resumo

Claudino, M.A. Uma Introdução ao Método dos Elementos Finitos e suas aplicações.

2012. 85 f. Trabalho de Conclusão de Curso - Instituto de Matemática e Estatística, Universidade

de São Paulo, São Paulo, 2012.

O Método dos Elementos Finitos (MEF) trata-se de um método numérico que busca obter uma

aproximação da equação diferencial sobre cada componente de uma malha de partições do domínio

(os chamados elementos) utilizando funções de interpolação construídas sobre cada um deles. Este

trabalho procura realizar uma primeira apresentação do MEF enfatizando alguns aspectos com-

putacionais do problema e apresenta algumas simulações computacionais de problemas envolvendo

distribuição de calor em regime estacionário e elasticidade linear. A construção do trabalho segue

de perto a obra [FB07] na apresentação do método e é complementada com aspectos computacionais.

Palavras-chave: Método dos Elementos Finitos, Equações Diferenciais Parciais, Modelagem Ma-

temática.

iii

iv

Sumário

Lista de Abreviaturas vii

Lista de Símbolos ix

Lista de Figuras xi

Lista de Tabelas xiii

1 Introdução 1

1.1 Pré-Requisitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Construção do método e denições iniciais . . . . . . . . . . . . . . . . . . . . . . . . 1

1.3 Organização do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 Introdução ao Método dos Elementos Finitos Unidimensional 5

2.1 Formulações do Problema Unidimensional . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Tipos de elementos e aproximações . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3 Formulação geral do método para problemas unidimensionais . . . . . . . . . . . . . 9

2.3.1 Elementos Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3.2 Elementos Quadraticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.3 Elementos de ordem k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Aspectos computacionais do problema unidimensional 19

3.1 Integração Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1.1 Quadratura de Gauss-Legendre . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Armazenamento do Sistema Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.3 Resolução de sistemas lineares com estrutura de banda . . . . . . . . . . . . . . . . . 24

3.3.1 Elementos com dois nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.3.2 Elementos com três nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3.3 Elementos com k nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.4 Convergência numérica do algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Simulações unidimensionais 29

4.1 Apresentação do Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.1.1 Resolução com o uso de Elementos Lineares . . . . . . . . . . . . . . . . . . . 30

4.1.2 Resolução com o uso de Elementos Quadraticos . . . . . . . . . . . . . . . . . 31

v

vi SUMÁRIO

5 Introdução ao Método dos Elementos Finitos bidimensional 35

5.1 Formulações do Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

5.2 Tipos de elementos e aproximações . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

5.2.1 Elemento Triangular com três nós . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.2.2 Coordenadas Triangulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.3 Formulação geral do método para problemas bidimensionais . . . . . . . . . . . . . . 43

5.3.1 Geração da Matriz Global de Rigidez e do Vetor Global de Forças . . . . . . 44

6 Aspectos computacionais do problema bidimensional 45

6.1 Geração da malha do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

6.1.1 Dados de Entrada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

6.1.2 Dados de Saída . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

6.2 Quadraturas numéricas para problemas bidimensionais . . . . . . . . . . . . . . . . . 47

6.2.1 Integração sobre elementos triangulares . . . . . . . . . . . . . . . . . . . . . 48

6.3 Armazenamento de matrizes esparsas . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

6.4 Método dos Gradientes Conjugados (MGC) . . . . . . . . . . . . . . . . . . . . . . . 49

6.4.1 Idéia do Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

6.4.2 Critério de Parada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

6.4.3 Propriedades do Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

6.4.4 Algoritmo do Método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

6.4.5 Pré-condicionamento do sistema . . . . . . . . . . . . . . . . . . . . . . . . . 52

6.5 Convergência numérica do algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

7 Simulações bidimensionais 53

7.1 Problema de Poisson no Quadrado . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

7.1.1 Simulação Computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

7.1.2 Análise de erro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

7.1.3 Ordem de convergência do Método . . . . . . . . . . . . . . . . . . . . . . . . 55

7.2 Distribuição de calor em uma chaminé . . . . . . . . . . . . . . . . . . . . . . . . . . 56

7.2.1 Modelagem do Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56


7.2.3 Análise dos resultados obtidos . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

7.3 Análise de deformações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

7.3.1 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

7.3.2 Tensões e trações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

7.3.3 Equilíbrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

7.3.4 Aplicação do Método dos Elementos Finitos . . . . . . . . . . . . . . . . . . . 62

7.3.5 Utilização de aproximações lineares em elementos triangulares com três nós . 63


Referências Bibliográcas 69

Lista de Abreviaturas

CSR Compactação esparsa por linhas (Compressed Sparse Row)

MDF Método das Diferenças Finitas (Finite Dierence Method)

MEF Método dos Elementos Finitos (Finite Element Method)

MGC Método dos Gradientes Conjugados (Conjugate Gradient Method)

vii

viii LISTA DE ABREVIATURAS

Lista de Símbolos

Ω Domínio do Problema

Γ Contorno do domínio Ω

Γe Trecho do contorno onde é conhecida a solução do problema (Condição de Dirichlet)

Γn Trecho do contorno onde é conhecida a componente normal da derivada da função

(Condição de Neumann)

Ki Matriz de rigidez do i-ésimo elemento

f iΩ Vetor de forças de campo do i-ésimo elemento

f iΓ Vetor de forças de contorno do i-ésimo elemento

K Matriz de rigidez global

fΩ Vetor global de forças de campo

fΓ Vetor global de forças de contorno

ix

x LISTA DE SÍMBOLOS

Lista de Figuras

2.1 Funções teste para elementos lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Funções teste para elementos quadráticos . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Partição do domínio para elementos lineares. . . . . . . . . . . . . . . . . . . . . . . 11

2.4 Partição do domínio para elementos quadrático. . . . . . . . . . . . . . . . . . . . . . 13

2.5 Partição do domínio para elementos de ordem k. . . . . . . . . . . . . . . . . . . . . 16

3.1 Matriz com propriedade de banda. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.1 Visualização do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2 Aproximações obtidas com elementos lineares para 2, 4, 8 e 16 elementos. A solução

exata é dada pela curva azul e as aproximações são dadas pelas curvas vermelhas. . . 30

4.3 A direita temos o gráco da aproximação obtida e a esquerda o gráco da solução

exata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.4 Erro cometido na aproximação com fator de escala 1e+12. . . . . . . . . . . . . . . . 32

4.5 As curvas em vermelho apresentam as aproximações numéricas e as curvas em azul

os valores da solução do problema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.1 Dois elementos quaisquer possuem em comum apenas o segmento que os une. . . . . 38

5.2 Elemento triangular de três nós . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.3 Continuidade entre elementos adjacentes . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.4 Coordenadas triangulares de um ponto P . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.5 Triângulo de referência em coordenadas triangulares . . . . . . . . . . . . . . . . . . 42

5.6 Montagem da matriz global de rigidez e do vetor global de forças. [Fla] . . . . . . . . 44

6.1 Arquivo de entrada para o EasyMesh do problema 7.8. . . . . . . . . . . . . . . . . . 46

6.2 Dados do arquivo entrada.e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

6.3 Dados do arquivo entrada.s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

7.1 Visualização das aproximações numéricas obtidas. . . . . . . . . . . . . . . . . . . . . 54

7.2 Renamento da malha de domínio do problema. . . . . . . . . . . . . . . . . . . . . . 55

7.3 Comportamento da ordem estimada do método em relação ao número de iterações

de renamento de malha do problema. . . . . . . . . . . . . . . . . . . . . . . . . . . 56

7.4 Modelagem matemática do problema [FB07]. . . . . . . . . . . . . . . . . . . . . . . 57

7.5 Modelagem do problema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

7.6 Renamento da malha de domínio do problema. . . . . . . . . . . . . . . . . . . . . . 59

7.7 Distribuição da temperatura com o MEF utilizando 1844 elementos triangulares. . . 60

xi

xii LISTA DE FIGURAS

7.8 Descrição do problema que será abordado: a chave inglesa encontra-se presa a um

parafuso que não se movimenta enquanto sofre a ação de uma força constante. . . . . 60

7.9 Visualização dos componentes do deslocamento [FB07]. . . . . . . . . . . . . . . . . . 61

7.10 Visualização dos componentes da tensão [FB07]. . . . . . . . . . . . . . . . . . . . . 62

7.11 Visualização dos componentes das forças em equilíbrio no quadrado unitário [FB07]. 62

7.12 Componentes do deslocamento em cada nó do elemento triangular. . . . . . . . . . . 64

7.13 Particionamento do domínio em elementos triangulares com constante de renamento

igual a 0.25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

7.14 Enumeração dos nós da malha com constante de renamento igual a 0.25. . . . . . . 66

7.15 Enumeração dos elementos da malha com constante de renamento igual a 0.25. . . . 66

7.16 Norma do deslocamento ocorrido em cada ponto do domínio utilizando-se da malha

com constante de renamento igual a 0.25. . . . . . . . . . . . . . . . . . . . . . . . . 67

7.17 Visualização do deslocamento da chave com um fator de escala 104 (malha com

constante de renamento igual a 0.25). . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Lista de Tabelas

3.1 Nós da quadratura de Gauss-Legendre e suas ponderações. . . . . . . . . . . . . . . . 22

4.1 Erros obtidos na simulação utilizando elementos lineares. . . . . . . . . . . . . . . . . 31

4.2 Erros obtidos na simulação utilizando elementos quadráticos. . . . . . . . . . . . . . 33

6.1 Nós da quadratura de Gauss-Legendre e suas ponderações para domínios triangulares 48

7.1 Dados obtidos com o renamento da malha. . . . . . . . . . . . . . . . . . . . . . . . 55

7.2 Dados obtidos com o renamento da malha do problema, onde n(i), e(i) e h(i)

representam a quantidade de nós, o erro máximo obtido e o tamanho h obtido na

malha i, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

7.3 Temperaturas (em C) obtidas na análise do problema com renamento da malha. . 59

7.4 Deslocamento obtido em alguns pontos da chave com fator de escala 104. . . . . . . . 66

xiii

xiv LISTA DE TABELAS

Capítulo 1

Introdução

"O chamado Método dos Elementos Finitos (MEF) consiste em diferentes métodos numéricosque aproximam a solução de problemas de valor de fronteira descritos tanto por equações diferenciaisordinárias quanto por equações diferenciais parciais através da subdivisão da geometria do problemaem elementos menores, chamados elementos nitos, nos quais a aproximação da solução exata podeser obtida por interpolação de uma solução aproximada. Atualmente o MEF encontra aplicação empraticamente todas as áreas de engenharia, como na análise de tensões e deformações, transferênciade calor, mecânica dos uidos e reologia, eletromagnetismo, etc..." [Shi]. Tal descrição é basicamenteencontrada em livros, artigos e outros materiais que abordam o Método dos Elementos Finitos.

De fato, obter a solução destes problemas sob a forma fechada costuma ser um trabalho bastantecomplicado (e algumas vezes impossível) de forma que faz-se necessário o uso de aproximaçõesnuméricas da solução do problema. Neste contexto, a primeira técnica normalmente apresentada éo Método das Diferenças Finitas (MDF).

Outra possível estratégia bastante utilizada é a aplicação do Método dos Elementos Finitos(MEF). Trata-se de um método numérico que busca obter uma boa aproximação da solução daequação diferencial sobre cada componente de uma malha de partições do domínio (os chamadoselementos) utilizando funções de interpolação construídas sobre cada um deles.

Uma das diculdades do estudo do MEF reside nas diferentes abordagens feitas ao método: sobo enfoque da engenharia o método é tratado como uma técnica de modelagem de problemas. Poroutro lado, sob o enfoque matemático, o método é apresentado como a minimização de um certofuncional em um determinado espaço de funções. Ambos os tratamentos apresentam complicacõese podem ser bastante trabalhosos para um aluno em início de graduação.

Desta forma, o objetivo deste trabalho é realizar uma primeira apresentação do MEF enfatizandoalguns aspectos computacionais do problema e apresentando algumas simulações computacionaisde problemas envolvendo distribuição de calor em regime estacionário e elasticidade linear. A cons-trução do trabalho segue de perto a obra [FB07] na apresentação do método e é complementadacom aspectos computacionais.

1.1 Pré-Requisitos

Durante a elaboração deste trabalho procuramos utilizar apenas de ferramentas normalmenteapresentadas nos cursos de cálculo vetorial e análise numérica, enfatizando as idéias utilizadas eos conceitos importantes. Desta forma, temos que tais pré-requisitos são aconselháveis, porém nãoindispensáveis para o entendimento do conteúdo apresentado.

1.2 Construção do método e denições iniciais

A construção do MEF divide-se basicamente em três partes: a obtenção da formulação fraca doproblema em estudo, a construção das funções teste e de ponderação e a resolução do sistema linearassociado.

1

2 INTRODUÇÃO 1.3

Iremos iniciar o estudo do método analisando sua aplicação em equações diferenciais ordinárias edepois, utilizando as idéias desenvolvidas, iremos generalizar a aplicação para equações diferenciaisparciais, onde o nosso foco serão as equações elípticas.

Para isto, apresentaremos aqui algumas denições e terminologias que serão utilizadas duranteo texto:

• Formulação forte: Dizemos que um dado problema encontra-se em sua formulação fortequando apresentado por uma equação diferencial denida sobre um determinado domínio Ω(Ω ⊂ R no caso unidimensional ou Ω ⊂ R2 no caso bidimensional) e são fornecidas tam-bém determinadas condições sobre a função na fronteira de Ω, denotada por Γ. O trecho dafronteira onde a condição fornecida é o próprio valor da função sera denominado condiçãode contorno essencial (Γe) do problema, enquanto que o trecho onde a condição fornecidaenvolve a derivada da função será denominado condição de contorno natural (Γn) do pro-blema. Tal terminologia vem do fato de que as condições naturais surgem automaticamenteda equivalência entre uma formulação fraca do problema e a formulação forte, enquanto quea condição essencial deve estar presente em ambas as formulações.

• Formulação fraca: Dizemos que um problema encontra-se em sua formulação fraca quandodescrito por uma formulação integral onde o integrando satisfaz a uma determinada proprie-dade.

Salientamos aqui que esta terminologia não é exclusiva: existem diversas formulações conside-radas forte na literatura, assim como existem outras consideradas formulações fracas. Porém, ouso de tais terminologias fez-se necessário, uma vez que o trabalho buscou acompanhar de perto aobra [FB07] no aspecto da apresentação da teoria do método, complementando-a apenas na exem-plicação do método das simulações computacionais.

1.3 Organização do Trabalho

Iniciamos nosso estudo do MEF apresentando a construção e a aplicação do método em equaçõesdiferenciais ordinárias. O Capítulo 2 apresenta como obter a formulação fraca associada ao problema,aborda a geração dos elementos unidimensionais e a construção das aproximações sobre cada umdeles. Ao nal, apresentamos como é gerado o sistema linear que fornece os coecientes a seremutilizados na aproximação global do problema.

No Capítulo 3 abordamos aspectos computacionais ligados ao problema unidimensional comoa solução numérica das integrais obtidas da formulação fraca, o armazenamento de forma ecienteda matriz de rigidez do sistema e as técnicas para resolução do sistema linear, utilizando-se daestrutura de banda do sistema. Ao nal, é apresentada uma estratégia numérica para avaliarmos ocomportamento do erro em função do renamento da malha de elementos.

Em seguida, o Capítulo 4 apresenta os resultados da simulação de um problema (extraidode [FB07]) envolvendo distribuição de calor em uma barra.

No Capítulo 5 é feita a introdução ao MEF bidimensional, com a construção da formulação fracaequivalente, a obtenção das funções teste sobre elementos triangulares e a construção do sistemalinear associado.

Já o Capítulo 6 aborda os aspectos computacionais relativos aos problemas bidimensionais,como a utilização de um gerador de malhas triangulares, quadraturas numéricas para a avaliaçãode integrais duplas sobre regiões triangulares e o Método dos Gradientes Conjugados (MGC) para aresolução do sistema linear associado ao problema bidimensional. Por m, generalizamos a estratégiaapresentada anteriormente para avaliarmos o comportamento do erro em relação ao renamento damalha de estudo.

Por m, no Capítulo 7 apresentamos três simulações utilizando o MEF: um problema de Pois-son em um quadrado unitário, o problema da condução de calor em uma chaminé e a análise dadeformação sofrida por uma chave inglesa após a ação de uma força constante.

1.3 ORGANIZAÇÃO DO TRABALHO 3

O uxograma a seguir descreve a construção do trabalho mostrando os passos a serem executa-dos para a aplicação do MEF.

Formulaçãoforte doproblema

Formulaçãofraca doproblema

Geração damalha deelementos

Geração dasmatrizes locais

de rigidez

Cálculo dasmatrizes derigidez locais

Geração dosvetores locaisde força decampo e decontorno

Cálculo dosvetores deforças locais

Geração damatriz globalde rigidez Geração do

vetor globalde forças

Resoluçãodo SistemaKd = f

Obtenção dovalor da apro-ximação porinterpolação

Ponto de interesse

Aproximação obtida

4 INTRODUÇÃO 1.3

Capítulo 2

Introdução ao Método dos Elementos

Finitos Unidimensional

2.1 Formulações do Problema Unidimensional

Para os modelos que serão abordados neste capítulo, dizemos que um problema encontra-se emsua formulação forte quando a equação diferencial ordinária é apresentada da seguinte forma:

- Obter u : [0, l]→ R de classe C2 que satisfaça as seguintes condições:

ddx

(p(x)dudx

)= f(x), ∀x ∈ (0, l) (2.1)

p(x)du

dx

∣∣∣x=0

= t0 (2.2)

u∣∣∣x=l

= u0 (2.3)

onde p(x) e f(x) são funções suaves, p(x) > 0, 0 ≤ x ≤ l, u0 e t0 ∈ R.Para desenvolvermos as equações do Método dos Elementos Finitos devemos inicialmente re-

formular o problema, de maneira a escrevê-lo como uma formulação integral, a qual denominamosformulação fraca.

A obtenção da formulação fraca inicia-se pela escolha de uma função w : [0, l] → R, de classeC1, denominada função peso (ou função de ponderação), que é escolhida de maneira arbitráriasatisfazendo a propriedade de que w(l) = 0.

Após a escolha da função peso, multiplicamos as equações (2.1) e (2.2) por w e integramos oprimeiro produto obtido no domínio do problema. Assim, obtemos:∫ l

0w(x)

d

dx

(p(x)

du

dx

)dx =

∫ l

0w(x)f(x)dx; (2.4)

(w(x)p(x)

du

dx

) ∣∣∣x=0

= w(0)t0 (2.5)

Utilizando integração por partes no termo a esquerda de (2.4), temos que:∫ l

0w(x)

d

dx

(p(x)

du

dx

)dx =

(w(x)p(x)

du

dx

) ∣∣∣l0−∫ l

0

(dw

dx

)(p(x)

du

dx

)dx (2.6)

5

6 INTRODUÇÃO AO MÉTODO DOS ELEMENTOS FINITOS UNIDIMENSIONAL 2.1

Como

(w(x)p(x)dudx

) ∣∣∣l0

=

=0, pois w(l)=0︷︸︸︷w(x)p(x)

du

dx

∣∣∣x=l−w(x)p(x)dudx

∣∣∣x=0

= −w(x)p(x)dudx

∣∣∣x=0

= −w(0)t0

(2.7)

e assim, das equações (2.4) e (2.7), temos que:∫ l

0

(dw

dx

)(p(x)

du

dx

)dx = −

∫ l

0w(x)f(x)dx− w(0)t0. (2.8)

Desta forma, dizemos que um problema encontra-se em sua formulação fraca quando é apresen-tado da seguinte forma:

- Obter u : [0, l]→ R , de classe C2 com u(l) = u0 que satisfaça:∫ l

0

(dw

dx

)(p(x)

du

dx

)dx = −

∫ l

0w(x)f(x)dx− w(0)t0, ∀w : w(l) = 0, (2.9)

onde w : [0, l]→ R de classe C1.Um questionamento natural a ser feito é sobre a equivalência entre as duas formulações. De

fato, provamos que a uma função que satisfaz a formulação forte satisfaz a formulação fraca. Vamosagora provar que uma solução da formulação fraca também satisfaz a formulação forte.

Para isso, simplesmente iremos inverter os passos realizados para a obtenção da formulaçãofraca. Ao invés de utilizarmos a integração por partes para eliminar a segunda derivada de u(x)invertemos a fórmula para obter uma integral com uma derivada maior e um termo de contorno.Assim, reorganizando os termos de (2.6) obtemos:∫ l

0

(dw

dx

)(p(x)

du

dx

)dx =

(w(x)p(x)

du

dx

) ∣∣∣l0−∫ l

0w(x)

d

dx

(p(x)

du

dx

)dx. (2.10)

Substituindo (2.10) em (2.9) e reorganizando os termos, obtemos:∫ l

0w

[(p(x)

du

dx

)− f(x)

]dx− w(x)(t0 − p(x)

du

dx)∣∣∣x=0

= 0,∀w : w(l) = 0. (2.11)

Devido a arbitrariedade da função peso, tomamos w(x) da forma

w(x) = ψ(x)

[d

dx

(p(x)

du

dx

)+ f(x)

], (2.12)

em que ψ(x) é uma função suave com ψ(x) > 0,∀x ∈ (0, l) e ψ(x) = 0 nos contornos. Uma funçãoψ(x) que tem essas propriedades, por exemplo, é a função ψ(x) = x(l − x). Desta forma, temosque a condição w(l) = 0 no contorno essencial é automaticamente satisfeita e, assim, w(x) pode serescolhida como uma função peso. Aplicando (2.12) em (2.11), obtemos:∫ l

0ψ(x)

[d

dx

(p(x)

du

dx

)+ f(x)

]2

dx = 0. (2.13)

Devido a construção da função peso, temos que o termo de contorno desaparece, pois w(0) = 0.Como o integrando acima é o produto de uma função positiva e o quadrado de uma função, segueque ele é positivo em cada ponto do domínio do problema. Desta forma, a igualdade (2.13) ocorre se,e somente se o integrando for igual a zero em cada um dos pontos do domínio, obtendo novamente(2.1). Por outro lado, sendo a igualdade (2.13) vericada, temos de (2.11) que:[

w(x)(t0 − p(x)du

dx)

] ∣∣∣x=0

= 0,∀w : w(l) = 0. (2.14)

2.2 TIPOS DE ELEMENTOS E APROXIMAÇÕES 7

Como a função peso é arbitrária, podemos tomá-la de maneira que w(0) = 1 obtendo então que acondição de contorno natural (2.2) é satisfeita.

A última equação da formulação forte (2.3) é satisfeita por todas as funções teste que satisfazema formulação fraca, pois exigimos que u(l) = u0.

Concluimos então que as funções que satisfazem a formulação fraca satisfazem também a for-mulação forte. Observe que as provas de equivalência estão fortemente baseadas no fato de que aformulação fraca é válida para qualquer função sucientemente diferenciável e que a função peso éarbitrariamente escolhida, porém, tomamos a liberdade de escolher uma que satiszesse as condiçõesque necessitavamos para a conclusão do resultado.

Um fato importante a ser enfatizado é a suavidade das funções envolvidas na solução da formu-lação fraca do problema. Se as funções teste e as funções peso forem muito irregulares, em relaçãoa descontinuidades, temos que a aplicação da formulação fraca torna-se bastante complicada. Poroutro lado, se uma dada função que satisfaz a formulação fraca não for duas vezes diferenciavelem algum ponto do domínio, temos que tal função não satisfaz a formulação forte. Desta forma,podemos observar que a formulação fraca permite relaxar as hipóteses sobre a solução do problemao que a torna, neste sentido, uma formulação mais geral do problema.

2.2 Tipos de elementos e aproximações

Iniciamos a aplicação do MEF por meio da colocação de nós no intervalo de estudo e pela divisãodo domínio do problema em elementos. Em cada um destes elementos, construimos funções espe-cícas, as quais denominamos funções teste, de forma que a aproximação obtida em cada elementoconsiste numa combinação linear das funções teste. Tais funções são escolhidas cuidadosamente, demaneira a obter a convergência do método.

Não apresentaremos aqui a análise matemática que aborda a convergência analítica do métodomas, seguindo a teoria do método, as funções escolhidas devem ser sucientemente diferenciáveis eapresentar boas propriedades operacionais. Como funções polinomiais são innitamente diferenciá-veis e sua a integração é facilmente calculável, iremos tomar tanto as funções peso como as funçõesde teste como polinômios.

Antes de analisarmos como obter os polinômios para as funções peso e as funções teste, vamosapresentar agora algumas notações que serão frequentemente utilizadas.

Denimos como domínio Ω do problema o intervalo [0, l]. Um elemento em Ω é um subintervalodenotado por Ωi, i = 1, 2, · · · , nel , formado por partições do domínio contendo uma quantidade nde pontos da malha, onde

Ωi = [xk, xk+n], xk ∈ Ω, k = 1, 2, · · · , (nt − n)

nel⋃i=1

Ωi = Ω = [0, l]

onde nt é o número total de pontos do domínio e nel é o número elementos em Ω.Sendo Ω um intervalo em R, temos que a fronteira do domínio é constituída por dois pontos.

Denotamos por Γn o ponto da fronteira onde é conhecida a condição natural e Γe o ponto onde éconhecida a condição essencial do problema. A nomenclatura dada vem do fato de que a condição decontorno essencial deve ser fornecida em ambas as formulações enquanto que a condição de contornonatural surge naturalmente na demonstração da equivalência entre as formulações.

Dada uma malha de pontos em Ω, a aproximação global por elementos nitos será indicada poruh(x), onde h refere-se a maior distância entre dois pontos da malha. A restrição da aproximaçãoao elemento e será indicada por ue(x), sendo que iremos impor a condição de ue(x) = 0 se x /∈ Ωe.Para a indicação dos nós pertencentes a cada um dos elementos, usaramos o índice inferior. Assim,por exemplo, xj1 representa o nó 1 do j-ésimo elemento.


Na construção de uma aproximação por elementos nitos, nosso objetivo é construir um po-linômo que aproxime bem a solução do problema em cada um dos elementos do domínio. Destaforma, temos que para cada elemento iremos construir um polinômio ue(x) da forma

ue(x) =

γe0 + γe1x+ γe2x

2 + · · · , se x ∈ Ωe

0, se x /∈ Ωe (2.15)

onde os γei ∈ R, i = 0, 1, 2, · · · são constantes escolhidas de maneira que não ocorram saltos naaproximação obtida sobre o ponto comum a dois elementos consecutivos. Desta forma, para umamalha contendo n pontos por elemento e nel elementos, a igualdade ue(xen) = ue+1(xe+1

1 ) deve servericada para e = 1, 2, · · · , (nel − 1).

É possível também desenvolver funções teste de forma que as derivadas nos contornos doselementos sejam iguais e a aproximação obtida seja de classe C1. No entanto, iremos enfatizar autilização de funções teste contínuas, isto é, de classe C0.

A construção das funções teste de classe C0 em uma dimensão são facilmente obtidas por meiodos Interpoladores de Lagrange que consistem em polinômios denidos em intervalos do domíniosobre um determinado número de pontos, neste caso, os nós do intervalo. Vamos exemplicar aobtenção das funções teste para elementos com dois e três pontos e, em seguida, iremos generalizara idéia utilizada para n pontos.1

• Funções teste lineares: Consideramos elementos que possuem dois pontos xe1 e xe2 e duasfunções teste, denidas por:

N e1 (x) =

(x− xe2)

(xe1 − xe2)

N e2 (x) =

(x− xe1)

(xe2 − xe1)

Figura 2.1: Funções teste para elementos lineares

• Funções teste quadráticas: Consideramos elementos que possuem três pontos xe1, xe2 e xe3 e

três funções teste, denidas por:

N e1 (x) =

(x− xe2)(x− xe3)

(xe1 − xe2)(xe1 − xe3)

N e2 (x) =

(x− xe1)(x− xe3)

(xe2 − xe1)(xe2 − xe3)

N e3 (x) =

(x− xe1)(x− xe2)

(xe3 − xe2)(xe3 − xe1)

1Em grande parte da literatura os elementos onde são denidas as funções teste são denotados pelas propriedadesdas funções teste, isto é, elementos cujas funções teste denidas sobre eles são lineares recebem o nome de elementoslineares, elementos cujas funções teste denidas sobre eles são quadraticas recebem o nome de elementos quadráticose assim por diante.

2.3 FORMULAÇÃO GERAL DO MÉTODO PARA PROBLEMAS UNIDIMENSIONAIS 9

Figura 2.2: Funções teste para elementos quadráticos

• Funções teste de ordem (n-1): De maneira geral, consideramos elementos que possuem npontos xei , i = 1, 2, · · · , n e n funções teste, denidas por:

N ei (x) =

n∏i,j=1;i 6=j

(x− xej)(xei − xej)

Observamos aqui que a construção das funções de teste desta forma apresentam a seguintepropriedade:

N ei (xej) =

1, se i = j0, se i 6= j

(2.16)

com e = 1, 2, · · · , nel.Desta forma, cada uma das funções teste serão polinômios de Lagrange não nulos apenas em

um determinado xei pertencente ao elemento e denido. Além disso, note que as funções teste paraelementos de ordem n− 1 formam uma base para os polinômios de grau n− 1.

O uso de tais polinômios fornecerá propriedades para a geração de algoritmos ecazes na reso-lução do problema.

2.3 Formulação geral do método para problemas unidimensionais

De posse da formulação fraca do problema (2.9), das funções de teste e de ponderação, partiremosagora para a formulação do Método dos Elementos Finitos por meio da discretizção do domínio emestudo, onde buscamos obter um número nito de equações algébricas oriundas da discretização daformulação fraca.

Inicialmente, particionamos o intervalo de domínio do problema Ω em nt pontos e denimoso tipo de elemento que será utilizado na formulação do problema. Em cada um dos elementos,denimos as funções teste obtidas pelos polinômios de Lagrange em relação aos pontos pertencentesa cada um dos elementos do domínio e tomaremos sobre este elemento uma aproximação ue(x) daforma:

ue(x) = αe1Ne1 + αe2N

e2 + αe3N

e3 + αe4N

e4 + · · ·+ αenN

en

onde cada uma das constantes αei , i = 1, 2, · · · , n serão escolhidas conforme vimos anteriormenteem (2.15). Observe que o número de coecientes αe é o mesmo número de pontos do elemento e N e

i

são as funções teste obtidas pelos interpoladores de Lagrange para n pontos por elemento.A aproximação global uh será obtida pela soma das aproximações locais, isto é

uh(x) =

nel∑e=1

ue(x)

onde cada ue(x) = 0 se x /∈ Ωe.


Como existem descontinuidades entre as funções teste entre cada um dos elementos, a integraçãoda formulação fraca deve ser avaliada como a somatória das integrações em cada um dos elementosda malha. Desta forma:

nel∑e=1

∫ xen

xe1

(dwe

dx

)(p(x)

due

dx

)dx+

∫ xen

xe1

wef(x)dx+

(wep(x)

due

dx

) ∣∣∣x=0

= 0 (2.17)

onde e representa o respectivo elemento onde a integral está sendo avaliada, n o número de pontoscontidos no elemento e, xe1 o primeiro ponto do elemento e xen o último ponto do elemento.

Denimos então os vetores das funções teste, das derivadas das funções teste e dos coecientesque satisfazem 2.17 restritos ao elemento e:

Ne(x) = [N e1 (x), N e

2 (x), · · · , N en(x)] (2.18)

Be(x) =

[dN e

1 (x)

dx,dN e

2 (x)

dx, · · · , dN

en(x)

dx

](2.19)

αe = [αe1, αe2, · · · , αen] (2.20)

Utilizando funções teste e funções peso da forma

ue(x) = Ne(x)αeT ⇒ due(x)dx = Be(x)αeT

we(x) = Ne(x)weT ⇒ dwe(x)dx = Be(x)weT

(2.21)

e substituindo em (2.17) temos

nel∑e=1

(we)T

(∫ xek

xe1

(Be)T p(x)Bedx

)︸︷︷︸

Ke

αe +

∫ xek

xe1

(Ne)T f(x)dx+

((Ne)T p(x)

due

dx

)︸︷︷︸

fe

∣∣∣x=0

= 0 (2.22)

Da equação anterior, denimos dois componentes que serão bastante úteis no Método dos Ele-mentos Finitos:

• A Matriz de Rigidez Ke do elemento e:

Ke =

∫ xen

xe1

(Be)T p(x)Bedx =

∫Ωe

(Be)T p(x)Bedx (2.23)

• O Vetor de Forças Externas fe do elemento e:

fe =

∫ xen

xe1

(Ne)T f(x)dx+ (Ne)Tw(0) =

∫Ωe

(Ne)T f(x)dx︸︷︷︸fΩe

+ (Ne)T t0

∣∣∣Γet︸︷︷︸

fΓe

(2.24)

Desta forma, escolhido o tipo de elemento que será utilizado, as matrizes de rigidez e os vetoresde forças externas de cada elemento podem ser calculados por meio das integrais descritas acima.Iremos apresentar, em detalhes, montagem do sistema para elementos lineares e quadráticos e, emseguida, generalizaremos para elementos de ordens superiores.

Um detalhe importante a ser observado é que os valores de we são arbitrários e, sendo o numéro


de graus de liberdade do problema igual ao número de constantes α a serem obtidas, podemosconsiderar que os valores de we na relação (2.22) não inuenciam na obtenção do sistema que geraa solução da formulação fraca.

2.3.1 Elementos Lineares

Numa primeira análise do método, iremos particionar Ω em nt pontos igualmente espaçados,denotando por le a distância entre dois pontos consecutivos da malha. Com o uso de elementoslineares temos que cada dois pontos da malha deve ser considerado um elemento Ωe. Temos entãoque:

Figura 2.3: Partição do domínio para elementos lineares.

Para cada um dos elementos, sendo xe1 e xe2 os pontos que compõem o elemento e, temos que asfunções teste são dadas por

Ne(x) = [N e1 (x), N e

2 (x)] =

[x− xe2xe1 − xe2

,x− xe1xe2 − xe1

]=

1

le[(xe2 − x), (x− xe1)] (2.25)

suas derivadas em relação a x serão dadas por

Be(x) =d

dxNe(x) =

[dN e

1 (x)

dx,dN e

2 (x)

dx

]=

1

le[−1, 1] (2.26)

e os seus coecientes αe na forma vetorial são dados por

αe = [αe1, αe2] (2.27)

Temos então que a matriz de rigidez de cada elemento é dada por:

Ke =

∫ xe2

xe1

(Be)T p(x)Bedx =

∫ xe2

xe1

1

le

[−11

]p(x)

1

le[−1 1

]dx

=1

(le)2

[−11

] [−1 1

] ∫ xe2

xe1

p(x)dx

=1

(le)2

[1 −1−1 1

] ∫ xe2

xe1

p(x)dx

Ke =1

(le)2

[1 −1−1 1

] ∫ xe2

xe1

p(x)dx (2.28)

Note que, se p(x) = p > 0, ∀x ∈ [0, l] temos que a matriz de rigidez do elemento é constante

Ke =p

le

[1 −1−1 1

](2.29)

Já o vetor de forças externas é composto por duas componentes: fΓe e fΩe . Na primeira com-


ponente, substituindo 2.25 em 2.24, temos que

fΓe =t0le

[(xe2 − x)(x− xe1)

] ∣∣∣∣∣Γn

(2.30)

Observe que se o ponto Γn não pertence ao elemento e, então o vetor fΓe é nulo pois as funçõesteste são denidas de forma que sejam identicamente nulas para todo ponto fora de Ωe.

Já a segunda componente fΩe é calculada por

fΩe =

∫Ωe

(Ne)T f(x)dx =

[ ∫ xe2xe1N1(x)f(x)dx∫ xe2

xe1N2(x)f(x)dx

]=

1

le

[ ∫ xe2xe1

(xe2 − x)f(x)dx∫ xe2xe1

(x− xe1)f(x)dx

](2.31)

Observe que se f(x) = f , ∀x ∈ [xe1, xe2] o vetor fΩe se reduz a:

fΩe =lef

2

[11

]Para obtermos os coecientes α que satisfazem (2.22) devemos construir uma matriz de rigidez

e um vetor de forças externas global, que associe os valores de cada uma das matrizes de rigidez edo vetor de forças externas de cada elemento as entradas da matriz de rigidez e do vetor de forçasexternas global. Para isso, tomando como exemplo os elementos lineares, temos que a matriz derigidez e o vetor de forças externas relativos ao elemento e são da forma:

Ke =

[ke11 ke12

ke21 ke22

]

f e =

[fe1fe2

]onde e = 1, 2, · · · , nel.

Observe que para cada dois elementos consecutivos da malha existe apenas um ponto quepertence a ambos os elementos. Desta forma, é necessário que o valor associado a este ponto namatriz de rigidez global seja igual ao valor associado a este ponto na matriz de rigidez associada aoelemento e somado ao valor obtido na matriz de rigidez (e+1). Neste caso, como estamos utilizando

elementos lineares, este processo é feito somando-se os elementos ke22 e k(e+1)11 na matriz global de

rigidez, com e = 0, 1, · · · , nel − 1. Portanto

K =

k111 k1

12 0 0 0 · · · 0

k121 (k1

22 + k211) k2

12 0 0. . . 0

0 k221 (k2

22 + k311) k3

12 0. . . 0

0 0 k321 (k3

22 + k411) k4

12. . . 0

0 0 0. . . . . . . . . 0

.... . . . . . . . . k

(nel−1)21 (k

(nel−1)22 + knel

11 ) knel12

0 0 0 · · · 0 knel21 knel

22

nt×nt

Observe também que na matriz de rigidez global a entrada kij = 0 se, e somente se, os pontosxi e xj da malha não formam um elemento da malha de estudo.

De maneira análoga, obtemos o vetor global de forças externas somando as entradas fe2 e f (e+1)1 ,


com e = 0, 1, · · · , nel − 1. Desta forma

f =

f11

f12 + f2

1

f22 + f3

1...

f(nel−1)2 + fnel

1

fnel2

nt×1

Desta forma, nosso problema torna-se obter a solução do sistema

Kα = f (2.32)

onde K é a matriz de rigidez global do problema, f é o vetor de forças externas e α é o vetor doscoecientes que satisfazem (2.22). Por outro lado, como conhecemos o valor da função em Γe, temosque o valor de αnt é conhecido. Desta forma, podemos reescrever o sistema linear como

Kα = f (2.33)

onde

K =

k11 · · · k1(nt−1)

k21 · · · k2(nt−1)...

......

k(nt−1)1 · · · k(nt−1)(nt−1)

(nt−1)×(nt−1)

α =

α1

α2...

α(nt−1)

(nt−1)×1

f =

f1 − k1ntαnt

f2 − k2ntαnt

...f(nt−1) − k(nt−1)nt

αnt

(nt−1)×1

2.3.2 Elementos Quadraticos

Assim como para os elementos lineares, iremos particionar Ω em nt pontos igualmente espaça-dos mas agora tomaremos como um elemento Ωe o intervalo obtido por três pontos consecutivos,conforme mostra a gura a seguir:

Figura 2.4: Partição do domínio para elementos quadrático.

Para cada um dos elementos, sendo xe1, xe2 e x

e3 os pontos que compõem o elemento e, temos que


as funções teste são dadas por

Ne(x) = [N e1 (x), N e

2 (x), N e3 (x)] =

[(x− xe2)(x− xe3)

(xe1 − xe2)(xe1 − xe3),

(x− xe1)(x− xe3)

(xe2 − xe1)(xe2 − xe3),

(x− xe1)(x− xe2)

(xe3 − xe2)(xe3 − xe1)

]e, sendo le a distância entre dois pontos, podemos escrever as funções teste como:

Ne(x) =1

2(le)2[(x− xe2)(x− xe3),−2(x− xe1)(x− xe3), (x− xe1)(x− xe2)] (2.34)

Suas derivadas em relação a x serão dadas por

Be(x) =d

dxNe(x) =

1

2(le)2[2x− (xe2 + xe3),−4x+ 2(xe1 + xe3), 2x− (xe1 + xe2)] (2.35)

e os seus coecientes αe, na forma vetorial, são dados por

αe = [αe1, αe2, α

e3] (2.36)

Com isto, temos que a matriz de rigidez do elemento é dada por:

Ke =

∫ xe3

xe1

(Be)T p(x)Bedx =

∫ xe3

xe1

Be1

Be2

Be3

p(x)[Be

1 Be2 Be

3

]dx

=

∫ xe3

xe1

Be1B

e1 Be

1Be2 Be

1Be3

Be2B

e1 Be

2Be2 Be

2Be3

Be3B

e1 Be

3Be2 Be

3Be3

p(x)dx

Assim

Ke =

∫ xe3

xe1

Be1B

e1 Be

1Be2 Be

1Be3

Be2B

e1 Be

2Be2 Be

2Be3

Be3B

e1 Be

3Be2 Be

3Be3

p(x)dx (2.37)

comBe

1 = 12(le)2 (2x− (xe2 + xe3))

Be2 = 1

2(le)2 (−4x+ 2(xe1 + xe3))

Be3 = 1

2(le)2 (2x− (xe1 + xe2))

Já o vetor de forças externas é composto por duas componentes: fΓe e fΩe . Na primeira com-ponente, substituindo (2.34) em (2.24), temos que

fΓe =t0

2(le)2

(x− xe2)(x− xe3)−2(x− xe1)(x− xe3)

(x− xe1)(x− xe2)

∣∣∣∣∣Γn

(2.38)

Novamente temos que o ponto se Γn não pertence ao elemento e o vetor fΓe é nulo pois as funçõesteste foram denidas de forma a serem identicamente nulas para todo ponto fora de Ωe.

A segunda componente fΩe é calculada por

fΩe =

∫Ωe

(Ne)T f(x)dx =

∫ xe3xe1N1(x)f(x)dx∫ xe3

xe1N2(x)f(x)dx∫ xe3

xe1N3(x)f(x)dx

=1

2(le)2

∫ xe3xe1

(x− xe2)(x− xe3)f(x)dx∫ xe3xe1−2(x− xe1)(x− xe3)f(x)dx∫ xe3xe1

(x− xe1)(x− xe2)f(x)dx

(2.39)

Para obtermos os coecientes α que satisfazem (2.22) devemos reunir as matrizes de rigidez eos vetores de forças externas de todos os elementos. Para elementos quadráticos, temos que cada


matriz de rigidez e cada vetor de forças externas são da forma:

Ke =

ke11 ke12 ke13

ke21 ke22 ke23

ke31 ke32 ke33

f e =

fe1fe2fe3

onde e = 1, 2, · · · , nel. Analogamente ao que zemos para elementos lineares, vamos associar osvalores obtidos ao ponto comum a dois elementos consecutivos somando as entradas ke33 e k(e+1)

11 ,para e = 1, 2, · · · , nel − 1 . Desta forma:

K =

k111 k1

12 k113 0 0 0 · · · 0 0

k121 k1

22 k123 0 0 0 · · · 0 0

k131 k1

32 k133 + k2

11 k212 k2

13 0 · · · 0 00 0 k2

21 k222 k2

23 0 · · · 0 00 0 k2

31 k232 k2

33 + k311 k3

22 · · · 0 0...

. . . . . . . . . . . . . . . . . . . . ....

0. . . . . . . . . . . . 0 k

(nel−1)33 + knel

11 knel12 knel

13

0 · · · · · · · · · · · · 0 knel21 knel

22 knel23

0 · · · · · · · · · · · · 0 knel31 knel

31 knel33

(nt×nt)

Novamente, temos que na matriz de rigidez global a entrada kij = 0 se, e somente se, os pontosxi e xj não pertencem a um mesmo elemento. O vetor global de forças é obtido de maneira análoga

ao que zemos para elementos lineares, somando as entradas fe3 e f (e+1)1 , para e = 1, 2, · · · , nel− 1.

Obtemos então o vetor global de forças externas:

f =

f11

f12

f13 + f2

1

f22

f23 + f3

1...

f(nel−1)3 + fnel

1

fnel2

fnel3

(nt×1)


Kα = f

onde K é a matriz de rigidez global do problema, f é o vetor de forças externas e α é o vetor doscoecientes que satisfazem (2.22). Por outro lado, como conhecemos o valor da função em Γe, temosque o valor de αnt é conhecido. Podemos então aplicar a mesma estratégia de resolução utilizadapara elementos lineares para obtermos o sistema equivalente como zemos em (2.33).

2.3.3 Elementos de ordem k

Vamos agora generalizar as idéias utilizadas nos elementos lineares e quadráticos para elementosde ordens maiores. Tomando novamente a malha de pontos igualmente espaçada utilizada nos


exemplos anteriores, consideramos agora como elemento o intervalo contendo k pontos da malha deestudo. Assim:

Figura 2.5: Partição do domínio para elementos de ordem k.

Para cada um dos elementos, sendo xe1, · · · , xek os pontos que compõem o elemento, temos queas funções teste são dadas por

Ne(x) = [N e1 (x), N e

2 (x), · · · , N ek(x)] (2.40)

onde

N ei (x) =

n∏i,j=1;i 6=j

(x− xej)(xei − xej)

e suas derivadas em relação a x serão indicadas por

Be(x) = [dN e

1 (x)

dx,dN e

2 (x)

dx, · · · ,

dN ek(x)

dx] (2.41)

e os coecientes αe, na forma vetorial, são dados por

de = [αe1, αe2, · · · , αek] (2.42)

A matriz de rigidez de cada elemento é dada por

Ke =

∫ xek

xe1

(Be)T p(x)Bedx =

∫ xek

xe1

Be1

· · ·Bek

p(x)[Be

1 · · · Bek

]dx

=

∫ xek

xe1

Be

1Be1 Be

1Be2 · · · Be

1Bek

Be2B

e1 Be

2Be2 · · · Be

2Bek

.... . . . . .

...BekB

e1 Be

kBe2 · · · Be

kBek

p(x)dx

Assim

Ke =

∫ xek

xe1

Be

1Be1 Be

1Be2 · · · Be

1Bek

Be2B

e1 Be

2Be2 · · · Be

2Bek

.... . . . . .

...BekB

e1 Be

kBe2 · · · Be

kBek

p(x)dx (2.43)

com Bej = d

dxNej , j = 1, 2, · · · , k

Assim como para os elementos lineares e quadráticos, o vetor de forças externas é composto por


duas componentes: fΓe e fΩe . Na primeira componente, substituindo (2.40) em (2.24), temos que

fΓe =

N1(x)t0...

Nk(x)t0

∣∣∣∣∣Γn

(2.44)

Novamente, temos que se o ponto Γn não pertence ao elemento e, então o vetor fΓe é nulo, pois asfunções teste são denidas identicamente nulas para todo ponto fora de Ωe.

A segunda componente fΩe é calculada por

fΩe =

∫Ωe

(Ne)T f(x)dx =

∫ xe3xe1N1(x)f(x)dx

...∫ xe3xe1Nk(x)f(x)dx

(2.45)

Para a obtenção dos coecientes αe que satisfazem (2.22), assim como zemos nos elementoslineares e quadráticos, vamos reunir as matrizes de rigidez locais e os vetores de forças externas decada elemento. Sendo as matrizes de rigidez locais e os vetores de forças externas da forma

Ke =

ke11 · · · ke1k...

. . ....

kek1 · · · kekk

f e =

fe1...fek

para e = 1, 2, · · · , nel e utilizando argumento análogo ao utilizado para elementos lineares e qua-dráticos, temos que os valores associados ao ponto pertencente a dois elementos consecutivos serãosomados. Desta forma, temos que:

K =

k111 · · · k1

1k 0 0 0 · · · 0...

. . .... 0 0 0 · · ·

...k1k1 · · · k1

kk + k211 · · · k2

1k 0 · · · 0

0 · · ·...

. . .... 0 · · ·

...0 · · · k2

k1 · · · k2kk + k3

11 0 · · · 0... · · · · · · · · · · · · · · · · · ·

...

0 · · · 0 0 0 k(nel−1)kk + knel

11 · · · k21k

0 · · · 0 0 0...

. . ....

0 · · · 0 0 0 knelk1 · · · knel

kk

(nt×nt)

Assim como nos casos anteriores a matriz de rigidez global a entrada kij = 0 se, e somente se, ospontos xi e xj não pertencem a um mesmo elemento. Assim como foi feito para elementos lineares


e quadráticos, o vetor global de forças é obtido por:

f =

f11...

f1k + f2

1...

f2k + f3

1...

f(nel−1)k + fnel

1...

fnelk

(nt×1)


Kα = f

onde K é a matriz de rigidez global do problema, f é o vetor de forças externas e α é o vetor doscoecientes que satisfazem (2.22). Por outro lado, como conhecemos o valor da função em Γe, temosque o valor de αnt é conhecido. Desta forma, podemos reescrever de maneira análoga ao que foifeito em (2.33).

Capítulo 3

Aspectos computacionais do problema

unidimensional

3.1 Integração Numérica

Na montagem do sistema linear obtido pela formulação do MEF unidimensional temos que cadaentrada da matriz de rigidez é dada pela integral de um produto dos polinômios usados como funçõesteste. Existem diversas técnicas que podem ser utilizadas para o cálculo destas integrais no entanto,devido a sua eciência quando utilizado para avaliar integrandos polinomiais, apresentaremos aquie utilizaremos em nossas simulações os algoritmos da quadratura de Gauss-Legendre.

Nosso objetivo consiste em obter uma aproximação da integral∫ ba f(x)dx da forma:∫ b

af(x)dx ≈

n∑i=1

Aif(xi)

onde xi, i = 1, · · · , n são os chamados nós da quadratura e os Ai, i = 1, · · · , n os respectivos pesosaplicados a cada nó. Nesta seção iremos mostrar como obter tais coecientes, de maneira que aaproximação obtida seja satisfatória.

3.1.1 Quadratura de Gauss-Legendre

Inicialmente, dada uma função f : [a, b] → R contínua, temos que a integral desejada pode serescrita como

I(f) =

∫ b

af(x)dx =

n∑i=1

Aif(xi)︸︷︷︸Q(f)

+E(f)

onde xi ∈ [a, b] são pontos a serem determinados, Ai ∈ R são pesos a serem obtidos e E(f) o errocometido na aproximação.

Nosso objetivo é construir uma fórmula que seja exata (E(f) = 0) para polinômios de maiorgrau possível utilizando o menor número de pontos possível. Observe que:

• Temos 2n incógnitas: x1, · · · , xn, A1, · · · , An

• Impondo que Q(f) seja exata para xk, com 0 ≤ k ≤ (2n − 1), como temos 2n equações,isto sugere que o maior grau a ser alcançado é (2n − 1) pois, para polinômios de grau 2n

19

20 ASPECTOS COMPUTACIONAIS DO PROBLEMA UNIDIMENSIONAL 3.1

considerando f(x) =n∏j=1

(x−xj)2 temos que Q(f) = 0 mas∫ ba f(x)dx > 0 (dado que f(x) > 0,

∀x ∈ R).

Sendo assim, para que exista uma quadradura Q(f) exata para polinômios de grau ≤ (2n− 1)é necessário que:

• Para todo polinômio q de grau menor ou igual a n− 1, vale que∫ b

a(x− x1) · · · (x− xn)︸︷︷︸

p(x)

q(x)dx = Q(f) = 0

Isto implica que p(x) é o polinômio de grau n da família dos polinômios ortogonais mônicosem relação ao produto interno dado por < u, v >=

∫ ba u(x)v(x)dx. É possível mostrar que o

polinômio de grau n da família dos polinômios ortogonais em relação a este produto internopossui exatamente n raízes distintas em (a, b) e que se houver outra família de polinômiosortogonais eles terão as mesmas raízes.

• Sendo Q(f) exata para os polinômios de Lagrange

Li(x) =

n∏j=1,j 6=i

(x− xj)(xi − xj)

e sendo Li(x) de grau (n− 1) < (2n− 1) segue que∫ b

aLi(x)dx = Q(Li) =

n∑j=1

AjLi(xj) = Ai

pois

Li(xj) =

1, se i = j0, se i 6= j

Desta forma, temos que se tal quadratura Q(f) existir, então ela consiste na integral do polinô-mio interpolador da função f em relação aos pontos xi, 1 ≤ i ≤ n raízes do polinômio ortogonalem relação ao produto interno denido pela integral. O teorema a seguir mostra que tal quadraturaexiste de fato e que as condições apresentadas são sucientes para obtê-la.

Teorema 1 (Quadratura Gaussiana). Seja q(x) um polinômio de grau n tal que∫ ba x

kq(x)dx = 0, ∀k ∈ N : 0 ≤ k ≤ n

Sejam xi, i = 1, · · · , n as raízes de q(x). Então, a fórmula∫ b

af(x)dx ≈

n∑i=1

Aif(xi) (3.1)

com

Ai =

∫ b

a

n∏j=1,j 6=i

(x− xjxi − xj

)dx

é exata para todos os polinômios de grau menor ou igual a 2n− 1.

Demonstração: Seja f(x) um polinômio de grau ≤ (2n− 1). Dividindo f(x) por q(x) temos umquociente p(x) e um resto r(x) com grau máximo igual a n− 1. Assim

f(x) = p(x)q(x) + r(x) (3.2)

3.1 INTEGRAÇÃO NUMÉRICA 21

Por hipótese, temos que∫ ba p(x)q(x)dx = 0. Além disso, como xi é raiz de q(x), temos que

f(xi) = p(xi)q(xi)︸︷︷︸q(xi)=0

+r(xi) = r(xi) (3.3)

∀xi, i = 1, · · · , n. Temos então que r(x) é o polinômio interpolador de f em relação aos pontos xi.Sendo r(x) um polinômio de grau menor ou igual a (n− 1), podemos escrevê-lo como

r(x) =

n∑i=1

n∏j=1,j 6=i

(x− xjxi − xj

) r(xi)e, desta forma, temos que

∫ b

ar(x)dx =

∫ b

a

n∑i=1

n∏j=1,j 6=i

(x− xjxi − xj

)r(xi)

︸︷︷︸

r(x)

dx

=n∑i=1

∫ b

a

n∏j=1,j 6=i

(x− xjxi − xj

)dx

︸︷︷︸

=Ai

r(xi) =n∑i=1

Air(xi)

∫ b

ar(x)dx =

n∑i=1

Air(xi) (3.4)

Logo, por 3.2, 3.3 e 3.4, temos que:∫ b

af(x)dx =

∫ b

ap(x)q(x)dx︸︷︷︸

=0

+

∫ b

ar(x)dx =

∫ b

ar(x)dx =

n∑i=1

Air(xi) =

n∑i=1

Aif(xi)

Desta forma, temos que 3.1 será exata para polinômios de grau menor ou igual a 2n− 1. A quadratura de Gauss-Legendre toma como nós do intervalo as raízes dos polinômios de Le-

gendre, os quais possuem a propriedade de serem ortogonais entre si em relação ao produto escalar

< f, g >=

∫ 1

−1f(x)g(x)dx

Tais polinômios podem ser obtidos recursivamente, através da relação:P0(x) = 1P1(x) = x

Pn+2(x) =[

2n+3n+2

]xPn+1(x)−

[n+1n+2

]Pn(x) ∀n ≥ 0

(3.5)

Escolhendo como nós da quadratura as raízes do polinômio de Legendre de graum e substituindona aproximação dada por (3.1) obtemos um sistema onde as m incógnitas serão os coecientes aserem aplicados nos valores da função calculados nos nós escolhidos.

Iremos exemplicar a obtenção dos coecientes para do método para a quadratura com 3 pontos.Os 4 primeiros polinômios de Legendre são dados por:

P0(x) = 1 P1(x) = x P2(x) = 32x

2 − 12 P3(x) = 5x3

2 −3x2

Utilizamos então as raízes do polinômio de grau 3: x0 = −√

0.6, x1 = 0 e x2 =√

0.6 e substituindo


em (3.1) temos:

f(x) = 1⇒∫ 1

−11dx = 2 = A0 +A1 +A2

f(x) = x⇒∫ 1

−1xdx = 0 = −

√0.6A0 + 0A1 +

√0.6A2

f(x) = x2 ⇒∫ 1

−1x2dx =

2

3= (−

√0.6)2A0 + 0A1 + (

√0.6)2A2

e, resolvendo tal sistema linear, obtemos as constantes

A0 = 59 A1 = 8

9 A2 = 59

e a fórmula de integração de Gauss-Legendre de três pontos é dada por∫ 1

−1f(x)dx =

1

9

(5f(−

√0.6) + 8f(0) + 5f(

√0.6)

)(3.6)

cujo valor obtido é exato se f(x) for uma função polinomial de grau ≤ 5, pelo teorema 1. Se afunção f(x) não é polinomial o valor obtido na aproximação representa a integral do polinômiointerpolador de grau ≤ 5 da função f(x) no intervalo [−1, 1].

Na tabela abaixo, apresentamos os valores das raízes dos polinômios de Legendre e suas res-pectivas ponderações, que podem ser utilizadas para obtermos uma maior precisão no cálculo dasintegrais. Para outros valores, o leitor pode consultar http://dlmf.nist.gov/3.5#v.1

n Nós xi Pesos Ai2 -0.5773502691896257645091488 1

0.5773502691896257645091488 13 -0.774596669241483377035853079 0.555555555555555555556

0 0.88888888888888888890.774596669241483377035853079 0.55555555555555555555556

4 -0.3399810435848562648026658 0.65214515486254614262693610.3399810435848562648026658 0.6521451548625461426269361-0.8611363115940525752239465 0.34785484513745385737306390.8611363115940525752239465 0.3478548451374538573730639

Tabela 3.1: Nós da quadratura de Gauss-Legendre e suas ponderações.

Embora as deduções tenham sido realizadas para a integração de uma dada função no intervalo[−1, 1], podemos utilizá-las em qualquer intervalo genérico [a, b] desde que realizemos uma mudançade variáveis. Seja

x =2t− (a+ b)

(b− a)⇒ t =

1

2[(b− a)x+ (b+ a)]

Desta forma, a integral de f(t) em [a, b] é transformada em uma integral sobre o domínio [−1, 1]da forma:∫ b

af(t)dt =

(b− a)

2

∫ 1

−1f

((b− a)

2x+

(b+ a)

2

)dx ≈ (b− a)

2

n∑i=0

Aif

((b− a)

2x+

(b+ a)

2

)(3.7)

com os valores dos Ai obtidos pela resolução do sistema linear gerado na integração dos polinômiosde Legendre.

1Os valores a serem utilizados na quadratura numérica são apresentados com um grande número de dígitos devidoao fato de que a quadratura é exata apenas para os valores exatos dos nós e dos pesos. Porém, como tais valores sãoirracionais, temos que o valor aproximado deve ser fornecido com precisão superior aquela utilizada nos cálculos.

http://dlmf.nist.gov/3.5#v

3.3 ARMAZENAMENTO DO SISTEMA LINEAR 23

3.2 Armazenamento do Sistema Linear

Vimos anteriormente que a construção da matriz global de rigidez do MEF unidimensional éfeita associando-se as matrizes de rigidez obtidas para cada um dos elementos do domínio numamatriz global de rigidez. Como as matrizes de cada elemento são simétricas, temos que a matrizglobal gerada também será simétrica. Por outro lado, a interseção de dois elementos diferentespossui apenas um ponto em comum, fazendo com que a matriz de rigidez global apresente estruturade banda, isto é, sendo K = kij a matriz de rigidez global, vale que ∀i, j ∈ [1, nt], se |i− j| > L⇒kij = 0, com L ∈ N dependendo do número de pontos n utilizados em cada elemento para deniras funções teste. A Figura 3.1 mostra uma matriz com a propriedade de estrutura de banda.

Figura 3.1: Matriz com propriedade de banda.

Além disso, devido a simetria da matriz, o armazenamento das entradas pode ser feito com o usode apenas n vetores. Tomando como exemplo a matriz de rigidez obtida de um domínio com doiselementos (nel = 2) cujas funções teste são obtidas com o uso de três pontos do elemento (n = 3,nt = 5), temos que

K =

k1

11 k112 k1

13 0 0k1

21 k122 k1

23 0 0k1

31 k132 k1

33 + k211 k2

12 k213

0 0 k221 k2

22 k223

0 0 k231 k2

32 k233

pode ser armazenada como

dj = [ki,i]j ; 1 ≤ i ≤ nt; 1 ≤ j ≤ ntaj = [ki,i+1]j ; 1 ≤ i ≤ (nt − 1); 1 ≤ j ≤ (nt − 1)

bj = [ki,i+2]j ; 1 ≤ i ≤ (nt − 2); 1 ≤ j ≤ (nt − 2)

onde os vetores d, a e b armazenam os valores da diagonal principal, da primeira diagonal acimada principal e da segunda diagonal acima da principal, respectivamente.

Observe que no caso geral para nel ≥ 2, tomando uma matriz com k pontos em Ω, para oarmazenamento da matriz completa seriam necessárias k2 posições de memória onde, utilizandoelementos quadráticos, apenas (3k − 3) necessitariam ser armazenadas.2

2Nesta estimativa ainda sim estamos armazenando alguns zeros nos vetores acima da diagonal principal. Isto éfeito para evitar problemas com os índices a serem utilizados nos cálculos a serem realizados com a matriz.


3.3 Resolução de sistemas lineares com estrutura de banda

Apresentaremos nesta seção os algoritmos para a resolução do sistema linear obtido para ele-mentos com dois e três nós. Antes de desenvolvermos os métodos de solução do sistema, observamosalguns fatos importantes sobre o sistema linear (2.33):

• A matriz do sistema possui a propriedade de ser diagonal dominante, isto é

|kii| ≥n∑j=1

|kij | (3.8)

• A matriz é inversível.

• A matriz do sistema possui estrutura de banda.

Devido as duas primeira propriedades, temos que a aplicação do processo de eliminação deGauss não requer troca de linhas. A partir deste fato e devido a terceira propriedade é possívelmostrar que, sendo U a matriz triangular superior obtida da eliminação de Gauss e L a matriztriangular inferior formada pelos multiplicadores, temos que U tem banda superior igual a bandasuperior da matriz do sistema e L tem banda inferior igual a banda inferior da matriz do sistema.Isto possibilita que obtenhamos algoritmos ecazes para resolver o sistema, conforme observamos aseguir.

3.3.1 Elementos com dois nós

Com o uso de elementos denidos sobre dois nós do intervalo, vimos anteriormente que a matrizde rigidez do sistema linear é dada por uma matriz tridiagonal. Utilizando -se do armazenamentoapresentado vetorial apresentado na seção anterior, temos que a resolução do sistema pode ser feitade maneira bastante eciente.

Para isto, partindo da matriz tridiagonal:

d1 a1 0 0 · · · 0a1 d2 a2 0 · · · 00 a2 d3 a3 · · · 0...

. . . . . . . . . . . ....

0 · · · 0 a(nt−2) d(nt−1) a(nt−1)

0 · · · 0 0 a(nt−1) dnt

x1

x2......

x(nt−1)

xnt

=

f1

f2......

f(nt−1)

fnt

(3.9)

eliminamos os valores abaixo da diagonal principal fazendo:

di ← di −(a(i−1)

d(i−1)

)ai

fi ← fi −(a(i−1)

d(i−1)

)fi

3.3 RESOLUÇÃO DE SISTEMAS LINEARES COM ESTRUTURA DE BANDA 25

para 2 ≤ i ≤ nt. Desta forma, o sistema resultante será da forma

d1 a1 0 0 · · · 00 d2 a2 0 · · · 00 0 d3 a3 · · · 0...

. . . . . . . . . . . ....

0 · · · 0 0 d(nt−1) a(nt−1)

0 · · · 0 0 0 dnt

x1

x2......

x(nt−1)

xnt

=

f1

f2......

f(nt−1)

fnt

A partir daí utilizamos a substituição inversa para obtermos a solução do sistema linear:

xnt =fnt

dnt

xi = (fi − cixi+1)/di i = (nt − 1), · · · , 1

obtendo assim os valores de x que satisfazem (3.9). Tal processo encontra-se descrito no Algoritmo 1:

Algoritmo 1: Resolução de sistema tridiagonal. [CK08]Entrada: Vetores d1×n,a1×(n−1), f1×nSaída: Vetor x1×n de soluções do sistemacte← 0;para i = 2 até nt faça

cte← ai−1/di−1;di ← di − (cte)ai−1;fi ← fi − (cte)fi−1;

m para

xnt ← fnt/dn;para i = (nt − 1) até 1 faça

xi ← (fi − cixi+1)/di;m para

3.3.2 Elementos com três nós

Com o uso de elementos denidos sobre tês nós do intervalo, vimos que a matriz de rigidezdo sistema linear é dada por uma matriz pentadiagonal. Utilizando o armazenamento apresentadovetorial apresentado anteriormente, temos que a resolução do sistema é ser feita de maneira bastanteeciente, utilizando-se de um método semelhante ao aplicado para elementos denidos sobre doisnós.

Partindo do sistema

d1 a1 b1 0 · · · · · · 0a1 d2 a2 b2 0 · · · 0b1 a2 d3 a3 b3 · · · 00 b2 a3 d4 a4 · · · 0...

. . . . . . . . . . . . . . ....

0 · · · b(nt−4) a(nt−3) d(nt−2) a(nt−1) b(nt−2)

0 · · · 0 b(nt−3) a(nt−2) d(nt−1) a(nt−1)

0 · · · 0 0 b(nt−2) a(nt−1) dnt

x1

x2.........

x(nt−1)

xnt

=

f1

f2.........

f(nt−1)

fnt

(3.10)

e utilizando uma idéia análoga á utilizada para elementos com dois nós, zeramos as entradas abaixo


da diagonal principal, obtendo:

d1 a1 b1 0 · · · · · · 00 d2 a2 b2 0 · · · 00 0 d3 a3 b3 · · · 00 0 0 d4 a4 · · · 0...

. . . . . . . . . . . . . . ....

0 · · · 0 0 d(nt−2) a(nt−1) b(nt−2)

0 · · · 0 0 0 d(nt−1) a(nt−1)

0 · · · 0 0 0 0 dnt

x1

x2.........

x(nt−1)

xnt

=

f1

f2.........

f(nt−1)

fnt

Em seguida, novamente utilizando a substituição inversa, obtemos a solução do sistema linear. Talprocesso encontra-se descrito no Algoritmo 2:

Algoritmo 2: Resolução de sistema pentadiagonal. [CK08]Entrada: Vetores d1×n,a1×(n−1),b1×(n−2), f1×nSaída: Vetor x1×n de soluções do sistemacte← 0;r ← a1;s← a2;t← b1;para i = 2 até nt − 1 faça

cte← r/di−1;di ← di − (cte)ai−1;ai ← ai − (cte)fi−1;fi ← fi − (cte)fi−1;cte← t/d(i−1);r ← s− (cte)a(i−1);di+1 ← d(i+1) − (cte)f(i−1);fi+1 ← f(i+1) − (cte)f(i−1);s← a(i+1);t← bi;

m para

cte← r/d(nt−1);dnt ← dnt − (cte)a(nt−1);xnt ← (fnt − (cte)f(nt−1))/dnt ;x(nt−1) ← (fnt−1 − a(nt−1)xnt)/dnt−1;para i = (nt − 2) até 1 faça

xi ← (fi − bix(i+2) − aix(i+1))/di;m para

3.3.3 Elementos com k nós

Podemos aplicar uma estratégia análoga para a resolução do sistema linear gerado por elementosutilizando k nós. Novamente, partindo da matriz inicial do sistema, eliminamos os valores abaixo dadiagonal principal e em seguiida utilizamos substituição inversa para obtermos a solução do sistemalinear. Não iremos apresentar aqui o algorimo de resolução para este tipo de sistema pois o mesmoé análogo aos Algoritmos 1 e 2.

3.4 CONVERGÊNCIA NUMÉRICA DO ALGORITMO 27

3.4 Convergência numérica do algoritmo

Após a implementação dos algoritmos, gostaríamos de obter uma estimativa do comportamentodo erro em relação ao espaçamento h entre os nós do domínio.

A teoria de convergência do método mostra que o erro obtido na aproximação com espaçamentoh é da ordem de Chp, com p dependendo do tipos de elementos escolhidos para a construção dasfunções teste e C uma constante que independe do renamento da malha. Para elementos lineares,temos que a ordem esperada do método é p = 2, enquanto que para elementos quadráticos, temosque a ordem esperada é p = 3.

Para analisarmos o comportamento do erro no algoritmo implementado iremos considerar duasmalhas igualmente espaçadas com 2p e 2p+1 elementos, com p ∈ N. Desta forma, temos que

erro(2p)

erro(2p+1)∼=

Chp

C(h2 )p+1= 2p

Chp

Chp⇒ p ∼= Log2

(erro(2p)

erro(2p+1)

)(3.11)

Assim, analisando o Log2 do quociente entre os erros obtidos para aproximações obtidas pormalhas renadas pela metade do espaçamento anterior é possível analisar se ocorre a convergênciados valores obtidos de acordo com a ordem esperada conforme aumentamos o número de elementos(e consequentemente o número de pontos na malha) do problema.


Capítulo 4

Simulações unidimensionais

Nesta seção apresentaremos a implementação computacional do Método dos Elementos Finitose os resultados obtidos na aplicação do método para obtenção da aproximação numéricas de umproblema estacionário de condução de calor apresentado em [FB07].

4.1 Apresentação do Problema

Considere uma barra com uma fonte de calor uniformemente distribuída de f = 5Wm−1. Abarra possui uma área de seção transversal uniforme de A = 0.1m2 e condutividade térmica dek = 2W C−1m−1. O comprimento da barra é de 4m. As condições de contorno são T (0) = 0C e

−p(x)dTdx

∣∣∣x=4

= q(x = 4) = 5Wm−1, conforme mostrado na Figura 4.1.

Figura 4.1: Visualização do problema

Com os dados acima, temos que a formulação forte do problema é dada por:

• Obter T : [0, 4]→ R de classe C2 que satisfaça as seguintes condições:

ddx

(0.2dTdx

)= −5, ∀x ∈ (0, 4)

0.2dTdx∣∣x=4

= −Aq(x = 4) = −0.5;

T (x = 0) = 0;

(4.1)

Observe que é possível encontrarmos a solução analítica do problema acima pela integração daequação diferencial e a aplicação das condições de contorno. Desta forma, obtemos que a soluçãoexata do problema é dada por:

T (x) = −12.5x2 + 97.5x (4.2)

Para aplicarmos o MEF inicialmente devemos obter a formulação fraca do problema em estudo.Seguindo os passos apresentados no capítulo 2, temos que a formulação fraca em estudo é dada por:

29

30 SIMULAÇÕES UNIDIMENSIONAIS 4.1

• Obter u : [0, 4]→ R , de classe C2 com u(0) = 0 que satisfaça:∫ 4

0

(dw

dx

)(0.2

dT

dx

)dx =

∫ 4

05w(x)dx+ 0.5w(4),∀w : w(0) = 0 (4.3)

4.1.1 Resolução com o uso de Elementos Lineares

Seguindo os passos apresentados no Capítulo 2, como p(x) e f(x) são constantes, utilizando amalha de pontos igualmente espaçada temos que para elementos lineares as matrizes de rigidez eos vetores de força em cada elemento são constantes e dados, respectivamente, por:

Ke =0.2

l

[1 −1−1 1

]

fΩe =5l

2

[11

]

fΓe = −0.5× (Ne)T∣∣∣Γn

=

00−0.5

se x = 4 ∈ Γe. 000

caso contrário.

para e = 1, 2, · · · , nel e l = 4/(nt − 1).Desta forma, utilizando a estratégia apresentada na seção (3.2) e o Algoritmo 1 obtemos os

valores em cada um dos nós da malha e, por interpolação, obtemos as aproximações nos pontos quenão pertencem à malha. A Figura 4.2 fornece a visualização da solução exata e das aproximaçõesobtidas para diferentes quantidades de elementos.

Figura 4.2: Aproximações obtidas com elementos lineares para 2, 4, 8 e 16 elementos. A solução exata édada pela curva azul e as aproximações são dadas pelas curvas vermelhas.

4.1 APRESENTAÇÃO DO PROBLEMA 31

Já a Tabela 4.1 apresenta a relação entre a quantidade de elementos no qual o domínio foidividido, a quantidade de pontos e o erro máximo obtido na comparação entre a aproximaçãoobtida e a solução exata em 50 pontos igualmente espaçados do domínio.

Qtd Elementos Qtd Pontos Erro máximo2 3 12.49474 5 3.12368 9 0.780916 17 0.195232 33 0.048864 65 0.0122128 129 0.00305256 257 0.0007

Tabela 4.1: Erros obtidos na simulação utilizando elementos lineares.

Observe que a cada iteração o erro obtido cai aproximadamente por um fator 0.25. Desta forma,temos que a ordem do método, conforme a estratégia apresentada em (3.11), é igual a 2 coincidindocom o valor apresentado em teoria para a aplicação do método com o uso de elementos lineares.

4.1.2 Resolução com o uso de Elementos Quadraticos

Para a resolução do problema em estudo com elementos quadráticos, é esperado que obtenhamosa solução exata do problema, a menos de pequenos erros de arredondamento. Isto ocorre devido aofato de que a solução do problema é um polinômio da mesma ordem das aproximações que serãoobtidas em cada um dos elementos (isto é ambos são polinômios quadráticos).

Após a realização das simulações e pela análise dos grácos da Figura (4.5) observamos que osvalores obtidos são exatos, a menos de um pequeno erro cometido durante os cálculos em aritméticade ponto utuante. Para a observação do erro cometido, temos que a Figura 4.4 apresenta um fatorde escala de 1e+ 12.

Figura 4.3: A direita temos o gráco da aproximação obtida e a esquerda o gráco da solução exata.


Figura 4.4: Erro cometido na aproximação com fator de escala 1e+12.

Optamos então por substituir a fonte de calor constante por uma fonte de calor dada porf(x) = 0.3125π2sen(πx8 )Wm−1 com as novas condições de contorno dadas por T (0) = 0C e

−p(x)dTdx

∣∣∣x=4

= q(x = 4) = 0Wm−1

Com os novos dados, temos que a formulação forte do problema é dada por:

• Obter T : [0, 4]→ R de classe C2 que satisfaça as seguintes condições:

ddx

(0.2dTdx

)= −0.3125π2sen

(πx8

), ∀x ∈ (0, 4)

0.2dTdx∣∣x=4

= −Aq(x = 4) = 0;

T (x = 0) = 0;

(4.4)

Novamente, temos que é possível obtermos a solução analítica do problema acima pela integraçãoda equação diferencial e a aplicação das condições de contorno. Para este caso a solução exata doproblema é dada por:

T (x) = 100sen(πx

8

)(4.5)

Com os novos dados a formulação fraca em estudo é dada por:

• Obter u : [0, 4]→ R , de classe C2 com u(0) = 0 que satisfaça:∫ 4

0

(dw

dx

)(0.2

dT

dx

)dx =

∫ 4

00.3125π2sen

(πx8

)w(x)dx, ∀w : w(0) = 0 (4.6)

Seguindo os passos apresentados no Capítulo 2 e utilizando uma malha igualmente espaçadageramos as matrizes de rigidez de cada elemento e seus respectivos vetores de forças externas, comodescrito em (2.37), (2.38) e (2.39). Em seguida, geramos o sistema global do problema (2.33) eutilizando a estratégia apresentada na seção (3.2) e o Algoritmo 2 obtivemos os valores em cadaum dos nós da malha e, por interpolação, obtemos as aproximações nos pontos que não pertencemà malha. A Figura 4.5 apresenta as aproximações obtidas.

4.1 APRESENTAÇÃO DO PROBLEMA 33

Figura 4.5: As curvas em vermelho apresentam as aproximações numéricas e as curvas em azul os valoresda solução do problema.

Por m, temos que a Tabela 4.2 apresenta a relação entre a quantidade de elementos no qualo domínio foi dividido, a quantidade de pontos e o erro máximo obtido na comparação entre aaproximação obtida e a solução exata em 50 pontos igualmente espaçados do domínio.

Qtd Elementos Qtd Pontos Erro máximo2 5 0.3585344 9 0.04667998 17 0.005946516 33 0.000679332 65 0.000091264 129 0.0000117128 257 0.0000015256 513 0.0000002

Tabela 4.2: Erros obtidos na simulação utilizando elementos quadráticos.

Neste caso a cada iteração o erro obtido cai aproximadamente por um fator 0.125. Utilizandonovamente a estratégia apresentada em (3.11), temos que a ordem do método utilizando-se deelementos quadráticos é igual a 3.


Capítulo 5

Introdução ao Método dos Elementos

Finitos bidimensional

Assim como na aplicação do MEF para problemas unidimensionais, inicialmente iremos desenvol-ver a formulação fraca associada ao problema, particionar o problema em uma malha de elementose, sob cada um destes elementos, iremos denir as funções teste que satisfaçam a formulação fraca.

Nesse contexto, diferentemente do que ocorria em problemas unidimensionais em que o domínioera um segmento de reta, temos agora que o domínio do problema é um subconjunto do R2 e, destaforma, a escolha dos elementos a serem utilizados depende da geometria do domínio do problema.Devido ao cárater introdutório deste trabalho, iremos enfatizar a discretização do domínio atravésde elementos triangulares de lados retos.

Neste capítulo iremos discutir como obter a formulação fraca do problema bidimensional asso-ciado, como obter as funções teste sobre elementos triangulares de três nós, como gerar as matrizesde rigidez, os vetores de força do problema e o sistema linear associado ao problema.

5.1 Formulações do Problema

Para os modelos bidimensionais analisados neste capítulo, diremos que um problema encontra-seem sua formulação forte quando a equação diferencial obtida é descrita pelo seguinte problema:

- Dado Ω ⊂ R2 um conjunto aberto e conexo, obter u : Ω → R, de classe C2, que satisfaça asseguintes condições:

∇ · (D∇u) = f(x, y), ∀(x, y) ∈ Ω (5.1)

D∇u · ~n = g(x, y), ∀(x, y) ∈ Γn (5.2)

u(x, y) = h(x, y), ∀(x, y) ∈ Γe (5.3)

onde D(x, y) uma matriz simétrica e positiva denida, f(x, y) uma função suave denida em Ω ∈R2, Γe e Γn são as condições essenciais e naturais do problema, respectivamente. Assim comopara problemas unidimensionais, temos que tais condições são subconjuntos de Γ = ∂Ω tais queΓ = Γe ∪ Γn e Γe ∩ Γn = ∅ .

Observe que com P = I2×2 temos que a equação (5.1) é dada por:

∇ · (∇u) =∂2u

∂x2+∂2u

∂y2= ∆u = f(x, y) (5.4)

que é a forma canônica de uma equação diferencial elíptica. Seguindo os objetivos deste trabalho,iremos enfatizar a aplicação do método para equações diferenciais elípticas. No entanto, ressaltamos

35

36 INTRODUÇÃO AO MÉTODO DOS ELEMENTOS FINITOS BIDIMENSIONAL 5.1

aqui que o método também é aplicável à equações diferenciais parabólicas e hiperbólicas. O leitorinteressado pode encontrar mais informações em [Goc06], [DH03] e [Bat96]. A obra [Bra06] apresentaa teoria com um enfoque matemático bastante interessante.

Assim como zemos nos problemas unidimensionais, a partir da formulação forte, iremos desen-volver a formulação fraca equivalente. Para isso, utilizaremos dois teoremas bastante conhecidos docálculo vetorial:

Teorema 2 (Teorema de Green). Seja f : Ω ⊂ R2 → R, contínua e integrável. Sendo Γ = ∂Ω umacurva seccionalmente de classe C1 orientada no sentido anti-horário, vale que∫ ∫

Ω∇fdΩ =

∮Γf~ndΓ (5.5)

onde ~n é a normal exterior unitária.

Teorema 3 (Teorema da divergência de Gauss). Seja f : Ω ⊂ R2 → R, de classe C2 e integrável.Sendo Γ = ∂Ω uma curva seccionalmente de classe C1 orientada no sentido anti-horário, vale que∫ ∫

Ω∇ · (∇f)dΩ =

∫ ∫ΩdivfdΩ =

∮Γ∇f · ~ndΓ (5.6)

onde ~n é a normal exterior unitária.

Não faremos aqui a demonstração destes resultados. Os leitores interessados podem encontrá-lasindicadas em [Gui02]. Observe que a equação (5.5) é similar ao resultado obtido na integração porpartes em (2.6): aqui o operador d

dx é equivalente ao operador gradiente.Utilizando os dois teoremas acima, iremos apresentar agora uma propriedade que será funda-

mental na obtenção da formulação fraca:

Teorema 4 (Fórmula de Green). Dadas f : Ω ⊂ R2 → R, de classe C2, w : Ω ∈ R2 → R de classeC1 temos que∫ ∫

Ωw∆fdΩ =

∫ ∫Ωw∇ · (∇f)dΩ =

∮w∇f · ~ndΓ−

∫ ∫Ω∇wT · ∇fdΩ (5.7)

onde ~n é a normal exterior unitária e Γ = ∂Ω é uma curva seccionalmente de classe C1 orientadano sentido anti-horário.

Demonstração: Para demonstrarmos a fórmula de Green, inicialmente iremos avaliar ∇· (w∇f)como uma derivada da regra do produto:

∇ · (w∇f) =∂

∂x

(w∂f

∂x

)+

∂

∂y

(w∂f

∂y

)=∂w

∂x

∂f

∂x+ w

∂2f

∂x2+∂w

∂y

∂w

∂y+ w

∂2f

∂y2

= w

(∂2f

∂x2+∂2f

∂y2

)︸︷︷︸∇·(∇f)=∆f

+

(∂w

∂x

∂f

∂x+∂w

∂y

∂w

∂y

)︸︷︷︸

∇wT ·∇f

= w∇ · (∇f) +∇wT · ∇f

⇒ ∇ · (w∇f) = w∇ · (∇f) +∇wT · ∇f (5.8)

Integrando (5.8) sobre o domínio, temos∫ ∫Ω∇ · (w∇f)dΩ =

∫ ∫Ωw∇ · (∇f)dΩ +

∫ ∫Ω∇wT · ∇fdΩ (5.9)

5.1 FORMULAÇÕES DO PROBLEMA 37

Observe que podemos aplicar o teorema da divergência no primeiro termo de (5.9):∫ ∫Ω∇(w∇f)dΩ =

∮Γw∇f · ~ndΓ

Desta forma, reorganizando os termos, obtemos a fórmula de Green (5.7).De posse da fórmula de Green, obtemos a formulação fraca do problema bidimensional seguindo

os mesmos passos utilizados na formulação fraca unidimensional. Escolhemos, de forma arbitrária,uma função de ponderação w : Ω ∈ R2 → R, de classe C1, tal que w(x, y) = 0, ∀(x, y) ∈ Γe (noteque, tal condição equivale a condição da função de ponderação unidimensional ser nula em x = l(2.3)). Para simplicarmos a notação, considere ~q(x, y) = D∇u. Desta forma, multiplicando asequações (5.1) e (5.2) pela função de ponderação e integrando sob os respectivos domínios, obtemos∫ ∫

Ωw∇~qdΩ =

∫ ∫ΩwfdΩ (5.10)

∫Γn

w~qdΓn =

∫Γn

t0~qdΓn (5.11)

Aplicando a fórmula de Green (5.7) ao termo da esquerda em (5.10) , temos∫ ∫Ωw∇~qdΩ =

∮Γ(w~q) · ~ndΓ−

∫ ∫Ω∇wT · ~qdΩ (5.12)

que, após substituída em (5.10) e reorganizados os termos, resulta em∫ ∫Ω∇wT · ~qdΩ =

∮Γ(w~q) · ~ndΓ−

∫ ∫ΩwfdΩ (5.13)

Como Γ = Γe ∪ Γn e Γe ∩ Γn = ∅ , podemos reescrever (5.13) como∫ ∫Ω∇wT · ~qdΩ =

∫Γn

(w~q) · ~ndΓ +

∫Γe

(w~q)~ndΓ−∫ ∫

ΩwfdΩ (5.14)

Porém, pela denição da função de ponderação, temos que w = 0 em Γe. Além disso, utilizando(5.11), temos que (5.14) pode ser reescrita como:∫ ∫

Ω∇wT · ~qdΩ =

∫Γn

wt0~ndΓ−∫ ∫

ΩwfdΩ (5.15)

Desta forma, dizemos que um problema bidimensional encontra-se em sua formulação fracaquando é apresentado da seguinte forma:

- Obter u : Ω ∈ R2 → R , de classe C2, com u(x, y) = h(x, y), ∀(x, y) ∈ Γe que satisfaça:∫ ∫Ω∇wT · ~qdΩ =

∫Γn

wgdΓ−∫ ∫

ΩwfdΩ (5.16)

∀w : w(x, y) = 0, ∀(x, y) ∈ Γe, onde w(x, y) : Ω→ R de classe C1.Para mostrarmos a equivalência entre as formulações do problema utilizaremos a mesma estra-

tégia utilizada na demonstração da equivalência dos problemas unidimensionais onde revertemos ospassos da fórmula de Green e usamos a arbitrariedade da função de ponderação para obtermos asequações da formulação forte.

Partindo de (5.16) e aplicando a fórmula de Green para o primeiro termo, temos∫ ∫Ωw(∇~q − f)dΩ +

∫Γn

w(t0 − ~q~n)dΓ−∫

Γe

w~q~ndΓ = 0 (5.17)


∀w(x, y) de classe C2 tal que w(x, y) = 0, ∀(x, y) ∈ Γe. Usando da arbitrariedade da função peso,podemos escolhê-la de forma que

w = ψ(x, y)(∇~q − f) (5.18)

onde ψ(x, y) é uma função sucientemente suave tal que

ψ(x, y) =

[0, ∀(x, y) ∈ Γe> 0, ∀(x, y) ∈ Ω

]Desta forma, substituindo (5.18) em (5.17) obtemos∫ ∫

Ωψ(∇~q − f)2 = 0 (5.19)

Tal igualdade ocorre somente se a condição (5.1) for satisfeita, devido a nossa escolha de w(x, y).Sendo satisfeita tal condição, segue de (5.17) que∫

Γn

ψ(t0 − ~q · ~n)2dΓ = 0 (5.20)

Sendo o integrando em (5.20) positivo, temos que a função entre parênteses deve ser nula emtodo ponto do contorno natural do problema, reestabelendo a condição (5.2). Desta forma, temosreestabelecida a formulação forte do problema.

Observamos então que a suavidade das funções envolvidas é extremamente importante tambémpara os problemas bidimensionais. Assim como em problemas unidimensionais, a formulação fracapermite relaxar as hipóteses sobre a solução do problema podendo novamente ser vista como umaformulação mais geral do problema.

5.2 Tipos de elementos e aproximações

Um fator importante a ser destacado neste momento é a geometria do problema em estudo.No caso unidimensional, o domínio Ω consiste em um intervalo de R, agora temos que o Ω é umasuperfície do R2. Desta forma, existe certa liberdade sobre a escolha dos elementos: triângulos,retângulos, trapézios e quaisquer outras formas geométricas simples (e até algumas cujos lados nãosejam segmentos de reta) podem ser escolhidas como elementos para a malha do problema.

Quando trabalhávamos com elementos unidimensionais, tínhamos que dois elementos consecu-tivos possuiam apenas um ponto em comum. Para problemas bidimensionais, vamos impor que oselementos adjacentes possuam, no máximo, uma aresta em comum (o segmento que os une). AFigura 5.1 exemplica isto usando dois elementos triangulares.

Figura 5.1: Dois elementos quaisquer possuem em comum apenas o segmento que os une.


Devido a liberdade na escolha dos elementos, o número de pontos em cada elemento variaconforme o elemento escolhido. Porém, assim como é feito nos problemas unidimensionais, nossoobjetivo é obter uma aproximação do problema (5.1)-(5.2)-(5.3) sobre cada um dos elementos pormeio da interpolação de funções teste denidas sobre cada um dos elementos.

Iremos apresentar em detalhes a aplicação do elemento triangular com três nós por se tratarde um dos mais versáteis e simples elementos em duas dimensões. Podemos representar facilmentequase toda geometria com elementos triangulares e, sem muitas diculdadades, construir malhasque tenham mais elementos em áreas onde estejamos interessados em uma maior precisão. Outravantagem é o fato de que os geradores de malhas triangulares são bastante robustos. Ao leitorinteressado, as referências [FB07], [Bat96], [Goc06] e [DH03] apresentam como utilizar outros tiposde elementos.

5.2.1 Elemento Triangular com três nós

Um modelo de elemento triangular com três nós é apresentado na Figura 5.2. Conforme podemosobservar, em cada um dos elementos, os nós do elemento são posicionados nas extremidades dotriângulo e denotados por (xei , y

ei ), com i = 1, 2, 3. A numeração dos nós é feita para cada elemento

da malha e, neste trabalho, será feita no sentido anti-horário, pois esta é a convenção adotada porgrande parte dos livros na literatura e por programas de elementos nitos.

Figura 5.2: Elemento triangular de três nós

Devido a natureza do elemento, temos que um número arbitrário de elementos podem ser conec-tados a um mesmo nó. Não há restrições sobre a topologia de uma malha de elementos triangulares,embora seja aconselhável que nenhum dos ângulos de qualquer elemento seja muito agudo.

Em elementos deste tipo, as aproximações a serem obtidas serão da forma

ue(x, y) = αe0 + αe1x+ αe2y (5.21)

com ue(x, y) = 0, se (x, y) /∈ Ωe e os coecientes αei são obtidos de maneira que a aproximaçãoobtida sobre as arestas comuns a dois elementos sejam iguais, isto é,

u1(x, y) = u2(x, y) (5.22)

∀(x, y) ∈ s, conforme observamos na Figura 5.3.A construção das funções teste sobre os elementos triangulares com três nós será feita utilizando

as mesmas idéias utilizadas nas funções teste unidimensionais: Vamos obter N1(x, y), N2(x, y) eN3(x, y) tais que

Nj(xek, y

ek) = δjk (5.23)

onde

δjk =

1, j = k0, j 6= k

com j, k = 1, 2, 3. Observe que isto é possível pois temos três nós conhecidos e três coecientes a


Figura 5.3: Continuidade entre elementos adjacentes

serem obtidos, de maneira que os constantes obtidas serão únicas. Temos então que: 100

=

1 xej yej1 xek yek1 xel yel

αe0αe1αe2

(5.24)

com k 6= l 6= k e j, k, l = 1, 2, 3. Resolvendo este sistema pela Regra de Cramer, temos que asfunções teste são dadas por

N ej (x, y) =

Dkl(x, y)

2Ae(5.25)

onde Ae é a área do elemento triangular e

Dkl(x, y) = det

1 x y1 xek yek1 xel yel

Sendo Ωe o conjunto dos pontos (x, y) ∈ R2 que compõem o elemento e, obtemos que as funções

teste para elementos triangulares com três nós são dadas por:

N e(x, y) =

[N e

1 (x, y) N e2 (x, y) N e

3 (x, y)], (x, y) ∈ Ωe

0, (x, y) /∈ Ωe (5.26)

onde

N e1 (x, y) =

1

2Ae((xe2y

e3 − xe3ye2) + (ye23)x+ (xe32)y)

N e2 (x, y) =

1

2Ae((xe3y

e1 − xe1ye3) + (ye31)x+ (xe13)y)

N e3 (x, y) =

1

2Ae((xe1y

e2 − xe2ye1) + (ye12)x+ (xe21)y)

Ae é a área do elemento em estudo, xeij = xei − xej e yeij = yei − yej .Podemos também representar as funções teste na forma matricial

Ne(x, y) =1

2Ae

(xe2ye3 − xe3ye2) (ye23) (xe32)

(xe3ye1 − xe1ye3) (ye31) (xe13)

(xe1ye2 − xe2ye1) (ye12) (xe21)

︸︷︷︸

Me

1xy

(5.27)

Temos então que a aproximação em cada um dos elementos é dada por

ue(x, y) = Ne(x, y)(αe)T


e o gradiente da aproximação em cada um dos elementos é dado por:

∇ue(x, y) =

[∂ue(x,y)

∂x∂ue(x,y)

∂y

]=

[∂Ne

1∂x α

e1 +

∂Ne2

∂x αe2 +

∂Ne3

∂x αe3

∂Ne1

∂y αe1 +

∂Ne2

∂y αe2 +

∂Ne3

∂y αe3

]=

[∂Ne

1∂x

∂Ne2

∂x∂Ne

3∂x

∂Ne1

∂y∂Ne

2∂y

∂Ne3

∂y

]︸︷︷︸

Be

αe1αe2αe3

︸︷︷︸

(αe)T

= Be(αe)T

Sendo as funções teste dadas por (5.26), segue que a matriz Be é dada por:

Be(x, y) =1

2Ae

[ye23 ye31 ye12

xe32 xe13 xe21

](5.28)

onde, novamente temos que xeij = xei − xej e yeij = yei − yej .Desta forma, a matriz Be é constante em cada elemento, dependendo apenas das coordenadas

dos nós em cada elemento. Isto implica que a aproximação obtida é linear e o gradiente de qual-quer combinação das funções teste será constante no interior do elemento triangular com três nós.Em uma dimensão tais características são as mesmas obtidas em elementos com dois nós, onde aaproximação possuia forma linear e sua derivada era constante.

Para simplicarmos nossos cálculos, ao invés de obtermos os coecientes para um elementotriangular no espaço, iremos mapear nosso elemento em um triangulo de referência utilizando ascoordenadas triangulares.

5.2.2 Coordenadas Triangulares

Para qualquer ponto P , temos que as coordenadas triangulares deste ponto são dadas por

ξj =AjA

(5.29)

onde Aj , j = 1, 2, 3 são as áreas dos triângulos gerados pelos dois nós diferentes de j e o ponto P,conforme podemos observar na Figura 5.4.

Figura 5.4: Coordenadas triangulares de um ponto P

Desta forma, quando o ponto P encontra-se sobre um dos nós, a coordenada triangular se tornaunitária e as outras coordenadas são nulas. Desta forma, temos que

ξi(xej , y

ej ) = δij (5.30)

onde

δij =

1, i = j0, i 6= j

Podemos observar, da denição (5.29), que

ξ1 + ξ2 + ξ3 =A1

A+A2

A+A3

A=A

A= 1

e, desta forma, pela propriedade (5.30) podemos escrever a relação entre as coordenadas (x, y) e as


coordenadas triangulares como

x =3∑i=1

xejξej , y =

3∑i=1

yejξej . (5.31)

Tal equação pode ser vista como um mapeamento entre um elemento de referência e o elementooriginal. Se olharmos o elemento no plano ξ1, ξ2, por 5.30, temos que ξ1(x1, y1) = 1, ξ2(x1, y1) =0, ξ2(x2, y2) = 1, ξ1(x2, y2) = 0 e ξ1(x3, y3) = ξ2(x3, y3) = 0, de forma que o triangulo de referênciaobtido é representado na Figura 5.5.

Figura 5.5: Triângulo de referência em coordenadas triangulares

Para concluirmos o desenvolvimento das coordenadas triangulares, devemos expressar as coor-denadas em termos de (x, y). A equação (5.31) fornece apenas duas equações, o que não é suciente.Porém, da denição (5.30) e da Figura 5.5, temos

ξ1 + ξ2 + ξ3 = 1

que, combinada a 5.29 fornece 1xy

=

1 1 1xe1 xe2 xe3ye1 ye2 ye3

ξ1

ξ2

ξ3

(5.32)

Resolvendo este sistema, obtemos ξ1

ξ2

ξ3

=1

2Ae

xe2ye3 − xe3ye2 ye23 xe32



︸︷︷︸

Me

1xy

(5.33)

onde Ae é a área do elemento em estudo, xeij = xei − xej e yeij = yei − yej .Concluímos assim que as coordenadas triangulares são lineares em x e y e, como satisfazem

(5.30), temos que as funções teste no plano ξ1, ξ2 devem ser idênticas ás funções teste lineares dotriângulo contido no plano x, y. Sendo ξ3 = 1 − ξ2 − ξ1, podemos então escrever a aproximaçãolinear como

ue =

3∑i=1

αei ξei (5.34)

Em resumo, temos que as funções teste dadas em (5.26) são idênticas aquelas obtidas nas coorde-nadas triângulares porém nas novas coordenadas a estrutura do triângulo de referência é conhecida(conforme visto na Figura 5.5).

5.3 FORMULAÇÃO GERAL DO MÉTODO PARA PROBLEMAS BIDIMENSIONAIS 43

5.3 Formulação geral do método para problemas bidimensionais

Após obtidas a formulação fraca, as funções teste e as funções de ponderação para elementostriangulares com três nós, partimos agora para a aplicação do MEF por meio da discretização dodomínio em estudo por elementos triangulares e da obtenção de um sistema de equações geradopela formulação fraca.

Para isto particionamos o domínio do problema em nel elementos triangulares de três nós e,conforme vimos anteriormente, as aproximações geradas sobre cada elemento serão da forma (5.21)onde cada um dos coecientes αei , i = 1, 2, 3 satisfazem a condição de continuidade (5.22) em relaçãoa aproximação obtida no elemento adjacente. Além disso, sendo Ωe um determinado elementotriangular e e usando que ue(x, y) = 0, se (x, y) /∈ Ωe, temos que a aproximação global obtida serádada por

u(x, y) =

nel∑e=1

ue(x, y) (5.35)

De maneira análoga ao que zemos para problemas unidimensionais, devido as descontinuidadesexistentes entre as funções teste em cada um dos elementos, avaliamos a integração da formulaçãofraca como a soma das integrais sobre cada um dos elementos. Desta forma, temos que

nel∑e=1

[∫ ∫Ω

(∇we)TD∇uedΩ +

∫ ∫ΩwefdΩ−

∫Γe

wgdΓ

]= 0 (5.36)

onde Ωe representa o conjunto de pontos que forma o elemento triangular e.Utilizando funções teste e funções de ponderação da forma

ue(x, y) = Ne(x, y)(αe)T ⇒ ∇ue(x, y) = Be(x, y)(αe)T

we(x, y) = Ne(x, y)(we)T ⇒ ∇we(x, y) = Be(x, y)(we)T

com Ne(x, y) e Be(x, y) dados por (5.26) e (5.28) respectivamente e substituindo em (5.36) temos:

nel∑e=1

(we)T

∫ ∫

Ωe

(Be)TDBedΩ︸︷︷︸Ke

αe +

∫ ∫Ωe

(Ne)T fdΩ−∫

Γe

(Ne)T gdΓ︸︷︷︸fe

= 0 (5.37)

Assim como zemos nos problemas unidimensionais, denimos

• A Matriz de Rigidez Ke do elemento e:

Ke =

∫ ∫Ωe

(Be)TDBedΩ (5.38)

• O Vetor de Forças Externas f e do elemento e:

f e =

∫ ∫Ωe

(Ne)T fdΩ︸︷︷︸fΩe

−∫

Γe

(Ne)T gdΓ︸︷︷︸fΓe

(5.39)

Observe que como a função peso é arbitrária e o número de graus de liberdade do problemaé igual ao número de constantes a serem obtidas podemos considerar que os valores de we nãoinuenciam na obtenção do sistema linear formado pelas matriz global de rigidez e o vetor globalde forças.


5.3.1 Geração da Matriz Global de Rigidez e do Vetor Global de Forças

Em problemas unidimensionais tinhamos que a geração da matriz global de rigidez e do vetorglobal de forças externas era feito pela soma das componentes relativas aos pontos que pertenciam adois elementos consecutivos. No caso bidimensional o procedimento realizado é equivalente porém,deve se levar em consideração que dois elementos adjacentes possuem uma aresta em comum e, destaforma, as componentes relativas aos dois pontos que geram tal aresta devem ser somadas na entradacorrespondente da matriz global. A Figura 5.6 fornece uma visualização interessante do processo aser realizado, onde um domínio arbitrário é discretizado por meio de elementos triangulares de trêsnós e é apresentada a montagem da matriz global de rigidez.

Figura 5.6: Montagem da matriz global de rigidez e do vetor global de forças. [Fla]

Assim como para a matriz de rigidez global obtida para problemas unidimensionais, temos quea entrada kij da matriz global de rigidez e nula se, e somente se, os pontos associados a ela nãoformam um elemento, isto é, se os pontos (xi, yi) e (xj , yj) não são aresta de nenhum elemento damalha.

Como conhecemos o valor da função nos ponto pertencentes a Γe, temos que os coecientesassociados a estes pontos são conhecidos. Podemos então reescrever o sistema linear, de forma aeliminar as linhas e as colunas associadas a essas variáveis, exatamente como fazíamos nos problemasunidimensionais. Temos então o sistema linear

Kα = f (5.40)

onde K é a matriz de rigidez, f o vetor global de forças e α o vetor contendo os coecientes a seremobtidos, já removidos os pontos pertencentes ao contorno essencial do problema.

Capítulo 6

Aspectos computacionais do problema

bidimensional

6.1 Geração da malha do problema

A geração de malhas bidimensionais por si só trata-se de um vasto campo de estudo, tanto paramatemáticos aplicados como para cientistas da computação. Como estamos interessados muito maisna aplicabilidade do gerador de malhas do que em sua eciência, não abordaremos aqui como é feitaa geração da malha ou qual estratégia deve ser tomada na enumeração dos nós da malha.

Em nossas simulações computacionais do MEF iremos utilizar um gerador de malhas já im-plementado, denominado EasyMesh. Trata-se de um gerador de malhas bastante simples e ecaz,que pode ser obtido livremente em http://www-dinma.univ.trieste.it/nirftc/research/easymesh/.Neste site também encontra-se disponível o visualizador de malhas ShowMesh que possibilita avisualização da malha gerada. As Figuras 7.13, 7.14 e 7.15 foram geradas pelo ShowMesh.

A seguir, iremos mostrar como é construído o arquivo de entrada para o gerador de malhas equais as informações e de que maneira são retornadas pelo gerador.

6.1.1 Dados de Entrada

O EasyMesh recebe os dados do domínio na forma de um arquivo simples de texto com aextensão .d. Arquivos desta forma podem ser facilmente gerados em compiladores C. O arquivodeve possuir a seguinte estrutura:

• Primeira linha: <Número de nós>

• Linhas seguintes: <Número do nó:> <x> <y> <Renamento> <Tipo de contorno>

• Próxima linha: <Número de segmentos>

• Linhas seguintes: <Número do segmento:> <Início> <Fim> <Tipo de contorno>

onde a constante renamento é um valor real maior que zero que representa o renamento damalha ao redor deste nó, ou seja, quanto menor o valor, mais na é a malha ao redor deste ponto.A constante tipo de contorno permite atribuir diferentes valores inteiros para diferentes trechos docontorno. Os valores de início e m são os pontos que formam o segmento.

É possível também realizar comentários no arquivo de entrada. Para isso, deve-se iniciar a linhacom o caracter '#' e, desta forma, tudo o que for escrito após não será interpretado pelo EasyMesh.A Figura 6.1 mostra o arquivo de entrada utilizado para gerar a malha da simulação da chaveinglesa.

45

http://www-dinma.univ.trieste.it/nirftc/research/easymesh/

46 ASPECTOS COMPUTACIONAIS DO PROBLEMA BIDIMENSIONAL 6.1

Figura 6.1: Arquivo de entrada para o EasyMesh do problema 7.8.

6.1.2 Dados de Saída

Sendo entrada.d o arquivo de entrada fornecido ao gerador, temos que são fornecidos como saidado programa três arquivos:

• entrada.n: Apresenta os dados sobre os nós da malha gerada. Fornece na primeira linha onúmero de nós gerado pela malha e as linhas seguintes são escritas da seguinte forma:

<Número do nó:> < x > < y > < marcador >

onde x e y são as respectivas coordenadas de cada nó e marcador indica a qual contorno onó pertence. Caso o nó esteja no interior do domínio, este valor é igual a zero.

• entrada.e: Apresenta os dados sobre os elementos gerados. Fornece na primeira linha o númerode elementos no qual o domínio foi particionado e as linhas seguintes são escritas da seguinteforma:

<Número do elemento:> < i > < j > < k > < ei > < ej > < ek > < si > < sj > < sk >< xV > < yV > < Tipo de contorno>

onde i, j, k são os nós que formam o elemento, ei, ej e ek são os vizinhos do elemento emquestão (caso o elemento esteja no contorno, o valor sera denido como −1), si, sj e sk sãoos segmentos que formam o domínio, xV e yV são as coordenadas do elemento circuncentrico

6.2 QUADRATURAS NUMÉRICAS PARA PROBLEMAS BIDIMENSIONAIS 47

(que não serão utilizadas em nossos estudos). A Figura 6.2 apresenta a visualização dos dadosfornecidos por este arquivo.

Figura 6.2: Dados do arquivo entrada.e

• entrada.s: Apresenta os dados sobre os segmentos de reta que formam os triângulos gerados.Fornece na primeira linha o número total de segmentos de retas da malha gerada e as linhasseguintes são escritas da forma:

<Número do segmento:> < c > < d > < ea > < eb > <marcador>

onde c e d são os pontos onde começa e termina o respectivo segmento, ea e eb são os elementosa esquerda e a direita do segmento, respectivamente. Se o segmento zer parte do contorno,temos que o elemento a direita do segmento não existe. Nesses casos, o valor denido é iguala −1. A Figura 6.3 apresenta a visualização dos dados fornecidos por este arquivo.

Figura 6.3: Dados do arquivo entrada.s

6.2 Quadraturas numéricas para problemas bidimensionais

Assim como no MEF unidimensional, a montagem das matrizes de rigidez de cada elementoé feita por integração numérica. Desta forma, apresentaremos aqui uma maneira de estender osconceitos vistos no Capítulo 3 relativos a integração númerica.


6.2.1 Integração sobre elementos triangulares

De maneira análoga á utilizada para problemas unidimensionais, temos que a fórmula de inte-gração para problemas bidimensionais será dada por∫ ∫

Ωe

f(x, y)dΩ =

ngp∑i=1

Wi|Je(ξi)|f((ξ1, ξ2)i) (6.1)

onde Ωe ∈ R2 é uma superfície triangular, ngp é o número de pontos utilizado na quadratura,Je éo jacobiano da transformação para o triângulo de referência, Wi as ponderações e ξi os pontos aserem utilizados na quadratura. Os pesos e os pontos da quadratura encontram-se na Tabela 6.1.

Ordem de Integração Grau de Precisão ξ1 ξ2 Pesos0.1666666666 0.1666666666 0.1666666666

Três pontos 2 0.6666666666 0.1666666666 0.16666666660.1666666666 0.6666666666 0.16666666660.1012865073 0.1012865073 0.06296959030.7974269853 0.1012865073 0.06296959030.1012865073 0.7974269853 0.0629695903

Sete pontos 5 0.4701420641 0.0597158717 0.06619707640.4701420641 0.4701420641 0.06619707640.0597158717 0.4701420641 0.06619707640.3333333333 0.3333333333 0.1125

Tabela 6.1: Nós da quadratura de Gauss-Legendre e suas ponderações para domínios triangulares

Observe que o jacobiano da transformação bidimensional é o equivalente ao fator (b−a)2 que

surge em (3.7). Como estamos utilizando triângulos com três nós e de lados retos, temos que estejacobiano é constante e dado por:

Je =

[xe1 − xe3 ye1 − ye3xe2 − xe3 ye2 − ye3

]ou seja, o jacobioano resultante é o dobro da área de Ωe. Temos então que a equação (6.1) é dadapor: ∫

Ωe

f(x, y)dΩ = 2A(Ωe)

ngp∑i=1

Wif((ξ1, ξ2)i) (6.2)

onde A(Ωe) é a área do elemento Ωe.

6.3 Armazenamento de matrizes esparsas

Da implementação do MEF utilizando elementos triangulares, assim como na implementaçãounidimensional, temos que a matriz do sistema normal (5.40) também apresenta estrutura esparsa.No entanto, diferentemente das matrizes obtidas nos problemas unidimensionais, não temos a es-trutura de banda xa pré determinada.

No entanto, podemos utilizar uma estrutura de armazenamento especial, de forma que não sejanecessário o armazenamento das entradas nulas do sistema. Uma estrutura de dados bem simplesutiliza apenas três vetores: um vetor de valores reais k contendo os elementos não nulos da matriz,um vetor de inteiros c que contém os índices das colunas dos elementos da matriz armaznados em

6.4 MÉTODO DOS GRADIENTES CONJUGADOS (MGC) 49

k e um vetor de inteiros p contendo os ponteiros para o início de cada linha nos vetores k e c.Especicamente, para uma matriz An×n com N elementos não nulos temos:

• Um vetor k, de tamanho N , contendo os elementos não nulos de A armazenados sequencial-mente por linhas.

• Um vetor c, de tamanho N , contendo os índices das colunas dos elementos de A armazenadosem k.

• Um vetor p, de tamanho n+ 1, contendo os ponteiuros para o início de cada linha nos vetoresk e c. Logo, temos que p[i] é a posição nos vetores k e c onde a linha i começa, com 1 ≤ i ≤ n.

Este formato de armazenamento recebe o nome de Compressed Sparse Row em inglês e é nor-malmente chamado pela sigla CSR. Existem outras estratégias de armazenamento para matrizesesparsas. Para o leitor interessado, [Saa03] trata-se de uma obra bastante interessante e acessível.

6.4 Método dos Gradientes Conjugados (MGC)

Para a resolução do sistema linear associado iremos utilizar um método iterativo de solução desistemas lineares, o chamado Método dos Gradientes Conjugados (MGC).

6.4.1 Idéia do Método

Sendo A uma matriz simétrica e positiva denida, o método dos Gradientes Conjugados paraa solução de um sistema linear Ax = b parte de uma aproximação inicial qualquer x0, a partir daqual denimos uma direção inicial e um resíduo inicial

d0 = r0 = b−Ax0

Se denotarmos a solução exata do sistema linear por x, temos que a relação entre o resíduo e oerro e0 = x− x0 é escrita por:

r0 = Ae0

Vamos então eliminar do erro qualquer componente na direção inicial d0. Para isso, ao invés deutilizarmos o produto interno usual do Rn, utilizaremos um produto interno induzido pela matrizA, uma vez que

< e0,d0 >A= < Ae0,d0 >= < r0,d0 >

que é possível de se obter, mesmo desconhecendo e0. Assim, o novo erro e1, ortogonal à direçãoinicial d0, é dado por

e1 = e0 −< r0,d0 >

< d0,d0 >d0

Note que a equação anterior equivale a obtermos uma nova aproximação

x1 = x0 +< r0,d0 >

< d0,d0 >Ad0

que podemos calcular apenas utilizando a matriz A e a aproximação inicial x0. O processo continuacom a adição de novas direções denidas através dos resíduos (ortogonalizados em relação as direçõesanteriores), sempre tornando o erro ortogonal aos espaços formados entre o resíduo e as direções d.

O processo terminará quando tivermos o erro ortogonal a um subespaço de dimensão n do próprioRn, ou seja, o próprio Rn. Porém, isto só é possível com o erro sendo nulo, ou seja, se chegarmos asolução exata do sistema. No entanto, devido aos erros de arredondamento, não conseguimos obter


a solução exata devendo assim denir um critério de parada de forma que o algoritmo execute asiterações do método até que a aproximação obtida seja satisfatória.

Observe que, em nossos problemas, a matriz A encontra-se armazenada no formato CSR. Destaforma, o produto matriz-vetor e é realizado conforme o Algoritmo 3:

Algoritmo 3: Produto entre uma matriz armazenada no formato CSR e um vetorEntrada: Vetores k1×N ∈ R, c1×N ∈ N,p1×n+1 ∈ N, f1×n ∈ RSaída: Vetor x1×n ∈ Rn contendo o vetor obtido pelo produto da matriz K pelo vetor fpara i = 0 até n− 1 faça

k1← p[i];k2← p[i+ 1]− 1;para j = k1 até k2 faça

x[i]← x[i] + k[j] ∗ f [c[j]];m para

m para

6.4.2 Critério de Parada

Como apresentado na seção anterior, devido aos erros de arredondamento temos que o MGCnão obtem a solução exata do sistema linear em estudo tornando-se necessário a denição de umcritério de parada para o algoritmo. Alguns critérios que podem ser utilizados são:

• Diferença relativa entre duas iterações

• Máxima diferença entre duas iterações

• Norma do resíduo

Um aspecto importante a ser observado é o fato de que, mesmo se o sistema linear fosse resolvidode maneira exata, a solução numérica obtida pelo MEF ainda sim é diferente do valor exato da solu-ção da equação diferencial. Sendo assim, temos que não é necessário que a solução do sistema linearseja exata, bastando que a aproximação obtida seja da mesma ordem de precisão da aproximaçãonumérica a ser gerada pelo MEF.

Em nossas simulações iremos tomar como critério a norma euclideana do resíduo menor que10−12 pois, desta forma as soluções numéricas obtidas na resolução do sistema linear serão sucientespara não afetar a ordem de precisão do MEF. Mais detalhes podem ser obtidos em [Saa03].

6.4.3 Propriedades do Método

Nesta seção iremos apresentar algumas propriedades do MGC. As demonstrações das proprie-dades a seguir não serão feitas aqui (são feitas por indução e razoávelmente trabalhosas). Ao leitorinteressado, as demonstrações podem ser obtidas em [SB02].

Propriedade 1 (Independencia Linear). Seja A ∈ Rn×n uma matriz simétrica positiva denida econsidere dj , j = 0, 1, 2, · · · , k um conjunto de vetores do Rn tal que

< Adi,dj >= 0 (6.3)

para todo i, j = 0, 1, 2, · · · , k. Temos então que os vetores dj , j = 0, 1, 2, · · · , k são linearmenteindependentes.

Propriedade 2 (Ortogonalidade entre os resíduos e as direções). Os vetores dk e rk, k = 0, 1, 2, · · · , ngerados pelo Algoritmo 4 são ortogonais em relação ao produto interno usual do Rn, isto é

(ri)Tdj = 0

6.4 MÉTODO DOS GRADIENTES CONJUGADOS (MGC) 51

para 0 ≤ j < i ≤ n.

Propriedade 3 (Igualdade entre produtos internos). Os vetores ri e di gerados pelo algoritmo 4satisfazem:

rTi di = rTi ri

para i ≤ n.

Propriedade 4 (Ortogonalidade entre os resíduos). Os vetores dk, k = 0, 1, 2, · · · , n gerados peloAlgoritmo 4 são ortogonais em relação ao produto interno usual do Rn, isto é

rTi rj = 0

para 0 ≤ j < i ≤ n.

Propriedade 5 (Solução em n passos). Considere o sistema linear

Ax = b (6.4)

com A ∈ Rn×n positiva denida, b ∈ Rn e x0 ∈ Rn uma aproximação inicial da solução do sistema.Desta forma, a solução do sistema utilizando o Algoritmo 4 é obtida em, no máximo, n passos.

6.4.4 Algoritmo do Método

Das idéias apresentadas na seção anterior, temos que o método dos gradientes conjugadosencontra-se descrito pelo Algoritmo 4:

Algoritmo 4: Método dos Gradientes Conjugados com Armazenamento CSR

Entrada: Matriz A ∈ Rn×n, vetor b ∈ Rn, a aproximação inicial x0 ∈ Rn e as constanteseps ∈ R e itmax ∈ N

Saída: Vetor x ∈ Rn de aproximações da solução do sistemaw← produto Ax0;para i = 0 até n− 1 faça

d[i] = r[i] = b[i]− w[i];m para

erro← 2 ∗ eps;k ← 0;enquanto k ≤ itmax e erro > eps faça

w← produto Adk;P1← produto interno < w,d >;P2← produto interno < r, r >;alfa← P2

P1 ;para i = 0 até n− 1 faça

x[i]← x[i] + alfa ∗ d[i];r[i]← r[i]− alfa ∗ w[i];

m para

P3← produto interno < r, r >;beta← P3

P2 ;para i = 0 até n− 1 faça

d[i]← r[i]− beta ∗ d[i];m para

eps← norma euclideana r;m enqto


6.4.5 Pré-condicionamento do sistema

A solução do sistema linear consiste em uma das etapas com maior custo computacional naimplementação do MEF. Uma estratégia bastante utilizada para melhorar o desempenho do MGCé o pré-condicionamento do sistema, onde algumas operações são realizadas visando acelerar oprocesso de convergência do algoritmo.

Existem diversas técnicas para o pré-condicionamento do sistema porém não existe um padrãoa ser seguido sobre qual delas é a mais ecaz para todos os problemas.

Seguindo os objetivos deste trabalho, não iremos nos aprofundar neste tópico. Ao leitor interes-sado, a obra [Saa03] faz um estudo interessante sobre os métodos de pré-condicionamento.

6.5 Convergência numérica do algoritmo

Durante o estudo teórico do MEF utilizando elementos triangulares, vimos que o erro obtidoentre a solução exata e a aproximação obtida utilizando-se de elementos triangulares é da forma

erro(k) ∼= Ch2 (6.5)

onde k é a malha de estudo e h é uma medida que indica o tamanho de cada elemento.Como a geração da malha é feita pelo EasyMesh de maneira a manter uma certa uniformidade

aos elementos do domínio, temos que as malhas renadas não são diretamente frações das anteriores.No entanto, podemos obter uma aproximação da ordem do método notando que:

erro(k1)

erro(k2)∼=Chn1Chn2

=hn1hn2

=

(h1

h2

)n⇒ n ∼=

log(erro(k1)erro(k2)

)log(h1h2

) (6.6)

Desta forma, tomando como medida do tamanho do elemento

h =√

2 · maxi∈[1,nel]

(Ai) (6.7)

onde nel é o número de elementos da malha e Ai é a área do i-ésimo elemento, obtemos uma maneirade avaliar o comportamento do erro em relação ao renamento da malha.

Capítulo 7

Simulações bidimensionais

Nesta seção faremos a modelagem e a aplicação do MEF para problemas envolvendo conduçãode calor e deformações em sólidos.

7.1 Problema de Poisson no Quadrado

Para iniciarmos o estudo da aplicação do MEF em problemas bidimensionais, iremos analisarnumericamente a solução de um problema modelo. Para isso, escolhemos o seguinte problema:

• Obter uma função u : [0, 1]× [0, 1]→ R, de classe C2, tal que:

∇ · (∇u) = 0 ∀(x, y) ∈ Ω

∇u · ~n =

[exsenyexcosy

]· ~n ∀(x, y) ∈ Γt

u(x, y) = exsen(y) ∀(x, y) ∈ Γe

ondeΩ = (x, y) ∈ R2 : (x, y) ∈ (0, 1)× (0, 1)

Γt = (x, y) ∈ R2 : (x = 1, y = t) ∪ (x = t, y = 1), t ∈ (0, 1]

Γe = (x, y) ∈ R2 : (x = t, y = 0) ∪ (x = 0, y = t), t ∈ [0, 1]

e ~n é a normal unitária exterior ao quadrado. Observe que ∂Ω = Γt ∪ Γe é a fronteira doquadrado e que Γt ∩ Γe = ∅.

Observe que o valor da função no contorno essencial satisfaz a equação diferencial e a condiçãode contorno natural do problema, pois

∇ · (∇u) = ∆u =∂2u

∂x2+∂2u

∂y2= exsen(y)− exsen(y) = 0

e

∇u =

[exsenyexcosy

]Temos então que a solução do problema é dada por:

u(x, y) = exsen(y)

53

54 SIMULAÇÕES BIDIMENSIONAIS 7.1

7.1.1 Simulação Computacional

Realizamos então a simulação do problema para malhas com 2, 14, 44, 150, 602 e 2390 (Fi-gura 7.2). Pela análise da Figura 7.1 podemos observar que as superfícies geradas nas simulaçõesconvergem para a superfície dada pela solução do problema conforme o aumento do número deelementos.

Figura 7.1: Visualização das aproximações numéricas obtidas.

7.1.2 Análise de erro

Para analisarmos o comportamento do erro calculamos, a partir da solução númerica obtida,o valor da aproximação em 676 pontos igualmente espaçados do quadrado unitário (grid 26x26) etomamos como medida do erro obtido na aproximação:

||u− v||∞ = max0≤i≤25,0≤j≤25

|u(xi, yj)− v(xi, yj)|

onde u é a solução do problema e v a aproximação obtida pelo algoritmo nas malhas contidas naFigura 7.2, ambas calculadas nos pontos (xi, yj) do grid. Usando como medida do tamanho doelemento (6.7), obtivemos os dados da Tabela 7.1.

7.1 PROBLEMA DE POISSON NO QUADRADO 55

Figura 7.2: Renamento da malha de domínio do problema.

Qtd Pontos Qtd Elementos h Erro máximoMalha 1 4 2 1.000 0.529211Malha 2 12 14 0.424 0.113973Malha 3 31 44 0.261 0.040709Malha 4 92 150 0.140 0.011695Malha 5 334 602 0.070 0.003124Malha 6 1260 2390 0.037 0.000903

Tabela 7.1: Dados obtidos com o renamento da malha.

7.1.3 Ordem de convergência do Método

Tomando as relações entre os erros obtidos e entre os tamanhos dos elementos para duas malhassucessivas podemos estimar a ordem do método, conforme observamos na Tabela 7.2 e no gráco


da Figura 7.3:

n(i)/n(j) e(i)/e(j) h(i)/h(j) Ordem4/12 4.6433 2.3570 1.790712/31 2.7997 1.6222 2.218031/92 3.4808 1.8681 1.995992/334 3.7435 2.0000 1.9044334/1260 3.4595 1.8708 1.9814

Tabela 7.2: Dados obtidos com o renamento da malha do problema, onde n(i), e(i) e h(i) representam aquantidade de nós, o erro máximo obtido e o tamanho h obtido na malha i, respectivamente.

Figura 7.3: Comportamento da ordem estimada do método em relação ao número de iterações de rena-mento de malha do problema.

Pela análise dos dados, podemos concluir que o erro em cada iteração do método cai aproxima-damente ao quadrado da medida h utilizada.

7.2 Distribuição de calor em uma chaminé

Considere agora o problema da distribuição de calor em uma chaminé cujo interior é mantidoa uma temperatura constante igual a 140C e a face exterior da chaminé encontra-se com umatemperatura constante igual a 10C.

Sabemos também que a chaminé é constituída de concreto em sua parte interior (cuja condu-tividade térmica é igual a 2.0W C−1 ) e de tijolos em parte externa (cuja condutividade térmicaé igual a 0.9W C−1). As medidas do problema encontram-se na Figura 7.5. Devido a simetria doproblema, iremos analisar a distribuição de calor em 1/4 da chaminé apenas.

7.2.1 Modelagem do Problema

Para obtermos as equações que regem o problema, considere uma chapa de espessura unitária,sujeita a uma fonte de calor s(x, y) (energia por unidade de área) e um volume de controle, de lados∆x e ∆y, conforme mostra a Figura 7.4:

7.2 DISTRIBUIÇÃO DE CALOR EM UMA CHAMINÉ 57

Figura 7.4: Modelagem matemática do problema [FB07].

Lembramos que em duas dimensões a lei de Fourier para materiais isotrópicos é dada por:

~q = −σ∇T

onde σ > 0 é a condutividade térmica do material e T a temperatura em cada ponto. Observe que,da simetria do problema, temos que a componente ~q · ~n é nula nos segmentos BC e DA da Figura7.5.

Estando o corpo em regime permanente, temos que o uxo de calor que sai pelas fronteiras dovolume de controle (~q = σ∇T ) é igual ao calor gerado s(x, y). Isto implica que o uxo de calor quesai tem de ser igual a energia caloríca gerada pela fonte.

O uxo de calor sobre o volume de controle pode ser expresso em termos de duas componen-tes: a componente tangencial (qt) ao contorno e a componente normal ao contorno (qn). Como acomponente tangencial ao contorno não contribui para a entrada ou saída de calor no volume decontrole, temos que o uxo depende apenas da componente normal, dada por:

qn = ~q · ~n = qxnx + qyny

Observe que na Figura 7.5, sobre AD, em que ~n = −~i, o uxo de calor entrando é −qn = −~q ·(−~i) =qx, enquanto sobre BC, em que ~n =~i, o uxo de calor entrando é −qn = −~q ·~i = −qx.

Desta forma, o balanço de energia no volume de controle é dado por:

qx(x− ∆x

2, y)∆y − qx(x+

∆x

2, y)∆y + qy(x, y −

∆y

2)∆x− qy(x, y +

∆y

2)∆x+ s(x, y)∆x∆y = 0;

Dividimos então a equação anterior por ∆x∆y e usando a denição de derivada parcial, obtemos:

∂qx∂x

= lim∆x→0

qx(x+ ∆x2 , y)− qx(x− ∆x

2 , y)

∆x

∂qy∂y

= lim∆y→0

qy(x, y + ∆y2 )− qy(x, y − ∆y

2 )

∆y

⇒ ∂qx∂x

+∂qy∂y− s(x, y) = 0

∇~q − s(x, y) = 0 (7.1)

Como não existe nenhuma fonte de calor externa temos que s(x, y) = 0, obtendo assim a equação:

∇ · (σ∇T ) = 0


Desta forma, a equação diferencial do problema é dada por:

∇ · (σ∇T ) = 0; ∀(x, y) ∈ Ωσ∇T · ~n = 0; ∀(x, y) ∈ ΓtT (x, y) = 140; ∀(x, y) ∈ Γe1T (x, y) = 10; ∀(x, y) ∈ Γe2

Figura 7.5: Modelagem do problema.

ondeΩ = (x, y) ∈ R2 : (x, y) ∈ int(ABCD),

Γt = (x, y) ∈ R2 : (x, y) ∈ (BC ∪AD),

Γe1 = (x, y) ∈ R2 : (x, y) ∈ (AB),

Γe2 = (x, y) ∈ R2 : (x, y) ∈ (CD),

e com σ = σ(x, y) assumindo os valores 0.9 ou 2.0 dependendo do material utilizado na região emestudo da chaminé.


Como não conhecemos a solução exata do problema iremos analisar o comportamento das solu-ções obtidas pelo renamento da malha de estudo. Geramos então malhas com 32 nós (40 elementos),73 nós (112 elementos), 284 nós (502 elementos) e 987 nós (1844 elementos) conforme mostra a Fi-gura 7.6, e calculamos aproximações da solução em pontos igualmente espaçados do quadrado delado 4.

Note que cada elemento gerado pertence a somente um dos materiais que compõem a chaminé.Na modelagem do nosso problema, temos que a condutividade térmica é constante nos trechos quecompõem a chaminé, mas apresenta uma descontinuidade no segmento que une os componentes.

Obtendo as aproximações em alguns pontos do dominio temos que os dados obtidos sugerem aconvergência do método pois os resultados gerados pelo algoritmo aplicado as malhas cada vez maisrenadas aparentemente estabilizam em um dado valor, conforme podemos observar na Tabela 7.3:

Utilizando os resultados do algoritmo gerado com o uso de 987 nós (1844 elementos), obtivemoso gráco da Figura 7.7 para a distribuição da temperatura.

7.3 DISTRIBUIÇÃO DE CALOR EM UMA CHAMINÉ 59

Figura 7.6: Renamento da malha de domínio do problema.

Ponto 32 nós 73 nós 284 nós 987 nós(0.5,3.5) 53.7737 53.6019 53.5061 53.4921(1.0,3.5) 53.131 52.5002 52.4911 52.4148(1.5,3.5) 50.9958 50.3440 50.1919 50.1291(2.0,3.5) 47.2271 46.4890 46.3779 46.1439(2.5,3.5) 43.7773 40.6480 40.1772 40.1287(0.5,2.5) 118.9358 118.5459 118.4578 118.4267(1.0,2.5) 118.1313 117.2222 117.17669 117.1001(1.5,2.5) 116.0670 114.15376 113.8418 113.7075(2.0,2.5) 112.8553 106.6802 105.3334 104.9242(2.5,2.5) 95.4843 90.5792 89.5212 89.0996

Tabela 7.3: Temperaturas (em C) obtidas na análise do problema com renamento da malha.

7.2.3 Análise dos resultados obtidos

Apesar de σ possuir uma descontinuidade no segmento que divide a composição da chaminéentre concreto e tijolos ainda sim é possível observar que a cada iteração de renamento da malhautilizada as soluções obtidas convergem para um determinado valor xo.

É possível observar também que, devido a descontinuidade de σ no segmento que divide osmateriais que compõem a chaminé, a solução obtida também será descontinua nesse trecho (noteque ocorre uma brusca queda de temperatura na faixa sobre o segmento (t, 3.0), (3.0, t),∀t ∈ [0, 3])porém é possível observar que o uxo de calor, dado por −σ∇T é continuo em todo domínio.


Figura 7.7: Distribuição da temperatura com o MEF utilizando 1844 elementos triangulares.

7.3 Análise de deformações

Considere o seguinte problema: uma chave inglesa encontra-se presa a um parafuso que não semovimenta. Sendo aplicada uma força constante no cabo da chave, temos que ocorre uma deforma-ção no corpo da chave. Vamos utilizar o Método dos Elementos Finitos para obtermos aproximaçõesdas deformações sofridas pelo sólido.

Figura 7.8: Descrição do problema que será abordado: a chave inglesa encontra-se presa a um parafuso quenão se movimenta enquanto sofre a ação de uma força constante.

Para obtermos as equações que regem o problema, iremos analisar a cinemática do problema, arelação constitutiva do material que relaciona as tensões e as deformações e as equações de equilíbriodo corpo.

7.3.1 Cinemática

Dado um ponto qualquer do corpo em estudo, o deslocamento ocorrido neste ponto pode serexpresso na forma vetorial por:

~u = ux~i+ uy~j (7.2)

7.3 ANÁLISE DE DEFORMAÇÕES 61

As variações nos comprimentos nos segmentos de linha innitesimais nas direções x e y, divididospelos comprimentos de linha nos fornecem as deformações extensoriais nos eixos x e y, respectiva-mente:

εx = lim∆x→0ux(x+∆x,y)−ux(x,y)

∆x = ∂ux∂x

εy = lim∆y→0uy(x,y+∆y)−uy(x,y)

∆y =∂uy∂y

(7.3)

Além das deformações extensoriais, temos também a deformação de cisalhamento, que indica avariação angular em relação aos vetores unitários nas direções x e y, dada por:

γxy = lim∆y→0ux(x,y+∆y)−ux(x,y)

∆y +

lim∆x→0uy(x+∆x,y)−uy(x,y)

∆x

⇒ γxy =∂uy∂x + ∂ux

∂y = α1 + α2

(7.4)

Figura 7.9: Visualização dos componentes do deslocamento [FB07].

Podemos então escrever um vetor de deformações contendo as deformações extensoriais nasdireções x e y e a deformação de cisalhamento na forma ε = [εx, εy, γxy] ou em forma matricialε = ∇S~u, onde ∇S é dado por:

∇S =

∂∂x 0

0 ∂∂y

∂∂y

∂∂x

Desta forma, obtemos uma primeira relação entre o deslocamento e as deformações sofridas pelosólido em estudo.

7.3.2 Tensões e trações

As tensões em duas dimensões que correspondem às forças por unidade de área agindo sobre osplanos normais aos eixos x ou y são chamadas trações. As trações sobre o plano com o vetor normal~n alinhado ao longo dos eixos x e y são denotadas por:

~σx = σxx~i+ σxy~j

~σy = σyx~i+ σyy~j(7.5)

A partir do equilíbrio do momento em um quadrado unitário, pode ser mostrada a igualdadeσxy = σyx. Uma possível forma de escrevermos as tensões é:

σT =[σxx σyx σxy

]Desta forma, é conveniente escrever as componentes da tensão em uma matriz simétrica τ2×2 dadapor:

τ =

[σxx σxyσyx σyy

]


Figura 7.10: Visualização dos componentes da tensão [FB07].

Sendo assim, especicada uma tração (~t0)T = [t0x, t0y] podemos reescrevê-la em função da tensãocomo:

~t0 = τ~n (7.6)

7.3.3 Equilíbrio

Considere um corpo de espessura unitária, sob a ação de uma força de tração ~t ao longo de seucontorno e uma força de campo ~b por unidade de volume, conforme mostra a Figura 7.7.

Figura 7.11: Visualização dos componentes das forças em equilíbrio no quadrado unitário [FB07].

Considerando o equilíbrio do dominio innitesimal de espessura unitária, para um problemaestático (isto é, sem efeitos dinâmicos) temos que a equação de equilíbrio é dada por:

∂ ~σx∂x

+∂ ~σy∂y

+~b = ~0 (7.7)

Multiplicando a equação anterior pela normal unitária e, reescrevendo na forma matricial, temos:

∇TSσ +~b = ~0 (7.8)

7.3.4 Aplicação do Método dos Elementos Finitos

Inicialmente, utilizando as equações da modelagem do problema iremos determinar as equaçõesdiferenciais que regem o problema e seus respectivos contornos. Tal formulação é denominada aformulação forte do problema.

Pelas relações vistas na seção anterior, temos o problema em estudo é dado por:

• Obter u(x, y) ∈ C2 que satisfaça∇TSσ +~b(x, y) = 0 (7.9)


denida no domínio Ω ∈ R2 sujeita a

σ = D∇Su (7.10)

~t0 = σ~n, ∀(x, y) ∈ Γt (7.11)

~u = ~u0, ∀(x, y) ∈ Γu (7.12)

onde a equação (7.9) é obtida da cinemática do problema relacionando as deformações extensoriaisnos eixos x e y. As equações (7.10) e (7.11) são obtidas pela relação constitutiva do elementocombinada com a condição de equiíbrio do corpo, onde ~n é o vetor normal unitário. O contorno dodomínio do problema Ω é dado por Γ = Γt ∪ Γu sendo Γt o trecho do contorno onde é conhecido atração e Γu o trecho do contorno onde é especicada a deformação.

Para a aplicação do Método dos Elementos Finitos inicialmente denimos as funções teste u eas funções de ponderação w, de modo que u ∈ U e w ∈ U0, onde

U = u(x, y)| u = u0, ∀(x, y) ∈ ΓeU0 = w(x, y)| w = 0, ∀(x, y) ∈ Γe

Multiplicamos então as equações (7.9) e (7.11) pela função de ponderação w e integramos sobre osrespectivos domínios obtendo: ∫ ∫

Ω(w∇TSσ + w~b)dΩ = 0 (7.13)

∫Γn

wτ~ndΓt =

∫Γt

w~t0dΓn (7.14)

Aplicando o teorema de Green em (7.13) , temos:∫ ∫Ωw∇TSσ =

∮Γwσ~ndΓ−

∫ ∫Ω

~∇w~σdΩ (7.15)

e, usando que a função de ponderação deve ser identicamente nula em Γe, substituimos (7.13) naigualdade anterior obtendo∫ ∫

Ω(∇Sw)TσdΩ =

∫Γn

wT~t0dΓ +

∫ ∫ΩwT~bdΩ (7.16)

Desta forma, podemos reescrever o problema como: obter u(x, y) ∈ U tal que∫Ω

(∇Sw)TD∇SudΩ︸︷︷︸K

=

∫Γt

wT ~t0dΓ︸︷︷︸feΓ

+

∫ΩwT~bdΩ︸︷︷︸feΩ

(7.17)

para todo w ∈ U0.

7.3.5 Utilização de aproximações lineares em elementos triangulares com três

nós

Na aplicação do MEF temos que as funções teste e as funções de ponderação em cada elementosão dadas por

u(x, y) ≈nel∑e=1

ue(x, y) =

nel∑e=1

N e(x, y)de (7.18)


e

w(x, y) ≈nel∑e=1

we(x, y) =

nel∑e=1

N e(x, y)we (7.19)

para todo (x, y) ∈ Ω, onde nel é o número de elementos triangulares no qual o Ω foi particionado.O elemento triangular de três nós é um elemento de deslocamento linear, onde as deformações

são constantes no elemento. Os nós, numerados no sentido anti-horário, possuem cada um delesduas componentes. Sendo assim, temos que o vetor N e(x, y) é dado por

(N e)T = [N e1 , N

e2 , N

e3 ]; (7.20)

onde

N e1 =

1

2Ae(xe2y

e3 − xe3ye2 + (y2

e − y3e)x+ (x3

e − x2e)y);

N e2 =

1

2Ae(xe3y

e1 − xe1ye3 + (y3

e − y1e)x+ (x1

e − x3e)y);

N e3 =

1

2Ae(xe1y

e2 − xe2ye1 + (y1

e − y2e)x+ (x2

e − x1e)y);

onde Ae é a área do elemento triangular. Desta forma, obtemos que o vetor e o campo de desloca-mentos são dados por:

de = [uex1, uey1, u

ex2, u

ey2, u

ex3, u

ey3]T

e

(uex, uey)T =

[N e

1 0 N e2 0 N e

3 00 N e

1 0 N e2 0 N e

3

]de

Figura 7.12: Componentes do deslocamento em cada nó do elemento triangular.

Utilizando o operador gradiente simétrico ∇S , obtemos a matriz de rigidez de cada elemento,dada por:

Ke =

∫Ωe

(Be)TDeBedΩ (7.21)

onde Be é dado por:

Be =1

2Ae

ye23 0 ye31 0 ye12 00 xe32 0 xe13 0 xe21

xe32 ye23 xe13 ye31 xe21 ye12

desendo xeIJ = xeI − xeJ e yeIJ = yeI − yeJ .

As matrizes de força de campo do elemento feΩ são calculadas numericamente por interpolação,utilizando-se das coordenadas triangulares. Porém, como estamos desprezando a ação gravitacionalsobre o corpo e nenhuma outra força de campo é aplicada sobre a chave, temos que em nosso modeloo vetor de força de campo é identicamente nulo em todos os elementos. Já os vetores de força de


contorno feΓ, para uma tração constante ~t0 = [tx, ty] são dados por:

feΓ =

∫ΓN eT~tdΓ =

l

2

txtytxty00

(7.22)

onde l é o comprimento da borda do contorno do elemento.Como as integrais são realizadas em cada um dos elementos, temos que a associação dos com-

ponentes de cada elemento em relação às componentes globais do problema gera o sistema linear:

Kd = f (7.23)

onde K é a matriz global de rigidez e f a soma das matrizes de força de contorno e de força decampo calculadas para cada um dos elementos.

Observando que a matriz de rigidez global é simétrica, positiva denida e esparsa podemosaplicar os algoritmos apresentados no Capítulo 6, obtendo então de maneira ecaz os valores dosdeslocamentos em cada um dos nós.


Dado que o material do qual a chave é composta possui módulo de Young E = 2× 107N/cm2 ecoeciente de Poisson ν = 0.3 (conforme dados obtidos em [Sch88]), iremos analisar a deformaçãoocorrida na chave quando aplicamos uma força constante de -200 N/cm na direção y. Considerandoque nenhuma carga é aplicada sobre a face relativa a espessura da chave, temos que a matriz D édada por:

D =E

1− ν2

1 ν 0ν 1 00 0 (1− ν)/2

(7.24)

Inicialmente particionamos o domínio em uma malha com 157 pontos e 248 elementos atravésdo EasyMesh, conforme mostram as Figuras 7.13, 7.14 e 7.15.

Figura 7.13: Particionamento do domínio em elementos triangulares com constante de renamento iguala 0.25.


Figura 7.14: Enumeração dos nós da malha com constante de renamento igual a 0.25.

Figura 7.15: Enumeração dos elementos da malha com constante de renamento igual a 0.25.

Resolvendo o sistema linear (7.23) obtivemos a distribuição do módulo do deslocamento emcada ponto da malha. Com isto, através de interpolação, obtemos aproximações do deslocamentoem todos os pontos do domínio Ω. Os resultados obtidos encontram-se representados nas Figuras7.16 e 7.17. Observe que o deslocamento é muito pequeno e por isso, na visualização dos resultadosé necessário o uso de um fator de escala de 104. Na tabela 7.4 temos os valores obtidos para omódulo do deslocamento em alguns pontos da chave e alguns valores de constantes de renamento.

Ponto (x,y) cte = 0.25 cte = 0.125 cte = 0.0625 cte = 0.03125(10.0,0.0) (-0.5588,-0.1718) (-0.589,-0.1876) (-0.6027,-0.195) (-0.6041,-0.1955)(10.0,30.0) (0.5552,-0.1721) (0.5888,-0.1876) (0.6027,-0.1951) (0.6048,-0.1961)(20.0,10.0) (-0.8023,-0.8298) (-0.8606,-0.8664) (-0.8908,-0.8755) (-0.8954,-0.8781)(20.0,20.0) (0.8001,-0.8311) (0.8613,-0.8674) (0.8875,-0.8746) (0.9005,-0.881)(30.0,10.0) (-2.3147,-4.5285) (-2.4826,-4.8552) (-2.5525,-4.9986) (-2.5571,-5.0165)(30.0,20.0) (2.3068,-4.5273) (2.4842,-4.8552) (2.5542,-4.9981) (2.561,-5.0158)(40.0,10.0) (-2.94,-10.2146) (-3.1465,-10.934) (-3.2263,-11.2297) (-3.23,-11.2577)(40.0,20.0) (2.9529,-10.2656) (3.1628,-10.9842) (3.2429,-11.2801) (3.2492,-11.3083)(50.0,10.0) (-3.0225,-16.3658) (-3.2306,-17.5046) (-3.3111,-17.9615) (-3.3155,-17.9998)(50.0,20.0) (3.0622,-16.4158) (3.2746,-17.5546) (3.3554,-18.0115) (3.3617,-18.0499)

Tabela 7.4: Deslocamento obtido em alguns pontos da chave com fator de escala 104.


Figura 7.16: Norma do deslocamento ocorrido em cada ponto do domínio utilizando-se da malha comconstante de renamento igual a 0.25.

Figura 7.17: Visualização do deslocamento da chave com um fator de escala 104 (malha com constante derenamento igual a 0.25).


Referências Bibliográcas

[Bat96] Klaus-Jurgen Bathe. Finite Element Procedures. Prentice-Hall, 1st edição, 1996. 36, 39

[Bra06] Dietrich Braess. Finite Elements - Theory, fast solvers and applications in solid mechanics.Cambridge, 3rd edição, 2006. 36

[CK08] Ward Cheney e David Kincaid. Numerical Mathematics and Computing. Thomson - BrooksCole, 6th edição, 2008. 25, 26

[DH03] Jean Donea e Antonio Huerta. Finite Element Methods for Flow Problems. Wiley, 1stedição, 2003. 36, 39

[FB07] Jacob Fish e Ted Belytscho. A rst course in nite elements. John Wiley & Sons, 1stedição, 2007. iii, xi, xii, 1, 2, 29, 39, 57, 61, 62

[Fla] Joseph E. Flaherty. Course notes - nite element analysis. http://www.cs.rpi.edu// aher-je/. xi, 44

[Goc06] Mark S. Gockenbach. Understanding and Implementing the Finite Element Method. SIAM,1st edição, 2006. 36, 39

[Gui02] Hamilton Luiz Guidorizzi. Um curso de Cálculo - Vol. 3. LTC, 5th edição, 2002. 36

[Saa03] Yousef Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2nd edição, 2003. 49,50, 52

[SB02] J Stoer e R. Bulirsch. Introduction to Numerical Analysis. Springer-Verlag, 3rd edição,2002. 50

[Sch88] Hans Rudolf Schwarz. Finite Element Methods. Academic Press, 1st edição, 1988. 65

[Shi] Carlos Y. Shigue. Material de aula - método dos elementos nitos.www.demar.eel.usp.br/metodos/mat_didatico.html. 1

69

Documents

Uma Introdução ao Método dos Elementos Finitos Marco