Apostila Métodos Numéricos - UTFPR

U NIVERSIDADE T ECNOLÓGICA F EDERAL DO P ARANÁ D EPARTAMENTO ACADÊMICO DE M ECÂNICA C URSO URSO DE E NGENHARIA M ECÂNICA

MÉTODOS NUMÉRICOS APLICADOS À ENGENHARIA ENGENHARIA INTRODUÇÃO AOS MÉTODOS DE DIFERENÇAS FINITAS E DE VOLUMES FINITOS

Prof. Dr. Admilson T. Franco

Agosto de 2011

Métodos Numéricos Aplicados Aplicados à Engenharia

CAPÍTULO 1 Introdução ao Método de Diferenças Finitas e Aplicações Nesse capítulo são apresentadas as motivações para o estudo de métodos numéricos em engenharia e exemplos de aplicações da solução numérica. 1.1 Introdução

Equações diferenciais ordinárias (EDOs) e parciais (EDPs) aparecem em inúmeros problemas da física-matemática. Em especial, na área de engenharia, todo cálculo um pouco mais elaborado normalmente recai em uma equação diferencial. Como poucas EDs têm solução analítica possível ou viável, os métodos numéricos aparecem como uma ferramenta extremamente eficiente para sua solução. O uso de métodos numéricos é extremamente difundido em engenharia, atualmente. A Tabela 1.1 mostra o vertiginoso desenvolvimento desenvolvimento do hardware nos últimos últimos 50 anos. Pode-se também afirmar que o desenvolvimento dos softwares seguiu a mesma proporção. Tabela 1.1 Comparação entre o Eniac e o Pentium. Pentium. Volume Peso Arquitetura do Sistema Operações por segundo Preço (US$)

Eniac (1.945) 558 m³ 30 ton. 17.468 válvulas 3,5 mil 260 milhões

Pentium (1.993) cabe em qualquer mesa 5 kg 3,1 milhões de transistores 100 milhões 3 mil

100 mil vezes mais rápido que o Eniac. Fonte: Revista ISTOÉ / 1.377 - 21/02/96 pág. 84. Pentium.Pro  100

Na Figura 1.1 é apresentada uma foto do primeiro computador, o ENIAC.

2


CAPÍTULO 1 Introdução ao Método de Diferenças Finitas e Aplicações Nesse capítulo são apresentadas as motivações para o estudo de métodos numéricos em engenharia e exemplos de aplicações da solução numérica. 1.1 Introdução

Equações diferenciais ordinárias (EDOs) e parciais (EDPs) aparecem em inúmeros problemas da física-matemática. Em especial, na área de engenharia, todo cálculo um pouco mais elaborado normalmente recai em uma equação diferencial. Como poucas EDs têm solução analítica possível ou viável, os métodos numéricos aparecem como uma ferramenta extremamente eficiente para sua solução. O uso de métodos numéricos é extremamente difundido em engenharia, atualmente. A Tabela 1.1 mostra o vertiginoso desenvolvimento desenvolvimento do hardware nos últimos últimos 50 anos. Pode-se também afirmar que o desenvolvimento dos softwares seguiu a mesma proporção. Tabela 1.1 Comparação entre o Eniac e o Pentium. Pentium. Volume Peso Arquitetura do Sistema Operações por segundo Preço (US$)

Eniac (1.945) 558 m³ 30 ton. 17.468 válvulas 3,5 mil 260 milhões

Pentium (1.993) cabe em qualquer mesa 5 kg 3,1 milhões de transistores 100 milhões 3 mil

100 mil vezes mais rápido que o Eniac. Fonte: Revista ISTOÉ / 1.377 - 21/02/96 pág. 84. Pentium.Pro  100

Na Figura 1.1 é apresentada uma foto do primeiro computador, o ENIAC.

2


Figura 1.1 Primeiro computador – ENIAC (1945) O uso de técnicas numéricas para solução de complexos problemas de engenharia engenharia tem permitido o projeto e otimização de equipamentos e sistemas de uma forma jamais imaginada a 30 ou 40 anos atrás. Isto foi possível principalmente ao fantástico desenvolvimento da capacidade computacional tanto em termos de velocidade e capacidade de armazenamento como ao aperfeiçoamento dos métodos numéricos. A facilidade de utilização dos métodos numéricos e a qualidade dos resultados obtidos tem sido um atrativo sempre crescente para aumentar sua utilização, aliada a economia de tempo de projeto e obviamente, o bviamente, do custo total do equipamento. Para se reforçar a atratividade dos métodos numéricos para simulação, atualmente um escoamento turbulento supersônico sobre um aerofólio requer alguns minutos de CPU com custo de centenas de dólares. A mesma simulação na década de 60, consumiria um tempo de computação de aproximadamente 30 anos, com custo de 10 milhões de dólares. A enorme difusão da simulação numérica se justifica, então, pelo barateamento dos equipamentos computacionais e desenvolvimento (aprimoramento) dos métodos numéricos utilizados na soluções das equações diferenciais que regem o fenômeno. 1.2 Métodos de Solução de Problemas

Em engenharia existem principalmente principalmente três métodos de solução de problemas: 

Métodos analíticos, os quais conduzem a resultados analíticos, (R. A.);



Métodos numéricos (experimentação (experimentação numérica), (R. N.); N.);



Experimentação em laboratório, (R.E.). (R.E.).

3

Métodos Numéricos Aplicados à Engenharia

Os métodos analíticos são a melhor forma de solucionar os problemas, pois fornecem um solução de forma fechada. Entretanto, poucos problemas de engenharia podem ser resolvidos dessa forma, devido as dificuldades impostas pelo conjunto de equações (normalmente EDPs de 2a ordem e não-lineares) que regem o fenômeno. Mesmo as soluções analíticas para alguns problemas, quando existem, normalmente contém séries infinitas, funções especiais ( erf,



), equações transcendentais para auto-

valores, etc. No entanto, as soluções analíticas, mesmo de alguns problemas simples e de pouco interesse prático, servem de base para a compreensão do comportamento do sistema de equações, para o desenvolvimento de métodos numéricos e validação de códigos computacionais. A experimentação em laboratório trata com a configuração real e é imprescindível quando estudando um novo fenômeno. Tem como desvantagem o custo muitas vezes proibitivo, problemas de segurança como transferência de calor no núcleo de reatores nucleares, ou impossibilidade de reprodução das condições reais, como vôos supersônicos em grandes altitudes ou simulação de reservatórios de petróleo. Na ausência de modelos matemáticos estabelecidos e em geometrias extremamente complexas, muitas vezes esta é a única alternativa viável. Os métodos numéricos praticamente não apresentam restrições e oferecem atrativas vantagens como: 



Baixo custo: a maior vantagem sobre os outros métodos; Velocidade: centenas de diferentes configurações podem ser testadas em

poucas horas; 

Informações Completas: fornece o valor das variáveis relevantes em

qualquer ponto de interesse; 

Facilidade de Simular Condições Realísticas : pode tratar qualquer

condição de contorno, velocidades altas ou baixas, temperaturas altas ou baixas, domínios pequenos ou amplos. Qualquer geometria arbitrária pode, a princípio, ser tratada. 1.3 Tipos de Erros

Dois tipos de erros que podem estar presentes na solução numérica quando os resultados são comparados com a realidade de um problema físico: 

erros numéricos: resultantes da má solução das equações diferenciais. É necessário

comparar com outras soluções analíticas ou numéricas. 4




erros: resultantes do uso de equações diferenciais que não representam

adequadamente o fenômeno. A comparação entre: (R.N.) e (R.A.) ou outro (R.N.)  validação numérica. (R.N.) e (R.E.)  validação física. 1.4 Solução Numérica

A solução numérica de um problema de engenharia é feita com as seguintes etapas: 1. Modelo matemático: construído a partir da observação do fenômeno, usando leis da física e da matemática. Julga-se que esse modelo represente muito apropriadamente o comportamento real do fenômeno físico. 2. Modelo numérico: obtido a partir do modelo matemático usando-se um método de aproximação para as equações diferenciais. Esse método consiste na discretização do domínio e na solução da equação em pontos específicos. Portanto, tem-se a seguinte sequência em uma solução numérica: 

Modelo Físico



(Processo iterativo)

Modelo Matemático  Resultados



Modelo Numérico



Solução Numérica

Numéricos

A Figura 1.2 ilustra de forma adequada essas etapas.

Modelagem Matemática

Ajuste do Modelo

Problema Físico

Equações Governantes

Discretização

Sistema de Equações Algébricas

Análise e Interpretação

Solução Aproximada

Resolução das Equações Algébricas

5


Figura 1.2 Etapas envolvidas na solução numérica de um problema físico Notar que, como desejamos que os resultados numéricos representem o fenômeno físico que estamos reproduzindo numericamente, deve-se ter o máximo cuidado em todas as etapas de solução de um problema, com o risco de se obter resultados que não possuem significado físico. Uma importante observação é: “A ferramenta numérica é a adequada e confiável quando se está de posse de um método numérico que resolva corretamente as equações diferenciais, e de um modelo matemático, que represente com fidelidade o fenômeno físico.” 1.5 Diferenças Finitas, Volumes Finitos e Elementos Finitos.

O Método das Diferenças Finitas (MDF) sempre foi empregado pelos analistas das área de escoamento de fluidos, enquanto o Método dos Elementos Finitos (MEF) o foi para a área estrutural, na solução de problemas de elasticidade. Os problemas de escoamento são altamente não-lineares, enquanto os da elasticidade não possuem os termos convectivos, não-lineares, e assemelham-se a problemas puramente difusivos de Transferência de Calor. O MDF, atualmente, pode ser utilizado em qualquer tipo de malha, mesmo a nãoestruturada usada em elementos finitos. O Método dos Volumes Finitos (MVF) garante a conservação da propriedade envolvida (massa, quantidade de movimento e entalpia) no volume elementar, e é preferível em relação ao MDF. Na atualidade, ambos os métodos (MVF e MEF) estão resolvendo problemas altamente convectivos, inclusive com ondas de choque e em geometrias arbitrárias. Isto é esperado, visto que, eles são derivados do mesmo princípio e diferem apenas na forma de minimização escolhida. O MEF, aplicado em níveis de volumes elementares produz o Volume Finite Element Method (CVFEM). O Método dos Elementos de Contorno (BEM - Boundary Element Method), possibilita tratar apenas com a discretização da fronteira, sem necessidade de discretizar o domínio interno. Na Figura 1.3 são mostrados alguns exemplos de aplicação dos métodos numéricos em problemas de engenharia.

6


AEROESPACIAL

Figura 1.3 Exemplos de Aplicação Aplicação dos Métodos Numéricos

7


Atividade 1

Realizar a simulação do processo de transferência de calor em uma placa com as condições de contorno mostradas. Determinar o campo de temperatura em regime permanente. Fazer a simulação numérica do problema de condução de calor proposto com o programa Trans-CalV1.1, usando as 11 etapas do Processo de Análise e explicando-as. Entregar na próxima semana. O programa Trans-Cal V1.1 pode ser baixado do link: http://www.sinmec.ufsc.br/sinmec/site/download.html

y

H =0,5 =0,5 m

T t = 100oC

T l = 400oC

o T r r=300 C

T b= 200oC L=1,0 m

8


Lista de Exercícios 1

1) Quais são os métodos usados em engenharia para solução de problemas? Quais as características de cada um e o que diferencia a escolha para determinada aplicação ? 2) No seu ponto de vista, quais as duas maiores vantagens da utilização de métodos numéricos para problemas de engenharia? 3) Explique as etapas para a solução numérica de um problema físico. 4) Foram obtidos dois conjuntos de dados sobre um certo problema físico. Um conjunto originou-se de uma simulação numérica, o outro a partir de um experimento físico. Qual dos conjuntos você considera mais representativo do fenômeno? Por quê? Discuta quais fatores podem levar um conjunto a ser mais representativo do que o outro, e vice-versa. 5) O que justifica o emprego em larga escala de métodos numéricos para problemas de engenharia atualmente? 6) Cite outros exemplos de aplicação de métodos numéricos não apresentados na Figura F igura 1.1.

7) Cite 2 programas comerciais de larga aplicação em dinâmica dos fluidos. Faça o mesmo para a análise estrutural ?

9


CA P ÍT U L O 2

A s p e ct o s M a t em át i c o s d a s E q u a çõe s, N ív ei s d e F o r m u l a ção d o s M o d elo s, E q u a ções D i f er en c i a i s P a r c i a i s Este capítulo apresenta os principais conceitos referentes aos aspectos matemáticos das equações diferenciais parciais (EDPs). Primeiramente, são definidos os níveis de formulação dos modelos para a solução dos problemas físicos. Na seqüência, discute-se os três tipos básicos de EDPs, salientando-se a diferença entre eles, tanto do ponto de vista físico como matemático.

2.1 Níveis de Formulação dos Modelos.

A obtenção da solução de qualquer problema físico requer a habilidade da criação do modelo matemático correspondente. O modelo matemático deve ser tal que possa ser resolvido com tempos de computação não proibitivos e que os resultados obtidos representem corretamente o fenômeno físico em questão. O sistema de equações gerado pelo modelo matemático pode ser obtido tanto a nível molecular, ou seja, uma equação para cada molécula, como a nível de v olumes de controle, onde se supõe um valor médio da propriedade dentro desse volume. Para o caso de escoamento de fluido newtoniano, as equações de conservação da massa, da quantidade de movimento e da energia, para sistema car tesiano de coordenadas são:

  t   t

  ui  

  t



 x j

  x   0

 u u    x

  x j

i

(2.1)

j

P

 x j

  T  





j



i



  u T    x j

j

   ui    x j   x j

 u   S 

 k   T  T    S  c p   x j 

i

(2.2)

(2.3)

10



sendo  a massa especifica, t a variável do tempo x j , u j representam a componente do sistema de coordenada e do vetor velocidade nas direções j  1,2,3 respectivamente,  a viscosidade dinâmica, T a temperatura, k o coeficiente de condutividade térmica do material e c p o calor especifico.

A Figura 2.1 apresenta de forma sucinta os três métodos de solução de problemas de engenharia, com os respectivos requisitos associados ao emprego de cada método.

11


Métodos Numéricos Aplicados à Engenharia e d a . c d i t t n e , a a u í q g r o e m n o e , c a , s o ã s ç a a m v , r o e t s n n e o m c i e v d o s m i e e L d 

e d s o l e d o m , s a v i t u t . i t c s t n e o , c i a s c e n ê õ l ç u a b l e r u R t 

S O C I R Ó E T S O D O T É M

s i s a r e o d p a i m d e r t a o e e ã n i ç l s i a l a o i a o c ã p h r a n l p s s a t e o s a t n e d n m i s o e a e d d e t õ n m a a ç e l z s e a e r m p r õ a o g t u ç . n c e a c a t t a u t r n e I T N F E

o n r o t n o c e d s e õ ç i d n o C























S S O O C I D R O É T É M U M N

N R

O C I O L T Á E M D E O T M A M

A R

O C I S Í F A M E L B O R P

S S O O C I D T Í O T L É A N M A

S I A S T N O E D M O T I É R E M P X E

s o d a s i c o a o s n i d v d o ê t d o o g a r o h h o r e n p e ã s n t v i a a ç e r n m u l a m m e o s o o e t a t a t c l s n c e e i i o l o e d d c d s d d o s s a l o o o a a h i d m h l a l a r s v h r é r o o e . t t c l o t e i e c s s a c é i r v t s t n C i d E M s E m E i

a t s a s i a x e i e õ c n ç o s a e ã d a u r ç e q f u l e i d o S

I E O M R D E Ó A T S D A E A T R C S O N B A E T B A L o o t t n n e e e o m d m s i ã r a a o ç e d s d p p e i s a e x d r a e d c e l i c s . s n i c m o o o e i r o t d s C L P d E 







S O D A T L U S E R

s o c i t í l a n a s o d a t l u s e R

s i a t n e m i r e p x e s o d a t l u s e R

s o c i r é m u n s o d a t l u s e R

A R

E R

N R

Figura 2.1 Métodos de Solução (Maliska,1995)

12



As equações (2.1) (2.2) e (2.3) podem ser escritas para um campo escalar geral f :

 t

    

 

  x



 u    

  y

  v  



 

 z

  w   

  

(2.4)

       S  x   x  y  y  z  z  







Os significados de f ,  e S  para cada equação são dados na Tabela 2.1. O fechamento do problema é obtido com a equação de estado:

    P , T 

(2.5)

O problema tridimensional representado pelas Eqs. (2.4) e (2.5), possui 6 equações e 6 incógnitas: r , u , v , w, P e T .

2.2 Equações Diferenciais Parciais (EDPs).

Os métodos numéricos tratam da obtenção de soluções numéricas para equações diferenciais parciais ou ordinárias. Para podermos aplicar técnicas computacionais a problemas que envolvem essas equações, é importante primeiro identificar suas características gerais. As EDPs que descrevem os fenômenos de interesse da dinâmica de fluidos computacional podem ser classificadas em três categorias básicas. 1. Elípticas. 2. Parabólicas. 3. Hiperbólicas. Essa classificação não é meramente acadêmica, uma vez que cada classe de equações está associada a uma categoria diferente de fenômenos físicos. Além disso, métodos numéricos que funcionam para uma categoria de equações podem n ão funcionar para as demais.

13



Tabela 2.1 Valores de  ,  e S 

Equação de conservação

Massa global

Quantidade de movimento em x

Quantidade de movimento em y

Quantidade de movimento em z

Energia *

Massa de um componente

*

é o termo de dissipação viscosa.

Uma equação diferencial pode ser classificada sob diferentes aspectos. A ordem de uma equação diferencial parcial é a ordem da maior derivada parcial presente na equação. Considerese a equação diferencial de segunda ordem em duas variáveis x e y , que não necessariamente representam coordenadas espaciais: A

 2  x

2

B

 2  x y

C

 2  y

2

D

 x

E



 F  G

 y

(2.6)

14



Quando A, B, C , D, E , F e G forem constantes ou funções de x , y apenas, a equação (2.6) é dita linear. Em função dos valores de A, B e C , pode-se classificar a equação (2.6) nos três tipos mencionados acima:



Elíptica, se   B 2  4 AC  0



Parabólica, se   B 2  4 AC  0



Hiperbólica, se   B 2  4 AC  0

Equações de ordem superior a dois podem ser classificadas nas três damílias anteriores, se forem reescritas como um sistema de equações de primeira ordem. Uma discussão mais detalhada sobre a classificação de EDPs pode-ser encontrada em Fletcher (1992).

15



Tabela 2.2 Níveis de Formulação dos Modelos (Maliska,1995)

Nivel em que os balanços de conservação são efetuados

Conservação para cada molécula

Balanços onde:

Informações necessárias

Tipo de equação resultante

Massa molecular, leis de troca de QM, campos de forças: eléricos, magnéticos, etc.

Equação para cada molécula

Propriedades refletindo o comportamento molecular

Conjunto de equações diferenciais parciais

etc.

Balanços onde:

Fornecer

etc. e as

tensões de Reynolds, relações de transferência de calor e massa turbulenta

Balanços onde o volume de controle coincide com o domínio de solução em alguma(s) direção(ões)

Fornecer as condições de contorno nas direções onde o volume de controle coincide com o domínio de solução

Conjunto de equações diferenciais parciais

Equações diferenciais, parciais, ordinárias ou algébricas

tempo médio sobre os quais os balanços de conservação são realizados tempo entre colisões moleculares escala de tempo para turbulência comprimento médio sobre os quais os balanços de conservação são realizados livre caminho médio entre as moléculas escala de comprimento para turbulência

16



2.2.1 Problemas de Equilíbrio

Problemas de equilíbrio são aqueles nos quais a propriedade de interesse não se altera como o passar do tempo. Matematicamente, esses problemas são, em geral, representados por equações diferenciais parciais elípticas, cuja equação modelo é a equação de Laplace. Em coordenadas cartesianas, essa equação pode ser escrita como:

   2

 2  x 2



 2  y 2

0

(2.7)

sendo  a variável dependente e 2 o operador laplaciano que, em coordenadas cartesianas bidimensionais, é dado por

  2



2

 x 2



 2  y 2

(2.8)

A solução única para problemas envolvendo equações diferenciais parciais é obtida especificando-se condições sobre a variável dependente na fronteira  R da região R , em que se quer resolver o problema. Problemas que exigem condições ao longo da fronteira (contorno)  R de toda a região são denominados de problemas de valor de co ntorno (PVC). Ver Figura 2.2.

 R

R

Figura 2.2 Região R com fronteira  R na qual resolvemos problemas de equilíbrio. As condições sobre a variável dependente são aplicadas ao longo de toda a fronteira da região.

Uma característica dos problemas regidos por equações elípticas é que toda a região R é imediatamente afetada por qualquer mudança no valor da variável dependente em um ponto P no interior de R , ou em sua fronteira  R (ver Figura 2.3). Isso equivale a dizer que perturbações deslocam-se em todas as direções dentro de R , afetando todos os demais pontos internos, embora essa influência diminua como o aumento da distância ao ponto P . Portanto, as soluções numéricas de problemas de equilíbrio variam suavemente em R . No caso da chapa, qualquer mudança de temperatura em P altera a distribuição estacionária de temperatura em todos os pontos da região R , sendo as mudanças menores conforme nos afastamos de P . 17



Fisicamente, observamos esse efeito quando colocamos ao fogo uma panela que está inicialmente em equilíbrio térmico. Após a panela atingir novamente o equilíbrio térmico, toda a panela está quente e não só a parte em contato com a ch ama.

Elíptico

P x

Parabólico

P x

Hiperbólico

P x

Figura 2.3 Problemas Elípticos, Parabólicos e Hiperbólicos O método numérico empregado para resolver o problema de equilíbrio deve levar em conta que cada ponto da fronteira  R afeta a solução global.

2.2.2 Problemas Transientes. Os problemas transientes, ou de propagação, envolvem a variação temporal das grandezas físicas de interesse. A partir dos valores iniciais dessas grandezas em um certo tempo t 0 , calculam-se pela solução numérica da EDP (quando não se dispõe da solução analítica), seus novos valores em sucessivos intervalos de tempo t , até alcançarmos o instante final t f . t0  t , t0  2t , t0  3t ,..., t f   t, t f

18



Note que caminha-se (marcha) na dimensão temporal em “passos” de comprimento t , desde t 0 até t f . Por isso, problemas transientes são também denominados “problemas em marcha”. Os fenômenos transientes são modelados por equações diferenciais parabólicas ou hiperbólicas. Quando apresentam mecanismos de dissipação de energia (por exemplo, na difusão de calor e no escoamento de fluidos viscosos), os fenômenos ditos dissipativos são descritos por equações parabólicas. Caso contrário, são representados por equações hiperbólicas.

Equações Parabólicas A equação modelo para problemas parabólicos é a equação transiente de difusão de calor:

1 T



 t

 2T

(2.9)

 x 2

na qual T é a temperatura e  é o coeficiente de difusividade térmica do material. Como exemplo, considere-se a barra delgada de comprimento 10 cm, mostrada na Figura 2.4, isolada termicamente ao longo de seu comprimento, mas que tem suas extremidades mantidas às temperaturas de T0  T1  50C .

T 0

T 1

Figura 2.4 Barra termicamente isolada

Supondo a barra inicialmente a temperatura de T  0C , deseja-se conhecer sua evolução temporal ao longo da barra. A Figura 2.5 mostra o gráfico da temperatura da barra em diferentes instantes de tempo. A partir do estado inicial, passa-se por uma sucessão de distribuições de temperaturas transientes. Após um tempo suficientemente longo, o equilíbrio térmico (estado estacionário) é alcançado e a temperatura da barra não mais varia. Obtém-se, então, uma distribuição uniforme de temperatura.

19


Métodos Numéricos Aplicados à Engenharia 60

Estacionária 50

)

40

C o (

a r u t a 30 r e p m e T 20

Transientes

10

0 0

1

2

3

4

5

6

7

8

9

10

x (cm)

Figura 2.5 Distribuição da temperatura em uma barra para diferentes instantes de tempo Para se poder estudar a evolução temporal da temperatura é necessário que, o valor inicial da temperatura ao longo de toda a barra seja especificado. Com essa informação, e sabendo que a temperatura nos extremos da barra é mantida a 50C (condição de fronteira), obtem-se a distribuição de temperatura ao longo da barra para diferentes instantes de tempo. As condições de fronteira, combinadas (quando apropriado) com condições iniciais, são denominadas condições auxiliares. Ao contrário dos problemas de equilíbrio, os problemas transientes necessitam de valores para a variável dependente  em t  0 (condições iniciais), além de condições de fronteira para t  0 . Problemas desse tipo são denominados problemas de valor inicial (PVI). A necessidade de condições iniciais fica clara quando observa-se que a equação (2.9) só relaciona variações espaciais e temporais de T , mas não especifica o valor de T . Deve-se especificar o valor inicial de T para, sabendo sua variação temporal e espacial (dadas pela equação (2.9)), obter os novos valores de T para t  0 . Pode-se, então escrever:

Valor inicial de T

+

Variação espacial e temporal de T dada pela equação diferencial ,

+

Condições de fronteira

=

Novo valor de T

Essa relação é válida para os problemas de valor inicial em geral. Finalmente, nota-se que, a princípio, é possível marchar indefinidamente no tempo, isto é, calcular a temperatura da barra em qualquer instante de tempo t  0 . Portanto, ao contrário 20



daquela dos problemas elípticos, a região R de solução de problemas parabólicos é aberta. (ver Figura 2.3)

Equações Hiperbólicas Equações hiperbólicas estão relacionadas a problemas de vibração ou de convecção, em que os fenômenos dissipativos são mínimos ou podem ser desprezados. Situações descritas por equações hiperbólicas necessitam tanto de condições iniciais como, em geral, de condições de fronteira (ver Figura 2.3). Portanto, são também problemas de valor inicial. E da mesma forma em que problemas parabólicos, a região R da solução é aberta. A equação modelo do problema hiperbólico é a equação de convecção que, para uma dimensão especial, é escrita como:

 t

 v

  x

(2.10)

Essa equação representa o transporte de  para a direita ao longo de x com velocidade 2 2 v  0 . Como essa equação não possui um termo dissipativo    x , o valor de  deve ser

apenas transportado ao longo de x , sem alteração, entre os instantes t 0 e t0  t (ver Figura 2.6). O produto v

 

(2.11)

 x

é denominado termo convectivo ou inercial .

v

Figura 2.6 Perfil senoidal transportado sem dissipação com velocidade v  0 por uma equação hiperbólica A ausência de mecanismos dissipativos faz com que quaisquer descontinuidades presentes nas condições iniciais se propaguem para a solução em t  0 . Isso implica que as equações hiperbólicas admitem soluções descontínuas, e que o método numérico utilizado deve ser capaz de lidar com elas. 21



Um processo de convecção que possua mecanismos dissipativos obedece à equação parabólica de convecção-difusão, também denominada equação de transporte.

 t

 v

  x



2  x 2

(2.12)

em que v é a velocidade de propagação de  e  é o coeficiente de difusão de  . Para ilustrar os efeitos que a presença ou ausência de mecanismos de dissipação têm sobre a solução numérica de problemas de convecção, a Figura 2.7 mostra a propagação com dissipação do perfil senoidal da Figura 2.6. Os efeitos dissipativos são claramente visíveis: o perfil senoidal foi “espalhado” (difundido) para pontos adjacentes e teve sua amplitude reduzida.

v

Figura 2.7 Perfil senoidal transportado com dissipação. Em relação ao perfil inicial, note a redução da amplitude e difusão do perfil para pontos adjacentes. Finalmente, convém mencionar outra equação hiperbólica importante, a equação da onda de segunda ordem:

 2  t2

c

2

2  x 2

(2.13)

Essa equação descreve, por exemplo, o deslocamento transversal de uma corda sob tensão. A constante c é a velocidade de propagação da onda na corda. Uma expressão semelhante modela a vibração de uma chapa de espessura desprezível, fixa pelas extremidades. A Tabela 2.3 adaptada de Versteeg & Malalasekera (1995), resume a classificação dos fenômenos físicos e as equações vistas até aqui.

Tabela 2.3 Sumário das características dos problemas físicos 22


Métodos Numéricos Aplicados à Engenharia Problema

Tipo de

Equação

equação

Modelo

Condições

Região

Admite solução descontínua?

Equilíbrio Transientes com dissipação Transientes sem dissipação

Elíptica

 2  0

Fronteira

Parabólica

   2 t

Fronteira e

Hiperbólica

   v t x

Fronteira e

iniciais

iniciais

Fechada

Não

Aberta

Não

Aberta

Sim

23




1) No contexto de métodos numéricos para a solução de EDPs (Equações Diferenciais Parciais), qual a importância de se conhecer a categoria (classificação) a que pertence uma dada EDP?

2) Explique as diferenças entre problemas elípticos, parabólicos e hiperbólicos. Exemplifique.

3) Verifique se as funções abaixo são soluções da equação de Laplace:

 2  2  0  x 2  y 2 (a)   xy (b)   x 2  y 2 (c)   ln( x 2  y 2 )

4) Suponha que as funções  1 e  2 sejam duas soluções da equação de Laplace. Mostre que, se c1 e c2 forem constantes, então: 3  c11  c2 1 também é solução.

u u  2u u   2 5) Considere a Equação de Burger dada por: t  x  x onde u é o componente x da velocidade, t é o tempo, x é a coordenadas espacial e  a viscosidade cinemática do fluido. Discuta a natureza desta equação.

6) Considere a seguinte equação diferencial:

 2T  2T 0 x 2   x  y 2 Esta equação pode ser elíptica, parabólica ou hiperbólica, dependendo se x  0 , x  0 ou x  0 . Explique quando está equação é elíptica, parabólica ou hiperbólica.

24

CA P ÍT U L O 3

OMé t o d o d a s D i f e r en ça s F i n i t a s Nesse capítulo é apresentado o Método de Diferenças Finitas e sua aplicação a problemas de engenharia, principalmente em de difusão de calor em regime transitório. 3.1 Processo de discretização

A solução analítica das equações diferenciais parciais envolve expressões na forma fechada as quais fornecem a variação das variáveis dependentes continuamente através do domínio. Em contraste, soluções numéricas podem dar respostas somente em pontos discretos do domínio, chamados pontos nodais ou pontos da malha. Para a solução numérica, as derivadas da função existentes na equação diferencial devem ser substituídas pelos valores discretos da função, como ilustrado na Figura 3.1. A maneira de obter essas equações algébricas é que caracteriza o tipo do método numérico. Lembramos que métodos como o BEM (Boundary Element Method), necessitam apenas da discretização da fronteira do domínio. Tais métodos são denominados meshless.

Equação diferencial

Sistema de equações

L(  ) = 0 e

algébricas [A] [  ] = [B]

condições de contorno Figura 3.1 Tarefa do método numérico. onde L(  ) é o operador diferencial. 25



Por conveniência, vamos assumir que o espaçamento dos pontos da grade, denominada de malha, na direção x é uniforme e igual a  x e que o espaçamento na direção y é também uniforme e dado por  y . Na figura abaixo, os pontos da malha são identificados na direção x pelo índice i e na direção y pelo índice j

y

x

Figura 3.2 Pontos discretos da malha.

3.2 Base do Método de Diferenças Finitas.

Representações das derivadas em diferenças finitas são baseadas na expansão em série de Taylor. Por exemplo, se ui,j denota o componente x da velocidade no ponto ( i, j ), então a velocidade ui+1,j no ponto (i+1,j ) pode ser expressa em termos de uma expansão em série de Taylor sobre o ponto (i, j ):

ui 1, j

  u    2 u   x 2  ui , j     x   2    x  i , j   x i , j 2

  3 u   x 3   3     x i, j 6

(3.1)

A equação (3.1) é matematicamente uma expressão exata se: (a) o número de termos é infinito e a série converge; (b) e / ou  x  0.

26



Para computação numérica, é impraticável calcular um número infinito de termos. Dessa maneira, a equação (3.1) é truncada. Por exemplo, se os termo s de magnitude (  x )3 e termos de ordens mais altas forem desprezados, a equação (3.1) reduz-se a: ui 1  ui , j

  u    2 u   x 2     x   2    x i , j   x i, j 2

(3.2)

Nós dizemos que a equação acima tem acurácia de segunda ordem, porque os termos iguais ou maiores que (  x )3 foram desprezados. Se os termos de ordem (  x )2 e de ordens maiores forem desprezados, nós obtemos: ui 1  ui , j

  u      x   x i , j

(3.3)

onde a equação (3.3) tem acurácia de primeira ordem. Nas equações (3.2) e (3.3), os termos desprezados representam o erro de truncamento na representação da série finita. Por exemplo, o erro de truncamento para a equação (3.2) é:

  n u   x n  n   n! n  3   x  i , j 

e o erro de trucamento para a equação (3.3) é:

  n u   x n  n   n! n  2   x  i , j 

O erro de trucamento pode ser reduzido por: (a) Tomando-se mais termos na série de Taylor, equação (3.1). Isto leva o aumento na ordem de acurácia na representação de ui+1, j ; (b) Reduzindo-se a magnitude de  x .

  u   :   x i, j

Vamos retornar a equação (3.1) e calcular  u u   u     i 1, j i, j  x   x i , j

  2u   x   3u   x 2   2    3   2 6  x  x  i , j  i, j erro de truncamento

ou u u   u     i 1, j i , j  O x   x   x i , j

(3.4)

27



 x ) - termos da ordem de  x .

O (

Derivada a jusante, com acurácia de primeira ordem: u u   u     i 1, j i , j  O x   x   x i , j

(3.5)

Vamos escrever a expansão em série de Taylor sobre o ponto (i-1,j) para u i-1,j ui 1, j

  u    2 u    x 2  ui , j      x    2    x  i , j   x  i , j 2

  u    u   x    u   x   u i , j     x    2    3     x  i , j   x  i , j 2   x  i, j 6 2

ui 1, j

  3 u    x 3   3   , ou 6  x   i , j 2

3

(3.6)

3

  u   , nós obtemos:  x   i, j

Resolvendo para 

u u   u     i , j i 1, j  O x   x   x i , j

(3.7)

Derivada a montante, com acurácia de primeira ordem. Vamos subtrair a equação (3.6) de (3.1): ui 1, j

 ui 1, j

  u    3 u   x 3  2   x   3   3  x  x  i , j  i , j

(3.8)

  u   , obtemos:  x   i, j

Resolvendo para 

u u   u     i 1, j i 1, j  O x 2 2 x   x i , j

(3.9)

Esta é a derivada central, com acurácia de segunda ordem, para o ponto ( i, j ). Para obter uma expressão em diferenças finitas para a derivada parcial de segunda

  2 u   , primeiro reescrevemos a equação (3.8). 2   x  

ordem 

28


Métodos Numéricos Aplicados à Engenharia u u   u     i 1, j i 1, j 2 x   x i , j

  3u   x 2   3     x i , j 6

(3.10)

Substituindo a equação (3.10) em (3.1)

 ui 1, j  ui 1, j ui 1, j  ui , j    2 x

   3u   x 2   2u   x 2   3   x   2   x x 6 2     i , j  i , j 

(3.11)

  u   x   u   x  3    4     x i, j 6   x i , j 24 3

3

4

4

  2 u  Resolvendo (3.11) para  2  , obtemos   x  i , j u  2ui , j  ui 1, j   2u   2   i 1, j  O x 2 2  x    x i , j

(3.12)

Derivada central de segunda ordem. Expressões em diferenças para as derivadas em y são obtidas da mesma maneira. u u   u     i , j 1 i , j  O y   y   y i , j

Derivada a jusante

u u   u     i , j i, j 1  O y   y   y i , j

Derivada a montante

u u   u     i, j 1 i , j 1  O y 2 2 y   y  i , j

Derivada Central

A derivada segunda também pode ser calculada de seguinte forma:

  2u   2    x i , j

   u    u      x     x       u   i 1, j  i 1, j  1         x  x   x   x  i , j    

 u  u   u  u  1   2u   2    i 1, j i, j    i , j i 1, j    x i , j   x    x   x u  2ui, j  ui 1, j   2u   2   i 1, j  x 2   x i, j

(3.13)

29



A equação (3.13) é análoga a equação (3.12).

  2u   no ponto (i, A mesma filosofia pode ser usada para se obter a derivada mista    x y i , j j ). Por exemplo:

  u    u           2u     u    y i 1, j   y i 1, j         x y 2 x  i , j  x   y 

(3.14)

 u u u   2u    u  1     i 1, j 1 i 1, j 1    i 1, j 1 i 1, j 1  2 y 2 y   x y i , j     2 x ou

  2u  1    ui 1, j 1  ui 1, j 1  ui 1, j 1  ui 1, j 1  O x2 ,  y 2    x y i, j 4 x y

(3.15)

A Figura 3 mostra a interpretação geométrica das derivadas.

Figura 3.3 Interpretação geométrica das fórmulas de diferenças finitas para as derivadas de 1a. ordem.

Exercício 3.1 – Dada a função f  x  =

1 4

x 2 , calcule a primeira derivada de f em x =2,0 usando

aproximação em diferenças finitas avançada e retrógrada. Compare esses resultados com uma 30



aproximação em diferenças centrais e o valor analítico exato. Use tamanhos de passos:  x = 0,1 e 0,4. Preencha o quadro abaixo.

 x = 0,1

Valor

 x = 0,4

 0,1

Analítico

 0,4

Der. Retrógrada Der. Avançada Der. Central

Exercício 3.2 – Dada a função f  x   sin2 x  , determine

 f  x

em x =0,375, usando

representação em diferença central de ordem (  x )2 e (  x )4 . Use passos 0,01, 0,1 e 0,25. Compare e discuta os resultados. Preencha o quadro abaixo.  f  x



f i 1  f i 1

2 x 

 O x 2

Aproximação por diferenças centrais

ordem (  x )2 :

 f

ordem (  x )4 :

 f

 x

 x

 

f i 1  f i 1

2 x 

 O x 2

 f i  2  8 f i 1  8 f i 1  f i  2  O x 4 12 x   Valor aproximado  Valor analitico   100 (%) Valor   analitico

determinação do erro percentual:   

 x = 0,01

 x = 0,1

 x = 0,25

Valor Analítico

 0,01

 0,1

 0,25

ordem (  x )2 ordem (  x )4

3.3 Tratamento da fronteira.

31



Muitas outras aproximações diferentes podem ser obtidas para as derivadas acima, como bem para as derivadas de ordens maiores. O que acontece com a fronteira? Que tipo de difereciamento é possível quando nós temos somente uma direção? Imagine o problema abaixo.

Figura 3.4 Pontos da malha na fronteira

  u 

 na Nós desejamos construir uma aproximação por diferenças finitas para   y  1 fronteira. A aproximação mais fácil é a derivada a jusante:

  u  u 2  u1     O y   y   y 1

(3.16)

a qual tem acurácia de primeira ordem. Como podemos obter uma acurácia de segunda ordem? A derivada central falha, porque não existe o ponto 2’. Vamos assumir que na fronteira u passa ser expressa pelo polinômio: u(y)= a +by +cy 2

(3.17)

Aplicando a equação (3.17) para os pontos da malha.

u1 = a  2 u2 = a +b y + c   y   2 u3 = a +b  2  y  + c  2  y 

(3.18)

Resolvendo este sistema para b : b

 3u1  4u2  u3 2 y

(3.19)

32



Retornando a equação (3.17) e diferenciando:  u  y

 b  2cy

(3.20)

Se (3.20) for avaliada na fronteira onde y =0,

  u     b   y 1

(3.21)

Combinando (3.19) e (3.20) temos:

  u   3u1  4u2  u3    2 y   y 1

(3.22)

Vamos verificar a ordem de aproximação da eq. (3.20) através da expansão em série de Taylor sobre o ponto 1.

  u    2 u  y 2   3u  y 3  y   2  u y   u1     3   y  y 2  1  1   y 1 6

(3.23)

Comparando as equações (3.23) e (3.17), percebemos que a expressão polinomial assumida é análoga aos primeiros três termos da série de Taylor. Portanto, a equação (3.17) é de O(  y )3 . Na formação da derivada (3.22) nós dividimos por

 y , o que então faz (3.22) da

ordem de O(  y )2 . Portanto:

  u   3u1  4u2  u3     O y 2 2 y   y 1

(3.24)

Exercício 3.3 – Considere a aproximação polinomial de 2 a. ordem.

Figura 3.5 Função polinomial passando por três pontos.

33



a) Calcular a função   x  pela expressão polinomial:   x   a  bx  cx 2

Solução de   x  :   x  

 x  x2  x  x3   x  x1  x  x3   x  x1  x  x2   1  2    x1  x2  x1  x3   x2  x1  x2  x3   x3  x1  x3  x2  3

Assuma que ( x 2 - x 1)  ( x 3 - x 2) b) Calcular:

        x  x

c) Calcular:

              ,   ,   e   x 1   x  2   x 3

  2    2  , assumindo   x  2

 x   x2  x1    x3  x2   cte

3.4 Formulações Explícita, Totalmente Implícita e Implícita.

Nessa seção serão expostos os princípios das formulações: explícita, totalmente implícita e implícita. Vamos exemplificar utilizando o problema de condução transiente unidimensional, que pode ser expresso por:  (  cPT )  t



  T  k   q  x   x  

(3.25)

sendo T a temperatura [ oC ],  a massa específica do material [ kg/m 3] , cP o calor específico à pressão constante [ J/(kg K)], k a condutividade térmica do material [ W/(m K)] e q o termo fonte de geração de calor [ W/m 3]. Para propriedades constantes a equação (3.25) se reduz para: 1T   t



 2 T  x

2



q k

(3.26)

sendo   k /(  cP ) a difusividade térmica do material [ m2/s]. A Figura 3.6 ilustra a notação empregada para a malha numérica em problemas unidimensionais. P denota o ponto nodal da malha sendo resolvido, e os seus respectivos nós vizinhos à direita, E , e à esquerda, W . 34



Figura 3.6 Malha numérica unidimensional. As aproximações espaciais são obtidas discretizando-se as derivadas com o Método de Diferenças Finitas,

 T  x

 x

 x

 x

2 P

Derivada de primeira ordem a montante

 T W 2  O x 2 x

Derivada de primeira ordem central

T E



T E



T E

P

 2 T

 T W  O x x

 P

 T

Derivada de primeira ordem a jusante

T E

P

 T

 T P  O   x x



 2TP  TW 2  O x 2   x 

Derivada de segunda ordem central

(3.27)

A aproximação do termo transitório pode ser feita como:

 T  t

 P

n 1

TP

 T Pn  O   t   t

(3.28)

sendo: T P  temperatura no ponto P no instante futuro ( n  1 ). 0 T P  temperatura no ponto P no instante atual ( n ).

As temperaturas no tempo anterior são conhecidas em todos os pontos do domínio. Uma questão importante surge agora com relação ao nível de tempo no qual será avaliado o lado direito da equação (3.28). Lembrando que este termo representa os fluxos difusivos de calor, e que estamos avançando a solução de um nível de tempo para outro, devemos decidir se vamos avaliar esses fluxos no início, no fim ou em uma posição intermediária do intervalo de tempo. Denotando por  a posição, no intervalo de tempo, de avaliação do termo difusivo, temos a seguinte aproximação numérica para (3.28).

35


Métodos Numéricos Aplicados à Engenharia n 1

TP

 TPn  t

T

 E

 2T P  T  W q     O[( t ,( x) 2 ] 2 ( x) k

(3.29)

3.4.1 Formulação Explícita 1D (  n) .

Escolhendo  = n, temos a formulação explícita, onde todas as temperaturas vizinhas a P são avaliadas no instante anterior e, portanto, já são co nhecidas. É possível explicitar a incógnita da equação ( T Pn1 ) em função de temperaturas vizinhas, todas conhecidas. A formulação explícita dá origem a um conjunto de equações algébricas que podem ser resolvidas uma a uma, obtendo-se a temperatura em cada ponto do espaço para o novo nível de tempo. Pelo fato das equações não serem acopladas entre si, não existe a necessidade de resolver uma matriz.

Figura 3.7 Malha regular unidimensional Para a formulação explícita, isto é,  = n, e a malha unidimensional da figura 3.7, a equação (3.29) pode ser reescrita como: n 1

TP

 2T Pn  T Wn q P   (  x) 2 k

(3.30)

 1  2rx  TPn  rx (TEn  TWn )  QP

(3.31)

 TPn  t

T

n E

ou, n 1

TP

onde r x

r x



 t

 x 2

QP



qP k

 t

(3.32)

 Fo =número de Fourier 36



Considere o problema unidimensional de transferência de calor por condução, sem termo-fonte, e com as seguintes condições: *

Condições iniciais: para t  0 s T2  0; T3  0; T 4  0

*

Condições de contorno: para t  0 s T1  100 ; T 5  0 Para a formulação explícita, a evolução no tempo exige a limitação do intervalo de

tempo, ou seja, um critério de estabilidade. Da análise de Von Neumann, na Eq. (3.31): r x 

1 2

,

para manter o coeficiente de T pn sempre positivo ou nulo. Se isto não for atendido, o coeficiente (1  2r x ) na equação (3.31) pode oscilar entre valores positivos e negativos, e isto ferir princípios termodinâmicos. Baseado na Eq, (3.31) e empregando r x 

1 , as expressões para T 2 , T 3 e 2

T 4 são

( QP  0 ):

T n 1  1 T n  T n  3 1 2 2  T n 1 1 T n T n  3  4 2 2  T n 1  1 T n  T n  5 3  4 2 Resultando, em regime permanente, no campo de temperatura: {T 1=100, T 2=75 , T 3=50, T 4=25 e T 5 =0 } A evolução do perfil de temperaturas com o tempo é mostrado na figura abaixo.

Figura 3.8 Evolução do perfil de temperatura. 37


Métodos Numéricos Aplicados à Engenharia Exemplo 3.1 – Condução de calor 1D – MDF – Formulação explícita (FE).

a) Determine o campo de temperatura em regime permanente para a barra isolada lateralmente e com termo fonte de geração uniforme de calor, q¢¢¢ = 5 ´ 103 W / m 3 . Inicialmente a barra está a uma temperatura isotérmica T i = 30oC , quando então são impostas as condições de contorno mostradas. Use o Método de Diferenças Finitas e 3 pontos internos, como mostrado na figura. O material da barra é o alumínio com propriedades: condutividade térmica térmica a = 1, 00 ´ 10k = 177 W / (m .K ) e difusividade

4

m2 / s .

b) Determine a solução analítica do problema na forma literal. Compare com os resultados do item anterior.

T 2

T 3

T 4

o

o

T 5 =20 C

T 1=80 C

 x

x L=0,26 m

3.4.2 Formulação Totalmente Implícita 1D (  n  1) .

Se na equação (3.29) fizemos  = n+1, teremos a formulação totalmente implícita.

 2T nP1 T Wn 1 qP   (  x )2 k

(3.33)

 r xTEn1  (1  2 rx ) TPn1  rxTWn 1  TPn  QP

(3.34)

n 1

TP

 TPn  t

T

n 1 E

ou,

onde as definições de r x e de QP são as mesmas apresentadas na Eq. (3.32). A Eq. (3.34) quando aplicada para cada ponto da malha, origina um sistema de equações algébricas, uma vez que elas estão acopladas entre si. Não existe mais possibilidade de coeficiente negativo para T P . A formulação é incondicionalmente estável e o intervalo de tempo é limitado por precisão. Para o problema discutido na seção 3.4.1 e com r x 

1 2

e QP  0 :

38



T n1  1 T n1  T n 1  1 T n 3 1  2 2 2 4  T n1  1 T n1  T n 1  1 T n 4  2 3 3 2 4   n1 1 n1 n 1 1 n T4  4 T5  T3   2 T4

(3.35)

As equações acima podem ser escritas na forma: aPT P

 a E T E  aW T W  b

(3.36)

Fazendo P =2, 3 e 4:

 a2 T2n 1  a3T3n1  b2  n 1 n 1 n1  a3 T3  a4 T4  a2 T2  b3  a T n 1  a T n1  b 3 3 4  4 4

(3.37)

onde as temperaturas T 1 e T 5 foram incluídas nos termos independentes b2 e b4, respectivamente. O sistema de equações (3.37) pode ser escrito na forma matricial como:

 AT    B

(3.38)

 1 T  T n  0  T  b2   2 1 2  rx (1  2r x )  r  T n1   b    T n  rx (1  2rx)  x    3n1   3   3    0 (1  2 rx )  T4  b4   1 r x n  T5  T 4  2  n 1 2

Sendo:

(3.39)

[A]  a matriz de coeficientes (constante). [T]  vetor de incógnitas. [B]  vetor lado direito (dependente do tempo). O sistema de equações (3.37) deve ser resolvido para cada intervalo de tempo, pois o

problema é transitório. Exemplo 3.2 – Condução de calor 1D – MDF – Formulação totalmente implícita (FTI).

Resolver o Exemplo 3.1, empregando agora na solução a formulação totalmente implícita.

39


Métodos Numéricos Aplicados à Engenharia 3.4.3 Formulação Implícita.

Na formulação implícita, os valores das temperaturas que entram no cálculo do fluxo de calor (difusivo) são tomados como uma ponderação dos valores dessas temperaturas no começo e no fim do intervalo de tempo. O método mais conhecido é o de Crank-Nicolson   1 / 2 , onde a temperatura é tomada como uma média aritmética entre as temperaturas T pn e T pn1 , como: n 

TP

n 1

TP

 1    T Pn

(3.40)

 Temperatura conhecida no instante presente (n).

T Pn T Pn



1

 Temperatura a ser determinada no instante futuro (n+ 1).

Figura 3.X Função de interpolação no tempo. É importante observar que basta  ser diferente de zero para que as equações fiquem acopladas, caracterizando a implicitude entre as mesmas. A figura 3.9 ilustra, para os três tipos de formulações, as conexões existentes entre o ponto P e seus vizinhos, no instante de tempo futuro e no instante de tempo presente.   n (Formulação Explícita)

n    n  1 (Formulação Implícita)

40



  n  1 (Formulação Totalmente Implícita)

Figura 3.9 Tipos de formulação Da equação (3.40), para  = n+1/2, obtemos: TPn1 2 n1 2

TP



 1

1 2

n1 P

T 2

 

TPn1  1 

1

n

 T p 2

 T Pn 

(3.41)

Na equação (3.29) : n 1

TP

 TPn  t

T

n 1/ 2 E

 2T nP1/ 2 T Wn1/ 2 q    ( x) 2 k

(3.42)

Substituindo na equação (3.42), expressões semelhantes a (3.41) para T nP1/ 2 , T E n1/ 2 e n1/ 2 T W , obtemos:

 r xTEn1  2(1  rx )TPn 1  rxTWn 1  2(1  rx )TPn  rx (TEn  TWn )  QP

(3.43)

Observe que a equação resultante para a formulação implícita é nesse caso uma ponderação entre as equações (3.31) para a formulação explícita e a (3.34) para a formulação totalmente implícita. As definições de r x e de QP são as mesmas apresentadas na equação (3.32). 41



Discuta qual o critério de estabilidade necessário para a equação (3.43). 3.5 Condições de contorno

As condições de contorno podem ser de três tipos: a. de 1ª espécie ou de Dirichlet; b. de 2ª espécie ou de Neuman; c. de 3ª espécie ou de Robin. 3.5.1 Condição de contorno de 1ª espécie ou de Dirichlet

Para essa condição de contorno, o valor da variável  deve ser especificado ao longo de uma fronteira  R , de um domínio R .  s

R

 R

Figura 3.10 Domínio de solução de uma equação diferencial. Para o caso de condução de calor, essa condição de contorno corresponde a ter a temperatura da fronteira especificada: T f . T f

T P

...  x

L

3.5.2 Condição de contorno de 2ª espécie ou de Neuman

Para essa condição de contorno, o gradiente normal:  /   f , ou tangencial:

 / s  g , deve ser especificado na fronteira. Quando

f  0 , dizemos que se trata de uma

condição de fronteira natural ou homogênea. No caso de condução de calor, o fluxo de calor é especificado na fronteira: q  k

T .  x

42


Métodos Numéricos Aplicados à Engenharia T f

q

T P

...  x

q   k T f

L

(T  T ) (T  T ) T   k P f  k f P  x x x qx   T P k

Uma aproximação de ordem O(x)2 , pode ser obtida para  T  x



3T1  4T2  T3 2 x

(3.44)

T , através da equação (3.22):  x

(3.45a)

Resultando para a temperatura na fronteira T f : T f



2 qx 4 1  TP  T P 1 3 k 3 3

(3.45b)

3.5.3 Condição de contorno de 3ª espécie ou de Robin

Nesse caso, a condição de contorno é uma combinação linear das condições de 1ª e de 2ª espécies. 

    h 

(3.45)

Sendo  e  valores constantes. Essa condição corresponde para condução de calor, a se ter os valores do coeficiente de película h e da temperatura T  da corrente convectiva, especificados. T f

T P

... h, T 

 x

L

O calor transportado por convecção tem que ser igual ao fluxo de calor que atravessa a

  h(T  Tf ) e o fluxo de calor fronteira. Igualando na fronteira o fluxo de calor convectivo: qconV TP  T f

  k difusivo: qconD

 x

: 43



T f



hT  (k / x)TP

(3.46)

h  ( k / x )

Exemplo 3.3 – Condução de calor 1D – MDF – FE – Condições de contorno de 2ª e 3ª

espécies. a) Determine o campo de temperatura em regime permanente para a barra isolada lateralmente e sem termo fonte de geração de calor. Inicialmente a barra está a uma temperatura isotérmica Ti

 80 C , quando então são impostas as condições de contorno mostradas. Use o Método de

Diferenças Finitas e a formulação explícita, e 3 pontos internos, como mostrado na figura abaixo. O material da barra é o alumínio com as propriedades: condutividade térmica k

 200 W / m. C e difusividade térmica   1,00  10 4 m 2 / s . 

b) Esse problema possui solução analítica? Justifique! Se sim, mostre como obtê-la e compare com os resultados do item a. c) Explique os significados físicos quando se faz o fluxo de calor igual a zero na fronteira direita,

T  0 ? Qual seria em regime permanente o provável campo de temperatura na barra  x

q   k

com essa condição de contorno? T 1 h  200W / ( m . C ) 2 

T1,  15 C

T 2

T 3

T 5

T 4

 x



x

q  7500W / m

2

L=0,20 m

Exemplo 3.4 – Condução de calor 1D – MDF – FTI – Condições de contorno de 2ª e 3ª

espécies. Resolver o Exemplo 3.3, empregando agora a formulação totalmente implícita.

3.6 Formulações em Diferenças Finitas para malhas cartesianas não - uniformes

44



A distribuição dos pontos nodais da malha, em geral, é feita de forma irregular, de maneira a se ter mais pontos próximos às regiões de maiores gradientes. A figura 3.11 ilustra em a) um malha irregular mas ortogonal e em b) uma malha irregular e não-ortogonal. a)

b)

Fig. 3. 11 Exemplos de malhas não-uniformes: a) ortogonais e b) não-ortogonais. Para malhas não-uniformes ou curvilíneas, a discretização da equação pode ser feita após a transformação do espaço físico (x, y, z) para um espaço cartesiano computacional ( ,,). As relações entre os dois espaços são definidas pelas fórmulas de transformação de coordenadas tais como: =(x,y,z) e fórmulas similares para  e . Todas as fórmulas derivadas anteriormente podem ser aplicadas no espaço ( ,,) das equações escritas em coordenadas curvilíneas. 45



Essas equações transformadas contêm coeficientes métricos que devem ser discretizados de uma maneira consistente. Eles introduzem a influência do tamanho da malha nas fórmulas de diferenças. Dessa maneira, é possível utilizar a técnica de diferenças finitas para geometria e malhas arbitrárias. A única restrição das malhas de diferenças finitas é que todos os pontos da malha devem estar posicionados em famílias de linha que não se cruzam.

Figura 3.12 Relação entre um sistema arbitrário de coordenadas ( ,) no plano físico e o plano computacional. A malha numérica 1D e irregular, mostrada na figura 3.13, é empregada para ilustrar a discretização para uma malha cartesiana não-uniforme.

Figura 3.13 Malha unidimensional arbitrária na direção x . 3.6.1 Discretização em Diferenças Finitas para malhas cartesianas não - uniformes

Será feita a expansão em série de Taylor baseado na malha mostrada na figura 3.13. Derivada a jusante:

  u    2 u   xi 1, j 2   xi 1, j   2  ui 1, j  ui , j    2  x  x  i , j  i, j

(3.47)

 ui 1, j  ui , j    2 u   xi 1, j   u       2   x x  i 1, j  i , j   x i , j 2

(3.48)

46



Derivada a montante: ui 1, j

  u    2 u    xi , j 2  ui, j      xi , j    2   2  x  x  i , j  i , j

 ui , j  ui 1, j    2 u   xi , j   u       2  x x    i , j i , j   x i , j 2

(3.49)

(3.50)

Diferença Central.

Combinando as equações (3.48) e (3.50) para eliminar os erros de truncamento de primeira ordem, obtemos a fórmula de segunda ordem.

  u   :   x i , j

Isolando das equações (3.48) e (3.50) as expressões para 

u  u    2 u   xi1, j   u     i 1, j i , j   2   xi 1, j   x i , j   x i , j 2

  3u   xi 1, j 2   3   x 6   i , j

u  u    2u   xi , j   u     i , j i 1, j   2   xi , j   x i , j   x i , j 2

  3u   xi , j 2   3   6  x  

(3.51)

(3.52)

Somando as equações (3.52) e 3.53): 2 2   u  ui 1  ui  ui  ui 1   2u    xi   xi 1    3u    xi 1  xi         2    2       3    x x x  x  x 2 6 6      i i 1 i  i   

(3.53)

2 2   u  1  ui 1, j  ui , j  ui , j  ui 1, j    xi 1, j   xi, j   2 u   xi 1, j    xi , j       2        2 4 12  x x x    x    i, j i 1, j i , j  i  

(3.54)

Observe que o erro de trucamento para as diferenças é proporcional a dois intervalos sucessivos da malha,  x i+1 e  x i . Se a malha varia abruptamente, então a fórmula terá acurácia de primeira ordem. Esta é uma propriedade geral da aproximação por diferenças finitas em malhas não-uniformes. Se a m a l h a n ão v a r i a su a v e m e n t e a p e r d a d e a c u r ác i a é i n e v i t áv e l .

A fórmula para a derivada central de segunda ordem é:

47


Métodos Numéricos Aplicados à Engenharia ui 1, j  ui , j ui , j  ui 1, j    2 u  1    3u  2 1    2       xi 1, j   xi , j  3    2 3  x x x x x       i i 1, j i , j i , j   x i   i 1, j

 x     x     u   12 x    x    x  3



3

i 1, j

(3.55)

4

i , j

4

i 1, j

i , j

3.6.2 Transformação geral das equações

Vamos considerar um escoamento bidimensional não-permanente, com as variáveis independentes x, y e t. Vamos transformar as variáveis no espaço físico ( x, y, t ) para o espaço transformado ( ,

, ) onde: = ( x,y,t ) = ( x,y,t )

(3.56)

= ( x,y,t ) Na transformação acima,  é considerada função de t somente, e é freqüentemente dada por =t Da regra da cadeia do cálculo diferencial:

=0

                                           x  y ,t     ,t   x  y ,t     ,t   x  y , t     ,   x  y ,t

(3.57)

Os subíndices indicam as variáveis que são mantidas constantes na diferenciação parcial.

                    x      x       x 

(3.58)

                     y      y       y 

(3.59)



Analogamente: 

Também temos:

       t  x, y

                                        ,t   t  x, y     ,   t  x, y     ,   t  x, y

(3.60)

ou

                  d                 t       t       t      d t

(3.61)

48



Os coeficientes:

       , , e são chamados fatores métricos.  x  y  x  y

A derivada segunda nas direções x e y são:  2  x 2

2

    2      2    2        2     2    2         x       x       x  2

(3.62)

                    2    2    x           x   x  2

 2  y 2

2

2

    2      2    2        2     2    2         y       y       y  (3.63) 2

  2       2             2    2    y           y   y  A aproximação da equação de Laplace para um sistema de coordenadas curvilíneas é mostrada abaixo.  2  x 2



 2  y 2

0

(3.64)

Substituindo (3.63) e (3.64) em (3.65), obtemos:  2    

2

            2   x    y  

2

   2     2     2     2               x    y  

2

  2       2                             2    2      x          x   x    y   y        

2

             x    y  

2

(3.65)

       2  2     2   2  2   0       x  y 

3.7 Aproximação da Equação de Condução de Calor Bidimensional. 49



Considere o problema bidimensional de condução de calor, com propriedades constantes, regido pela equação (3.67). 1T   t





2

T

 x

2



 2 T 2

 y



q k

(3.66)

A equação discretizada pelo MDF possui a seguinte forma: n 1

TP

 T E  2T P T W T N  2T P  T S  q    TPn       O[(t,( x) 2 ,( y ) 2 ]  2 2 ( x) (  y) t   k

(3.67)

A escolha do valor de n    n  1 determina o tipo de formulação. As tempreaturas do fluxo de calor são obtidas da Eq. (3.40). Na Fig. 3.14, as dimensões L e W correspondem ao comprimento e altura da placa bidimensional, respectivamente. h, T  ; q”; T

h, T 

h, T 

q”

q”

T

T

W

h, T  ; q”; T

L

Figura 3.14 Malha 2D para discretização da Eq. (3.67). Nas próximas seções serão apresentadas as formulações explícita e totalmente implícita para os domínio 2D da Fig. 3.14. 3.7.1 Formulação Explícita 2D (  n) .

Fazendo   n na equação (3.68), obtemos: TP 1  TP n

n

t TPn

1

 T n  2T nP T Wn T Nn  2T nP T nS  q     E    k  2 2 ( ) ( ) x y    

 1  2(rx  ry )  TPn  rx TEn  TWn   ry T Nn  TSn   Q P

(3.68)

(3.69) 50



Sendo: r x 

  t

  x 

2

e r y 

  t

  y 

2

. A definição de QP é dada na Eq. (3.41).

O critério de estabilidade para a Eq. (3.70) é obtido fazendo o coeficiente de T Pn sempre maior ou igual a zero. (1  2r x  2 r y )  0 

t 

 t

  t

  x

 y 

 2

2



1 2

1

 1 1   2 2    x   y  

2 

(3.70)

A Eq. (3.71) determina o valor máximo do passo de tempo admissível para a formulação explícita. Qualquer valor menor de t pode ser utilizado. A Eq. (3.70) em geral é escrita da forma: n 1

TP

 aW TWn  a ETEn  a N TNn  a STSn  b

(3.71a)

Onde:

aPn  (1  2( rx  r y ))  a E  aW  r x a  a  r  N S y b  a n T n  Q P P P 

(3.71b)

Se  x  y  r x  ry  Fo , então a Eq. (3.70) pode ser reescrita como: n 1

TP

 1  4 Fo  TPn  Fo TEn  TWn   Fo TNn  TSn   QP

(3.71)

O conjunto de equações algébricas resultantes da aplicação da Eq. (3.70) a um domínio 2D, pode ser resolvido, por exemplo, pelo método de Jacobi quando desejamos conhecer o processo transitório; ou pelo método de Gauss-Seidel desejando apenas o regime permanente. Exemplo 3.5 – Condução de calor 2D – MDF – FE – Condições de contorno de 2ª e 3ª

espécies. a) Determine o campo de temperatura em regime permanente, para a placa retangular mostrada na figura abaixo. Use o método de Diferenças Finitas com 3 pontos internos na direção x e 2 pontos internos na direção

y ,

como mostrado na figura. Use aproximação polinomial para

tratamento das condições de contorno de 2ª espécie (fluxo de calor especificado). Inicialmente a 51



placa está a uma temperatura

T ij = 20oC .

O material da placa é o aço AISI 304 com

condutividade térmica k = 14, 2W / (m K ) e difusividade térmica

a = 3, 047 ´ 10- 6 m 2 / s .

b) Trace o gráfico de uma das temperaturas internas versus o tempo ( T ij ´ tempo ). Comente o comportamento da curva.

qt ¢¢ = 1200W / m 2

y

H = 0,27 m

T r ,¥ = 5oC

T l , ¥ = 5oC h = 400W / (m 2K )

h = 400W / (m 2K )

x Tb = 0oC L = 0, 60 m

3.7.2 Formulação Totalmente Implícita 2D (  n  1) .

Fazendo   n  1 na equação (3.68), obtemos: n 1

TP

 T nE1 2T nP1 T Wn1 T nN1 2T Pn1 T nS1  q   TPn      k  (  x) 2 ( y) 2 t  

r x TEn1  TWn1   1  2(rx  ry )  TPn 1  ry TNn 1  TSn 1   TPn Q P Sendo: r x 

  t

  x 

2

e r y 

  t

  y 

2

(3.72)

(3.73)

. A definição de QP é dada na Eq. (3.41).

A Eq. (3.74) não necessita de um critério de estabilidade, sendo o valor do  t normalmente restringido por questões de precisão. Se  x  y  r x  ry  Fo , então a Eq. (3.74) pode ser reescrita como:

 Fo T En 1  TWn 1   1  4 Fo  TPn 1  Fo TNn 1  TSn 1   TPn  QP

(3.74)

52



A Eq. (3.74) pode ser escrita na forma: n 1

a PTP

 aW TWn1  aETEn 1  a STSn 1  a N TNn 1  b

(3.75)

onde:

a P  (1  2( rx  r y ))  aW  aE  r x  a S  aN  r y b  T n  Q P P 

(3.76)

O sistema resultante da aplicação da Eq. (3.74) a um domínio 2D, tem a forma: AT

B

(3.77)

Devido ao acoplamento das equações, o sistema acima pode ser resolvido pelos métodos de inversão de matrizes, eliminação de Gauss, Cholesky, Gradiente Conjugado, discutidos nas próximas seções. Exemplo 3.6 – Condução de calor 2D – MDF – FTI – Condições de contorno de 2ª e 3ª

espécies. Resolver o Exemplo 3.5 empregando agora a formulação totalmente implícita. 3.8 Estrutura da Matriz de Coeficientes.

A estrutura da matriz de coeficientes obtida na aproximação numérica é de fundamental importância na escolha do método de solução do sistema linear. A matriz de coeficientes do problema para os casos uni, bi, ou tridimensional é mostrada abaixo:  X X  X X   X       0    

0

X X X X X X X X X

X X X X X X X X X

(a)

           X  X 

 X X  X X   X  X 0   X     0   

X X

0 X

X X

0 X

X X X

X

X X X X

0 X

X X X X

0 X X X X

X X X

X

        X    X  X 

X  X X  X X X 0   X X X   X 0 X X  X X  X 0   X 0 X   X 0 X 

(b)

X

0

  X     X  0 X   X  X X   X X 

0 X

X X

X 0 X X X X X X X X X 0 X X

(c)

Figura 3.15 Estrutura da matriz de coeficientes

A

.

53



a) problemas 1 - D, matriz tridiagonal b) problemas 2 - D, matriz pentadiagonal c) problemas 3 - D, matriz heptadiagonal. A estrutura da matriz é função de no cálculo do valor de um nó estarem envolvidos os valores dos nós vizinhos. Também de termos utilizado na discretização da equação derivadas centrais. É possível envolver mais pontos do que somente os vizinhos adjacentes para construir a matriz do problema. Se isto for feito, o número de zeros naquela linha se altera. Pode-se envolver todos os pontos do domínio para a discretização de cada equação, acarretando uma linha cheia (sem zeros), no sistema de equações. A maioria das metodologias emprega o esquema de 5 pontos (célula de 5 pontos) já visto. Assim, a matriz de rigidez do problema,

A

, é a mais simples possível.

Uma característica importante das matrizes obtidas das aproximações numéricas é o alto índice de esparsidade, como mostrado na Tab. 3.1. O índice de esparsidade é definido pela porcentagem de nulos comparados com os não-nulos da matriz. As matrizes até agora estudadas são originárias de aproximações numéricas estruturadas, cujos volumes possuem sempre o mesmo número de vizinhos. Em discretizações não-estruturadas, podemos ter diferentes números de vizinhos para cada volume, originando matrizes com banda variável. Tabela 3.1 Indice de esparsidade em problemas bidimensionais. Malha Número de Elementos da Matriz Não-zeros % de não-zeros Índice de Esparsidade Memória requerida* [Kbytes]

10x10

20x20

40x40

104

16x104

256x104

500 5% 0,95

2000 1,25% 0,9875

8000 0,31% 0,9968

40/2

640/8

10240/32

* O número superior indica a quantidade de Kbytes necessária para armazenar em simples precisão a matriz completa. O número inferior indica quantos Kbytes serão necessários para armazenar somente os elementos não-nulos desta matriz.

3.8 Solução do Sistema Linear de Equações.

Nesta seção serão apresentados alguns métodos de solução de sistemas lineares comumente utilizados para resolver o sistema de equações gerado pelas discretizações. O sistema de equações a ser resolvido tem a seguinte forma:

54



 E1 :  E :  2    En :

a11x1 a21x1

 a12 x2   a1n xn  b1  a22 x2  a2 n xn  b2





a n1x1



(3.78)



an 2 x2  ann xn  bn

O qual pode ser escrito na forma matricial como: AX

B

(3.79)

Em geral há dois métodos de solução do sistema linear de equações algébricas,

X,

a

saber, os métodos Direto e Iterativo. Métodos Diretos: são aqueles que necessitam da inversão da matriz completa, incluindo os não - zeros. Há alguns métodos diretos, como o método de Cholesky, que necessitam apenas do armazenamento da banda da matriz, mas exige que a matriz seja de banda, simétrica e positiva definida. Os métodos direitos, além da enorme memória requerida, exigem o manuseio dos nãonulos, influindo grandemente na taxa de convergência do método. Métodos Iterativos: os métodos iterativos podem ser classificados em iterativos ponto a ponto, linha a linha ou plano a plano.

3.8.1 Métodos Diretos.

a) Inversão de Matriz O método de inversão de matrizes consiste em encontrar o valor de pelo vetor do lado direito,

B ,

A

1

que multiplicado

do sistema de equações, conduza diretamente a solução.

Normalmente, o método de inversão de matrizes aplica-se quando os sistemas lineares são pequenos, por limitações de memória e necessidade de envolver todos os elementos da matriz no processo de solução. Após a inversão da matriz, a solução pode ser obtida como: X

 A 1 B

(3.80)

b) Método de eliminação de Gauss

55



O método de eliminação de Gauss consiste em triangularizar a matriz, de modo a obterse o valor de xn e executar a substituição retrógrada. A matriz de coeficientes

A e

a matriz lado direito

B do

sistema de Eqs. (3.79), após a

triangularização pode ser reescrita como:

a11 x1  0   0   0

a12 x2



a1n xn

a22 x2



a2 n xn

 b1  b2

0







0

0

ann xn

 bn

(3.81)

Os valores de aij e b j da matriz não são necessariamente os mesmos valores iniciais. O valor de xn pode ser determinado como: xn



bn ann

O valor de xn1 pode ser determinado em função de xn : xn 1 

bn 1  an 1, n xn an 1, n 1

Continuando o processo para as outras equações, chegamos ao seguinte algoritmo: xi



bi

  j 1 ai , j x j n

(3.82)

ai ,i

para i = n-1, n-2, ..., 2, 1. Exemplo 3.7 – Resolver o sistema abaixo utilizando o método de eliminação de Gauss.

0 0  2 x1  x2  x  2 x  x 0  1 2 3  0  x2  2 x3  x4   0  0  x3  2 x4

1 0 0 1

(3.83)

O sistema acima pode ser reescrito na forma:

1  E 2 :  1 2  E 3 :  0  1  0 E 4 :  0 E 1 :  2

0

0

1 

 1 0 0 2  1 0   1 2 1

Fazendo (E 1 + 2E 2 ) E 2: 3x 2 - 2x 3 = 1 Fazendo (E 2 + 3E 3 ) E 3: 4x 3 - 3x 4 = 1 56



Fazendo (E 3 + 4E 4 ) E 4:

5 x 4 = 5

Reescrevendo o sistema: E 1 : 2

1

0

E 2 : 0

3

2

0

4

0

0

  E 3 : 0  E 4 : 0

0

1

0

1

   3 1  5 5

Após a triangularização do sistema, a aplicação da fórmula de recorrência (3.83) conduz a solução:

1 1 1 1

 x1  x  2  x  3  x4

3.8.2 Métodos Iterativos.

c) Método de Jacobi Neste método, a variável dependente de cada ponto da grade é resolvida, usando-se valores de um “chute” inicial ou valores previamente calculados. O novo valor de uij na nova iteração n  1 é: 1

u in, j



1 4

u 

n i 1, j

 u in1, j  u in, j 1 u in, j 1 

(3.84)

onde n corresponde aos valores previamente estimados ou calculados. O processo iterativo é repetido até que o critério de convergência especificado seja atendido. Esse método é raramente usado (ou nunca) para a solução de equações elípticas. De forma generalizada: AP TP

 AE TE  AW TW  AN TN  AS TS  B

(3.85)

A equação (3.61) pode ser reescrita como: AP1 TP 1 n

n

  Anbn Tnbn  B

(3.86)

O ciclo iterativo é: 57



 Estimar: campo inicial da variável;  Iterar em n;  Calcular T p usando a Eq. (3.86) para todos os pontos do domínio;  Checar a convergência;  Retornar se o critério não foi satisfeito. O método de Jacobi é de lenta convergência e requer uma matriz com dominância diagonal. d) Método de Gauss-Seidel Neste método, o valor da variável dependente é utilizado para calcular os pontos vizinhos tão logo ele seja disponível. Isto certamente aumenta a taxa de convergência dramaticamente sobre o método de Jacobi (~ 100%). O método é convergente se os maiores elementos são localizados na diagonal principal da matriz coeficiente (dominância diagonal). O processo de iterativo pode ser resumido como:

 Estimar campo inicial da variável  Iterar em n 

Calcular T P pela Eq. (3.86) onde foi considerada uma varredura de sul para norte e de oeste para leste, podendo-se considerar conhecidas, na mesma varredura, T W e T S .



Checar convergência

 Retornar se o critério não foi atendido n 1

AP TP

  Anb Tnbn  AEn TEn  ANn TNn  B

(3.87)

Para um algoritmo vetorizado, o método de Jacobi pode ser mais eficiente que o de Gauss-Seidel. e) Método das Sobre - Relaxações Sucessivas - S.O.R O Método de iteração Gauss-Seidel pode ser acelerado introduzindo-se um parâmetro de relaxação. A estrutura do processo iterativo é:

 Estimar campo inicial da variável;  Iterar em n; 58





Calcular os valores de T usando a Eq. (3.88);

 Sobre - relaxar usando a Eq. (3.88) ; 

Checar convergência;

 Retornar se o critério não foi satisfeito. T Pn

1

 w T Pn 1 GS 1  w T Pn

(3.88)

onde T Pn 1 GS representa o valor calculado com o método de Gauss-Seidel e

w o

coeficiente de

relaxação. O coeficiente de relaxação serve para avançar mais rapidamente a solução, quando o processo esta lento, ou “segurar” a variável quando a mesma avança em demasia e pode causar divergência. Os valores recomendados de

w para

avançar mais rapidamente a solução variam

de 1.5 a 1.7, apesar deste valor depender do tamanho da malha. Na prática, o método de tentativa e erro é usado para determinar

w ótimo para um problema em particular.

f) Método Linha a Linha Os métodos linha a linha resolvem diretamente uma linha. O método mais conhecido é o algoritmo de Thomas ou (TDMA). Para problemas 2-D e 3-D são iterativos, com a varredura se processando linha a linha e coluna por coluna. Considere a malha unidimensional abaixo:

1

2

3

...

...

...

N-1

N

Figura 3.16 Malha unidimensional, sendo 1 e N são pontos sobre a fronteira. O valor da temperatura no ponto i pode ser determinado pela seguinte equação algébrica: aiT i  biT i 1  ciT i 1  d i

(3.89)

Quando as condições de contorno são dadas, as equações para o ponto próximo a fronteira tomam a forma trivial. Por exemplo, se T 1 é dado, nós temos a1=1, b1=0, c1=0 e d 1= T 1. Essas condições implicam que T 1 é conhecido em termos de T 2. A equação para i =2 é uma relação entre T 1, T 2 e T 3. Como T 1 pode ser expresso em termos de T 2, esta relação reduz-se a uma relação entre T 2 e T 3. Em outras palavras, T 2 pode ser expresso em termos de T 3. Este processo de substituição pode ser continuado até que T N seja formalmente expresso em termos de T N+1. Como T N+1 não tem significado, nós obtemos o valor de T N neste estágio. Isto permite 59



que façamos agora a substituição retrógrada na qual T N-1 é obtida de T N , T N-2 de T N-1, T 2 de T 3 e T 1 de T 2.

Suponha que na substituição avançada, nós tenhamos uma relação: T i  PiT i 1  Qi

(3.90)

Antes obtivemos a seguinte relação: T i 1  Pi 1T i  Qi 1

(3.91)

Substituindo a Eq. (3.92) em (3.91): aiT i  biT i 1  ci ( Pi 1T i  Qi 1 )  d i

(3.92)

Rearranjando a equação acima para parecer-se com a Eq. (3.91)



  c Q d  T i 1   i i 1  i   ai  ci Pi 1   ai  ci Pi 1 

T i  

bi

(3.93)

Comparando a Eq. (3.94) com a Eq. (3.91): Pi

 bi     a c P   i i i 1 

(3.94)

Qi

 c Q  d    i i 1 i   ai  ci Pi 1 

(3.95)

Estas são relações de recorrência, desde que elas dão o valor de P i e Qi em termos de P i-1 e Qi-1. Para começar o processo de avanço, escrevemos a Eq. (3.91) para o ponto i =1. a1T 1  b1T 2  d 1

(3.96)

Dividindo por a1: T 1 

b1 a1

T 2 

d 1

(3.97)

a1

Onde os valores de P 1 e Q1 são: P1 

b1 a1

,

Q1



d 1 a1

(3.98)

Para o ponto N, temos b N =0  P N =0. 60



Portanto: T N  Q N

(3.99)

Podemos agora executar a substituição retrógrada via a Eq . (3.91). O algoritmo para o problema pode ser resumido da seguinte forma: Algoritmo TDMA Calcule P 1 e Q1 da Eq. (3.96) Use as relações de recorrência (3.95) e (3.96) para obter P i e Qi para i =2, 3,...,N. Faça T N =Q N . Use a Eq. (3.91) para i =N-1, N-2,..., 3, 2, 1 para obter T N-1, T N-2, ...,T 3, T 2, T 1. Exemplo 3.8 – Resolver o sistema abaixo utilizando o método TDMA.

0 0  2 x1  x2  x  2 x  x 0  1 2 3  0  x  2 x  x 2 3 4   0  0  x3  2 x4

1 0 0 1

(3.100)

A nossa equação deve estar na forma: aiT i  biT i 1  ciT i 1  d i

que foi utilizada para deduzir as relacões de recorrência (3.95) e (3.96). As equações no sistema (3.96) estão escritas na forma: ciT i 1  aiT i  biT i 1  d i

Portanto, as Eqs. (3.95) e (3.96) são reescritas como: Pi

  bi      a c P  i i i 1 

(3.95)

Qi

 d  c Q    i i i 1   ai  ci Pi 1 

(3.96)

Os valores de P 1 e Q1 são reescritos como: P1 

 b1 a1

, Q1 

d 1 a1

61



Para a primeira equação do sistema (3.96), utilizando a equação acima obtemos:



d 1  b1 1  , Q1  1 

n=1

P1

n=2

1 3 (a2  c2 P1 )  (2  1. )  2 2

a1

2

2

a1

  b2  1 2    ; P2    a c P 2 2 1   3 3

Q2

 d  c Q    2 2 1    a2  c2 P1 

2

n=3

3

1 2

2 6

 

1 3

2

2 4 (a3  c3 P2 )  ( 2  1. )  3 3

  b3  1 3    ; P3    a3  c3 P2  4 4 3

n=4

0  1.

( a4

3 4

 c4 P3 )  (2  1. ) 

1

 d 3  c3Q2  0  1. 3 1   Q3    4 4 a c P  3 3 2   3

5 4

P4  0

1

5

 d 4  c4Q3  1  1. 4 4   Q4    1 5 5 a c P   4 4 3  4

4

Da Eq. (3.96): x4

 Q4  1

Utilizando a Eq. (3.91) para executar a substituição retrógrada: x N  P N x N 1  Q N

3

1

4

4

2

1

3

3

x3

 P3 x4  Q3  .1   1

x2

 P2 x3  Q2  .1   1

x1  P1 x2  Q1 

1 2

.1 

1 2

1

Finalmente, temos a solução do sistema linear:

 x1  x  2   x3  x4

1 1 1 1 62



No método TDMA é importante observar as condições de contorno dominantes para realizar o processo de varredura, principalmente nesta direção. e) Outros métodos utilizados na solução de sistemas lineares: MSI - Modified Strongly Implicit; ADI - Alternating Direction Implicit; CGM - Conjugate Gradient Method; ICCG - Incomplete Cholesky Conjugated Gradient. 3.9 Consistência, Estabilidade e Convergência.

Os problemas da engenharia e da física dão origem a sistemas de equações complexas sobre cujos comportamentos matemáticos pouco se conhece. A dificuldade está em se obter as condições (tamanho de malha, tamanho de intervalo de tempo, coeficientes de relaxação, etc..) para que as aproximações numéricas sejam estáveis e convergentes. A terminologia para esses conceitos importantes em simulação numérica são dados abaixo. Consistência: a aproximação numérica deve reproduzir a equação diferencial quando os tamanhos da malha espacial e temporal tenderem a zero. Os erros de trucamento na série de Taylor tendem a zero quando ( x , t )0. Em resumo, as equações discretizadas devem tender às equações diferenciais, quando a malha tender a zero. Estabilidade: um esquema numérico é estável se a solução numérica obtida é a solução exata das equações discretizadas. Erros de arredondamento da máquina e dificuldades no tratamento de acoplamento entre variáveis pode causar instabilidade. Convergência: consistência e estabilidade são condições necessárias e suficientes para a convergência. Portanto, a solução numérica é convergente quando é estável e tende para a solução das equações diferenciais quando a malha é refinada.

63




1) a) Determine o campo de temperatura em regime permanente, para a placa retangular mostrada na figura abaixo, transcorridos 5000 s do instante inicial. Use o método de Diferenças Finitas e cinco pontos internos em cada direção. Inicialmente a placa está a uma temperatura Tij

 20 C . O material da placa é o cobre com condutividade térmica

k

 250 W / m. C e 

difusividade térmica   1,0 10 5 m 2 / s . b) Trace o gráfico de uma das temperaturas internas versus o tempo ( T ij  tempo ). Através desse gráfico é possível dizer se estamos ou não próximos do regime permanente? Justifique sua resposta. c) Estime a temperatura do centro da placa em regime permanente.

y

H=0,9 m

o

T t = 200 C

o

T e = 100 C

o

T d = 100 C

x o

T b = 200 C

L=0,6 m

2) a) Determine o campo de temperatura em regime permanente para a barra com termo fonte uniforme de geração de calor. Inicialmente a barra está a uma temperatura isotérmica Ti

 50 C , quando então são impostas as condições de contorno mostradas. Use o Método de

Diferenças Finitas e 3 pontos internos para obter o campo transitório de temperatura. O material da barra é o alumínio com propriedades: condutividade térmica k  200 W / m. C e

64



difusividade térmica   1,00  10 4 m2 / s . O termo fonte é: q  1, 0 105 W / m3 . T1  30 C , r1

 20 mm , T5  80 C e 

r2

 40 mm .

b) Compare os resultados com a solução analítica em regime permanente.

T5 r 2

r 1

1 ¶ T 1 ¶ = a ¶ t r ¶r

T1

r

θ

æ ¶T ççr çç ¶ r è

ö q ¢¢¢ ÷ ÷ + ÷ ÷ k ø

3) Um reservatório tem forma de um cilindro de raio R t = 0, 5 m e contem água a uma altura ho = 2, 0 m . Abre-se então o orifício de raio R 2 = 0, 02 m na base do reservatório, drenando a

água com velocidade V 2 =

2gh . A EDO que rege a altura instantânea de água h = h (t ) no

reservatório é dada por: dh = dt

cuja solução analítica é dada por: t =

2ho g

2

æR ö ÷ 2gh çç 2 ÷ è R t ø

2 æR t ö æ ÷ ç ÷ çç1 ÷ çè çè R 2 ø

1/ 2

æh ö ç ÷ ÷ ÷ èç ho ø

ö ÷. ÷ ÷ ÷ ø

a) Determine o tempo necessário para esvaziamento total do reservatório, t esv , resolvendo o problema numericamente com o método de Diferenças Finitas. Pretende-se conhecer o processo transitório com boa precisão. b) Trace o gráfico ( h ´ t ) e discuta-o. Compare com os resultados da solução analítica.

65



Rt

ho h(t)

R2 V 2

4) Resolver o sistema linear abaixo usando os seguintes métodos: a) Eliminação de Gauss b) Jacobi c) Gauss-Seidel d) Sobre-Relaxações Sucessivas ( w = 1,6) e) TDMA

 3 x1  x2  0  x4  E 2 :   x1  7 x2  2 x3  2 x4  E 3 :   3 x1  x2  6 x3  x4 E 4 :   x2  x3  3 x4  x1 E 1 :

3 2 1 6

66

Apostila Métodos Numéricos - UTFPR

Recommend Documents