Métodos Numéricos em Fluido Térmica - feb.unesp.br · a) E.D.O. - Sistema Massa Mola b) E.D.P. - Eq. da quantidade de movimento Figura 1.1: Exemplos de E.D. ordinárias e parciais

UNIVERSIDADE ESTADUAL PAULISTA

Campus Universitario de Bauru

FACULDADE DE ENGENHARIA DE BAURU

DEPARTAMENTO DE ENGENHARIA MECANICA

Topicos Especiais em Fluido-Termica

Metodos Numericos em Fluido Termica

Build:08-2007 Autor: VICENTE LUIZ SCALON

SUMARIO

1 Equacoes Diferenciais e Metodos Numericos 1

2 Introducao aos Sistemas Matlab/GNU Octave 7

2.1 Operacoes Fundamentais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Definicoes e operacoes com matrizes e vetores . . . . . . . . . . . . . . . . . . . 8

2.2.1 Funcoes e operacoes especiais . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3 Definicao de funcoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4 Montagem de graficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.5 Operacoes logicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.6 Diferencas basicas entre o Matlab e o Octave . . . . . . . . . . . . . . . . . . . . 18

3 Metodos para a solucao de Eq. Diferenciais Ordinarias 19

3.1 Metodo de Runge Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.1.1 Implementacao da solucao no Gnu-Otave . . . . . . . . . . . . . . . . . . 26

i

SUMARIO ii

3.2 Solucoes de sistemas de equacoes diferenciais ou equacoes diferenciais de ordem

superior utilizando Runge Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.2.1 Procedimento de solucao usando o octave . . . . . . . . . . . . . . . . . . 34

3.3 Condicao de Contorno Deslocada . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.3.1 Procedimento de solucao usando o octave . . . . . . . . . . . . . . . . . . 40

4 Solucao de Sistema de Equacoes Lineares 42

4.1 Metodos iterativos de solucao . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.2 Metodos de inversao direta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5 Equacoes Diferenciais Parciais - Caracterizacao e Metodos de Discretizacao

e Solucao 48

5.1 Subdivisao das equacoes diferenciais parciais . . . . . . . . . . . . . . . . . . . . 48

5.2 Formas de solucao de equacoes diferenciais parciais . . . . . . . . . . . . . . . . 50

5.3 Discretizacao do domınio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

6 Metodo dos Elementos Finitos 55

6.1 Princıpios gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

6.1.1 Aproximacao por funcoes . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

6.1.2 Aproximacao nodal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

6.2 Aproximacao por Elementos Finitos . . . . . . . . . . . . . . . . . . . . . . . . . 61

6.2.1 Definicao geometrica dos elementos. . . . . . . . . . . . . . . . . . . . . . 62

6.2.2 Regras para a discretizacao de domınios atraves de elementos finitos . . . 62

6.3 Elementos de referencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

SUMARIO iii

6.4 Propriedades das funcoes de aproximacao u(x) . . . . . . . . . . . . . . . . . . . 66

6.5 Construcao das funcoes de interpolacao . . . . . . . . . . . . . . . . . . . . . . . 67

6.6 Transformacao de operadores diferenciais . . . . . . . . . . . . . . . . . . . . . . 70

6.6.1 Montagem do Jacobiano . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.6.2 Transformacao de uma integral . . . . . . . . . . . . . . . . . . . . . . . 72

6.7 Coordenadas nodais e conectividade . . . . . . . . . . . . . . . . . . . . . . . . . 73

6.8 O metodo dos resıduos ponderados . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.8.1 Transformacao integral: a integral por partes . . . . . . . . . . . . . . . . 78

6.8.2 Metodo de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6.9 Tratamento das condicoes de contorno . . . . . . . . . . . . . . . . . . . . . . . 80

6.9.1 Implementacao da condicao de contorno generalizada . . . . . . . . . . . 81

6.10 Integracao Numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.11 Exemplos de Aplicacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7 Metodo de Diferencas Finitas 94

7.1 Aplicacao da formulacao de diferencas finitas em sistemas de equacoes diferenciais 97

8 Aplicacao do Metodo de Diferencas Finitas em Problemas Bidimensionais e

Transientes 102

8.1 Solucao de problemas bidimensionais . . . . . . . . . . . . . . . . . . . . . . . . 102

8.2 Solucao de problemas transientes . . . . . . . . . . . . . . . . . . . . . . . . . . 106

8.2.1 Formulacao Explıcita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

8.2.2 Formulacao Totalmente Implıcita . . . . . . . . . . . . . . . . . . . . . . 111

SUMARIO iv

8.2.3 Formulacao Implıcita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

CAPITULO 1

Equacoes Diferenciais e Metodos Numericos

Uma grande parcela dos problemas de engenharia depende, para a obtencao de resulta-

dos, da solucao de uma unica ou ate de um sistema de equacoes diferenciais. Ate a primeira

metade deste seculo buscou-se, de forma intensa, a solucao analıtica destas equacoes utilizando

uma grande gama de ferramentas matematicas como transformadas, solucao por series, etc.

Com este esforco concentrado foram obtidos alguns resultados e muitas equacoes diferenciais

puderam ser resolvidas. No entanto, quase todas estas respondem por problemas fısicos simples

e nao representam uma amostra significativa dos problemas de engenharia que normalmente

tem geometria e condicoes de contorno complexas.

Desta forma, a partir da segunda metade do seculo XX, este panorama foi completamente

modificado. Deixou-se de buscar a solucao puramente analıtica para estes problemas e passou-

se a trabalhar com os metodos numericos na tentativa da obtencao de solucoes aproximadas.

Embora a grande maioria dos metodos numericos sejam conhecidos a bastante tempo, a sua

utilizacao em massa so ocorreu gracas a um acontecimento: o grande avanco dos computadores.

No passado quando se tentava resolver um problema, mesmo um extremamente simples, atraves

de metodos numericos este processo demandava um enorme tempo em calculos e, muitas vezes,

envolvia uma equipe de pessoas. Esta dificuldade praticamente inviabilizava a sua utilizacao ate

o aparecimento dos computadores de grande porte e, mais recentemente, os super computadores.

Esta situacao se inverteu completamente e nos dias de hoje sao os metodos numericos que

1

Top. Esp. Fluido-Temica 2

respondem praticamente pela totalidade dos problemas complexos em engenharia. Alem disto

o barateamento ocorrido em termos de custo/hora do tempo de CPU tem tornado-a, cada vez

mais, acessıvel a um numero maior de pessoas. Um exemplo disto e citado por (Maliska, 1995)

em seu livro: a solucao de um escoamento supersonico usando computadores do tipo IBM 704,

existente na decada de 60, consumiria um tempo de computacao de 30 anos a um custo de

alguns milhoes de dolares enquanto o mesmo problema e resolvido hoje com alguns minutos de

CPU a um custo de algumas centenas de dolares.

Apesar de tudo isto que foi dito aqui os metodos analıticos ainda tem grandes utilidades.

Em primeiro lugar, nao tem muito sentido resolver um problema numericamente se a sua

solucao analıtica e conhecida e pode fornecer a valores para todo o domınio, ao contrario

da solucao numerica que so so fornece a solucao para os pontos considerados. Alem disto

a solucao analıtica e, muitas vezes, utilizada como padrao e ate mesmo fonte de inspiracao

na resolucao de problemas com metodos numericos. Em outras vezes ainda, um modelo e

testado em determinadas condicoes nas quais existe solucao analıtica para poder confrontar os

resultados e, so depois do sucesso nesta consideracao, e que o mesmo e estendido para casos

onde a solucao analıtica nao e conhecida.

Quando nao existe solucao analıtica que permita esta verificacao do metodo numerico outra

ferramenta e entao utilizada: os ensaios experimentais. Da mesma forma que os metodos

analıticos os metodos experimentais foram, e diria ate que ainda sao, muito utilizados na

solucao de problemas de engenharia. Isto se deve ao fato que a simulacao pode ser feita sobre

as condicoes desejadas (ou proximas dela) sem se preocupar sequer com as condicoes contorno ou

tipo de equacoes diferenciais. Embora os metodos experimentais tenham os seus atrativos, tem

tambem um grande problema: o seu elevado custo. Tanto a instrumentacao necessaria, como

a montagem do experimento e a mao de obra necessaria necessitam de um alto investimento,

o que dificulta grandes desenvolvimentos nesta area. Alem do mais a presenca de sondas para

tomada de medidas ou quaisquer outros objetos fısicos, tendem a alterar as condicoes ideais para

a retirada de medidas de um experimento, tornando muitas vezes o experimento impraticavel.

Um grande exemplo deste fato e a industria aeronautica que tem substituıdo grande parte dos

seus ensaios em tunel de vento por simulacoes numericas com excelentes resultados. Da mesma

forma que a solucao analıtica muitas vezes sao utilizados resultados experimentais para validar

resultados de metodos numericos, dispensando a partir de entao a repeticao da experiencia para

casos similares.

Para a solucao de equacoes diferenciais por metodos numericos, da mesma forma que quando

sao resolvidas por metodos analıticos, e conveniente o conhecimento de alguns preceitos sobre

equacoes diferenciais:


a) E.D.O. - Sistema Massa Mola b) E.D.P. - Eq. da quantidade de movimento

Figura 1.1: Exemplos de E.D. ordinarias e parciais.

Equacao Diferencial Ordinaria: onde a funcao que e a solucao da equacao diferencial e

dependente de uma unica variavel. O exemplo classico deste tipo de equacao e a primeira

lei de Newton, comumente aplicada em problemas envolvendo vibracoes, como mostrado

na figura (??a):

mu′′(t) = F [t, u(t), u′(t)]

onde todas as variaveis sao direta ou indiretamente funcao de t. Um grande numero de

problemas recai em equacoes diferenciais deste tipo.

Equacao Diferencial Parcial: e aquela em que as funcoes solucao para a equacao diferencial

sao dependentes de mais de uma variavel. A maior parte dos problemas envolvendo

geometrias bi e tri dimensionais e analise de transientes em corpos com gradientes internos

recai neste tipo de equacao diferencial. Tomando como exemplo utilizar-se-a a equacao

da quantidade de movimento linear, bidimensional e em coordenadas cartesianas, para

um fluido newtoniano na direcao x:

ρ

(∂u

∂t+ u

∂u

∂x+ v

∂u

∂y

)= −∂P

∂x+ µ

(∂2u

∂x2+∂2u

∂y2

)Um exemplo de aplicacao desta equacao seria para a obtencao do perfil de velocidades

em torno de uma asa, como mostrado na figura (??b).

Equacao Diferencial Linear: e o tipo de equacao diferencial onde aparece uma composicao

linear da funcao e suas derivadas. Exemplo disto e a equacao da conducao de calor, para

substancias isotropicas e com propriedades fısicas admitidas independentes da tempera-


a) E. D. Linear - Conduca numa placa b) E.D. nao linear - Eq. do pendulo

Figura 1.2: Exemplos de E.D. linear e nao linear.

tura:∂2T

∂x2+∂2T

∂y2= 0

A distribuicao de temperaturas em uma placa plana mostrada na figura (??a), e um

exemplo tıpico deste tipo de aplicacao.

Equacao Diferencial Nao Linear: e aquela a funcao nao aparece de forma linear na equacao

diferencial, mas sim com termos quadraticos ou atraves de outras funcoes nao lineares. Um

exemplo deste tipo de equacao e a equacao de um pendulo que, com base na figura (??b),

pode ser dada como:d2θ

dt2+g

lsin θ = 0

que tem formas de solucao bem complexa. Prova disto e a propria expressao acima, que

e quase sempre resolvida atraves da linearizacao da equacao (trabalhando-se com angulos

pequenos sin θ ≈ θ):d2θ

dt2+g

lθ = 0

Equacao Diferencial Homogenea: e aquela em que o termo independente nao se apresenta

como funcao de nenhuma variavel, ou melhor:

y′(x) = f(x, y)

e homogenea se f(x, y) for uma constante ou uma expressao que permita uma trans-

formacao de variavel que a deixe independente da variavel transformada (normalmente

um f(y/x)). Um exemplo de equacao diferencial homogenea:

y′′(x) = y′(x) + y(x) =⇒ y′′(x)− y′(x)− y(x) = 0


Equacao Diferencial Nao Homogenea: e a contraposicao da da definicao acima ou seja,

quando o termo independente e necessariamente funcao de uma das variaveis do problema

e nao se conhece transformacao de variaveis que o deixe homogeneo. Exemplo de equacao

diferencial nao homogenea:

y′′(x)− y′(x)− y(x) = cos(x)

Ordem de uma Equacao Diferencial: representa o ındice de derivacao da maior derivada

das que compoe a equacao diferencial. Veja por exemplo a expressao:

y′′′ + 2 y′ − 3y = 0

que e uma equacao de terceira ordem e que, consequentemente, precisa de pelo menos

tres condicoes de contorno para ser resolvida. Esta e na realidade a grande vantagem

desta metodologia de classificacao de equacoes diferenciais: indicar instantaneamente o

numero de condicoes de contorno necessarias para resolver o problema.

Estas classificacoes indicam caracterısticas principais da equacoes diferenciais e fornecem

informacoes sobre como estas podem ser resolvidas. Nao existe uma forma imediata para

afirmar que uma equacao diferencial tem solucao exata pois ate mesmo o tipo de condicao

de contorno pode influir sobre este aspecto. Um exemplo disto e a equacao de conducao bi-

dimensional que tem solucao analıtica para temperaturas de parede dadas, mas nao tem se

uma das paredes esta sujeita a um processo de conveccao, por exemplo. Assim sendo a unica

maneira de descobrir sobre a existencia da solucao analıtica para uma dada equacao diferencial,

caso voce nao saiba de antemao, e consultando um livro sobre o assunto.

Mas e importante ressaltar que as caracterısticas acima sao importantes tambem para iden-

tificar qual o metodo numerico de solucao que deve ser utilizado na solucao de um problema.

Vejamos por exemplo o tipo de solucao mais utilizados segundo a natureza das equacoes dife-

renciais:

Eq. Diferenciais Ordinarias

Metodo de Euler

Metodo de Passos Multiplos (ou Preditor-Corretor)

Metodo de Runge Kutta1

Eq. Diferenciais Parciais

Metodo de Diferencas Finitas

Metodo de Volumes Finitos

Metodo de Elementos Finitos

1Muitas vezes este metodo tambem e utilizado na solucao da parcela envolvendo o tempo de Eq. DiferenciaisParciais


E importante ressaltar que nada impede que se utilize esquemas de solucao para Eq. Dife-

renciais Parciais em Eq. Diferenciais Ordinarias, uma vez que o seu procedimento e generico,

no entanto a recıproca nao e verdadeira. Normalmente se adotam esquemas diferenciados para

solucao de equacoes diferenciais ordinarias, pois estes podem ser um pouco mais simples ou

entao mais precisos permitindo uma maior acuracidade na solucao.

CAPITULO 2

Introducao aos Sistemas Matlab/GNU Octave

Existem uma serie de ambientes matematicos propıcios para a solucao de algumas tarefas a

serem realizadas cotidianamente em calculos da Engenharia: Matlab, Mathemathica, GNU Oc-

tave, SciLab, Maxima, etc. Alguns destes sao capazes, inclusive, de trabalhar com manipulacao

simbolica como o caso do Maxima, Mathemathica, Matlab (versoes posteriores a 5.0), SAGE

e o proprio octave se utilizando de pacotes adicionais.Entretanto, para o caso de utilizacao em

simulacao numerica a manipulacao simbolica nao represbta um fator decisivo.

Este capıtulo, basicamente, ficara restrito ao uso dos sistemas Matlab/GNU Octave sendo

o primeiro um sistema licensiado e o segundo uma alternativa livre de ambientes matematicos.

Embora similares em grande numero de comandos existem algumas diferencas entre os co-

mandos em cada um dos sistemas. Na maioria das vezes octave suporta tanto a sua sintaxe

especıfica como aquela que seria utilizada pelo Matlab. O SciLab tambem e considerado uma

boa alternativa livre ao uso do Matlab, mas o seu uso nao sera abordado neste material.

Existem uma serie de referencias que podem complementar as informacoes aqui fornecidas,

dentre as quais destaco os materiais de (Domingues & Mendes, 2002) e (Eaton, 2006). Outra

importante fonte de ajuda e o proprio programa, onde uma serie de informacoes a respeito de

um comando podem ser obtidas utilizando-se help -i nome_do_comando.

Inicialmente, sera visto simplesmente algumas operacoes fundamentais com matrizes e ve-

7


tores que nao apresentam variacao entre estes sistemas. Com estas informacoes ja sao possıveis

realizar uma serie de procedimentos do nosso curso.

2.1 Operacoes Fundamentais

Neste tipo de plataformas estao contemplados todos os tipos de operadores, tanto para

operacao com reais com inteiros. Assim sao possıveis a soma(+), subtracao(-), divisao(/),

multiplicacao (*) divisao reversa (\) e exponencial (ˆ). Operacoes com inteiros sao tambem

possıveis como a divisao, utilizando o truncamento dos decimais (floor), e resto (mod ou rem).

Assim:

octave>mod(5,2)ans =

1

octave> rem(5,2)ans =

1

octave> floor(5/2)ans =

2

octave> disp(5**2), disp(’ ou ’),disp(5^2)25ou25

O disp e um comando utilizado para escrever na tela e converte a saida para caracteres.

Comandos para arredontamento como round ou ceil tambem estao disponıveis no Octave.

Alem disto, existe uma extensa biblioteca matematica pre-implementada que permite o

calculo de uma serie de funcoes hiperbolicas (exp, log , sinh, etc.), trigonometricas (sin, cos ,

tan, etc.), de Bessel (besselj , besselk , besseli , etc.) e uma infinidade de outras.

2.2 Definicoes e operacoes com matrizes e vetores

Antes de mais nada e possıvel criar vetores e matrizes atraves de um valor inicial, um valor

final e incrementos constantes do tipo:


octave> 1:10ans =

1 2 3 4 5 6 7 8 9 10

octave> 1:2:10ans =

1 3 5 7 9

ou entao se estabelecendo nao o incremento, mas sim o numero de componentes da matriz:

octave> linspace(1,10,5)ans =

1.0000 3.2500 5.5000 7.7500 10.0000

Para criar uma matriz ou um vetor incluindoos valores de cada posicao e armazena-lo numa

variavel, o procedimento tambem e simples, basta inseri-lo da maneira mostrada abaixo:

octave> a=[1 2; 4 7]a =1 24 7

ou ainda utilizando um <enter>, ao inves do ”; ”, para indicar mudanca de linha:

octave> b=[3 6> 9 4]b =3 69 4

Definidas as matrizes pode-se realizar operacoes entre elas. Veja por exemplo como realizar

uma adicao entre as matrizes a e b, definidas anteriormente.

octave> a+bans =

4 813 11

Da mesma maneira pode-se utilizar uma resposta anterior, mesmo que nao armazenada em

variavel nenhuma utilizando da variavel ans. Como exemplo disto, veja como apresentar a

segunda coluna da matriz resposta anterior:


octave> c=ans(:,2)c =

811

sendo que para isto e bastante util o ”:”da maneira apresentada. Ele pode representar, quando

usado desta maneira, todas as linhas ou colunas de uma matriz. Caso desejasse mostrar apenas

um componente da matriz, bastaria colocar o seu endereco ente parenteses:

octave> a(2,1)ans = 4

Da mesma maneira que a adicao, outras operacoes entre as matrizes poderiam ser realizadas,

como por exemplo a multiplicacao:

octave> a*bans =21 1475 52

Outra forma desta operacao, a chamada multiplicacao termo a termo, pode tambem ser

necessaria e neste caso ela pode ser realizada atraves da forma:

octave> a.*bans =

3 1236 28

sendo ainda existente uma operacao equivalente a esta para a divisao termo a termo, represen-

tada pelo operador ”.”.

2.2.1 Funcoes e operacoes especiais

Sao ainda possıveis uma serie de outras operacoes com matrizes, sendo destacadas aqui:

• Determinante (det):


octave> det(a)ans = -1

• Matriz Inversa (inv):

octave> inv(b)ans =-0.095238 0.1428570.214286 -0.071429

• Matriz Transposta (’):

octave> b’ans =3 96 4

• Matriz nula de qualquer tamanho (zeros):

octave> zeros(4)ans =0 0 0 00 0 0 00 0 0 00 0 0 0

ou ainda para qualquer matriz nao quadrada definindo-se o numero de linhas e colunas:

octave> zeros(1,7)ans =0 0 0 0 0 0 0

• Matriz Unitaria tambem pode ser montada de dorma analoga (ones):

octave:1> ones(3,2)ans =1 11 11 1

• Matriz de numeros aleatorios (rand): com todos os numeros aleatorios variando entre 0

e 1.

octave> rand(2,4)ans =0.539648 0.061666 0.070065 0.3248830.569649 0.023215 0.673922 0.419023


Em funcao do exposto se o interesse e por uma matriz cujo o valor maximo e 10, basta

mutiplicar o resultado anterior pelo valor maximo.

• Matriz Identidade de qualquer tamanho (eye):

octave> eye(4)ans =1 0 0 00 1 0 00 0 1 00 0 0 1

que e uma operacao bastante util se voce estiver interassado em montar uma linha qual-

quer com um valor ”1”na posicao da diagonal principal e o restante zeros:

octave> eye(10)(5,:)ans =0 0 0 0 1 0 0 0 0 0

• Matriz Diagonal generica a partir de um vetor (d iag):

octave> a=[1 2 3]a =1 2 3

octave> diag(a)ans =1 0 00 2 00 0 3

O vetor diagonal tambem pode ser usada para montar uma diagonal secundaria da matriz,

para isto basta fornecer como segundo argumento inteiro que representa a sua posicao

na matriz. Numeros negativos podem ser usados para representar diagonais secundarias

abaixo da posicao atual:

octave> diag(a,-2)ans =0 0 0 0 00 0 0 0 01 0 0 0 00 2 0 0 00 0 3 0 0

Se aplicado em uma matriz bidimensional, o comando diag retorna a respectiva diago-

nal indicada na forma de vetor, como se fizesse uma operacao inversa da anteriormente

demonstrada:


octave:18> b=[8 5 5 6 7; 0 3 4 1 5; 1 2 2 4 0;3 9 7 8 6 ]b =8 5 5 6 70 3 4 1 51 2 2 4 03 9 7 8 6

octave:19> diag(b,1)’ans =5 4 4 6

• Operacoes com as colunas de componentes de uma matriz: no caso da soma (sum)

octave> sum(a)5 9

e ainda existem outros comandos que permitem a obtencao da media(mean), o produto

dos termos(prod), o valor maximo (max ), o valor mınimo (min) e a ordenacao de matrizes

(sort). Todos estes comandos realizam estas operacoes entre os elementos pertencentes a

mesma coluna.

Deve-se lembrar ainda que mesmo nos ambientes deste tipo nao existe a comutatividade em

operacoes com matrizes assim:

octave> c*c’ans =64 8888 121

e diferente de:

octave> c’*cans = 185

como era de se esperar.

Bem este texto serve como uma referencia basica para o tratamento de matrizes e vetores

nos referidos sistemas entretanto existem ainda uma serie de diferentes comandos relacionados

a este que podem ser encontrados em documentacoes mais aprofundadas e atraves do Help dos

programas.

Existem uma serie de outras operacoes que permitem operacoes basicas com vetores, prin-

cipalmente com relacao a uniao de vetores (union) e a idendificacao de posicoes que obedeccam

a caracterısticas definidas (find).


2.3 Definicao de funcoes

Para definir funcoes no octave nomalmente e indicado criar um arquivo com extensao .m no

diretorio corrente obedecendo a uma estrutura basica:

i. a primeira linha deve conter a palavra chave function, em seguida a variavel que armazena

o valor a ser retornado que, por sua vez, e igualada ao nome da funcao seguida da sequencia

de parametros de entrada. E fundamental que o nome da funcao seja identico ao fornecido

ao arquivo .m.

ii. na linha a seguir sao definidas as variaveis globais, se existirem.

iii. depois vem o corpo da funcao com a sua sequencia de comandos.

iv. o procedimento e finalizado com a palavra end .

Veja por exemplo a criacao de uma funcao do tipo sinal de um numero. Desta forma sera

editado um arquivo sinal.m do tipo:

# func~ao sinalfunction ret=sinal(x)if (x!=0) ret=x/abs(x);else ret=0;endifend

A partir deste ponto existe uma funcao pronta no octave de nome sinal que pode ser chamada

em qualquer instante. Cabe ressaltar entretanto que esta funcao deve estar no diretorio corrente

ou no diretorio de funcoes do octave. Assim:

octave> sinal(100)ans= 1

octave> sinal(-10)ans=- 1

octave> sinal(0)ans = 0

Esta mesma funcao poderia ser criada simplesmente digitando a sequencia de comandos apre-

sentada no octave dispensando, assim, a necessidade da criacao de um novo arquivo. O inco-

veniente desta forma e que a mesma so estaria disponıvel depois de carregada para a memoria

do octave em cada secao.


Quando se trata de funcoes mais simples, que envolve o seu calculo diretamente a partir

de parametros fornecidos o comandoinline pode ser uma boa alternativa. Sofre das mesmas

limitacoes de quando se define a funcao no interior de um script, entretanto sua definicao e bem

mais simples:

octave> f=inline("2*x.^2-3*x+4")f =f(x) = 2*x.^2-3*x+4

octave> f(2)ans = 6

sendo que neste caso todos os parametros envolvidos na funcao seriam tambem argumentos da

mesma dificultando, assim, definicoes mais complexas. Existem alternativas para personalizar

esta definicao uma vez que este comando pode ser utilizado com maior numero de parametros.

Maiores detalhes podem ser encontrados com a utilizacao da ajuda da funcao.

2.4 Montagem de graficos

Para elaboracao de graficos o octave se utiliza de um programa externo denominado ”GNU-

PLOT”. Existem comandos internos do proprio gnuplot que muitas vezes sao utilizados para

definir parametros preliminares dos graficos. Para um bom conhecimento destas funcoes sugere-

se uma leitura do manual do proprio programa.

Com relacao ao comando para plotagem plot ele pode ser utilizado com a entrada de pelo

menos dois vetores (x, y), mas sua forma geral permite a utilizacao de um formato em sequencia

identificando como vai ser a linha

Apenas para ilustrar, foi feito um grafico personalizado alterando algus parametros mais

importantes do gnuplot via gset e utilizando-se de um script do octave:

octave> x1=0:0.1:pi; %define vetor xoctave> a=cos(x1); %define o primeiro vetor yoctave> b=sin(x1); %define o segundo vetor yoctave> __gnuplot_set__ xlabel "x" % define nome do eixo xoctave> __gnuplot_set__ ylabel "y" % define nome do eixo yoctave> __gnuplot_set__ key outside box % define legenda

% do lado de fora do grafico e com bordaoctave> plot(x1,a,"-;cos(x);",x1,b,"-;sin(x);")


Figura 2.1: Grafico gerado no octave

e com isto foi criado o grafico mostrado na figura (2.1). Deve-se ressaltar que em versoes antigas

do programa utilizava-se substiturir o comando __gnuplot_set__ por gset.

Um dos aspectos mais complexos e a utilizacao de estilos de linhas e pontos neste tipo

de plotagem. Alem de escolher o tıtulo da legenda da curva e ainda possıvel nestes graficos,

escolher tanto a cor como a forma das linhas ou pontos da curva. Para tanto e interessante

conhecer os esquemas a serem utilizados:

• - define a curva na forma de linhas;.

• ., +, *, o e x define a curva na forma de diferentes estilos de pontos;

• ^ define grafico de impulsos;

• L define grafico de ”steps”;

• ‘n’ ou ‘c’, definem a cor a ser utilizada de forma:

Num. Letra Cor

0 k preto

1 r vermelho

2 g green

3 b azul

4 m magenta

5 c cyan

Alem destas cores pode-se ainda utilizar o w para o branco e os numeros maiores que

cinco para outras variacoes.


Cabe ressaltar ainda que o comando fplot pode ser utilizado diretamente para a elaboracao

de graficos a partir de funcoes diretamente. Ele pode ser utilizada de maneira analoga ao plot ,

excetuando-se pelas mudancas de formatos anteriormente demonstradas.

octave> __gnuplot_set__ xlabel "x" % define nome do eixo xoctave> __gnuplot_set__ ylabel "y" % define nome do eixo yoctave> fplot("[cos(x), sin(x)]", [0,pi])

2.5 Operacoes logicas

E possıvel realizar uma serie de operacoes logicas e testes usando o Octave. As operacoes

mais usuais sao maior (>) ou maior ou igual (>=), menor(<) ou menor ou igual (<=), igual

(==) e diferente (! = ou ˜=). E conveniente notar que o teste de igualdade (==) e diferente

da atribuicao (=).

• if e utilizado para realizacao de comparacoes diretas e direcionar o fluxo do programa em

funcao de seu resultado.

octave> a=22

octave> if (mod(a,2)==0) disp(\"Numero par\") else disp(\"Numero impar\") endifNumero par

octave> a=3;3

octave> if (mod(a,2)==0) disp(\"Numero par\") else disp(\"Numero impar\") endifNumero impar

• while utilizado para o caso de repeticoes onde o teste e feito por diversas vezes a cada

interacao do problema.

octave> z=1;octave> while (z<5) disp(z); z+=2; endwhile13

• for no caso de operacoes que usam um contador com incrementos constantes o comando

for e o mais indicado.

octave> for z=1:2:4 disp(z); endfor13


• switch permite a selecao de uma alternativa entre diversas. Pode ser substituido por um

conjunto de if s em cascata.

octave> nlados=3;octave> switch (nlados)> case (3) disp("Triangulo")> case (4) disp("Quadrado")> case (5) disp("Pentagono")> otherwise disp("Figura n~ao classificada")> endswitchTriangulo

Em todos os comandos acima o final endif , endwhile, endfor e endswitch pode ser substi-

tuido por end sem comprometer o funcionamento do script (e mantendo compatibilidade com

o Matlab)

2.6 Diferencas basicas entre o Matlab e o Octave

Algumas diferencas basicas que podem afetar a compatibilidade entre ambos sao:

• o nome de algumas funcoes sao diferentes

• o comentario no Matlab e % e no Octave tambem e aceito o #

• no Matlab, os blocos formados por while, if e for e as “functions“ necessariamente ter-

minam com end . No octave pode-se usar, opcionalmente, endwhile, endif , endfor e

endfunction.

• no Matlab a unica forma aceita para a desigualdade e o ˜=. O != e aceito apenas no

Octave.

• operadores incrementais ++ e - - nao sao aceitos no Matlab.

CAPITULO 3

Metodos para a solucao de Eq. Diferenciais

Ordinarias

Os problemas envolvendo as Equacoes Diferenciais Ordinarias sao normalmente subdi-

vididos em dois grandes grupos:

Problemas de Valor Inicial: que sao aqueles onde as condicoes de contorno sao estabeleci-

das em um ponto inicial, e partir destas vao sendo calculados os valores para posicoes

subsequentes. Este procedimento e normalmente caracterizado por procedimento de mar-

cha (os transientes sao os melhores exemplos deste tipo de problemas).

Problemas de Valor de Contorno: sao aqueles onde as condicoes de contorno sao estabe-

lecidas em posicoes diferentes dentro de um problema, o que implica em uma solucao que

envolva todo o domınio.

Para exemplificar melhor cada um dos problemas, imagine um caso de conducao unidimen-

sional em um solido de condutividade termica independente da temperatura com geracao de

calor. A equacao diferencial para este problema e:

T ′′(x) = − qk

19


Imaginando que a uma das temperaturas das faces e conhecida e T (x = 0) = T0. Se alem

dela, for conhecido o fluxo de calor que esta entrando por esta face ou seja, na mesma posicao

q(x = 0) = kT ′(0) = q0 este problema sera um problema de valor inicial e o procedimento de

marcha podera ser adotado.

Por sua vez se, ao inves disso, forem conhecidas as duas temperaturas de parede em ambas

as faces T (x = 0) = T0 e T (x = L) = TL este problema passa a ser um problema de valor de

contorno.

3.1 Metodo de Runge Kutta

Existem diversos graus para o metodo de Runge Kutta, sera mostrada aqui a deducao

para o metodo de segunda ordem, para os metodos de ordem superior o procedimento e analogo.

Sabe-se que por serie de Taylor e possıvel expressar o valor qualquer ponto a partir de um

ponto conhecido:

y(x+ ∆x) = y(x) + ∆x y′(x) +∆x2

2!y′′(x) + · · ·

que pode ser escrito para um sistema discreto com espacamento ∆x = h na forma:

yn+1 = yn + h y′n +h2

2!y′′n + · · ·

Suponhamos que uma dada equacao diferencial pode ser manipulada de forma a obter uma

expressao para a derivada primeira de forma que:

y = f ′(x, y) =⇒ y′n = f(xn, yn)

No caso da derivada de uma funcao qualquer em relacao a x e conveniente lembrar da regra

da cadeia:d

dxf(x, y) = fx(x, y) + fy(x, y)

dy

dx

sendo que os ındices x ou y da funcao representam a derivada parcial da funcao em relacao a

tal derivada.

Reescrevendo a expansao por Taylor, mas levando em conta esta informacao tem-se agora

uma expansao para funcoes de duas variaveis :

yn+1 = yn + hf(xn, yn) +h2

2![fx(xn, yn) + fy(xn, yn)y′n] + · · ·


ou ainda:

yn+1 = yn + hf(xn, yn) +h2

2![fx(xn, yn) + f(xn, yn)fy(xn, yn)] + · · · (3.1)

O metodo de Runge Kutta se baseia no princıpio de que existe um ponto em um determinado

intervalo cuja a derivada fornece exatamente a tangente para o calculo do ponto no extremo

deste intervalo. Esta derivada pode ser expressa como uma composicao linear na forma:

yn+1 = yn + (a1y′n + a2y

′∗)h (3.2)

sendo y′∗ um valor previsto para a derivada para o ponto generico (xn + b1h, yn + b2y′nh).

Assim o valor para y′∗:

y′∗ = f(xn + b1h, yn + b2y′nh)

que pode ser ainda expandida numa serie de Taylor de primeira ordem (envolvendo os termos

de h)

f(xn + b1h, yn + b2y′nh) = f(xn, yn) + b1fx(xn, yn)h+ b2y

′nhfy(xn, yn)

= f(xn, yn) + b1fx(xn, yn)h+ b2hf(xn, yn)fy(xn, yn)

Substituindo o valor da expansao de y′∗ na equacao (3.2) tem-se que:

yn+1 = yn + {a1f(xn, yn) + a2[f(xn, yn) + b1fx(xn, yn)h+ b2hf(xn, yn)fy(xn, yn)]}h

Que rearranjada resulta em:

yn+1 = yn + (a1 + a2)f(xn, yn)h+ a2b1fx(xn, yn)h2 + a2b2f(xn, yn)fy(xn, yn)h2 (3.3)

Finalmente, comparando-se as equacoes (3.1) e (3.3) obtem- se o seguinte sistema de

equacoes:

a1 + a2 = 1

a2b1 =1

2

a2b2 =1

2

Este sistema tem mais incognitas que equacoes e a sua famılia de solucoes e dada por:

a2 = 1− a1

b1 =1

2− 2a1

b2 =1

2− 2a1


e

e

e

eK1

2K2

2

K3

xn

h

h/2

e

xn+1

K4

uyn+1

Figura 3.1: Esquema de funcionamento do Runge Kutta de 4a ordem

Uma solucao comum para o sistema seria dada se a1 = a2 que implica em a1 = a2 = 1/2 e

b1 = b2 = 1 que igualaria ao metodo de Runge Kutta ao Metodo de Euler. Maiores detalhes

desta deducao e do metodo de Euler podem ser encontrados em (Ruggiero & Lopes, 1988).

Exercıcio: Resolva a equacao diferencial para um problema de conducao onde a geracao de

energia e funcao da temperatura dada por:

dT

dx= x+ 2T

utilizando o Runge Kutta de 2a ordem numa placa de 2 m (0 ≤ x ≤ 2) e usando um incremento

h = 0.4. A temperatura numa face da parede e dada e e igual 10◦C. (x = 0 → T = 10).

A figura (3.1) mostra graficamente como funciona a aproximacao mais utilizada deste tipo

de esquema: O Runge Kutta de 4a ordem. Sao aproximadas tanto as retas tangentes como os

pontos utilizados para os calculos intermediarios. E apresentada tambem na figura o valor do

ponto calculado numericamente utilizando as escalas dos Ks.

As expressoes que sao utilizadas para o Runge Kutta de 4a ordem na aproximacao de um

problema y′(x) = f(x, y) sao:


K1 = h f(xn, yn)

K2 = h f(xn + h/2, yn + K1/2)

K3 = h f(xn + h/2, yn + K2/2)

K4 = h f(xn + h, yn +K3)

sendo que o valor posterior e calculado na forma:

yn+1 = yn +1

6(K1 + 2K2 + 2K3 +K4) (3.4)

Exemplo de Aplicacao: Determinar o tempo de resposta de um termometro de mercurio

de vidro quando este instrumento e utilizado na leitura de um ambiente cuja temperatura varia

de forma senoidal com o tempo: T∞ = 40 + 22.48 sin(2πt) sendo o tempo t dado em horas.

Admitir que o coeficiente global de transferencia de calor entre o termometro e o fluido e de

25 kcal/h m2◦C. O termometro pode ser idealizado como um cilindro de mercurio de 25 mm

de comprimento e 6 mm de diametro e sua temperatura inicial de 15◦C.

Dados:

ρmerc = 13600 kg/m3

cmerc = 0.0325 kcal/kg ◦C

Volmerc = 7.07× 10−7 m3

Solucao:

Pelo balanco termico tem-se que:

[Qtde de calor

armazenada notermometro

]=

Qtde de calortransferida pelo

fluido

ρmerc cmerc VolmercdT

dt= hA (T∞ − T ) (3.5)

sendo que a area de transferencia de calor por conveccao e: A = πdL = π(0.006)(0.025) =

4.71× 10−4 m2.

Substituindo os valores na equacao (3.5) tem-se:

13600× 0.0325× (7.07× 10−7)dT

dt= 25× (4.71× 10−4) (T∞ − T )


ou de forma rearranjada tem-se:

dT

dt= 37.681× (T∞ − T )

Como o valor de T∞ = 40 + 22.48 sin(2πt):

dT

dt= 37.681[40 + 22.48 sin(2πt)− T ] (3.6)

que e uma expressao na forma T ′ = F (t, T ).

A condicao inicial do problema tambem e conhecida e e dada por t = 0 ⇒ T = 15◦C.

Para resolver este problema basta aplicar o metodo de Runge-Kutta ou entao resolver esta

expressao analiticamente. A solucao exata para este problema pode ser encontrada em (Kreith,

1977) e e dada por:

T = 40 + 22.17 sin(2πt− 0.165)− 21.36 exp(−37.17 t) (3.7)

Aplicando-se Runge Kutta para o problema em intervalos de tempo de 2 min (h = 0.03333 h)

tem-se que:

• t = 0 : T=15◦C

• t = 0.03333

K1 = hF (0, 15) = 0.03333{37.681[40 + 22.48 sin(2π0)− 15]} = 31.401

K2 = hF (0 + 0.03333/2, 15 + 31.401/2) = 0.03333F (0.01667, 30.700) = 14.632

K3 = hF (0 + 0.03333/2, 15 + 14.632/2) = 0.03333F (0.01667, 22.316) = 25.163

K4 = hF (0 + 0.03333, 15 + 25.163) = 0.03333F (0.03333, 40.163) = 5.666

T2 = T1 +1

6(K1 + 2K2 + 2K3 +K4) ⇒

T2 = 15 +1

6(31.401 + 2.14.632 + 2.25.163 + 5.66)

isto implica em que T2 = 34.44◦C.

Utilizando a solucao analıtica para o problema a temperatura indicada pelo termometro,

para o mesmo tempo, e de 34.80◦C. Isto equivale a um erro de cerca de 1%, indicando a

boa precisao do metodo.


A tabela (3.1) mostra a solucao para este problema para varios valores de tempo. Ela

fornece os valores da temperatura calculada pela solucao analıtica e por Runge Kutta, alem de

fornecer o valores para os Ks em cada tempo.

A figura (3.2) mostra o comportamento das temperaturas lidas pelo termometro, exata e

numerica, e a temperatura real do banho. Pela figura e possıvel mostrar a boa concordancia

entre o resultado numerico e o resultado exato.

Tabela 3.1: Tabela de resultados para o exemplo 1

Tempo (h) Texata Tnumerica K1 K2 K3 K4

0.00000 14.999 15.0000.03333 34.797 34.443 31.401 14.632 25.163 5.6660.06667 43.776 43.560 12.851 7.635 10.910 4.7610.10000 49.389 49.285 7.013 5.242 6.354 4.1440.13333 53.665 53.613 4.934 4.132 4.636 3.4980.16667 57.075 57.043 3.885 3.305 3.669 2.7460.20000 59.661 59.635 3.046 2.475 2.834 1.8880.23333 61.365 61.340 2.192 1.580 1.964 0.9520.26667 62.129 62.104 1.277 0.630 1.036 −0.0250.30000 61.924 61.900 0.317 −0.344 0.071 −0.9990.33333 60.760 60.737 −0.653 −1.302 −0.894 −1.9300.36667 58.688 58.668 −1.593 −2.202 −1.820 −2.7770.40000 55.800 55.783 −2.464 −3.006 −2.666 −3.5020.43333 52.222 52.208 −3.227 −3.679 −3.395 −4.0750.46667 48.109 48.099 −3.849 −4.191 −3.976 −4.4690.50000 43.641 43.637 −4.302 −4.519 −4.383 −4.6680.53333 39.015 39.015 −4.568 −4.651 −4.599 −4.6620.56667 34.432 34.437 −4.634 −4.578 −4.613 −4.4530.60000 30.092 30.102 −4.497 −4.306 −4.426 −4.0500.63333 26.185 26.199 −4.164 −3.846 −4.046 −3.4690.66667 22.882 22.900 −3.649 −3.217 −3.488 −2.7370.70000 20.327 20.348 −2.974 −2.448 −2.778 −1.8850.73333 18.632 18.654 −2.170 −1.572 −1.947 −0.9510.76667 17.870 17.894 −1.270 −0.627 −1.031 0.0250.80000 18.076 18.100 −0.315 0.345 −0.070 0.9990.83333 19.240 19.263 0.653 1.302 0.895 1.9300.86667 21.312 21.332 1.594 2.202 1.820 2.7770.90000 24.200 24.217 2.464 3.006 2.666 3.5020.93333 27.778 27.792 3.227 3.679 3.395 4.0750.96667 31.891 31.901 3.849 4.191 3.976 4.4691.00000 36.359 36.363 4.302 4.519 4.383 4.668


10

15

20

25

30

35

40

45

50

55

60

65

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Tem

pera

tura

[oC

]

Tempo[h]

Sol. AnaliticaSol. Numerica

Temperatura do banho

Figura 3.2: Temperaturas do termometro e banho no exemplo 1

3.1.1 Implementacao da solucao no GNU-Otave

Para arepresentar o problema sera escolhida a solucao da relacao entre a temperatura de

um termometro e de seu banho apresentada na apostila. Neste caso a equacao do problema e

dada por:

T ′(t, T ) = 37.681(40 + 22.48 sin(2πt)− T )

sendo a condicao de partida a temperatura inicial do termometro T (0) = 15.

Para definir esta equacao do problema e preciso criar um arquivo representando esta funcao,

sendo que ambos (a funcao e arquivo) devem ter o mesmo nome. No caso sera criado o arquivo

dert.m, composto por:

function dr=dert(tempo,temper)dr=37.681*(40+22.48*sin(2*pi*tempo)-temper);endfunction

e a partir dele e possıvel calcular o valor da funcao para qualquer par ordenado T ′(t, T ). Por

exemplo:

octave> dert(0.15)ans = 942.02


Cabe ressaltar aqui que no Matlab nao existe o comando endfunction, sendo que o mesmo

nao precisa ser incluido no caso de arquivos deste tipo.

Entretanto o objetivo principal e a solucao desta equacao diferencial, portanto e preciso

resolve-la. Para tanto sera montado uma outra funcao runge.m na qual sera implementado o

metodo de Runge-Kutta. A funcao referida que utilizara a dert.m, anteriormente definida.

function resp=runge(tempo,temper,h)k1=h*dert(tempo,temper);k2=h*dert(tempo+h/2, temper+k1/2);k3=h*dert(tempo+h/2, temper+k2/2);k4=h*dert(tempo+h, temper+k3);resp=temper+(k1+2*k2+2*k3+k4)/6;endfunction

Com esta funcao agora e possıvel realizar a marcha do processo de solucao. Caso se desejasse

a temperatura apos decorridos 2 min, e possıvel obte-la usando a funcao recem definida:

octave> runge(0.15,2/60)ans = 34.443

Para a obtencao da solucao completa e preciso repetir este procedimento por diversas vezes

este procedimento e armazenar a solucao em vetores. Para isto, a sequencia de comandos abaixo

pode ser implementada diretamente ou via arquivo e a solucao armazenada nos valores de te

e TT:

dt=2/60;te=0:dt:1;TT=zeros(1,31);TT(1)=15;for i=2:31

TT(i)=runge(te(i-1),TT(i-1),dt);endfor

que se executado no octave. Deve-se observar que o comando te0:2/60:1;= cria um vetor com

todos os termos de 0 a 1, incrementados em intervalos de 2/60.

Feita esta analise implementar ainda solucao analıtica do problema e a evolucao da tem-

peratura do banho em funcoes distintas (t exata.m e tbanho.m). O resultado obtido destas

funcoes sao armazenadas em variaveis:


function tex=t_exata(t)tex=40+22.17*sin(2*pi*t-0.165)-21.36*exp(-37.681*t);endfunction

e

function tba=tbanho(t)tba=40+22.48*sin(2*pi*t);endfunction

Definidas as funcoes pode-se avaliar as solucoes e armazena-los nas variaveis texv e tba,

a partir das quais, pode-se plotar os resultados. Estes valores estao mostrados no grafico em

funcao do tempo armazenado na variavel te, e o resultado esta mostrado na figura (3.3). Para

comparacao uma tabela de saıda de dados foi montada com base nos resultados obtidos. Desta

forma e possıvel uma analise dos valores numericos de cada caso:

octave> [te’ TT’ texv’ tba’]ans =

0.00000 15.00000 14.99853 40.000000.03333 34.44279 34.90197 44.673850.06667 43.55979 43.83592 49.143440.10000 49.28485 49.41487 53.213410.13333 53.61306 53.67464 56.705900.16667 57.04284 57.07827 59.468250.20000 59.63482 59.66188 61.379750.23333 61.34022 61.36521 62.356850.26667 62.10423 62.12881 62.356850.30000 61.89951 61.92357 61.379750.33333 60.73686 60.75968 59.468250.36667 58.66766 58.68836 56.705900.40000 55.78252 55.80022 53.213410.43333 52.20758 52.22153 49.143440.46667 48.09910 48.10870 44.673850.50000 43.63665 43.64147 40.000000.53333 39.01527 39.01510 35.326150.56667 34.43692 34.43177 30.856560.60000 30.10170 30.09180 26.786590.63333 26.19909 26.18487 23.294100.66667 22.89964 22.88172 20.531750.70000 20.34756 20.32672 18.620250.73333 18.65439 18.63154 17.643150.76667 17.89412 17.87027 17.643150.80000 18.09998 18.07617 18.620250.83333 19.26298 19.24024 20.531750.86667 21.33229 21.31162 23.294100.90000 24.21747 24.19977 26.786590.93333 27.79242 27.77847 30.856560.96667 31.90090 31.89130 35.326151.00000 36.36335 36.35853 40.00000


Figura 3.3: Grafico gerado no GNU-Octave

Neste grafico estao mostradas a solucao numerica e exata do problema, alem da temperatura

do banho. O comando utilizado para plotagem, considerando-se as avriaveis anteriormente

definidas e dado por:

xlabel "Tempo [horas]"ylabel "Tempertura [C]"plot(te,texv,"-;Sol Analitica;",te,TT,"*;Sol Numerica;",te,tban, "-;Temp. Banho;");

Solucao usando comandos preexistentes no GNU-Octave

Embora o procedimento acima possa ser realizado sem maiores problemas ele depende,

como foi mostrado, da elaboracao de uma rotina para a solucao do problema, no caso usando

o procedimento de Runge-Kutta. Existe uma alternativa um pouco mais simples que consiste

na utilizacao do procedimento de solucao ja implementado no octave usando o comando lsode.

Este comando consiste na utilizacao do algoritmo de Hindmarsh, um algoritmo um pouco

mais recente que o de Runge-Kutta e otimizado para sistemas de Equacoes Diferenciais Or-

dinarias. Para a solucao do problema anteriormente proposto e preciso conhecer o uso de:

lsode(”nome da funcao”, condicao, domınio)

sendo que:

nome da funcao e o nome do arquivo que contem a expressao da funcao a ser integrada e cujo

nome vem entre” . Embora esta seja basicamente igual a definida anteriormente deve-se


tomar o cuidado de que a funcao deve ser chamada sempre com o primeiro argumento

sendo o vetor da grandeza a ser calculada e o segundo o do parametro da solucao.

Nestas condicoes a funcao utilizada na solucao do problema do termometro, chamada de

dert2.m deve ser redefinida como:

function dr=dert2(temper, tempo)dr=37.681*(40+22.48*sin(2*pi*tempo)-temper);endfunction

condicao sao a(s) condicao(oes) de partida necessarias para a solucao do problema.

domınio representa a faixa de valores para os quais a solucao vai ser obtida. Neste caso deve-se

estabelecer a forma de um vetor do tipo inıcio: passo: fim ou forma equivalente.

Para o caso do exemplo anterior pode-se utilizar:

octave> sol2=lsode("dert2",15,te);

considerando o dominio do problema anteriormente definido por te. Depois disto a solucao

poderia ser plotada e comparada com as anteriores ou mesmo verificado a diferenca entre a

nova solucao e a anterior.

Um procedimento similar a este poderia ser elaborado usando-se o Matlab para obter a

solucao deste problema com as seguintes diferencas:

• existem uma serie de funcoes que permitem a solucao de odes no Matlab, sendo a mais

utilizada a ode45, que normalmente substituir o lsode.

• a chamada da funcao dentro do ode45 se da com um @ na frente do nome e nao entre ”.

• o domınio pode ser estabelecido atraves de um vetor na forma [inıcio final ], alem das

formas anteriormente apresentadas.

EXERCICIO: Repita o procedimento anteriormente apresentado agora para um termopar,

cuja capacidade termica e, por consequencia, o tempo de resposta sao bem menores. Supondo

que as condicoes do banho sao as mesmas propostas para o caso do termometro e que as

propriedades fısicas do termopar sao: ρ = 7600 kg/m3; c = 0.12 kcal/kg K; D = 0.0008

m.Considere neste caso que o volume submerso e igual ao diametro total do par termoeletrico.


Solucao: Neste caso a EDO do problema seria dada por:

T ′(t, T ) = 137.06(40 + 22.48 sin(2πt)− T )

e a solucao analıtica do problema seria dada por:

T (t) = 40 + 22.46 sin(2πt− 0.0458)− 23.97 exp(−137.06t)

Agora apresente a solucao numerica deste problema e compare-a com a resposta obtida a partir

da solucao analıtica.

3.2 Solucoes de sistemas de equacoes diferenciais ou equa-

coes diferenciais de ordem superior utilizando Runge

Kutta

O mesmo procedimento mostrado anteriormente pode ser estendido para equacoes di-

ferenciais de ordem superior, mas para isto e necessario algumas adaptacoes no esquema. A

adaptacao utilizada e converter a E.D. (Equacao Diferencial) em um sistema. Este fato e melhor

explicado tomando-se uma E.D. como exemplo:

y′′ + y′ − y + 2x = 0 =⇒ y′′ = −y′ + y − 2x

Esta equacao diferencial e analoga ao sistema abaixo:

y′ = Z

y′′ = −Z + y − 2x

O procedimento para a solucao de sistemas e identico ao metodo de solucao de uma unica

equacao desde que seja tomado o cuidado de calcular os Ks em paralelo para todas as equacoes.

Isto e muito importante e deve ser respeitado e implica que antes de calcular qualquer K2 os

valores de K1 devem ter sido calculados para todas as equacoes.

Para ilustrar a forma de utilizacao deste procedimento sera apresentado um exemplo a

seguir.

Exemplo de Aplicacao: Considere uma aleta circular sobre um duto de de 12 cm de

diametro e 0.5 cm de espessura e feita de alumınio (k = 215 W/m2◦C). A aleta troca ca-

lor com um ambiente a 25◦C e com coeficiente de pelıcula de 50 W/m2◦C. Sabendo que a base


da aleta trabalhara ha uma temperatura de 50◦C e que deve dissipar uma quantidade de calor

de 120 W pergunta-se qual o comprimento mınimo que esta aleta devera ter?

Obs: O calor trocado pela ponta da aleta pode ser desprezado.

Solucao:

Fazendo um balanco termico em um anel de espessura dr nesta aleta tem-se que:

qr = qr+dr + qconv

−k(2π r t) dTdr

∣∣∣∣r

= −k(2π (r + dr) t)dT

dr

∣∣∣∣r+dr

+ h 2 (2πr dr) (T − T∞)

Fazendo a a aproximacao da derivada (dT/dr)r+dr por serie de Taylor na forma:

dT

dr

∣∣∣∣r+dr

=dT

dr

∣∣∣∣r

+d2T

dr2

∣∣∣∣r

dr

e ainda simplificando a equacao obtem-se:

rdT

dr= (r + dr)

(dT

dr+d2T

dr2dr

)− 2h r dr

k t(T − T∞)

Note que os valores de avaliacao da derivada deixaram de aparecer uma vez que todas agora

passam a ser avaliadas na posicao r. Expandindo a expressao:

rdT

dr= r

dT

dr+ r

d2T

dr2dr + dr

dT

dr+d2T

dr2dr2 − 2h r dr

k t(T − T∞)

Desprezando-se o termo de ordem O(dr2) e simplificando a expressao:

d2T

dr2+

1

r

dT

dr− 2h

k t(T − T∞)

ou rearranjadad2T

dr2=

2h

k t(T − T∞)− 1

r

dT

dr

Substituindo pelos valores numericos:

d2T

dr2= 93.02 (T − 25)− 1

r

dT

dr(3.8)

Esta equacao esta sujeita as condicoes de contorno:

r = R, T = Tb na base da aleta a temperatura e conhecida. Substituindo os valores tem-se

que neste caso: T=50◦C.


r = R, q = qb o fluxo de calor na base tambem e conhecido. No entanto e importante ressaltar

que o fluxo de calor tem que ser adaptado para se tornar uma condicao de contorno:

qb = −k 2π R tdT

dr

∣∣∣∣r=R

=⇒ dT

dr

∣∣∣∣r=R

= − qb2π k R t

Substituindo os valores:

dT

dr

∣∣∣∣r=R

= − 120

2π 215 0.06 0.005= −296.10

Feito isto pode-se partir para a solucao da equacao propriamente dita. Para questao de

nomenclatura utilizar-se-a o sistema de equacoes na forma:

G(r, T,G) = G

F (r, T,G) = 93.02 (T − 25)− 1

rG

sendo que T ′′ = F (r, T,G) e T ′ = G(r, T,G)

• r = 0.06 T1 = 50; e G1 = −296.1.

• r = 0.07 adotar-se a o segundo sub-ındice 1 para os valores da integracao da funcao e o

ındice 2 para a integracao da sua derivada, assim

K1,1 = h ∗G(0.06, 50,−296.1) = 0.01 ∗ −296.1 = −2.96

K1,2 = h ∗ F (0.06, 50,−296.1) = 0.01 ∗ [93.02 (50− 25)− 1

0.06− 296.1] = 72.61

K2,1 = h ∗G(0.06 + 0.01/2, 50− 2.96/2,−296.1 + 72.61/2) = −2.60

K2,2 = h ∗ F (0.065, 48.52,−259.8) = 61.851

K3,1 = h ∗G(0.06 + 0.01/2, 50− 2.60/2,−296.1 + 61.85/2) = −2.65

K3,2 = h ∗ F (0.065, 48.70,−265.17) = 62.84

K4,1 = h ∗G(0.06 + 0.01, 50− 2.65,−296.1 + 62.84) = −2.33

K4,2 = h ∗ F (0.07, 47.35,−233.26) = 54.11

1Os valores dos parametros de F e G sao os mesmos para a mesma posicao e foram somente apresentadosde maneira diferente: sendo indicados em G e calculados em F .


E agora e possıvel calcular os valores desejados:

T2 = T1 +1

6(K1,1 + 2K2,1 + 2K3,1 +K4,1) =

50 +1

6(−2.96− 2 2.60− 2 2.65− 2.33) = 47.37

G2 = G1 +1

6(K1,2 + 2K2,2 + 2K3,2 +K4,2) =

−296.1 +1

6(72.61 + 2 61.85 + 2 62.64 + 54.11) = −233.41

Este mesmo procedimento pode ser repetido por diversas vezes ate que o fluxo de calor G se

anule, que e quando a aleta deixa de transferir calor. Os resultados para este problema podem

ser encontrados na tabela (3.2).

A solucao analıtica para este problema tambem pode ser obtida no entanto ele envolve

funcoes de Bessel e pode ser dada na forma:

T (r) = 25 + 5.472I0(9.645 r) + 23.648K0(9.645 r)

T ′(r) = 26.388[I−1(9.645 r) + I1(9.645 r)]− 114.04[K−1(9.645 r) +K1(9.645 r)]

sendo ainda que os valores calculados por estas funcoes tambem estao presentes na tabela (3.2).

Um ponto importante no que tange ao calculo da eficiencia da aleta e descobrir o ponto

onde a derivada se anula. Utilizando a solucao analıtica obtem-se r = 0.16408, e fazendo uma

regressao linear na tabela para os resultados numericos obtem-se r = 0.16413 que equivale a

um erro de 0.03%.

Deste resultado e possıvel avaliar a eficiencia da aleta:

η =qrealqideal

=120

50 2π (0.16412 − 0.062) (50− 25)=

120

183.22= 65.49%

3.2.1 Procedimento de solucao usando o GNU-Octave

Tambem no caso da EDO de segunda ordem e possıvel utilizar estes metodos. Como foi

visto a saıda e converter a solucao para um sistema de equacoes diferenciais em que cada

uma representa uma ordem diferente. Uma solucao implementando a tecnica de Runge Kutta

apresentada anteriormente poderia ser vista, entretanto, aqui so sera apresentada a tecnica que

se utiliza da funcao pre-implementada no octave lsode.


Tabela 3.2: Tabela de solucao para o exemplo resolvido 2

rT

Tex

ata

G=dT/dr

T′ ex

ata

k1,1

k1,2

k2.1

k2,2

k3,1

k3,2

k4,1

k4,2

0.06

50.0

050.0

0−

296.

10−

296.

10

0.07

47.3

747.3

7−

233.

42−

233.

42−

2.96

72.6

1−

2.60

61.8

5−

2.65

62.8

4−

2.33

54.1

1

0.08

45.2

845.2

8−

185.

70−

185.

70−

2.33

54.1

5−

2.06

47.2

3−

2.10

47.8

2−

1.86

42.0

6

0.09

43.6

243.6

2−

148.

02−

148.

02−

1.86

42.0

8−

1.65

37.3

8−

1.67

37.7

5−

1.48

33.7

5

0.10

42.3

042.3

0−

117.

38−

117.

38−

1.48

33.7

7−

1.31

30.4

4−

1.33

30.6

9−

1.17

27.8

2

0.11

41.2

641.2

6−

91.8

4−

91.8

4−

1.17

27.8

3−

1.03

25.4

0−

1.05

25.5

8−

0.92

23.4

6

0.12

40.4

540.4

5−

70.0

7−

70.0

7−

0.92

23.4

7−

0.80

21.6

6−

0.81

21.7

9−

0.70

20.2

1

0.13

39.8

539.8

5−

51.1

5−

51.1

5−

0.70

20.2

1−

0.60

18.8

4−

0.61

18.9

4−

0.51

17.7

4

0.14

39.4

239.4

2−

34.3

9−

34.3

9−

0.51

17.7

4−

0.42

16.7

0−

0.43

16.7

8−

0.34

15.8

7

0.15

39.1

539.1

5−

19.2

6−

19.2

6−

0.34

15.8

7−

0.26

15.0

8−

0.27

15.1

4−

0.19

14.4

5

0.16

39.0

339.0

3−

5.37

−5.

37−

0.19

14.4

5−

0.12

13.8

5−

0.12

13.9

0−

0.05

13.3

8

0.17

39.0

439.0

47.

617.

604

−0.

0513.3

90.

0112.9

50.

0112.9

90.

0812.6

1


Bem para isto e necessario definir uma funcao que represente a equacao diferencial do

problema que foi apresentada anteriormente, sendo:

T ′′(r, T, T ′) = 93.02(T − 25)− 1

rT ′

que sera representada atraves da funcao sol ale1.m, que tem como parametros de entrada um

vetor que armazena T e suas derivadas e uma variavel para armazenar a posicao radial r.

function temr=sol_ale1(temper, raio)temr=zeros(2,1);temr(1)=temper(2);temr(2)=93.02*(temper(1)-25)-1/raio*temper(2);endfunction

onde a variavel temr (1) e (2) representam as expressoes para a primeira e segundas derivadas

de T nos pontos considerados, respectivamente.

Feito isto a solucao pode ser obtida estabelecendo a regiao do domınio para a qual seria

solucionada (raios rr), as condicoes iniciais) e executando a chamada do lsode, como mostra

o script a seguir.

octave> saida=lsode("sol_ale1", [50 -296.1], rr=[0.06:0.01:0.2]);octave> xlabel "Tempo [horas]"octave> ylabel "Temperatura [C]"octave> plot(rr,saida(:,1),"-*; Sol. Numerica;")octave> figure(2)octave> xlabel "Tempo [horas]"octave> ylabel "Grad. Temper [C/h]"octave> gset key left topoctave> plot(rr,saida(:,2),"-*; Sol. Numerica;")

que apresentaria a solucao do problema na tela e o armazenaria na variavel saida. A partir

destes dados, poderia ser tracar os graficos da forma que se mostrasse adequada e tambem

identificar o tamanho real da aleta, a partir do ponto em que a derivada se anula. A figura (3.4),

mostra o comportamento da solucao.

3.3 Condicao de Contorno Deslocada

Ate agora foi visto apenas casos onde as condicoes iniciais eram conhecidas de antemao.

No entanto existem casos onde se deseja satisfazer condicoes que estejam deslocadas em relacao

ao ponto inicial. O procedimento de Runge Kutta funciona bem tambem neste caso, no entanto


Figura 3.4: Temperatura e gradientes no caso da aleta

e necessario utilizar-se de um procedimento iterativo repetindo-o para diferentes condicoes

iniciais ate obter a condicao deslocada desejada. Este procedimento e usado como uma maneira

de contornar as limitacoes do metodo utilizando-o na solucao de um problema de contorno.

Embora este acerto na condicao de contorno possa ser feito por tentativa e erro existem

procedimentos mais otimizados para a busca dos valores iniciais: sao eles os metodos de busca

de zero de funcoes (metodo da biparticao, secante e Newton-Raphson). O metodo que tende a

ser mais eficiente para este procedimento e o Metodo de da Secante, que e uma adaptacao do

metodo de Newton-Raphson. Para tanto e apresentada a seguir uma breve revisao a respeito

dos metodos.

Metodo de Newton Raphson: e um metodo que permite achar a raiz de expressoes. Sua

forma basica e:

xn+1 = xn −f(xn)

f ′(xn)

Metodo da Secante: e uma adaptacao do Metodo de Newton-Raphson que procura expres-

sar a derivada da funcao de forma discreta. Relembrando a definicao de derivada e:

f ′(x) = lim∆x→0

f(x+ ∆x)− f(x)

∆x

No caso da expressao pela Metodo da Secante e tomada a expressao discreta para a derivada


utilizando-se dos dois pontos anteriormente calculados, ou seja, fazendo ∆x = xn−xn−1 e assim:

xn+1 = xn −f(xn)

f(xn)− f(xn−1)

(xn − xn−1)

ou xn+1 = xn −f(xn)

f(xn)− f(xn−1)(xn − xn−1)

Esta mesma expressao pode ser reagrupada na sua forma mais usual que e:

xn+1 =f(xn)xn−1 − f(xn−1)xn

f(xn)− f(xn−1)(3.9)

Esta expressao permite o calculo de forma a obter valores estimados para as condicoes

iniciais no caso de condicao de contorno deslocada.

Exemplo de Aplicacao: Considere uma aleta circular sobre um duto de de 12 cm de

diametro e 0.5 cm de espessura e feita de alumınio (k = 215 W/m2◦C). A aleta troca ca-

lor com um ambiente a 25◦C e com coeficiente de pelıcula de 50 W/m2◦C. Sabendo que a base

da aleta trabalhara ha uma temperatura de 50◦C pergunta-se qual maximo fluxo de calor que

podera ser dissipado pela aleta, independente do seu comprimento?

Solucao:

A aleta tera a sua maxima dissipacao de calor se o seu comprimento tender a infinito. Mas

para que o seu comprimento tenda a infinito e preciso respeitar a condicao de contorno:

dT

dr

∣∣∣∣r→∞

= 0

No caso da solucao numerica, e preciso estipular um valor para este infinito, neste caso

estipulou-se como infinito o valor de 2 m.

A outra condicao de contorno e dada pela temperatura na base da aleta, ou seja T (r =

0.06) = 50◦C.

Como este problema e similar ao exemplo anterior tomou-se inicialmente a condicao de

contorno para a derivada, ja calculada anteriormente, T ′(r = 0.06) = −296.1. Depois dele

adotou-se um valor arredondado para o T ′(r = 0.06) = −300. Os valores obtidos para a

derivada em x = 2 m e mostrada a seguir:


dT/dx|x=0dT/dx|x=∞

2

-296.1 307617.567

-300 296887.049

Tomando por base estes dois valores e possıvel aplicar a expressao do metodo da secante,

equacao (3.9), e obter uma estimativa do valor que vai anular o fluxo afastado da base da aleta.

Assim:

xn+1 =(−296.1)(296887.049)− (−300)(307617.567)

296887.049− 307617.567=

4377014.78

−10730.52= −407.903415

Este e o novo valor que deve ser utilizado de forma a buscar a condicao de contorno desejada.

Executando o Runge Kutta com este valor inicial obtem-se o resultado mostrado abaixo:

dT/dr|r=0.06dT/dr|r→∞

2

-407.903415 -2.688E-09

.

Isto indica que o maximo fluxo de calor que podera ser dissipado por uma aleta deste tipo

e sujeita as condicoes pre-estabelecidas e:

q = −k AbdT

dr

∣∣∣∣r=0.06

= −215 (2π 0.06 0.005) (−407.903415) = 165.31 W

A solucao exata para este problema tambem e obtida em termos de funcao de Bessel e e

dada por:

T (r) = 25 + 31.019K0(9.645 r)

dT

dr= −149.589 (K−1(9.645 r) +K1(9.645 r))

Desta forma pode-se calcular o fluxo de calor maximo transmitido:

dT

dr

∣∣∣∣r=0.06

= −409.23 =⇒ qb = 165.86 W

Desta forma e possıvel verificar que o valor calculado atraves de Runge Kutta, praticamente

coincide com o exato.2O valor para o ∞ foi adotado como sendo 2 m


3.3.1 Procedimento de solucao usando o octave

Os problemas de ordem superior nem sempre fornecem as condicoes necessarias para a

partida do processo de marcha, ou seja, o valor da funcao e de sua derivada no ponto. Em um

grande numero de situacoes serao conhecidos o valor da funcao, ou mesmo da funcao e derivadas,

em dois pontos distintos. Imagine o caso da aleta anterior, se conhece o valor da temperatura

na base e a da temperatura na ponta, ou mesmo a temperatura na base e as condicoes de

conveccao na ponta, o problema nao poderia ser resolvido da maneira apresentada.

Para obter-se a solucao neste caso e necessario o procedimento iterativo, onde admite-se os

valores da funcao e suas derivadas no ponto e busca-se os outros valores nos pontos desejados.

Desta forma imagine o caso resolvido na apostila onde deseja-se saber o fluxo maximo de

calor que pode ser dissipado por uma aleta com estas caracterısticas. Admitindo-se que o

comprimento infinito seja uma aleta de 2 m, busca-se o valor da derivada na base que obtenha

fluxo de calor nulo nesta ponta.

Calcula-se os valores inicialmente para duas ”estimativas”da derivada na base da aleta,

usando-se a seguinte sequencia de comandos:

octave> rr=0.06:0.01:2;octave> fx1=-296.1;octave> fx2=-300;octave> saida=lsode("sol_ale1", [50 fx1], rr);octave> y1=saida(end,2)y1 = 1.1202e+09octave> saida=lsode("sol_ale1", [50 fx2], rr);octave> y2=saida(end,2)y2 = 1.0816e+09

onde os valores fx1 e fx2 sao os valores da derivada da temperatura no ponto r = 2 m.

A partir destes valores utiliza-se o metodo da reta tangente mostrado no texto da apostila

para encontrar a nova (neste caso, terceira) estimativa para o dT/dx|r=6 cm.

octave> fx3=(y2*fx1-y1*fx2)/(y2-y1)fx3 = -409.26

Esta e a nova estimativa da derivada na base da aleta, com a qual podem ser reefetuados

os calculos sucessivas vezes ate que seja encontrado um valor para a derivada da temperatura

nula na posicao desejada. Seguindo o procedimento abaixo:


octave> fx1=fx2;octave> fx2=fx3;octave> saida=lsode("sol_ale1", [50 fx2], rr);octave> y2=saida(end,2)y2 = -5744.3octave> fx3=(y2*fx1-y1*fx2)/(y2-y1)fx3 = -409.25

que resulta num valor para a derivada na base de −409.2547 K/m, e que pode tranquilamente

ser convertido para o valor do fluxo de calor na base q = 164.23 W. O comportamento do

grafico da derivada da temperatura ao longo da posicao pode ser visto utilizando um script de

plotagem similar ao proposto anteriormente na pagina (36).

Um procedimento similar a este, embora utilizando-se do comando pre-implementado no

GNU-Octave, o fsolve (equivalente ao fzero no Matlab), para encontrar zero de funcoes pode

tambem ser utilizado. Para tanto sera necessaria a definicao de uma funcao que se anule quando

o resultado esperado de derivada nula em x = 2 m ocorra. A implementacao foi realizada na

funcao ale inf.m, onde o parametro de entrada x e o vetor com o chute inicial da derivada e

o valor retornado e a derivada da temperatura no ponto afastado da base, como mostrado a

seguir:

function rr=ale_inf(x)raio=0.06:0.01:2;y=lsode("sol_ale1",[50 x], raio);rr=y(length(y),2);endfunction

Definida esta funcao a mesma pode ser facilmente solucionada buscando-se o valor de x que

zera a funcao acima:

octave:65> solu=fsolve("ale_inf", -300)solu = -409.25octave:66> printf("%15.8f\n",solu);-409.25472132

obtendo-se com isto um resultado similar ao anterior e o grafico dos valores de tempertura e

sua derivada tambem podem ser montados a partir do valor correto encontrado em solu.

CAPITULO 4

Solucao de Sistema de Equacoes Lineares

Diversos metodos se apresentam como eficientes na solucao de sistemas de equacoes

lineares. Destacam-se dois grandes grupos que, neste caso, sao os metodos iterativos de solucao

e os metodos diretos de solucao. Neste capıtulo apresentaremos os dois grupos com dois dos

metodos mais usuais, sendo um representante de cada grupo.

4.1 Metodos iterativos de solucao

O grupo dos metodos iterativos de solucao tem como caracterıstica comum o uso de

procedimentos que vao sucessivamente se repetindo, cada vez com novos argumentos, ate a

obtencao da solucao final (dentro de uma certa tolerancia).

Sao representantes destes metodos: Metodo da Sucessiva Sobre-Relaxacao, Metodo de

Gauss-Jacobi, Metodo de Gauss-Siedel, etc. O mais utilizado dentre todos e o metodo de

Gauss-Siedel, cujo o procedimento pode ser descrito em alguns passos:

1. isolar em todo o sistema de equacao, em cada linha a variavel correspondente a diagonal

da matriz (elemento central).

2. admitir uma solucao de partida, em funcao da qual serao calculadas as solucoes posteri-

42


ores.

3. substituir em cada uma das equacoes o valor atual para cada uma das variaveis obtendo,

desta forma, um novo valor para cada variavel.

4. repetir o passo 3 ate que a oscilacao da solucao atinja um valor inferior a tolerancia

adotada.

5. adotar a solucao obtida como solucao do sistema de equacoes.

Tanto o metodo de Gauss-Jacobi como o metodo da sucessiva sobre-relaxacao (S.O.R.) tem

procedimento similar ao descrito. No entanto convem ressaltar que o metodo de Gauss-Jacobi

tem uma menor velocidade de convergencia para a solucao do sistema pois atualiza os valores das

variaveis somente ao final de cada iteracao. O metodo S.O.R. tambem e um metodo eficiente,

no entanto o seu maior complicador e o calculo do valor da sobre-relaxacao, que depende das

condicoes dos sistema linear. Na maior parte dos casos a sua convergencia e ate mais rapida

que a de Gauss-Seidel. Maiores detalhes a respeito destes metodos podem ser encontrados em

(Ruggiero & Lopes, 1988) e (Maliska, 1995).

Exemplo: Resolva o sistema linear mostrado abaixo:

4x1 + 3x2 = 10

2x1 + 4x2 + 3x3 = 19

3x2 + 5x3 + 2x4 = 29

2x3 + 8x4 + 2x5 = 48

3x4 + 5x5 + 4x6 = 61

x5 + 4x6 + 2x7 = 43

2x6 + 3x7 + x8 = 41

2x7 + 6x8 = 62

Solucao:


Se isolarmos o elemento central de cada no obtem-se:

x1 = 14(−3x2 + 10)

x2 = 14(−3x3 − 2x1 + 19)

x3 = 15(−2x4 − 3x2 + 29)

x4 = 18(−2x5 − 2x3 + 48)

x5 = 15(−4x6 − 3x4 + 61)

x6 = 14(−2x7 − x5 + 43)

x7 = 13(−x8 − 2x6 + 41)

x8 = 16(−2x7 + 62)

Partindo de um chute inicial onde xi = 0 para 1 < i < 8 tem-se, aplicando os conceitos de

Gauss-Seidel, que:

x1 = 14(10) = 2, 5

x2 = 14(−2 2.5 + 19) = 3, 5

x3 = 15(−3 3, 5 + 29) = 3, 7

x4 = 18(−2 3, 7 + 48) = 5, 075

x5 = 15(−3 5, 075 + 61) = 9, 155

x6 = 14(−9, 155 + 43) = 8, 461

x7 = 13(−2 8, 461 + 41) = 8, 026

x8 = 16(−2 8, 026 + 62) = 7, 658

Repetido o procedimento por diversas vezes encontra-se as solucoes mostradas na tabela (4.1),

e a solucao foi buscada ate 30 iteracoes, mas notem que para parar neste ponto foi necessario

uma tolerancia superior a 10−3.

Chute Numero de Iteracoes Solucao

inicial 1 2 3 4 5 10 20 30 Exata

0.000 2.500 -0.125 0.972 0.735 0.990 0.953 0.985 0.995 1

0.000 3.500 2.037 2.353 2.014 2.106 2.055 2.018 2.007 2

0.000 3.700 2.548 3.158 2.866 2.917 2.963 2.987 2.995 3

0.000 5.075 3.074 4.314 4.049 4.009 4.010 4.004 4.002 4

0.000 9.155 3.586 4.939 5.047 5.039 4.996 4.996 4.998 5

0.000 8.461 5.840 5.905 5.944 5.967 6.000 6.002 6.001 6

0.000 8.026 7.220 7.088 7.047 7.027 7.000 6.999 6.999 7

0.000 7.658 7.927 7.971 7.984 7.991 8.000 8.000 8.000 8

Tabela 4.1: Resultados obtidos com o metodo de Gauss-Siedel


A1 B1 0 0 · · · 0 0 0 · · · 0 0

C2 A2 B2 0 · · · 0 0 0 · · · 0 0

0 C3 A3 B3 · · · 0 0 0 · · · 0 0

· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·0 0 0 0 · · · Ci Ai Bi · · · 0 0

· · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·0 0 0 0 · · · 0 0 0 · · · Cn An

Figura 4.1: Esquema de uma matriz triangular

4.2 Metodos de inversao direta

Estes metodos diferem dos demais por nao necessitar de procedimento iterativo para

a busca da solucao, eles buscam uma forma de solucionar diretamente o sistema. Dentre

estes metodos destacam-se Metodo de Triangularizacao de Gauss, Metodo de Fatoracao LU,

etc. Existem ainda alguns metodos especıficos para matrizes especiais, o caso tıpico (e muito

utilizado pela peculiaridade de que na matriz cada ponto esta relacionado somente aos seu dois

vizinhos) e o TDMA (TriDiagonal Matrix Algorithm) ou Metodo Linha a Linha. Este metodo

e extremamente util principalmente quando se trabalha com Diferencas Finitas, como veremos

posteriormente, que geram matrizes como estas.

O metodo TDMA se baseia na triangularizacao de Gauss, mas feita para uma matriz Tri-

diagonal. A figura (4.1) mostra a forma geral deste tipo de matriz.

O procedimento para utilizacao do TDMA e um pouco mais complexo que o de Gauss Siedel,

imaginemos que uma linha qualquer do sistema pode ser expressa por:

Ai xi +Bi xi+1 + Ci xi−1 = Di

Repare que para obedecer as condicoes do sistema e necessario que C1 e Bn sejam nulos

resultando em um sistema composto por n incognitas expresso atraves de uma matriz quadrada.

O sistema pode, com isto, ser resolvido atraves da expressao:

xi = Pi xi+1 +Qi

onde:

Pi =−Bi

Ai + CiPi−1

Qi =Di − CiQi−1

Ai + Ci Pi−1


De acordo com as condicoes ja descritas C1 = 0 logo:

P1 =−B1

A1

Q1 =D1

A1

e a partir deles e possıvel calcular todos os outros ate os valores Pn e Qn.

No entanto para obter o campo de temperaturas e preciso comecar pela ultima linha e ir

regredindo ate a primeira. Sabe-se que Bn = 0 e portanto Pn = 0. Desta forma:

xn = Qn

donde e possıvel calcular todos os demais valores para a solucao. Uma descricao mais detalhada

do metodo e sua utilizacao pode ser encontrada em (Maliska, 1995).

Exemplo: Resolva o sistema linear, identico ao solucionado por Gauss-Seidel, mostrado

abaixo:4x1 + 3x2 = 10

2x1 + 4x2 + 3x3 = 19

3x2 + 5x3 + 2x4 = 29

2x3 + 8x4 + 2x5 = 48

3x4 + 5x5 + 4x6 = 61

x5 + 4x6 + 2x7 = 43

2x6 + 3x7 + x8 = 41

2x7 + 6x8 = 62

Solucao:

Este mesmo sistema pode ser expresso matricialmente na forma:

4 3 0 0 0 0 0 0

2 4 3 0 0 0 0 0

0 3 5 2 0 0 0 0

0 0 2 8 2 0 0 0

0 0 0 3 5 4 0 0

0 0 0 0 1 4 2 0

0 0 0 0 0 2 3 1

0 0 0 0 0 0 2 6

x1

x2

x3

x4

x5

x6

x7

x8

=

10

19

29

48

61

43

41

62

e repare na sua formacao Tridiagonal.


Os valores para A, B, C e D sao retirados diretamente de cada linha do sistema de equacoes

e podem ser vistos na tabela (4.2), assim sendo os valores de P1 e Q1 sao:

P1 =−B1

A1

= −3

4= −0, 75 Q1 =

D1

A1

=10

4= 2, 5

E os proximos valores de podem ser calculados atraves da expressao:

P2 =−B2

A2 + C2P1

=−3

4 + 2 (−0, 75)= −1, 2

Q2 =Di − CiQi−1

Ai + Ci Pi−1

=19− 2 2, 5

4 + 2 (−0, 75)= 5, 6

sendo todos os demais P e Q calculados da mesma forma. O resultado para os mesmos neste

caso pode ser encontrado na tabela (4.2).

Calculados todos os valores de P e Q e possıvel voltar calculando a solucao para o sistema

de equacoes. Para a ultima linha sabe-se que:

x8 = Q8 = 8

Com este resultado pode-se calcular o valor de x7 atraves da expressao:

x7 = P7 x8 +Q7 = −0, 607 8 + 11.857 = 7

e assim sucessivamente ate obter a solucao para todas as variaveis. A solucao completa esta a

disposicao na tabela (4.2), exceto pela solucao exata que e identica a obtida pelo TDMA.

Linha A B C D P Q x

1 4 3 0 10 -0.750 2.500 1

2 4 3 2 19 -1.200 5.600 2

3 5 2 3 29 -1.429 8.714 3

4 8 2 2 48 -0.389 5.944 4

5 5 4 3 61 -1.043 11.261 5

6 4 2 1 43 -0.676 10.735 6

7 3 1 2 41 -0.607 11.857 7

8 6 0 2 62 -0.000 8.000 8

Tabela 4.2: Resultados obtidos atraves do metodo TDMA

CAPITULO 5

Equacoes Diferenciais Parciais - Caracterizacao

e Metodos de Discretizacao e Solucao

Ate agora foi visto metodos de solucao de equacoes diferenciais ordinarias, a partir de

agora comecaremos a ver as formas de solucao numerica das equacoes diferenciais parciais. Na

introducao desta apostila ja foram discutidos alguns aspectos deste tipo de solucao e a partir

de agora sera feito um aprofundamento nestes aspectos com a introducao de alguns conceitos

novos.

5.1 Subdivisao das equacoes diferenciais parciais

Da mesma forma que as equacoes diferenciais ordinarias, as E.D. parciais tambem po-

dem ser subdivididas em subgrupos. Embora estes grupos sejam de interesse eminentemente

matematico e importante conhecer as suas caracterısticas principais e assim poder identificar

as caracterısticas de suas solucoes. No caso de problemas fısicos nao se deve preocupar-se tanto

com a natureza de uma equacao diferencial como um todo, que e calculada a partir dos valores

dos coeficientes dos termos da equacao diferencial, mas sim como se comportam cada uma das

coordenadas da mesma.

48


Ponto P -

---

c -

� c

c

(c)

(b)

(a)

Figura 5.1: Figuras ilustrativas dos tipos de equacoes diferenciais

Quanto a classificacao pode-se dizer que as equacoes diferenciais se dividem em:

Equacao Diferencial Parabolica: e um tipo de equacao diferencial em que existe uma

direcao definida para a a propagacao de uma dada perturbacao, vide figura (5.1a). Este tipo

de equacao e a mais simples de ser solucionada numericamente pois a solucao de cada passo de-

pende basicamente do conhecimento do passo anterior. Um bom exemplo da parcela parabolica

de uma equacao diferencial e a coordenada que envolve o tempo, cujo o comportamento e pu-

ramente parabolico. Veja por exemplo que para se conhecer a situacao em um determinado

tempo basta conhecer como este se encontrava no instante anterior. Embora este seja o melhor

exemplo de comportamento deste tipo existem uma serie de outros problemas que podem ter

parcelas deste tipo (equacoes de camada limite, por exemplo).

Equacao Diferencial Elıptica e a mais comuns das formas das equacoes diferenciais, e

neste tipo cada ponto esta ligado a todo o domınio e so se pode resolver o problema a partir

do conhecimento da solucao de todo o domınio, vide figura (5.1b). Todos os problemas que

necessitam de condicoes de contorno em dois pontos recaem em situacoes deste tipo, assim o

grande exemplo de equacao diferencial elıptica sao as coordenadas envolvendo o espaco.

Equacao Diferencial Hiperbolica e aquela em que a perturbacao tambem tem uma direcao

preferencial, no entanto esta direcao preferencial nao e conhecida, depende dos conhecimento,

muitas vezes da propria solucao do problema, vide figura (5.1c). Desta a solucao deste tipo

e praticamente um misto das outras, quando a direcao preferencial e conhecida ela pode ser


resolvida como uma equacao parabolica, quando nao ela e resolvida como uma elıptica. Exem-

plos de equacoes deste tipo sao a parcela convectiva das equacoes de transporte, equacoes de

onda, etc.

5.2 Formas de solucao de equacoes diferenciais parciais

As formas basicas de solucao de equacoes diferenciais parciais ja foram discutidas na

introducao. A respeito do seu uso nao existe propriamente um consenso e varias discussoes se

apresentam. Embora nao deva ficar apenas nisto vou transcrever aqui um trecho da referencia

(Maliska, 1995), cujo o autor e um dos mais respeitados da area numerica de todo o paıs. E

importante notar que se trata apenas da opiniao do autor e nao propriamente de um consenso.

”Com o grande desenvolvimento experimentado pelos metodos numericos e a consequente

penetracao dos mesmos na engenharia, nao raramente se travam discussoes acaloradas a res-

peito da eficiencia do metodo das diferencas finitas (MDF) e elementos finitos (MEF). Minha

intervencao neste ponto polemico deve-se ao fato de ter observado, ao longo dos ultimos 10 anos,

que muitas afirmacoes acerca desses metodos sao oriundas do desconhecimento desta natureza.

Um breve historico e importante para o entendimento. 0 MDF sempre foi empregado pelos

especialistas da area de escoamento de fluidos, enquanto o MEF o foi para area estrutural, na

solucao de problemas da elasticidade. Os problemas, do ponto de vista fısico, sao completa-

mente diferentes. Os de escoamento sao altamente nao-lineares (equacoes de Navier-Stokes),

enquanto os de elasticidade nao possuem os termos convectivos, nao-lineares, e assemelham-se

a problemas puramente difusivos de transferencia de calor. Foi natural, portanto, o fato de os

pesquisadores do MDF terem se concentrado na tentativa de dominar as nao-linearidades dos

termos convectivos e no problema do difıcil acoplamento entre as equacoes, dificuldades nao

encontradas em problemas de elasticidade. Por muito tempo foi deixado para segundo plano o

problema do tratamento de geometrias complexas, e o MDF teve todo o seu desenvolvimento

baseado nos sistemas coordenados ortogonais, como o cartesiano, o cilındrico e o esferico. Por

esta razao, muitas pessoas ainda vinculam o MDF com malhas cartesianas, equivocadamente,

uma vez que ele pode ser aplicado a qualquer tipo de malha, mesmo a nao-estruturada usada

em elementos finitos.

Por outro lado, o MEF sempre teve a vantagem de usar malhas nao-estruturadas, o que

permite que problemas em geometrias complexas possam ser resolvidos. 0 MEF nao teve

penetracao forte na area de fluidos por muito tempo, porque se acreditava que a equacao

diferencial a ser resolvida necessitava de um princıpio variacional para que o metodo pudesse

ser aplicado. Como a equacao de Navier-Stokes nao tem esta propriedade, a aplicacao do MEF


em fluidos foi retardada.

Ate o inicio da decada de 70, tinha-se, portanto, o MDF com grande experiencia na area de

fluidos, mas sem habilidades para tratar geometrias complexas; e o MEF, habil no tratamento

da geometria, mas sem ferramentas para tratar os termos convectivos presentes nas equacoes do

movimento. Mesmo suplantando a questao do princıpio variacional, atraves do uso do metodo

de Galerkin e outras variantes, o MEF nao teve sucesso imediato em problemas de fluidos,

uma vez que o metodo de Galerkin (que e equivalente ao uso de diferencas centrais no MDF)

e adequado apenas para problemas puramente difusivos.

0 uso do metodo de Galerkin em elementos finitos e equivalente ao uso de diferencas cen-

trais em diferencas finitas, ambos produzindo instabilidade em problemas de conveccao domi-

nante. Este e outros problemas similares, que possuem a adequada interpretacao fısica pelo

nao-funcionamento, motivaram pesquisas para o aprimoramento do metodo dos volumes finitos

(MVF), no qual as equacoes aproximadas sao obtidas atraves de balancos de conservacao da

propriedade envolvida (massa, quantidade de movimento, entalpia, etc.) no volume elementar.

A observacao do carater fısico de cada termo da equacao diferencial permitiu que metodos

mais robustos fossem desenvolvidos. A possibilidade de associar a interpretacao fısica com a

matematica influiu de modo consideravel para que praticamente todos os analistas envolvidos

com o MDF passassem a usar o MVF, visto que ambos, por serem equivalentes para uma serie

de problemas, levam muitas pessoas a confundi-los. Importantes desenvolvimentos foram entao

realizados no MVF, mas ainda em coordenadas ortogonais, principalmente cartesianas.

Uma grande transformacao, motivada pelo aparecimento de equipamentos mais velozes,

processou-se na decada de 70. Em meados dessa decada, os sistemas coordenados ortogonais

convencionais comecaram a ceder espaco para os sistemas coordenados generalizados coinciden-

tes com a fronteira do domınio, e o MVF passou a resolver problemas de fluidos em geometria,

irregulares. Nos ultimos 15 anos, foi espantoso o crescimento experimentado pelo MVF em co-

ordenadas coincidentes com a fronteira. Praticamente todos os grandes pacotes hoje disponıveis

no mercado para a solucao de problemas de escoamento de fluidos com transferencia de calor

empregam coordenadas generalizadas no ambito do MVF.

Paralelamente, o MEF passou a empregar outras funcoes de interpolacao para permitir o

tratamento adequado dos termos convectivos nao-lineares. As funcoes do tipo Petrov-Galerkin,

que nada mais sao do que a ponderacao entre os efeitos difusivos e convectivos, semelhantes

aos esquemas hıbridos empregados em volumes finitos, possibilitaram um expressiva avaliacao

do MEF na area de escoamento de fluidos. Recentes formulacoes, onde estas funcoes sao

desenvolvidas ao longo da linha de corrente, tambem equivalentes aos esquemas skew usados

em volumes finitos, permitiram que o MEF passasse, tambem, a tratar problemas de fluidos


minimizando os efeitos de difusao numerica.

Atualmente, um grande esforco de pesquisa esta sendo dedicado ao desenvolvimento de

metodos em volumes finitos, usando malhas nao-estruturadas, semelhantes, portanto, aquelas

usadas em elementos finitos.

No panorama atual, observa-se que ambos os metodos (MVF e MEF) estao resolvendo

problemas altamente convectivos, inclusive com ondas de choque, em geometrias arbitrarias,

mostrando que existe entre eles uma forte semelhanca em termos de generalidade. Se olhar-

mos do ponto de vista matematico, isto poderia ser diferente, uma vez que todos os metodos

numericos podem ser derivados do metodo dos resıduos ponderados, empregando-se diferentes

funcoes peso. Por exemplo, o MDF surge quando a funcao peso e feita igual a funcao delta no

ponto em consideracao; o MVF aparece quando esta funcao peso e feita igual a 1 no volume

elementar, e a zero em todos os outros volumes elementares, ja o MEF-Galerkin surge quando

estas funcoes peso sao feitas iguais as funcoes tentativas. Logo, nao existe sentido em argu-

mentar que um determinado metodo e sempre superior a outro, visto que eles sao derivados

do mesmo principio e diferem apenas na forma de minimizacao escolhida. 0 que se tem, na

pratica, sao diferentes graus de experiencia dos diversos metodos para diferentes problemas.

A preferencia pessoal deste autor pelo metodo dos volumes finitos (MVF) para problemas

de escoamento de fluidos e justificada primeiro pela escola seguida na sua formacao e, segundo,

pelo fato de o MVF, ao criar suas equacoes aproximadas, estar realizando um balanco da propri-

edade em nıvel de volumes elementares. Se o que se busca com o metodo numerico e a solucao

da equacao diferencial que representa a conservacao da propriedade em inves de ponto (infinite-

simal), parece logico que as equacoes aproximadas (que formam o sistema linear) representem

a conservacao em nıvel de volumes elementares (discreto). A depuracao de um programa com-

putacional tambem fica mais facil quando o analista tem etapas a serem conferidas. Como no

MVF os balanco de conservacao devem ser satisfeitos em nıvel de volumes elementares, para

qualquer tamanho de malha, todos os princıpios de conservacao podem ser checados em uma

malha bastante grosseira. Ou seja, quase tudo pode ser feito manuseando-se poucos resultados

em execucoes rapidas no computador. Em outros metodos, pode-se apenas conferir a solucao

com uma malha refinada.

Recentes desenvolvimentos mostram tambem o MEF aplicado em nıvel de volumes ele-

mentares, sendo denominado metodo dos elementos finitos baseado no volume de controle,

conhecido na literatura internacional como CVFEM - Control Volume Finite Element Method,

cujo objetivo e obter as equacoes aproximadas em nıvel de volumes elementares em uma base

de elementos finitos. Muitos autores, principalmente aqueles ligados ao MEF classico, nao con-

sideram o CVFEM como um MEF. Entretanto, foge do nosso escopo aprofundar esta e outras


questoes especıficas.

Um outro metodo que vem ganhando destaque e espaco e o metodo dos elementos no

contorno (Boundary Element Method - BEM da literatura internacional). Sua vantagem e a

possibilidade de tratar apenas com a discretizacao da fronteira, sem necessidade de discretizar

o domınio interno. 0 metodo e aplicado quando e possıvel transferir a influencia do operador do

domınio para a fronteira. Apesar de atraente, e um metodo que ainda esta longe de responder

as solicitacoes dos problemas complexos resolvidos pelos outros metodos. Sem duvida e uma

area de pesquisa que merece esforcos.”

Depois desta discussao toda o metodo que sera estudado mais a fundo e o metodo de

diferencas finitas (MDF) por se tratar de um metodo tradicional, ainda amplamente utilizado

e, que segundo o proprio texto mostra, nao difere muito dos demais.

5.3 Discretizacao do domınio

Como ja foi discutido a solucao numerica nao gera solucoes contınuas como as analıticas

mas sim valores para pontos determinados. Trata-se portanto de uma etapa fundamental na

solucao numerica de um problema a escolha dos pontos para os quais vai se obter solucao, o

que e determinado pela discretizacao do domınio. Muitas vezes se obtem valores para pontos

fora dos pontos da malha, mas estes valores sao obtidos por interpolacao da solucao original, o

que nao garante a correcao do seu valor, embora este procedimento forneca bons resultados na

maior parte dos casos.

Um dos grandes exemplos onde este procedimento falha e quando se tem um problema

de solidificacao ou fusao. Entre o ponto um ponto da malha que e solido e outro lıquido ha

uma forte descontinuidade, e uma simples interpolacao linear neste caso nao fornecera bons

resultados, embora a solucao para cada ponto esteja correta.

A escolha da malha tambem pode influir decisivamente na solucao de um problema numerico.

Um grande exemplo disto e que a utilizacao de dois pontos muito proximos numa malha onde a

distancia media entre os pontos e significativamente maior pode induzir em erros significativos

na solucao final do problema. Isto e explicavel pela inclusao de uma quase-singularidade na

matriz (dois pontos que respondem pela mesma equacao), dificultando a solucao do sistema.

Recentemente tem sido desenvolvidas muitas tecnicas de geracao automatica de malhas, em

(Maliska, 1995) pode ser encontrada descricao de uma serie destes metodos. Este procedimento

e um conforto necessario pois, a medida que as geometrias dos problemas vao se tornando mais


complexas, o tempo necessario para gerar uma malha manualmente aumenta numa proporcao

muito maior. E importante ressaltar que para se trabalhar com a formulacao de diferencas a

ser apresentada nesta apostila e necessaria que a malha seja ortogonal1.

Existem formulacoes alternativas para diferencas finitas e volumes finitos que trabalham

com malhas ortogonais. O metodo de elementos finitos nao requer malhas ortogonais.

1Malha ortogonal e o tipo de malha em que cada ponto tem seus vizinhos diretos numa direcao normal outangencial.

CAPITULO 6

Metodo dos Elementos Finitos

Sera mostrado neste capıtulo todos os princıpios da formulacao de um problema atraves

do metodo dos elementos finitos. Serao apresentados, de forma sucinta, os princıpios da meto-

dologia de elementos finitos, todos os passos a serem cumpridos na formulacao de um problema

desde a discretizacao ate solucao final do problema.

6.1 Princıpios gerais

6.1.1 Aproximacao por funcoes

O princıpio fundamental de elementos finitos consiste em utilizar funcoes, de diferentes

ordem, para aproximar a solucao dentro do domınio do elemento (subdomınio do problema).

Normalmente quando se tem uma determinada quantidade expressa em termos de uma variacao,

espacial ou temporal, pode-se determinar uma equacao de aproximacao. A qualidade desta

aproximacao pode ser determinada atraves de um desvio expresso por:

e(x) = u(x)− uex(x)

sendo e(x) o desvio, u(x) a aproximacao e uex o valor exato desta aproximacao para este mesmo

ponto.

55


Para construir uma aproximacao exata e necessario escrever a aproximacao de u(x) em

funcao da posicao e dos parametros de aproximacao:

• escrever uma funcao contendo n parametros ai como, por exemplo u(x) = a1 + a2 x +

a3x2 + · · ·+ an x

n − 1, ou qualquer outra funcao onde:

u = f(x, a1, a2, a3, · · · an)

• determinar os parametros ai da aproximacao usando uma determinada regra. Esta regra

pode ser a de mınimos quadrados, ou ainda, a mais comum para determinacao das funcoes

de elementos finitos, a funcao em que o desvio e(x) = 0, nos pontos considerados.

Desta forma obtem-se uma funcao que simplesmente aproxima o comportamento da variavel

ao longo do domınio. Com esta aproximacao podem ser obtidas, dentre outras coisas,:

• um expressao simples para uma funcao complexa ou difıcil de ser manipulada, valida para

um certo numero de pontos que se deseje ou ainda com boa aproximacao dentro de uma

certa regiao;

• solucao para equacoes diferenciais parciais ou ordinarias (normalmente associadas a pro-

blemas fısicos).

Exemplo Aproximacao da quantidade fısica: Suponha que deseja-se obter a distribuicao de

temperatura numa dada regiao e sabe-se que esta e uma funcao contınua (sem ressaltos) e sao

conhecidos os seus valores em tres pontos:

x uex(x)

0,0 20◦C

0,5 25◦C

1,0 22◦C

Pode-se determinar uma expressao para as temperaturas entre os pontos 0 e 1, utilizando a

o desvio nulo (e(x) = 0) nos pontos conhecidos e uma aproximacao polinomial quadratica:

uex(x) ≈ u(x) = a1 + a2 x+ a3 x2


assim:

uex(x = 0) = u(x = 0) = a1 = 20

uex(x = 0, 5) = u(x = 0, 5) a1 + 0, 5 a2 + (0, 5)2 a3 = 25

uex(x = 1) = u(x = 1) a1 + 1 a2 + (1)2 a3 = 22

Se resolvido o sistema com tres equacoes e tres incognitas obtem-se que:

a1 = 20; a2 = 18; a3 = −16;

que resulta em u(x) = 20+18x−16x2 e representa uma aproximacao do uex(x). Por exemplo,

para uma aproximacao da temperatura na posicao x = 0, 7 pode ser usado:

u(x = 0, 7) = 20 + 18× 0, 7− 16× (0, 7)2 = 20 + 12.6− 7, 84 = 24.76

Exemplo: Solucao aproximada de uma equacao diferencial: Supondo que se tem um pro-

blema de conducao com a geracao de energia variando com a posicao. Neste caso, a equacao

diferencial que rege o problema seria do tipo:

d2u

dx2= −f(x) para 0 ≤ x ≤ 1

e sabendo que:

x f(x)

0,25 1

0,75 0,25

Tem-se, ainda, que as temperaturas nas extremidades do problema seriam conhecidas. Neste

caso, as condicoes de contorno adotadas sao:

? x = 0 : a temperatura u = 0

? x = 1 : a temperatura u = 0

Para a obtencao da solucao tambem se faz necessario a escolha de uma funcao que obedeca

as condicoes de contorno. Assim, foi escolhida a seguinte aproximacao:

uex(x) ≈ u(x) = a1 sin(π x) + a2 sin(2π x)


Normalmente em solucoes deste tipo e mais difıcil encontrar uma solucao que satisfaca as

condicoes de contorno do que partir dela, e encontrar a solucao do problema, propriamente

dita. Desta forma tem-se que para este caso, nos dois pontos considerados:

d2u

dx2

∣∣∣∣x=xi

= −a1 π2 sin(π xi)− 4a2 π

2 sin(2π xi)

e para os pontos analisados:

d2u

dx2

∣∣∣∣x=0,25

= a1 π2 sin(0.25π)− 4a2 π

2 sin(0.5π) = −1

d2u

dx2

∣∣∣∣x=0,75

= a1 π2 sin(0.75π)− 4a2 π

2 sin(1.5π) = −0.25

Resolvendo o sistema composto por duas equacoes e duas incognitas, obtem-se:

a1 =5

4√

2

1

π2; a2 =

3

32

1

π2;

que implica numa solucao:

u(x) =5

4√

2

1

π2sin(π x) +

3

32

1

π2sin(2π x)

Pode-se utilizar ainda a expressao para avaliar a temperatura, ou outra variavel de interesse,

em qualquer ponto. Por exemplo para a posicao x = 0, 25, tem-se:

u(x = 0, 25) =23

32

1

π2= 0, 0728

6.1.2 Aproximacao nodal

Este mesmo procedimento pode ser utilizado para expressar uma funcao para uma

variavel generica em funcao dos pontos nodais. Para isto, note primeiro que nos casos an-

teriores poderia se fazer uma aproximacao em funcao das variaveis a1, · · · , an, bastando para

isto:

u(x) = P1(x) a1 + P2(x) a2 + · · ·+ Pn(x) an

sendo Pi chamada de funcao basica de interpolacao (basis function).

Imagine que se tenha um domınio qualquer e se conheca os valores de uma funcao qualquer

em pontos determinados. Pode-se montar uma funcao que represente o comportamento da

variavel em funcao do seu valor nestes pontos. Para uma funcao u(x) tem-se que para uma


serie de pontos nodais x1, x2, · · · , xn os valores da funcao sao, respectivamente, u1, u2, · · · , un.

Pode-se desta forma fazer a aproximacao nodal:

u(x) = N1(x)u1 +N2(x)u2 + · · ·+Nn(x)un

ou ainda numa forma matricial:

u(x) = [N1(x)N2(x) · · · Nn(x)] =

u1

u2

...

un

= [N ] · u

sendo Ni chamada de funcao de interpolacao (interpolation function)

Assim sendo, e possıvel notar que

(a) como u(xi) = ui, as funcoes de interpolacao assumem os valores::

Nj(ui) =

{0 sei 6= j

1 sei = j

(b) o erro da aproximacao nos pontos nodais e nulo.

e(x) = 0 se x = xi

Exemplo Aproximacao para quatro pontos: Considere a situacao onde sao escolhidos quatro

pontos alinhados para solucionar o problema unidimensional. Desta forma a aproximacao nodal

e dada por:

u(x) = N1(x)u1 +N2(x)u2 +N3(x)u3 +N4(x)u4

Neste caso e adotada como funcao de interpolacao que respeite as condicoes anteriormente

expostas: a funcao de Lagrange de terceiro grau,

Ni(x) =4∏

j=1j 6=i

(x− xj)

(xi − xj)

(6.1)

Tomando o exemplo de N1(x) tem-se:

N1(x) =(x− x2)(x− x3)(x− x4)

(x1 − x2)(x1 − x3)(x1 − x4)


1 2 3 4 5 6 7-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

N1(

x)x

Figura 6.1: Comportamento da funcao N1(x) na regiao considerada

Assim se x1 = 1, x2 = 2, x3 = 5 e x4 = 7 tem-se que:

N1(x) = − 1

24(x− 2) (x− 5) (x− 7)

que se substituıdo

x 1 1,5 2 3 4 5 6 7

N1(x) 1 − 77192

0 −13

−14

0 16

0

O comportamento desta funcao pode ser visto para todos os pontos na figura (6.1).

Raciocınio analogo e usado para as demais funcoes:

N2(x) =(x− x1)(x− x3)(x− x4)

(x2 − x1)(x2 − x3)(x2 − x4)=

1

15(x− 1) (x− 5) (x− 7)

N3(x) =(x− x1)(x− x2)(x− x4)

(x3 − x1)(x3 − x2)(x3 − x4)= − 1

24(x− 1) (x− 2) (x− 7)

N4(x) =(x− x1)(x− x2)(x− x3)

(x4 − x1)(x4 − x2)(x4 − x3)=

1

60(x− 1) (x− 2) (x− 5)

Assim, considerando a interpolacao da propria variavel espacial, ou seja u1 = 1, u2 = 2,

u3 = 5 e u4 = 7: tem se que a funcao de interpolacao u(x) e dada por:


u(x) = N1(x) u1 +N2(x) u2 +N3(x) u3 +N4(x) u4

= − 1

24(x− 2) (x− 5) (x− 7) +

2

15(x− 1) (x− 5) (x− 7)−

5

24(x− 1) (x− 2) (x− 7) +

7

60(x− 1) (x− 2) (x− 5)

= x

que resulta com isto na solucao exata: u(x) = x.

Este e apenas um exemplo unidimensional, no entanto na maior parte dos problemas que sao

tratados em elementos finitos, tem-se duas ou mais dimensoes, e as novas funcoes de interpolacao

passam a ser funcao de todas as variaveis espaciais:

uex(x, y, z) = uex(~x)

ou ainda

u(~x) = [N1(~x)N2(~x) · · · Nn(~x)] =

u1

u2

...

un

= [N ] · {u}

6.2 Aproximacao por Elementos Finitos

A construcao das funcoes de aproximacao u(x) vao se tornando mais difıceis a medida

que o numero de nos aumenta. Uma maior complexidade ainda e verificada quando o domınio

V e irregular ou possui condicoes de contorno mais complexas. Por outro lado, as aproximacoes

nodais de sub-domınio simplificam a obtencao da solucao u(x) e sao extremamente faceis de

serem implementadas em um computador. Este procedimento consiste, basicamente, dos se-

guintes passos:

• divisao do domınio principal V em sub-domınios V e;

• a escolha de uma aproximacao nodal adequada para cada subdomınio. De maneira geral,

ela depende das dos pontos nodais e aproximacoes utilizadas nas vizinhancas. O metodo

dos elementos finitos e apenas um tipo de aproximacao nodal por sub-domınio, sendo

suas caracterısticas principais:

– a aproximacao nodal dentro do sub-domınio depende apenas dos nos do proprio

elemento;


– a aproximacao elementar ue(x) deve garantir um mınimo de continuidade entre dois

elementos assim como nas suas fronteiras.

Algumas definicoes importantes:

Nos: sao os pontos do sub-domınio onde a funcao e avaliada;

Coordenadas nodais: sao as coordenadas geometricas dos pontos avaliados.

Variaveis nodais: sao valores da funcao de interesse, u(x), nos nos.

6.2.1 Definicao geometrica dos elementos.

Imagine o problema generico aproximado por elementos compostos pelos nos 1, 2, 3 e 4 e

compreendendo o domınio total entre x1 e x4. Os tres elementos lineares seriam responsaveis

pelo seguintes domınios:

Elemento 1 V 1 x1 < x < x2



Utilizando para cada um destes elementos de dois nos a funcao Lagrange tem-se que para

cada elemento:

ui(x) = N1 u1 +N2 u2

sendo Ni dado pela expressao de Lagrange, equacao (6.1). Deve-se notar que no exemplo

anteriormente resolvido, utilizava-se um polinomio de Lagrange de quarta ordem, enquanto no

caso de elementos de dois nos utiliza-se apenas o polinomio de segunda ordem.

6.2.2 Regras para a discretizacao de domınios atraves de elementos

finitos

A subdivisao de um domınio V em sub-domınios V e deve obedecer a algumas regras.

(a) deve haver sempre uma fronteira comum entre os elementos adjacentes onde estarao os

unicos pontos comuns entre os elementos. Estas fronteiras podem ser compostas por

pontos, linhas ou areas.


(a) Overlapping (b) Hole

Figura 6.2: Anomalias comuns na discretizacao de um domınio atraves de malhas.

(b) nao e permitida a existencia de regioes comuns a mais de um elemento (overlapping) e nem

regioes dentro do domınio que nao pertencam a regiao alguma (holes). Estas anomalias

estao mostradas na figura (6.2).

(c) quando a fronteira do domınio do problema nao coincidem inclui-se uma anomalia (hole)

que acarreta em um erro que nao pode ser mensurado. Entretanto, estes erros, deno-

minados por erros geometricos, podem ser minimizados utilizando-se elementos menores

ou elementos de maior ordem, que melhor se adequam a fronteira. Estes procedimen-

tos sao melhores compreendidos quando discretiza-se uma fronteira do tipo mostrado na

figura (6.3).

Note que na fig. (6.3a) a discretizacao e feita com elementos grosseiros e nao e possıvel

representar a fronteira de forma adequada, apresentando maiores vaos (holes), nestas

regioes. Na fig. (6.3b) utiliza-se elementos menores e a fronteira ja pode ser melhor repre-

sentada. E, finalmente, na figura (6.3c), nota-se que foi obtida uma boa representacao da

fronteira mesmo com elementos mais grosseiros, desde que estes sejam de ordem superior

e suas fronteiras possam ser deformadas para se ajustar ao domınio do problema.

(a) Malha Linear Grosseira (b) Malha mais refinada (c) Malha com elementos

de maior ordem

Figura 6.3: Erro comum na discretizacao do domınio.


Figura 6.4: Transformacao do domınio elementar.

6.3 Elementos de referencia

Quando se trabalha com elementos finitos existem dois tipos basicos de sistemas de

coordenadas. A metodologia de elementos finitos usa ferramentas matematicas para efetuar a

transformacao (τ) de um sistema de coordenadas para outro. A figura (6.4) mostra o sistemas

de referencia para um elemento triangular. Os sistemas de coordenadas que coexistem neste

caso sao:

Sistema de coordenadas local: e aquele que existe dentro do elemento ideal e cujas coor-

denadas sao expressas em termos de ξ e η.

Sistema de coordenadas global: que existe no sistema fısico real e expresso em termos de

coordenadas globais para todos os elementos.

Pode-se expressar a transformacao de um sistema de coordenadas em outro atraves da

expressao:

τ e : ξ =⇒ xe(ξ)

onde o ponto ξ representa um ponto no sistema de coordenadas local e x o ponto no sistema

global.

Esta mesma expressao pode ser escrita em termos de variaveis nodais sendo expressa na

forma:

τ e : ξ =⇒ xe (ξ) = N(ξ) · {xe}

Para que esta transformacao seja feita de maneira adequada sao necessarias algumas condicoes:


• o sistema de transformacao deve ser biunıvoco (um ponto em cada conjunto tem corres-

pondencia a outro unico ponto do outro sistema).

• cada no local do sistema de coordenadas corresponde a um no do sistema generalizado e

vice-versa.

• da mesma forma que os nos cada fronteira do elemento corresponde a uma fronteira global.

Neste sistema de coordenadas pode-se utilizar as proprias coordenadas para verificar as

funcoes de interpolacao. Desta forma:

x(ξ, η) = N1(ξ, η) · xi +N2(ξ, η) · xj +N3(ξ, η) · xk = [N ] · {x}

y(ξ, η) = N1(ξ, η) · yi +N2(ξ, η) · yj +N3(ξ, η) · yk = [N ] · {y}

sendo N as funcoes de interpolacao.

Estes conceitos de transformacao de coordenadas podem ser vistos mais detalhadamente

quando se realiza uma transformacao para o elemento triangular, originalmente proposta na

figura (6.4), assim:

~P | = ~R + ~r

~r = ξ · iξ + η · iη~R = Xi

iξ = Xj −Xi

iη = Xk −Xi

P = Xi + ξ · (Xj −Xi) + η · (Xk −Xi)

P = (1− ξ − η) ·Xi + ξ ·Xj + η ·Xk

P = {(1− ξ − η) ξ η}

Xi

Xj

Xk

E possıvel notar que para que a transformacao de coordenadas seja adequada e preciso

considerar as amplitudes da transformacao, assim como a multipla composicao de dimensoes

(η, por exemplo, apresenta derivadas tanto na direcao x como y). Para corrigir este fator

e necessario conhecer o Jacobiano [J ] da transformacao. O Jacobiano pode ser dado pela

expressao:

[J ] =

[∂x∂ξ

∂y∂ξ

∂x∂η

∂y∂η

]=

[xj − xi yj − yi

xk − xi yk − yi

]det[J ] = (xj − xi)(yk − yi)− (xk − xi)(yj − yi)


E importante notar que o valor do determinante do Jacobiano resulta em duas vezes a area

do mesmo.

6.4 Propriedades das funcoes de aproximacao u(x)

Propriedade Fundamental: o valor da funcao deve coincidir com os respectivos valores no-

dais. Com base na expressao geral:

u(x) = N1(x)x1 +N2(x)x2 + · · ·Nn(x)xn =n∑

i=1

Ni(x)xi

tem-se que a propriedade e atendida se sobre o ponto nodal considerado j tem-se:

Nj(xi) =

{0 se i 6= j

1 se i = j

Continuidade no elemento: as funcoes Ni(x) devem ser contınuas em todo elemento, assim

como as derivadas ate a ordem s considerada.

Continuidade entre elementos: tanto os valores da funcao como de suas derivadas devem

ser os mesmos nas fronteiras de elementos adjacentes, sejam eles calculados por um ou

outro elemento.

Funcao polinomial completa O erro de truncamento e minimizado com a diminuicao do

elemento. Por muitas outras razoes e necessario diminuir o erro das derivadas e da funcao

de aproximacao. Alem do mais, e necessario que para que se resolva um problema ele

possua pelo menos a mesma ordem (s) de derivadas contınuas na funcao de interpolacao.

Caso contrario a funcao seria automaticamente anulada e nao se conseguiria os resultados

desejados para o problema.

Alem disto existem algumas definicoes importantes que devem ser destacadas:

• se uma funcao e contınua pelas suas fronteiras ele e classificada com C0, se a funcao e a

primeira derivada sao contınuas, ela e classificada de C1, e assim sucessivamente.

• se a transformacao de coordenadas e as funcoes de aproximacao se utilizam das mesmas

funcoes de interpolacao, o elemento e chamado de isoparametrico.

• se as funcoes de interpolacao sao diferentes podem haver dois tipos de transformacao:

pseudo-parametrico: se as funcoes sao diferentes mas utilizam a mesma base.


subparametrica: quando as funcoes de interpolacao da geometria sao de ordem inferior

a da variavel de interesse.

• o numero de variaveis nodais associados a cada um dos nos de cada elemento e denominado

por grau de liberdade do sistema.

6.5 Construcao das funcoes de interpolacao

Escolha da base polinomial: esta diretamente associada ao numero de pontos nodais do

elemento, sendo que quanto maior o numero de pontos, maior a base polinomial. A tabela (6.1)

mostra a base associada ao numero de pontos do elemento.

Tabela 6.1: Base polinomial para elementos de ate duas dimensoes.

Dimensao Grau do Base Polinomial Nos

Polinomio nd

Base completa

1 1 < 1 ξ > (linear) 2

1 2 < 1 ξ ξ2 > (quadratico) 3

2 1 < 1 ξ η > (linear) 3

2 2 < 1 ξ η ξ2 ξ η η2 > (quadratico) 3

Base incompleta

2 2 < 1 ξ ξ η η > (quadratico) 3

Para a construcao do polinomio e utilizada a propriedade basica da funcao de interpolacao,

ou seja, o fato de que no ponto nodal o valor resultante e igual ao proprio valor da funcao.

Sabendo que a transformacao de sistemas e dada na forma:

u(ξ) =< a(x) > (6.2)

sendo a base polinomial e < a > as variaveis generalizadas para montagem da funcao

de interpolacao. A expressao apresentada na eq. (6.2) representa a aproximacao generalizada,

diferenciada da aproximacao nodal.


Relacao entre o sistema generalizado e nodal: e dada pela aproximacao avaliada no

ponto nodal:

{un} =

< P1(ξ1) P2(ξ1) · · · Pn(ξ1) >

< P1(ξ2) P2(ξ2) · · · Pn(ξ2) >

· · · · · ·< P1(ξn) P2(ξn) · · · Pn(ξn) >

{a} = [Pn]{a}

ou ainda invertendo a matriz dos polinomios:

{a} = [Pn]−1{un} (6.3)

Expressoes analıticas para as funcoes de interpolacao: podem ser obtidas a partir dos

resultados anteriores com respeito da base polinomial e o vetor a. Se substituir-se a eq. (6.3)

na eq. (6.2), obtem-se:

u(ξ) = [Pn]−1{un}

que implica em dizer que as funcoes de interpolacao para a aproximacao nodal:

u(ξ) =< N(ξ) > {un} = [Pn]−1{un}

ou ainda as funcoes de aproximacao nodal sao dados por:

< N(ξ) >= [Pn]−1 (6.4)

Derivadas de u(ξ): podem ser obtidas diferenciando em relacao a cada uma das variaveis

geometricas, como no caso unidimensional:{uξ

uη

}=

[< Pξ >

< Pη >

][Pn]−1 {un} =

[< Nξ >

< Nη >

]{un} = [Bξ] {un}

Exemplo: Construcao das funcoes de interpolacao para um elemento quadrilateral de quatro

nos: Serao utilizados agora os preceitos fornecidos para a obtencao das funcoes de interpolacao

de um elemento linear de quatro nos, como mostrado na figura (6.5).

Escolha da base polinomial: ja que o elemento possui quatro nos, e preciso escolher uma

base de acordo com este tamanho. Ja que nao existe nenhuma base completa para quatro

nos a melhor escolha, que respeita as condicoes de de simetria e inter-continuidade dos

elementos e:

〈P 〉 = 〈 1 ξ η ξ · η 〉


Figura 6.5: Elemento bidimensional, linear e isoparametrico.

Montagem da matriz Pn: e montada com base na avaliacao da base polinomial nos quatro

nos possıveis (com seus respectivos valores de ξ e η:

ξ =

−1

1

1

−1

, η =

−1

−1

1

1

[Pn] =

1 −1 −1 1

1 1 −1 −1

1 1 1 1

1 −1 1 −1

Inversao de Pn: neste caso a matriz e inversıvel:

[Pn]−1 =1

4[Pn]T =

1

4

1 1 1 1

−1 1 1 −1

−1 −1 1 1

1 −1 1 −1

Obtencao das expressoes de 〈N〉: a partir da multiplicacao da base polinomial pela inversa

de [Pn]


〈N〉 = 〈P 〉 · [Pn]−1 = { 1 ξ η ξ · η } ·1

4

1 1 1 1

−1 1 1 −1

−1 −1 1 1

1 −1 1 −1

〈N〉 =

⟨N1(ξ, η) N2(ξ, η) N3(ξ, η) N4(ξ, η)

⟩〈N〉 =

1

4

⟨1− ξ − η + ξη 1 + ξ − η − ξη 1 + ξ + η + ξη 1− ξ + η − ξη

⟩〈N〉 =

1

4

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩

O elemento determinado e isoparametrico e portanto as mesmas funcoes utilizadas para a

interpolacao de uma quantidade generica podem ser utilizadas para as variaveis espaciais (x e

y no caso deste elemento bidimensional):

x(ξ, η) = 〈N〉 ·

x1

x2

x3

x4

, y(ξ, η) = 〈N〉 ·

y1

y2

y3

y4

6.6 Transformacao de operadores diferenciais

As equacoes que governam os fenomenos fısicos envolvem normalmente nao somente as

funcoes u, mas tambem as suas derivadas. Como ja foi destacado, a aproximacao no espaco

real e sempre complexa devendo ser dada preferencia para se trabalhar no domınio elementar.

Isto implica na utilizacao das funcoes de aproximacao no espaco ξ:

uex ≈ u(ξ) = 〈< N(ξ)〉{un}

sendo possıvel a obtencao da solucao atraves da transformacao de coordenadas. Normalmente a

transformacao e complexa, como ja foi enfatizado. De qualquer forma, quando e preciso avaliar

uma derivada da referida funcao, seja qual for a transformacao desejada, esta e feita atraves do

princıpio da regra da cadeia, que numa forma matricial poderia ser expressa por:∂∂ξ∂∂η∂∂ζ

=

∂x∂ξ

∂y∂ξ

∂z∂ξ

∂x∂η

∂y∂η

∂z∂η

∂x∂ζ

∂y∂ζ

∂z∂ζ

·

∂∂x∂∂y∂∂z

=⇒ {∂ξ} = [J ] · {∂x}

sendo [J ] o Jacobiano da matriz de transformacao.


De forma analoga pode-se enfatizar que transformacao inversa respeita a mesma regra:∂∂x∂∂y∂∂z

=

∂ξ∂x

∂η∂x

∂ζ∂x

∂ξ∂y

∂η∂y

∂ζ∂y

∂ξ∂z

∂η∂z

∂ζ∂z

·

∂∂ξ∂∂η∂∂ζ

=⇒ {∂x} = [j] · {∂ξ}

e que nos permite afirmar que [j] = [J ]−1.

Com isto e possıvel obter qualquer matriz de transformacao, desde que esta seja biunıvoca

e, portanto, gere matrizes inversıveis.

Alguns Jacobianos e seus inversos

• Unidimensional:

[J ] = J1,1 [j] = [J ]−1 =1

J1,1

• Bi-dimensional:

[J ] =

[J1,1 J1,2

J2,1 J2,2

][j] = [J ]−1 =

1

det[J ]

[J2,2 −J1,2

−J2,1 J1,1

]det[J ] = J1,1 · J2,2 − J2,1 · J1,2

• Tri-dimensional:

[J ] =

J1,1 J1,2 J1,3

J2,1 J2,2 J2,3

J3,1 J3,2 J3,3

[j] = [J ]−1 =1

det[J ]

J2,2J3,3 − J2,3J3,2 −J1,2J3,3 + J1,3J3,2 J1,2J2,3 − J1,3J2,2

−J2,1J3,3 + J2,3J3,1 J1,1J3,3 − J1,3J3,1 −J1,1J2,3 + J1,3J2,1

J2,1J3,2 − J2,2J3,1 −J1,1J3,2 + J1,2J3,1 J1,1J2,2 − J1,2J2,1

det[J ] = J1,1(J2,2J3,3 − J2,3J3,2) + J1,2(J3,1J2,3 − J2,1J3,3) + J1,3(J2,1J3,2 +−J3,1J2,2)

6.6.1 Montagem do Jacobiano

Considerando a matriz de transformacao de variaveis entre um espaco real x e um espaco

elementar ξ tem-se que a aproximacao nodal se apresenta na forma:

{ x y z } = 〈N〉[ {x}n {yn} {zn} ]


E o Jacobiano pode ser montado a partir das proprias derivadas da funcao:

[J ] =

∂∂ξ∂∂η∂∂ζ

{x y z

}=

∂〈N〉∂ξ

∂〈N〉∂η

∂〈N〉∂ζ

[ {xn} {yn} {zn}]

Exemplo: Jacobiano para o elemento de quatro nos:

[J ] =1

4

[−(1− η) (1− η) (1 + η) −(1 + η)

−(1− ξ) −(1 + ξ) (1 + ξ) (1− ξ)

]×

x1 y1

x2 y2

x3 y3

x4 y4

que apos realizada a multiplicacao resulta em:

[J ] =1

4

−x1 + x2 + x3 − x4 + η (x1 − x2 + x3 − x4) −y1 + y2 + y3 − y4 + η (y1 − y2 + y3 − y4)

−x1 − x2 + x3 + x4 + ξ (x1 − x2 + x3 − x4) −y1 − y2 + y3 + y4 + ξ (y1 − y2 + y3 − y4)

sendo que o determinante da matriz pode ser calculado atraves da expressao:

det[J ] = A0 + A1ξ + A2η

e onde:

A0 =1

8[(y4 − y2)(x3 − x1)− (y3 − y1)(x4 − x2)]

A1 =1

8[(y4 − y3) (x1 − x2)− (y2 − y1)(x3 − x4)]

A2 =1

8[(y1 − y4)(x2 − x3)− (y2 − y3)(x1 − x4)]

6.6.2 Transformacao de uma integral

Supondo agora que se deseje realizar uma integracao de uma funcao generica f no

domınio elementar, muito mais simples, e transforma-la para o domınio real. Para realizar este

objetivo e preciso analisar a natureza da integracao. Para uma integracao sobre um determinado

volume no espaco real tem-se que:

dV = (dx× dy) · dz

sendo

d~x = J1,1 · d~ξ + J2,1 · d~η + J3,1 · d~ζ

d~y = J1,2 · d~ξ + J2,2 · d~η + J3,2 · d~ζ

d~z = J1,3 · d~ξ + J2,3 · d~η + J3,3 · d~ζ


que implica que esta mesma funcao dV pode ser representada na forma:

dV = det[J ] · dξ · dη · dζ

6.7 Coordenadas nodais e conectividade

Com os nos numerados de forma sequencial de 1 ate n e expresso num sistema de

coordenadas global (generalizado), pode-se montar um descricao com todos os pontos em uma

tabela. No caso bidimensional devem constar na tabela o numero do no e sua respectiva

coordenada x e y.

A uniao de alguns nos faz surgir os elementos que tambem podem ser numerados sequenci-

almente de 1 ate m (m < n). Estes elementos podem ser descritos atraves do numero dos nos

que o compoem, e que podem ser associados ainda a tabela anterior, que possui as coordenadas

dos pontos. Esta tabela que apresenta as conexoes entre os nos de cada elemento e chamada

de conectividade.

A formulacao por elementos finitos exige alguns cuidados sendo o principal deles o fato

que a mesma deve obedecer uma sequencia de conexao entre os elementos, nao podendo ser

apresentada em uma ordem aleatoria. Normalmente e utilizado um ponto como referencia

inicial e um sentido de numeracao, horario ou anti-horario, para expressar a conectividade.

Considerando todo o procedimento nao faz diferenca qual o ponto que se adota como origem

para o elemento e nem o sentido de rotacao, no entanto, para todos os elementos deve ser

adotado o MESMO sentido de rotacao.

Um exemplo da montagem da tabela de coordenadas dos nos e de conectividade, pode ser

visto no exemplo que se segue.

Exemplo Numerico: Considere o volume total de chuva que incide sobre uma determinada

regiao de area A e na qual sao lidos uma serie de pontos ui. Os pontos numerados de 1 a 10

estao mostrados na tabela (6.2) e figura (6.6) sendo apresentados suas coordenadas x e y e os

nıveis de chuva medidos sobre o ponto.

O total de chuva Q pode ser definido como sendo:

Q =

∫A

u(x, y) · dA

sendo que u(x, y) representa o nıvel de precipitacao.


Figura 6.6: Discretizacao do problema em elementos

Com estes dados determine o valor total da precipitacao sobre a regiao considerada e o

ındice medio de precipitacao (umed).

Tabela 6.2: Nıvel de precipitacao em diversos pontos selecionados

No x[km] y[km] un[cm]

1 0 33, 3 4, 62

2 13, 2 62, 3 3, 81

3 39, 3 84, 5 4, 76

4 22, 2 30, 1 5, 45

5 49, 9 57, 6 4, 90

6 78, 8 78, 2 10, 35

7 39, 3 10, 0 4, 96

8 59, 7 34, 3 4, 26

9 73, 9 36, 2 18, 36

10 69, 8 5, 1 15, 69

Solucao: A primeira etapa consiste na montagem da tabela de conectividade para represen-

tar como os pontos estao interconectados. Esta tabela esta representada na tabela (6.3).

Aproximacao nodal:


Tabela 6.3: Tabela de Conectividade

Elemento 1 2 3 4

1 1 4 5 2

2 2 5 6 3

3 4 7 8 5

4 5 8 9 6

5 7 10 9 8

u(ξ, η) =⟨N1(ξ, η) N2(ξ, η) N3(ξ, η) N4(ξ, η)

⟩·

u1

u2

u3

u4

=

1

4

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩

u1

u2

u3

u4

Integracao:

Q =

∫A

u(x, y) · dx · dy =

∫A

u(ξ, η) · det [J ] · dξ · dη

ou, na forma matricial:

Q =1

4

∫ 1

−1

∫ 1

−1

(A0 + A1ξ + A2η)

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩dξ · dη

u1

u2

u3

u4

que depois de integrada:

Q =⟨−1

3A1 + A0 − 1

3A2

13A1 + A0 − 1

3A2

13A1 + A0 + 1

3A2 −1

3A1 + A0 + 1

3A2

⟩

u1

u2

u3

u4


que na forma polinomial pode ser escrita como:

Q =

(−1

3A1 + A0 −

1

3A2

)u1+

(1

3A1 + A0 −

1

3A2

)u2+(

1

3A1 + A0 +

1

3A2

)u3 +

(−1

3A1 + A0 +

1

3A2

)u4

e depois de reagrupada resulta em:

Q = A0 (u2 + u1 + u4 + u3) +A1

3(−u1 − u4 + u3 + u2) +

A2

3(u3 − u2 − u1 + u4)

ou ainda representado na forma:

Q = A0 U1 +A1

3U2 +

A2

3U3

sendo:

U1 = u2 + u1 + u4 + u3

U2 = −u1 − u4 + u3 + u2

U3 = u3 − u2 − u1 + u4

Para o primeiro elemento tem-se que:

No 1 ⇒ x1 = 0, y1 = 33, 3, u1 = 4, 62

No 4 ⇒ x2 = 22, 2, y2 = 30, 1, u2 = 5, 45

No 5 ⇒ x3 = 49, 9, y3 = 57, 6, u3 = 4, 90

No 2 ⇒ x4 = 13, 2, y4 = 62, 3, u4 = 3, 81

Utilizando das definicoes anteriormente obtidas para elementos de quatro nos:

A0 =1

8[(y4 − y2)(x3 − x1)− (y3 − y1)(x4 − x2)] ⇒ A0 = 228, 19

A1 =1

8[(y4 − y3) (x1 − x2)− (y2 − y1)(x3 − x4)] ⇒ A1 = 1, 638

A2 =1

8[(y1 − y4)(x2 − x3)− (y2 − y3)(x1 − x4)] ⇒ A2 = 55, 04

Consequentemente a vazao para o elemento 1 pode ser representada por:

Q1 = A0 (u2 + u1 + u4 + u3) +A1

3(−u1 − u4 + u3 + u2) +

A2

3(u3 − u2 − u1 + u4) ⇒

Q1 = 4261, 5 cm km2


Quanto ao calculo da area, tem-se que:

A =

∫A

dx · dy =

∫A

det [J ] · dξ · dη

A =

∫ 1

−1

∫ 1

−1

(A0 + A1ξ + A2η) dξ · dη =

∫ 1

−1

(2 · A0 + 2 · A2η)dη

A = 4 · A0

que implica em que:

A1 = 4 · A0 = 4 · 228, 19 = 912, 76 km2

Tabela 6.4: Resumo dos resultados do exemplo.

x1 y1 u1 x2 y2 u2 x3 y3 u3 x4 y4 u4

1 0,0 33,3 4,62 22,2 30,1 5,45 49,9 57,6 4,90 13,2 62,3 3,81

2 13,2 62,3 3,81 49,9 57,6 4,90 78,8 78,2 10,35 39,3 84,5 4,76

3 22,2 30,1 5,45 39,3 10,0 4,96 59,7 34,3 4,26 49,9 57,6 4,90

4 49,9 57,6 4,90 59,7 34,3 4,26 73,9 36,2 18,36 78,8 78,2 10,35

5 39,3 10,0 4,96 69,8 5,1 15,69 73,9 36,2 18,36 59,7 34,3 4,26

A0 A1 A2 U1 U2 U3 Q Area

1 228,19 1,64 55,04 18,78 1,92 -1,36 4261,41 912,74

2 241,65 -5,70 12,99 23,82 6,68 6,40 5771,07 966,59

3 217,56 -25,18 -14,01 19,57 -1,13 -1,25 4272,97 870,24

4 182,79 -65,72 29,70 37,87 7,37 19,55 6954,45 731,17

5 159,37 15,94 -66,85 43,27 24,83 1,97 6983,87 637,47

Soma 28243,78 4118,21

Media =

∑Qi∑Ai

= 6, 86 cm

6.8 O metodo dos resıduos ponderados

Considere a resolucao de um sistema fısico qualquer, do qual se conhece a equacao

diferencial. Trata-se de uma equacao em termos de derivadas parciais ou totais, linear ou nao

linear e de ordem m e pode ser expressa na forma:

L(u) + fv = 0 em todo o domınio V.

e sujeita as condicoes de contorno:

C(u) = fs no contorno S.


sendo u a variavel principal do problema e varia de acordo com a posicao no espaco. Define-se

a funcao residual como sendo:

R(u) = L(u) + fv

O metodo dos resıduos ponderados e normalmente expresso em sua forma integral:

W (u) =

∫V

〈ψ〉{R(u)}dV =

∫V

〈ψ〉{L(u) + fv}dV (6.5)

sendo ψ uma funcao peso dentre as possıveis e u a solucao para o problema que satisfaz as

condicoes de contorno C(u).

6.8.1 Transformacao integral: a integral por partes

As transformacoes integrais sao utilizadas para reduzir a ordem das equacoes diferenciais.

Relembrando ainda que para uma integracao generica:

u · v =

∫V

u · dv +

∫V

v · du =⇒∫V

u · dV = −∫V

v · du+ u · v

Caso se aplique este mesmo preceito para a equacao integral de resıduos ponderados,

eq. (6.5), este mesmo preceito quando esta acompanha uma derivada tem-se a forma fraca

da equacao, que e dada por:

• Unidimensional ∫ x2

x1

ψdu

dx· dx =−

∫ x2

x1

dψ

dxu · dx+ (ψ u)|x2

x1(6.6)∫ x2

x1

ψd2u

dx2· dx =−

∫ x2

x1

dψ

dx

du

dx· dx+

(ψdu

dx

)∣∣∣∣x2

x1

(6.7)

• Bidimensional (vide figura 6.7)

∫A

ψ∂u

∂x· dx dy =−

∫A

∂ψ

∂xu · dx dy +

∮S

(ψ u) · dy∫A

ψ∂u

∂y· dx dy =−

∫A

∂ψ

∂yu · dx dy +

∮S

(ψ u) · dx∫A

ψ

(∂u

∂x+∂u

∂y

)· dx dy =−

∫A

(∂ψ

∂x+∂ψ

∂y

)u · dx dy +

∮S

(ψ un) · dS (6.8)


∫A

ψ

(∂2u

∂x2+∂2u

∂y2

)· dx dy =−

∫A

(∂ψ

∂x

∂u

∂x+∂ψ

∂y

∂u

∂y

)· dx dy +

∮S

(ψ∂u

∂~n

)· dS (6.9)

sendo que

∮S

(ψ∂u

∂~n

)· dS =

∫S

(ψ∂u

∂y

)· dx+

∫S

(ψ∂u

∂x

)· dy

A forma fraca da equacao, obtida a partir da integracao por partes mostrada anteriormente,

apresenta algumas caracterısticas e requisitos diferenciados:

(a) a ordem da maior derivada de variavel de interesse u e reduzida, o que relaxa a condicao

de continuidade necessaria para a convergencia;

(b) algumas das condicoes de contorno aparecem diretamente na expressao geral;

(c) reduz a ordem necessaria para a funcao de interpolacao entre os nos;

(d) aumenta a ordem necessaria para a funcao peso da solucao.

6.8.2 Metodo de Galerkin

Existem varias opcoes de escolha para a funcao peso, sendo que cada uma pode apre-

sentar melhores ou piores resultados de acordo com o tipo de problema a ser solucionado. Esta

escolha pode ser auxiliada com a utilizacao de referencias basicas de elementos finitos, como

(Dhatt & Touzot, 1984).

No entanto existe um esquema de escolha de funcao peso que se destaca pela grande uti-

lizacao nos mais diversos tipos de problemas. Este esquema de solucao, conhecido como o

Figura 6.7: Composicao vetorial para integral sobre a superfıcie


metodo de Galerkin, se utiliza da mesma funcao usada na interpolacao da variavel principal

(u) para a funcao de interpolacao (ψ), ou seja:

ψ = 〈N〉

6.9 Tratamento das condicoes de contorno

O tratamento das condicoes de contorno atraves do metodo dos elementos finitos nao e

uma tarefa simples, ele vem acompanhado de uma serie de operacoes que nao sao imediatas.

Para simplificar esta tarefa sera utilizado neste estudo um tratamento simplificando, vide (Pa-

tankar, 1980), das condicoes de contorno onde e definida uma condicao de contorno generica

na forma:∂u

∂n= αc u+ βc (6.10)

Esta equacao e capaz de representar qualquer condicao de contorno usual (das tres especies)

e outras ainda podem ser aproximadas. Para esta representacao e necessario que:

Condicao de contorno da primeira especie onde o valor da variavel u e especificado na

fronteira. Isto e possıvel de se obter substituindo os valores constantes na eq. (6.10) por:

αc → −∞ e βc →∞×ue sendo ue o valor especificado para a variavel u na fronteira. Este

esquema pode ser facilmente explicado indicando uma tendencia de para qualquer valor

resultante da derivada (∂u/∂~n),o valor de u estara muito proximo do valor especificado

ue.

Uma outra possibilidade para este tipo de condicao de contorno, e ate mesmo mais comum

no tratamento de condicoes deste tipo e a substituicao da linha da matriz equivalente ao

no pela respectiva igualdade, ou sejam, ui = ue onde i e o numero do no considerado.

Condicao de contorno de segunda especie onde o valor conhecido e a taxa de variacao,

ou fluxo, da variavel considerada. Neste caso substitui-se os valores da eq. (6.10) por:

αc = 0 e βc = u′e onde u′e e o valor da derivada na fronteira considerada, nao se esquecendo

de considerar o sinal da derivada, que deve estar de acordo com o sentido de ~n.

Condicao de contorno de terceira especie e uma condicao tıpica de transmissao de calor

com a condicao de conveccao e deve expressar a igualdade (caso em que se trata de uma

fronteira final de x):

−k∂T∂x

= h(T − T∞) =⇒ ∂T

∂x= −h

kT +

hT∞k


Figura 6.8: Tipos de fronteira para a discretizacao da C.C. de terceira especie.

e se igualarmos esta expressao a eq. (6.10) obtem-se que:

αc = −hk

e βc =hT∞k

E importante ressaltar que caso se trate de uma fronteira inicial de x, como mostrada na

figura (6.8), ~n e ~i tem direcoes opostas, no entanto o sinal tambem Da expressao para

conveccao tambem se inverte pela posicao da fronteira. Desta forma, o resultado final e

sempre o mesmo e independe de se tratar de uma fronteira inicial ou final.

6.9.1 Implementacao da condicao de contorno generalizada

Vejam o caso de implementacao desta condicao de contorno num elemento quadrilateral

simples, figura (6.5), e cuja funcao de interpolacao ja foi obtida sendo:

〈N〉 =1

4

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩Desta forma, a variavel de interesse u pode ser expressa em termos da funcao de interpolacao

atraves da aproximacao nodal:

u(ξ, η) =⟨N1(ξ, η) N2(ξ, η) N3(ξ, η) N4(ξ, η)

⟩·

u1

u2

u3

u4

=

1

4

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩

u1

u2

u3

u4


Pode-se escolher qualquer uma das fronteiras para realizar a integracao sobre a fronteira,

mudando apenas o sentido da derivada e os pontos considerados e a direcao sobre a qual e feita

a derivada (ξ ou η).

∂ 〈N〉∂ξ

=1

4

∂

∂ξ

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩=

1

4

⟨−1 + η 1− η 1 + η −1− η

⟩∂ 〈N〉∂η

=1

4

∂

∂η

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩=

1

4

⟨−1 + ξ −1− ξ 1 + ξ 1− ξ

⟩

Para a face superior do elemento quadrado (η = 1) unindo os nos 3 e 4, a derivada normal

passa a ser ∂u/∂η e a integral sobre a superfıcie passa a ser:

∮S

(ψ∂u

∂~n

)dS =

1∫ξ=−1

ψ · (αu+ β) · dξ =

1∫ξ=−1

1

4

⟨ 0

0

2 · (1 + ξ)

2 · (1− ξ)

⟩(αu+ β) · dξ

=α

16

1∫ξ=−1

⟨ 0

0

2 (1 + ξ)

2 (1− ξ)

⟩⟨0 0 2(1 + ξ) 2(1− ξ)

⟩dξ

u1

u2

u3

u4

+β

4

1∫ξ=−1

⟨ 0

0

2 (1 + ξ)

2 (1− ξ)

⟩dξ

=α

16

1∫ξ=−1

0 0 0 0

0 0 0 0

0 0 (2 + 2ξ)2 (2 + 2ξ) (2− 2ξ)

0 0 (2 + 2ξ) (2− 2ξ) (2− 2ξ)2

dξ

u1

u2

u3

u4

+β

4

0

0

2ξ + ξ2

2ξ − ξ2

1

−1


∮S

(ψ∂u

∂~n

)dS =

α

16

0 0 0 0

0 0 0 0

0 0 16(2 + 2ξ)3 −4

3ξ3 + 4ξ

0 0 −43ξ3 + 4ξ −1

6(2− 2ξ)3

1

−1

+β

4

0

0

4

4

=

α

16

0 0 0 0

0 0 0 0

0 0 323

163

0 0 163

323

u1

u2

u3

u4

+ β

0

0

1

1

=

0 0 0 0

0 0 0 0

0 0 2α3

α3

0 0 α3

2α3

u1

u2

u3

u4

+

0

0

β

β

=

0

023αu3 + 1

3αu4

13αu3 + 2

3αu4

+

0

0

β

β

Considerando o comprimento total ` sendo representado por ξ variando entre -1 e 1 (2

unidades), tem-se:

∮S

(ψ∂u

∂~n

)dS =

0

023αu3 + 1

3αu4

13αu3 + 2

3αu4

+

0

0

β

β

`

2=`

6

0

0

2αu3 + αu4 + 3β

αu3 + 2αu4 + 3β

6.10 Integracao Numerica

Ate agora todas as integrais a serem resolvidas o foram analiticamente, no entanto nao

pode deixar de se considerar a integral numerica como ferramenta importante neste tipo de

solucao. Da forma ja vista anteriormente pode-se dizer que a integral a ser solucionada e:∫e

g(x) · dx =

∫ 1

−1

G(ξ, η) · dξdη ≈ns∑

j=1

ωj

ns∑k=1

ωk G(ξj, ηk) =ns∑

j=1

ns∑k=1

ωj ωk G(ξj, ηk)

Os pontos de integracao sao obtidos assim como a funcao peso para diversas geometrias sendo as

mais comuns a quadratura de Gauss e o metodo de Newton-Cotes para elementos quadrilaterais


Esquema Num. pontos Posicoes Peso

4 ±1/√

3, ±1/√

3 1

3

1/6, 1/62/3, 1/61/6, 2/3 1/6

Figura 6.9: Pontos de integracao para elementos triangulares e quadrangulares.

e o metodo de Gauss-Radau para elementos triangulares. Existem uma serie de pontos de

integracao escolhidos com funcoes peso tambem determinadas. Dentre estas possibilidades

estao mostradas na figura (6.9) as mais utilizadas que sao as de 4 pontos para elementos

quadrangulares e 3 pontos para elementos triangulares.

6.11 Exemplos de Aplicacao

Um grande forno industrial e suportado por uma longa coluna de tijolos refratarios com

1 m por 1 m de lado. Durante a operacao em regime estacionario, as condicoes sao tais que que

tres superfıcies da coluna sao mantidas a 500 K e a outra superfıcie e exposta ha uma corrente

de ar com T∞ = 300 K e o coeficiente de pelıcula h = 10 W/m2K. Usando uma malha com

δx = δy = 0.25 m, determine a distribuicao bidimensional de temperatura bem como a taxa

de calor transmitida ao ar por unidade de comprimento da coluna. (problema apresentado por

(Incropera & DeWitt, 1999), pg 101)

Solucao: Para a solucao deste problema sera utilizada uma malha relativamente grosseira

com 16 nos e 8 elementos, mostrada a seguir:


Se trata de um problema de conducao simples em um meio isotropico e a equacao geral que

rege o problema e dada por:∂2T

∂x2+∂2T

∂y2= 0

Considerando a sua forma de resıduos ponderados, tem-se que:∫e

Ψ

(∂2T

∂x2+∂2T

∂y2

)dVe = 0

e reduzida a sua forma fraca resulta em:

−∫

e

(∂Ψ

∂x

∂T

∂x+∂Ψ

∂y

∂T

∂y

)dVe +

∮e

(Ψ∂T

∂n

)dSe = 0

Usando a aproximacao nodal e multiplicando a expressao por (-1):∫e

(∂Ψ

∂x

∂

∂x〈N〉+

∂Ψ

∂y

∂

∂y〈N〉

)dVe · {T} −

∮e

(Ψ∂T

∂n

)dSe = 0

e depois, aplicando Galerkin:∫e

(∂

∂x{N} ∂

∂x〈N〉+

∂

∂y{N} ∂

∂y〈N〉

)dVe · {T} −

∮e

({N}∂T

∂n

)dSe = 0

que e a expressao geral para o problema e que pode ser representada, em termos matriciais, na

forma: ∫e

({N} [ ∂

∂x∂∂y

] ·

[∂∂x∂∂y

]〈N〉

)dVe · {T} −

∮e

({N} ∂T

∂n

)dSe = 0

Pelas propriedades matriz transposta, sabe-se que:[∂∂x∂∂y

]〈N〉 = [B] =⇒ [B]T = {N}

[∂∂x

∂∂y

]


e portanto a representacao matricial pode ser feita atraves da expressao:∫e

[B]T · [B] dVe · {T}+

∮e

({N} ∂T

∂n

)dSe = 0

Sendo que para os elementos internos a integral de superfıcie e desprezada resultando

em: ∫e

[B]T · [B] dVe · {T} = 0

que em termos de variaveis locais se reduz a:∫ 1

−1

∫ 1

−1

[B]T · [B] det[J ] · dξ · dη · {T} = 0

onde a matriz [B] e escrita como:

[B] =

[ ⟨∂N∂x

⟩⟨∂N∂y

⟩ ] =

[∂ξ∂x

∂η∂x

∂ξ∂y

∂η∂y

] ⟨∂N∂ξ

⟩⟨∂N∂η

⟩ = [J ]−1 · [Bξ] = [Q] · [Bξ]

Reescrevendo a expressao acima em termos de variaveis locais tem-se:∫ 1

−1

∫ 1

−1

([Q] [Bξ])T [Q] [Bξ] det[J ]·dξdη·{T} = 0 ⇒

∫ 1

−1

∫ 1

−1

[Bξ]T [Q]T [Q] [Bξ] det[J ]·dξdη·{T} = 0

Para elementos quadrilateros as funcoes de aproximacao sao dadas por:

〈N〉 =1

4

⟨(1− ξ)(1− η) (1 + ξ)(1− η) (1 + ξ)(1 + η) (1− ξ)(1 + η)

⟩⟨∂N

∂ξ

⟩=

1

4

⟨−1 + η 1− η 1 + η −1− η

⟩⟨∂N

∂η

⟩=

1

4

⟨−1 + ξ −1− ξ 1 + ξ 1− ξ

⟩[Bξ] =

1

4

[−1 + η 1− η 1 + η −1− η

−1 + ξ −1− ξ 1 + ξ 1− ξ

]

Quanto ao Jacobiano, ele pode ser representado a partir da matriz [B], na forma:

[J ] = [Bξ] ·

x1 y1

x2 y2

x3 y3

x4 y4

e [Q] = [J ]−1


O procedimento descrito ate aqui e generico e e identico para qualquer problema de conducao

simples. A partir de agora o procedimento passa a representar nao so as caracterısticas do

problema, como tambem as caracterısticas de malha e discretizacao representados pelo numero

de nos e tipo e numero de elementos.

Para o problema considerado, e necessario montar uma tabela com o numero de nos, que

pode ser dada por:

No x y

1 0 0

2 0,25 0

3 0,5 0

4 0 0,25

5 0,25 0,25

6 0,5 0,25

7 0 0,5

8 0,25 0,5

9 0,5 0,5

10 0 0,75

11 0,25 0,75

12 0,5 0,75

13 0 1

14 0,25 1

15 0,5 1

E tambem e necessario montar uma tabela que represente como estes nos estao associados

para formar os elementos, e a chamada tabela de conectividade:

elem. 1 2 3 4

1 1 2 5 4

2 2 3 6 5

3 4 5 8 7

4 5 6 9 8

5 7 8 11 10

6 8 9 12 11

7 10 11 14 13

8 11 12 15 14


Feita a conectividade e possıvel notar que todos os elementos sao iguais e a solucao de um

sera identica a todos os demais. Assim, tomando o elemento 1 tem-se:

1. Integracao Numerica: e feita sobre a equacao anteriormente obtida, no domınio elementar

transformando-se a integral em uma soma, avaliada nos pontos determinados e com a

funcao peso adequada:

∫ 1

−1

∫ 1

−1

[Bξ]T [Q]T [Q] [Bξ] det[J ]·dξdη·{T} =

np∑i=1

np∑j=1

(ωiωj [Bξ]

T [Q]T [Q] [Bξ] det[J ])·{T}

sendo que neste caso foram escolhidos quatro pontos de integracao:

• Ponto 1: ωi = 1, ωj = 1, ξ = 1/√

3, η = 1/√

3

[Bξ] =1

4

[−1 + 1

3

√3 1− 1

3

√3 1 + 1

3

√3 −1− 1

3

√3

−1 + 13

√3 −1− 1

3

√3 1 + 1

3

√3 1− 1

3

√3

]

=

[−0. 105 66 0. 105 66 0. 394 35 −0. 394 35

−0. 105 66 −0. 394 35 0. 394 35 0. 105 66

]e o Jacobiano para este elemento:

[J ] = [Bξ] ·

x1 y1

x2 y2

x3 y3

x4 y4

=

[−0. 105 66 0. 105 66 0. 394 35 −0. 394 35

−0. 105 66 −0. 394 35 0. 394 35 0. 105 66

]0 0

0.25 0

0.25 0.25

0 0.25

=

0. 125 0

0 0. 125sendo o valor do det [J ] = 1. 562 5× 10−2 e o [J ]−1 = Q =

[8.0 0

0 8.0

]e assim

[Bξ]T [Q]T [Q] [Bξ] =

−0. 105 66 −0. 105 66

0. 105 66 −0. 394 35

0. 394 35 0. 394 35

−0. 394 35 0. 105 66

[

8.0 0

0 8.0

]·

·

[8.0 0

0 8.0

][−0. 105 66 0. 105 66 0. 394 35 −0. 394 35

−0. 105 66 −0. 394 35 0. 394 35 0. 105 66

]

=

1. 429 1. 952 2 −5. 333 4 1. 952 2

1. 952 2 10. 667 −7. 286 1 −5. 333 4

−5. 333 4 −7. 286 1 19. 906 −7. 286 1

1. 952 2 −5. 333 4 −7. 286 1 10. 667


Que se usada toda a expressao para o primeiro ponto de integracao resulta em:

[Vi] = ωiωj [Bξ]T [Q]T [Q] [Bξ] det[J ]

sendo

[V1] = 1 · 1 ·

1. 429 1. 952 2 −5. 333 4 1. 952 2

1. 952 2 10. 667 −7. 286 1 −5. 333 4

−5. 333 4 −7. 286 1 19. 906 −7. 286 1

1. 952 2 −5. 333 4 −7. 286 1 10. 667

· 1. 562 5× 10−2

=

2. 232 8× 10−2 3. 050 3× 10−2 −8. 333 4× 10−2 3. 050 3× 10−2

3. 050 3× 10−2 0. 166 67 −0. 113 85 −8. 333 4× 10−2

−8. 333 4× 10−2 −0. 113 85 0. 311 03 −0. 113 85

3. 050 3× 10−2 −8. 333 4× 10−2 −0. 113 85 0. 166 67

• Ponto 2: ωi = 1, ωj = 1, ξ = 1/

√3, η = −1/

√3 obtem-se det [J ] = 0.015625

[V2] =

0. 166 67 −0. 113 84 −8. 333 1× 10−2 3. 050 2× 10−2

−0. 113 84 0. 311 02 −0. 113 84 −8. 333 1× 10−2

−8. 333 1× 10−2 −0. 113 84 0. 166 67 3. 050 2× 10−2

3. 050 2× 10−2 −8. 333 1× 10−2 3. 050 2× 10−2 2. 232 8× 10−2

• Ponto 3: ωi = 1, ωj = 1, ξ = −1/

√3, η = −1/

√3

[V3] =

0. 311 02 −0. 113 84 −8. 333 1× 10−2 −0. 113 84

−0. 113 84 0. 166 67 3. 050 2× 10−2 −8. 333 1× 10−2

−8. 333 1× 10−2 3. 050 2× 10−2 2. 232 8× 10−2 3. 050 2× 10−2

−0. 113 84 −8. 333 1× 10−2 3. 050 2× 10−2 0. 166 67

• Ponto 4: ωi = 1, ωj = 1, ξ = −1/

√3, η = 1/

√3

[V4] =

0. 166 67 3. 050 2× 10−2 −8. 333 1× 10−2 −0. 113 84

3. 050 2× 10−2 2. 232 8× 10−2 3. 050 2× 10−2 −8. 333 1× 10−2

−8. 333 1× 10−2 3. 050 2× 10−2 0. 166 67 −0. 113 84

−0. 113 84 −8. 333 1× 10−2 −0. 113 84 0. 311 02

E se feita a somatoria resulta na matriz elementar:

4∑i=1

[Vi] =

23−1

6−1

3−1

6

−16

23−1

6−1

3

−13−1

623−1

6

−16−1

3−1

623


2. O mesmo procedimento deve ser repetido para todos os elementos. Neste caso, entretanto,

tem-se o fato conveniente de que todos os elementos se encontram na mesma disposicao e

tem as mesmas caracterısticas geometricas e, portanto, esta matriz encontrada vale para

todos os elementos.

3. A proxima etapa consiste na montagem da matriz global a partir das matrizes elementares.

Neste caso utiliza-se a tabela de conectividade para converter da numeracao local para a

numeracao global, e cada linha e coluna da numeracao local passam a responder por uma

outra na numeracao global. Vejamos para como exemplo para o primeiro elemento (onde

a matriz global A ainda esta vazia):

numeracao local 1 2 3 4

numeracao global 1 2 5 4

A =

0 + 0. 666 69 0− 0. 166 68 0 0− 0. 166 68 0− 0. 333 33 0 0 0 0 0 0 0 0 0 00− 0. 166 68 0 + 0. 666 69 0 0− 0. 333 33 0− 0. 166 68 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 00− 0. 166 68 0− 0. 333 33 0 0 + 0. 666 69 0− 0. 166 68 0 0 0 0 0 0 0 0 0 00− 0. 333 33 0− 0. 166 68 0 0− 0. 166 68 0 + 0. 666 69 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0


que depois de incluıdos todos os elementos resulta em:

[A] =

23

−16

0 −16−1

30 0 0 0 0 0 0 0 0 0

−16

43

−16−1

3−1

3−1

30 0 0 0 0 0 0 0 0

0 −16

23

0 −13−1

60 0 0 0 0 0 0 0 0

−16−1

30 4

3−1

30 −1

6−1

30 0 0 0 0 0 0

−13−1

3−1

3−1

383

−13−1

3−1

3−1

30 0 0 0 0 0

0 −13−1

60 −1

343

0 −13−1

60 0 0 0 0 0

0 0 0 −16−1

30 4

3−1

30 −1

6−1

30 0 0 0

0 0 0 −13−1

3−1

3−1

383

−13−1

3−1

3−1

30 0 0

0 0 0 0 −13−1

60 −1

343

0 −13−1

60 0 0

0 0 0 0 0 0 −16−1

30 4

3−1

30 −1

6−1

30

0 0 0 0 0 0 −13−1

3−1

3−1

383

−13−1

3−1

3−1

3

0 0 0 0 0 0 0 −13−1

60 −1

343

0 −13−1

6

0 0 0 0 0 0 0 0 0 −16−1

30 2

3−1

60

0 0 0 0 0 0 0 0 0 −13−1

3−1

3−1

643

−16

0 0 0 0 0 0 0 0 0 0 −13−1

60 −1

623

4. Deve-se entao incluir as condicoes de contorno que de acordo com a proposicao geral pode

ser expressa na forma da eq.() e respeitando a tabela abaixo(lembrando que l = 0.25 m):

No 1 No 2 α β −α·l3

−α·l6

β·l2

1 2 −101

10·3001

0. 833 33 0. 416 67 375

2 3 −101

10·3001

0. 833 33 0. 416 67 375

4 7 106 500× 106 8. 333 3× 104 4. 166 7× 104 6. 25× 107

7 10 106 500× 106 8. 333 3× 104 4. 166 7× 104 6. 25× 107

10 13 106 500× 106 8. 333 3× 104 4. 166 7× 104 6. 25× 107

13 14 106 500× 106 8. 333 3× 104 4. 166 7× 104 6. 25× 107

14 15 106 500× 106 8. 333 3× 104 4. 166 7× 104 6. 25× 107

sendo que na linha relativa a cada condicao devem ser somados os respectivos valores.

Veja o exemplo para uma condicao i-j:

i j

i −α·l3

−α·l6

j −α·l6

−α·l3

= β·l2

β·l2

Implementando numa matriz terıamos:


−[C] =

266666666666666666666666666666664

0.833333 0.416667 0 0 0 0 0 0 0 0 0 0 0 0 0

0.416667 1.66667 0.416667 0 0 0 0 0 0 0 0 0 0 0 0

0 0.416667 0.833333 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 83333.3 0 0 41666.7 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 41666.7 0 0 166667. 0 0 41666.7 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 41666.7 0 0 166667. 0 0 41666.7 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 41666.7 0 0 166667. 41666.7 0

0 0 0 0 0 0 0 0 0 0 0 0 41666.7 166667. 41666.7

0 0 0 0 0 0 0 0 0 0 0 0 0 41666.7 83333.3

377777777777777777777777777777775

{Bc} =

375

750

375

6.25× 107

0

0

1.25× 108

0

0

1.25× 108

0

0

1.25× 108

1.25× 108

6.25× 107

E a matriz final:

[M ] = [A] + [C]

=

266666666666666666666666666666664

1.5 0.25 0 − 16

− 13

0 0 0 0 0 0 0 0 0 0

0.25 3. 0.25 − 13

− 13

− 13

0 0 0 0 0 0 0 0 0

0 0.25 1.5 0 − 13

− 16

0 0 0 0 0 0 0 0 0

− 16

− 13

0 83334.7 − 13

0 41666.5 − 13

0 0 0 0 0 0 0

− 13

− 13

− 13

− 13

83

− 13

− 13

− 13

− 13

0 0 0 0 0 0

0 − 13

− 16

0 − 13

43

0 − 13

− 16

0 0 0 0 0 0

0 0 0 41666.5 − 13

0 166668. − 13

0 41666.5 − 13

0 0 0 0

0 0 0 − 13

− 13

− 13

− 13

83

− 13

− 13

− 13

− 13

0 0 0

0 0 0 0 − 13

− 16

0 − 13

43

0 − 13

− 16

0 0 0

0 0 0 0 0 0 41666.5 − 13

0 166668. − 13

0 41666.5 − 13

0

0 0 0 0 0 0 − 13

− 13

− 13

− 13

83

− 13

− 13

− 13

− 13

0 0 0 0 0 0 0 − 13

− 16

0 − 13

43

0 − 13

− 16

0 0 0 0 0 0 0 0 0 41666.5 − 13

0 166667. 41666.5 0

0 0 0 0 0 0 0 0 0 − 13

− 13

− 13

41666.5 166668. 41666.5

0 0 0 0 0 0 0 0 0 0 − 13

− 16

0 41666.5 83334.

377777777777777777777777777777775


Tabela 6.5: Solucoes diversas obtidas para o exemplo considerado.

Num. Solucao T2 T3

1 Diferencas finitas ((Incropera & DeWitt, 1999)) 357,0 339,1

2 Elementos finitos (15 pontos) 340,9 331,0



5 Idem a 2, mas utilizando T1 = 500 330,2 338,4

A matriz de carga {B}, nao sofre alteracoes neste caso ja que ela e alterada apenas pelascondicoes de contorno (sendo nula para a matriz elementar). A solucao do sistema ficasendo:

T1 → 341.57, T2 → 340.933, T3 → 331.023, T4 → 499.999, T5 → 417.765, T6 → 405.08, T7 → 500,

T8 → 468.708, T9 → 454.805, T10 → 500, T11 → 488.392, T12 → 483.626, T13 → 500, T14 → 500, T15 → 500

Nao existe solucao analıtica para este problema e nao existe forma de determinar o valor

exato da temperatura para estes pontos, no entanto o refinamento da malha tende a fornecer

solucoes mais proximas da solucao real. Sao mostradas na tabela (6.5) uma serie de solucoes

obtidas para este problema.

CAPITULO 7

Metodo de Diferencas Finitas

O metodo de diferencas finitas foi, durante muito tempo, o unico metodo numerico

utilizado para problemas fluido-termicos. Embora recentemente este venha perdendo espaco

para outros metodos ainda e muito utilizado nesta area.

A solucao de qualquer problema em diferencas finitas se baseia no princıpio da serie de

Taylor em que um ponto qualquer pode ser expresso em funcao dos seu pontos proximos, ou

seja:

f(x+ ∆x) = f(x) +df(x)

dx∆x+

d2f(x)

dx2

(∆x)2

2!+d3f(x)

dx3

(∆x)3

3!+ · · ·

que e na realidade uma serie infinita.

Para exemplificar esta aproximacao imaginemos que a aproximacao de pontos Ti e Ti+1 a

partir de ponto Ti. Imaginando que eles sao igualmente espacados e a distancia entre os pontos

e ∆x. Neste caso:

Ti+1 = Ti +∂T

∂x

∣∣∣∣i

∆x+∂2T

∂x2

∣∣∣∣i

(∆x)2

2!+∂3T

∂x3

∣∣∣∣i

(∆x)3

3!+O(∆x)4 (7.1)

e o outro ponto, considerado na direcao negativa de x, pode ser expresso por:

Ti−1 = Ti −∂T

∂x

∣∣∣∣i

∆x+∂2T

∂x2

∣∣∣∣i

(∆x)2

2!− ∂3T

∂x3

∣∣∣∣i

(∆x)3

3!+O(∆x)4 (7.2)

Como se tratam de series infinitas, pode-se fazer a aproximacao com a ordem do erro que se

desejar. Grande numero de trabalhos envolvendo a primeira derivada foram resolvidos atraves

94


de esquemas de primeira ordem (O(∆x)2), mas atualmente sao muito utilizados os esquemas

de segunda ordem1. A obtencao das expressoes discretas e facil a partir dai.

No caso das aproximacoes acima, se utilizarmos ate o termo de primeira ordem somente,

cada uma pode resultar em aproximacoes para a primeira derivada. No caso da equacao (7.1):

∂T

∂x

∣∣∣∣i

=Ti+1 − Ti

∆x(7.3)

que e conhecida como a derivada avancada. Da mesma forma utilizando-se da expressao (7.2)

pode-se chegar numa outra expressao, denominada por aproximacao recuada:

∂T

∂x

∣∣∣∣i

=Ti − Ti−1

∆x(7.4)

Estas duas expressoes sao aproximacoes de primeira ordem e portanto, em condicoes nor-

mais, apresentam piores resultados que as expressoes de segunda ordem, as mais utilizadas. As

expressoes de segunda ordem sao as mais utilizadas, principalmente em problemas lineares, por

fornecerem melhores resultados. A expressao para a primeira derivada considerando os termos

de ate segunda ordem pode ser obtida trabalhando-se com as expressoes, basicamente fazendo

a subtracao da equacao (7.2) da (7.1). Desta forma obtem-se:

Ti+1 − Ti−1 = 2∂T

∂x

∣∣∣∣i

∆x

que rearranjada resulta em:∂T

∂x

∣∣∣∣i

=Ti+1 − Ti−1

2 ∆x(7.5)

Como se esta trabalhando com os termos de ate segunda ordem pode-se obter tambem uma

expressao para a segunda derivada. Isto e possıvel fazendo a soma das equacoes (7.1) e (7.2),

donde obtem-se:

Ti+1 + Ti−1 = 2Ti + 2∂2T

∂x2

∣∣∣∣i

(∆x)2

2!

que rearranjada na forma desejada resulta em:

∂2T

∂x2

∣∣∣∣i

=Ti+1 − 2Ti + Ti−1

(∆x)2(7.6)

Como pode ser visto este procedimento e extremamente generico e pode ser utilizado em

esquemas de ordem superiores sem maiores problemas, dependendo apenas da utilizacao de um

maior numero de pontos. Repare que para chegar nas expressoes de primeira ordem, apenas uma

das expansoes em Serie de Taylor precisava ser conhecida, a partir das expressoes de segunda

ordem ja sao necessarias duas equacoes para resolver o problema (tanto as aproximacoes de

Ti+1 como Ti−1). Este mesmo raciocınio continua valendo para expressoes de ordem superior.

1esquemas de ordem superiores tambem tem sido motivos de estudos recentemente.


Exemplo: A partir de uma malha irregular como a mostrada a seguir obtenha as expressoes

de segunda ordem para a primeira e segunda derivadas.

e e eTi−1 Ti Ti+1

∆x2 ∆x1

Solucao: Expandindo cada um dos pontos por serie de Taylor, a partir do ponto i:

Ti+1 = Ti +∂T

∂x

∣∣∣∣i

∆x1 +∂2T

∂x2

∣∣∣∣i

(∆x1)2

2!+O(∆x)3

Ti−1 = Ti −∂T

∂x

∣∣∣∣i

∆x2 +∂2T

∂x2

∣∣∣∣i

(∆x2)2

2!+O(∆x)3

A partir desta expansao e possıvel obter a expressao para a primeira derivada fazendo-se a

soma da primeira equacao multiplicada por (∆x2)2 e da segunda multiplicada por −(∆x1)

2:

Ti+1 (∆x2)2 = Ti (∆x2)

2 +∂T

∂x

∣∣∣∣i

∆x1(∆x2)2 +

∂2T

∂x2

∣∣∣∣i

(∆x1 ∆x2)2

2

−Ti−1 (∆x1)2 = −Ti(∆x1)

2 +∂T

∂x

∣∣∣∣i

∆x2 (∆x1)2 − ∂2T

∂x2

∣∣∣∣i

(∆x2 ∆x1)2

2

Ti+1 (∆x2)2 − Ti−1 (∆x1)

2 = Ti (∆x22 −∆x2

1) +∂T

∂x

∣∣∣∣i

(∆x1∆x22 + ∆x2 ∆x2

1)

Rearranjando a expressao de forma a isolar o valor da primeira derivada obtem-se:

∂T

∂x

∣∣∣∣i

=Ti+1 (∆x2)

2 − Ti−1 (∆x1)2 + Ti (∆x

21 −∆x2

2)

∆x1∆x2(∆x2 + ∆x1)

Da mesma forma pode se obter a expressao para a segunda derivada, mas agora pena

multiplicando a expansao em Serie de Taylor de Ti+1 por ∆x2 e a de Ti−1 por ∆x1. Assim:

Ti+1 ∆x2 = Ti ∆x2 +∂T

∂x

∣∣∣∣i

∆x1 ∆x2 +∂2T

∂x2

∣∣∣∣i

∆x21 ∆x2

2

Ti−1 ∆x1 = Ti ∆x1 −∂T

∂x

∣∣∣∣i

∆x2 ∆x1 +∂2T

∂x2

∣∣∣∣i

∆x1 ∆x22

2

Ti+1 ∆x2 + Ti−1 ∆x1 = Ti (∆x1 + ∆x2) +∂2T

∂x2

∣∣∣∣i

∆x1 ∆x2 (∆x2 + ∆x1)

2

E desta forma a expressao final para a segunda derivada fica:

∂2T

∂x2

∣∣∣∣i

= 2Ti+1 ∆x2 − Ti (∆x1 + ∆x2) + Ti−1 ∆x1

∆x1 ∆x2 (∆x2 + ∆x1)


Repare que estas duas expressoes obtidas se reduzem as expressoes (7.5) e (7.6) respectiva-

mente quando ∆x1 = ∆x2 = ∆x.

7.1 Aplicacao da formulacao de diferencas finitas em sis-

temas de equacoes diferenciais

A utilizacao da metodologia de diferencas finitas em equacoes diferenciais em um dados

problema e bastante simples, consistindo basicamente das seguintes etapas:

1. Discretizacao do domınio, que implica na escolha dos pontos onde serao obtidos valores

para a solucao.

2. Substituicao das derivadas por sua aproximacao calculada a partir de uma serie de Taylor

para um ponto generico i. A aproximacao central de segunda ordem para uma derivada

de primeira ordem e dada pela equacao (7.5), por exemplo. Outras aproximacoes ja foram

mostradas.

3. Transformacao das equacoes diferenciais em sistema de equacoes baseados nos valores dos

pontos discretos. Para isto e necessario a utilizacao da equacao generica para o ponto i,

obtida no passo anterior para cada um dos pontos discretos.

4. Discretizacao e aplicacao das condicoes de contorno nos pontos da fronteira do domınio

e, com isto, obter equacoes discretizadas neste ponto.

5. Solucao do sistema de equacoes obtido.

Exemplo: Obtenha a equacao geral discretizada para a equacao de conducao bidimensional:

∂2T

∂x2+∂2T

∂y2= 0

Solucao: Utilizando as aproximacoes centrais para a segunda derivada em cada uma das

direcoes, equacao (7.6) tem-se que:

∂2T

∂x2=

Ti+1,j − 2Ti,j + Ti−1,j

(∆x)2

∂2T

∂y2=

Ti,j+1 − 2Ti,j + Ti,j−1

(∆y)2


Desta forma a equacao geral fica:

Ti+1,j − 2Ti,j + Ti−1,j

(∆x)2+Ti,j+1 − 2Ti,j + Ti,j−1

(∆y)2= 0

ou rearranjada na forma de mostrar os coeficientes de cada termo:(2

(∆x)2+

2

(∆y)2

)Ti,j =

1

(∆x)2Ti+1,j +

1

(∆x)2Ti−1,j +

1

(∆y)2Ti,j+1 +

1

(∆y)2Ti,j−1

Por curiosidade repare que se ∆x = ∆y, a expressao anterior fica:

Ti,j =Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1

4

Exemplo: Resolva o problema proposto anteriormente para a aleta com fluxo de calor na

base conhecido, resolvido por Runge Kutta, e compare os resultados com os obtidos atraves do

metodo de diferencas finitas (problema da pagina 32).

Solucao: A equacao diferencial da aleta foi obtida no exemplo anterior e e dada por:

d2T

dr2=

2h

k t(T − T∞)− 1

r

dT

dr

que discretizada pela aproximacao central:

Ti+1 − 2Ti + Ti−1

(∆r)2=

2h

k t(Ti − T∞)− 1

ri

Ti+1 − Ti−1

2 ∆r

Rearranjando a expressao tem-se:(1

(∆r)2+

1

2 ri ∆r

)Ti+1 −

(2

(∆r)2+

2h

k t

)Ti +

(1

(∆r)2− 1

2 ri ∆r

)Ti−1 = −2h

k tT∞

Considerando os valores de incremento de raio constante ∆r iguais ao usados no Runge

Kutta de 0,01 m e que o r0 = 0, 06 m e os valores sao identicos ao do exercıcio anterior:(10000 +

50

ri

)Ti+1 − 20093, 02Ti +

(10000− 50

ri

)Ti−1 = −2325, 6 (7.7)

Tem-se tambem que as condicoes de contorno sao conhecidas:

T = 50◦C em r = 6 cm =⇒ T1 = 50

q = 120 W em r = 6 cm =⇒ ∂T

∂r= − 120

π.0, 12.0, 005.215= −296, 10

Discretizando a condicao de contorno do fluxo de calor especificado pela discretizacao

avancada da primeira derivada tem-se que:

dT

dr

∣∣∣∣r=6

=T2 − T1

∆r= −296, 10 =⇒ T2 − T1 = −2, 96 =⇒ T2 = 47, 04


Tabela 7.1: Resultados para o problema da aleta finita com fluxo de calor na base conhecido

r T Texata dT/dr2 dT/drexata

0.06 50.00 50.00 -296.10 -296.10

0.07 47.04 47.37 -237.49 -233.42

0.08 44.66 45.28 -192.33 -185.70

0.09 42.74 43.62 -156.45 -148.02

0.10 41.18 42.30 -127.22 -117.38

0.11 39.90 41.26 -102.90 -91.84

0.12 38.88 40.45 -82.27 -70.07

0.13 38.05 39.85 -64.49 -51.15

0.14 37.41 39.42 -48.90 -34.39

0.15 36.92 39.15 -35.01 -19.26

0.16 36.57 39.03 -22.46 -5.37

Utilizando agora a equacao (7.7) para o ponto 2 (r = 0, 07) tem-se:(10000 +

50

0, 07

)T3 − 20093, 02 47, 04 +

(10000− 50

0, 07

)50 = −2325, 6

que implica em:

10714.3T3 = 945175, 7− 464285, 7− 2325, 6 =⇒ T3 = 44, 66

Aplicando-se esta temperatura na equacao do ponto 3 e possıvel calcular a temperatura do

ponto 4 e assim sucessivamente. A tabela (7.1) mostra os resultados para este caso.

Comparando estes resultados com os obtidos pela solucao exata percebe-se um desvio. Este

desvio tende a diminuir com a reducao do espacamento da malha. No entanto e conveniente

ressaltar que atraves de Runge Kutta foram obtidos resultados muito mais precisos. Isto ja

era esperado uma vez que foi utilizado um procedimento de 4a ordem contra o esquema de 2a

usado no caso de diferencas finitas.

Na tabela (7.2) sao mostrados resultados do valor obtido para a temperatura quando r =

0.15 m em diferentes malhas e o ponto nodal onde a derivada se inverte. Obviamente que ha

uma melhora na precisao quando se utiliza malhas mais densas.

Exemplo: Resolva o problema anterior para o caso onde o fluxo de calor que passa pela base

e desconhecido, mas a aleta e infinita.


Solucao: E o mesmo caso do exemplo resolvido para Runge-Kutta na pagina 38. A equacao

diferencial, e consequentemente a sua discretizacao, nao mudam em nada do obtido no exemplo

anterior, assim como a condicao de temperatura prescrita na base. A unica condicao modificada

foi a de fluxo de calor na base que deixou de existir, passando agora para:

∂T

∂r= 0 quando r →∞

Estipulando o valor em que o fluxo de calor como nulo a 1 m da base (r=1,06 m). Desta

forma utilizando ∆r = 0, 1, sera obtido um sistema com 100 incognitas.

Desta forma a nova equacao de contorno fica:

dT

dr

∣∣∣∣r=1.06

=T101 − T100

∆r= 0 =⇒ T101 − T100 = 0

Utilizando-se a equacao (7.7) o sistema de equacoes fica sendo:

T1 = 50 (Condicao de Contorno)

10714, 3T3 − 20093, 02T2 + 9285, 7T1 = −2325, 6

10625T4 − 20093, 02T3 + 9375T2 = −2325, 6

10555, 6T5 − 20093, 02T4 + 9444, 4T3 = −2325, 6...

10047, 6T101 − 20093, 02T100 + 9952, 4T99 = −2325, 6

T101 − T100 = 0 (Condicao de Contorno)

Este sistema de equacoes deve ser resolvido por TDMA ou por Gauss-Seidel. Neste caso

o sistema foi resolvido por TDMA e o resultado para os primeiro 20 pontos e os ultimos 10

encontra-se na tabela (7.3).

Comparando o fluxo de calor dissipado pela base calculado por diferencas finitas com o

calculado por Ruge Kutta, percebe-se novamente uma diferenca significativa (qb = 175, 96 W

Tabela 7.2: Comparacao do caso da aleta com fluxo de calor na base conhecido para diversas

malhas

∆r 0,01 0,005 0.001 0.0001 Exato

T (r = 0.15) 36,92 38,03 38,93 39,13 39,15

x∗i 0.19 0.175 0.1660 0.1640 0.16408

∗ Posicao nodal onde o sinal da derivada se inverte.


Tabela 7.3: Resultados para o problema da aleta infinitar A B C D P Q T dT/dx −kAdT/dx

0.06 1.00 50.00 -0.00000 50.000 50.00 -367.51 175.960.07 -20093 10714.29 9285.71 -2325.50 0.53323 23.223 46.32 -300.00 167.570.08 -20093 10625.00 9375.00 -2325.50 0.70392 14.578 43.32 -248.66 158.740.09 -20093 10555.56 9444.44 -2325.50 0.78510 10.413 40.84 -208.53 149.760.10 -20093 10500.00 9500.00 -2325.50 0.83105 8.014 38.75 -176.48 140.830.11 -20093 10454.55 9545.45 -2325.50 0.85973 6.482 36.99 -150.47 132.080.12 -20093 10416.67 9583.33 -2325.50 0.87875 5.436 35.48 -129.07 123.590.13 -20093 10384.62 9615.38 -2325.50 0.89188 4.689 34.19 -111.28 115.430.14 -20093 10357.14 9642.86 -2325.50 0.90119 4.137 33.08 -96.35 107.630.15 -20093 10333.33 9666.67 -2325.50 0.90791 3.718 32.12 -83.72 100.210.16 -20093 10312.50 9687.50 -2325.50 0.91280 3.394 31.28 -72.99 93.180.17 -20093 10294.12 9705.88 -2325.50 0.91638 3.139 30.55 -63.80 86.550.18 -20093 10277.78 9722.22 -2325.50 0.91899 2.937 29.91 -55.91 80.300.19 -20093 10263.16 9736.84 -2325.50 0.92088 2.775 29.35 -49.09 74.430.20 -20093 10250.00 9750.00 -2325.50 0.92222 2.643 28.86 -43.20 68.940.21 -20093 10238.10 9761.90 -2325.50 0.92315 2.536 28.43 -38.07 63.800.22 -20093 10227.27 9772.73 -2325.50 0.92376 2.449 28.05 -33.61 58.990.23 -20093 10217.39 9782.61 -2325.50 0.92413 2.377 27.71 -29.71 54.520.24 -20093 10208.33 9791.67 -2325.50 0.92431 2.318 27.42 -26.29 50.350.25 -20093 10200.00 9800.00 -2325.50 0.92435 2.269 27.15 -23.30 46.48

......

......

......

......

......

0.97 -20093 10051.55 9948.45 -2325.50 0.91290 2.178 25.00 -0.01 0.060.98 -20093 10051.02 9948.98 -2325.50 0.91285 2.179 25.00 -0.01 0.060.99 -20093 10050.51 9949.49 -2325.50 0.91280 2.180 25.00 -0.01 0.051.00 -20093 10050.00 9950.00 -2325.50 0.91275 2.181 25.00 -0.00 0.041.01 -20093 10049.50 9950.50 -2325.50 0.91271 2.182 25.00 -0.00 0.031.02 -20093 10049.02 9950.98 -2325.50 0.91266 2.184 25.00 -0.00 0.021.03 -20093 10048.54 9951.46 -2325.50 0.91261 2.185 25.00 -0.00 0.021.04 -20093 10048.08 9951.92 -2325.50 0.91257 2.186 25.00 -0.00 0.011.05 -20093 10047.62 9952.38 -2325.50 0.91253 2.187 25.00 0.00 -0.001.06 -1 0.00 1.00 0 0.00000 25.001 25.00

contra q = 165, 31 W, no caso da aproximacao de 4a ordem e 165,86 W da solucao exata. A

unica forma de contornar esta diferenca na precisao e utilizando malhas mais refinadas.

CAPITULO 8

Aplicacao do Metodo de Diferencas Finitas em

Problemas Bidimensionais e Transientes

Foi visto no capıtulo anterior a forma de discretizacao geral atraves de diferencas finitas.

Embora o procedimento em si nao mude tanto os problemas bidimensionais como transientes

apresentam certos detalhes que serao tratados separadamente.

Alem disto estes sao os grupos onde ha uma maior utilizacao deste tipo de discretizacao

uma vez que, nos casos vistos ate agora, o metodo de Runge Kutta poderia ser aplicado com

sucesso. De agora para frente isto nao mais acontece pois as equacoes diferenciais tanto num

caso como no outro sao equacoes diferenciais parciais.

8.1 Solucao de problemas bidimensionais

A solucao de problemas tridimensionais nao envolve grandes dificuldades em termos de

formulacao, o seu unico inconveniente e o fato de que o sistema de equacoes obtidos nao e

Tridiagonal e portanto nao pode ser resolvido por TDMA. Desta forma a unica possibilidade

de solucao para o sistema e o metodo de Gauss Seidel, dentre os vistos ate agora.

No entanto sera apresentado a seguir um procedimento denominado por ADI, sigla em ingles

102


para direcao alternadamente implıcita. Este metodo nada mais e do que uma forma de converter

o sistema de equacoes para um sistema tridiagonal e, portanto, possıvel de ser solucionado pelo

TDMA.

A metodologia e simples e consiste em tratar uma das direcoes implicitamente, enquanto a

outra e tratada de forma explıcita. E mostrado a seguir um exemplo desta forma de solucao.

Exemplo: Resolva o problema o problema de conducao em uma placa plana de alumınio

onde os seus lados estao submetidos a temperaturas especificadas de 100, 200, 100, 200◦C,

rotacionalmente. Considere a placa quadrada com lados de 2 m. Utilizando a discretizacao

mostrada na figura abaixo obtenha o valor das temperaturas nos pontos.

Dados: kalumınio = 215 W/m◦C.

6

-e e e ee e e ee e e ee e e e

1,1 1,2 1,3 1,4

2,1 2,2 2,3 2,4

3,1 3,2 3,3 3,4

4,1 4,2 4,3 4,4

T = 200◦C

T = 100◦C

∂T/∂y = 0

∂T/∂x = 0

Solucao: A equacao de conducao bidimensional ja foi discretizada em exemplo anterior

(pagina 97), e neste caso onde a malha e uniforme foi mostrado que a expressao para a tempe-

ratura generica i, j e dada por:

Ti,j =Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1

4

ou na forma matricial:Ti+1,j

4+Ti,j+1

4− Ti,j

Ti−1,j

4+Ti,j−1

4= 0

Considerando o procedimento ADI e preciso adotar valores iniciais para a solucao Ti,j =

150◦C, para todos os valores de i e j fora das condicoes de contorno.


E importante ressaltar que a presenca de condicoes de simetria no interior da placa e equi-

valente a existencia de uma linha ou coluna identica a do lado oposto. No caso, por exemplo,

da linha de simetria inferior e equivalente a existencia de uma linha identica a linha 2 do lado

inferior a linha 1. Fato similar ocorre na superfıcie direita onde e equivalente a existir do

lado esquerdo da coluna 1 uma coluna identica a 2. Quando nos utilizamos desta condicao de

simetria a expressao geral do problema passa a valer tambem para estes pontos no contorno.

Adotando como primeira direcao a ser resolvida implicitamente a direcao x tem-se que a

equacao para o ponto i, j fica sendo:

Ti+1,j

4− Ti,j +

Ti−1,j

4= −

T ∗i,j+1

4−T ∗

i,j−1

4

onde ∗ indica a utilizacao dos valores admitidos, ou calculados na iteracao anterior. Repare que

neste ponto a matriz ja se tornou uma tridiagonal e pode ser resolvida atraves do algoritmo

TDMA.

No caso deste problema o sistema de equacoes fica:

2T2,1 + 2T1,2 − 4T1,1 = 0

T1,3 + 2T2,2 − 4T1,2 + T1,1 = 0

T1,4 + 2T2,3 − 4T1,3 + T1,2 = 0

T1,4 = 100

2T2,2 + T3,1 − 4T2,1 + T1,1 = 0

T2,3 + T3,2 − 4T2,2 + T2,1 + T1,2 = 0

T2,4 + T3,3 − 4T2,3 + T2,2 + T1,3 = 0

T2,4 = 100

2T3,2 + T4,1 − 4T3,1 + T2,1 = 0

T3,3 + T4,2 − 4T3,2 + T3,1 + T2,2 = 0

T3,4 + T4,3 − 4T3,3 + T3,2 + T2,3 = 0

T3,4 = 100

T4,1 = 200

T4,2 = 200

T4,3 = 200

T4,4 = 200

Adotando a direcao x como a primeira a ser resolvida implicitamente, o primeiro sistema


tridiagonal a ser resolvido fica:

2T1,2 − 4T1,1 = −300

T1,3 − 4T1,2 + T1,1 = −300

T1,4 − 4T1,3 + T1,2 = −300

T1,4 = 100

2T2,2 − 4T2,1 = −300

T2,3 − 4T2,2 + T2,1 = −300

T2,4 − 4T2,3 + T2,2 = −300

T2,4 = 100

2T3,2 − 4T3,1 = −350

T3,3 − 4T3,2 + T3,1 = −350

T3,4 − 4T3,3 + T3,2 = −350

T3,4 = 100

T4,1 = 200

T4,2 = 200

T4,3 = 200

T4,4 = 200

Resolvendo por TDMA obtem-se a seguinte solucao para o campo de temperaturas:

T1,1 T1,2 T1,3 T1,4 T2,1 T2,2 T2,3 T2,4

148,07 146,15 136,53 100 148,07 146,15 136,53 100

T3,1 T3,2 T3,3 T3,4 T4,1 T4,2 T4,3 T4,4

172,12 169,23 154,81 100 200 200 200 200

De posse desta nova solucao e possıvel partir para a proxima etapa tratando explicitamente

os termos na direcao x e implicitamente os da direcao y. Desta forma o novo sistema de equacoes

fica:2T2,1 − 4T1,1 = −292, 3

T3,1 − 4T2,1 + T1,1 = −292.3

T4,1 − 4T3,1 + T2,1 = −338, 5

T4,1 = 200


2T2,2 − 4T1,2 = −284, 6

T3,2 − 4T2,2 + T1,2 = −284, 6

T4,2 − 4T3,2 + T2,2 = −326, 9

T4,2 = 200

2T2,3 − 4T1,3 = −246, 2

T3,3 − 4T2,3 + T1,3 = −246, 2

T4,3 − 4T3,3 + T2,3 = −269, 2

T4,3 = 200

T1,4 = 100

T2,4 = 100

T3,4 = 100

T4,4 = 200

Resolvendo por TDMA obtem-se a seguinte solucao para o campo de temperaturas:

T1,1 T2,1 T3,1 T4,1 T1,2 T2,2 T3,2 T4,2

150 143,84 173,07 200 146,15 150 169,23 200

T1,3 T2,3 T3,3 T4,3 T1,4 T2,4 T3,4 T4,4

126,92 130,76 150 200 100 100 100 200

Resolve-se diversas vezes este mesmo procedimento, e sempre utilizando-se dos valores ob-

tidos do passo anterior, a solucao final sera obtida quando os valores do campo de temperatura

deixem de variar (ou sua variacao atinja um limite admitido como toleravel). As diversas

solucoes obtidas para este caso estao mostradas na tabela (8.1) ate a convergencia.

8.2 Solucao de problemas transientes

Ate agora tratamos de problemas onde o buscado era simplesmente o estado estacionario

(ou regime permanente), mas existe uma parcela muito grande de problemas em engenharia em

que o interesse maior recai sobre o tempo necessario ate que o regime estacionario (ou algum

outro fenomeno) ocorra. Nestes casos e preciso considerar o regime transiente do problema.

O estudo do transiente em si nao apresenta grandes problemas, ele apenas inclui novos

termos relacionados a variacao da variavel em relacao ao tempo na equacao diferencial. O

problema passa a ser, desta forma, uma composicao de diversas solucoes intermediarias que

funcionam como fotos do problema em cada instante estudado.


Tabela 8.1: Valores das iteracoes para o caso de conducao numa placa plana

Posi- Valores Iteracao 1 Iteracao 2 Iteracao 3 Iteracao 4

cao Iniciais Linha Coluna Linha Coluna Linha Coluna Linha Coluna

1,1 150.00 148.08 150.00 148.82 150.00 150.27 150.00 149.94 150.00

1,2 150.00 146.15 146.15 143.79 143.79 144.33 144.33 144.21 144.21

1,3 150.00 136.54 126.92 126.33 126.92 127.06 126.92 126.89 126.92

1,4 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

2,1 150.00 148.08 153.85 156.21 156.21 155.67 155.67 155.79 155.79

2,2 150.00 146.15 150.00 150.89 150.00 149.80 150.00 150.05 150.00

2,3 150.00 136.54 130.77 131.95 131.95 131.68 131.68 131.74 131.74

2,4 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

3,1 150.00 172.12 173.08 172.49 173.08 173.21 173.08 173.05 173.08

3,2 150.00 169.23 169.23 168.05 168.05 168.32 168.32 168.26 168.26

3,3 150.00 154.81 150.00 149.70 150.00 150.07 150.00 149.98 150.00

3,4 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00

4,1 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00

4,2 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00

4,3 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00

4,4 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00 200.00

Alem de tudo isto os problemas transientes trazem junto consigo uma serie de novas de-

finicoes na bagagem, de acordo com a forma que e feita a sua discretizacao. A seguir apresen-

taremos algumas das mais importantes.

8.2.1 Formulacao Explıcita

A formulacao explıcita e aquela que mais se aproveita do carater parabolico da solucao

transiente, como ja foi discutido em capıtulos anteriores. Ela se utiliza de uma solucao num

instante anterior ja conhecida para calcular a solucao posterior desejada. Utilizando-se deste

procedimento ela nao depende da inversao de matrizes tornando-se extremamente facil e rapida

para ser aplicada, pois o a obtencao dos resultados depende unica e exclusivamente das tempe-

raturas calculadas para o instante anterior.

No entanto, esta facilidade tem o seu custo: este tipo de formulacao nao e incondicionalmente

estavel. Isto implica em que ha determinadas condicoes que devem ser satisfeitas para que se

consiga obter uma solucao para problema. No caso, por exemplo da equacao de calor, a condicao


a ser satisfeita e:

α∆t

∆x2≤ 1

2

onde ∆t e o incremento de tempo, ∆x o menor valor de espacamento da malha e α a difusividade

termica.

Exemplo: Considere uma placa de aco carbono com 1,6 cm de espessura que inicialmente se

encontra ha uma temperatura uniforme de 600◦C. A placa num instante t = 0, num banho de

15◦C. O coeficiente transmissao de calor e adotado igual a 10.000 W/m2◦C. As propriedades do

aco carbono podem ser adotadas como constantes (k = 40 W/m◦C e α = 1×10−5m2/s). Calcule

o tempo necessario para que a temperatura do centro da placa atinja 100◦C e a temperatura

de um plano situado a 0,2 cm de uma das superfıcies da chapa neste instante.

Solucao: A equacao diferencial para o problema transiente de conducao unidimensional e:

∂T

∂t= α

∂2T

∂x2

que se discretizada na forma explıcita fica:

Tm+1i − Tm

i

∆t= α

Tmi+1 − 2Tm

i + Tmi−1

∆x2

No caso da formulacao explıcita o unico termo desconhecido e o do instante m + 1 e desta

forma:

Tm+1i =

(1− 2

α∆t

∆x2

)Tm

i +α∆t

∆x2Tm

i+1 +α∆t

∆x2Tm

i−1

Adotando um espacamento de malha regular e igual a 0,2 cm obtem-se 5 pontos no domınio

como mostra a figura (8.1), e neste caso o valor maximo para o incremento de tempo e:

α∆t

∆x2= 1× 10−5 ∆t

(2× 10−3)2=

∆t

0, 4≤ 1

2=⇒ ∆t ≤ 0, 2

Adotando o valor maximo possıvel para o ∆t = 0, 2 s a equacao geral fica sendo:

Tm+1i =

1

2Tm

i+1 +1

2Tm

i−1

No caso deste problema temos condicao de simetria no centro da placa. Podemos adotar

neste ponto tambem a origem do eixo de coordenadas, vide figura (8.1). Isto implica que a

condicao de contorno nesta face e dada por:

∂T

∂x= 0 =⇒ Tm+1

2 − Tm+11

∆x= 0 =⇒ Tm+1

1 = Tm+12


-u u u u ue1 2 3 4 52

Figura 8.1: Malha na placa considerando a simetria.

A outra condicao e dada na superfıcie e implica que o fluxo de calor que sai por conveccao

e igual ao que chega por conducao:

−k A ∂T

∂x

∣∣∣∣p

= hA (Tp − T∞)

que pode ser divida por A e discretizada na forma:

Tm+15 − Tm+1

4

∆x=h

k(T∞ − Tm+1

5 ) →(

1 +h∆x

k

)Tm+1

5 = Tm+14 +

h∆x

kT∞

ou substituindo os valores (h = 1× 104 W/m2◦C, ∆x = 2× 10−3 m e k = 40 W/m◦C):

1, 5Tm+15 = Tm+1

4 + 0, 5T∞ −→ Tm+15 =

Tm+14 + 0, 5T∞

1, 5

Desta forma de posse da equacao geral para os pontos internos e das duas condicoes de

contorno, pode-se partir dos valores iniciais e marchar rumo a solucao desejada. Neste caso a

temperatura da placa e inicialmente uniforme em 600◦C, podemos calcular o campo de tempe-

raturas em t = 0, 2 s:

T 12 =

T 01 + T 0

3

2=

600 + 600

2= 600

T 13 =

T 02 + T 0

4

2=

600 + 600

2= 600

T 14 =

T 03 + T 0

5

2=

600 + 600

2= 600

e aplicando as condicoes de contorno:

T 11 = T 1

2 = 600

T 15 =

T 14 + 0, 5T∞

1, 5=

600 + 0, 5 15

1, 5= 405

A tabela (8.2) mostra os resultados obtidos para as iteracoes. Para economia de espaco a

partir de t = 2.0 s sao mostrados apenas resultados em intervalos de 1 s.


Tabela 8.2: Transiente numa placa plana.

t T1 T2 T3 T4 T5

0.0 600.00 600.00 600.00 600.00 600.00

0.2 600.00 600.00 600.00 600.00 405.00

0.4 600.00 600.00 600.00 502.50 340.00

0.6 600.00 600.00 551.25 470.00 318.33

0.8 575.63 575.63 535.00 434.79 294.86

1.0 555.31 555.31 505.21 414.93 281.62

1.2 530.26 530.26 485.12 393.41 267.28

1.4 507.69 507.69 461.84 376.20 255.80

1.6 484.76 484.76 441.94 358.82 244.21

1.8 463.35 463.35 421.79 343.08 233.72

2.0 442.57 442.57 403.22 327.76 223.50

3.0 352.59 352.59 321.40 261.91 179.61

4.0 281.50 281.50 256.88 209.91 144.94

5.0 225.38 225.38 205.94 168.87 117.58

6.0 181.08 181.08 165.73 136.47 95.98

7.0 146.10 146.10 133.99 110.89 78.92

8.0 118.50 118.50 108.93 90.70 65.46

8.8 100.66 100.66 92.74 77.65 56.77

9.0 96.70 96.70 89.15 74.76 54.84

10.0 79.50 79.50 73.54 62.17 46.45

11.0 65.91 65.91 61.21 52.24 39.83

12.0 55.19 55.19 51.48 44.40 34.60

Muito embora tenha sido obtido o resultado ele esta ainda um pouco afastado da solucao

analıtica que preve que a temperatura do centro alcancara os 100◦C em 11,6 s (contra os 8,8 s

da solucao numerica). Neste mesmo instante esperava-se obter uma temperatura de 73,6◦C a

0,2 cm da superfıcie contra os 77,65◦C obtidos numericamente. Detalhes da solucao analıtica

podem ser vistos em (Bejan, 1996), pagina 137. Embora este erro tenha valores significativos,

principalmente no que se refere ao tempo e possıvel minimizar os seus valores com o refinamento

da malha. Vejamos por exemplo as solucoes obtidas para diversas malhas.


Malha Tempo para Temperatura a 0,2 cm

(∆t×∆x) o centro atingir 100◦C da superfıcie da placa

0, 2× 0, 2 8,8 s 77,65 ◦C

0, 05× 0, 2 8,85 s 76,98◦C

0, 05× 0, 1 10,2 s 75,10◦C

0, 0125× 0, 05 10,86 s 74,53◦C

8.2.2 Formulacao Totalmente Implıcita

Ao contrario da formulacao explıcita, este tipo se utiliza das temperaturas no proprio

instante de tempo a ser calculado para avaliar as derivadas. Desta forma o grau de complexidade

do problema aumenta consideravelmente, uma vez que a solucao do campo de temperaturas so

pode ser obtida atraves da solucao de um sistema linear de equacoes.

Este aumento na complexidade da solucao do problema, no entanto, apresenta em con-

traponto uma vantagem, o esquema torna-se totalmente estavel. Desta forma o incremento

de tempo utilizado se torna independente da malha, permitindo uma melhor otimizacao dos

parametros de solucao.

8.2.3 Formulacao Implıcita

Este tipo de formulacao e qualquer uma em que estejam envolvidos termos na posicao

de tempo atual. Vejamos por exemplo se utilizamos uma aproximacao generica para o campo

de temperaturas onde:

T θi = θ Tm+1

i + (1− θ)Tmi (8.1)

temos que atraves desta expressao generica que se:

θ = 0 o esquema recai na formulacao explıcita;

θ = 1 o esquema recai na formulacao totalmente implıcita;

0 < θ < 1 o esquema recai na formulacao implıcita.

E obvio que a escolha do valor de θ influira nas caracterısticas do sistema a ser solucionado.

Um grande exemplo disto e a estabilidade, a formulacao implıcita so e incondicionalmente

estavel para valores de θ ≥ 1/2.


ee eee e m+ 1

m

e

e

e

e

e

ee

e

e

e

e

e

m+ 1

m+ 1

m

m

Esquema Totalmente Implıcito

Esquema Implıcito

Esquema Explıcito

6k3

�-

6

6k3�-

Figura 8.2: Interrelacao entre os pontos em cada esquema de discretizacao de transientes.

Existem diversos valores de θ utilizados comumente na solucao de problemas transientes

e, de acordo com a escolha do valor de θ o esquema recebe um nome diferente. De todos

estes esquemas o mais difundido e tambem, que apresenta melhores resultados, e o de Crank-

Nicholson, obtido quando se toma θ = 1/2. Este esquema e muito utilizado pois funciona como

um tipo de aproximacao de segunda ordem em relacao ao tempo, e e altamente recomendado

quando se esta estudando problemas transientes onde se deseja acompanhar o tempo de evolucao

de um dado fenomeno.

A figura (8.2) mostra esquemas dos valores utilizados em cada iteracao pelos metodos

explıcito, implıcito e totalmente implıcito. Esta figura e bastante ilustrativa e permite en-

tender o relacionamento existente entre os pontos em um dado instante para cada um dos

esquemas.

Exemplo: Considere uma placa de aco carbono com 1,6 cm de espessura que inicialmente se

encontra ha uma temperatura uniforme de 600◦C. A placa num instante t = 0, num banho de

15◦C. O coeficiente transmissao de calor e adotado igual a 10.000 W/m2◦C. As propriedades do

aco carbono podem ser adotadas como constantes (k = 40 W/m◦C e α = 1×10−5m2/s. Calcule

o tempo necessario para que a temperatura do centro da placa atinja 100◦C e a temperatura

de um plano situado a 0,2 cm de uma das superfıcies da chapa neste instante.


Solucao: Este problema e identico ao resolvido anteriormente utilizando-se da formulacao

explıcita, no entanto agora nos utilizaremos da discretizacao generica, equacao (8.1):

∂T

∂t= α

∂2T

∂x2

que se discretizada na forma generica fica:

Tm+1i − Tm

i

∆t= α

Tm+θi+1 − 2Tm+θ

i + Tm+θi−1

∆x2

Adotando a formulacao de Crank-Nicholson, temos que θ = 1/2:

Tm+1i − Tm

i

∆t= α

Tm+1/2i+1 − 2T

m+1/2i + T

m+1/2i−1

∆x2

e as temperaturas podem ser admitidas como:

Tm+1/2 =1

2Tm +

1

2Tm+1

Substituindo os valores na expressao e reagrupando os termos de forma a deixar os no tempo

m+ 1 no primeiro termo da equacao:(1 +

α∆t

∆x2

)Tm+1

i −− α∆t

2 ∆x2Tm+1

i+1 − α∆t

2 ∆x2Tm+1

i−1 =(1− α∆t

∆x2

)Tm

i +α∆t

2 ∆x2Tm

i+1 +α∆t

2 ∆x2Tm

i−1

Adotando um espacamento de malha regular e igual a 0,2 cm obtem-se 5 pontos no domınio

e um incremento de tempo de 0,2 s a equacao geral fica sendo:

1, 5Tm+1i − 0, 25Tm+1

i+1 − 0, 25Tm+1i−1 = 0, 5Tm

i + 0, 25Tmi+1 + 0, 25Tm

i−1

No caso deste problema temos condicao de simetria no centro da placa. Podemos adotar

neste ponto tambem a origem do eixo de coordenadas, vide figura (8.1). Isto implica que a

condicao de contorno nesta face e dada pela existencia de um ponto identico ao ponto 4:

Tm+11 = Tm+1

2

que se substituirmos na equacao geral para o problema:

1, 5Tm+11 − 0, 5Tm+1

2 = 0, 5Tm1 + 0, 5Tm

2


A outra condicao e dada na superfıcie e implica que o fluxo de calor que sai por conveccao

e igual ao que chega por conducao:

−k A ∂T

∂x

∣∣∣∣p

= hA (Tp − T∞)

que pode ser divida por A e discretizada na forma:

Tm+15 − Tm+1

4

∆x=h

k(T∞ − Tm+1

5 ) →(

1 +h∆x

k

)Tm+1

5 = Tm+14 +

h∆x

kT∞

ou substituindo os valores (h = 1× 104 W/m2◦C, ∆x = 2× 10−3 m e k = 40 W/m◦C):

1, 5Tm+15 − Tm+1

4 = 0, 5T∞

Desta forma e possıvel montar o sistema de equacoes para o primeiro passo no tempo

considerando-se as duas condicoes de contorno e o valor inicial uniforme para a temperatura

de 600◦C:1, 5Tm+1

1 − 0, 5Tm+12 = 600

1, 5Tm+12 − 0, 25Tm+1

3 − 0, 25Tm+11 = 600

1, 5Tm+13 − 0, 25Tm+1

4 − 0, 25Tm+12 = 600

1, 5Tm+14 − 0, 25Tm+1

5 − 0, 25Tm+13 = 600

1.5Tm+15 − Tm+1

4 = 7, 5

A tabela (8.3) mostra os resultados obtidos para as iteracoes. Para economia de espaco a

partir de t = 0, 1 s sao mostrados apenas alguns resultados.

Relembrando do exercıcio anterior onde se viu que o tempo para que a temperatura do

centro atinja 100◦C e a temperatura a 0,2 cm da superfıcie no mesmo instante, obtidos atraves

de metodos analıticos, sao respectivamente 11,6 s e 73,6 s. Comparando-se com os resultados

obtidos atraves deste metodo percebe-se que estes se aproximam bem mais da solucao exata

que o metodo explıcito (11 s e 73,1◦C). Estes resultados podem ser ainda melhorados com o

refinamento da malha, veja por exemplo se reduzirmo o ∆t para 0,1 s e a malha para ∆x =

0, 4 mm obtem-se o tempo para chegar a 100◦Cigual a 11,5 s e a a 0,2 cm da superfıcie de

73.3◦C. Refinamentos maiores fornecerao resultados melhores.

A tabela (8.4) mostra a evolucao transiente da temperatura na posicao central da placa

utilizando-se ∆t = 0, 2 s e ∆x = 2 mm utilizando-se dos tres esquemas: explıcito, totalmente

implıcito e Crank-Nicholson. Esta tabela permite uma comparacao de como funciona cada um

dos esquemas.


Tabela 8.3: Solucao por Crank-Nicholson do transiente numa placa plana.

t T1 T2 T3 T4 T5

0.0 600.00 600.00 600.00 600.00 600.00

0.2 599.62 598.86 593.51 562.22 379.81

0.4 597.39 593.70 573.95 501.88 339.59

0.6 591.38 583.07 548.30 463.24 313.83

0.8 580.99 568.53 524.34 434.61 294.74

1.0 567.18 552.03 502.60 411.73 279.49

1.2 551.28 534.63 482.68 392.48 266.65

1.4 534.29 516.95 464.18 375.65 255.43

1.6 516.86 499.34 446.80 360.51 245.34

1.8 499.40 482.01 430.34 346.59 236.06

2.0 482.17 465.10 414.67 333.61 227.41

3.0 402.81 388.23 345.54 277.86 190.24

4.0 336.21 324.10 288.66 232.57 160.05

5.0 280.99 270.95 241.60 195.15 135.10

6.0 235.25 226.94 202.63 164.17 114.45

7.0 197.38 190.50 170.37 138.52 97.35

8.0 166.02 160.32 143.65 117.28 83.19

9.0 140.05 135.33 121.53 99.69 71.46

10.0 118.54 114.64 103.21 85.13 61.75

11.0 100.74 97.50 88.04 73.07 53.71

12.0 86.00 83.32 75.48 63.08 47.06


Tabela 8.4: Comparacao da solucao do transiente para uma placa plana utilizando-se de tres

metodos para discretizacao no tempo.

Esquemas

Tempo Explıcito Tot. Implıcito Crank-Nicholson

(s) θ = 0 θ = 1 θ = 1/2

0.00 600.00 600.00 600.00

0.20 600.00 597.55 599.62

0.40 600.00 591.59 597.39

0.60 600.00 582.24 591.38

0.80 575.63 570.15 580.99

1.00 555.31 556.08 567.18

1.20 530.26 540.74 551.28

1.40 507.69 524.64 534.29

1.60 484.76 508.21 516.86

1.80 463.35 491.73 499.40

2.00 442.57 475.40 482.17

3.00 352.59 399.25 402.81

4.00 281.50 334.55 336.21

5.00 225.38 280.56 280.99

6.00 181.08 235.67 235.25

7.00 146.10 198.37 197.38

8.00 118.50 167.37 166.02

9.00 96.70 141.61 140.05

10.00 79.50 120.20 118.54

11.00 65.91 102.42 100.74

12.00 55.19 87.64 86.00

REFERENCIAS BIBLIOGRAFICAS

Bejan, A., Transferencia de Calor, Editora Edgard Blucher, 1996.

Dhatt, G. e Touzot, G., The Finite Element Method Displayed, Ed. J. Wiley & Sons, 1984.

Domingues, M. O. and Mendes, Jr, O., Introducao aos Programas Ci-

entıficos de Distribuicao Gratuita: GNU/Octave, GNU:Maxima, LaTeX

e GNU/Rcs, Instituto Nacional de Pesquisas Espaciais - INPE, 2002,

http://mtc-m16.sid.inpe.br/col/sid.inpe.br/iris\%401916/2005/08.16.14.21/

doc/Ermac2002.pdf, Acesso em 28/09/2006.

Eaton. J.W., Octave Manual, http://www.gnu.org/software/octave/doc/interpreter/,

2006, Acesso 25/09/2006.

Holman, J.P., Transferencia de Calor, Editora McGraw Hill, 1983.

Incropera, F.P. e DeWitt,D.P., Fundamentos da Transferencia de Calor e Massa, 4a ed., LTC

Ed. 1998.

Kreith, F. Princıpios da Transmissao de Calor, Editora Edgard Blucher, 1977.

Lewis, R.W., Morgan,K., Thomas, H.R. e Seetharamu, K.N.,The Finite Element Method in

Heat Transfer Analysis, Ed. J. Wiley & Sons, 1996.

Maliska, C. R., Transferencia de Calor e Mecanica dos Fluidos Computacional, LTC Editora,

1995.

Patankar, S., Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980.

117


Ruggiero, M. A. G. e Lopes, V.L.R., Calculo Numerico: Aspectos Teoricos e Computacionais,

McGraw Hill, 1988.

Apostila de Transferencia de Calor Computacional, UNICAMP, 1989.

Documents

Métodos Numéricos em Fluido Térmica - feb.unesp.br · a) E.D.O. - Sistema Massa Mola b) E.D.P. - Eq. da quantidade de movimento Figura 1.1: Exemplos de E.D. ordinárias e parciais