Dinâmica dos Fluidos Computacional

Prévia do material em texto

Universidade Federal de Campina Grande
Centro de Ciência e Tecnologia
Unidade Acadêmica de Engenharia Mecânica
Departamento de Engenharia Mecânica
Agenda Semanal 07
Dinâmica dos Fluidos Computacional
(Capítulo 15, Mecânica dos Fluidos, Merlee C. Potter, 3ª edição)
Disciplina
Mecânica dos Fluidos II
Aluna
Joyce Ingrid Venceslau de Souto
Campina Grande - PB
Setembro de 2021
Sumário
1.	INTRODUÇÃO 	 	 3
2. 	VISÃO GERAL DOS MÉTODOS DE DIFERENÇAS FINITAS E DE VOLUMES FINITOS 4
3.	EXEMPLOS DOS MÉTODOS SIMPLES DE DIFERENÇAS FINITAS	 5
	3.1 Discretização do Domínio 5
	3.2 Discretização da Equação Governante 6
	3.3 Definição do Algoritmo de Solução 11
	3.4 Comentários sobre a Escolha dos Operadores de Diferenças 14
4.	EXEMPLOS DOS MÉTODOS SIMPLES DE VOLUMES FINITOS 21
1. INTRODUÇÃO 
Entende-se que a Dinâmica dos Fluidos Computacional (DFC) – ou, do inglês, Computational Fluid Dynamics (CFD) – começou a ser desenvolvida tendo em vista a necessidade de se resolver as equações que regem o escoamento fluido de forma que estas não permitiam o emprego dos métodos analíticos convencionais. Sabe-se que essas equações, das quais se destacam: a da continuidade, a de Navier-Stokes e a da conservação da energia. O conjunto formado por elas compõe um sistema de EDPs (Equações Diferenciais Parciais) acopladas e não-lineares. 
Frente a isso, infere-se que a impossibilidade da obtenção de soluções analíticas para problemas de escoamentos fluidos com esse sistema se deve à presença dos termos não-lineares. Para um punhado de problemas, é razoável realizar considerações acerca desses termos de forma a torná-los nulos ou desprezíveis, seja devido à própria física do problema, ou porque são pequenos em comparação aos demais termos do sistema. Contudo, para a maior parte das aplicações em engenharia, isso não acontece e, sendo assim, é necessário recorrer à métodos numéricos assistidos por computador, para encontrar as soluções das EDPs, sendo este o principal objetivo do estudo da DFC. Já que se tratam de métodos que se utilizam de computador, o grande desenvolvimento em termos de memória e processamento de dados deste permitiu aos métodos numéricos uma ampla gama de soluções quanto ao escoamento em fluidos, podendo ser compressível/incompressível, laminares/turbulentos, reagentes quimicamente ou não, com solução homogênea ou não.
Figura 1 – Análise Numérica de Escoamentos utilizando CFD.
Fonte: https://www.esss.co/iesss/especializacao/analise-numerica-escoamentos-utilizando-cfd-online/
2. VISÃO GERAL DOS MÉTODOS DE DIFERENÇAS FINITAS E DE VOLUMES FINITOS
Em contrapartida ao que foi descrito no final da seção anterior, observa-se que existem certas limitações computacionais, que cercam às soluções encontradas e que norteiam o processo de elaboração do algoritmo, e das quais se podem destacar quatro, sendo elas: computadores são capazes de efetuar apenas operações aritméticas e lógicas, logo, as derivadas e integrais devem ser apresentadas em termos daquelas; computadores representam números por uma quantidade finita de dígitos, logo, sempre existirá um erro de arredondamento associado a essa transcrição e que deve ser controlado; computadores apresentam memória de armazenamento finita e limitada, então as soluções serão aproximadas a uma quantidade limitada de pontos espaciais e temporais; e os computadores excutam uma quantidade limitada de operações por unidade de tempo. 
A partir das quatro limitações associadas ao computador descritas anteriormente, é possível inferir que o processo de solução não deve apenas objetivar alcançá-la, mas o fazer otimizando ambos o tempo de processamento e a quantidade de operações necessários. Com base nisso, os métodos das Diferenças Finitas (DF) e dos Volumes Finitos (VF) solucionam as EDPs através das seguintes etapas:
1. Discretização do Domínio: como o domínio dos problemas reais é contínuo em tempo e espaço, deve-se substituí-lo por um segmentado em pontos de uma grade e níveis temporais. Essa etapa realizada idealmente envolve utilizar o menor número de pontos e níveis temporais necessário para obter a solução com a exatidão desejada.
2. Discretização das EDPs: frente à limitação computacional e a mudança no domínio, deve-se substituir as EDPs por um sistema de equações algébricas para o domínio discretizado. Se ele for representado de forma idealizada, esse sistema possui a mesma descrição física que as EDPs regentes do escoamento fluido.
3. Especificação do Algoritmo: A partir dos passos anteriores, pode-se estabelecer o procedimento por etapas através do qual as equações algébricas, seja de DF ou de VF, irão resultar em soluções a medida em que se avança a cada ponto na grade e nível temporal.
3. EXEMPLOS DOS MÉTODOS SIMPLES DE DIFERENÇAS FINITAS 
Para visualizar o caso do MDF, em primeiro lugar, é preciso considerar o escoamento como sendo não-permanente, incompressível e laminar cujo fluido com uma dada viscosidade escoa entre duas placas paralelas. No instante inicial, as placas e o fluido estão estacionários, contudo, após um dado instante temporal t, a placa inferior começa a se mover para a direção do eixo axial positiva em um movimento retilíneo uniforme, de acordo com a Fig. 2. 
Figura 2 – Escoamento entre placas paralelas, com movimento relativo.
Fonte: http://mecflu2.usuarios.rdc.puc-rio.br/Grad_Eng1707/10-MecanicaFluidosII-IntroducaoCompressivel.pdf
Assim sendo, pode-se inferir que as equações governantes desse caso de escoamento fluido são as da conservação da massa e a de Navier-Stokes. Tendo em vista as considerações já estabelecidas, tem-se que essas equações resultam numa EDP linear, dada por: 
(1)
Além disso, as condições iniciais e de contorno da Eq. 1 são 
(2)
3.1 Discretização do Domínio
Primeiramente, a partir da Fig. 2, é notório destacar que o domínio espacial é um segmento de linha partindo de 0 até H e que o domínio temporal é um raio contínuo cuja fonte é o instante inicial mas, apesar disso, a duração de interesse do escoamento é finita e o instante final será aquele no qual se atinge o estado estacionário. No ponto associado a esse instante, enquanto o domínio espacial é discretizado através da sua substituição por JL pontos de grade fixos uniformemente distribuídos, o domínio temporal é substituído por níveis temporais incrementados de forma constante. Sendo este apenas mais um tipo de discretização, ambos os domínios não precisam necessariamente sofrem acréscimo uniforme ou serem estacionários em todos os casos e métodos. 
Figura 3 – Sistema de grade e níveis temporais para o escoamento da Fig. 2.
Fonte: http://mecflu2.usuarios.rdc.puc-rio.br/Grad_Eng1707/10-MecanicaFluidosII-IntroducaoCompressivel.pdf
De forma geral, o tamanho desses espaçamentos da grade entre pontos e do passo entre os níveis temporais é função das escalas de comprimento e tempo a serem resolvidas e das propriedades das EDPs discretizadas. Assim sendo, tem-se que a solução obtida se restringe somente aos pontos de grade e níveis temporais declarados, sendo indicada por
 A partir do qual os sobrescritos indicam os níveis temporais e os subscritos sinalizam os locais dos pontos de grade.
3.2 Discretização da Equação Governante
Uma vez que se tenha realizado o processo de discretização do domínio, a presente etapa consiste em converter as EDPs governantes do problema em equações algébricas cujas entradas serão os pontos de grade e níveis temporais discretizados na etapa anterior. No caso do MDF, substitui-se as derivadas por operadores de diferenças,de forma que, a discretização das equações envolve, primeiramente, obter esses operadores e, em seguida, selecioná-los. 
1. Obtenção dos Operadores de Diferenças
Dentre uma vasta gama de métodos para obtenção dos operadores de diferenças, tem-se que, para problemas cujas soluções não envolvem descontinuidades, um método baseado no seguinte teorema que diz que, se uma função e suas derivadas são contínuas, a um valor único e finitas, logo o valor dessa função pode ser expresso em função de e suas derivadas, em qualquer ponto, através de uma expansão em série de Taylor, desde que o outro ponto em análise esteja no campo de convergência da série. Com esse teorema, obtém-se os quatro operadores de diferenças a seguir:
(3)
(4)
(5)
(6)
Nestes, e indicam erros de truncamento, associados aos termos suprimidos da expansão em séries de Taylor. Enquanto nas Eqs. 3 e 4 a potência de em é 1, pois o termo à esquerda dos erros de truncamento para esses dois operadores está sendo multiplicado por aquele elevado à primeira potência, nas Eqs. 5 e 6 a potência de em é 2 por causalidade análoga. Para esse último caso, diz-se que esses operadores de diferenças possuem precisão de segunda ordem, o que é um valor adequado de forma geral, visto que essa precisão indica a quantidade de termos truncados na expansão e, além disso, expressa adequadamente uma função suave (ou sem descontinuidades).
Além disso, os pontos de grade formam uma estrutura subjacente ao operador de diferenças, de maneira que se ela envolver pontos com índices maiores ou iguais a então o operador será denominado de operador de diferenças à frente, como na Eq. 3; caso envolva índices menores ou iguais a , ele será denominado de operador de diferenças a trás, assim como na Eq. 4; se o número de pontos atrás e a frente de for igual, ele será chamado de operador de diferenças central como nas Eqs. 5 e 6. Ademais, se o número de pontos atrás e a frente de for diferente, ele será chamado de operador de diferenças que tende para trás ou à frente.
De maneira geral, as séries de Taylor podem obter operadores de diferenças de qualquer derivada de qualquer ordem de precisão, usando qualquer estrutura. Como exemplo, pode-se considerar um operador de diferenças para uma derivada de primeira ordem, com precisão de terceira ordem e uma estrutura que possui apenas um ponto da grade a jusante. A obtenção desse operador segue a seguinte sequência: determinar o número de termos a serem mantidos na série de Taylor, sendo este dado pela soma da ordem da derivada para qual se busca o operador com a ordem de precisão requerida (neste caso, ); selecionar uma estrutura com pontos para operador no ponto (neste caso, ); construir séries de Taylor em torno de para todos os pontos do modelo, a exceção do ponto ; resolver as equações lineares acopladas para o operador de diferenças procurado. 
Na Tab. 1 abaixo, encontra-se um resumo de operadores de diferenças comumente utilizados. Nota-se que se e forem substituídos por e , os operadores também poderão ser utilizados para derivadas temporais.
Tabela 1 – Escoamento entre placas paralelas, com movimento relativo.
Fonte: http://mecflu2.usuarios.rdc.puc-rio.br/Grad_Eng1707/10-MecanicaFluidosII-IntroducaoCompressivel.pdf
2. Seleção dos Operadores de Diferenças 
Para selecionar o operador de diferenças, aplica-se a Eq. 1 a um ponto arbitrário interno da grade e em um tempo entre e , de forma que 
(7)
Onde algum valor no intervalo, ou igual a, e . Utilizando a Eq. 7, supõe-se que a solução é procurada no tempo e conhecida em . Essa consideração geralmente é aplicável pois sempre se conhece a solução no nível temporal anterior e, no caso do nível temporal zero, a solução é a condição inicial. 
3. Método Explícito
Para a Eq. 7, se , logo todas as equações são calculadas no tempo (nível temporal anterior). Supondo que, num processo de tomada de decisão, escolhe-se o operador de diferenças à frente para e o operador central para , a Eq. 7 se torna em
(8a)
A EDF acima tem precisão espacial de segunda ordem e temporal de primeira ordem. A partir da Eq. 8a, resolve-se para , obtendo
(8b)
A EDF numa das formas das Eqs. 8 podem ser aplicadas em qualquer ponto interno arbitrário da grade, enquanto que as EDFs nos pontos de contorno da grade são obtidas através das condições de contorno (CCs). Para este caso, as condições de contorno dadas pela Eq. 2 produzem as EDFs a seguir:
(9a)
(9b)
É notório que as Eqs 8 e 9 apresentam apenas uma incógnita em consequência de se considerar , logo, um método que apresente apenas essas EDFs é chamado de explícito. Caso particular, denominado de esquema explícito de Euler, consiste em uma diferenciação à frente no tempo com precisão de primeira ordem. 
4. Método Implícito
Já se na Eq. 7, as derivadas espaciais serão calculadas no nível temporal cuja solução é desconhecida (). Ao escolher o operador de diferenças para trás para e o operador central para , a Eq. 7 resulta em
(10)
Podendo ser rearranjada em outra equação e possui a mesma analogia com relação a ordem de precisão da Eq. 8ª. A diferença entre ambas se constitui no fato de que a Eq. 10 possui 3 incógnitas ao invés de apenas uma EDF, e ela é contextualizada pelo cálculo das derivadas espaciais em . A plicando a outra forma da Eq. 10 em cada ponto interno da grade, obtém-se o seguinte sistema de equações:
(11a)
(11b)
Em resumo, o método implícito utiliza uma EDF com mais de 1 incógnita, inferindo a necessidade de solução das equações de forma simultânea. Esse método particular descrito, que consiste na diferenciação no tempo para trás com precisão de primeira ordem é denominado de esquema implícito de Euler.
5. Método Generalizado
Agora, supondo que , a Eq. 7 pode resultar em
(12)
Onde as derivadas espaciais são iguais às temporais. Caso se escolha o operador central, dado pela Eq. 5, para , a Eq. 12 se torna em
(13)
Para , a Eq. 13 se reduz ao esquema obtido pelo método explícito e, para , a equação acima se reduz ao esquema implícito. Tem-se um caso particular, cuja fórmula de diferenciação temporal resultante é denominada de método trapezoidal ou de Crank-Nicolson, aplicado a , da onde a Eq. 13 terá precisão de segunda ordem no espaço e tempo. 
3.3 Definição do Algoritmo de Solução
A partir do domínio e equações governantes discretizados, resta especificar processo de obtenção das soluções, ou seja, o seu algoritmo de solução. De acordo com os conjuntos de EDFs obtidos anteriormente, sabe-se que há três algoritmos distintos de solução. Tendo em vista o método explícito, expresso pelas Eqs. 8 e 9, o algoritmo de solução é obtido ao realizar as seguintes etapas:
1. Especificar o problema, inserindo os valores para , ;
2. Especificar o número de pontos desejados, inserindo JL;
3. Especificar o passo temporal e a duração desejada ;
4. Calcular o espaçamento de grade com ;
5. Calcular as constantes e ;
6. Ajustar o contador de nível de tempo para 0;
7. Especificar a solução em cada ponto da grade no tempo utilizando a condição inicial dada pela Eq. 2;
8. Calcular em cada ponto interno da grade usando a Eq. 8b;
9. Calcular em cada ponto do contorno da grade usando as Eqs. 9;
10. Gravar a solução no nível de tempo para análise de dados;
11. Incrementar o contador de nível temporal a cada laço no valor da unidade;
12. Se , repetir os passos 7-12.
Considerando o método implícito, basta repetir as etapas descritas anteriormente, a exceção do passo 8 que deve ser substituído pelo seguinte:
8. Calcular em cada ponto interno da grade resolvendo o sistema linear de equações, usando as Eqs. 11;
Para as EDFs dos pontos de contorno e do método generalizado, o algoritmo depende do valor de . Caso seja igual a 0 ou a 1, os algoritmos são idênticos aos obtidos através das etapas descritas anteriormente; caso , o algoritmo será aquele referente ao método implícito, com o detalhe da substituição dos elementos na matriz A e no vetor b na Eqs. 11 por aqueles referentes à Eq. 13. Além disso, a solução exata pode ser utilizada tantopara o cálculo da precisão dos métodos explícitos e implícito de Euler e o de Crank-Nicolson, quanto para fornecer indícios para o espaçamento da grade e do tamanho de nível temporal necessário para encontrar uma solução com dada exatidão. 
Ademais, é do interesse do presente apresenta rum método particular bastante eficiente, denominado de algoritmo de Thomas, no que tange à solução de sistemas lineares de equações com matrizes de coeficientes tridiagonais, sendo estes comumente encontrados na DFC. Para descrever este, tem-se um sistema de equações lineares independentes a seguir:
(14a)
Onde,
(14b)
Nas Eqs. acima, os elementos de A e b são conhecidos e se busca os elementos de x. Primeiramente, fatora-se a matriz de coeficientes A em duas matrizes bidiagonais, assim
(15)
(16)
Agora, para determinar os elementos de L e U, iguala-se a Eq. 16 com a 14b e se comparam os termos. Esse procedimento fornece uma fórmula recursiva para o cálculo dos elementos de L e U. O próximo passo é a substituição da fatoração de LU de A na Eq. 14a, seguindo com a solução para z pela substituição à frente:
(17)
E a solução para x pela substituição para trás:
(18)
Logo, conclui-se o algoritmo de Thomas, que demanda operações aritméticas para obter uma solução.
3.4 Comentários sobre a Escolha dos Operadores de Diferenças
A partir dos detalhamentos da presente seção, pode-se perceber que existe uma variedade de operadores de diferenças que podem ser utilizadas para cada derivada de modo que a EDF não é unicamente associada a uma dada EDP. Caso a maior ordem de derivada temporal for 1 (caso das equações da continuidade, Navier-Stokes e da energia), as derivadas espaciais de segunda ordem representam a difusão e as de primeira ordem representam a convecção: a difusão espalha a perturbação em todas as direções, de modo que as derivadas espaciais de segunda ordem são substituídas pelos operadores de diferenças centrais, enquanto que a convecção é direcional então é conveniente o uso de operadores ascendente ou com essa tendência para a substituição das derivadas espaciais de primeira ordem.
Sabe-se que o operador que substitui a derivada temporal determina se as EDFs serão acopladas entre si, ou não, para uma EDP em pontos diferentes. Sendo assim, um operador de diferenças temporais generalizado com 3 níveis é dado por:
(19)
Onde e são constantes determinadas pelo usuário. Além de conter a Eq. 13 no caso especial , ela tem precisão de segunda ordem quando e de primeira ordem em caso contrário. Outros operadores desse tipo podem ser obtidos quando se escolhe valores convenientes para e : para , tem-se Euler explícito; para , tem-se Euler implícito; para , tem-se Crank-Nicolson; para , tem-se para trás três pontos.
De uma forma geral, os operadores presentes nas Eqs. 12 e 19 s classificados como operadores de estágio único, sendo utilizados em problemas de PVI (Problemas de Valor Inicial) com a integração das EDOs (Equação Diferencial Ordinária), com o detalhe de que qualquer método pode ser utilizado na solução de PVIs, podendo ser usado ou generalizado como operador de diferenças para as derivadas temporais nas EDPs. Desse odo, pode-se utilizar qualquer um dos operadores de estágio único dados pelas fórmulas de Adams-Bashfort ou de Adams-Moulton. As fórmulas de Adams-Bashfort são explícitas e dadas por:
(20)
 Os coeficientes , para uma precisão de até quarta ordem, são: para precisão de primeira ordem, e (fórmula explícita de Euler); para precisão de segunda ordem, , e ; para precisão de terceira ordem, , , e ; para precisão de quarta ordem, , , , e .
As fórmulas de Adams-Moulton são implícitas e podem ser expressas por:
(21)
Os coeficientes , para uma precisão de até quarta ordem, são: para precisão de primeira ordem, e (fórmula implícita de Euler); para precisão de segunda ordem, , e (fórmula implícita de Crank-Nicolson); para precisão de terceira ordem, , , e ; para precisão de quarta ordem, , , , e .
Ao fazer uso das fórmulas descritas anteriormente na presente seção, nota-se que as fórmulas de ordem mais alta não são “auto-iniciáveis”, ou seja, elas precisam de fórmulas de baixa ordem para iniciar os cálculos. Além disso, deve-se perceber que quanto mais alto for a ordem de precisão, maior será o número de níveis temporais envolvidos, o que pode aumentar notoriamente a demanda de memória e processamento de dados do computador. 
 Paralelamente a isso, existem operadores alternativos de integração de EDO e PVIs aos mostrados há pouco, que são os de múltiplos estágios, como os de Runge-Kutta, sendo as suas fórmulas construídas a partir da interpolação e extrapolação das derivadas calculadas previamente. Vale destacar que todas elas são explícitas e auto-iniciáveis, independentemente da ordem de precisão. Um desses operadores é dado por:
(22a
(22b)
Em síntese, demonstrou-se que, além da grande variedade, o operador de diferenças deve ser escolhido com base nas soluções a serem procuradas: se são, ou não, estacionárias. No caso de solução estacionária, não existe preocupação quanto a precisão temporal, então, o operador deve ser escolhido para acelerar a convergência para o estado desejado. Caso contrário, pode-se utilizar operadores implícitos ou explícitos, sendo a ordem de precisão temporal normalmente igual a 2 para operadores de único estágio e esse número pode ser diferente para operadores de diversos estágios. 
4. EXEMPLOS DOS MÉTODOS SIMPLES DE VOLUMES FINITOS
Para entender sobre o método dos VF (Volumes Finitos), é mais conveniente entender que a sua diferença, comparado ao método das DF, está na interpretação das soluções nos pontos da grade: no MDF, a solução é visualizada como uma função pontual, então pode ser interpolada em qualquer ponto a partir dos pontos na grade e dos níveis temporais; no MVF, a solução é visualizada como o valor médio da função em uma célula. Além disso, existem duas opções de abordagem para seleção das células, destacando que, independente desta, as células adjacentes dividem os mesmos contornos de células sem sobreposição ou lacunas entre elas:
1. Abordagem de vértice centrado: Seleciona-se o contorno da célula como bissetor de uma linha que liga pontos adjacente da grade;
2. Abordagem de célula centrada: Escolhe-se os contornos da célula como os pontos da grade em 1D, as linhas que conectam esses pontos em 3D, e os planos constituídos por essas linhas em 3D. 
Figura 4 – Célula de vértice centrado e célula centrada
Fonte: http://mecflu2.usuarios.rdc.puc-rio.br/Grad_Eng1707/10-MecanicaFluidosII-IntroducaoCompressivel.pdf
A partir da Fig. 4, consegue-se visualizar a diferença entre as abordagens e, na Fig. 5, é ilustrado a diferença entre as interpretações de DF e VF.
2