25

Tutorial Elementos Finitos

Embed Size (px)

Citation preview

Departamento de Matemática | Instituto de Ciências Exatas | Universidade Federal de Minas Gerais

Uma breve introdução ao

Método dos Elementos Finitos

Breno Loureiro Giacchini

Janeiro de 2012

Conteúdo

Prefácio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1 Introdução 2

2 O problema unidimensional 32.1 Formulação fraca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Discretização do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Existência e unicidade da solução do problema aproximado . . . . . . . . . . . . . . . 62.4 Um caso particular: partição regular do intervalo . . . . . . . . . . . . . . . . . . . . . 7

3 O problema bidimensional 83.1 Formulação fraca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2 Problema aproximado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.3 Uma base para Vd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.3.1 Enumerações dos vértices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.3.2 Funções da base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.3.3 Cálculo dos gradientes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.4 Alguns exemplos de malhas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.4.1 Exemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.4.2 Exemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.5 Outros casos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Bibliograa 23

Prefácio

Em meio a tantos bons livros e apostilas sobre o Método dos Elementos Finitos, o questionamento doporquê da escrita deste texto não é de todo descabido. O que nos motivou a escrevê-lo é a diculdadede se encontrar um texto, em português, que apresente o Método de forma simples e direta, que lheforneça uma idéia geral e ao mesmo tempo permita sua implementação em casos simples, mas semgrandes delongas em formalismos matemáticos e preâmbulos sobre análise funcional, por exemplo.

Admitimos, pois, que nossa exposição do Método não é feita com todo rigor matemático nem emtoda sua generalidade, mas cremos que essa opção satisfaz ao estudante que deseja entender sua essênciae aplicá-lo, de forma rápida, em alguns casos; ou ao interessado em ter um primeiro contato com essatécnica de resolver numericamente problemas de valores de contorno. Por ser um texto introdutório,que dá apenas um sabor do Método, nos limitamos contemplar o problema de Dirichlet homogêneouni e bidimensional e a utilizar, neste caso, apenas elementos triangulares.

Deixamos expresso nosso agradecimento ao apoio da Fundação de Amparo à Pesquisa do Estadode Minas Gerais FAPEMIG , que nanciou o projeto de pesquisa Obtenção dos autovalores e

autofunções do laplaciano via o quociente de Rayleigh, do qual o estudo do Método dos ElementosFinitos e a escrita deste texto foram partes integrantes. Também agradecemos ao professor RodneyJosué Biezuner, nosso orientador neste projeto.

1

Capítulo 1

Introdução

Diversos problemas da Física, Engenharia e outras ciências aparecem sob a forma de uma equaçãode Poisson

−∆u = f (x) em Ω (1.1)

com condição de fronteira de Dirichlet u = c sobre ∂Ω, sendo c uma função constante por partes.Aqui ∆ é o operador laplaciano1, Ω representa o aberto limitado no qual o problema está denido e∂Ω, sua fronteira. Quando c = 0 temos a condição de Dirichlet homogênea. Ao conjunto de umaequação de Poisson com uma condição de Dirichlet homogênea chamamos um problema de Dirichlethomogêneo:

−∆u = f(x) em Ω,u = 0 sobre ∂Ω.

(1.2)

Dependendo da geometria do domínio Ω a solução do problema pode ser obtida analiticamentena forma de séries de Fourier. Exemplos clássicos normalmente estudados num curso de equaçõesdiferenciais parciais são o caso de retângulos, semiplanos, discos e paralelepípedos. No entanto, épreciso recorrer a métodos numéricos caso o domínio se torne mais elaborado. O método dos elementosnitos (MEF) é conhecido por ser robusto e aplicável em domínios deveras elaborados. Essas tambémsão algumas de suas vantagens sobre o método das Diferenças Finitas, também bastante popular.

A idéia central do MEF é discretizar o domínio, representando-o, ainda que de forma aproximada,por uma reunião de um número nito de elementos; e resolver não o problema original (1.2), mas sim umque lhe é associado sua forma fraca. No caso de um domínio plano, os elementos podem ser triângulosou quadriláteros. O método pode ser utilizado para resolver não só problemas elípticos, como o hápouco mencionado; e as condições não necessitam ser de Dirichlet: o MEF também é aplicável no casode condição de Neumann ou Robin. Optamos por explorar neste texto apenas elementos triangularese considerar somente o problema (1.2), já que nosso objetivo é propiciar um primeiro contato com oMEF.

Analisaremos, primeiramente, o caso do problema unidimensional, que é bastante simples e útilcomo introdução ao método. Em seguida, passaremos ao problema bidimensional, apresentando eexemplicando como o MEF se lhe aplica. Cremos que a partir daí o leitor ou a leitora já estarão aptosa utilizar dessa ferramenta na resolução de alguns problemas de interesse.

1Se Ω ⊂ Rn e u função u : Ω→ R, ∆u =

n∑i=1

d2u

dxi2.

2

Capítulo 2

O problema unidimensional

2.1 Formulação fraca

O problema de Dirichlet homogêneo unidimensional se escreve−d

2u

dx2= f (x) em [0,1],

u(0) = u(1) = 0.

(2.1)

Assumiremos que a função f : [0, 1]→ R é limitada e contínua por partes. Isso é necessário porqueo MEF supõe que f é integrável. Notamos que aqui Ω = [0, 1].

Ao invés de resolver o problema (2.1) da forma como está escrito, o MEF se propõe a solucionarum problema equivalente, chamado formulação fraca do original. Para escrever (2.1) dessa forma,

principiamos denindo o espaço de funções V = v; v é função contínua em [0, 1],dv

dxé limitada e

contínua por partes, e v(0) = v(1) = 0.Em seguida, multiplicamos a primeira equação de (2.1) por uma função qualquer de V e integramos

a equação resultante em Ω:

−1∫

0

d2u

dx2· v dx =

1∫0

f(x) · v dx .

Integrando por partes e lembrando que v sastifaz a condição de Dirichlet homogênea:

−(vdu

dx

) ∣∣∣∣∣1

0

+

1∫0

du

dx

dv

dxdx =

1∫0

f v dx

⇒1∫

0

du

dx

dv

dxdx =

1∫0

f v dx (2.2)

para todo v ∈ V . A equação (2.2) juntamente com a condição de Dirichlet homogênea é a formulaçãofraca do problema (2.1).

Mostraremos agora que a existência de uma solução de (2.1), problema original, implica na equi-valência entre os problemas de formulação forte e fraca. Vimos logo acima que formulação forte ⇒formução fraca. Resta, pois, vericar a recíproca.

Já que supomos que o problema original tem solução, sabemos qued2u

dx2existe e é contínua por

partes. Podemos, então, integrar por partes a formulação fraca (2.2). O que obtemos é justamente aformulação forte (original):

3

1∫0

du

dx

dv

dxdx =

(vdu

dx

) ∣∣∣∣∣1

0

−1∫

0

vd2u

dx2=

1∫0

f(x) · v(x) dx

⇒1∫

0

(f(x) +

d2u

dx2

)v(x) dx = 0 , ∀ v ∈ V.

Como a igualdade acima se verica para qualquer função v em V , o termo do integrando que está entreparêntesis deve ser nulo; e assim chegamos ao problema original:

f(x) +d2u

dx2= 0 , 0 < x < 1.

Com isso provamos que uma função que resolve o problema forte também é solução do fraco; e que,se a solução do problema fraco for sucientemente regular, ela também resolverá o problema forte. Nométodo dos elementos nitos resolveremos o problema fraco, (2.2).

2.2 Discretização do problema

Na formulação original, o problema de Dirichlet é contínuo e seu espaço de soluções pode ter dimensãoinnita1. Aproximaremos o problema contínuo por outro discreto, cuja solução está em um espaço dedimensão nita. Isso é feito dividindo o domínio (o intervalo [0, 1]) em um número nito de subintervalosIj = [xj−1, xj ], 1 6 j 6 N + 1, N ∈ N, com 0 = x0 < x1 < · · · < xN < xN+1 = 1. Cada subintervalotem comprimento hj = xj − xj−1.

Essa discretização é uma partição [do intervalo], cuja norma denimos h = maxjhj, o compri-mento do maior dos subintervalos.

O termo discretização é usado justamente porque passamos de um contínuo (a função originalestá denida num domínio que é uma reunião não-enumerável de pontos) para um conjunto discreto: odomínio passa a ser uma reunião nita de intervalos2. Em cada um desses intervalos Ij , aproximamosa função original u por um segmento de reta de extremos u(xj−1) e u(xj) (Figura 1a). Evidentemente,quanto menor o comprimento dos subintervalos, ou seja, quanto menor a norma da partição, maisa função discretizada ud se aproximará da original u (Figura 1a-b). Notemos, ainda, que ud comodenida é contínua.

Figura 1: Aproximação de uma função suave por outra linear por partes. Quanto menor a norma dapartição, melhor a aproximação.

1Exemplos de problemas cujas soluções estão num espaço de dimensão innita são aqueles tradicionalmente estudadosnum curso de equações diferenciais parciais que têm como resultado séries de Fourier innitas. Cada uma daquelas funções

sen(nπxL

)ou cos

(nπxL

)é uma função da base do espaço que contém a solução do problema. Como existem innitos

n ∈ N, a base é um conjunto innito.2Uma reunião nita de intervalos ainda é um conjunto não-enumerável de pontos. A idéia aqui é que, ao invés de

buscarmos u com denição ponto a ponto, vamos aproximá-la por uma função que é denida subintervalo por subintervaloe, nesse sentido, ela será denida discretamente. Mais adiante no texto cará clara essa idéia.

4

Para discretizar o problema na forma fraca, devemos também aproximar o espaço V por um dedimensão nita, Vd = v; v é contínua em [0, 1], v é linear em cada Ij e v(0) = v(1) = 0. Notemosque Vd ⊂ V de sorte que ao tomarmos uma função v ∈ Vd não ferimos a condição v ∈ V da formulaçãofraca.

Nosso problema discretizado (ou aproximado) é, então, encontrar ud ∈ Vd tal que

1∫0

duddx

dv

dxdx =

1∫0

f(x)v(x) dx ∀ vd ∈ Vd . (2.3)

Esse ud será a aproximação para a função u desejada.Observação 1: note que a condição de fronteira u(0) = u(1) = 0 está contida no enunciado do pro-blema discretizado (2.3) já que ud ∈ Vd implica na condição de Dirichlet homogênea.Observação 2: talvez tenha parecido ao leitor que nos precipitamos ao declarar que a solução apro-ximada ud que buscamos está em Vd. Veremos porque podemos assumir isso. A primeira condiçãoque ud deve satisfazer para pertencer a Vd é ser contínua. Já vimos que como foi denida ud ≈ u écontínua. A segunda condição (linearidade em cada subintervalo) vem também da discretização doproblema e está relacionada com a qualidade dessa aproximação - aproximamos uma curva suave poroutra poligonal. Ao assumirmos que buscamos uma solução com essa aproximação, ud satisfaz, porconseguinte, a segunda condição. Por m, a terceira é justamente a condição de Dirichlet que tantou quando ud satisfazem por hipótese. Concluímos assim que ud ∈ Vd.

Ao discretizarmos o espaço V , o aproximamos por um de dimensão nita. Como v ∈ Vd é linearem cada Ij , as funções

φj(x) =

x− xj−1

hjse x ∈ [xj−1, xj ],

xj+1 − xhj+1

se x ∈ [xj , xj+1],

0 caso contrário.

(2.4)

formam uma base, B, de Vd.

Figura 2: Gráco da função φj de B.

Mostramos na Figura 2 um gráco de uma dessas funções-base de Vd. É fácil ver que se i 6= j, asfunções φi e φj são linearmente independentes. Um momento de reexão bastará para que a leitora seconvença de que qualquer função de Vd se escreve em termos das φj acima denidas.

Ora, como ud pertence a Vd, será da forma

ud(x) =

N∑j=1

αj φj(x) , x ∈ [0, 1], (2.5)

e o nosso problema (2.3) se escreverá

1∫0

d

dx

(∑j

αjφj

)dv

dxdx =

1∫0

f(x)v(x) dx ∀ vd ∈ Vd . (2.6)

5

Recordemos que a função v é uma função qualquer de Vd. Escolhemos, então, v como sendo umadas funções da base: v = φi para algum i 6 N . Para esse i, (2.6) implica em

1∫0

(N∑j=1

αjdφjdx

)dφidx

dx =

1∫0

f(x)φi(x) dx

⇒N∑j=1

(αj

1∫0

dφjdx

dφidx

dx

)=

1∫0

f(x)φi(x) dx (2.7)

Variando i de 1 a N, (2.7) resulta em um sistema de N equações e N incógnitas αj :

1∫0

dφ1

dx

dφ1

dxdx

1∫0

dφ2

dx

dφ1

dxdx · · ·

1∫0

dφNdx

dφ1

dxdx

1∫0

dφ1

dx

dφ2

dxdx

1∫0

dφ2

dx

dφ2

dxdx · · ·

1∫0

dφNdx

dφ2

dxdx

......

. . ....

1∫0

dφ1

dx

dφNdx

dx1∫0

dφ2

dx

dφNdx

dx · · ·1∫0

dφNdx

dφNdx

dx

α1

α2...αN

=

1∫0

f(x)φ1(x) dx

1∫0

f(x)φ2(x) dx

...1∫0

f(x)φN (x) dx

(2.8)

Chamaremos a matriz do sistema acima deM e seus elementos denotaremos pormij . M é amatrizde rigidez, enquanto que o vetor que aparece no membro direito de (2.8) é denominado vetor decarga.

Vimos, portanto, que o problema de achar a função ud ∈ Vd que satisfaz (2.3) se reduz à resoluçãode um sistema linear. Resolvendo-o, determinamos os coecientes αj e podemos construir, usando (2.5)e (2.4), a função ud ≈ u. Antes de darmos o assunto por encerrado, vericaremos algumas propriedadesda matriz M e provaremos que sempre existirá uma (única) solução para o sistema (2.8) - e, portanto,para o problema aproximado (2.3).

2.3 Existência e unicidade da solução do problema aproximado

Proposição 1. A matriz M goza das seguintes propriedades:

R1) é simétrica;

R2) é tridiagonal;

R3) é positiva denida - isto é, wTMw > 0 ∀ w não-nulo em RN .

Demonstração.R1) É conseqüência da comutatividade do produto de funções:

mij =

1∫0

dφjdx

dφidx

dx =

1∫0

dφidx

dφjdx

dx = mji .

R2) Calcularemos os elementos mij para mostrar que M é tridiagonal. Para tanto, usaremos asderivadas das funções φj denidas por (2.4).

Se i = j, mii =

xi∫xi−1

1

h2i

dx+

xi+1∫xi

1

h2i+1

dx =1

hi+

1

hi+1.

Se i e j diferem por apenas 1 unidade, mi,i−1 = mi−1,i =

xi∫xi−1

−1

hi

1

hidx = −xi − xi−1

h2i

= − 1

hi.

Finalmente, se i e j diferem por mais de 1 unidade, mij = 0.Conclui-se, pois, que apenas os termos da forma mi,i e mi,i±1 são não-nulos.

6

R3) Ora, wTMw =

N∑i=1

N∑j=1

wj

( 1∫0

dφjdx

dφidx

dx

)wi =

1∫0

[(N∑j=1

wjdφjdx

)(N∑i=1

widφidx

)]dx =

=

1∫0

(N∑j=1

wjdφjdx

)2

dx > 0.

Como essa relação vale para qualquer vetor w em RN − 0, a igualdade só se vericaria casodφjdx

= 0 para cada j e em todo o intervalo [0, 1]. Como tal situação não ocorre, temos a desigualdade

estrita wTMw > 0 ∀ w ∈ RN − 0.

Um conhecido teorema da Álgebra Linear garante que uma matriz positiva denida tem deter-minante não-nulo3. Outro teorema reza que que se a matriz de um sistema linear tem determinantenão-nulo, o sistema tem solução única. Esses teoremas, juntamente com o terceiro item da Proposição1, nos asseguram que (2.8) tem solução - e ela é única.

2.4 Um caso particular: partição regular do intervalo

Concluimos o estudo do caso unidimensional escrevendo o sistema (2.8) num caso particular de partiçãodo intervalo, a saber, considerando que todos os subintervalos Ij têm mesmo comprimento h. A umapartição deste tipo dá-se o nome de regular.

Utilizando os cálculos realizados na prova do segundo item da Proposição 1, temos que uma partiçãoregular do domínio fornece a matriz de rigidez:

M =1

h

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

E o sistema (2.8) pode ser escrito como:

2 −1−1 2 −1

−1. . . . . .. . . . . . −1

−1 2 −1−1 2

α1

α2

...

αN

= h

1∫0

f(x)φ1(x) dx

1∫0

f(x)φ2(x) dx

...1∫0

f(x)φN (x) dx

(2.9)

O leitor que já estudou o Método das Diferenças Finitas notará que a matriz M obtida para umapartição regular é bastante semelhante à encontrada naquele método - elas só diferem por um termo 1/hmultiplicando. Dependendo da maneira de discretizar as integrais do vetor de carga, os dois métodoscoincidirão.

3Mais precisamente, seu determinante é maior que zero.

7

Capítulo 3

O problema bidimensional

3.1 Formulação fraca

Sejam Ω ⊂ R2 um aberto limitado e f uma função real contínua por partes e limitada, em Ω. Oproblema de Dirichlet homogêneo bidimensional se escreve

−∆u = f(x, y) em Ω,u = 0 sobre ∂Ω.

(3.1)

Assim como zemos no caso unidimensional, escreveremos o problema (3.1) na forma fraca.

Denimos o espaço de funções V = v : R2 → R; v é função contínua em Ω,∂v

∂xe∂v

∂ysão contínuas

por partes em Ω, e v = 0 sobre ∂Ω. Multiplicando a equação de Poisson do problema (3.1) por umafunção qualquer de V e integrando sobre Ω temos:

−v ·∆u = v · f ⇒ −∫Ω

v ·∆u dV =

∫Ω

v · f dV. (3.2)

Podemos reescrever a equação acima de forma mais conveniente usando a fórmula de Green, quese baseia noTeorema 1. [Teorema do divergente] Seja Ω ⊂ Rn compacto e com fronteira suave por partes.

Se w é um campo de vetores diferenciável denido em Ω, então:∫Ω

div w dV =

∫∂Ω

〈w , n〉 ds,

onde n representa o vetor unitário normal à ∂Ω.A notação1 〈w , n〉 indica o produto escalar dos vetores w e n.Acreditamos que esse teorema, cuja prova omitiremos, já foi estudado pela leitora em algum curso

de Cálculo, pelo menos para n = 2 e n = 3. O caso n = 2 que nos interessa é às vezes chamadoTeorema de Green. Para obter a fórmula de Green, aplicamos o Teorema para os campos de

vetores a(x, y) =

(g · ∂h

∂x, 0

)e b(x, y) =

(0 , g · ∂h

∂y

), sendo as funções g,h : R2 → R. Considerando

que o vetor normal unitário é n = (n1 , n2), temos, para a,∫Ω

(g · ∂

2h

∂x2+∂g

∂x

∂h

∂x

)dV =

∫∂Ω

g · ∂h∂x· n1 ds (3.3)

1Outras notações e denições que utilizamos ao longo deste texto:o divergente de um campo de vetores w: Ω→ Rn: div w =

∑nk=1 ∂wk/∂xk;

o gradiente de uma função f : Rn → R: gradf =

(∂f

∂x1, · · · , ∂f

∂xn

).

8

e, para b: ∫Ω

(g · ∂

2h

∂y2+∂g

∂y

∂h

∂y

)dV =

∫∂Ω

g · ∂h∂y· n2 ds. (3.4)

Somando membro a membro as equações (3.3) e (3.4) e reagrupando:∫Ω

[g ·(∂2h

∂x2+∂2h

∂y2

)+∂g

∂x

∂h

∂x+∂g

∂y

∂h

∂y

]dV =

∫∂Ω

g ·(n1 ·

∂h

∂x+ n2 ·

∂h

∂y

)ds

⇒∫Ω

(g ·∆h+ 〈gradg , gradh〉

)dV =

∫∂Ω

g · 〈n , gradh〉 ds. (Fórmula de Green) (3.5)

Notemos que se a função g(x, y) acima satisfaz a condição de Dirichlet homogênea, a integral sobre∂Ω em (3.5) é nula e a fórmula de Green implica em

−∫Ω

g ·∆h dV =

∫Ω

〈gradg , gradh〉 dV.

Comparando (3.2) com a equação acima, vemos que os membros esquerdos são iguais se zermosg = v e h = u. Temos, portanto, que∫

Ω

v · f dV =

∫Ω

〈gradv , gradu〉 dV ∀ v ∈ V. (3.6)

Esta equação acrescida da condição de Dirichlet homogênea formam a formulação fraca do problemabidimensional (3.1). É possível mostrar, como o zemos no caso unidimensional, que as formas fracae forte são equivalentes e que uma solução da forma fraca, se sucientemente regular, também serásolução da forma forte.

3.2 Problema aproximado

Uma vez compreendida a essência do método no caso unidimensional, o caso dos domínios planosnão apresenta maiores diculdades no que tange essa essência. A diculdade surge no momento dediscretizar Ω e trabalhar com a malha resultante, como veremos em seções seguintes.

Começaremos a discretização do problema dividindo o domínio Ω em triângulos. Obviamente, nãoé qualquer domínio que aceita essa divisão perto de sua borda. Neste caso, aproximamos Ω por Ωd cujafronteira é uma curva poligonal (formada por uniões nitas de segmentos de retas). A cada um dessestriângulos chamamos elemento. A discretização em triângulos (ou triangulação) deve cumprir asseguintes condições:D1) A reunião de todos os elementos forma Ωd, que aproxima Ω;

D2) Os elementos não se sobrepõem;

D3) Os vértices de um elemento nunca ocorrem no lado de outro elemento.

A Figura 3 mostra exemplos de triangulações permitidas e não permitidas no método dos elementosnitos, além de ilustrar como podemos fazer a aproximação da fronteira.

Figura 3: (a) exemplo de triangulação permitida. A partição em (b) não é permitida pois o traço em azuldene um vértice que ocorre em um lado de outro elemento.

9

No problema discretizado, buscamos uma função ud que aproxima u. Aproveitamos nossa malha(domínio discretizado) para impor uma condição sobre ud: que ela seja contínua e linear em cadaelemento. Esta última imposição signica que o gráco de ud em cada elemento é um pedaço de planocontido no R3. Ao fazermos essas exigências estamos apenas escolhendo como iremos aproximar u.

Chamamos atenção para o fato de que ao aproximar Ω por Ωd nos propomos a resolver o problema(3.1) com o domínio Ωd. Portanto, quanto mais parecido for ∂Ωd de ∂Ω, mais a função encontrada udserá parecida com a função real u.

Para discretizar o problema na forma fraca, devemos ainda aproximar o espaço V por um nitoVd = v; v é contínua em Ωd, v é linear em cada elemento e v = 0 sobre ∂Ωd. Como no casounidimensional, v ∈ Vd ⇒ v ∈ V , implicando que v satisfaz a formulação fraca contínua.

O problema aproximado é, então: achar ud ∈ Vd tal que∫Ωd

〈gradv , gradud〉 dV =

∫Ωd

v · f dV ∀ v ∈ Vd. (3.7)

A maior diculdade que surge no problema bidimensional é a manipulação numérica das funçõesda base de Vd. Por isso, optamos por primeiro expor a teoria supondo que temos uma base - mas semescrevê-la - e obter o sistema linear resultante da discretização. Mostraremos ainda que, assim comono caso unidimensional, o sistema tem única solução. Na seção seguinte escreveremos explicitamenteuma base e faremos algumas contas, já com vistas à implementação de um algoritmo de MEF.

Seja, pois, B uma base do espaço Vd. Sabemos que ud, por estar em Vd, tem (única) representaçãocomo combinação linear das funções de B. Denotando essas funções por φj , escrevemos

ud(x, y) =

N∑j=1

αj φj(x, y) , (x, y) ∈ Ωd,

onde N é a dimensão de Vd.Substituindo ud acima no problema discretizado (3.7) obtemos

N∑j=1

αj

∫Ωd

〈gradv , gradφj〉 dV =

∫Ωd

v · f dV ∀ v ∈ Vd.

(Note que o integrando está dentro do somatório.)Em particular, para v = φi qualquer da base:

N∑j=1

αj

∫Ωd

〈gradφi , gradφj〉 dV =

∫Ωd

φi · f dV (3.8)

Variando i de 1 a N , (3.8) se mostra um sistema linear N × N de incógnitas αj . DenindoMij =

∫Ωd

〈gradφi , gradφj〉 dV , o nosso problema equivale ao sistema

M11 · · · M1N...

. . ....

MN1 · · · MNN

α1

...αN

=

Ωd

φ1 · f dV

...∫Ωd

φN · f dV

(3.9)

.Vejamos algumas propiedades da matriz M do sistema, a matriz de rigidez.

Proposição 2. A matriz M é simétrica e positiva denida.

Demonstração.Como o produto interno é simétrico por denição,Mij =

∫Ωd

〈gradφi , gradφj〉 dV =∫

Ωd

〈gradφj , gradφi〉 dV = Mji ⇒ M é simétrica.

10

Para mostrar que M é positiva denida temos que provar que ∀ w ∈ RN −~0 se tem wTMw > 0.Ora,

wTMw =

N∑i=1

wi

N∑j=1

wj

∫Ωd

〈gradφi , gradφj〉 dV =

∫Ωd

⟨ N∑i=1

wi · gradφi ,N∑j=1

wj · gradφj⟩dV > 0

pois o produto escalar de um vetor por ele mesmo é sempre > 0, sendo que a igualdade só se verica casoo vetor seja nulo. Como aqui w é um vetor qualquer, isso ocorreria somente se tivéssemos gradφi = 0 ∀ ie em todo Ωd. Como por hipótese φi é contínua e cumpre a condição de fronteira de Dirichlet, issoimplicaria que as funções da base são identicamente nulas, um absurdo. Logo, wTMw é estritamentepositivo, denindo M postivamente.

Como conseqüência do fato deM ser positiva denida, temos que o sistema (3.9) sempre admite umaúnica solução, os coecientes αj que, juntamente com as funções da base, determinam a aproximaçãoud de u.

Agora que vimos como o problema de Dirichlet se escreve na forma discreta usando o método doselementos nitos, usaremos uma base de Vd para estudar como fazer os cálculos e efetivamente resolvero problema.

3.3 Uma base para Vd

3.3.1 Enumerações dos vértices

Antes de buscarmos uma base para o espaço Vd, introduziremos alguns conceitos que se mostrarãoimportantes naquela tarefa. Denimos o número de vértices da malha como sendo o número totalde vértices dos elementos, com a condição de que mesmo se determinado vértice é comum a várioselementos, o contamos apenas uma vez. Chamamos de vértices interiores aqueles que não estão sobrea fronteira de Ωd (Figura 4).

Figura 4: Uma malha com 16 elementos e 13 vértices, sendo que destes, 5 são interiores.

Vamos supor que o nosso domínio Ω foi dividido em m triângulos, resultando em uma malha de Nvértices, sendo que N são interiores. Faremos algumas denições.Denição 1. Chamamos de enumeração dos elementos a uma bijeção que associa a cada elemento

triangular da malha um número natural entre 1 e m. Representamos, pois, cada elemento pela letra Tseguida de seu número como sub-índice. Por exemplo, Tk é o k-ésimo elemento da malha.

Denição 2. Chamamos de enumeração global dos vértices interiores a uma bijeção que associa a

cada vértice interior da malha um número natural entre 1 e N . Representamos, pois, cada vértice

interior pela letra p seguida de seu número como sub-índice. Por exemplo, pi é o vértice interior i damalha.

Denição 3. Chamamos de enumeração global dos vértices a uma bijeção que associa a cada vértice

da malha um número natural entre 1 e N , respeitando a enumeração global dos vértices interiores. Esta

enumeração consiste em adotar a enumeração da Denição 2 e ainda atribuir números entre N + 1 e

N aos vértices da fronteira de Ωd.

Denição 4. Chamamos de enumeração local dos vértices a uma bijeção que

i) associa a cada vértice de um elemento Tk um número do conjunto 1, 2, 3;

11

ii) percorre o elemento em sentido anti-horário. Isto é, denido o vértice número 1 do elemento Tk,percorre-se a fronteira do elemento em sentido anti-horário a partir desse vértice 1. O próximo vértice

será o de número 2 e o último será o número 3. O vértice s, s em 1, 2, 3, do elemento Tk tem

coordenadas (x(k)s , y

(k)s ).

Notemos que um mesmo vértice comum a dois elementos pode ter numeração local diferente emcada elemento. Por exemplo, pode ser o vértice 1 do elemento Tk e o vértice 3 do Tl. Neste caso(x

(k)1 , y

(k)1 ) e (x

(l)3 , y

(l)3 ) representam o mesmo ponto da malha. Supondo ainda que esse vértice seja o

vértice interior h global, então (x(k)1 , y

(k)1 ) = (x

(l)3 , y

(l)3 ) = ph.

Na Figura 5 mostramos um exemplo de uma malha e de uma possível enumeração (global) doselementos e dos vértices. A Tabela 1 complementa a Figura 5 exemplicando uma enumeração localdos vértices. Note que, para cada elemento, um dos vértices globais assume a posição local 1, 2 ou 3.Na Figura 6 mostramos alguns elementos e a enumeração local de seus vértices.

Figura 5: Exemplo de enumeração dos elementos (a), e enumeração global dos vértices (b).

Tabela 1: Exemplo de enumeração local dos vérticesElemento Vértice 1 Vértice 2 Vértice 3

1 16 5 12 1 5 63 1 6 24 6 7 25 2 7 96 7 8 97 15 16 38 16 1 39 1 4 310 4 1 211 4 2 1012 10 2 913 13 14 1514 13 15 315 13 3 1216 4 12 317 4 11 1218 10 11 4

Figura 6: Alguns elementos da malha daFigura 5 e exemplo de enumeração localde seus vértices. O único cuidado nessaenumeração é que seu sentido seja anti-horário.

3.3.2 Funções da base

Com os conceitos de enumerações globais e local dos vértices bem estabelecidos, podemos principiarnossa busca por uma base de Vd. Como discretizamos o problema de modo que v ∈ Vd fosse linear em

12

cada elemento, as funções φj : R2 → R tais que

φj(pi) =

1 se i = j,0 se i 6= j

(3.10)

e

gráco de φj no elemento Tk =

plano 6= 0 se Tk tem o vértice pj ,0 caso contrário

(3.11)

formam uma base B de Vd. Lembramos a notação: os vértices interiores estamos representando pelospontos pi. Tanto i quanto j acima assumem valores em 1, 2, 3, ..., N. Essas funções têm formatospiramidais, como se vê na Figura 7. Diremos que a função φj e o vértice pj são associados. Note queuma função é associada a um único ponto, e vice-versa.

Figura 7: Função chapéu.

As únicas funções de B que assumem valores não-nulos em um dado elemento são aquelas três asso-ciadas aos seus vértices. Diremos que essas funções são associadas ao elemento, que reciprocamentelhes é associado. Note que a cada elemento podem existir no máximo três funções associadas, mas queuma função pode ser associada a um número qualquer de elementos, dependendo da triangulação damalha.

Por (3.10) já sabemos quanto vale φj(x, y) se o ponto (x, y) for um vértice de um elemento. Nossatarefa agora é determinar o valor que a função assume num ponto no interior de um triângulo. Isto é,determinar φj(x, y) para qualquer (x, y) ∈ Ωd.

Já sabemos por (3.11) que se o ponto (x, y) não pertence a nenhum elemento que tenha o vérticepj , então φj(x, y) = 0. Vejamos, então, o que ocorre se (x, y) ∈ Tk e o triângulo Tk for associado àfunção φj .

Como vimos, existe uma enumeração local dos vértices de Tk: (x(k)1 , y

(k)1 ), (x

(k)2 , y

(k)2 ) e (x

(k)3 , y

(k)3 ).

Vamos supor, então, que a função φj ∈ B é a associada ao vértice 1, (x(k)1 , y

(k)1 ), do elemento

Tk. Determinaremos o valor de φj(x, y) se (x, y) ∈ Tk. Como sabemos que estamos no elementoTk, podemos dispensar os índices superiores nas coordenadas dos vértices, escrevendo simplesmente,(x1, y1). (Faremos isso apenas para deixar a notação mais limpa durante a dedução da fórmula; ao nalrestituiremos os índices superiores.) Por (3.10) sabemos que φj(x1, y1) = 1 e φj(x2, y2) = φj(x3, y3) =0. Por sua vez, (3.11) implica que se (x, y) ∈ Tk, então (x, y, φj(x, y)) está no plano determinadopelos pontos (x1, y1, φj(x1, y1)), (x2, y2, φj(x2, y2)) e (x3, y3, φj(x3, y3)), ou, substituindo os valores dafunção nos vértices, (x1, y1, 1), (x2, y2, 0) e (x3, y3, 0).

Para que o gráco de φj(x, y) esteja nesse plano, os três vetores que ligam (x, y, φj(x, y)) a cadaum dos pontos (x1, y1, 1), (x2, y2, 0) e (x3, y3, 0) devem ser coplanares ou, equivalentemente, o produto

13

misto2 dos três deve ser nulo. Esses vetores são:

(x, y, φj)− (x1, y1, 1) = (x− x1 , y − y1 , φj − 1)

(x, y, φj)− (x2, y2, 0) = (x− x2 , y − y2 , φj)

(x, y, φj)− (x3, y3, 0) = (x− x3 , y − y3 , φj),

onde, φj = φj(x, y). Igualando o produto misto a zero,∣∣∣∣∣∣x− x1 y − y1 φj − 1x− x2 y − y2 φjx− x3 y − y3 φj

∣∣∣∣∣∣ = 0.

Desenvolvendo o determinante em cofatores com relação à terceira coluna:

(φj − 1)

∣∣∣∣x− x2 y − y2

x− x3 y − y3

∣∣∣∣− φj ∣∣∣∣x− x1 y − y1

x− x3 y − y3

∣∣∣∣+ φj

∣∣∣∣x− x1 y − y1

x− x2 y − y2

∣∣∣∣ = 0

⇒ φj(x, y) =

∣∣∣∣x− x2 y − y2

x− x3 y − y3

∣∣∣∣∣∣∣∣x− x2 y − y2

x− x3 y − y3

∣∣∣∣− ∣∣∣∣x− x1 y − y1

x− x3 y − y3

∣∣∣∣+

∣∣∣∣x− x1 y − y1

x− x2 y − y2

∣∣∣∣ . (3.12)

Faremos algumas manipulações algébricas usando propriedades dos determinantes para escrever aequação acima de maneira mais conveniente. Por exemplo, no numerador:

∣∣∣∣x− x2 y − y2

x− x3 y − y3

∣∣∣∣ =

∣∣∣∣x− x2 y − y2

x y

∣∣∣∣−∣∣∣∣x− x2 y − y2

x3 y3

∣∣∣∣ =

∣∣∣∣x yx y

∣∣∣∣−∣∣∣∣x2 y2

x y

∣∣∣∣−∣∣∣∣ x yx3 y3

∣∣∣∣+∣∣∣∣x2 y2

x3 y3

∣∣∣∣ =

∣∣∣∣∣∣1 x y1 x2 y2

1 x3 y3

∣∣∣∣∣∣ .A última igualdade pode ser facilmente vericada desenvolvendo o determinante 3 × 3 em cofatoresrelativos à primeira coluna.

Realizando os mesmos passos que zemos com o determinante do numerador nos outros dois de-terminantes do denominador, a leitora é convidada a mostrar que (3.12) equivale à:

φj(x, y) =

∣∣∣∣∣∣∣1 x y

1 x(k)2 y

(k)2

1 x(k)3 y

(k)3

∣∣∣∣∣∣∣∣∣∣∣∣∣∣1 x

(k)1 y

(k)1

1 x(k)2 y

(k)2

1 x(k)3 y

(k)3

∣∣∣∣∣∣∣, (3.13)

onde foram restituídos os índices superiores.Um resultado da Geometria Analítica informa que o determinante do denominador acima é justa-

mente o dobro da área do triângulo de vértices (x(k)1 , y

(k)1 ), (x

(k)2 , y

(k)2 ) e (x

(k)3 , y

(k)3 ), ou seja, é o dobro

da área Ak do elemento Tk. Reescrevemos, pois, (3.13) como

φj(x, y) =1

2Ak

∣∣∣∣∣∣∣1 x y

1 x(k)2 y

(k)2

1 x(k)3 y

(k)3

∣∣∣∣∣∣∣ . (3.14)

2Denotando o produto vetorial entre dois vetores por × e o escalar por 〈 , 〉, o produto misto de três vetores a, b ec (nesta ordem) pertencentes a R3 é denido por 〈a , b× c〉 e pode ser calculado como o determinante da matriz cujaslinhas são a, b, e c, nesta ordem. A interpretação geométrica desse produto é o volume do paralelepípedo determinadopelos três vetores. Caso o resultado seja nulo, os três vetores não determinam volume algum, estando, pois, num mesmoplano.

14

A expressão acima fornece o valor de φj num ponto (x, y) qualquer do elemento Tk. Em outroselementos associados, a função pode não ser dada por (3.14). Relembremos as suposições feitas queresultaram em (3.14): consideramos que φj era associada à Tk, mais precisamente, era associada ao

vértice 1 de Tk. Sob essas hipóteses, encontramos a fórmula acima para φj neste elemento.Neste ponto deverá o leitor estar se perguntando o que ocorreria se a função φj fosse associada a

outro vértice r, (x(k)r , y

(k)r ), local que não o de número 1. Neste caso, o vértice r é o que cumpre o

papel de vértice 1 na equação acima. O vértice seguinte a r cumprirá o papel de vértice 2 e o último,o de vértice 3. É imporante que nos lembremos que a numeração local dos vértices sempre é feita emsentido anti-horário: a ordem local dos vértices é sempre 1 − 2 − 3 − 1 − 2 − 3 − 1 − · · · (o vértice 1sempre segue ao 3; o 3 segue ao 2, que segue ao 1).

Por exemplo, se φj é associada ao vértice 2 do elemento Tk, entãovértice 2 → vértice 1 em (3.14);vértice 3 → vértice 2 em (3.14) evértice 1 → vértice 3 em (3.14),onde → signica cumpre o papel de ou corresponde ao. Isso resulta, pois, em

φj(x, y) =1

2Ak

∣∣∣∣∣∣∣1 x y

1 x(k)3 y

(k)3

1 x(k)1 y

(k)1

∣∣∣∣∣∣∣ ∀(x, y) ∈ Tk.

Pode-se fazer o mesmo procedimento para a associação ao vértice 3 e assim chegamos numa ex-pressão geral para o valor de uma função qualquer φj de B em um elemento qualquer Tk:

φj(x, y) =

1

2Ak

∣∣∣∣∣∣∣1 x y

1 x(k)2 y

(k)2

1 x(k)3 y

(k)3

∣∣∣∣∣∣∣ se (x, y) ∈ Tk e φj for associada ao vértice 1 de Tk,

1

2Ak

∣∣∣∣∣∣∣1 x y

1 x(k)3 y

(k)3

1 x(k)1 y

(k)1

∣∣∣∣∣∣∣ se (x, y) ∈ Tk e φj for associada ao vértice 2 de Tk,

1

2Ak

∣∣∣∣∣∣∣1 x y

1 x(k)1 y

(k)1

1 x(k)2 y

(k)2

∣∣∣∣∣∣∣ se (x, y) ∈ Tk e φj for associada ao vértice 3 de Tk,

0 se (x, y) ∈ Tk mas φj não for associada a Tk.

(3.15)

Note o caráter local da expressão acima: φj é denida elemento por elemento. Às vezes, para

deixar bem claro que estamos calculando a função restrita ao elemento Tk, escreveremos φ(k)j . Em

cada elemento φj poderá será dada por uma expressão diferente, dependendo das coordenadas dos seusvértices e do número local do vértice associado. Da mesma forma, o gradiente de φj dependerá doelemento considerado.

Lembremos que as entradas da matriz de rigidez dependem dos gradientes das funções da base.Veremos agora como, a partir de (3.15), podemos calculá-los.

3.3.3 Cálculo dos gradientes

Supondo que φj é associada ao vértice 1 de Tk, utilizaremos o primeiro caso de (3.15). Desenvolvendoo determinante em termos dos cofatores da primeira linha, temos que

15

φ(k)j (x, y) =

1

2Ak

(∣∣∣∣∣x(k)2 y

(k)2

x(k)3 y

(k)3

∣∣∣∣∣− x∣∣∣∣∣1 y

(k)2

1 y(k)3

∣∣∣∣∣+ y

∣∣∣∣∣1 x(k)2

1 x(k)3

∣∣∣∣∣)

⇒∂φ

(k)j

∂x(x, y) = − 1

2Ak

∣∣∣∣∣1 y(k)2

1 y(k)3

∣∣∣∣∣ e∂φ

(k)j

∂y(x, y) =

1

2Ak

∣∣∣∣∣1 x(k)2

1 x(k)3

∣∣∣∣∣⇒ gradφ(k)

j (x, y) =1

2Ak

(y

(k)2 − y(k)

3 , x(k)3 − x

(k)2

). (3.16)

Caso a função φj seja a associada ao vértice 2 de Tk, basta fazer a troca de índices na expressãoacima. O vértice seguinte ao associado será o 3 (posição em (3.16) ocupada pelo vértice 2) e o que lhesegue será o 1 (no lugar do 3 em (3.16)). Obtemos:

gradφ(k)j (x, y) =

1

2Ak

(y

(k)3 − y(k)

1 , x(k)1 − x

(k)3

). (3.17)

Para uma função associada ao vértice 3:

gradφ(k)j (x, y) =

1

2Ak

(y

(k)1 − y(k)

2 , x(k)2 − x

(k)1

). (3.18)

Obviamente, se φj não é associada a nenhum vértice do triângulo Tk, então gradφ(k)j (x, y) = 0.

As expressões (3.16), (3.17) e (3.18) dão o valor do gradiente de uma função da base em um elementose ela lhe for associada ao vértice 1, 2 ou 3, respectivamente. Poderá a leitora se perguntar: digamosque φj é associada ao vértice 1 do elemento Tk mas também é associada ao vértice 3 de outro elemento,

Tl. Neste caso, quanto vale gradφ(l)j ? Ora, como essa função é associada ao vértice 3 de Tl, usamos

(3.18) com as coordenadas de Tl: gradφ(l)j (x, y) =

1

2Al

(y

(l)1 − y

(l)2 , x

(l)2 − x

(l)1

).

Como ao discretizar a malha conhecemos as coordenadas dos vértices dos elementos, podemoscalcular os gradientes de todas as N funções de B em todos os elementos da malha. A partir daí, épossível calcular os elementos Mij da matriz de rigidez e, usando (3.15), o vetor de carga. Escrevemosassim o sistema (3.9).

Antes de darmos o assunto por encerrado, veremos alguns detalhes do cálculo deM, cujas entradassão Mij =

∫Ωd

〈gradφi , gradφj〉 dV . Reparemos que a integral sobre Ωd se decompõe em uma soma de

integrais, cada uma sobre um elemento. Isto é:

Mij =

∫Ωd

〈gradφi , gradφj〉 dV =

m∑k=1

∫Tk

〈gradφ(k)i , gradφ(k)

j 〉 dV.

Sabemos por (3.16)-(3.18) que o gradiente de uma função da base é constante em cada elemento.Por isso, podemos passar o produto dos gradientes para fora das integrais, obtendo

Mij =m∑k=1

(〈gradφ(k)

i , gradφ(k)j 〉

∫Tk

dV

)=

m∑k=1

〈gradφ(k)i , gradφ(k)

j 〉Ak. (3.19)

Da forma como denimos as funções de B, elas só têm valores e gradientes não-nulos nos seuselementos associados. Como conseqüência disso, o produto 〈gradφ(k)

i , gradφ(k)j 〉 só será diferente de

zero se ambas funções φi e φj forem associadas à Tk. Isso faz com que a maior parte dos termos dosomatório (3.19) sejam nulos. Mais ainda, um grande número de entradas de M são zeros: a matrizde rigidez é esparsa.

Concluímos que para calcular os Mij basta considerar os elementos associados ao mesmo tempo aφi e a φj .

16

Da mesma forma como transformamos uma integral sobre Ωd em uma soma de integrais sobre oselementos para calcular os termos da matriz de rigidez, podemos fazê-lo também para o vetor de carga.Seus elementos são do tipo ∫

Ωd

φi · f dV =m∑k=1

∫Tk

φ(k)i · f dV.

Novamente, apenas os termos do somatório com Tk associado a φi serão não-nulos.O que queríamos mostrar é como transformamos a integral sobre Ωd em uma soma de integrais

sobre os elementos relevantes no cálculo. De certo modo nossa exposição do problema de Dirichletbidimensional está terminada. Apenas a título de ilustração dos procedimentos, calcularemos a matrizde rigidez em dois exemplos de malhas.

3.4 Alguns exemplos de malhas

Nesta seção consideraremos alguns exemplos de triangulações para mostrar como é feita a construçãodo sistema (3.9). As malhas que mostraremos são bastante simples, com poucos elementos, já quedesejamos apenas ilustrar o método. Em aplicações práticas um número bem superior de elementosdeve ser utilizado. Nossas malhas podem ser consideradas células de malhas maiores. Se for mantidasua regularidade, os resultados aqui obtidos podem ser muito facilmente adaptados para aquelas.

3.4.1 Exemplo 1

Considere Ω = [0, L] × [0, L], o quadrado de lado L. Dividimos nosso domínio em nove quadrados delado L/3 e, em seguida, traçamos uma diagonal em cada quadrado, à maneira da Figura 8a.

Figura 8: Malha e enumeração global dos elementos (a) e dos vértices (b). (c) mostra a enumeração localdos vértices dos elementos que serão utilizados neste exemplo.

Nossa malha tem 18 elementos e 16 vértices, sendo que apenas 4 são interiores. Enumeraremos osvértices globalmente como mostra a Figura 8b, e localmente consoante a Figura 8c. As áreas de todoselementos são iguais, e representaremos simplesmente por A.

Notemos que cada função da base (associadas aos vértices 1, 2, 3 e 4) é associada a seis elementos.Mais ainda, esse conjunto função-base + os seis elementos associados forma uma espécie de célula,sendo transladado equivale aos outros conjuntos semelhantes. Isso é conseqüência da regularidadedesta malha.

A dependência de φ1 com os seis elementos vizinhos ao vértice p1 é a mesma de φi qualquer comseus seis elementos associados. Assim, calculando o gradiente de φ1 em T1, T2, T3, T10, T9 e T8, játeremos os gradientes das outras funções.

Começemos, pois, pelo elemento T1. φ1 é associada ao seu vértice 1, logo, por (3.16),

gradφ(1)1 (x, y) =

1

2A

(y

(1)2 − y

(1)3 , x

(1)3 − x

(1)2

)=

1

2A(L/3− 0 , 0− L/3) =

L

6A(1,−1).

17

Em T2, φ1 também é associada ao vértice (local) 1. Então, usando novamente (3.16),

gradφ(2)1 (x, y) =

1

2A

(y

(2)2 − y

(2)3 , x

(2)3 − x

(2)2

)=

1

2A(0− 0 , 2L/3− L/3) =

L

6A(0, 1).

Em T3, φ1 é associada ao vértice (local) 2. Usamos, pois, (3.17):

gradφ(3)1 (x, y) =

1

2A

(y

(3)3 − y

(3)1 , x

(3)1 − x

(3)3

)=

1

2A(0− L/3 , 2L/3− L/3) =

L

6A(−1, 1).

Para calcular o gradiente de φ1 em T10 não faremos conta alguma: descobriremos seu valor pelageometria da malha. Notemos que o lado 1, 2 do elemento T10 é paralelo ao lado 2, 3 de T1. Assim,como o gráco de φ1 é uma pirâmide, é fácil ver que a direção de crescimento de φ1 é a mesma nessesdois elementos. Como o gradiente tem justamente essa direção, nesses elementos eles são paralelos. Jáque T10 é simétrico com relação à p1 ao elemento T1, seus gradientes têm mesmo módulo e sentidosopostos. Logo,

gradφ(10)1 = −gradφ(1)

1 =L

6A(−1, 1).

Da mesma forma, T9 é simétrico com T2 e T8 o é com T3. Portanto:

gradφ(9)1 = −gradφ(2)

1 =L

6A(0,−1),

gradφ(8)1 = −gradφ(3)

1 =L

6A(1,−1).

Já temos, então, os gradientes de φ1. Podemos calcular o primeiro termo da matriz de rigidez:

M11 =

∫Ω

〈gradφ1 , gradφ1〉 dV =

∫T1

〈 L6A

(1,−1) ,L

6A(1,−1)〉 dV +

∫T2

〈 L6A

(0, 1) ,L

6A(0, 1)〉 dV+

+

∫T3

〈 L6A

(−1, 1) ,L

6A(−1, 1)〉 dV +

∫T10

〈 L6A

(−1, 1) ,L

6A(−1, 1)〉 dV +

∫T9

〈 L6A

(0,−1) ,L

6A(0,−1)〉 dV+

+

∫T8

〈 L6A

(1,−1) ,L

6A(1,−1)〉 dV =

L2

36A(2 + 1 + 2 + 2 + 1 + 2) =

5L2

18A.

Como a área total do domínio é L2 e ele está dividido em 18 triângulos de mesma área A, temosque A = L2/18. Substiuindo isso no resultado acima, M11 = 5.

Por simples inspeção da malha, e nos baseando nas considerações já feitas sobre a simetria datriangulação utilizada, vemos que

gradφ(1)1 = gradφ(3)

2 = gradφ(7)3 = gradφ(9)

4 ,

gradφ(2)1 = gradφ(4)

2 = gradφ(8)3 = gradφ(10)

4 ,

gradφ(3)1 = gradφ(5)

2 = gradφ(9)3 = gradφ(11)

4 ,

gradφ(10)1 = gradφ(12)

2 = gradφ(16)3 = gradφ(18)

4 ,

gradφ(9)1 = gradφ(11)

2 = gradφ(15)3 = gradφ(17)

4 ,

gradφ(8)1 = gradφ(10)

2 = gradφ(14)3 = gradφ(16)

4 .

Já conhecemos então todos os gradientes. Nos elementos que não estão relacionados na lista acima, osgradientes são nulos.

Daí então, M11 = M22 = M33 = M44 = 5 .Por inspeção da malha, vemos que as funções

(a) φ1 e φ2 têm os elementos associados T3 e T10 em comum;(b) φ1 e φ3 têm os elementos associados T8 e T9 em comum;

18

(c) φ1 e φ4 têm os elementos associados T9 e T10 em comum;

(d) φ2 e φ3 não têm elementos associados em comum;(e) φ2 e φ4 têm os elementos associados T10 e T11 em comum;

(f) φ3 e φ4 têm os elementos associados T9 e T16 em comum.Notemos, ainda com base na malha, que as relações (e),(b) e (f),(a) têm mesma geometria. Isso

mostra que basta analisar uma célula da malha e como esta se relaciona com suas vizinhas para entendero comportamento de toda a malha - claro, no caso de uma triangulação regular.

Feitas essas considerações, podemos calcular os demais elementos de M usando (3.19):

(a) ⇒ M12 = M21 = A〈gradφ(3)1 , gradφ(3)

2 〉 + A〈gradφ(10)1 , gradφ(10)

2 〉 = A〈gradφ(3)1 , gradφ(1)

1 〉 +

A〈gradφ(10)1 , gradφ(8)

1 〉 = 2A(−4)L2

36A2= −2L2

9A= −4.

(b) ⇒ M13 = M31 = A〈gradφ(8)1 , gradφ(8)

3 〉 + A〈gradφ(9)1 , gradφ(9)

3 〉 = A〈gradφ(8)1 , gradφ(2)

1 〉 +

A〈gradφ(9)1 , gradφ(3)

1 〉 = 2A(−2)L2

36A2= −L

2

9A= −2.

(c) ⇒ M14 = M41 = A〈gradφ(9)1 , gradφ(9)

4 〉 + A〈gradφ(10)1 , gradφ(10)

4 〉 = A〈gradφ(9)1 , gradφ(1)

1 〉 +

A〈gradφ(10)1 , gradφ(2)

1 〉 = 2A(2)L2

36A2=L2

9A= 2.

(d) ⇒ M23 = M32 = 0.

(e), (b) ⇒ M24 = M42 = −2.

(f), (a) ⇒ M34 = M43 = −4.

Logo,

M =

5 −4 −2 2−4 5 0 −2−2 0 5 −42 −2 −4 5

.Vemos queM acima tem poucos zeros, em contradição com o que há pouco armamos, que a matriz

de rigidez é esparsa. Essa aparente incoerência ocorre devido ao tamanho da malha considerada. Osúnicos zeros de esparsidade que ocorrem são devidos às funções φ2 e φ3 (que não têm elementosassociados em comum). No entanto, se considerássemos uma malha formada com o mesmo padrão,porém com nove pontos interiores, o número de funções sem elementos associados em comum aumentarábastante. Um pouco de reexão bastará para que o leitor se convença de queM para a malha mostradana Figura 9 será dada por

M =

5 −4 0 −2 2 0 0 0 0−4 5 −4 0 −2 2 0 0 00 −4 5 0 0 −2 0 0 0−2 0 0 5 −4 0 −2 2 02 −2 0 −4 5 −4 0 −2 20 2 −2 0 −4 5 0 0 −20 0 0 −2 0 0 5 −4 00 0 0 2 −2 0 −4 5 −40 0 0 0 2 −2 0 −4 5

,

que tem aproximadamente metade dos seus elementos nulos. À medida que o número de vérticesinteriores da malha aumentar, essa proporção também o fará.

19

Figura 9: Malha de 32 elementos e 9 vértices interiores (enumerados).

Incentivamos o interessado a sempre buscar compreender a simetria de uma malha regular, comozemos neste exemplo. Isso torna fácil a tarefa de escrever a matriz de rigidez para malhas maioresque seguem o mesmo padrão.

Por m, chamamos atenção para o fato de que uma mudança no sistema de enumeração dos vérticesinteriores causa alteração na matriz de rigidez (pois ocorre uma reordenação da base). Por exemplo, seao escrever a matriz para o caso de quatro vértices interiores tivéssemos usado a numeração da Figura10 ao invés da Figura 8, chegaríamos na matriz M 6= M (verique):

M =

5 −4 2 −2−4 5 −2 02 −2 5 −4−2 0 −4 5

.

Figura 10: Malha de 18 elementos com outra enumeração dos vértices interiores.

3.4.2 Exemplo 2

Considere Ω = [0, L]× [0, L], o quadrado de lado L. Dividimos nosso domínio em quatro quadrados delado L/2 e, em seguida, traçamos as duas diagonais em cada quadrado, à maneira da Figura 11a. Amalha formada tem 16 elementos e 13 vértices, sendo que 5 são interiores. Enumeramo-los globalmenteconsoante a Figura 11b. Na Figura 11c mostramos enumerações locais de vértices em alguns triângulos.Veremos que só precisaremos desses elementos para escrever a matriz de rigidez.

20

Figura 11: Enumeração dos elementos (a), enumeração global dos vértices (b) e enumeração local dosvértices de alguns elementos (c).

A grande diferença dessa malha para a do Exemplo 1 é que, enquanto lá cada função era associadaa seis elementos, aqui existem funções associadas a quatro e a oito elementos. Por inspeção da malha,vemos que as funções que são associadas a 4 elementos estão, por assim dizer, encerradas; que elasnão têm nenhum elemento associado em comum. A única função que tem elementos em comum comoutras é a associada a oito elementos, φ3. A partir dessa análise simples, já podemos garantir queM12 = M21 = M14 = M41 = M15 = M51 = M24 = M42 = M45 = M54 = 0. Para os outrostermos teremos que fazer contas, mas, na medida do possível, utilizaremos da simetria da malha parasimplicá-las.

Consideremos primeiro φ1. Ela é não-nula apenas em T1, T2, T3 e T4. Pela geometria da malhavemos que T1 e T3 têm gradientes de mesmo módulo e sentidos opostos, da mesma forma que T2 e T4.Então, como T1 e T2 têm p1 como vértice local 1, (3.16) implica em:

gradφ(1)1 =

1

2A

(y

(1)2 − y

(1)3 , x

(1)3 − x

(1)2

)=

1

2A(L/2− 0 , 0− 0) =

L

4A(1, 0)

e

gradφ(2)1 =

1

2A

(y

(2)2 − y

(2)3 , x

(2)3 − x

(2)2

)=

1

2A(0− 0 , L/2− 0) =

L

4A(0, 1).

Ainda por inspeção da malha, vemos que essa estrutura de função da base com quatro elementosao redor se repete pela malha, por uma simples translação. Portanto:

gradφ(1)1 = gradφ(5)

2 = gradφ(9)4 = gradφ(13)

5 = −gradφ(3)1 = −gradφ(7)

2 = −gradφ(11)4 = −gradφ(15)

5 ;

gradφ(2)1 = gradφ(6)

2 = gradφ(10)4 = gradφ(14)

5 = −gradφ(4)1 = −gradφ(8)

2 = −gradφ(12)4 = −gradφ(16)

5 ,

e, ainda, M11 = M22 = M44 = M55. Logo,

Mii(i 6=3) = 4L2

16A2= 4,

21

pois a área de cada elemento é A = L2/16.Resta, agora, calcular o gradiente de φ3, aquela que é associada a 8 elementos. Novamente usaremos

o argumento de que os vértices simétricos com relação p3 têm gradientes de mesmo módulo e sentidosopostos.

Podemos numerar, localmente, os vértices de T3, T4, T5 e T8 de modo que o vértice 1 seja semprep3. Isso nos faz usar apenas a fórmula (3.16) para cálculo dos gradientes. Temos, pois,

gradφ(4)3 =

1

2A

(y

(4)2 − y

(4)3 , x

(4)3 − x

(4)2

)=

1

2A(L/2− L/4 , L/4− 0) =

L

8A(1, 1) = −gradφ(14)

3

gradφ(3)3 =

1

2A

(y

(3)2 − y

(3)3 , x

(3)3 − x

(3)2

)=

1

2A(L/4− 0 , L/2− L/4) =

L

8A(1, 1) = −gradφ(13)

3

Repare que encontramos gradφ(4)3 = gradφ(3)

3 . Observando a malha, já poderíamos esperar isso, pois ográco de φ3 é uma pirâmide de base quadrada e T3 e T4 formam um mesmo lado desse quadrado.

Usaremos esse argumento para armar que

gradφ(5)3 = gradφ(8)

3 =1

2A

(y

(8)2 − y

(8)3 , x

(8)3 − x

(8)2

)=

1

2A(L/4− L/2 , L− 3L/4) =

L

8A(−1, 1) =

= −gradφ(10)3 = −gradφ(11)

3 .Pronto: já conhecemos os gradientes das funções φi em todos os elementos da malha. Calculemos,

pois, o restante das entradas de M. Usando (3.19):

M13 = M31 = 〈gradφ(3)1 , gradφ(3)

3 〉A+ 〈gradφ(4)1 , gradφ(4)

3 〉A = − L2

32A2= −2.

Usando as considerações feitas, é possível mostrar que M23 = M32 = M34 = M43 = M35 = M53 =M13 = −2.

Só resta, pois, M33 = 8L2

64A22 = 4. Portanto,

M =

4 0 −2 0 00 4 −2 0 0−2 −2 4 −2 −20 0 −2 4 00 0 −2 0 4

. (3.20)

Um pouco de reexão mostrará que para a malha da Figura 12 (verique):

M =

4 0 0 −2 0 0 0 0 0 0 0 0 00 4 0 −2 −2 0 0 0 0 0 0 0 00 0 4 0 −2 0 0 0 0 0 0 0 0−2 −2 0 4 0 −2 −2 0 0 0 0 0 0

0 −2 −2 0 4 0 −2 −2 0 0 0 0 00 0 0 −2 0 4 0 0 −2 0 0 0 00 0 0 −2 −2 0 4 0 −2 −2 0 0 00 0 0 0 −2 0 0 4 0 −2 0 0 00 0 0 0 0 −2 −2 0 4 0 −2 −2 00 0 0 0 0 0 −2 −2 0 4 0 −2 −20 0 0 0 0 0 0 0 −2 0 4 0 00 0 0 0 0 0 0 0 −2 −2 0 4 00 0 0 0 0 0 0 0 0 −2 0 0 4

.

Chamamos atenção para o fato de que a matriz de rigidez depende da numeração global dos vérticesinteriores, mas independe da enumeração local. Em casos de malhas simétricas podemos escolher estaúltima de modo a usar apenas uma fórmula para o gradiente, como zemos neste exemplo.

22

Figura 12: Malha com 36 elementos e 13 vértices interiores (enumerados).

3.5 Outros casos

Como dissemos logo no início deste texto, nosso objetivo é apenas transmitir a essência do método doselementos nitos, por isso este material é bastante simples. Mencionamos aqui, brevemente, algumasoutras possibilidades que o Método permite.

Contemplamos apenas os casos de uma e duas dimensões. A formulação do caso tridimensionalpode ser deduzida sem grandes diculdades a partir da dedução feita neste capítulo. Foi, inclusive,com esse intuito que deixamos o Teorema do Divergente enunciado em sua forma geral. Novamente, adiculdade irá surgir ao discretizar o domínio (agora em tetraedros) e buscar escrever uma base parao espaço de funções Vd.

Ao longo deste texto, sempre aproximamos a função u por outra que tinha a propriedade deser linear em cada elemento da malha. Existem outras possibilidades: podemos desejar que ud sejaquadrática por partes, ou mesmo polinomial por partes, fornecendo aproximações mais suaves.

Por m, o Método não se aplica apenas ao problema de Dirichlet. Condições de contorno deNeumann e Robin também são aceitas, com algumas alterações no exposto neste texto. Por exemplo,no caso de condição de Neumann, não podemos, nas regiões de ∂Ω com essa condição, igualar o membroà direita de (3.5), fórmula de Green, a zero.

23

Bibliograa

[1] Jochen ALBERTY, Carsten CARSTENSEN, Stefan A. FUNKEN, Remarks around 50 lines of

Matlab: short nite element implementation. Numerical Algorithms 20 (1999), 117-137.

[2] Rodney Josué BIEZUNER, Notas de aula: Autovalores do Laplaciano. UFMG, 2006.

[3] Giovanni CALDERÓN, Rodolfo GALLO, Introducción al Método de los Elementos Finitos: un

enfoque matemático. Caracas: Ediciones IVIC, 2011.

[4] Jichun LI, Yi-Tung CHEN, Computational partial dierential equations using MATLAB. BocaRaton: CRC Press, 2009.

[5] James STEWART, Cálculo: volume 2. São Paulo: Cengage Learning, 2009.

24