Universidade Federal de Santa Catarina –UFSC
Departamento de Engenharia Química e Alimentos
Volumes Finitos
Por
Ricardo Vicente de Paula Rezende
Elaine Vosniak Takeshita
Gabriela Iris da Silva
Florianópolis, Janeiro de 2003
Prepare-se para as provas
Estude fácil! Tem muito documento disponível na Docsity
Ganhe pontos para baixar
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Guias e Dicas
Prepare-se para as provas
Estude fácil! Tem muito documento disponível na Docsity
Encontrar documentos
Prepare-se para as provas com trabalhos de outros alunos como você, aqui na Docsity
Encontra documentos específicos para os exames da tua universidade
Encontra documentos específicos para os exames da tua universidade
Videoaulas
Prepare-se com as videoaulas e exercícios resolvidos criados a partir da grade da sua Universidade
Quiz
Responda perguntas de provas passadas e avalie sua preparação.
Ganhe pontos para baixar
Ganhe pontos ajudando outros esrudantes ou compre um plano Premium
Volumes Finitos2, Notas de estudo de Engenharia Mecânica
Apostilha sobre o metodo dos volumes finitos
Tipologia: Notas de estudo
Antes de 2010
1 / 32
Esta página não é visível na pré-visualização
Não perca as partes importantes!
Documentos relacionados
Pré-visualização parcial do texto
Baixe Volumes Finitos2 e outras Notas de estudo em PDF para Engenharia Mecânica, somente na Docsity!
Universidade Federal de Santa Catarina –UFSC
Departamento de Engenharia Química e Alimentos
Volumes Finitos
Por
Ricardo Vicente de Paula Rezende
Elaine Vosniak Takeshita
Gabriela Iris da Silva
Florianópolis, Janeiro de 2003
Índice
- I NTRODUÇÃO
- M ODELOS M ATEMÁTICOS
- O M ÉTODO
- Condições de Contorno
- Métodos de solução de um sistema linear
- Método Gauss-Seidel
- Método TDMA (TriDiagonal Matrix Algorithm)
- Funções de Interpolação
- Esquema Power-Law
- Acoplamento Pressão-Velocidade
- Método SIMPLE
- CONSIDERAÇÕES FINAIS
- NOMENCLATURA
- BIBLIOGRAFIA
Apesar de tudo os métodos numéricos apresentam pontos delicados em sua utilização, as principais são:
¾ Formulação do Modelo Matemático , o qual seve ser condizente com o fenômeno em estudo; ¾ A escolha do Método Numérico adequado, pois constitui uma solução dos modelos matemáticos, e sua escolha deve atender à complexidade do fenômeno; ¾ Implementação Computacional, que deve minimizar os erros no processo de cálculo e permitir o menor esforço computacional e de programação.
O ideal é aliar as vantagens dos métodos, utilizando o método analítico no auxílio dos modelos numéricos e o experimental em sua validação. Quando se opta por um Método Numérico tem-se em mente, segundo Maliska (1995), o seguinte objetivo:
“Resolver uma ou mais equações diferenciais, substituindo as derivadas existentes na equação por expressões algébricas que envolvem a função incógnita.”
Entre os diversos métodos disponíveis, podemos citar os mais utilizados em problemas de engenharia: ( A ordem não implica em uma hierarquia).
¾ Elementos Finitos; ¾ Elementos de Contorno; ¾ Diferenças Finitas; ¾ Volumes Finitos.
Vamos focar a nossa atenção no Método de Volumes Finitos ( MVF ).
Copyright 2003. All rights reserved.
Modelos Matemáticos
Um grande número de problemas de Engenharia se concentram em fenômenos de transferência de massa e energia, sendo que este último, em geral, refere-se ao calor e à quantidade de movimento. Outros fenômenos podem estar envolvidos mas, quase sempre, estes são os de maior interesse para a Engenharia. Genericamente, a formulação do modelo deve atender aos seguintes critérios (Smye e Clayton 2002):
¾ Acuracidade : O modelo deve confirmar, com o mínimo de erros, os dados experimentais. O que requer uma compreensão profunda do fenômeno; ¾ Predição : O modelo deve ser preditivo, deve poder avaliar se um evento é ou não possível, viável ou não, e em que condições ele pode ocorrer.; ¾ Economia : Deve saber avaliar quais os parâmetros essenciais do fenômeno. Isto requer habilidade por parte do proponente do modelo; ¾ Estabilidade : não ser sensível a “ruídos numéricos” e/ou experimentais; ¾ Utilidade : deve ser útil à prática acima de tudo.
As formulações matemáticas – em coordenadas cartesianas - em geral tomam a forma da Eq.(1):
( ) ( ) ( ) ( ) φ^ φ φ + φ
∂ S
t x u y v z w x x y y z z (1)
onde ρ : densidade; u,v,w: componentes da velocidade nas direções x,y e z respectivamente; φ : propriedade em estudo; Γ^ φ : coeficiente que depende da propriedade φ. Ver Tabela 1.
No lado esquerdo da igualdade, temos o termo transiente e os termos convectivos nas três direções. No lado direito, os termos referentes aos processos de difusão nas três direções e o termo fonte da
Copyright 2003. All rights reserved.
Existem diversas maneiras de discretizar o domínio. A discretização gera Malhas que podem ou não ser uniformes e estruturadas. A Figura 2 ilustra o conceito. O número de pontos e o tipo de malha adequados dependerão tanto do tipo de fenômeno quanto do nível de precisão requerida.
Figura 1- Domínio contínuo Figura 2 – Domínio discretizado ( Malha)
A unidade fundamental do domínio é o Volume de Controle (VC), mostrado na Figura 3, e a sua geometria é dependente do tipo de malha utilizada. Como o domínio possui um número finito de regiões, o número de VC gerados, logicamente, deverá também ser finito.
Figura 3 – Volume de Controle
No MVF existem duas formas das equações aproximadas serem obtidas: a) realizando balanços da propriedade no VC; ou b) integrando a equação diferencial em sua forma conservativa no tempo e no espaço. Ambas representam balanços de tudo que atravessa as fronteiras e é gerado ou consumido no VC,
Copyright 2003. All rights reserved.
e, fundamentalmente, são a mes ma coisa. O processo de obtenção das equações aproximadas é o processo de discretização do modelo. Para uma melhor compreensão, suponha um balanço de massa em um volume elementar (VC) bidimensional em estado estacionário representado na Figura 4. A coordenada em z foi igualada a unidade, assim:
Figura 4- Balanço diferencial em um volume de controle bidimensional
ρ u ∆ ye −ρ u ∆ yw +ρ v ∆ xn −ρ v ∆ xs = 0 (2)
Dividindo-se a equação por ∆ y ∆ x e aplicando o limite chegamos a equação diferencial escrita na
forma conservativa :
∂∂ x^^ (ρ^ u )^ +∂∂ y (^ ρ v )^ =^0 (3)
Integrando-se a Eq. (3) somente no espaço ( o regime é estacionário, independente do tempo portanto) tem-se
(^) ρ ∂ ρ + ∂ ∂
∫ ∫ (^) x u y v dydx
e w
n s
Copyright 2003. All rights reserved.
Figura 5 – Discretização unidimensional do domínio de cálculo
Deve-se agora integrar a Eq.(7) no tempo e no espaço:
( ) dxdt
x
T
C
k t dxdt x
T t t t
e w p
t t t
e ∫ ∫ w ∫ ∫
+∆ +∆
∂ ρ (8)
( T T ) dx Ck Tx Ck Tx dt
t t t p e p w
e w ∫^ ∫
+∆
ρ −ρ^0 0 = ∂ (9)
Precisamos escolher uma função de interpolação espacial para temperatura, pois devemos avaliar a sua derivada nas faces dos VC’s. Este é um ponto muito delicado do método como será visto mais adiante. Podemos para isso usar Diferenças Centrais, que nos apresenta uma aproximação linear para as derivadas, como mostrado nas Equações 10 e 11:
e
E P e x
T T
x
T
∂ θ (10)
w
P W w x
T T
x
T
∂ θ (11)
Copyright 2003. All rights reserved.
De posse destas equações e integrando a Eq. 9 no tempo, temos:
w
P W e p
E P p
P P P P x
T T
C
k x
T T
C
k t
M T M T
− 0 0 θ θ θ θ (12)
onde M (^) P^0^ = ρ^0 P ∆ x e M (^) p = ρ P ∆ x.
Agrupando os coeficientes e resolvendo para Tp temos:
AP TP = AeTE θ^ + AwTW θ+ AP^0 T P^0 (13)
Onde os termos Ai representam os coeficientes agrupados de cada termo da temperatura em sua respectiva posição no domínio de cálculo. O termo T (^) P é o que queremos avaliar no instante seguinte t+ ∆ t. O termo TP^0 , representa o valor da propriedade avaliado em relação ao próprio ponto P no instante anterior [ t ] ao que se quer calcular [ t+ ∆ t ], por isso o índice 0. Os pontos TE^ θ^ e TW^ θ^ , os pontos a leste [E] e a oeste [W] do ponto P respectivamente. O índice θ refere-se a que posição na dimensão temporal a temperatura em E e W está sendo avaliada. A temperatura pode ser avaliada em t , em t+ ∆ t ou em uma posição intermediária. Como foi feito para o espaço, precisamos de uma função de interpolação da temperatura no tempo. Necessitamos de uma função que descreva seu comportamento no intervalo. Uma escolha aceitável é uma função linear que pode ser dada como a Eq. (14).:
T θ^ =θ T + ( 1 −θ) T 0 (14)
De acordo com a escolha para θ temos três formulações possíveis:
¾ Explícita , onde θ =0:
Na formulação explícita o ponto P em t+ ∆ t é avaliado em relação a todos os pontos vizinhos a P em t , e estes valores são conhecidos gerando um conjunto de equações não acopladas; o processo anda em “marcha” no tempo. Esta formulação possui limitação quanto ao uso de ∆ t pois valores muito pequenos
Copyright 2003. All rights reserved.
Escolhido o tipo de formulação, substitui-se a função de interpolação nos pontos E e W na Eq. (13), e resolve-se a equação resultante para Tp reagrupando os coeficientes. Apesar disto, a estrutura da Eq.(13) não será modificada. A diferença reside nos termos que compõem os coeficientes e além disto, ao termo AP^0 Tp^0 poderão somar-se, dependendo da formulação, os termos Ae^0 TE^0 e Aw^0 TW^0 que são avaliados no instante t , e são, portanto, valores conhecidos. Poderão também aparecer alguns termos referentes ao termo fonte ( caso este seja considerado no modelo ). Todos estes termos podem ser agrupados em um novo coeficiente que chamaremos de B. Assim, para este nosso exemplo a equação, independente da formulação usada, ficaria na forma:
A (^) P TP = AeTE θ^ + AwTW θ+ B (15)
Este foi um exemplo unidimensional onde apenas um tipo de fenômeno foi considerado. No entanto, o mesmo raciocínio pode ser empregado na formulação de um fenômeno tridimensional e composto por todos os termos da Eq.(1), sem acarretar grandes alterações na Eq.(15). Generalizando a Eq.(15) e a Eq.(14) em coordenadas cartesianas e com base na Eq.(1), temos:
AP φ (^) P = Ae φθ E^ + Aw φ W θ + An φ N θ + As φ S θ+ At φ T θ+ Ab φθ B + B (16)
φ θ^ =θφ+ ( 1 −θ)φ 0 (17)
onde os índices N e S representam os pontos ao Norte e ao Sul de P, e os índices T e B os pontos Acima (Top) e Abaixo ( Botton) de P, respectivamente.
A Eq(16) pode ainda ser escrita de uma forma mais sintética:
AP φ (^) P = ∑ AnbTNB θ^ + B (18)
onde o índice NB refere-se aos pontos vizinhos a P.
Copyright 2003. All rights reserved.
Condições de Contorno
Um ponto ainda não discutido é o que diz respeito aos VC’s nas fronteiras do domínio de cálculo, ou seja, as condições de contorno. Todas as equações até aqui deduzidas foram para os volumes internos. Há algumas maneiras de se discretizar a fronteira do domínio. Uma delas é colocar o ponto central do VC exatamente na fronteira do domínio (Patankar, 1980), mas este método gera volumes não-inteiros e viola a conservação da propriedade quando o valor desta é conhecida ( prescrita ) na fronteira. Podem-se ainda imaginar volumes fictícios, entretanto, se a malha é muito refinada e possui mais de uma dimensão, o número de VC’s gerados aumenta o esforço computacional sobremaneira. O método indicado por Maliska (1995), tem uma melhor solidez física e pode ser aplicado a qualquer tipo de sistemas coordenados. O procedimento indicado é o mesmo para os volumes internos: fazer um balanço da propriedade para os volumes de fronteira se utilizando das condições de contorno. Este procedimento não aumenta o número de VC’s e é fisicamente mais coerente. Ainda no modelo anterior, nós poderemos ter na fronteira:
¾ Propriedade prescrita:
f f f f P x q k T T ∆
onde Tf , é a temperatura na fronteira e tem um valor conhecido.
¾ Fluxo prescrito "
q f = valor conhecido (20)
¾ Convecção: ( ) f f f f f P x q hT T k T T ∆
( (^) f ) f
f (^) f T T k
h x q h − = (^1) + ∆ ∞
Copyright 2003. All rights reserved.
Método Gauss-Seidel
O método Gauss-Seidel pertence ao grupo de métodos ponto a ponto, onde se resolve o sistema linear através da resolução de equação por equação utilizando-se, dentro de um mesmo ciclo iterativo, os valores das variáveis já calculadas neste ciclo. Isto acelera a convergência. O algoritmo deste método, é o seguinte: ¾ Estimar o campo inicial da variável; ¾ Calcular φ (^) P pela Eq.(26). A varredura é feita de baixo para cima, do sul para o norte e de oeste para leste. Os valores da propriedade φ W , φ S e φ B podem ser consideradas conhecidas para a mesma varredura; ¾ Checar se houve convergência; ¾ Caso o critério não tenha sido satisfeito, retornar.
A (^) P φ (^) Pk +^1 =∑ Anb φ NBk +^1 + Ae φ Ek + An φ Nk + At φ Tk + B (26)
Esse método é de lenta convergência e requer uma matriz com dominância diagonal. Em contrapartida é de fácil implementação e é totalmente vetorizável.
Método TDMA (TriDiagonal Matrix Algorithm)
Esse método é conhecido como método linha por linha, resolvendo diretamente uma linha, ou seja, um problema unidimensional. A resolução de problemas bi e tridimensionais são iterativas, com a varredura se processando linha por linha e coluna por coluna. Como exemplo, considera-se um problema bidimensional no qual a equação a ser resolvida será:
A (^) P TP = AeTE + AwTW + AnTN + AsTS + B (27)
Copyright 2003. All rights reserved.
Escrevendo a Eq. (27) de forma mais conveniente para procedimento recursivos, temos:
Am Tm + BmTm + 1 + CmTm − 1 = D m (28)
Deseja-se determinar uma relação recursiva da forma:
Tm = PmTm + 1 + Q m (29)
para que, usando as condições de contorno, possa-se varrer uma linha em um sentido, determinando os coeficientes P e Q , e voltar, determinando os valores da variável T. Baixando um índice da Eq. (29) temos:
Tm (^) − 1 = Pm − 1 Tm + Qm − 1 (30)
Substituindo a Eq.(30) na Eq. (28) e comparando o resultado com a Eq. (29), encontramos as seguintes expressões para os coeficientes P e Q :
m m m m m A C P
P B (31)
1
1 −
−
m m m m m m m A C P
Q D C Q (32)
As Eqs. (31) e (32) permitem, depois de conhecidos P 1 e Q 1 , determinar todos os valores de P e Q. A determinação de P 1 e Q 1 é fácil de inferir, inspecionando as Eqs. (31) e (32) e considerando que C 1 é zero. Neste caso temos as seguintes equações.:
1 1 1 A
P = − B (33)
1 1 1 A
Q = D (34)
Copyright 2003. All rights reserved.
Quando α>1 temos sobre-relaxação e quando α<1temos sub-relaxações. A sua escolha só obedece a um critério: Tente alguns valores e verifique os que respondem melhor ao caso.
A Eq.(35) apresenta certa similaridade com a Eq.(17) mas não são iguais! A Eq.(35) é um artifício numérico puramente iterativo, enquanto a Eq.(17) descreve o comportamento de uma propriedade dentro de um intervalo de tempo.
Funções de Interpolação
Quando discretizamos o nosso modelo, tivemos a necessidade de usar funções de interpolação que descrevessem o comportamento da propriedade e de suas n derivadas nas faces dos VC’s localizadas entre os nós que formam a malha, tanto para o tempo como para o espaço. E, para o espaço, utilizamos uma função linear dada por diferenças centrais mas poderíamos ter usado qualquer outra função, linear ou não. O uso de funções de interpolação adequada varia de acordo com o fenômeno que se estuda. Às vezes funções lineares não são boas aproximações para descrever fenômenos que se comportam não- linearmente. Esse tipo de problema ocorre principalmente quando se tem a presença dos fenômenos de convecção e de difusão concomitantemente. Quando a velocidade de um escoamento é muito alta, a influência da difusão não é percebida, e a propriedade é “arrastada” pelo domínio quase sem sofrer alterações. Por outro lado, quando o inverso se dá, a propriedade varia linearmente pelo domínio. O problema está quando os dois fenômenos se somam: a aproximação linear para a função de interpolação não será a melhor opção. Pois o perfil da propriedade se encontra em algum lugar entre os dois extremos. A Figura 7 esclarece:
Copyright 2003. All rights reserved.
A : Convecção B : Difusão C : Ambos os fenômenos
Figura 7- Comportamento idealizado da propriedade φ perante os fenômenos de convecção e difusão.
Os problemas que podem aparecer com a interpolação inadequada são a Oscilação Numérica e a Difusão Numérica, apresentadas na Figura 8.
Figura 8 – Tipos de desvios numéricos para uma variação em degrau da propriedade
Algumas funções de interpolação são apresentadas:
Copyright 2003. All rights reserved.