Departamento de Matemtica | Instituto de Cincias Exatas | Universidade Federal de Minas Gerais

Uma breve introduo ao

Mtodo dos Elementos Finitos

Breno Loureiro Giacchini

Janeiro de 2012

Contedo

Prefcio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1 Introduo 2

2 O problema unidimensional 32.1 Formulao fraca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Discretizao do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Existncia e unicidade da soluo do problema aproximado . . . . . . . . . . . . . . . 62.4 Um caso particular: partio regular do intervalo . . . . . . . . . . . . . . . . . . . . . 7

3 O problema bidimensional 83.1 Formulao fraca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2 Problema aproximado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.3 Uma base para Vd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.3.1 Enumeraes dos vrtices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.3.2 Funes da base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.3.3 Clculo 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

Prefcio

Em meio a tantos bons livros e apostilas sobre o Mtodo dos Elementos Finitos, o questionamento doporqu da escrita deste texto no de todo descabido. O que nos motivou a escrev-lo a diculdadede se encontrar um texto, em portugus, que apresente o Mtodo de forma simples e direta, que lhefornea uma idia geral e ao mesmo tempo permita sua implementao em casos simples, mas semgrandes delongas em formalismos matemticos e prembulos sobre anlise funcional, por exemplo.

Admitimos, pois, que nossa exposio do Mtodo no feita com todo rigor matemtico nem emtoda sua generalidade, mas cremos que essa opo satisfaz ao estudante que deseja entender sua essnciae aplic-lo, de forma rpida, em alguns casos; ou ao interessado em ter um primeiro contato com essatcnica de resolver numericamente problemas de valores de contorno. Por ser um texto introdutrio,que d apenas um sabor do Mtodo, nos limitamos contemplar o problema de Dirichlet homogneouni e bidimensional e a utilizar, neste caso, apenas elementos triangulares.

Deixamos expresso nosso agradecimento ao apoio da Fundao de Amparo Pesquisa do Estadode Minas Gerais FAPEMIG , que nanciou o projeto de pesquisa Obteno dos autovalores eautofunes do laplaciano via o quociente de Rayleigh, do qual o estudo do Mtodo dos ElementosFinitos e a escrita deste texto foram partes integrantes. Tambm agradecemos ao professor RodneyJosu Biezuner, nosso orientador neste projeto.

1

Captulo 1

Introduo

Diversos problemas da Fsica, Engenharia e outras cincias aparecem sob a forma de uma equaode Poisson

u = f (x) em (1.1)

com condio de fronteira de Dirichlet u = c sobre , sendo c uma funo 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 condio de Dirichlet homognea. Ao conjunto de umaequao de Poisson com uma condio de Dirichlet homognea chamamos um problema de Dirichlethomogneo: {

u = f(x) em ,u = 0 sobre .

(1.2)

Dependendo da geometria do domnio a soluo do problema pode ser obtida analiticamentena forma de sries de Fourier. Exemplos clssicos normalmente estudados num curso de equaesdiferenciais parciais so o caso de retngulos, semiplanos, discos e paraleleppedos. No entanto, preciso recorrer a mtodos numricos caso o domnio se torne mais elaborado. O mtodo dos elementosnitos (MEF) conhecido por ser robusto e aplicvel em domnios deveras elaborados. Essas tambmso algumas de suas vantagens sobre o mtodo das Diferenas Finitas, tambm bastante popular.

A idia central do MEF discretizar o domnio, representando-o, ainda que de forma aproximada,por uma reunio de um nmero nito de elementos; e resolver no o problema original (1.2), mas sim umque lhe associado sua forma fraca. No caso de um domnio plano, os elementos podem ser tringulosou quadrilteros. O mtodo pode ser utilizado para resolver no s problemas elpticos, como o hpouco mencionado; e as condies no necessitam ser de Dirichlet: o MEF tambm aplicvel no casode condio 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 introduo ao mtodo. 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 estaro aptosa utilizar dessa ferramenta na resoluo de alguns problemas de interesse.

1Se Rn e u funo u : R, u =n

i=1

d2u

dxi2.

2

Captulo 2

O problema unidimensional

2.1 Formulao fraca

O problema de Dirichlet homogneo unidimensional se escreved

2u

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

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

(2.1)

Assumiremos que a funo f : [0, 1] R limitada e contnua por partes. Isso necessrio porqueo MEF supe que f integrvel. Notamos que aqui = [0, 1].

Ao invs de resolver o problema (2.1) da forma como est escrito, o MEF se prope a solucionarum problema equivalente, chamado formulao fraca do original. Para escrever (2.1) dessa forma,

principiamos denindo o espao de funes V = {v; v funo contnua em [0, 1], dvdx

limitada e

contnua por partes, e v(0) = v(1) = 0}.Em seguida, multiplicamos a primeira equao de (2.1) por uma funo qualquer de V e integramos

a equao resultante em :

1

0

d2u

dx2 v dx =

10

f(x) v dx .

Integrando por partes e lembrando que v sastifaz a condio de Dirichlet homognea:

(vdu

dx

) 1

0

+

10

du

dx

dv

dxdx =

10

f v dx

1

0

du

dx

dv

dxdx =

10

f v dx (2.2)

para todo v V . A equao (2.2) juntamente com a condio de Dirichlet homognea a formulaofraca do problema (2.1).

Mostraremos agora que a existncia de uma soluo de (2.1), problema original, implica na equi-valncia entre os problemas de formulao forte e fraca. Vimos logo acima que formulao forte formuo fraca. Resta, pois, vericar a recproca.

J que supomos que o problema original tem soluo, sabemos qued2u

dx2existe e contnua por

partes. Podemos, ento, integrar por partes a formulao fraca (2.2). O que obtemos justamente aformulao forte (original):

3

10

du

dx

dv

dxdx =

(vdu

dx

) 1

0

1

0

vd2u

dx2=

10

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 funo v em V , o termo do integrando que est entreparntesis deve ser nulo; e assim chegamos ao problema original:

f(x) +d2u

dx2= 0 , 0 < x < 1.

Com isso provamos que uma funo que resolve o problema forte tambm soluo do fraco; e que,se a soluo do problema fraco for sucientemente regular, ela tambm resolver o problema forte. Nomtodo dos elementos nitos resolveremos o problema fraco, (2.2).

2.2 Discretizao do problema

Na formulao original, o problema de Dirichlet contnuo e seu espao de solues pode ter dimensoinnita1. Aproximaremos o problema contnuo por outro discreto, cuja soluo est em um espao dedimenso nita. Isso feito dividindo o domnio (o intervalo [0, 1]) em um nmero nito de subintervalosIj = [xj1, xj ], 1 6 j 6 N + 1, N N, com 0 = x0 < x1 < < xN < xN+1 = 1. Cada subintervalotem comprimento hj = xj xj1.

Essa discretizao uma partio [do intervalo], cuja norma denimos h = maxj{hj}, o compri-mento do maior dos subintervalos.

O termo discretizao usado justamente porque passamos de um contnuo (a funo originalest denida num domnio que uma reunio no-enumervel de pontos) para um conjunto discreto: odomnio passa a ser uma reunio nita de intervalos2. Em cada um desses intervalos Ij , aproximamosa funo original u por um segmento de reta de extremos u(xj1) e u(xj) (Figura 1a). Evidentemente,quanto menor o comprimento dos subintervalos, ou seja, quanto menor a norma da partio, maisa funo discretizada ud se aproximar da original u (Figura 1a-b). Notemos, ainda, que ud comodenida contnua.

Figura 1: Aproximao de uma funo suave por outra linear por partes. Quanto menor a norma dapartio, melhor a aproximao.

1Exemplos de problemas cujas solues esto num espao de dimenso innita so aqueles tradicionalmente estudadosnum curso de equaes diferenciais parciais que tm como resultado sries de Fourier innitas. Cada uma daquelas funes

sen(nxL

)ou cos

(nxL

) uma funo da base do espao que contm a soluo do problema. Como existem innitos

n N, a base um conjunto innito.2Uma reunio nita de intervalos ainda um conjunto no-enumervel de pontos. A idia aqui que, ao invs de

buscarmos u com denio ponto a ponto, vamos aproxim-la por uma funo que denida subintervalo por subintervaloe, nesse sentido, ela ser denida discretamente. Mais adiante no texto car clara essa idia.

4

Para discretizar o problema na forma fraca, devemos tambm aproximar o espao V por um dedimenso nita, Vd = {v; v contnua em [0, 1], v linear em cada Ij e v(0) = v(1) = 0}. Notemosque Vd V de sorte que ao tomarmos uma funo v Vd no ferimos a condio v V da formulaofraca.

Nosso problema discretizado (ou aproximado) , ento, encontrar ud Vd tal que

10

duddx

dv

dxdx =

10

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

Esse ud ser a aproximao para a funo u desejada.Observao 1: note que a condio de fronteira u(0) = u(1) = 0 est contida no enunciado do pro-blema discretizado (2.3) j que ud Vd implica na condio de Dirichlet homognea.Observao 2: talvez tenha parecido ao leitor que nos precipitamos ao declarar que a soluo apro-ximada ud que buscamos est em Vd. Veremos porque podemos assumir isso. A primeira condioque ud deve satisfazer para pertencer a Vd ser contnua. J vimos que como foi denida ud u contnua. A segunda condio (linearidade em cada subintervalo) vem tambm da discretizao doproblema e est relacionada com a qualidade dessa aproximao - aproximamos uma curva suave poroutra poligonal. Ao assumirmos que buscamos uma soluo com essa aproximao, ud satisfaz, porconseguinte, a segunda condio. Por m, a terceira justamente a condio de Dirichlet que tantou quando ud satisfazem por hiptese. Conclumos assim que ud Vd.

Ao discretizarmos o espao V , o aproximamos por um de dimenso nita. Como v Vd linearem cada Ij , as funes

j(x) =

x xj1

hjse x [xj1, xj ],

xj+1 xhj+1

se x [xj , xj+1],

0 caso contrrio.

(2.4)

formam uma base, B, de Vd.

Figura 2: Grco da funo j de B.

Mostramos na Figura 2 um grco de uma dessas funes-base de Vd. fcil ver que se i 6= j, asfunes i e j so linearmente independentes. Um momento de reexo bastar para que a leitora seconvena de que qualquer funo de Vd se escreve em termos das j acima denidas.

Ora, como ud pertence a Vd, ser da forma

ud(x) =

Nj=1

j j(x) , x [0, 1], (2.5)

e o nosso problema (2.3) se escrever

10

d

dx

(j

jj

)dv

dxdx =

10

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

5

Recordemos que a funo v uma funo qualquer de Vd. Escolhemos, ento, v como sendo umadas funes da base: v = i para algum i 6 N . Para esse i, (2.6) implica em

10

(Nj=1

jdjdx

)didx

dx =

10

f(x)i(x) dx

Nj=1

(j

10

djdx

didx

dx

)=

10

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

Variando i de 1 a N, (2.7) resulta em um sistema de N equaes e N incgnitas j :

10

d1dx

d1dx

dx10

d2dx

d1dx

dx 10

dNdx

d1dx

dx

10

d1dx

d2dx

dx10

d2dx

d2dx

dx 10

dNdx

d2dx

dx

......

. . ....

10

d1dx

dNdx

dx10

d2dx

dNdx

dx 10

dNdx

dNdx

dx

12...N

=

10

f(x)1(x) dx

10

f(x)2(x) dx

...10

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 funo ud Vd que satisfaz (2.3) se reduz resoluode um sistema linear. Resolvendo-o, determinamos os coecientes j e podemos construir, usando (2.5)e (2.4), a funo ud u. Antes de darmos o assunto por encerrado, vericaremos algumas propriedadesda matriz M e provaremos que sempre existir uma (nica) soluo para o sistema (2.8) - e, portanto,para o problema aproximado (2.3).

2.3 Existncia e unicidade da soluo do problema aproximado

Proposio 1. A matriz M goza das seguintes propriedades:R1) simtrica;

R2) tridiagonal;

R3) positiva denida - isto , wTMw > 0 w no-nulo em RN .

Demonstrao.R1) conseqncia da comutatividade do produto de funes:

mij =

10

djdx

didx

dx =

10

didx

djdx

dx = mji .

R2) Calcularemos os elementos mij para mostrar que M tridiagonal. Para tanto, usaremos asderivadas das funes j denidas por (2.4).

Se i = j, mii =

xixi1

1

h2idx+

xi+1xi

1

h2i+1dx =

1

hi+

1

hi+1.

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

xixi1

1hi

1

hidx = xi xi1

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,i1 so no-nulos.

6

R3) Ora, wTMw =Ni=1

Nj=1

wj

( 10

djdx

didx

dx

)wi =

10

[(Nj=1

wjdjdx

)(Ni=1

wididx

)]dx =

=

10

(Nj=1

wjdjdx

)2dx > 0.

Como essa relao vale para qualquer vetor w em RN {0}, a igualdade s se vericaria casodjdx

= 0 para cada j e em todo o intervalo [0, 1]. Como tal situao no 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 no-nulo3. Outro teorema reza que que se a matriz de um sistema linear tem determinanteno-nulo, o sistema tem soluo nica. Esses teoremas, juntamente com o terceiro item da Proposio1, nos asseguram que (2.8) tem soluo - e ela nica.

2.4 Um caso particular: partio regular do intervalo

Concluimos o estudo do caso unidimensional escrevendo o sistema (2.8) num caso particular de partiodo intervalo, a saber, considerando que todos os subintervalos Ij tm mesmo comprimento h. A umapartio deste tipo d-se o nome de regular.

Utilizando os clculos realizados na prova do segundo item da Proposio 1, temos que uma partioregular do domnio fornece a matriz de rigidez:

M =1

h

2 11 2 1

1 . . . . . .. . . . . . 1

1 2 11 2

E o sistema (2.8) pode ser escrito como:

2 11 2 1

1 . . . . . .. . . . . . 1

1 2 11 2

12

...

N

= h

10

f(x)1(x) dx

10

f(x)2(x) dx

...10

f(x)N (x) dx

(2.9)

O leitor que j estudou o Mtodo das Diferenas Finitas notar que a matriz M obtida para umapartio regular bastante semelhante encontrada naquele mtodo - elas s diferem por um termo 1/hmultiplicando. Dependendo da maneira de discretizar as integrais do vetor de carga, os dois mtodoscoincidiro.

3Mais precisamente, seu determinante maior que zero.

7

Captulo 3

O problema bidimensional

3.1 Formulao fraca

Sejam R2 um aberto limitado e f uma funo real contnua por partes e limitada, em . Oproblema de Dirichlet homogneo 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 espao de funes V = {v : R2 R; v funo contnua em , vx

ev

yso contnuas

por partes em , e v = 0 sobre }. Multiplicando a equao de Poisson do problema (3.1) por umafuno qualquer de V e integrando sobre temos:

v u = v f

v u dV =

v f dV. (3.2)

Podemos reescrever a equao acima de forma mais conveniente usando a frmula de Green, quese baseia noTeorema 1. [Teorema do divergente] Seja Rn compacto e com fronteira suave por partes.Se w um campo de vetores diferencivel denido em , ento:

div w dV =

w , n ds,

onde n representa o vetor unitrio normal .A notao1 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 Clculo, pelo menos para n = 2 e n = 3. O caso n = 2 que nos interessa s vezes chamadoTeorema de Green. Para obter a frmula 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 funes g,h : R2 R. Considerando

que o vetor normal unitrio n = (n1 , n2), temos, para a,

(g

2h

x2+g

x

h

x

)dV =

g hx n1 ds (3.3)

1Outras notaes e denies 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 funo f : Rn R: gradf =(f

x1, , f

xn

).

8

e, para b:

(g

2h

y2+g

y

h

y

)dV =

g hy n2 ds. (3.4)

Somando membro a membro as equaes (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. (Frmula de Green) (3.5)

Notemos que se a funo g(x, y) acima satisfaz a condio de Dirichlet homognea, a integral sobre em (3.5) nula e a frmula de Green implica em

g h dV =

gradg , gradh dV.

Comparando (3.2) com a equao acima, vemos que os membros esquerdos so iguais se zermosg = v e h = u. Temos, portanto, que

v f dV =

gradv , gradu dV v V. (3.6)

Esta equao acrescida da condio de Dirichlet homognea formam a formulao fraca do problemabidimensional (3.1). possvel mostrar, como o zemos no caso unidimensional, que as formas fracae forte so equivalentes e que uma soluo da forma fraca, se sucientemente regular, tambm sersoluo da forma forte.

3.2 Problema aproximado

Uma vez compreendida a essncia do mtodo no caso unidimensional, o caso dos domnios planosno apresenta maiores diculdades no que tange essa essncia. A diculdade surge no momento dediscretizar e trabalhar com a malha resultante, como veremos em sees seguintes.

Comearemos a discretizao do problema dividindo o domnio em tringulos. Obviamente, no qualquer domnio que aceita essa diviso perto de sua borda. Neste caso, aproximamos por d cujafronteira uma curva poligonal (formada por unies nitas de segmentos de retas). A cada um dessestringulos chamamos elemento. A discretizao em tringulos (ou triangulao) deve cumprir asseguintes condies:D1) A reunio de todos os elementos forma d, que aproxima ;D2) Os elementos no se sobrepem;

D3) Os vrtices de um elemento nunca ocorrem no lado de outro elemento.

A Figura 3 mostra exemplos de triangulaes permitidas e no permitidas no mtodo dos elementosnitos, alm de ilustrar como podemos fazer a aproximao da fronteira.

Figura 3: (a) exemplo de triangulao permitida. A partio em (b) no permitida pois o trao em azuldene um vrtice que ocorre em um lado de outro elemento.

9

No problema discretizado, buscamos uma funo ud que aproxima u. Aproveitamos nossa malha(domnio discretizado) para impor uma condio sobre ud: que ela seja contnua e linear em cadaelemento. Esta ltima imposio signica que o grco de ud em cada elemento um pedao de planocontido no R3. Ao fazermos essas exigncias estamos apenas escolhendo como iremos aproximar u.

Chamamos ateno para o fato de que ao aproximar por d nos propomos a resolver o problema(3.1) com o domnio d. Portanto, quanto mais parecido for d de , mais a funo encontrada udser parecida com a funo real u.

Para discretizar o problema na forma fraca, devemos ainda aproximar o espao V por um nitoVd = {v; v contnua em d, v linear em cada elemento e v = 0 sobre d}. Como no casounidimensional, v Vd v V , implicando que v satisfaz a formulao fraca contnua.

O problema aproximado , ento: achar ud Vd tal qued

gradv , gradud dV =d

v f dV v Vd. (3.7)

A maior diculdade que surge no problema bidimensional a manipulao numrica das funesda 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 discretizao. Mostraremos ainda que, assim comono caso unidimensional, o sistema tem nica soluo. Na seo seguinte escreveremos explicitamenteuma base e faremos algumas contas, j com vistas implementao de um algoritmo de MEF.

Seja, pois, B uma base do espao Vd. Sabemos que ud, por estar em Vd, tem (nica) representaocomo combinao linear das funes de B. Denotando essas funes por j , escrevemos

ud(x, y) =

Nj=1

j j(x, y) , (x, y) d,

onde N a dimenso de Vd.Substituindo ud acima no problema discretizado (3.7) obtemos

Nj=1

j

d

gradv , gradj dV =d

v f dV v Vd.

(Note que o integrando est dentro do somatrio.)Em particular, para v = i qualquer da base:

Nj=1

j

d

gradi , gradj dV =d

i f dV (3.8)

Variando i de 1 a N , (3.8) se mostra um sistema linear N N de incgnitas j . DenindoMij =

d

gradi , gradj 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.Proposio 2. A matriz M simtrica e positiva denida.Demonstrao.Como o produto interno simtrico por denio,Mij =

d

gradi , gradj dV =

d

gradj , gradi dV = Mji M simtrica.

10

Para mostrar que M positiva denida temos que provar que w RN {~0} se tem wTMw > 0.Ora,

wTMw =Ni=1

wi

Nj=1

wj

d

gradi , gradj dV =d

Ni=1

wi gradi ,Nj=1

wj gradjdV > 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 tivssemos gradi = 0 ie em todo d. Como por hiptese i contnua e cumpre a condio de fronteira de Dirichlet, issoimplicaria que as funes da base so identicamente nulas, um absurdo. Logo, wTMw estritamentepositivo, denindo M postivamente.

Como conseqncia do fato deM ser positiva denida, temos que o sistema (3.9) sempre admite umanica soluo, os coecientes j que, juntamente com as funes da base, determinam a aproximaoud de u.

Agora que vimos como o problema de Dirichlet se escreve na forma discreta usando o mtodo doselementos nitos, usaremos uma base de Vd para estudar como fazer os clculos e efetivamente resolvero problema.

3.3 Uma base para Vd

3.3.1 Enumeraes dos vrtices

Antes de buscarmos uma base para o espao Vd, introduziremos alguns conceitos que se mostraroimportantes naquela tarefa. Denimos o nmero de vrtices da malha como sendo o nmero totalde vrtices dos elementos, com a condio de que mesmo se determinado vrtice comum a vrioselementos, o contamos apenas uma vez. Chamamos de vrtices interiores aqueles que no esto sobrea fronteira de d (Figura 4).

Figura 4: Uma malha com 16 elementos e 13 vrtices, sendo que destes, 5 so interiores.

Vamos supor que o nosso domnio foi dividido em m tringulos, resultando em uma malha de Nvrtices, sendo que N so interiores. Faremos algumas denies.Denio 1. Chamamos de enumerao dos elementos a uma bijeo que associa a cada elementotriangular da malha um nmero natural entre 1 e m. Representamos, pois, cada elemento pela letra Tseguida de seu nmero como sub-ndice. Por exemplo, Tk o k-simo elemento da malha.Denio 2. Chamamos de enumerao global dos vrtices interiores a uma bijeo que associa acada vrtice interior da malha um nmero natural entre 1 e N . Representamos, pois, cada vrticeinterior pela letra p seguida de seu nmero como sub-ndice. Por exemplo, pi o vrtice interior i damalha.

Denio 3. Chamamos de enumerao global dos vrtices a uma bijeo que associa a cada vrticeda malha um nmero natural entre 1 e N , respeitando a enumerao global dos vrtices interiores. Estaenumerao consiste em adotar a enumerao da Denio 2 e ainda atribuir nmeros entre N + 1 eN aos vrtices da fronteira de d.Denio 4. Chamamos de enumerao local dos vrtices a uma bijeo quei) associa a cada vrtice de um elemento Tk um nmero do conjunto {1, 2, 3};

11

ii) percorre o elemento em sentido anti-horrio. Isto , denido o vrtice nmero 1 do elemento Tk,percorre-se a fronteira do elemento em sentido anti-horrio a partir desse vrtice 1. O prximo vrtice

ser o de nmero 2 e o ltimo ser o nmero 3. O vrtice s, s em {1, 2, 3}, do elemento Tk temcoordenadas (x

(k)s , y

(k)s ).

Notemos que um mesmo vrtice comum a dois elementos pode ter numerao local diferente emcada elemento. Por exemplo, pode ser o vrtice 1 do elemento Tk e o vrtice 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 vrtice seja o

vrtice interior h global, ento (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 possvel enumerao (global) doselementos e dos vrtices. A Tabela 1 complementa a Figura 5 exemplicando uma enumerao localdos vrtices. Note que, para cada elemento, um dos vrtices globais assume a posio local 1, 2 ou 3.Na Figura 6 mostramos alguns elementos e a enumerao local de seus vrtices.

Figura 5: Exemplo de enumerao dos elementos (a), e enumerao global dos vrtices (b).

Tabela 1: Exemplo de enumerao local dos vrticesElemento Vrtice 1 Vrtice 2 Vrtice 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 enumerao localde seus vrtices. O nico cuidado nessaenumerao que seu sentido seja anti-horrio.

3.3.2 Funes da base

Com os conceitos de enumeraes globais e local dos vrtices 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 funes j : R2 R tais que

j(pi) =

{1 se i = j,0 se i 6= j (3.10)

e

grco de j no elemento Tk =

{plano 6= 0 se Tk tem o vrtice pj ,0 caso contrrio

(3.11)

formam uma base B de Vd. Lembramos a notao: os vrtices interiores estamos representando pelospontos pi. Tanto i quanto j acima assumem valores em {1, 2, 3, ..., N}. Essas funes tm formatospiramidais, como se v na Figura 7. Diremos que a funo j e o vrtice pj so associados. Note queuma funo associada a um nico ponto, e vice-versa.

Figura 7: Funo chapu.

As nicas funes de B que assumem valores no-nulos em um dado elemento so aquelas trs asso-ciadas aos seus vrtices. Diremos que essas funes so associadas ao elemento, que reciprocamentelhes associado. Note que a cada elemento podem existir no mximo trs funes associadas, mas queuma funo pode ser associada a um nmero qualquer de elementos, dependendo da triangulao damalha.

Por (3.10) j sabemos quanto vale j(x, y) se o ponto (x, y) for um vrtice de um elemento. Nossatarefa agora determinar o valor que a funo assume num ponto no interior de um tringulo. Isto ,determinar j(x, y) para qualquer (x, y) d.

J sabemos por (3.11) que se o ponto (x, y) no pertence a nenhum elemento que tenha o vrticepj , ento j(x, y) = 0. Vejamos, ento, o que ocorre se (x, y) Tk e o tringulo Tk for associado funo j .

Como vimos, existe uma enumerao local dos vrtices de Tk: (x(k)1 , y

(k)1 ), (x

(k)2 , y

(k)2 ) e (x

(k)3 , y

(k)3 ).

Vamos supor, ento, que a funo j B a associada ao vrtice 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 vrtices, escrevendo simplesmente,(x1, y1). (Faremos isso apenas para deixar a notao mais limpa durante a deduo da frmula; 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, ento (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 dafuno nos vrtices, (x1, y1, 1), (x2, y2, 0) e (x3, y3, 0).

Para que o grco de j(x, y) esteja nesse plano, os trs 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 trs deve ser nulo. Esses vetores so:

(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 relao terceira coluna:

(j 1)x x2 y y2x x3 y y3

j x x1 y y1x x3 y y3+ j x x1 y y1x x2 y y2

= 0 j(x, y) =

x x2 y y2x x3 y y3x x2 y y2x x3 y y3

x x1 y y1x x3 y y3+ x x1 y y1x x2 y y2

. (3.12)Faremos algumas manipulaes algbricas usando propriedades dos determinantes para escrever a

equao acima de maneira mais conveniente. Por exemplo, no numerador:

x x2 y y2x x3 y y3 = x x2 y y2x y

x x2 y y2x3 y3 = x yx y

x2 y2x y x yx3 y3

+x2 y2x3 y3 =

1 x y1 x2 y21 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 restitudos os ndices superiores.Um resultado da Geometria Analtica informa que o determinante do denominador acima justa-

mente o dobro da rea do tringulo de vrtices (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 trs vetores a, b e

c (nesta ordem) pertencentes a R3 denido por a , b c e pode ser calculado como o determinante da matriz cujaslinhas so a, b, e c, nesta ordem. A interpretao geomtrica desse produto o volume do paraleleppedo determinadopelos trs vetores. Caso o resultado seja nulo, os trs vetores no determinam volume algum, estando, pois, num mesmoplano.

14

A expresso acima fornece o valor de j num ponto (x, y) qualquer do elemento Tk. Em outroselementos associados, a funo pode no ser dada por (3.14). Relembremos as suposies feitas queresultaram em (3.14): consideramos que j era associada Tk, mais precisamente, era associada aovrtice 1 de Tk. Sob essas hipteses, encontramos a frmula acima para j neste elemento.

Neste ponto dever o leitor estar se perguntando o que ocorreria se a funo j fosse associada a

outro vrtice r, (x(k)r , y(k)r ), local que no o de nmero 1. Neste caso, o vrtice r o que cumpre o

papel de vrtice 1 na equao acima. O vrtice seguinte a r cumprir o papel de vrtice 2 e o ltimo,o de vrtice 3. imporante que nos lembremos que a numerao local dos vrtices sempre feita emsentido anti-horrio: a ordem local dos vrtices sempre 1 2 3 1 2 3 1 (o vrtice 1sempre segue ao 3; o 3 segue ao 2, que segue ao 1).

Por exemplo, se j associada ao vrtice 2 do elemento Tk, entovrtice 2 vrtice 1 em (3.14);vrtice 3 vrtice 2 em (3.14) evrtice 1 vrtice 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 associao ao vrtice 3 e assim chegamos numa ex-

presso geral para o valor de uma funo 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 vrtice 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 vrtice 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 vrtice 3 de Tk,0 se (x, y) Tk mas j no for associada a Tk.

(3.15)

Note o carter local da expresso acima: j denida elemento por elemento. s vezes, para

deixar bem claro que estamos calculando a funo restrita ao elemento Tk, escreveremos (k)j . Em

cada elemento j poder ser dada por uma expresso diferente, dependendo das coordenadas dos seusvrtices e do nmero local do vrtice associado. Da mesma forma, o gradiente de j depender doelemento considerado.

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

3.3.3 Clculo dos gradientes

Supondo que j associada ao vrtice 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)2x(k)3 y(k)3 x

1 y(k)21 y(k)3+ y

1 x(k)21 x(k)3)

(k)j

x(x, y) = 1

2Ak

1 y(k)21 y(k)3 e

(k)j

y(x, y) =

1

2Ak

1 x(k)21 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 funo j seja a associada ao vrtice 2 de Tk, basta fazer a troca de ndices na expressoacima. O vrtice seguinte ao associado ser o 3 (posio em (3.16) ocupada pelo vrtice 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 funo associada ao vrtice 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 no associada a nenhum vrtice do tringulo Tk, ento grad(k)j (x, y) = 0.

As expresses (3.16), (3.17) e (3.18) do o valor do gradiente de uma funo da base em um elementose ela lhe for associada ao vrtice 1, 2 ou 3, respectivamente. Poder a leitora se perguntar: digamosque j associada ao vrtice 1 do elemento Tk mas tambm associada ao vrtice 3 de outro elemento,

Tl. Neste caso, quanto vale grad(l)j ? Ora, como essa funo associada ao vrtice 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 vrtices dos elementos, podemoscalcular os gradientes de todas as N funes de B em todos os elementos da malha. A partir da, possvel 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 clculo deM, cujas entradasso Mij =

d

gradi , gradj dV . Reparemos que a integral sobre d se decompe em uma soma de

integrais, cada uma sobre um elemento. Isto :

Mij =

d

gradi , gradj dV =mk=1

Tk

grad(k)i , grad(k)j dV.

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

Mij =mk=1

(grad(k)i , grad

(k)j

Tk

dV

)=

mk=1

grad(k)i , grad(k)j Ak. (3.19)

Da forma como denimos as funes de B, elas s tm valores e gradientes no-nulos nos seuselementos associados. Como conseqncia disso, o produto grad(k)i , grad

(k)j s ser diferente de

zero se ambas funes i e j forem associadas Tk. Isso faz com que a maior parte dos termos dosomatrio (3.19) sejam nulos. Mais ainda, um grande nmero de entradas de M so zeros: a matrizde rigidez esparsa.

Conclumos que para calcular os Mij basta considerar os elementos associados ao mesmo tempo ai 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 tambm para o vetor de carga.Seus elementos so do tipo

d

i f dV =mk=1

Tk

(k)i f dV.

Novamente, apenas os termos do somatrio com Tk associado a i sero no-nulos.O que queramos mostrar como transformamos a integral sobre d em uma soma de integrais

sobre os elementos relevantes no clculo. De certo modo nossa exposio do problema de Dirichletbidimensional est terminada. Apenas a ttulo de ilustrao dos procedimentos, calcularemos a matrizde rigidez em dois exemplos de malhas.

3.4 Alguns exemplos de malhas

Nesta seo consideraremos alguns exemplos de triangulaes para mostrar como feita a construodo sistema (3.9). As malhas que mostraremos so bastante simples, com poucos elementos, j quedesejamos apenas ilustrar o mtodo. Em aplicaes prticas um nmero bem superior de elementosdeve ser utilizado. Nossas malhas podem ser consideradas clulas 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 domnio em nove quadrados delado L/3 e, em seguida, traamos uma diagonal em cada quadrado, maneira da Figura 8a.

Figura 8: Malha e enumerao global dos elementos (a) e dos vrtices (b). (c) mostra a enumerao localdos vrtices dos elementos que sero utilizados neste exemplo.

Nossa malha tem 18 elementos e 16 vrtices, sendo que apenas 4 so interiores. Enumeraremos osvrtices globalmente como mostra a Figura 8b, e localmente consoante a Figura 8c. As reas de todoselementos so iguais, e representaremos simplesmente por A.

Notemos que cada funo da base (associadas aos vrtices 1, 2, 3 e 4) associada a seis elementos.Mais ainda, esse conjunto funo-base + os seis elementos associados forma uma espcie de clula,sendo transladado equivale aos outros conjuntos semelhantes. Isso conseqncia da regularidadedesta malha.

A dependncia de 1 com os seis elementos vizinhos ao vrtice p1 a mesma de i qualquer comseus seis elementos associados. Assim, calculando o gradiente de 1 em T1, T2, T3, T10, T9 e T8, jteremos os gradientes das outras funes.

Comeemos, pois, pelo elemento T1. 1 associada ao seu vrtice 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 tambm associada ao vrtice (local) 1. Ento, 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 vrtice (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 no 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 grco de 1 uma pirmide, fcil ver que a direo de crescimento de 1 a mesma nessesdois elementos. Como o gradiente tem justamente essa direo, nesses elementos eles so paralelos. Jque T10 simtrico com relao p1 ao elemento T1, seus gradientes tm mesmo mdulo e sentidosopostos. Logo,

grad(10)1 = grad(1)1 =

L

6A(1, 1).

Da mesma forma, T9 simtrico 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, ento, os gradientes de 1. Podemos calcular o primeiro termo da matriz de rigidez:

M11 =

grad1 , grad1 dV =T1

L6A

(1,1) , L6A

(1,1) dV +T2

L6A

(0, 1) ,L

6A(0, 1) dV+

+

T3

L6A

(1, 1) , L6A

(1, 1) dV +T10

L6A

(1, 1) , L6A

(1, 1) dV +T9

L6A

(0,1) , L6A

(0,1) dV+

+

T8

L6A

(1,1) , L6A

(1,1) dV = L2

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

5L2

18A.

Como a rea total do domnio L2 e ele est dividido em 18 tringulos de mesma rea A, temosque A = L2/18. Substiuindo isso no resultado acima, M11 = 5.

Por simples inspeo da malha, e nos baseando nas consideraes j feitas sobre a simetria datriangulao 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 ento todos os gradientes. Nos elementos que no esto relacionados na lista acima, osgradientes so nulos.

Da ento, M11 = M22 = M33 = M44 = 5 .Por inspeo da malha, vemos que as funes

(a) 1 e 2 tm os elementos associados T3 e T10 em comum;(b) 1 e 3 tm os elementos associados T8 e T9 em comum;

18

(c) 1 e 4 tm os elementos associados T9 e T10 em comum;

(d) 2 e 3 no tm elementos associados em comum;(e) 2 e 4 tm os elementos associados T10 e T11 em comum;

(f) 3 e 4 tm os elementos associados T9 e T16 em comum.Notemos, ainda com base na malha, que as relaes (e),(b) e (f),(a) tm mesma geometria. Isso

mostra que basta analisar uma clula da malha e como esta se relaciona com suas vizinhas para entendero comportamento de toda a malha - claro, no caso de uma triangulao regular.

Feitas essas consideraes, podemos calcular os demais elementos de M usando (3.19):

(a) M12 = M21 = Agrad(3)1 , grad(3)2 + Agrad

(10)1 , grad

(10)2 = Agrad

(3)1 , grad

(1)1 +

Agrad(10)1 , grad(8)1 = 2A(4)

L2

36A2= 2L

2

9A= 4.

(b) M13 = M31 = Agrad(8)1 , grad(8)3 + Agrad

(9)1 , grad

(9)3 = Agrad

(8)1 , grad

(2)1 +

Agrad(9)1 , grad(3)1 = 2A(2)

L2

36A2= L

2

9A= 2.

(c) M14 = M41 = Agrad(9)1 , grad(9)4 + Agrad

(10)1 , grad

(10)4 = Agrad

(9)1 , grad

(1)1 +

Agrad(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 24 5 0 22 0 5 42 2 4 5

.Vemos queM acima tem poucos zeros, em contradio com o que h pouco armamos, que a matriz

de rigidez esparsa. Essa aparente incoerncia ocorre devido ao tamanho da malha considerada. Osnicos zeros de esparsidade que ocorrem so devidos s funes 2 e 3 (que no tm elementosassociados em comum). No entanto, se considerssemos uma malha formada com o mesmo padro,porm com nove pontos interiores, o nmero de funes sem elementos associados em comum aumentarbastante. Um pouco de reexo bastar para que o leitor se convena de queM para a malha mostradana Figura 9 ser dada por

M =

5 4 0 2 2 0 0 0 04 5 4 0 2 2 0 0 00 4 5 0 0 2 0 0 02 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 nmero de vrticesinteriores da malha aumentar, essa proporo tambm o far.

19

Figura 9: Malha de 32 elementos e 9 vrtices interiores (enumerados).

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

Por m, chamamos ateno para o fato de que uma mudana no sistema de enumerao dos vrticesinteriores causa alterao na matriz de rigidez (pois ocorre uma reordenao da base). Por exemplo, seao escrever a matriz para o caso de quatro vrtices interiores tivssemos usado a numerao da Figura10 ao invs da Figura 8, chegaramos na matriz M 6= M (verique):

M =

5 4 2 24 5 2 02 2 5 42 0 4 5

.

Figura 10: Malha de 18 elementos com outra enumerao dos vrtices interiores.

3.4.2 Exemplo 2

Considere = [0, L] [0, L], o quadrado de lado L. Dividimos nosso domnio em quatro quadrados delado L/2 e, em seguida, traamos as duas diagonais em cada quadrado, maneira da Figura 11a. Amalha formada tem 16 elementos e 13 vrtices, sendo que 5 so interiores. Enumeramo-los globalmenteconsoante a Figura 11b. Na Figura 11c mostramos enumeraes locais de vrtices em alguns tringulos.Veremos que s precisaremos desses elementos para escrever a matriz de rigidez.

20

Figura 11: Enumerao dos elementos (a), enumerao global dos vrtices (b) e enumerao local dosvrtices de alguns elementos (c).

A grande diferena dessa malha para a do Exemplo 1 que, enquanto l cada funo era associadaa seis elementos, aqui existem funes associadas a quatro e a oito elementos. Por inspeo da malha,vemos que as funes que so associadas a 4 elementos esto, por assim dizer, encerradas; que elasno tm nenhum elemento associado em comum. A nica funo que tem elementos em comum comoutras a associada a oito elementos, 3. A partir dessa anlise 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 possvel, utilizaremos da simetria da malha parasimplic-las.

Consideremos primeiro 1. Ela no-nula apenas em T1, T2, T3 e T4. Pela geometria da malhavemos que T1 e T3 tm gradientes de mesmo mdulo e sentidos opostos, da mesma forma que T2 e T4.Ento, como T1 e T2 tm p1 como vrtice 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 inspeo da malha, vemos que essa estrutura de funo da base com quatro elementosao redor se repete pela malha, por uma simples translao. 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 vrtices simtricos com relao p3 tm gradientes de mesmo mdulo e sentidosopostos.

Podemos numerar, localmente, os vrtices de T3, T4, T5 e T8 de modo que o vrtice 1 seja semprep3. Isso nos faz usar apenas a frmula (3.16) para clculo 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 poderamos esperar isso, pois o

grco de 3 uma pirmide 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 funes 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 consideraes feitas, possvel 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 02 2 4 2 20 0 2 4 00 0 2 0 4

. (3.20)

Um pouco de reexo 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 02 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 ateno para o fato de que a matriz de rigidez depende da numerao global dos vrticesinteriores, mas independe da enumerao local. Em casos de malhas simtricas podemos escolher estaltima de modo a usar apenas uma frmula para o gradiente, como zemos neste exemplo.

22

Figura 12: Malha com 36 elementos e 13 vrtices interiores (enumerados).

3.5 Outros casos

Como dissemos logo no incio deste texto, nosso objetivo apenas transmitir a essncia do mtodo doselementos nitos, por isso este material bastante simples. Mencionamos aqui, brevemente, algumasoutras possibilidades que o Mtodo permite.

Contemplamos apenas os casos de uma e duas dimenses. A formulao do caso tridimensionalpode ser deduzida sem grandes diculdades a partir da deduo feita neste captulo. Foi, inclusive,com esse intuito que deixamos o Teorema do Divergente enunciado em sua forma geral. Novamente, adiculdade ir surgir ao discretizar o domnio (agora em tetraedros) e buscar escrever uma base parao espao de funes Vd.

Ao longo deste texto, sempre aproximamos a funo u por outra que tinha a propriedade deser linear em cada elemento da malha. Existem outras possibilidades: podemos desejar que ud sejaquadrtica por partes, ou mesmo polinomial por partes, fornecendo aproximaes mais suaves.

Por m, o Mtodo no se aplica apenas ao problema de Dirichlet. Condies de contorno deNeumann e Robin tambm so aceitas, com algumas alteraes no exposto neste texto. Por exemplo,no caso de condio de Neumann, no podemos, nas regies de com essa condio, igualar o membro direita de (3.5), frmula de Green, a zero.

23

Bibliograa

[1] Jochen ALBERTY, Carsten CARSTENSEN, Stefan A. FUNKEN, Remarks around 50 lines ofMatlab: short nite element implementation. Numerical Algorithms 20 (1999), 117-137.

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

[3] Giovanni CALDERN, Rodolfo GALLO, Introduccin al Mtodo de los Elementos Finitos: unenfoque matemtico. Caracas: Ediciones IVIC, 2011.

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

[5] James STEWART, Clculo: volume 2. So Paulo: Cengage Learning, 2009.

24