GUSTAVO TAIJI NAOZUKA

SIMULAÇÃO NUMÉRICA DO PROBLEMA DE DISPERSÃODE POLUENTES NA ATMOSFERA

LONDRINA–PR

2016

GUSTAVO TAIJI NAOZUKA

SIMULAÇÃO NUMÉRICA DO PROBLEMA DE DISPERSÃODE POLUENTES NA ATMOSFERA

Trabalho de Conclusão de Curso apresentadoao curso de Bacharelado em Ciência da Com-putação da Universidade Estadual de Lon-drina para obtenção do título de Bacharel emCiência da Computação.

Orientador: Prof. Dr. Alan Salvany FelintoCoorientador: Profa. Dra. Neyva Maria LopesRomeiro

LONDRINA–PR

2016

Gustavo Taiji NaozukaSimulação numérica do problema de dispersão de poluentes na atmosfera/

Gustavo Taiji Naozuka. – Londrina–PR, 2016-92 p. : il. (algumas color.) ; 30 cm.

Orientador: Prof. Dr. Alan Salvany Felinto

– Universidade Estadual de Londrina, 2016.

1. Dispersão de poluentes. 2. Poluição atmosférica. 3. Simulação numérica. 4.Equações diferenciais. I. Prof(a). Dr(a). Alan Salvany Felinto. II. UniversidadeEstadual de Londrina. III. Faculdade de Ciência da Computação. IV. Simulaçãonumérica do problema de dispersão de poluentes na atmosfera

CDU 02:141:005.7

GUSTAVO TAIJI NAOZUKA

SIMULAÇÃO NUMÉRICA DO PROBLEMA DE DISPERSÃODE POLUENTES NA ATMOSFERA

Trabalho de Conclusão de Curso apresentadoao curso de Bacharelado em Ciência da Com-putação da Universidade Estadual de Lon-drina para obtenção do título de Bacharel emCiência da Computação.

BANCA EXAMINADORA

Prof. Dr. Alan Salvany FelintoUniversidade Estadual de Londrina

Orientador

Prof. Dr. Cosmo Damião SantiagoUniversidade Tecnológica Federal do Paraná

Prof. Dr. Elieser Botelho Manhas Jr.Universidade Estadual de Londrina

Londrina–PR, 26 de Fevereiro de 2016

Dedico este trabalho à minha família e aos meus amigos.

AGRADECIMENTOS

Em um primeiro momento, gostaria de agradecer à minha família, em especial aosmeus pais, Marcelo e Julia, e às minhas irmãs, Elaine e Juliana, por me auxiliarem emalguns trabalhos e estudos e sempre me apoiarem nas escolhas que fiz desde o início daminha vida acadêmica. Agradeço a eles por terem me ajudado a conquistar tudo o que jáconquistei, além de outras conquistas que ainda estão por vir.

Agradeço também aos meus amigos, principalmente àqueles que conheci duranteo período da graduação, por terem encarado as dificuldades desse período, como provas etrabalhos, juntos a mim. Apesar de algumas discussões ocorridas, todos acabaram se aju-dando e cooperando para enfrentar os problemas. Agradeço a eles por terem me ensinadoalgo muito importante, o trabalho em grupo, o qual não se constituía ainda em uma dasminhas qualidades.

Além disso, gostaria de agradecer aos meus orientadores do Departamento de Ma-temática, Neyva e Eliandro, que me convidaram a trabalhar com eles enquanto aindaera calouro e me forneceram importantes ensinamentos desde então. Com certeza, grandeparte do que continuarei trabalhando estará embasado nesses ensinamentos. Agradeçotambém ao meu orientador do Departamento de Computação, Alan, o qual, apesar de teriniciado um trabalho com ele somente este ano, me concedeu seus conhecimentos sobre oTrabalho de Conclusão de Curso e o Estágio Curricular Obrigatório.

Finalmente, agradeço à Fundação Araucária, por meio do convênio 24/2012 - Pes-quisa Básica e Aplicada - em relação ao projeto contemplado de número 39225, pelosuporte financeiro.

“Ciência é uma equação diferencial.Religião é a condição de contorno.”

(Alan Turing)

NAOZUKA, G. T.. Simulação numérica do problema de dispersão de poluen-tes na atmosfera. 92 p. Trabalho de Conclusão de Curso (Bacharelado em Ciência daComputação) – Universidade Estadual de Londrina, Londrina–PR, 2016.

RESUMO

Os impactos ambientais e na saúde humana provocados pela emissão de poluentes naatmosfera tornaram-se uma temática bastante relevante no contexto atual. A preservaçãoe a melhoria da qualidade do ar consistem em uma questão de extrema importância ea ser levada em consideração. Por esse motivo, a realização de simulações numéricas deequações diferenciais proporcionou a criação de uma técnica interessante e eficiente dese avaliar e prever tais impactos. Nesse sentido, este trabalho tem como objetivo sugerirum modelo numérico capaz de solucionar a equação que governa o fluxo de poluentesna atmosfera. O modelo proposto apresenta, como inovação, a representação no sistemade coordenadas generalizadas bidimensionais combinada à discretização por meio do mé-todo de diferenças finitas, utilizando o esquema upwind de primeira ordem para a suaresolução. Para a comprovação teórica do funcionamento do modelo, três problemas, comdiferentes configurações geométricas, são simulados. De acordo com os experimentos, re-sultados bastante satisfatórios foram encontrados, condizendo com a explicação teóricados processos dispersivos.

Palavras-chave: Dispersão de poluentes. Poluição atmosférica. Simulação numérica.Equações diferenciais.

NAOZUKA, G. T.. Numerical simulation of the pollutants dispersion problemin the atmosphere. 92 p. Final Project (Bachelor of Science in Computer Science) –State University of Londrina, Londrina–PR, 2016.

ABSTRACT

The environmental and human health impacts caused by the emission of pollutants in theatmosphere have become a very relevant issue in the current context. The preservationand the improvement of the air quality are an extremely important and to be taken intoaccount matter. For this reason, the realization of numerical simulations of differentialequations provided creating an interesting and efficient technique to assess and predictsuch impacts. In this sense, this paper aims to suggest a numerical model able to solvethe equation that governs the flow of pollutants in the atmosphere. The proposed modelpresents, as innovation, the representation in the system of two-dimensional generalizedcoordinates combined with the discretization using the finite difference method, apply-ing the first order upwind scheme for its resolution. For the theoretical verification ofthe operation of the model, three problems, with different geometric configurations, aresimulated. According to the experiments, satisfactory results were found, matching thetheoretical explanation of dispersive processes.

Keywords: Pollutant dispersion. Atmospheric pollution. Numerical simulation. Differen-tial equations.

LISTA DE ILUSTRAÇÕES

Figura 1 – Mecanismos em um processo de dispersão de poluentes - retirada deStockie [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

Figura 2 – Classificação de malhas quanto à distribuição espacial dos pontos. . . . 36Figura 3 – Sistemas de coordenadas . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figura 4 – Rotulagem de índices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figura 5 – Discretização por Diferenças Finitas . . . . . . . . . . . . . . . . . . . 42Figura 6 – Disposição dos nós 𝐷, 𝑈 e 𝑅 . . . . . . . . . . . . . . . . . . . . . . . 43Figura 7 – Elementos de malha - retirada de Barba [2] . . . . . . . . . . . . . . . 56Figura 8 – Fluxograma do modelo de dispersão de poluentes na atmosfera . . . . . 62Figura 9 – Malha fantasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figura 10 – Malha deslocada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figura 11 – Malha do duto simples . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figura 12 – Campo de velocidade do duto simples . . . . . . . . . . . . . . . . . . . 76Figura 13 – Concentrações do poluente no duto simples (caso 1) . . . . . . . . . . . 76Figura 14 – Concentrações do poluente no duto simples (caso 2) . . . . . . . . . . . 77Figura 15 – Concentrações do poluente no duto simples (caso 3) . . . . . . . . . . . 78Figura 16 – Concentrações do poluente no duto simples (caso 4) . . . . . . . . . . . 79Figura 17 – Malha do duto com um obstáculo . . . . . . . . . . . . . . . . . . . . . 80Figura 18 – Campo de velocidade do duto com um obstáculo . . . . . . . . . . . . 80Figura 19 – Concentrações do poluente no duto com um obstáculo . . . . . . . . . 81Figura 20 – Malha do duto com dois obstáculos . . . . . . . . . . . . . . . . . . . . 82Figura 21 – Campo de velocidade do duto com dois obstáculos . . . . . . . . . . . . 83Figura 22 – Concentrações do poluente no duto com dois obstáculos . . . . . . . . 84

LISTA DE TABELAS

Tabela 1 – Parâmetros para a geração da malha computacional . . . . . . . . . . . 73Tabela 2 – Parâmetros para a resolução numérica da Equação de Navier-Stokes . . 73Tabela 3 – Parâmetros para a resolução numérica da equação de transporte . . . . 74Tabela 4 – Tempos de processamento dos experimentos realizados . . . . . . . . . 85

LISTA DE ALGORITMOS

4.1 Modelo de dispersão de poluentes na atmosfera . . . . . . . . . . . . . . . 634.2 Geração da malha fantasma . . . . . . . . . . . . . . . . . . . . . . . . . . 654.3 Inicialização do campo de velocidade na malha fantasma . . . . . . . . . . 664.4 Atribuição das condições iniciais . . . . . . . . . . . . . . . . . . . . . . . . 664.5 Atribuição das condições de fronteira . . . . . . . . . . . . . . . . . . . . . 674.6 Geração da malha deslocada . . . . . . . . . . . . . . . . . . . . . . . . . . 674.7 Inicialização da matriz de Jacobianos . . . . . . . . . . . . . . . . . . . . . 684.8 Resolução numérica da equação de transporte . . . . . . . . . . . . . . . . 694.9 Funções upwind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.10 Arquivo Make . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

LISTA DE ABREVIATURAS E SIGLAS

VOC Compostos Orgânicos Voláteis

CFD Computational Fluid Dynamics

EDO Equação Diferencial Ordinária

EDP Equação Diferencial Parcial

ET Erro de Truncamento

FOU First Order Upwinding

SOU Second Order Upwinding

QUICK Quadratic Upstream Interpolation for Convective Kinematics

CUBISTA Convergent and Universally Bounded Interpolation Scheme for the Tre-atment of Advection

CTAG Comprehensive Turbulent Aerosol Dynamics and Gas Chemistry

RANS Reynolds-averaged Navier-Stokes

LES Large-Eddy Simulation

DES Detached-Eddy Simulation

CTM Chemistry Transport Model

CNTP Condições Normais de Temperatura e Pressão

LISTA DE SÍMBOLOS

𝜌 Massa específica do fluido

𝑢 = 𝑣𝑥 Componente do vetor velocidade em 𝑥

𝑣 = 𝑣𝑦 Componente do vetor velocidade em 𝑦

𝑡 Tempo em coordenadas cartesianas

𝑥, 𝑦 Posicionamento em coordenadas cartesianas

𝑝 Pressão em 𝑥 e 𝑦

𝜇 Viscosidade dinâmica do fluido

𝑅1 Termo fonte na direção 𝑥 da Equação de Navier-Stokes

𝑅2 Termo fonte na direção 𝑦 da Equação de Navier-Stokes

𝐶 Concentração do poluente

−→𝑣 Vetor velocidade−→𝐾 Vetor difusão

𝑅 Termo fonte da Equação de Transporte

𝐾𝑥 Coeficiente de difusão em 𝑥

𝐾𝑦 Coeficiente de difusão em 𝑦

𝜅 Coeficiente de decaimento

𝜉, 𝜂 Posicionamento em coordenadas generalizadas

𝐽 Jacobiano da transformação

𝑃 , 𝑄 Termos fonte das equações governantes para a geração de malhas

𝑛𝑖 Número total de linhas na direção 𝜉

𝑛𝑗 Número total de linhas na direção 𝜂

𝑃 Ponto cardeal central

𝐸 Ponto cardeal leste

𝑊 Ponto cardeal oeste

𝑁 Ponto cardeal norte

𝑆 Ponto cardeal sul

𝑁𝑊 Ponto cardeal noroeste

𝑆𝑊 Ponto cardeal sudoeste

𝑆𝐸 Ponto cardeal sudeste

𝑁𝐸 Ponto cardeal nordeste

𝐷 Downstream

𝑈 Upstream

𝑅 Remote-upstream

𝑣𝑓 Velocidade na face 𝑓

𝜏 Tempo em coordenadas generalizadas

𝒟(𝐶) Termo difusivo da equacão de transporte em coordenadas cartesianas

𝒜(𝐶) Termo advectivo da equação de transporte em coordenadas cartesianas

𝑈, 𝑉 Componentes contravariantes para o vetor velocidade

D(𝐶) Termo difusivo da equação de transporte em coordenadas generalizadas

A (𝐶) Termo advectivo da equação de transporte em coordenadas generaliza-das

𝑡𝑓 Tempo final

𝑣𝑖 Velocidade de injeção prescrita

𝐿 Comprimento característico

𝑅𝑒 Número de Reynolds

𝐶𝑖 Concentração de injeção prescrita

SUMÁRIO

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271.1 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 281.2 Organização do Trabalho . . . . . . . . . . . . . . . . . . . . . . . 29

2 FUNDAMENTAÇÃO TEÓRICO-METODOLÓGICA . . . . . 312.1 Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . . . . . 312.1.1 Equação de Navier-Stokes . . . . . . . . . . . . . . . . . . . . . . 332.1.2 Equação de Transporte . . . . . . . . . . . . . . . . . . . . . . . . 342.2 Malhas Computacionais . . . . . . . . . . . . . . . . . . . . . . . . 352.2.1 Sistema de Coordenadas Generalizadas . . . . . . . . . . . . . . 372.3 Condições Auxiliares . . . . . . . . . . . . . . . . . . . . . . . . . . 392.4 Métodos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . 402.4.1 Método de Diferenças Finitas . . . . . . . . . . . . . . . . . . . . 412.4.2 Esquema Upwind . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3 TRABALHOS CORRELATOS . . . . . . . . . . . . . . . . . . . 45

4 PROCEDIMENTOS METODOLÓGICOS . . . . . . . . . . . . 494.1 Desenvolvimento Matemático . . . . . . . . . . . . . . . . . . . . 494.1.1 Transformação de Coordenadas da Equação Governante . . . 504.1.2 Discretização dos Termos da Equação Governante . . . . . . . 564.1.2.1 Discretização em Coordenadas Cartesianas . . . . . . . . . . . . . . . . 57

4.1.2.2 Discretização em Coordenadas Generalizadas . . . . . . . . . . . . . . 58

4.2 Implementação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5 RESULTADOS NUMÉRICOS E DISCUSSÃO . . . . . . . . . 735.1 Duto Simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755.2 Duto com um Obstáculo . . . . . . . . . . . . . . . . . . . . . . . 795.3 Duto com dois Obstáculos . . . . . . . . . . . . . . . . . . . . . . 815.4 Análise dos Tempos de Processamento . . . . . . . . . . . . . . 84

6 CONSIDERAÇÕES FINAIS . . . . . . . . . . . . . . . . . . . . 87

REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

27

1 INTRODUÇÃO

Inseridas no contexto atual, questões ambientais relacionadas à poluição ganharamdestaque nos últimos anos, pelo fato de que o problema tem provocado consequênciasimpactantes nos ecossistemas do planeta, assim como na saúde humana. Dentre os diversostipos de poluição, a poluição atmosférica constitui-se de grande relevância, devido àsseveras influências nas alterações climáticas da Terra e na composição química do ar[3]. Um agravante do problema está na característica dinâmica dos poluentes durante oprocesso de dispersão [3], a qual amplifica o seu raio de alcance.

Algumas das maiores causas para o aumento substancial dos níveis de poluiçãono meio ambiente estão concentradas na crescente demanda mundial por energia e noaumento da produção industrial [3, 4], impulsionada pela Revolução Industrial ocorridano século XVIII. Atualmente, existem diversas fontes de poluição atmosférica, e essasfontes podem ser classificadas de acordo com a sua origem, da seguinte maneira [4]:

• Antropogênicas: emissões originadas por atividades humanas, tais como emissõesindustriais, de veículos automotores e de usinas;

• Geogênicas: emissões originadas pelo próprio planeta, tais como emissões vulcânicas,de águas salgadas e incêndios naturais;

• Biogênicas: emissões originadas por organismos vivos, tais como emissões de com-postos orgânicos voláteis (VOC) em florestas e de metano (𝐶𝐻4) em pântanos.

As fontes geogênicas e biogênicas também são denominadas de fontes naturais. No quediz respeito às fontes antropogênicas, ainda existe a classificação quanto à mobilidade,podendo, dessa forma, serem consideradas fontes móveis ou fixas [3].

De modo formal, um poluente atmosférico pode ser definido como qualquer subs-tância emitida no ar por meio de alguma das fontes mencionadas anteriormente, podendoprovocar efeitos adversos [4]. Os poluentes podem ser categorizados [4] como segue:

• Primários: substâncias emitidas diretamente na atmosfera, a partir das fontes;

• Secundários: substâncias emitidas indiretamente na atmosfera, formadas a partir dereações químicas entre os poluentes primários.

Existe uma ampla variedade de compostos químicos poluentes, tanto primários quantosecundários, nos estados gasosos, líquidos e sólidos. Uma lista desses poluentes pode serconsultada no livro de Daly et al. [4].

28

Com relação aos efeitos adversos provocados pela poluição atmosférica, acredita-se, de acordo com grande parte dos pesquisadores, que o problema se constitui a principalcausa dos fenômenos do efeito estufa e das chuvas ácidas [5], ocorridos no meio ambiente.Nos seres humanos, distúrbios respiratórios, alergia, lesões degenerativas no sistema ner-voso e em órgãos vitais e câncer são exemplos de doenças ou sintomas que podem serprovocados pela poluição do ar [5].

Por essas razões, o desenvolvimento e a aplicação de estratégias de controle paraa preservação e melhoria da qualidade do ar tornou-se imprescindível [6]. Algumas dasprincipais medidas que podem ser tomadas para amenizar o problema são: a redução naquantidade de emissões, o tratamento dos poluentes antes da emissão ou o reposiciona-mento das fontes emissoras [5]. Entretanto, um estudo preliminar faz-se necessário paracada situação, com a finalidade de se encontrar a solução mais eficiente e viável [5].

Nesse sentido, a modelagem da dispersão de poluentes na atmosfera constitui-seuma forma bastante interessante de se avaliar e prever o seu comportamento [7]. Essamodelagem consiste basicamente na descrição matemática do transporte de poluentes noar [1] e é realizada através da simulação numérica de uma equação diferencial: a equação detransporte ou também denominada equação de advecção-difusão. Essa equação possibilitaestimar os níveis de concentração de um determinado poluente conforme o tempo, dadaa região que está sendo investigada.

A área de Dinâmica dos Fluidos Computacional (CFD - Computational Fluid Dy-namics) é responsável por proporcionar a simulação real de escoamentos, em escalas re-duzidas ou, até mesmo, em largas escalas, por meio da solução numérica das equaçõesgovernantes do problema [8]. Em grande parte dos modelos matemáticos, a solução analí-tica das equações diferenciais para problemas de escoamento não pode ser obtida [8]. Poresse motivo, técnicas computacionais possibilitam transformar essas equações em sistemasde equações algébricas [8], tornando possível, no âmbito computacional, a sua solução.

1.1 Objetivos

Este trabalho tem como objetivo geral, além de estudar e analisar as atuais me-todologias empregadas para a modelagem da dispersão de poluentes na atmosfera, apre-sentar uma nova proposta de modelagem numérica, aplicando um conjunto de técnicasmatemáticas e computacionais do campo de CFD para a simulação numérica do problema.

Nessa proposta, a equação de transporte é representada no sistema de coordenadasgeneralizadas bidimensionais e discretizada por meio do método de diferenças finitas, uti-lizando o esquema upwind de primeira ordem para a resolução do seu termo convectivo.Mais detalhes acerca da equação governante e das teorias envolvidas nas técnicas em-pregadas podem ser encontrados na fundamentação teórico-metodológica deste trabalho,

29

presente no Capítulo 2.

É importante levar em consideração que a equação de transporte é dependentede dois outros algoritmos, a geração de malhas no sistema de coordenadas generalizadase a resolução numérica da Equação de Navier-Stokes, as quais foram provenientes dostrabalhos de Cirilo et al. [9] e de Barba [2], respectivamente.

1.2 Organização do Trabalho

Este trabalho encontra-se organizado da seguinte forma:

• No Capítulo 2, esclarecem-se os principais conceitos atrelados ao modelo de disper-são de poluentes proposto;

• No Capítulo 3, realiza-se uma revisão de trabalhos relacionados, apontando as res-pectivas metodologias aplicadas;

• No Capítulo 4, demonstram-se os procedimentos metodológicos empregados para aformulação do modelo;

• No Capítulo 5, comprova-se, por meio de experimentos numéricos, a aplicabilidadedo modelo;

• No Capítulo 6, identificam-se as considerações finais, incluindo as contribuições e aspróximas etapas deste trabalho.

31

2 FUNDAMENTAÇÃO TEÓRICO-METODOLÓGICA

Neste capítulo, apresentam-se os conceitos mais relevantes para facilitar o enten-dimento da metodologia adotada no modelo proposto. Expondo uma visão geral dessesconceitos, na Seção 2.1, define-se o conceito de equações diferenciais e suas classificações,enfatizando as equações aplicadas neste trabalho; na Seção 2.2, define-se a concepção demalhas computacionais e suas classificações, destacando o sistema de coordenadas gene-ralizadas; na Seção 2.3, define-se a noção de condições auxiliares, associadas às soluçõesde equações diferenciais; finalmente, na Seção 2.4, define-se a ideia de métodos numéricos,salientando o método de diferenças finitas e o esquema upwind.

2.1 Equações Diferenciais

No contexto da matemática, equação diferencial é a denominação dada às equa-ções que contêm derivadas [10]. As equações diferenciais são extremamente importantesna área de Ciências Exatas pelo fato de que possibilitam determinar o comportamentode diversos fenômenos físicos e/ou químicos. Em outras palavras, uma grande variedadede aplicações é proporcionada pelo estudo e resolução das equações diferenciais [10, 11].Alguns exemplos de aplicações modeladas por essas equações são: problemas envolvendoo movimento de fluidos, o fluxo de corrente elétrica em circuitos, a dissipação de ca-lor em objetos sólidos, a propagação e detecção de ondas sísmicas e a variação da taxapopulacional [10].

Nas equações diferenciais, uma função constitui-se a incógnita do problema, o quepermite classificá-las da seguinte maneira [10]:

• Equações Diferenciais Ordinárias (EDOs): equações compostas por derivadas sim-ples, isto é, a função desconhecida depende de uma única variável independente;

• Equações Diferenciais Parciais (EDPs): equações compostas por derivadas parciais,ou seja, a função desconhecida depende de duas ou mais variáveis independentes.

Apesar de possuírem soluções que descrevem comportamentos bastante diferentes, umaEDO consiste basicamente em um caso especial de EDP [11].

Em muitas situações, outra questão que também deve ser levada em considera-ção está no fato de que pode existir mais de uma função desconhecida. Nesses casos, énecessário um sistema de equações diferenciais para a sua resolução [10].

32

Um conceito muito importante quando se trata de equações diferenciais é a ordem,a qual corresponde à ordem da maior derivada presente na equação [12, 10]. Esse conceitopossibilita formular, de maneira matemática, as EDOs e as EDPs.

Uma EDO de ordem 𝑛 para uma função desconhecida 𝑦(𝑥) é descrita pela Equação(2.1) [11]:

𝐹 (𝑥, 𝑦(𝑥), 𝑦′(𝑥), . . . , 𝑦(𝑛)(𝑥)) = 0. (2.1)

Uma EDP para uma função desconhecida 𝑢(𝑥, 𝑦), isto é, dependente de duas variáveis, éexpressa, por exemplo, pela Equação (2.2) de segunda ordem [11]:

𝐹 (𝑥, 𝑦, 𝑢, 𝑢𝑥, 𝑢𝑦, 𝑢𝑥𝑥, 𝑢𝑥𝑦, 𝑢𝑦𝑦) = 0. (2.2)

As EDPs de ordem 𝑛 com mais variáveis independentes possuem formulações análogas.Ainda, nas Equações (2.1) e (2.2), 𝐹 é uma função conhecida, a qual representa a equaçãodiferencial, dado o problema.

Outra classificação possível para as equações diferenciais é realizada segundo àlinearidade, ou seja, uma equação diferencial é dita linear se a função 𝐹 é linear emrelação à função desconhecida e suas derivadas [10, 11]. Atualmente, os métodos existentespara a resolução de equações lineares são bastante satisfatórios, porém ainda há grandedificuldade em se resolver equações não lineares [10]. Por esse motivo, quando o fenômenopode ser representado por uma equação linear, é comum realizar um processo denominadode linearização, consistindo na aproximação de uma equação não linear por uma linear[10].

De modo particular, as EDPs possuem também uma classificação específica. Con-sidere reescrever a Equação (2.2) como combinação linear da função desconhecida 𝑢(𝑥, 𝑦)e de suas derivadas, de acordo com a Equação (2.3) [13]:

𝑎𝜕2𝑢

𝜕𝑥2 + 𝑏𝜕2𝑢

𝜕𝑥𝜕𝑦+ 𝑐

𝜕2𝑢

𝜕𝑦2 + 𝑑𝜕𝑢

𝜕𝑥+ 𝑒

𝜕𝑢

𝜕𝑦+ 𝑓𝑢 + 𝑔 = 0, (2.3)

em que 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, 𝑓 e 𝑔 podem ser funções dependentes de 𝑥, 𝑦 e 𝑢. Sabendo-se queΔ = 𝑏2 − 4𝑎𝑐, a EDP, representada na Equação (2.3), pode ser classificada como segue[13]:

• Elíptica: quando Δ < 0;

• Parabólica: quando Δ = 0;

• Hiperbólica: quando Δ > 0.

Tendo em vista que boa parte das equações diferenciais não possuem soluçõesanalíticas, isto é, soluções exatas, a obtenção de uma solução numérica apropriada éfundamental. Para tal fim, é necessário o desenvolvimento de um modelo matemático

33

capaz de representar o problema em questão [2]. Esse modelo é constituído pelas equaçõese pelo domínio físico do problema, os quais devem estar adequadamente expressos [2].

O domínio físico, representado na forma de malhas computacionais, é explicado demaneira mais detalhada na Seção 2.2. As equações diferenciais aproximadas são encon-tradas através de métodos de resolução numérica, os quais são especificados na Seção 2.4.Além disso, os métodos numéricos, bem como as equações diferenciais, estão vinculadosa condições auxiliares, descritas na Seção 2.3.

Apresentada uma visão geral acerca das equações diferenciais, nas Subseções 2.1.1e 2.1.2, descrevem-se duas equações diferenciais indispensáveis para a proposta destetrabalho.

2.1.1 Equação de Navier-Stokes

A partir dos estudos de Leonard Euler, Claude Navier, Simeon Poisson e GeorgeStokes, foram deduzidas as equações que governam o movimento de fluidos, as quaispassaram a ser denominadas como equações de Navier-Stokes (Cenedese[14], 2005 apudBARBA[2], 2015). Classificadas como EDPs não lineares, essas equações podem ser apre-sentadas de diferentes formas, dependendo das características do problema investigado[2].

De maneira genérica, não é possível encontrar soluções analíticas para as equaçõesde Navier-Stokes [15]. Contudo, para versões simplificadas dessas equações, soluções ana-líticas são encontradas. No trabalho de Barba [2], cujo modelo matemático proposto foiutilizado também neste trabalho para se obter os campos de velocidade necessários paraa equação de transporte, as equações de Navier-Stokes são descritas pelas Equações (2.4)e (2.5), em duas dimensões.

𝜕

𝜕𝑡(𝜌𝑢)⏟ ⏞

Termo temporal

+ 𝜕

𝜕𝑥(𝜌𝑢𝑢) + 𝜕

𝜕𝑦(𝜌𝑢𝑣)⏟ ⏞

Termo convectivo

= −𝜕𝑝

𝜕𝑥⏟ ⏞ Termo de pressão

+ 𝜇

(𝜕2𝑢

𝜕𝑥2 + 𝜕2𝑢

𝜕𝑦2

)⏟ ⏞

Termo difusivo

+𝑅1, (2.4)

𝜕

𝜕𝑡(𝜌𝑣)⏟ ⏞

Termo temporal

+ 𝜕

𝜕𝑥(𝜌𝑢𝑣) + 𝜕

𝜕𝑦(𝜌𝑣𝑣)⏟ ⏞

Termo convectivo

= −𝜕𝑝

𝜕𝑦⏟ ⏞ Termo de pressão

+ 𝜇

(𝜕2𝑣

𝜕𝑥2 + 𝜕2𝑣

𝜕𝑦2

)⏟ ⏞

Termo difusivo

+𝑅2, (2.5)

nas quais 𝑢 e 𝑣 representam as componentes da velocidade, respectivamente, em 𝑥 e 𝑦,e 𝑝 corresponde à pressão em 𝑥 e 𝑦. Nas derivadas parciais, 𝑡 refere-se ao tempo, e 𝑥, 𝑦,ao posicionamento no plano cartesiano. Nos termos temporais e convectivos, 𝜌 indica amassa específica. Nos termos difusivos, 𝜇 descreve a viscosidade dinâmica. 𝑅1 e 𝑅2 são ostermos fonte1, respectivamente, nas direções 𝑥 e 𝑦.

1 Os termos fonte comumente informam a perda ou o ganho de uma determinada quantidade (Costa eRibeiro[16], 2010 apud BARBA[2], 2015).

34

Dessa forma, as Equações de Navier-Stokes fornecem os módulos dos vetores develocidade, assim como o valor da pressão em cada célula da malha computacional.

2.1.2 Equação de Transporte

Como já foi mencionado na introdução deste trabalho, a equação diferencial res-ponsável pela modelagem da dispersão de poluentes na atmosfera é denominada equaçãode transporte, ou equação de advecção-difusão. Assim como a Equação de Navier-Stokes,essa equação possui várias formulações diferentes, dependendo das características da situ-ação estudada. Quando a equação possui termo fonte não nulo, pode-se ter a denominaçãode equação de advecção-difusão-reação, uma vez que o termo fonte, nesse caso, representao termo reativo da equação diferencial.

Nessa equação, dois fenômenos físicos são extremamente relevantes: a advecção ea difusão. O processo de advecção, correspondente ao processo de convecção da Equaçãode Navier-Stokes, consiste naquele em que a nuvem de poluente é transportada pelodeslocamento do ar sem alterar sua forma. Já o processo de difusão constitui-se naqueleem que a nuvem se espalha à medida que sua concentração diminui [17]. O processo dereação, indicando a perda ou produção de um determinado poluente, pode ser consideradoou não na formulação.

De maneira aplicada, a Figura 1 ilustra os três mecanismos presentes em um pro-cesso de dispersão de poluentes. A imagem representa as emissões da Inco Superstack, umachaminé localizada em Sudbury, Canadá. A deposição apresentada, no caso, é referenteao processo de reação.

Figura 1 – Mecanismos em um processo de dispersão de poluentes - retirada de Stockie[1]

Calculando os níveis de concentração do contaminante em cada instante de tempo,

35

a equação de transporte genérica para uma, duas ou três dimensões pode ser expressa pelaEquação (2.6) [18, 19, 1].

𝜕𝐶

𝜕𝑡+ −→

∇ · (𝐶−→𝑣 ) −−→∇ · (−→𝐾−→

∇ · 𝐶) = 𝑅, (2.6)

em que −→∇ é o vetor gradiente. Já os demais termos são explicados a seguir, na formulação

utilizada nesta proposta de modelagem.

A representação matemática da equação de transporte empregada neste trabalhoé descrita pela Equação (2.7), em duas dimensões.

𝜕𝐶

𝜕𝑡⏟ ⏞ Termo temporal

+ 𝜕

𝜕𝑥(𝐶𝑣𝑥) + 𝜕

𝜕𝑦(𝐶𝑣𝑦)⏟ ⏞

Termo convectivo

− 𝜕

𝜕𝑥

(𝐾𝑥

𝜕𝐶

𝜕𝑥

)− 𝜕

𝜕𝑦

(𝐾𝑦

𝜕𝐶

𝜕𝑦

)⏟ ⏞

Termo difusivo

= −𝜅𝐶⏟ ⏞ Termo reativo

, (2.7)

em que 𝐶 indica a concentração do poluente. Nas derivadas parciais, 𝑡 refere-se ao tempo,e 𝑥, 𝑦, ao posicionamento no plano cartesiano. No termo convectivo, também chamadocomo termo advectivo, 𝑣𝑥 e 𝑣𝑦 correspondem às velocidades dos fluxos de vento da regiãoanalisada, obtidas por meio da solução da Equação de Navier-Stokes e representando,respectivamente, os valores de 𝑢 e 𝑣 das Equações (2.4) e (2.5). No termo difusivo, asconstantes 𝐾𝑥 e 𝐾𝑦 descrevem os coeficientes de difusão, os quais são dependentes domeio em que o poluente se propaga. No termo reativo 𝑅 = −𝜅𝐶, a constante 𝜅 equivaleao coeficiente de decaimento.

Ainda na Equação (2.7), vale ressaltar que, quando 𝑣𝑥 e 𝑣𝑦 são nulos, elimina-se oprocesso de advecção; assim como, quando 𝐾𝑥 e 𝐾𝑦 são nulos, elimina-se o processo dedifusão do poluente [17].

Definido o conceito de equações diferenciais e apresentadas as duas principais EDPsaplicadas neste trabalho, detalha-se o conceito de malhas computacionais, bem como osistema de coordenadas generalizadas na Seção 2.2.

2.2 Malhas Computacionais

Uma das características mais importantes dos métodos numéricos é a presença deuma malha para descrever a geometria do meio que está sendo investigado. Essa malhaconsiste no domínio físico discretizado, ou seja, no domínio físico representado atravésde um conjunto de pontos [15, 20]. O processo de discretização constitui-se um passoindispensável, pelo fato de que esses métodos de resolução não permitem a obtenção deuma solução para todos os pontos do domínio [2].

Uma malha computacional é formada por um agrupamento de células, limitadaspor arestas denominadas de faces e contendo vértices chamados de nós (Marques[21],2005 apud BARBA[2], 2015). Dependendo do tipo de malha, essas células podem pos-

36

suir diferentes formas, tais como triângulos ou quadriláteros quando se consideram duasdimensões, e tetraedros ou hexaedros, em três dimensões [20].

As malhas possuem certas classificações dependendo do parâmetro empregado,como segue [2]:

• De acordo com o espaçamento entre os pontos: malhas uniformes ou não uniformes;

• De acordo com a distribuição espacial dos pontos: malhas estruturadas ou não es-truturadas.

A Figura 2 demonstra, por meio de um exemplo, as principais diferenças entre as malhasestruturadas, Figura 2a, e as não estruturadas, Figura 2b.

(a) Malha estruturada (b) Malha não estruturada

Figura 2 – Classificação de malhas quanto à distribuição espacial dos pontos.

Existe também o conceito de malhas híbridas ou estruturadas em blocos, as quaisconsistem em malhas constituídas por uma pequena quantidade de elementos estrutura-dos combinados através de um padrão não estruturado [20]. Contudo, sua geração paraquaisquer geometrias se constitui ainda um obstáculo [20].

As malhas não estruturadas possuem a característica de serem mais versáteis que asestruturadas [22, 20], o que permite a sua adaptabilidade a geometrias consideradas com-plexas. Em contrapartida, apresentam dificuldade de ordenação dos volumes elementares,ou seja, as células da malha, implicando diretamente no consumo adicional de memória,devido ao aumento das dimensões da matriz computacional. Em virtude da facilidade demanipulação dos elementos da malha, considera-se, neste trabalho, a utilização de malhasestruturadas.

No que diz respeito ao sistema de coordenadas, o sistema de coordenadas cartesi-anas é comumente adotado para a representação do domínio físico do problema, devidoa sua simplicidade. Entretanto, para problemas que envolvem geometrias complexas, essesistema de coordenadas acarreta uma má adequação da fronteira do problema, pelo fatode que o domínio físico não coincide com o domínio da malha computacional.

A partir dessas motivações, surge então a ideia do sistema de coordenadas gene-ralizadas, utilizado no modelo sugerido e melhor esclarecido na Subseção 2.2.1.

37

2.2.1 Sistema de Coordenadas Generalizadas

O sistema de coordenadas generalizadas, também denominadas como coordenadascurvilíneas, cujas malhas produzidas são do tipo estruturadas, possibilitam também aadaptabilidade presente em malhas não estruturadas, permitindo, dessa maneira, a repre-sentação computacional se adequar à geometria do problema, a qual pode ser complexaou não [22, 23, 24, 25, 9]. Outro fator que se deve levar em consideração é a possibilidadede se concentrar pontos da malha onde for necessário, proporcionando um estudo maisaprofundado em uma certa localização dentro da geometria.

Para que se possa descrever um objeto em estudo no sistema de coordenadasgeneralizadas, com coordenadas dadas pelas variáveis (𝜉, 𝜂), é indispensável realizar umatransformação, a qual se fundamenta no sistema de coordenadas cartesianas (𝑥, 𝑦). Coma finalidade de ilustrar as diferenças entre os dois sistemas de coordenadas, apresentam-seduas malhas na Figura 3.

(a) Coordenadas cartesianas (b) Coordenadas generalizadas

Figura 3 – Sistemas de coordenadas

As métricas da transformação são responsáveis por efetuar essa mudança do domí-nio físico (𝑥, 𝑦) ao domínio computacional, ou transformado, (𝜉, 𝜂) e realizar compensaçõesnecessárias [22], uma vez que o plano computacional adotado é do tipo retangular, e seusvolumes elementares, com dimensões unitárias, isto é, os espaçamentos das células Δ𝜉 = 1e Δ𝜂 = 1. Em outras palavras, ainda que as células possuam espaçamentos arbitrários nodomínio físico, suas dimensões são unitárias no domínio computacional, possibilitando asimplificação de certas formulações no decorrer da transformação de coordenadas.

As métricas de transformação são descritas matematicamente por:

𝜕𝜉

𝜕𝑥= 𝐽

𝜕𝑦

𝜕𝜂,

𝜕𝜉

𝜕𝑦= −𝐽

𝜕𝑥

𝜕𝜂,

𝜕𝜂

𝜕𝑥= −𝐽

𝜕𝑦

𝜕𝜉,

𝜕𝜂

𝜕𝑦= 𝐽

𝜕𝑥

𝜕𝜉, (2.8)

onde 𝐽 =(

𝜕𝑥𝜕𝜉

𝜕𝑦𝜕𝜂

− 𝜕𝑥𝜕𝜂

𝜕𝑦𝜕𝜉

)−1é o Jacobiano da transformação [22, 25]. O Jacobiano é uma

variável bastante importante pelo fato de que, quando a malha é mal construída, este tendea ser nulo, não sendo possível a representação do plano físico [22, 24]. O problema ocorre,

38

uma vez que a área do elemento no plano físico corresponde ao inverso do Jacobiano, ouseja, 𝑆 = 1

𝐽, quando Δ𝜉 = Δ𝜂 = 1 é assumido.

Para se realizar a geração de malhas computacionais, existem algumas metodo-logias na literatura, dentre as quais se emprega, neste trabalho, a abordagem utilizandoequações diferenciais elípticas. Esse tipo de equações possibilita obter soluções que nãogeram Jacobiano nulo no domínio, e ainda soluções nas quais duas linhas coordenadas, 𝜉

ou 𝜂, nunca se interceptam [22], constituindo, dessa forma, as principais motivações paratal uso.

Vale salientar que o sistema de coordenadas generalizadas parte do pressuposto deque o número de linhas que partem da fronteira esquerda da geometria deve ser igual aonúmero de linhas que chegam na fronteira direita, e, de maneira análoga, o mesmo ocorreentre as fronteiras superior e inferior.

Assim sendo, as equações governantes para a geração de malhas bidimensionaissão expressas por:

−→∇2𝜉 = 𝑃 (𝜉, 𝜂), (2.9)−→∇2𝜂 = 𝑄(𝜉, 𝜂), (2.10)

em que −→∇ é o vetor gradiente. 𝑃 e 𝑄 são os termos fonte, os quais, nesse caso, possibilitam

determinar a concentração de pontos da malha conforme a necessidade [22, 9]. Essasfunções são dadas respectivamente pelas Equações (2.11) e (2.12) [22, 9, 2].

𝑃 (𝜉, 𝜂) = −𝑛𝑗∑

𝑗=1𝑎𝑗sign(𝜉 − 𝜉𝑗)𝑒−𝑐𝑗 |𝜉−𝜉𝑗 | −

𝑛𝑖∑𝑖=1

𝑏𝑖sign(𝜉 − 𝜉𝑖)𝑒−𝑑𝑖

√(𝜉−𝜉𝑖)2+(𝜂−𝜂𝑖)2

, (2.11)

𝑄(𝜉, 𝜂) = −𝑛𝑗∑

𝑗=1𝑎𝑗sign(𝜂 − 𝜂𝑗)𝑒−𝑐𝑗 |𝜂−𝜂𝑗 | −

𝑛𝑖∑𝑖=1

𝑏𝑖sign(𝜂 − 𝜂𝑖)𝑒−𝑑𝑖

√(𝜉−𝜉𝑖)2+(𝜂−𝜂𝑖)2

, (2.12)

em que 𝑛𝑖 e 𝑛𝑗, existentes nos somatórios, indicam o número total de linhas nas direções 𝜉

e 𝜂 respectivamente. Os termos 𝑎𝑗, 𝑏𝑖, 𝑐𝑗 e 𝑑𝑖 correspondem a números reais ajustados pormeio de experimentação numérica, com a finalidade de atrair as linhas 𝜉 e 𝜂 em direçãoa 𝜉𝑖 e 𝜂𝑖 [9].

Aplicando as métricas de transformação, mostradas na Equação (2.8), as equa-ções governantes devem ser adaptadas para o sistema de coordenadas generalizadas. Apartir das Equações (2.9) e (2.10), essas equações, após a transformação, são descritas,respectivamente, pelas Equações (2.13) e (2.14).

𝛼𝜕2𝑥

𝜕𝜉2 + 𝛾𝜕2𝑥

𝜕𝜂2 − 2𝛽𝜕2𝑥

𝜕𝜉𝜂+ 1

𝐽2

(𝑃

𝜕𝑥

𝜕𝜉+ 𝑄

𝜕𝑥

𝜕𝜂

)= 0, (2.13)

𝛼𝜕2𝑦

𝜕𝜉2 + 𝛾𝜕2𝑦

𝜕𝜂2 − 2𝛽𝜕2𝑦

𝜕𝜉𝜂+ 1

𝐽2

(𝑃

𝜕𝑦

𝜕𝜉+ 𝑄

𝜕𝑦

𝜕𝜂

)= 0, (2.14)

39

em que 𝛼 =(

𝜕𝑥𝜕𝜂

)2+(

𝜕𝑦𝜕𝜂

)2, 𝛾 =

(𝜕𝑥𝜕𝜉

𝜕𝑦𝜕𝜉

)2e 𝛽 = 𝜕𝑥

𝜕𝜉𝜕𝑥𝜕𝜂

+ 𝜕𝑦𝜕𝜉

𝜕𝑦𝜕𝜂

.

Partindo das Equações (2.13) e (2.14), deve-se aplicar então algum método nu-mérico de resolução. Em um primeiro momento, para a aproximação numérica dessasequações, generalizam-se as derivadas parciais em 𝑥 e 𝑦, utilizando uma única variável 𝜑

para a sua representação [22, 9], como demonstra a Equação (2.15).

𝛼𝜕2𝜑

𝜕𝜉2 + 𝛾𝜕2𝜑

𝜕𝜂2 − 2𝛽𝜕2𝜑

𝜕𝜉𝜂+ 1

𝐽2

(𝑃

𝜕𝜑

𝜕𝜉+ 𝑄

𝜕𝜑

𝜕𝜂

)= 0. (2.15)

Visando melhorar também a representação da solução numérica, utiliza-se uma rotulagemde índices nos nós da malha, ilustrada pela Figura 4, através dos pontos cardeais 𝑃

(centro), 𝐸 (leste), 𝑊 (oeste), 𝑁 (norte), 𝑆 (sul), 𝑁𝑊 (noroeste), 𝑆𝑊 (sudoeste), 𝑆𝐸

(sudeste) e 𝑁𝐸 (nordeste). No âmbito computacional, esses pontos cardeais podem serinterpretados como sendo posições em uma matriz bidimensional.

Figura 4 – Rotulagem de índices

Portanto, a solução numérica da equação genérica para a geração da malha, Equa-ção (2.15), aproximando-se os termos das derivadas para o ponto 𝑃 pelo método de dife-renças finitas centrais, é indicada pela Equação (2.16). O método numérico empregado éexplicado detalhadamente na Subseção 2.4.1.

𝜑𝑃 = 1𝐴𝑃

(𝐴𝐸𝜑𝐸 + 𝐴𝑊 𝜑𝑊 + 𝐴𝑁𝜑𝑁 + 𝐴𝑆𝜑𝑆+𝐴𝑁𝐸𝜑𝑁𝐸 + 𝐴𝑆𝐸𝜑𝑆𝐸

+𝐴𝑁𝑊 𝜑𝑁𝑊 + 𝐴𝑆𝑊 𝜑𝑆𝑊 ), (2.16)

em que 𝐴𝑃 = 2𝛼 + 2𝛾, 𝐴𝑁 = 𝛾 + 𝑄2𝐽2 , 𝐴𝑆𝐸 = 𝛽

2 , 𝐴𝐸 = 𝛼 + 𝑃2𝐽2 , 𝐴𝑆 = 𝛾 − 𝑄

2𝐽2 , 𝐴𝑁𝑊 = 𝛽2 ,

𝐴𝑊 = 𝛼 − 𝑃2𝐽2 , 𝐴𝑁𝐸 = −𝛽

2 e 𝐴𝑆𝑊 = −𝛽2 , com as derivadas presentes em 𝐽 também

aproximadas por diferenças centrais.

Apontadas as principais teorias associadas às malhas computacionais e ao sistemade coordenadas generalizadas, apresenta-se, na Seção 2.3, o conceito de condições auxili-ares para a solução numérica de equações diferenciais.

2.3 Condições Auxiliares

Equações diferenciais podem possuir soluções gerais e particulares. De maneiraexemplificada, a solução geral 𝑢(𝑥, 𝑦) de uma EDP, como apresentada na Equação (2.2),

40

sobre um conjunto tal contido em R2 constitui-se a solução que engloba todas as soluçõespossíveis sobre esse conjunto [26]. Já uma solução particular da mesma EDP consiste emuma função específica que a satisfaça, sob condições auxiliares [26]. Definições análogassão válidas para EDOs e outras EDPs.

Para se resolver numericamente equações diferenciais, a determinação adequadadas condições auxiliares é fundamental, uma vez que, quando especificadas de formaincorreta, geralmente produzem soluções não únicas ou fisicamente incompatíveis com oproblema [15]. Existem duas classificações de condições auxiliares: as condições iniciaise de fronteira. As condições iniciais descrevem o estado inicial do problema, informandoo valor inicial da solução [2]. Já as condições de fronteira, também denominadas comocondições de contorno, demonstram o comportamento da solução nas fronteiras da regiãoanalisada [2].

Ainda, as condições de fronteira podem ser subdivididas em três tipos básicos [15],sobre a variável dependente 𝜑 na fronteira da região 𝑅:

• Dirichlet: o valor de 𝜑 é especificado em 𝛿𝑅;

• Neumann: o gradiente normal 𝜕𝜑𝜕𝑛

= 𝑓 ou tangencial 𝜕𝜑𝜕𝑠

= 𝑔 é especificado em 𝛿𝑅;

• Robin: uma combinação linear das condições de fronteira tipo Dirichlet e tipo Neu-mann é especificada em 𝛿𝑅.

Vale ressaltar que um mesmo problema pode apresentar diferentes condições de fronteira,significando que ele possui condições de fronteira mistas [15], como é o caso deste trabalho,no qual se empregam as condições de fronteira tipo Dirichlet e tipo Neumann.

Definido o conceito de condições auxiliares e sua importância, apresentam-se osmétodos numéricos de resolução de equações diferenciais na Seção 2.4, destacando-se ométodo de diferenças finitas aplicado nesta proposta.

2.4 Métodos Numéricos

Para se solucionar numericamente uma equação diferencial, assim como o do-mínio físico, também é necessário discretizá-la, realizando manipulações algébricas paraexpressá-la na forma de operações aritméticas facilmente interpretadas pelo computador[2, 13]. Essa discretização é realizada através de alguma metodologia, dentre as quais asmais conhecidas atualmente são [27]: diferenças finitas, elementos finitos, volumes finitose os métodos espectrais.

Na proposta de modelagem deste trabalho, aplica-se o método de diferenças fi-nitas como estratégia de discretização das Equações de Navier-Stokes e de transporte,

41

explicado com mais detalhes na Subseção 2.4.1. Devido a problemas numéricos decorren-tes da discretização do termo convectivo, emprega-se o esquema upwind como técnica deaproximação, apresentado na Subseção 2.4.2.

2.4.1 Método de Diferenças Finitas

O método numérico por diferenças finitas realiza o processo de discretização dasequações diferenciais em função dos pontos (nós) da malha computacional, nos quaisse deve representar a sua solução exata [13]. Essa aproximação pode ser encontrada devárias maneiras, das quais as mais empregadas são [15]: a interpolação polinomial e aexpansão em série de Taylor. Neste trabalho, adota-se a expansão em série de Taylor paraa aproximação.

Dessa maneira, para desenvolver as equações de aproximação, considere uma fun-ção 𝑢 = 𝑢(𝑥) contínua e com derivadas contínuas, além de um espaçamento ℎ. A expansãode 𝑢 em série de Taylor é descrita por [13]:

𝑢(𝑥 + ℎ) = 𝑢(𝑥) + ℎ𝑢′(𝑥) + ℎ2

2 𝑢′′(𝑥) + ℎ3

6 𝑢′′′(𝑥) + . . . (2.17)

Da mesma forma,

𝑢(𝑥 − ℎ) = 𝑢(𝑥) − ℎ𝑢′(𝑥) + ℎ2

2 𝑢′′(𝑥) − ℎ3

6 𝑢′′′(𝑥) + . . . (2.18)

A partir da Equação (2.17), desprezando-se a expressão ℎ2

2 𝑢′′(𝑥) + ℎ3

6 𝑢′′′(𝑥) + . . . ,tem-se que 𝑢(𝑥+ℎ) ≈ 𝑢(𝑥)+ℎ𝑢′(𝑥). Dessa forma, a aproximação progressiva (ou avançada)para a primeira derivada é denotada por [13]:

𝑢′(𝑥) = 𝑢(𝑥 + ℎ) − 𝑢(𝑥)ℎ

+ 𝑂(ℎ), (2.19)

em que 𝑂(ℎ) = ℎ2 𝑢′′(𝑥) + ℎ2

6 𝑢′′′(𝑥) + . . . é o erro de truncamento da série (ET).

Analogamente, realizando o mesmo processo para a Equação (2.18), a aproximaçãoregressiva (ou atrasada) para a primeira derivada é expressa por [13]:

𝑢′(𝑥) = 𝑢(𝑥) − 𝑢(𝑥 − ℎ)ℎ

+ 𝑂(ℎ), (2.20)

em que ET é 𝑂(ℎ) = ℎ2 𝑢′′(𝑥) − ℎ2

6 𝑢′′′(𝑥) + . . . .

Além disso, subtraindo a Equação (2.17) da Equação (2.18) e realizando algumasmanipulações matemáticas, obtém-se a aproximação de diferença central para a primeiraderivada, dada por [13]:

𝑢′(𝑥) = 𝑢(𝑥 + ℎ) − 𝑢(𝑥 − ℎ)2ℎ

+ 𝑂(ℎ2), (2.21)

em que 𝑂(ℎ2) = ℎ2

6 𝑢′′′(𝑥) + ℎ4

120𝑢𝑣(𝑥) + . . . .

42

De maneira análoga, é possível calcular aproximações para derivadas de ordenssuperiores. Adicionando as Equações (2.17) e (2.18), obtém-se, por meio de alguns desen-volvimentos matemáticos, a aproximação de diferença central para a segunda derivada,expressada por [13]:

𝑢′′(𝑥) = 𝑢(𝑥 + ℎ) − 2𝑢(𝑥) + 𝑢(𝑥 − ℎ)ℎ2 + 𝑂(ℎ2), (2.22)

em que a notação 𝑂(ℎ2) = ℎ2

12 𝑢𝑖𝑣(𝑥)+. . . indica que o erro é de segunda ordem de acurácia.

Vale ressaltar ainda que as aproximações demonstradas pelas Equações (2.19 -2.22) podem ser utilizadas para funções de várias variáveis [13], alterando apenas suanotação. Essa notação faz uso da rotulagem de índices, ilustrada pela Figura 4, e é em-pregada de forma idêntica nas discretizações realizadas neste trabalho. Considere umafunção 𝑢 = 𝑢(𝑥, 𝑡), um espaçamento ℎ na direção 𝑥 e um espaçamento 𝑘 na direção 𝑡, aaproximação de diferença central para 𝜕2𝑢

𝜕𝑥2 no ponto 𝑃 , por exemplo, é descrita por [13]:

𝜕2𝑢

𝜕𝑥2

𝑃

= 𝑢|𝐸 − 2𝑢|𝑃 + 𝑢|𝑊ℎ2 + 𝑂(ℎ2). (2.23)

Do mesmo modo, a aproximação de diferença central para 𝜕2𝑢𝜕𝑡2 no ponto 𝑃 é denotada por

[13]:𝜕2𝑢

𝜕𝑡2

𝑃

= 𝑢|𝑁 − 2𝑢|𝑃 + 𝑢|𝑆𝑘2 + 𝑂(𝑘2). (2.24)

Assim sendo, as Equações (2.19 - 2.22) e suas formulações derivadas englobamgrande parte das formas de discretização que são utilizadas pelo método de diferençasfinitas. Para que o computador possa interpretá-las, a implementação das expressões écomumente realizada pela manipulação de matrizes ou vetores. A fim de facilitar a com-preensão, considere, por exemplo, um vetor contendo os pontos da malha pertencentesa um meio unidimensional e a função 𝑢 = 𝑢(𝑥) previamente definida. Dessa forma, osvalores de 𝑢(𝑥 − ℎ), 𝑢(𝑥) e 𝑢(𝑥 + ℎ) presentes nas aproximações correspondem respecti-vamente a 𝑢𝑖−1, 𝑢𝑖 e 𝑢𝑖+1, onde 𝑖 equivale à posição no vetor sobre a qual a derivada de𝑢 está sendo aproximada, como se pode verificar na Figura 5.

Figura 5 – Discretização por Diferenças Finitas

O termo convectivo presente em equações que transportam propriedades físicas,como é o caso das Equações (2.4 - 2.7), provoca, na maioria das vezes, dificuldades de seobter soluções numéricas de boa qualidade (Queiroz e Ferreira[28], 2007 apud BARBA[2],

43

2015). Esse problema depende significativamente do esquema de discretização escolhido,devido à aproximação imprecisa das derivadas advectivas nesse tipo de equações [13].

Grandes esforços estão direcionados à criação de esquemas capazes de resolvero problema, dentre os quais, a estratégia upwind, descrita na Subseção 2.4.2, tem sedestacado [13].

2.4.2 Esquema Upwind

De acordo com o esquema upwind, o termo convectivo é discretizado conforme osinal do vetor velocidade no local, ou seja, no mesmo sentido do fluxo. A técnica se baseiana disposição de três nós computacionais adjacentes ao ponto de discretização [13]: ajusante 𝐷 (Downstream) e as montantes 𝑈 (Upstream) e 𝑅 (Remote-upstream). Pode-seperceber pela Figura 6 que essa disposição está relacionada com o sentido da velocidade𝑣𝑓 , a qual se refere a uma variável convectada 𝜑𝑓 na face 𝑓 de uma célula da malhacomputacional.

Figura 6 – Disposição dos nós 𝐷, 𝑈 e 𝑅

Dessa forma, para avaliar 𝜑𝑓 , basta escrever a variável em função dos três nóscomputacionais [13], por meio da relação:

𝜑𝑓 = 𝜑𝑓 (𝜑𝐷, 𝜑𝑈 , 𝜑𝑅). (2.25)

No que diz respeito às suas classificações, a estratégia upwind pode ser categori-zada segundo a ordem do esquema, do seguinte modo (Queiroz e Ferreira[28], 2007 apudBARBA[2], 2015):

• Primeira ordem: é estável incondicionalmente e produz um caráter difusivo, o qualgeralmente suaviza a solução. Exemplo: esquema FOU (First Order Upwinding);

• Alta ordem: aumenta a acurácia da aproximação numérica, porém pode introduziroscilações não físicas que são capazes de influenciar na convergência da solução.Exemplos: esquemas SOU (Second Order Upwinding), QUICK (Quadratic UpstreamInterpolation for Convective Kinematics) e CUBISTA (Convergent and UniversallyBounded Interpolation Scheme for the Treatment of Advection).

44

Ainda, outros esquemas também são mencionados no trabalho de Corrêa et al. [13]. Nestetrabalho, emprega-se o esquema FOU para a discretização do termo convectivo da equaçãode transporte.

Fundamentado em Corrêa et al. [13], para entender o processo de discretização doesquema FOU, considere um exemplo de aproximação de um dos termos convectivos daEquação (2.7), utilizando diferenças centrais para a primeira derivada, como segue:

𝜕

𝜕𝑥(𝐶𝑣𝑥)

(𝑖,𝑗)

≈(𝐶𝑣𝑥)|(𝑖+ 1

2 ,𝑗) − (𝐶𝑣𝑥)|(𝑖− 12 ,𝑗)

2(

Δ𝑥2

)=

𝐶|(𝑖+ 12 ,𝑗)𝑣𝑥|(𝑖+ 1

2 ,𝑗) − 𝐶|(𝑖− 12 ,𝑗)𝑣𝑥|(𝑖− 1

2 ,𝑗)Δ𝑥

.

(2.26)

Com a finalidade de manter consistência e estabilidade nas propriedades físicas da equaçãode transporte, observe, na Equação (2.26), que a aproximação realizada sobre a posição(𝑖, 𝑗) utilizou os pontos

(𝑖 + 1

2 , 𝑗)

e(𝑖 − 1

2 , 𝑗). Note também que o espaçamento ℎ = Δ𝑥

foi dividido por 2.

As concentrações 𝐶 nos pontos indicados (nó computacional 𝑈) são então cal-culadas dependendo do sinal da velocidade 𝑣𝑥 correspondente e, por consequência, dalocalização dos nós 𝐷 e 𝑅. De forma equacionada, o esquema FOU para 𝐶|(𝑖+ 1

2 ,𝑗) é ex-pressado por:

𝐶|(𝑖+ 12 ,𝑗) =

⎛⎝1 + 𝛼|(𝑖+ 12 ,𝑗)

2

⎞⎠𝐶|(𝑖,𝑗) +⎛⎝1 − 𝛼|(𝑖+ 1

2 ,𝑗)2

⎞⎠𝐶|(𝑖+1,𝑗), (2.27)

em que 𝛼|(𝑖+ 12 ,𝑗) =

⎧⎪⎨⎪⎩1, se 𝑣𝑥|(𝑖+ 1

2 ,𝑗) ≥ 0−1, se 𝑣𝑥|(𝑖+ 1

2 ,𝑗) < 0.

Analogamente, o esquema FOU para 𝐶|(𝑖− 12 ,𝑗) é descrito por:

𝐶|(𝑖− 12 ,𝑗) =

⎛⎝1 + 𝛼|(𝑖− 12 ,𝑗)

2

⎞⎠𝐶|(𝑖−1,𝑗) +⎛⎝1 − 𝛼|(𝑖− 1

2 ,𝑗)2

⎞⎠𝐶|(𝑖,𝑗), (2.28)

em que 𝛼|(𝑖− 12 ,𝑗) =

⎧⎪⎨⎪⎩1, se 𝑣𝑥|(𝑖− 1

2 ,𝑗) ≥ 0−1, se 𝑣𝑥|(𝑖− 1

2 ,𝑗) < 0.

Dessa maneira, tendo em vista que a concentração 𝐶 é equivalente à variávelconvectada 𝜑𝑓 na Equação (2.25), o valor da concentração em 𝑈 (𝜑𝑈) é dependente apenasdo valor da concentração em 𝑅 (𝜑𝑅), uma vez que o valor de 𝛼, influenciado pelo sentidodo vetor velocidade 𝑣𝑥, elimina o valor da concentração em 𝐷 (𝜑𝐷) durante o processo demultiplicação.

Apresentada a fundamentação teórico-metodológica deste trabalho, no Capítulo3, é realizada uma revisão bibliográfica dos modelos para o tratamento da dispersãode poluentes encontrados na literatura, acompanhados com as respectivas metodologiasutilizadas.

45

3 TRABALHOS CORRELATOS

Neste capítulo, um levantamento de material bibliográfico foi realizado em bi-bliotecas digitais como ScienceDirect1, ACM 2 e IEEE Xplore3, utilizando palavras-chaverelevantes e associadas ao assunto. No entanto, obtiveram-se resultados pertinentes apenasna biblioteca ScienceDirect. Para cada referência levantada, buscaram-se os procedimen-tos metodológicos aplicados para a modelagem do problema, com a finalidade de analisarcriticamente as vantagens e desvantagens de cada estratégia adotada.

Em seu artigo, Fang et al. [29] apresenta um novo modelo de ordem reduzida paraa simulação da poluição no ar, no qual se utiliza uma malha adaptativa não estruturadae a técnica de discretização de elementos finitos. Como demonstração do funcionamentodo modelo, realiza-se a aplicação sobre cânions urbanos em ambientes bidimensionais etridimensionais. Comparado a um modelo de elementos finitos completo, os resultadosmostraram que a precisão é mantida, com ricos detalhes dos escoamentos de ar, e osrequisitos computacionais são reduzidos.

Ferragut et al. [30] resolve numericamente a equação diferencial de transporte,empregando também o método de elementos finitos de maneira adaptativa, contudo comcaracterísticas apenas na direção horizontal; na direção vertical, a técnica de diferençasfinitas é empregada. Os poluentes simulados são passivos, isto é, não reativos entre si oucom o ar atmosférico, e avaliados sobre a escala urbana. Ainda, um modelo de campode velocidade do fluxo de vento desenvolvido pelos próprios autores é incluído. Foramaplicadas técnicas de paralelização do algoritmo, melhorando a precisão da solução esimultaneamente reduzindo o tempo computacional.

Analisando a dispersão de poluentes em cânions urbanos, Madalozzo et al. [31]investiga numericamente três casos, os quais se diferenciam apenas nas configurações dageometria. O modelo foi baseado na formulação de elementos finitos utilizada por algunsdos autores, com a finalidade de solucionar problemas de transporte de massa e calor emmicroescala urbana. Nos resultados, foi constatado que o padrão de escoamento do argoverna a dispersão e que houve grande influência da temperatura, quando comparadoscom resultados de simulações isotérmicas, nas quais não há variação na temperatura.

Fazendo uso do modelo Comprehensive Turbulent Aerosol Dynamics and Gas Che-mistry (CTAG), pertencente à área de CFD, Steffens et al. [32] calculou os gradientesespaciais das concentrações do poluente hexafluoreto de enxofre (𝑆𝐹6) atrás de uma bar-reira sólida, nas proximidades de uma estrada. Duas técnicas de CFD foram empregadas,

1 <http://www.sciencedirect.com/>2 <http://dl.acm.org/>3 <http://ieeexplore.ieee.org/Xplore/home.jsp>

46

Reynolds-averaged Navier-Stokes (RANS) e Large-Eddy Simulation (LES), ambas utili-zadas para a simulação de escoamentos turbulentos. Apesar da menor complexidade, atécnica RANS foi considerada insuficiente para se prever a dispersão de poluentes na re-gião analisada. Já a abordagem LES apresentou-se como melhor opção para a simulaçãonessas condições, onde uma zona de recirculação formada atrás da barreira é considerá-vel. Finalmente, o trabalho concluiu que pode ser bastante útil como forma de auxílio àprojeção inteligente de estradas.

Empregando também a técnica LES como metodologia, Zhong et al. [33] estudaa dispersão de poluentes reativos em um profundo cânion urbano, sob condições meteo-rológicas neutras. Por causa da formação de dois vórtices instáveis dentro do cânion, avariação espacial dos poluentes é bastante significante. Como conclusão, comprovou-se queo modelo desenvolvido poderia auxiliar na criação de estratégias de planejamento urbanoe administração de tráfego, assim como na avaliação da exposição humana à poluição.

Gousseau et al. [34] aplica o método LES para a simulação da dispersão de polu-entes ao redor de uma construção cúbica isolada. Esse tipo de construção provoca grandesvórtices e, por esse motivo, exerce papel fundamental no transporte de massa turbulento.Três anos depois, Gousseau et al. [35], trabalhando sobre uma área urbana atual, realizaum estudo semelhante no centro de Montreal, no Canadá, utilizando a mesma metodolo-gia e considerando duas direções de vento. A meta desse estudo foi averiguar e entenderos processos ambientais como eles ocorrem em situações reais e mais complexas. Comoresultado, destaca-se a influência do posicionamento das construções em relação ao fluxode vento no processo de transporte.

Outra aplicação de LES foi desenvolvida por Kikumoto et al. [36] para realizar adispersão turbulenta de poluentes do ar quimicamente reativos em um cânion urbano bi-dimensional. As reações tiveram influências significantes nas concentrações dos poluentesno espaço do cânion. Além disso, a estrutura da geometria teve grande impacto sobre otransporte dos poluentes.

Aplicando um outro modelo de escoamento turbulento, denominado Detached-Eddy Simulation (DES), Lateb et al. [37] simula numericamente, sobre uma malha estru-turada altamente refinada, o campo de velocidade do vento e o campo de dispersão deum poluente emitido através de uma fonte pontual rodeada de duas construções. Compa-rando os resultados com um outro modelo de CFD, foi evidenciada a mesma média de erroaproximadamente das concentrações calculadas, entretanto o tempo de processamento domodelo adotado foi significativamente maior.

Diferente da maioria das metodologias comumente encontradas, Feng et al. [38]realiza a previsão das concentrações médias diárias de material particulado (𝑃𝑀2.5), pormeio da criação de um modelo híbrido, composto por análises da trajetória da massa dear e pela técnica de transformação wavelet. A trajetória da massa de ar foi empregada

47

para distinguir corredores de transporte de ar limpo e ar poluído em algumas estaçõesde monitoramento selecionadas. Para cada corredor, uma rede de estação triangular foiconstruída. Em seguida, a sequência de tempo original das concentrações do poluente ava-liado foi decomposta pela transformação wavelet em pequenas subsequências com menorvariabilidade, e então, para cada uma delas, a estratégia de previsão foi aplicada. Com opropósito de melhorar a precisão de uma rede neural artificial, a combinação de técnicasdemonstrou ser bastante eficaz.

Como proposta para a diminuição da poluição no ar, Schaap et al. [39] desenvolvecinco modelos de transporte químico (CTMs - Chemistry Transport Models) em três di-mensões, estudados sobre a Europa. Os CTMs são tipicamente aplicados com a finalidadede realizar simulações da química atmosférica. O objetivo principal do trabalho era deter-minar a resolução ótima de malha horizontal, na qual a adição de esforços computacionaisnão proporcionam um desempenho melhor do modelo. O estudo mostrou que o aumentona resolução é algo vantajoso, e as respostas dos CTMs foram bastante consistentes, sendoque a melhor foi encontrada para os poluentes: óxido nítrico (𝑁𝑂2), material particulado(𝑃𝑀10) e ozônio (𝑂3).

Portanto, pode-se perceber que diversas metodologias são aplicáveis neste contextode pesquisa relacionado ao problema de dispersão de poluentes no meio atmosférico. Emdestaque, uma vasta área de conhecimento, denominada como CFD, desempenha um pa-pel fundamental na modelagem desse tipo problema, sendo amplamente empregada devidoao agrupamento de um conjunto de métodos e técnicas que a constituem. Nesse sentido, autilização de um sistema de coordenadas generalizadas como forma de modelagem para aresolução numérica da equação de transporte poderá oferecer um grande diferencial nesseramo, tendo em vista sua capacidade de adaptação à geometria do problema.

Identificados os principais modelos de dispersão de poluentes na atmosfera, noCapítulo 4, apresentam-se os procedimentos metodológicos empregados para o desenvol-vimento do modelo proposto neste trabalho, incluindo sua formulação matemática e suaimplementação.

49

4 PROCEDIMENTOS METODOLÓGICOS

Neste capítulo, aspectos relacionados ao desenvolvimento matemático para a ob-tenção do modelo numérico da equação de transporte e sua implementação são explicadosde maneira detalhada, respectivamente, nas Seções 4.1 e 4.2. É importante salientar queo algoritmo de geração de malhas bidimensionais em coordenadas generalizadas e de re-solução numérica da Equação de Navier-Stokes pertencem aos trabalhos de Cirilo et al.[9] e de Barba [2], fornecendo importantes dados de entrada para o algoritmo de resolu-ção numérica da equação de transporte. Por esse motivo, não se apresentam os demaisdesenvolvimentos nem suas implementações.

4.1 Desenvolvimento Matemático

Nesta seção, desenvolvem-se os cálculos para a obtenção da equação de transporteem sua forma discretizada. Vale ressaltar que o desenvolvimento matemático realizadoneste trabalho foi embasado na dissertação de Barba [2], na qual um processo de de-senvolvimento semelhante é empregado para a Equação de Navier-Stokes. Ambos os de-senvolvimentos englobam a transformação do sistema de coordenadas cartesianas parao de coordenadas generalizadas, a qual é descrita na Subseção 4.1.1, acompanhada peladiscretização dos termos, especificada na Subseção 4.1.2.

Como já mencionado no Capítulo 2, a equação governante para a dispersão depoluentes no meio atmosférico, conhecida como equação de transporte, é representadapela Equação (4.1). Apenas para fins didáticos, essa equação está sendo redescrita nestaseção.

𝜕𝐶

𝜕𝑡⏟ ⏞ Termo temporal

+ 𝜕

𝜕𝑥(𝐶𝑣𝑥) + 𝜕

𝜕𝑦(𝐶𝑣𝑦)⏟ ⏞

Termo convectivo

− 𝜕

𝜕𝑥

(𝐾𝑥

𝜕𝐶

𝜕𝑥

)− 𝜕

𝜕𝑦

(𝐾𝑦

𝜕𝐶

𝜕𝑦

)⏟ ⏞

Termo difusivo

= −𝜅𝐶⏟ ⏞ Termo reativo

, (4.1)

em que:

• 𝐶 é a concentração do poluente;

• 𝑡 refere-se ao tempo;

• 𝑥 e 𝑦 indicam o posicionamento no plano cartesiano;

• 𝑣𝑥 e 𝑣𝑦 são as velocidades dos fluxos de vento da região analisada, respectivamente,nas direções 𝑥 e 𝑦;

• 𝐾𝑥 e 𝐾𝑦 são os coeficientes de difusão, respectivamente, nas direções 𝑥 e 𝑦;

• 𝜅 é o coeficiente de decaimento.

50

4.1.1 Transformação de Coordenadas da Equação Governante

Para que seja possível estudar a dispersão de poluentes sobre geometrias comple-xas, isto é, mais realísticas, é necessário expressar as equações governantes no sistema decoordenadas generalizadas. A princípio, a Equação (4.1) está representada no sistema decoordenadas cartesianas.

Um ponto arbitrário (𝜉, 𝜂) no sistema de coordenadas generalizadas está relacio-nado a um ponto (𝑥, 𝑦) no sistema de coordenadas cartesianas através das equações detransformação, dadas na Equação (4.2).⎧⎪⎪⎪⎨⎪⎪⎪⎩

𝜉 = 𝜉(𝑥, 𝑦, 𝑡),𝜂 = 𝜂(𝑥, 𝑦, 𝑡),𝜏 = 𝜏(𝑡).

(4.2)

Observa-se, na Equação (4.2), que a variável 𝜏 se refere ao tempo e corresponde a umafunção dependente apenas de 𝑡, significando a não existência de malhas móveis, ou seja,malhas cujos pontos se movimentam conforme o tempo (OLIVEIRA et al.[40], 2013 apudBARBA[2], 2015), alterando sua localização.

Para iniciar o processo de transformação de coordenadas, alguns ajustes na Equa-ção (4.1) são realizados, agrupando-se alguns termos de mesma derivada parcial, comomostra a Equação (4.3).

𝜕𝐶

𝜕𝑡+ 𝜕

𝜕𝑥

(𝐶𝑣𝑥 − 𝐾𝑥

𝜕𝐶

𝜕𝑥

)+ 𝜕

𝜕𝑦

(𝐶𝑣𝑦 − 𝐾𝑦

𝜕𝐶

𝜕𝑦

)= −𝜅𝐶. (4.3)

Substituindo-se alguns termos da Equação (4.3) por variáveis, determina-se aEquação (4.4).

𝜕𝐴1

𝜕𝑡+ 𝜕𝐴2

𝜕𝑥+ 𝜕𝐴3

𝜕𝑦= 𝑆, (4.4)

em que

𝐴1 = 𝐶, 𝐴2 = 𝐶𝑣𝑥 − 𝐾𝑥𝜕𝐶

𝜕𝑥, 𝐴3 = 𝐶𝑣𝑦 − 𝐾𝑦

𝜕𝐶

𝜕𝑦, 𝑆 = −𝜅𝐶.

A Equação (4.4) é denominada forma conservativa da equação de transporte em coorde-nadas cartesianas.

Aplicando-se a regra da cadeia nas derivadas parciais da Equação (4.4), tem-seque:

𝜕𝐴1

𝜕𝑡= 𝜕𝐴1

𝜕𝜉

𝜕𝜉

𝜕𝑡+ 𝜕𝐴1

𝜕𝜂

𝜕𝜂

𝜕𝑡+ 𝜕𝐴1

𝜕𝜏, (4.5)

𝜕𝐴2

𝜕𝑥= 𝜕𝐴2

𝜕𝜉

𝜕𝜉

𝜕𝑥+ 𝜕𝐴2

𝜕𝜂

𝜕𝜂

𝜕𝑥, (4.6)

𝜕𝐴3

𝜕𝑦= 𝜕𝐴3

𝜕𝜉

𝜕𝜉

𝜕𝑦+ 𝜕𝐴3

𝜕𝜂

𝜕𝜂

𝜕𝑦. (4.7)

51

Percebe-se que, nas Equações (4.5) - (4.7), algumas simplificações implícitas foram reali-zadas, eliminando-se alguns termos do resultado da regra da cadeia, uma vez que 𝜕𝜏

𝜕𝑡= 1

e 𝜕𝜏𝜕𝑥

= 𝜕𝜏𝜕𝑦

= 0.

Substituindo-se as Equações (4.5) - (4.7) na Equação (4.4), obtém-se a Equação(4.8):(

𝜕𝐴1

𝜕𝜉

𝜕𝜉

𝜕𝑡+ 𝜕𝐴1

𝜕𝜂

𝜕𝜂

𝜕𝑡+ 𝜕𝐴1

𝜕𝜏

)+(

𝜕𝐴2

𝜕𝜉

𝜕𝜉

𝜕𝑥+ 𝜕𝐴2

𝜕𝜂

𝜕𝜂

𝜕𝑥

)+(

𝜕𝐴3

𝜕𝜉

𝜕𝜉

𝜕𝑦+ 𝜕𝐴3

𝜕𝜂

𝜕𝜂

𝜕𝑦

)= 𝑆. (4.8)

Posteriormente, multiplicando-se a Equação (4.8) pelo inverso do Jacobiano 1𝐽, em

que 𝐽 é definido na Subseção 2.2.1, encontra-se a Equação (4.9):

𝜕𝐴1

𝜕𝜉

𝜕𝜉

𝜕𝑡

1𝐽

+ 𝜕𝐴1

𝜕𝜂

𝜕𝜂

𝜕𝑡

1𝐽

+ 𝜕𝐴1

𝜕𝜏

1𝐽

+ 𝜕𝐴2

𝜕𝜉

𝜕𝜉

𝜕𝑥

1𝐽

+ 𝜕𝐴2

𝜕𝜂

𝜕𝜂

𝜕𝑥

1𝐽

+ 𝜕𝐴3

𝜕𝜉

𝜕𝜉

𝜕𝑦

1𝐽

+ 𝜕𝐴3

𝜕𝜂

𝜕𝜂

𝜕𝑦

1𝐽

= 𝑆

𝐽. (4.9)

A partir da Equação (4.9), nota-se que alguns termos se constituem parte doresultado de uma regra do produto. Observe os desenvolvimentos das Equações (4.10) -(4.16):

𝜕

𝜕𝜉

(𝐴1

𝜕𝜉

𝜕𝑡

1𝐽

)= 𝜕𝐴1

𝜕𝜉

𝜕𝜉

𝜕𝑡

1𝐽

+ 𝐴1

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)]

⇒ 𝜕𝐴1

𝜕𝜉

𝜕𝜉

𝜕𝑡

1𝐽

= 𝜕

𝜕𝜉

(𝐴1

𝜕𝜉

𝜕𝑡

1𝐽

)− 𝐴1

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)], (4.10)

𝜕

𝜕𝜂

(𝐴1

𝜕𝜂

𝜕𝑡

1𝐽

)= 𝜕𝐴1

𝜕𝜂

𝜕𝜂

𝜕𝑡

1𝐽

+ 𝐴1

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)]

⇒ 𝜕𝐴1

𝜕𝜂

𝜕𝜂

𝜕𝑡

1𝐽

= 𝜕

𝜕𝜂

(𝐴1

𝜕𝜂

𝜕𝑡

1𝐽

)− 𝐴1

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)], (4.11)

𝜕

𝜕𝜏

(𝐴1

1𝐽

)= 𝜕𝐴1

𝜕𝜏

1𝐽

+ 𝐴1

[𝜕

𝜕𝜏

( 1𝐽

)]

⇒ 𝜕𝐴1

𝜕𝜏

1𝐽

= 𝜕

𝜕𝜏

(𝐴1

1𝐽

)− 𝐴1

[𝜕

𝜕𝜏

( 1𝐽

)], (4.12)

𝜕

𝜕𝜉

(𝐴2

𝜕𝜉

𝜕𝑥

1𝐽

)= 𝜕𝐴2

𝜕𝜉

𝜕𝜉

𝜕𝑥

1𝐽

+ 𝐴2

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)]

⇒ 𝜕𝐴2

𝜕𝜉

𝜕𝜉

𝜕𝑥

1𝐽

= 𝜕

𝜕𝜉

(𝐴2

𝜕𝜉

𝜕𝑥

1𝐽

)− 𝐴2

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)], (4.13)

𝜕

𝜕𝜂

(𝐴2

𝜕𝜂

𝜕𝑥

1𝐽

)= 𝜕𝐴2

𝜕𝜂

𝜕𝜂

𝜕𝑥

1𝐽

+ 𝐴2

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)]

⇒ 𝜕𝐴2

𝜕𝜂

𝜕𝜂

𝜕𝑥

1𝐽

= 𝜕

𝜕𝜂

(𝐴2

𝜕𝜂

𝜕𝑥

1𝐽

)− 𝐴2

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)], (4.14)

𝜕

𝜕𝜉

(𝐴3

𝜕𝜉

𝜕𝑦

1𝐽

)= 𝜕𝐴3

𝜕𝜉

𝜕𝜉

𝜕𝑦

1𝐽

+ 𝐴3

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)]

⇒ 𝜕𝐴3

𝜕𝜉

𝜕𝜉

𝜕𝑦

1𝐽

= 𝜕

𝜕𝜉

(𝐴3

𝜕𝜉

𝜕𝑦

1𝐽

)− 𝐴3

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)], (4.15)

52

𝜕

𝜕𝜂

(𝐴3

𝜕𝜂

𝜕𝑦

1𝐽

)= 𝜕𝐴3

𝜕𝜂

𝜕𝜂

𝜕𝑦

1𝐽

+ 𝐴3

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)]

⇒ 𝜕𝐴3

𝜕𝜂

𝜕𝜂

𝜕𝑦

1𝐽

= 𝜕

𝜕𝜂

(𝐴3

𝜕𝜂

𝜕𝑦

1𝐽

)− 𝐴3

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)]. (4.16)

Dessa forma, substituindo-se as Equações (4.10) - (4.16) na Equação (4.9), define-se a Equação (4.17):

𝜕

𝜕𝜉

(𝐴1

𝜕𝜉

𝜕𝑡

1𝐽

)− 𝐴1

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)]+ 𝜕

𝜕𝜂

(𝐴1

𝜕𝜂

𝜕𝑡

1𝐽

)− 𝐴1

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)]

+ 𝜕

𝜕𝜏

(𝐴1

1𝐽

)− 𝐴1

[𝜕

𝜕𝜏

( 1𝐽

)]+ 𝜕

𝜕𝜉

(𝐴2

𝜕𝜉

𝜕𝑥

1𝐽

)− 𝐴2

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)]

+ 𝜕

𝜕𝜂

(𝐴2

𝜕𝜂

𝜕𝑥

1𝐽

)− 𝐴2

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)]+ 𝜕

𝜕𝜉

(𝐴3

𝜕𝜉

𝜕𝑦

1𝐽

)− 𝐴3

[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)]

+ 𝜕

𝜕𝜂

(𝐴3

𝜕𝜂

𝜕𝑦

1𝐽

)− 𝐴3

[𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)]= 𝑆

𝐽. (4.17)

Separando-se então alguns termos em diferentes lados da equação e realizando-sealguns agrupamentos na Equação (4.17), tem-se que:

𝜕

𝜕𝜉

(𝐴1

𝜕𝜉

𝜕𝑡

1𝐽

)+ 𝜕

𝜕𝜂

(𝐴1

𝜕𝜂

𝜕𝑡

1𝐽

)+ 𝜕

𝜕𝜏

(𝐴1

1𝐽

)

+ 𝜕

𝜕𝜉

(𝐴2

𝜕𝜉

𝜕𝑥

1𝐽

)+ 𝜕

𝜕𝜂

(𝐴2

𝜕𝜂

𝜕𝑥

1𝐽

)+ 𝜕

𝜕𝜉

(𝐴3

𝜕𝜉

𝜕𝑦

1𝐽

)+ 𝜕

𝜕𝜂

(𝐴3

𝜕𝜂

𝜕𝑦

1𝐽

)= 𝑆

𝐽

+𝐴1

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)]+[

𝜕

𝜕𝜏

( 1𝐽

)]}

+𝐴2

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)]}

+𝐴3

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)]}, (4.18)

ou ainda,𝜕

𝜕𝜉

(𝐴1

𝐽

𝜕𝜉

𝜕𝑡+ 𝐴2

𝐽

𝜕𝜉

𝜕𝑥+ 𝐴3

𝐽

𝜕𝜉

𝜕𝑦

)+ 𝜕

𝜕𝜂

(𝐴1

𝐽

𝜕𝜂

𝜕𝑡+ 𝐴2

𝐽

𝜕𝜂

𝜕𝑥+ 𝐴3

𝐽

𝜕𝜂

𝜕𝑦

)+ 𝜕

𝜕𝜏

(𝐴1

𝐽

)= 𝑆

𝐽

+𝐴1

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)]+[

𝜕

𝜕𝜏

( 1𝐽

)]}

+𝐴2

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)]}

+𝐴3

{[𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)]+[

𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)]}. (4.19)

Para continuar o processo de desenvolvimento matemático, alguns cálculos com-plementares são efetuados, com a finalidade de realizar simplificações sobre a Equação(4.19). Aplicando-se a regra da cadeia sobre os termos 𝜕𝜉

𝜕𝜏e 𝜕𝜂

𝜕𝜏, obtém-se que:

𝜕𝜉

𝜕𝜏= 𝜕𝜉

𝜕𝑥

𝜕𝑥

𝜕𝜏+ 𝜕𝜉

𝜕𝑦

𝜕𝑦

𝜕𝜏+ 𝜕𝜉

𝜕𝑡

𝜕𝑡

𝜕𝜏, (4.20)

53

𝜕𝜂

𝜕𝜏= 𝜕𝜂

𝜕𝑥

𝜕𝑥

𝜕𝜏+ 𝜕𝜂

𝜕𝑦

𝜕𝑦

𝜕𝜏+ 𝜕𝜂

𝜕𝑡

𝜕𝑡

𝜕𝜏. (4.21)

Tendo em vista que 𝜕𝜉𝜕𝜏

= 𝜕𝜂𝜕𝜏

= 0 e 𝜕𝑡𝜕𝜏

= 1, as Equações (4.20) e (4.21) implicam que:

𝜕𝜉

𝜕𝑥

𝜕𝑥

𝜕𝜏+ 𝜕𝜉

𝜕𝑦

𝜕𝑦

𝜕𝜏+ 𝜕𝜉

𝜕𝑡= 0, (4.22)

𝜕𝜂

𝜕𝑥

𝜕𝑥

𝜕𝜏+ 𝜕𝜂

𝜕𝑦

𝜕𝑦

𝜕𝜏+ 𝜕𝜂

𝜕𝑡= 0. (4.23)

Finalmente, isolando-se 𝜕𝜉𝜕𝑡

e 𝜕𝜂𝜕𝑡

, respectivamente, das Equações (4.22) e (4.23) e, emseguida, aplicando-se as métricas de transformação, Equação (2.8) mostrada na Subseção2.2.1, encontra-se:

𝜕𝜉

𝜕𝑡= −𝐽

𝜕𝑦

𝜕𝜂

𝜕𝑥

𝜕𝜏+ 𝐽

𝜕𝑥

𝜕𝜂

𝜕𝑦

𝜕𝜏, (4.24)

𝜕𝜂

𝜕𝑡= 𝐽

𝜕𝑦

𝜕𝜉

𝜕𝑥

𝜕𝜏− 𝐽

𝜕𝑥

𝜕𝜉

𝜕𝑦

𝜕𝜏. (4.25)

Retornando-se ao processo de desenvolvimento principal, observa-se, na Equação(4.19), que alguns termos podem ser anulados, conforme segue nas Equações (4.26), (4.27)e (4.28). Nas simplificações, as métricas de transformação, Equação (2.8), e os resultadosobtidos nas Equações (4.24) e (4.25) são utilizados:

𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑡

1𝐽

)+ 𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑡

1𝐽

)+ 𝜕

𝜕𝜏

( 1𝐽

)

= 𝜕

𝜕𝜉

(−𝜕𝑦

𝜕𝜂

𝜕𝑥

𝜕𝜏+ 𝜕𝑥

𝜕𝜂

𝜕𝑦

𝜕𝜏

)+ 𝜕

𝜕𝜂

(𝜕𝑦

𝜕𝜉

𝜕𝑥

𝜕𝜏− 𝜕𝑥

𝜕𝜉

𝜕𝑦

𝜕𝜏

)+ 𝜕

𝜕𝜏

(𝜕𝑥

𝜕𝜉

𝜕𝑦

𝜕𝜂− 𝜕𝑥

𝜕𝜂

𝜕𝑦

𝜕𝜉

)

= − 𝜕2𝑦

𝜕𝜉𝜕𝜂

𝜕𝑥

𝜕𝜏− 𝜕2𝑥

𝜕𝜉𝜕𝜏

𝜕𝑦

𝜕𝜂+ 𝜕2𝑥

𝜕𝜉𝜕𝜂

𝜕𝑦

𝜕𝜏+ 𝜕2𝑦

𝜕𝜉𝜕𝜏

𝜕𝑥

𝜕𝜂+ 𝜕2𝑦

𝜕𝜂𝜕𝜉

𝜕𝑥

𝜕𝜏+ 𝜕2𝑥

𝜕𝜂𝜕𝜏

𝜕𝑦

𝜕𝜉− 𝜕2𝑥

𝜕𝜂𝜕𝜉

𝜕𝑦

𝜕𝜏

− 𝜕2𝑦

𝜕𝜂𝜕𝜏

𝜕𝑥

𝜕𝜉+ 𝜕2𝑥

𝜕𝜏𝜕𝜉

𝜕𝑦

𝜕𝜂+ 𝜕2𝑦

𝜕𝜏𝜕𝜂

𝜕𝑥

𝜕𝜉− 𝜕2𝑥

𝜕𝜏𝜕𝜂

𝜕𝑦

𝜕𝜉− 𝜕2𝑦

𝜕𝜏𝜕𝜉

𝜕𝑥

𝜕𝜂= 0, (4.26)

𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑥

1𝐽

)+ 𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑥

1𝐽

)= 𝜕

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)+ 𝜕

𝜕𝜂

(−𝜕𝑦

𝜕𝜉

)= 𝜕2𝑦

𝜕𝜉𝜕𝜂− 𝜕2𝑦

𝜕𝜉𝜕𝜂= 0, (4.27)

𝜕

𝜕𝜉

(𝜕𝜉

𝜕𝑦

1𝐽

)+ 𝜕

𝜕𝜂

(𝜕𝜂

𝜕𝑦

1𝐽

)= 𝜕

𝜕𝜉

(−𝜕𝑥

𝜕𝜂

)+ 𝜕

𝜕𝜂

(𝜕𝑥

𝜕𝜉

)= − 𝜕2𝑥

𝜕𝜉𝜕𝜂+ 𝜕2𝑥

𝜕𝜉𝜕𝜂= 0. (4.28)

Assim, a partir das Equações (4.26) - (4.28), simplifica-se a Equação (4.19), comoé expresso na Equação (4.29).

𝜕

𝜕𝜏

(𝐴1

𝐽

)+ 𝜕

𝜕𝜉

(𝐴1

𝐽

𝜕𝜉

𝜕𝑡+ 𝐴2

𝐽

𝜕𝜉

𝜕𝑥+ 𝐴3

𝐽

𝜕𝜉

𝜕𝑦

)+ 𝜕

𝜕𝜂

(𝐴1

𝐽

𝜕𝜂

𝜕𝑡+ 𝐴2

𝐽

𝜕𝜂

𝜕𝑥+ 𝐴3

𝐽

𝜕𝜂

𝜕𝑦

)= 𝑆

𝐽. (4.29)

Substituindo-se alguns termos da Equação (4.29) por variáveis, determina-se aEquação (4.30).

𝜕𝐵1

𝜕𝜏+ 𝜕𝐵2

𝜕𝜉+ 𝜕𝐵3

𝜕𝜂= 𝑇, (4.30)

54

na qual

𝐵1 = 𝐴1

𝐽, 𝐵2 = 𝐴1

𝐽

𝜕𝜉

𝜕𝑡+𝐴2

𝐽

𝜕𝜉

𝜕𝑥+𝐴3

𝐽

𝜕𝜉

𝜕𝑦, 𝐵3 = 𝐴1

𝐽

𝜕𝜂

𝜕𝑡+𝐴2

𝐽

𝜕𝜂

𝜕𝑥+𝐴3

𝐽

𝜕𝜂

𝜕𝑦, 𝑇 = 𝑆

𝐽.

A Equação (4.30) é denominada forma conservativa da equação de transporte em coorde-nadas generalizadas.

Substituindo-se as variáveis definidas na Equação (4.4) sobre a Equação (4.29),obtém-se que:

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉

(𝐶

𝐽

𝜕𝜉

𝜕𝑡+(

𝐶𝑣𝑥

𝐽− 𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

)𝜕𝜉

𝜕𝑥+(

𝐶𝑣𝑦

𝐽− 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

)𝜕𝜉

𝜕𝑦

)

+ 𝜕

𝜕𝜂

(𝐶

𝐽

𝜕𝜂

𝜕𝑡+(

𝐶𝑣𝑥

𝐽− 𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

)𝜕𝜂

𝜕𝑥+(

𝐶𝑣𝑦

𝐽− 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

)𝜕𝜂

𝜕𝑦

)= −𝜅𝐶

𝐽, (4.31)

ou ainda,

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉

(𝐶

𝐽

(𝜕𝜉

𝜕𝑡+ 𝑣𝑥

𝜕𝜉

𝜕𝑥+ 𝑣𝑦

𝜕𝜉

𝜕𝑦

)− 𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

𝜕𝜉

𝜕𝑥− 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

𝜕𝜉

𝜕𝑦

)

+ 𝜕

𝜕𝜂

(𝐶

𝐽

(𝜕𝜂

𝜕𝑡+ 𝑣𝑥

𝜕𝜂

𝜕𝑥+ 𝑣𝑦

𝜕𝜂

𝜕𝑦

)− 𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

𝜕𝜂

𝜕𝑥− 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

𝜕𝜂

𝜕𝑦

)= −𝜅𝐶

𝐽. (4.32)

Separando-se alguns termos das derivadas parciais na Equação (4.32), encontra-sea Equação (4.33).

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉

(𝐶

𝐽

(𝜕𝜉

𝜕𝑡+ 𝑣𝑥

𝜕𝜉

𝜕𝑥+ 𝑣𝑦

𝜕𝜉

𝜕𝑦

))+ 𝜕

𝜕𝜂

(𝐶

𝐽

(𝜕𝜂

𝜕𝑡+ 𝑣𝑥

𝜕𝜂

𝜕𝑥+ 𝑣𝑦

𝜕𝜂

𝜕𝑦

))

=−𝜅𝐶

𝐽+ 𝜕

𝜕𝜉

(𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

𝜕𝜉

𝜕𝑥+ 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

𝜕𝜉

𝜕𝑦

)+ 𝜕

𝜕𝜂

(𝐾𝑥

𝐽

𝜕𝐶

𝜕𝑥

𝜕𝜂

𝜕𝑥+ 𝐾𝑦

𝐽

𝜕𝐶

𝜕𝑦

𝜕𝜂

𝜕𝑦

). (4.33)

Sabendo-se que a concentração 𝐶 = 𝐶(𝑥, 𝑦), isto é, uma função dependente de 𝑥

e 𝑦, aplica-se a regra da cadeia nos termos 𝜕𝐶𝜕𝑥

e 𝜕𝐶𝜕𝑦

da Equação (4.33), obtendo-se:

𝜕𝐶

𝜕𝑥= 𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑥+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑥, (4.34)

𝜕𝐶

𝜕𝑦= 𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑦+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑦. (4.35)

Substituindo-se as Equações (4.34) e (4.35) em alguns termos da Equação (4.33), segueque:

𝜕𝐶

𝜕𝑥

𝜕𝜉

𝜕𝑥=(

𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑥+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑥

)𝜕𝜉

𝜕𝑥= 𝜕𝐶

𝜕𝜉

(𝜕𝜉

𝜕𝑥

)2

+ 𝜕𝐶

𝜕𝜂

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥, (4.36)

𝜕𝐶

𝜕𝑦

𝜕𝜉

𝜕𝑦=(

𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑦+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑦

)𝜕𝜉

𝜕𝑦= 𝜕𝐶

𝜕𝜉

(𝜕𝜉

𝜕𝑦

)2

+ 𝜕𝐶

𝜕𝜂

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦, (4.37)

𝜕𝐶

𝜕𝑥

𝜕𝜂

𝜕𝑥=(

𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑥+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑥

)𝜕𝜂

𝜕𝑥= 𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥+ 𝜕𝐶

𝜕𝜂

(𝜕𝜂

𝜕𝑥

)2

, (4.38)

55

𝜕𝐶

𝜕𝑦

𝜕𝜂

𝜕𝑦=(

𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑦+ 𝜕𝐶

𝜕𝜂

𝜕𝜂

𝜕𝑦

)𝜕𝜂

𝜕𝑦= 𝜕𝐶

𝜕𝜉

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦+ 𝜕𝐶

𝜕𝜂

(𝜕𝜂

𝜕𝑦

)2

. (4.39)

Ainda, dadas as Equações (4.36) - (4.39), expandem-se alguns termos utilizando as mé-tricas de transformação, Equação (2.8), como demonstram as Equações (4.40) - (4.45).(

𝜕𝜉

𝜕𝑥

)2

=(

𝐽𝜕𝑦

𝜕𝜂

)2

= 𝐽2(

𝜕𝑦

𝜕𝜂

)2

, (4.40)

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥= 𝐽

𝜕𝑦

𝜕𝜂

(−𝐽

𝜕𝑦

𝜕𝜉

)= −𝐽2 𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉, (4.41)(

𝜕𝜉

𝜕𝑦

)2

=(

−𝐽𝜕𝑥

𝜕𝜂

)2

= 𝐽2(

𝜕𝑥

𝜕𝜂

)2

, (4.42)

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦= −𝐽

𝜕𝑥

𝜕𝜂𝐽

𝜕𝑥

𝜕𝜉= −𝐽2 𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉, (4.43)(

𝜕𝜂

𝜕𝑥

)2

=(

−𝐽𝜕𝑦

𝜕𝜉

)2

= 𝐽2(

𝜕𝑦

𝜕𝜉

)2

, (4.44)(𝜕𝜂

𝜕𝑦

)2

=(

𝐽𝜕𝑥

𝜕𝜉

)2

= 𝐽2(

𝜕𝑥

𝜕𝜉

)2

. (4.45)

Substituindo-se as Equações (4.40) - (4.45) em (4.36) - (4.39) e, em seguida,aplicando-se os resultados obtidos sobre a Equação (4.33), encontra-se a Equação (4.46).

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉

(𝐶

𝐽

(𝜕𝜉

𝜕𝑡+ 𝑣𝑥

𝜕𝜉

𝜕𝑥+ 𝑣𝑦

𝜕𝜉

𝜕𝑦

))+ 𝜕

𝜕𝜂

(𝐶

𝐽

(𝜕𝜂

𝜕𝑡+ 𝑣𝑥

𝜕𝜂

𝜕𝑥+ 𝑣𝑦

𝜕𝜂

𝜕𝑦

))= −𝜅𝐶

𝐽

+ 𝜕

𝜕𝜉

⎛⎝𝐾𝑥

𝐽

⎛⎝𝜕𝐶

𝜕𝜉𝐽2(

𝜕𝑦

𝜕𝜂

)2

+ 𝜕𝐶

𝜕𝜂

(−𝐽2 𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)⎞⎠+𝐾𝑦

𝐽

⎛⎝𝜕𝐶

𝜕𝜉𝐽2(

𝜕𝑥

𝜕𝜂

)2

+ 𝜕𝐶

𝜕𝜂

(−𝐽2 𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)⎞⎠⎞⎠+ 𝜕

𝜕𝜂

⎛⎝𝐾𝑥

𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(−𝐽2 𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)+ 𝜕𝐶

𝜕𝜂𝐽2(

𝜕𝑦

𝜕𝜉

)2⎞⎠

+𝐾𝑦

𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(−𝐽2 𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)+ 𝜕𝐶

𝜕𝜂𝐽2(

𝜕𝑥

𝜕𝜉

)2⎞⎠⎞⎠ . (4.46)

Finalmente, realizando alguns ajustes na Equação (4.46), determina-se a equaçãoque governa a dispersão de poluentes na atmosfera, descrita no sistema de coordenadasgeneralizadas e representada pela Equação (4.47).

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉

(𝐶

𝐽

(𝜕𝜉

𝜕𝑡+ 𝑣𝑥

𝜕𝜉

𝜕𝑥+ 𝑣𝑦

𝜕𝜉

𝜕𝑦

))+ 𝜕

𝜕𝜂

(𝐶

𝐽

(𝜕𝜂

𝜕𝑡+ 𝑣𝑥

𝜕𝜂

𝜕𝑥+ 𝑣𝑦

𝜕𝜂

𝜕𝑦

))= −𝜅𝐶

𝐽

+𝐾𝑥

⎡⎣ 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠+ 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑦

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠⎤⎦+𝐾𝑦

⎡⎣ 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑥

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠+ 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑥

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠⎤⎦ .

(4.47)

56

Com base na Equação (4.47), o próximo passo consiste na discretização dos termosda equação de transporte, para que se torne possível de ser manipulada pelo computador.Na Subseção 4.1.2, efetua-se essa discretização tanto no sistema de coordenadas cartesi-anas quanto no sistema de coordenadas generalizadas.

4.1.2 Discretização dos Termos da Equação Governante

Antes de se iniciar o processo de discretização, é importante levar em consideraçãoque, neste trabalho, malhas do tipo deslocada são aplicadas, isto é, aquelas em que asvariáveis do problema, no caso a concentração 𝐶 do poluente, são armazenadas em po-sições diferentes das posições comumente utilizadas [15]. A aplicação dessa metodologiaé necessária pelo fato de que possibilita a redução da instabilidade da solução numérica(FERREIRA[41], 2001 apud BARBA[2], 2015).

Na Figura 7, ilustra-se, à esquerda, a nomenclatura das posições de uma malhahipotética, derivada dos pontos cardeais em inglês. À direita, apresenta-se o armazena-mento da variável de pressão 𝑝 no centro da célula e das velocidades 𝑢 e 𝑣 nos centros dasfaces para o caso da equação de Navier-Stokes. Para a equação de transporte, a concen-tração 𝐶 ocupa o lugar da pressão 𝑝, e a posição das velocidades, definidas como 𝑣𝑥 e 𝑣𝑦

respectivamente, não se altera.

Figura 7 – Elementos de malha - retirada de Barba [2]

A apresentação das nomenclaturas das localizações na malha, de acordo com aFigura 7, constitui-se bastante importante para as discretizações realizadas nas próximasduas subseções. Em ambas as discretizações, efetuam-se as aproximações a partir doelemento de centro 𝑃 , separando-se os termos: temporal, difusivo e advectivo (convectivo).

57

4.1.2.1 Discretização em Coordenadas Cartesianas

Nesta subseção, realiza-se a discretização da equação governante em coordenadascartesianas apenas para fins didáticos. Em um primeiro momento, é necessário efetuaralguns ajustes na Equação (4.1), obtendo-se a Equação (4.48).

𝜕𝐶

𝜕𝑡+ 𝜕

𝜕𝑥(𝐶𝑣𝑥) + 𝜕

𝜕𝑦(𝐶𝑣𝑦) = −𝜅𝐶 + 𝐾𝑥

𝜕2𝐶

𝜕𝑥2 + 𝐾𝑦𝜕2𝐶

𝜕𝑦2 . (4.48)

A aproximação para o termo temporal 𝜕𝐶𝜕𝑡

é realizada por diferenças finitas pro-gressivas para a primeira derivada, como mostra a Equação (4.49).

𝜕𝐶

𝜕𝑡

𝑘+1

𝑃

≈ 𝐶|𝑘+1𝑃 − 𝐶|𝑘𝑃

Δ𝑡. (4.49)

O termo difusivo, definido como 𝒟(𝐶), é discretizado por diferenças finitas centraispara a segunda derivada, de acordo com a Equação (4.50).

𝒟(𝐶)|𝑘𝑃 = 𝐾𝑥𝜕2𝐶

𝜕𝑥2

𝑘

𝑃

+ 𝐾𝑦𝜕2𝐶

𝜕𝑦2

𝑘

𝑃

≈ 𝐾𝑥

(𝐶|𝑘𝑊 − 2𝐶|𝑘𝑃 + 𝐶|𝑘𝐸

(Δ𝑥)2

)+ 𝐾𝑦

(𝐶|𝑘𝑆 − 2𝐶|𝑘𝑃 + 𝐶|𝑘𝑁

(Δ𝑦)2

). (4.50)

O termo advectivo, nomeado como 𝒜(𝐶), é aproximado por diferenças finitas cen-trais para a primeira derivada, segundo a Equação (4.51).

𝒜(𝐶)|𝑘𝑃 = 𝜕

𝜕𝑥(𝐶𝑣𝑥)

𝑘

𝑃

+ 𝜕

𝜕𝑦(𝐶𝑣𝑦)

𝑘

𝑃

≈ 𝑣𝑥|𝑒𝐶|𝑘𝑒 − 𝑣𝑥|𝑤𝐶|𝑘𝑤2(

Δ𝑥2

) + 𝑣𝑦|𝑛𝐶|𝑘𝑛 − 𝑣𝑦|𝑠𝐶|𝑘𝑠2(

Δ𝑦2

)= 𝑣𝑥|𝑒𝐶|𝑘𝑒 − 𝑣𝑥|𝑤𝐶|𝑘𝑤

Δ𝑥+ 𝑣𝑦|𝑛𝐶|𝑘𝑛 − 𝑣𝑦|𝑠𝐶|𝑘𝑠

Δ𝑦. (4.51)

Percebe-se na Equação (4.51) que os termos 𝐶|𝑘𝑒 , 𝐶|𝑘𝑤, 𝐶|𝑘𝑛 e 𝐶|𝑘𝑠 não estão localizados noscentros da célula da malha computacional. Por esse motivo, o esquema FOU é aplicado,com a finalidade de evitar instabilidades numéricas, conforme apresentam as Equações(4.52) - (4.55).

𝐶|𝑘𝑒 =(

1 + ℛ|𝑘𝑒2

)𝐶|𝑘𝑃 +

(1 − ℛ|𝑘𝑒

2

)𝐶|𝑘𝐸, ℛ|𝑘𝑒 =

⎧⎨⎩ 1, 𝑣𝑥|𝑒 ≥ 0−1, 𝑣𝑥|𝑒 < 0

, (4.52)

𝐶|𝑘𝑤 =(

1 + ℛ|𝑘𝑤2

)𝐶|𝑘𝑊 +

(1 − ℛ|𝑘𝑤

2

)𝐶|𝑘𝑃 , ℛ|𝑘𝑤 =

⎧⎨⎩ 1, 𝑣𝑥|𝑤 ≥ 0−1, 𝑣𝑥|𝑤 < 0

, (4.53)

𝐶|𝑘𝑛 =(

1 + ℛ|𝑘𝑛2

)𝐶|𝑘𝑃 +

(1 − ℛ|𝑘𝑛

2

)𝐶|𝑘𝑁 , ℛ|𝑘𝑛 =

⎧⎨⎩ 1, 𝑣𝑦|𝑛 ≥ 0−1, 𝑣𝑦|𝑛 < 0

, (4.54)

𝐶|𝑘𝑠 =(

1 + ℛ|𝑘𝑠2

)𝐶|𝑘𝑆 +

(1 − ℛ|𝑘𝑠

2

)𝐶|𝑘𝑃 , ℛ|𝑘𝑠 =

⎧⎨⎩ 1, 𝑣𝑦|𝑠 ≥ 0−1, 𝑣𝑦|𝑠 < 0

. (4.55)

58

Portanto, dos resultados presentes nas Equações (4.49) - (4.51), substituem-se ostermos na Equação (4.48), obtendo-se:

𝜕𝐶

𝜕𝑡

𝑘+1

𝑃

+ 𝒜(𝐶)|𝑘𝑃 = −𝜅𝐶|𝑘𝑃 + 𝒟(𝐶)|𝑘𝑃 , (4.56)

ou ainda,𝐶|𝑘+1

𝑃 − 𝐶|𝑘𝑃Δ𝑡

+ 𝒜(𝐶)|𝑘𝑃 = −𝜅𝐶|𝑘𝑃 + 𝒟(𝐶)|𝑘𝑃 . (4.57)

Isolando-se o termo 𝐶|𝑘+1𝑃 na Equação (4.57), encontra-se:

𝐶|𝑘+1𝑃 = Δ𝑡(−𝜅𝐶|𝑘𝑃 + 𝒟(𝐶)|𝑘𝑃 − 𝒜(𝐶)|𝑘𝑃 ) + 𝐶|𝑘𝑃 . (4.58)

A Equação (4.58) representa a equação resultante da discretização da equação de trans-porte no sistema de coordenadas cartesianas.

Na Subseção 4.1.2.2, demonstra-se então a discretização, em coordenadas genera-lizadas, da equação de transporte.

4.1.2.2 Discretização em Coordenadas Generalizadas

Nesta subseção, realiza-se a discretização da equação governante em coordenadasgeneralizadas. A princípio, nomeiam-se as variáveis 𝑈 e 𝑉 , conhecidas como componentescontravariantes para o vetor velocidade, aplicando-as sobre a Equação (4.47):

𝜕

𝜕𝜏

(𝐶

𝐽

)+ 𝜕

𝜕𝜉(𝑈𝐶) + 𝜕

𝜕𝜂(𝑉 𝐶) = −𝜅𝐶

𝐽

+𝐾𝑥

⎡⎣ 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠+ 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑦

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠⎤⎦+𝐾𝑦

⎡⎣ 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑥

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠+ 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑥

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠⎤⎦ ,

(4.59)

em que

𝑈 = 1𝐽

(𝜕𝜉

𝜕𝑡+ 𝑣𝑥

𝜕𝜉

𝜕𝑥+ 𝑣𝑦

𝜕𝜉

𝜕𝑦

),

𝑉 = 1𝐽

(𝜕𝜂

𝜕𝑡+ 𝑣𝑥

𝜕𝜂

𝜕𝑥+ 𝑣𝑦

𝜕𝜂

𝜕𝑦

).

A aproximação para o termo temporal 𝜕𝜕𝜏

(𝐶𝐽

)é realizada por diferenças finitas

progressivas para a primeira derivada, como mostra a Equação (4.60):

𝜕

𝜕𝜏

(𝐶

𝐽

)𝑘+1

𝑃

≈

(𝐶𝐽

)𝑘+1

𝑃−(

𝐶𝐽

)𝑘𝑃

Δ𝜏=

(1

𝐽 |𝑃

)𝐶|𝑘+1

𝑃 −(

1𝐽 |𝑃

)𝐶|𝑘𝑃

Δ𝜏. (4.60)

59

O termo difusivo, definido como D(𝐶), é subdividido em quatro termos – D1(𝐶),D2(𝐶), D3(𝐶) e D4(𝐶) – para facilitar os cálculos e melhorar a visualização, conforme aEquação (4.61):

D(𝐶)|𝑘𝑃 = 𝐾𝑥(D1(𝐶)|𝑘𝑃 + D2(𝐶)|𝑘𝑃 ) + 𝐾𝑦(D3(𝐶)|𝑘𝑃 + D4(𝐶)|𝑘𝑃 ). (4.61)

Os termos D1(𝐶), D2(𝐶), D3(𝐶) e D4(𝐶) são discretizados por diferenças finitas centraispara a primeira derivada, de acordo com as Equações (4.62) - (4.65):

D1(𝐶)|𝑘𝑃 = 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑃

≈ 12Δ𝜉

⎡⎢⎣⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝐸

−

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑊

⎤⎥⎦=1

2

⎡⎢⎣𝐽 |𝐸

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠𝑘

𝐸

− 𝐽 |𝑊

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑦

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠𝑘

𝑊

⎤⎥⎦=1

2

⎡⎣⎛⎝𝐽

(𝜕𝑦

𝜕𝜂

)2⎞⎠

𝐸

(𝜕𝐶

𝜕𝜉

)𝑘

𝐸

−(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝐸

(𝜕𝐶

𝜕𝜂

)𝑘

𝐸

−

⎛⎝𝐽

(𝜕𝑦

𝜕𝜂

)2⎞⎠

𝑊

(𝜕𝐶

𝜕𝜉

)𝑘

𝑊

+(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑊

(𝜕𝐶

𝜕𝜂

)𝑘

𝑊

⎤⎦=1

2

⎡⎣⎛⎝𝐽

(𝜕𝑦

𝜕𝜂

)2⎞⎠

𝐸

𝛿1 −(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝐸

𝛿2 −

⎛⎝𝐽

(𝜕𝑦

𝜕𝜂

)2⎞⎠

𝑊

𝛿3 +(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑊

𝛿4

⎤⎦ ,

(4.62)

D2(𝐶)|𝑘𝑃 = 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑦

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑃

≈ 12Δ𝜂

⎡⎢⎣⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑦

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑁

−

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑦

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑆

⎤⎥⎦=1

2

⎡⎣⎛⎝𝐽

(𝜕𝑦

𝜕𝜉

)2⎞⎠

𝑁

(𝜕𝐶

𝜕𝜂

)𝑘

𝑁

−(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑁

(𝜕𝐶

𝜕𝜉

)𝑘

𝑁

−

⎛⎝𝐽

(𝜕𝑦

𝜕𝜉

)2⎞⎠

𝑆

(𝜕𝐶

𝜕𝜂

)𝑘

𝑆

+(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑆

(𝜕𝐶

𝜕𝜉

)𝑘

𝑆

⎤⎦=1

2

⎡⎣⎛⎝𝐽

(𝜕𝑦

𝜕𝜉

)2⎞⎠

𝑁

𝛿5 −(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑁

𝛿6 −

⎛⎝𝐽

(𝜕𝑦

𝜕𝜉

)2⎞⎠

𝑆

𝛿7 +(

𝐽𝜕𝑦

𝜕𝜂

𝜕𝑦

𝜕𝜉

)𝑆

𝛿8

⎤⎦ , (4.63)

D3(𝐶)|𝑘𝑃 = 𝜕

𝜕𝜉

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜉

(𝜕𝑥

𝜕𝜂

)2

− 𝜕𝐶

𝜕𝜂

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑃

≈12

⎡⎣⎛⎝𝐽

(𝜕𝑥

𝜕𝜂

)2⎞⎠

𝐸

𝛿1 −(

𝐽𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)𝐸

𝛿2 −

⎛⎝𝐽

(𝜕𝑥

𝜕𝜂

)2⎞⎠

𝑊

𝛿3 +(

𝐽𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)𝑊

𝛿4

⎤⎦ ,

(4.64)

60

D4(𝐶)|𝑘𝑃 = 𝜕

𝜕𝜂

⎛⎝𝐽

⎛⎝𝜕𝐶

𝜕𝜂

(𝜕𝑥

𝜕𝜉

)2

− 𝜕𝐶

𝜕𝜉

𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

⎞⎠⎞⎠𝑘

𝑃

≈12

⎡⎣⎛⎝𝐽

(𝜕𝑥

𝜕𝜉

)2⎞⎠

𝑁

𝛿5 −(

𝐽𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)𝑁

𝛿6 −

⎛⎝𝐽

(𝜕𝑥

𝜕𝜉

)2⎞⎠

𝑆

𝛿7 +(

𝐽𝜕𝑥

𝜕𝜂

𝜕𝑥

𝜕𝜉

)𝑆

𝛿8

⎤⎦ ,

(4.65)

em que

𝛿1 =(

𝜕𝐶

𝜕𝜉

)𝑘

𝐸

≈

⎧⎨⎩𝐶|𝑘𝐸𝐸𝐸−𝐶|𝑘𝑃

2 (diferenças finitas centrais)𝐶|𝑘𝐸 − 𝐶|𝑘𝑃 (diferenças finitas regressivas)

,

𝛿2 =(

𝜕𝐶

𝜕𝜂

)𝑘

𝐸

≈

⎧⎪⎪⎪⎨⎪⎪⎪⎩𝐶|𝑘𝑁𝐸−𝐶|𝑘𝑆𝐸

2 (diferenças finitas centrais)𝐶|𝑘𝑁𝐸 − 𝐶|𝑘𝐸 (diferenças finitas progressivas)𝐶|𝑘𝐸 − 𝐶|𝑘𝑆𝐸 (diferenças finitas regressivas)

,

𝛿3 =(

𝜕𝐶

𝜕𝜉

)𝑘

𝑊

≈

⎧⎨⎩𝐶|𝑘𝑃 −𝐶|𝑘𝑊 𝑊 𝑊

2 (diferenças finitas centrais)𝐶|𝑘𝑃 − 𝐶|𝑘𝑊 (diferenças finitas progressivas)

,

𝛿4 =(

𝜕𝐶

𝜕𝜂

)𝑘

𝑊

≈

⎧⎪⎪⎪⎨⎪⎪⎪⎩𝐶|𝑘𝑁𝑊 −𝐶|𝑘𝑆𝑊

2 (diferenças finitas centrais)𝐶|𝑘𝑁𝑊 − 𝐶|𝑘𝑊 (diferenças finitas progressivas)𝐶|𝑘𝑊 − 𝐶|𝑘𝑆𝑊 (diferenças finitas regressivas)

,

𝛿5 =(

𝜕𝐶

𝜕𝜂

)𝑘

𝑁

≈

⎧⎨⎩𝐶|𝑘𝑁𝑁𝑁 −𝐶|𝑘𝑃

2 (diferenças finitas centrais)𝐶|𝑘𝑁 − 𝐶|𝑘𝑃 (diferenças finitas regressivas)

,

𝛿6 =(

𝜕𝐶

𝜕𝜉

)𝑘

𝑁

≈

⎧⎪⎪⎪⎨⎪⎪⎪⎩𝐶|𝑘𝑁𝐸−𝐶|𝑘𝑁𝑊

2 (diferenças finitas centrais)𝐶|𝑘𝑁𝐸 − 𝐶|𝑘𝑁 (diferenças finitas progressivas)𝐶|𝑘𝑁 − 𝐶|𝑘𝑁𝑊 (diferenças finitas regressivas)

,

𝛿7 =(

𝜕𝐶

𝜕𝜂

)𝑘

𝑆

≈

⎧⎨⎩𝐶|𝑘𝑃 −𝐶|𝑘𝑆𝑆𝑆

2 (diferenças finitas centrais)𝐶|𝑘𝑃 − 𝐶|𝑘𝑆 (diferenças finitas progressivas)

,

𝛿8 =(

𝜕𝐶

𝜕𝜉

)𝑘

𝑆

≈

⎧⎪⎪⎪⎨⎪⎪⎪⎩𝐶|𝑘𝑆𝐸−𝐶|𝑘𝑆𝑊

2 (diferenças finitas centrais)𝐶|𝑘𝑆𝐸 − 𝐶|𝑘𝑆 (diferenças finitas progressivas)𝐶|𝑘𝑆 − 𝐶|𝑘𝑆𝑊 (diferenças finitas regressivas)

.

Pode-se notar, nas Equações (4.62) - (4.65), que as variáveis 𝛿𝑖 definidas, onde 𝑖 = 1, . . . , 8,podem ser discretizadas de diferentes maneiras. O modo de escolha de qual discretizaçãoé efetuada depende diretamente da localização da célula na malha, na qual a concentra-ção está sendo calculada. Vale ressaltar que se tem como prioridade a aproximação pordiferenças finitas centrais, quando possível de ser aplicada.

O termo advectivo, nomeado como A (𝐶), é aproximado por diferenças finitascentrais para a primeira derivada, segundo a Equação (4.66):

A (𝐶)|𝑘𝑃 = 𝜕

𝜕𝜉(𝑈𝐶)

𝑘

𝑃

+ 𝜕

𝜕𝜂(𝑉 𝐶)

𝑘

𝑃

≈ 𝑈 |𝑒𝐶|𝑘𝑒 − 𝑈 |𝑤𝐶|𝑘𝑤2(

Δ𝜉2

) + 𝑉 |𝑛𝐶|𝑘𝑛 − 𝑉 |𝑠𝐶|𝑘𝑠2(

Δ𝜂2

)

61

= 𝑈 |𝑒𝐶|𝑘𝑒 − 𝑈 |𝑤𝐶|𝑘𝑤 + 𝑉 |𝑛𝐶|𝑘𝑛 − 𝑉 |𝑠𝐶|𝑘𝑠 . (4.66)

Observa-se que os termos 𝐶|𝑘𝑒 , 𝐶|𝑘𝑤, 𝐶|𝑘𝑛 e 𝐶|𝑘𝑠 não estão localizados nos centros da célulada malha computacional. Por essa razão, o esquema FOU é aplicado, conforme apresentamas Equações (4.67) - (4.70):

𝐶|𝑘𝑒 =(

1 + S |𝑘𝑒2

)𝐶|𝑘𝑃 +

(1 − S |𝑘𝑒

2

)𝐶|𝑘𝐸, S |𝑘𝑒 =

⎧⎨⎩ 1, 𝑈 |𝑒 ≥ 0−1, 𝑈 |𝑒 < 0

, (4.67)

𝐶|𝑘𝑤 =(

1 + S |𝑘𝑤2

)𝐶|𝑘𝑊 +

(1 − S |𝑘𝑤

2

)𝐶|𝑘𝑃 , S |𝑘𝑤 =

⎧⎨⎩ 1, 𝑈 |𝑤 ≥ 0−1, 𝑈 |𝑤 < 0

, (4.68)

𝐶|𝑘𝑛 =(

1 + S |𝑘𝑛2

)𝐶|𝑘𝑃 +

(1 − S |𝑘𝑛

2

)𝐶|𝑘𝑁 , S |𝑘𝑛 =

⎧⎨⎩ 1, 𝑉 |𝑛 ≥ 0−1, 𝑉 |𝑛 < 0

, (4.69)

𝐶|𝑘𝑠 =(

1 + S |𝑘𝑠2

)𝐶|𝑘𝑆 +

(1 − S |𝑘𝑠

2

)𝐶|𝑘𝑃 , S |𝑘𝑠 =

⎧⎨⎩ 1, 𝑉 |𝑠 ≥ 0−1, 𝑉 |𝑠 < 0

. (4.70)

Portanto, das Equações (4.60), (4.61) e (4.66), substituem-se os termos na Equação(4.59), obtendo-se:

𝜕

𝜕𝜏

(𝐶

𝐽

)𝑘+1

𝑃

+ A (𝐶)|𝑘𝑃 = −𝜅𝐶|𝑘𝑃𝐽 |𝑃

+ D(𝐶)|𝑘𝑃 , (4.71)

ou ainda, (1

𝐽 |𝑃

)𝐶|𝑘+1

𝑃 −(

1𝐽 |𝑃

)𝐶|𝑘𝑃

Δ𝜏+ A (𝐶)|𝑘𝑃 = −𝜅𝐶|𝑘𝑃

𝐽 |𝑃+ D(𝐶)|𝑘𝑃 . (4.72)

Isolando-se o termo 𝐶|𝑘+1𝑃 na Equação (4.72), encontra-se:

𝐶|𝑘+1𝑃 = 𝐽 |𝑃

[Δ𝜏

(−𝜅𝐶|𝑘𝑃

𝐽 |𝑃+ D(𝐶)|𝑘𝑃 − A (𝐶)|𝑘𝑃

)+(

1𝐽 |𝑃

)𝐶|𝑘𝑃

], (4.73)

ou ainda,

𝐶|𝑘+1𝑃 = 𝐽 |𝑃 Δ𝜏

(−𝜅𝐶|𝑘𝑃

𝐽 |𝑃+ D(𝐶)|𝑘𝑃 − A (𝐶)|𝑘𝑃

)+ 𝐶|𝑘𝑃 . (4.74)

A Equação (4.74) representa a equação resultante da discretização da equação de trans-porte no sistema de coordenadas generalizadas. Essa equação é empregada na implemen-tação do modelo numérico para a dispersão de poluentes na atmosfera, descrita na Seção4.2.

4.2 Implementação

Nesta seção, a implementação da Equação (4.74) é detalhada, na qual se empregoua linguagem Fortran 901 como forma de desenvolvimento, uma vez que o foco principal1 <http://www.fortran.com/>

62

desta proposta está embasada no alto desempenho e na precisão dos resultados, buscandocoerência com a situação real.

Como visão geral da implementação, representa-se, na Figura 8, um fluxograma domodelo adotado neste trabalho, com todos os passos do algoritmo. A separação em corestem apenas fins didáticos, na qual a cor verde retrata os módulos de leitura e alocação dememória; a cor vermelha, os módulos principais deste modelo; e a cor azul, os módulosde escrita e desalocação de memória. Atenção especial é dada aos módulos descritos pelacor vermelha, explicados minuciosamente nos próximos parágrafos.

Figura 8 – Fluxograma do modelo de dispersão de poluentes na atmosfera

De maneira codificada, o fluxograma da Figura 8 é representado pelo trecho decódigo no Algoritmo 4.1. Cada chamada (CALL) executa uma subrotina passando os pa-râmetros especificados e indicando cada passo do fluxograma. Os passos de alocação edesalocação de memória não podem ser executados em subrotinas e, por essa razão, sãoexecutados no próprio programa principal.

63

1 ! LEITURA DOS NOMES DOS ARQUIVOS PASSADOS VIA LINHA DE COMANDO2 CALL READ_FNAMES (FNAME)34 ! LEITURA DE INFORMACOES PREVIAS5 CALL READ_INFO (FNAME (1) , KAPPA , KX , KY , C_PR , CSI , ETA , TAU , TF)67 ! ALOCACAO DE MEMORIA8 ALLOCATE (X(CSI + 3, ETA + 3))9 ALLOCATE (Y(CSI + 3, ETA + 3))

10 ALLOCATE (U(CSI + 2, ETA + 2))11 ALLOCATE (V(CSI + 2, ETA + 2))12 ALLOCATE (XM(CSI + 2, ETA + 2))13 ALLOCATE (YM(CSI + 2, ETA + 2))14 ALLOCATE (X_CSI(CSI + 2, ETA + 2))15 ALLOCATE (X_ETA(CSI + 2, ETA + 2))16 ALLOCATE (Y_CSI(CSI + 2, ETA + 2))17 ALLOCATE (Y_ETA(CSI + 2, ETA + 2))18 ALLOCATE (JACOB(CSI + 2, ETA + 2))19 ALLOCATE (C(CSI + 2, ETA + 2, TAU + 1))2021 ! LEITURA DOS PONTOS DA MALHA22 CALL READ_GRID (FNAME (2) , CSI , ETA , X, Y)2324 ! LEITURA DO CAMPO DE VELOCIDADE25 CALL READ_VEL_FIELD (FNAME (3) , CSI , ETA , U, V)2627 ! GERACAO DA MALHA FANTASMA28 CALL GHOST_MESH_GENERATION (CSI , ETA , X, Y)2930 ! INICIALIZACAO DO CAMPO DE VELOCIDADE NA MALHA FANTASMA31 CALL INIT_VEL_FIELD (CSI , ETA , U, V)3233 ! ATRIBUICAO DAS CONDICOES INICIAIS34 CALL INITIAL_CONDITIONS (CSI , ETA , TAU , C)3536 ! ATRIBUICAO DAS CONDICOES DE FRONTEIRA37 CALL BOUNDARY_CONDITIONS (CSI , ETA , TAU , C_PR , C)3839 ! GERACAO DA MALHA DESLOCADA40 CALL DISPLACED_MESH_GENERATION (CSI , ETA , X, Y, XM , YM)4142 ! INICIALIZACAO DA MATRIZ DE JACOBIANOS43 CALL INIT_JACOB_MATRIX (CSI , ETA , X, Y, X_CSI , X_ETA , Y_CSI , Y_ETA , JACOB)4445 ! RESOLUCAO NUMERICA DA EQUACAO DE TRANSPORTE46 DELTA_T = TF / REAL(TAU)47 CALL SOLVER (TF , KX , KY , DELTA_T , KAPPA , TAU , CSI , ETA , C, X_CSI , X_ETA , Y_CSI ,

Y_ETA , JACOB , U, V)4849 ! PREPARACAO DA ESCRITA DAS SAIDAS50 CALL PRE_WRITE ( FOLDER )5152 ! ESCRITA DOS PONTOS DA MALHA53 CALL WRITE_GRID (FOLDER , CSI , ETA , X, Y)5455 ! ESCRITA DO CAMPO DE VELOCIDADE

64

56 CALL WRITE_VEL_FIELD (FOLDER , DELTA_T , CSI , ETA , XM , YM , U, V)5758 ! ESCRITA DAS CONCENTRACOES A CADA NIVEL DE TEMPO59 CALL WRITE_CONC (TAU , CSI , ETA , FOLDER , DELTA_T , X, Y, C)6061 ! DESALOCACAO DE MEMORIA62 DEALLOCATE (X)63 DEALLOCATE (Y)64 DEALLOCATE (U)65 DEALLOCATE (V)66 DEALLOCATE (XM)67 DEALLOCATE (YM)68 DEALLOCATE (X_CSI)69 DEALLOCATE (X_ETA)70 DEALLOCATE (Y_CSI)71 DEALLOCATE (Y_ETA)72 DEALLOCATE (JACOB)73 DEALLOCATE (C)

Algoritmo 4.1 – Modelo de dispersão de poluentes na atmosfera

No Algoritmo 4.1, a subrotina READ_FNAMES obtém os nomes de três arquivospassados via linha de comando: um arquivo de informações gerais, um arquivo contendoos pontos da malha e um arquivo contendo o campo de velocidade para a malha analisada.

Em um primeiro momento, leem-se as informações gerais na subrotina READ_INFO,as quais se constituem: o coeficiente de decaimento (KAPPA), os coeficientes de difusão (KXe KY), uma concentração prescrita (C_PR), o número de partições na direção 𝜉 (CSI), 𝜂

(ETA) e 𝜏 (TAU), além de um tempo final (TF). A concentração prescrita é utilizada noAlgoritmo 4.5 e indica a concentração do poluente que estaria sendo lançada na região. Onúmero de partições representam a quantidade de espaçamentos nas respectivas direções.O tempo final determina qual o instante em que se finaliza o lançamento de poluente naregião. As demais informações fazem parte da equação de transporte.

Em seguida, é necessário alocar memória para as matrizes utilizadas no algoritmo,tendo em vista que já se tem conhecimento do número de partições em 𝜉, 𝜂 e 𝜏 . Alocadasas estruturas de dados, leem-se os pontos da malha na subrotina READ_GRID, especifica-das pelas localizações em 𝑥 e 𝑦, atribuindo-os às variáveis X e Y. Lê-se também o campode velocidade na subrotina READ_VEL_FIELD proveniente dos resultados da Equação deNavier-Stokes, dada a malha previamente lida, atribuindo-o às variáveis U e V. É impor-tante destacar que o campo de velocidade lido se encontra no seu regime permanente,sendo empregado em todos os passos de tempo da resolução numérica da equação detransporte.

No Algoritmo 4.2 da subrotina GHOST_MESH_GENERATION, gera-se uma malha fan-tasma, a qual é composta por um conjunto de células-fantasma que envolvem a malhaoriginal, com o objetivo de atribuir as condições de contorno do problema. Dado um ponto𝐴 adjacente a um ponto 𝐵 pertencente a uma das fronteiras da malha original, o ponto

65

𝐶 pertencente à malha fantasma é encontrado pela relação −→𝐴𝐵 = −−→

𝐵𝐶. Esse cálculo podeser melhor visualizado pela Figura 9 e é realizado para todas as quatro fronteiras da ge-ometria. Na Figura 9, as células em vermelho constituem a malha fantasma, e as célulasem preto se referem à malha original. Ainda, pela Figura 9, pode-se perceber que, nospontos que interligam duas fronteiras, atribui-se o próprio ponto da malha original queinterliga as mesmas fronteiras.

1 DO I = 2, CSI + 22 X(I, 1) = 2 * X(I, 2) - X(I, 3)3 Y(I, 1) = 2 * Y(I, 2) - Y(I, 3)45 X(I, ETA + 3) = 2 * X(I, ETA + 2) - X(I, ETA + 1)6 Y(I, ETA + 3) = 2 * Y(I, ETA + 2) - Y(I, ETA + 1)7 END DO89 DO J = 2, ETA + 2

10 X(1, J) = 2 * X(2, J) - X(3, J)11 Y(1, J) = 2 * Y(2, J) - Y(3, J)1213 X(CSI + 3, J) = 2 * X(CSI + 2, J) - X(CSI + 1, J)14 Y(CSI + 3, J) = 2 * Y(CSI + 2, J) - Y(CSI + 1, J)15 END DO1617 X(1, 1) = X(2, 2)18 Y(1, 1) = Y(2, 2)1920 X(1, ETA + 3) = X(2, ETA + 2)21 Y(1, ETA + 3) = Y(2, ETA + 2)2223 X(CSI + 3, 1) = X(CSI + 2, 2)24 Y(CSI + 3, 1) = Y(CSI + 2, 2)2526 X(CSI + 3, ETA + 3) = X(CSI + 2, ETA + 2)27 Y(CSI + 3, ETA + 3) = Y(CSI + 2, ETA + 2)

Algoritmo 4.2 – Geração da malha fantasma

Figura 9 – Malha fantasma

Assim, com a malha fantasma gerada, é preciso inicializar o campo de velocidadenessa malha, uma vez que seus valores estão definidos apenas para a malha original. NoAlgoritmo 4.3 da subrotina INIT_VEL_FIELD, atribuem-se as velocidades U e V nas células-fantasma, a partir das respectivas velocidades das células adjacentes pertencentes à malhaoriginal. O valor da velocidade nos pontos que interligam duas fronteiras é nulo, uma vezque não interfere na solução do problema.

66

1 DO I = 2, CSI + 12 U(I, 1) = U(I, 2)3 V(I, 1) = V(I, 2)45 U(I, ETA + 2) = U(I, ETA + 1)6 V(I, ETA + 2) = V(I, ETA + 1)7 END DO89 DO J = 2, ETA + 1

10 U(1, J) = U(2, J)11 V(1, J) = V(2, J)1213 U(CSI + 2, J) = U(CSI + 1, J)14 V(CSI + 2, J) = V(CSI + 1, J)15 END DO

Algoritmo 4.3 – Inicialização do campo de velocidade na malha fantasma

Nos Algoritmos 4.4 e 4.5 das correspondentes subrotinas INITIAL_CONDITIONS eBOUNDARY_CONDITIONS, determinam-se, respectivamente, as condições iniciais e de fron-teira do problema. As condições iniciais são definidas nas células da malha original parao primeiro passo de tempo, e, no caso, atribuiu-se um valor nulo às concentrações C. OAlgoritmo 4.4 foi implementado apenas para ser facilmente alterado, caso as condiçõesiniciais do problema sejam estabelecidas por valores não nulos ou, até mesmo, funções.

1 DO I = 1, CSI + 22 DO J = 1, ETA + 23 C(I, J, 1) = 0.04 END DO5 END DO

Algoritmo 4.4 – Atribuição das condições iniciais

As condições de fronteira, por sua vez, são definidas nas células-fantasma geradaspara todos os passos de tempo. Foram empregados dois tipos de condições: de Dirichletpara a fronteira esquerda e de Neumann para as demais fronteiras. A condição de contornode Dirichlet, também denominada, no caso, como condição de injeção prescrita, consiste naatribuição de um valor pré-estabelecido às concentrações C da fronteira, dado pela variávelC_PR lida no arquivo de informações gerais. Já a condição de contorno de Neumann utilizaas respectivas concentrações das células adjacentes pertencentes à malha original para aatribuição. Assim como para o campo de velocidade, o valor da concentração nos pontosque interligam duas fronteiras é nulo.

67

1 DO K = 1, TAU + 12 DO I = 2, CSI + 13 C(I, 1, K) = C(I, 2, K)4 C(I, ETA + 2, K) = C(I, ETA + 1, K)5 END DO67 DO J = 2, ETA + 18 C(1, J, K) = C_PR9 C(CSI + 2, J, K) = C(CSI + 1, J, K)

10 END DO11 END DO

Algoritmo 4.5 – Atribuição das condições de fronteira

No Algoritmo 4.6 da subrotina DISPLACED_MESH_GENERATION, constrói-se umamalha deslocada, alocada nas matrizes XM e YM, na qual seus nós computacionais cor-respondem aos centros das células da malha original. Esse centro é calculado pela médiaaritmética dos quatro vértices da célula, e esse passo é realizado para todas as células damalha original. Na Figura 10, ilustra-se a malha original na cor preta e a malha deslocada,após construída, na cor vermelha. A construção dessa malha tem como finalidade apenasa escrita da saída do campo de velocidade, uma vez que os vetores de velocidade nasdireções 𝑥 e 𝑦 são calculados sobre a face da célula, mas representados no seu centro.

1 DO I = 1, CSI + 22 DO J = 1, ETA + 23 XM(I, J) = 0.25 * (X(I, J) + X(I + 1, J) + X(I, J + 1) + X(I + 1, J + 1))4 YM(I, J) = 0.25 * (Y(I, J) + Y(I + 1, J) + Y(I, J + 1) + Y(I + 1, J + 1))5 END DO6 END DO

Algoritmo 4.6 – Geração da malha deslocada

Figura 10 – Malha deslocada

Lembrando que o Jacobiano constitui uma variável bastante importante no con-texto do sistema de coordenadas generalizadas e que seu valor é dependente da malha com-putacional, calculam-se os Jacobianos no Algoritmo 4.7 da subrotina INIT_JACOB_MATRIXpreviamente, atribuindo-os à matriz JACOB. Conforme os Jacobianos são calculados, asmatrizes auxiliares X_CSI, X_ETA, Y_CSI e Y_ETA também são construídas, indicando res-pectivamente as derivadas parciais 𝑥𝜉, 𝑥𝜂, 𝑦𝜉 e 𝑦𝜂. Esse passo é realizado antes da resoluçãonumérica da equação de transporte, a fim de evitar processamentos desnecessários e re-duzir o tempo computacional.

68

1 DO I = 1, CSI + 22 DO J = 1, ETA + 23 X_CSI(I, J) = 0.5 * (X(I + 1, J + 1) + X(I + 1, J) - X(I, J + 1) - X(I, J)

)4 X_ETA(I, J) = 0.5 * (X(I + 1, J + 1) + X(I, J + 1) - X(I + 1, J) - X(I, J)

)5 Y_CSI(I, J) = 0.5 * (Y(I + 1, J + 1) + Y(I + 1, J) - Y(I, J + 1) - Y(I, J)

)6 Y_ETA(I, J) = 0.5 * (Y(I + 1, J + 1) + Y(I, J + 1) - Y(I + 1, J) - Y(I, J)

)78 JACOB(I, J) = 1.0 / (X_CSI(I, J) * Y_ETA(I, J) - X_ETA(I, J) * Y_CSI(I, J)

)9 END DO

10 END DO

Algoritmo 4.7 – Inicialização da matriz de Jacobianos

No programa principal, mostrado no Algoritmo 4.1, calcula-se o DELTA_T, referindo-se à variação no tempo. Posteriormente, no Algoritmo 4.8 da subrotina SOLVER, apresenta-se a resolução do modelo numérico da equação de transporte. Esse algoritmo constitui-seno módulo fundamental de toda a implementação. Para cada nível de tempo, calcula-sea concentração do poluente no próximo nível para todas as células da malha original.

Em um primeiro momento, verifica-se se os coeficientes de difusão KX e KY sãonulos ou não, para evitar cálculos desnecessários. Caso os valores de ambos sejam nulos,atribui-se 0.0 à variável DIFFUSION, correspondente à D(𝐶)|𝑘𝑃 no modelo numérico. Casocontrário, é preciso calcular os valores de 𝛿𝑖, atribuídos à variável D_SUBPART, onde 𝑖 =1, . . . , 8. Em seguida, calculam-se os valores de D𝑗(𝐶)|𝑘𝑃 , atribuídos à variável D_PART, emque 𝑗 = 1, . . . , 4. Finalmente, calcula-se o valor de D(𝐶)|𝑘𝑃 e atribui-se à DIFFUSION. Valelembrar que, dependendo da localização na malha no processo iterativo, diferentes técnicasde discretização em diferenças finitas são empregadas para o cálculo de D_SUBPART.

A próxima etapa consiste em determinar o valor de A (𝐶)|𝑘𝑃 do modelo numérico,o qual é atribuído à variável ADVECTION. Para tal, calculam-se, em um primeiro instante,os valores de 𝐶|𝑘𝑒 , 𝐶|𝑘𝑤, 𝐶|𝑘𝑛 e 𝐶|𝑘𝑠 , as quais realizam chamadas às funções upwind em suasoperações, descritas no Algoritmo 4.9, e são atribuídas à variável A_PART.

Com os valores de DIFFUSION e ADVECTION obtidos, assim como outras informaçõesdefinidas anteriormente, é possível encontrar o valor da concentração do poluente nalocalização da malha atual durante o processo iterativo para o próximo nível de tempo.

69

1 DO K = 1, TAU2 DO I = 2, CSI + 13 DO J = 2, ETA + 14 IF (KX /= 0.0 .OR. KY /= 0.0) THEN5 IF (I == 2) THEN6 D_SUBPART (1) = (C(I + 2, J, K) - C(I, J, K)) * 0.57 D_SUBPART (2) = (C(I + 1, J + 1, K) - C(I + 1, J - 1, K)) * 0.58 D_SUBPART (3) = C(I, J, K) - C(I - 1, J, K)9 IF (J == 2) THEN

10 D_SUBPART (4) = C(I - 1, J + 1, K) - C(I - 1, J, K)11 ELSE IF (J == ETA + 1) THEN12 D_SUBPART (4) = C(I - 1, J, K) - C(I - 1, J - 1, K)13 ELSE14 D_SUBPART (4) = (C(I - 1, J + 1, K) - C(I - 1, J - 1, K)) * 0.515 END IF16 ELSE IF (I == CSI + 1) THEN17 D_SUBPART (1) = (C(I + 1, J, K) - C(I, J, K))18 IF (J == 2) THEN19 D_SUBPART (2) = C(I + 1, J + 1, K) - C(I + 1, J, K)20 ELSE IF (J == ETA + 1) THEN21 D_SUBPART (2) = C(I + 1, J, K) - C(I + 1, J - 1, K)22 ELSE23 D_SUBPART (2) = (C(I + 1, J + 1, K) - C(I + 1, J - 1, K)) * 0.524 END IF25 D_SUBPART (3) = (C(I, J, K) - C(I - 2, J, K)) * 0.526 D_SUBPART (4) = (C(I - 1, J + 1, K) - C(I - 1, J - 1, K)) * 0.527 ELSE28 D_SUBPART (1) = (C(I + 2, J, K) - C(I, J, K)) * 0.529 D_SUBPART (2) = (C(I + 1, J + 1, K) - C(I + 1, J - 1, K)) * 0.530 D_SUBPART (3) = (C(I, J, K) - C(I - 2, J, K)) * 0.531 D_SUBPART (4) = (C(I - 1, J + 1, K) - C(I - 1, J - 1, K)) * 0.532 END IF3334 IF (J == 2) THEN35 D_SUBPART (5) = (C(I, J + 2, K) - C(I, J, K)) * 0.536 D_SUBPART (6) = (C(I + 1, J + 1, K) - C(I - 1, J + 1, K)) * 0.537 D_SUBPART (7) = C(I, J, K) - C(I, J - 1, K)38 IF (I == 2) THEN39 D_SUBPART (8) = C(I + 1, J - 1, K) - C(I, J - 1, K)40 ELSE IF (I == CSI + 1) THEN41 D_SUBPART (8) = C(I, J - 1, K) - C(I - 1, J - 1, K)42 ELSE43 D_SUBPART (8) = (C(I + 1, J - 1, K) - C(I - 1, J - 1, K)) * 0.544 END IF45 ELSE IF (J == ETA + 1) THEN46 D_SUBPART (5) = C(I, J + 1, K) - C(I, J, K)47 IF (I == 2) THEN48 D_SUBPART (6) = C(I + 1, J + 1, K) - C(I, J + 1, K)49 ELSE IF (I == CSI + 1) THEN50 D_SUBPART (6) = C(I, J + 1, K) - C(I - 1, J + 1, K)51 ELSE52 D_SUBPART (6) = (C(I + 1, J + 1, K) - C(I - 1, J + 1, K)) * 0.553 END IF54 D_SUBPART (7) = (C(I, J, K) - C(I, J - 2, K)) * 0.555 D_SUBPART (8) = (C(I + 1, J - 1, K) - C(I - 1, J - 1, K)) * 0.556 ELSE

70

57 D_SUBPART (5) = (C(I, J + 2, K) - C(I, J, K)) * 0.558 D_SUBPART (6) = (C(I + 1, J + 1, K) - C(I - 1, J + 1, K)) * 0.559 D_SUBPART (7) = (C(I, J, K) - C(I, J - 2, K)) * 0.560 D_SUBPART (8) = (C(I + 1, J - 1, K) - C(I - 1, J - 1, K)) * 0.561 END IF6263 D_PART (1) = 0.5 * (JACOB(I + 1, J) * (Y_ETA(I + 1, J) ** 2.0) *

D_SUBPART (1) &64 - JACOB(I + 1, J) * Y_ETA(I + 1, J) * Y_CSI(I + 1, J) * D_SUBPART

(2) &65 - JACOB(I - 1, J) * (Y_ETA(I - 1, J) ** 2.0) * D_SUBPART (3) &66 + JACOB(I - 1, J) * Y_ETA(I - 1, J) * Y_CSI(I - 1, J) * D_SUBPART

(4))67 D_PART (2) = 0.5 * (JACOB(I, J + 1) * (Y_CSI(I, J + 1) ** 2.0) *

D_SUBPART (5) &68 - JACOB(I, J + 1) * Y_ETA(I, J + 1) * Y_CSI(I, J + 1) * D_SUBPART

(6) &69 - JACOB(I, J - 1) * (Y_CSI(I, J - 1) ** 2.0) * D_SUBPART (7) &70 + JACOB(I, J - 1) * Y_ETA(I, J - 1) * Y_CSI(I, J - 1) * D_SUBPART

(8))71 D_PART (3) = 0.5 * (JACOB(I + 1, J) * (X_ETA(I + 1, J) ** 2.0) *

D_SUBPART (1) &72 - JACOB(I + 1, J) * X_ETA(I + 1, J) * X_CSI(I + 1, J) * D_SUBPART

(2) &73 - JACOB(I - 1, J) * (X_ETA(I - 1, J) ** 2.0) * D_SUBPART (3) &74 + JACOB(I - 1, J) * X_ETA(I - 1, J) * X_CSI(I - 1, J) * D_SUBPART

(4))75 D_PART (4) = 0.5 * (JACOB(I, J + 1) * (X_CSI(I, J + 1) ** 2.0) *

D_SUBPART (5) &76 - JACOB(I, J + 1) * X_ETA(I, J + 1) * X_CSI(I, J + 1) * D_SUBPART

(6) &77 - JACOB(I, J - 1) * (X_CSI(I, J - 1) ** 2.0) * D_SUBPART (7) &78 + JACOB(I, J - 1) * X_ETA(I, J - 1) * X_CSI(I, J - 1) * D_SUBPART

(8))7980 DIFFUSION = KX * ( D_PART (1) + D_PART (2)) + KY * ( D_PART (3) + D_PART

(4))81 ELSE82 DIFFUSION = 0.083 END IF8485 A_PART (1) = (1 + UPWIND_U (I, J)) * C(I, J, K) * 0.5 + (1 - UPWIND_U (I,

J)) * C(I + 1, J, K) * 0.586 A_PART (2) = (1 + UPWIND_U (I - 1, J)) * C(I - 1, J, K) * 0.5 + (1 -

UPWIND_U (I - 1, J)) * C(I, J, K) * 0.587 A_PART (3) = (1 + UPWIND_V (I, J)) * C(I, J, K) * 0.5 + (1 - UPWIND_V (I,

J)) * C(I, J + 1, K) * 0.588 A_PART (4) = (1 + UPWIND_V (I, J - 1)) * C(I, J - 1, K) * 0.5 + (1 -

UPWIND_V (I, J - 1)) * C(I, J, K) * 0.58990 ADVECTION = U(I, J) * A_PART (1) - U(I - 1, J) * A_PART (2) + V(I, J) *

A_PART (3) - V(I, J - 1) * A_PART (4)9192 C(I, J, K + 1) = DELTA_T * (-KAPPA * C(I, J, K) + JACOB(I, J) * (

DIFFUSION - ADVECTION )) + C(I, J, K)93 END DO94 END DO

71

95 END DO

Algoritmo 4.8 – Resolução numérica da equação de transporte

O Algoritmo 4.9 demonstra duas funções upwind: uma para o fluxo de vento nadireção 𝑥 (UPWIND_U) e outro para a direção 𝑦 (UPWIND_V). Essas funções também estãocontidas na subrotina SOLVER e retornam -1.0 se o sentido do fluxo de vento for negativo;ou 1.0, caso contrário.

1 REAL FUNCTION UPWIND_U (I, J)23 IMPLICIT NONE4 INTEGER , INTENT (IN) :: I, J56 IF (U(I, J) >= 0.0) THEN7 UPWIND_U = 1.08 ELSE9 UPWIND_U = -1.0

10 END IF1112 END FUNCTION UPWIND_U1314 REAL FUNCTION UPWIND_V (I, J)1516 IMPLICIT NONE17 INTEGER , INTENT (IN) :: I, J1819 IF (V(I, J) >= 0.0) THEN20 UPWIND_V = 1.021 ELSE22 UPWIND_V = -1.023 END IF2425 END FUNCTION UPWIND_V

Algoritmo 4.9 – Funções upwind

Resolvida numericamente a equação de transporte, o Algoritmo 4.1 se direci-ona à escrita dos resultados. Uma pequena subrotina foi desenvolvida, nomeada comoPRE_WRITE, para a remoção e recriação de uma pasta para o armazenamento das saídas,caso ela já exista; caso contrário, apenas cria-se a nova pasta. O nome dessa pasta estácontido em FOLDER e é determinada como parâmetro em sua declaração, se constituindouma constante no programa.

As subrotinas WRITE_GRID, WRITE_VEL_FIELD e WRITE_CONC são então executa-das, criando, respectivamente, um arquivo contendo a malha computacional, um arquivocontendo o campo de velocidade e vários arquivos contendo as concentrações do poluente,um para cada nível de tempo. É importante levar em consideração que todos os arquivosgerados na saída do algoritmo estão no formato .vtk, tendo em vista a visualização gráficados resultados no software empregado. Ressalta-se também que os arquivos contendo os

72

pontos da malha e o campo de velocidade na entrada do algoritmo possuem formatosdiferentes da saída, e, por esse motivo, é necessário a criação de novos arquivos.

Assim, o Algoritmo 4.1 realiza a desalocação de memória das matrizes alocadase finaliza sua execução. No Algoritmo 4.10, identifica-se o arquivo Make, contendo asinstruções para a compilação do código implementado em Fortran 90 e gerando o arquivoexecutável EXEC.

1 GF = gfortran2 FLAGS = -o3 FILES = EQ_TRANSPORTE .f90 READ_FNAMES .f90 READ_INFO .f90 READ_GRID .f90

READ_VEL_FIELD .f90 \4 GHOST_MESH_GENERATION .f90 INIT_VEL_FIELD .f90 INITIAL_CONDITIONS .f90

BOUNDARY_CONDITIONS .f90 DISPLACED_MESH_GENERATION .f90 \5 INIT_JACOB_MATRIX .f90 SOLVER .f90 PRE_WRITE .f90 WRITE_GRID .f90 WRITE_VEL_FIELD

.f90 \6 WRITE_CONC .f9078 PROG: $(FILES)9 $(GF) $(FLAGS) EXEC $(FILES)

Algoritmo 4.10 – Arquivo Make

Apontados os procedimentos metodológicos para a realização deste trabalho, noCapítulo 5, mostram-se os experimentos executados e discutem-se os resultados encontra-dos.

73

5 RESULTADOS NUMÉRICOS E DISCUSSÃO

Neste capítulo, apresentam-se os experimentos efetuados, assim como os resultadosobtidos. Três testes são realizados: um duto simples, um duto com um obstáculo e umduto com dois obstáculos. Para o teste do duto simples, quatro casos são tratados, coma finalidade de mostrar os processos atrelados à dispersão de poluentes. Além disso, paraa visualização gráfica de todos os testes executados, o software Paraview1 é empregado,cujos arquivos de entrada estão no formato .vtk, como relatado na Seção 4.2.

A Tabela 1 contém as principais informações necessárias para a geração da ma-lha computacional: o número de partições nas direções 𝜉 e 𝜂. As localizações dos pontospertencentes às fronteiras são definidas manualmente, tendo, em todas as malhas, o com-primento máximo de 20 metros e a altura máxima de 3 metros.

Tabela 1 – Parâmetros para a geração da malha computacional

Duto Partições em 𝜉 Partições em 𝜂Simples 100 15Com um obstáculo 110 18Com dois obstáculos 120 18

De acordo com a Tabela 1, é possível calcular a quantidade de células da malha,multiplicando-se os dois números de partições. As quantidades de partições em ambas asdireções foram estipuladas visualmente, buscando uma malha refinada e, em alguns casos,a concentração de linhas em uma determinada localização da geometria para uma melhoranálise dos resultados na região.

A Tabela 2 informa os parâmetros necessários para o modelo numérico da Equaçãode Navier-Stokes: o tempo final 𝑡𝑓 (em s), o número de partições na direção 𝜏 , a massaespecífica 𝜌 (em kg/m3), a viscosidade dinâmica 𝜇 (em N.s/m2) e a velocidade de injeçãoprescrita 𝑣𝑖 (em m/s).

Tabela 2 – Parâmetros para a resolução numérica da Equação de Navier-Stokes

Duto 𝑡𝑓 Partições em 𝜏 𝜌 𝜇 𝑣𝑖

Simples 100.0 100000 1.293 0.029093 0.75Com um obstáculo 100.0 100000 1.293 0.11637 0.3Com dois obstáculos 100.0 100000 1.293 0.11637 0.3

O valor da densidade 𝜌 do ar é definido com base no livro de Taylor [42], nascondições normais de temperatura e pressão (CNTP), em que a temperatura do ambi-1 <http://www.paraview.org/>

74

ente corresponde a 0°C e a pressão equivale a 101325 Pa. A velocidade de injeção 𝑣𝑖 éencontrada empiricamente, pois é influenciada pela geometria do problema. Finalmente,a viscosidade dinâmica 𝜇 é calculada através da Equação (5.1) [43].

𝑅𝑒 = 𝜌𝐿𝑣𝑖

𝜇, (5.1)

em que 𝑅𝑒 é o número de Reynolds2, e 𝐿, o comprimento característico, o qual, no caso,corresponde a 3 metros. Como na dissertação de Barba [2], a resolução da Equação deNavier-Stokes está direcionada ao tratamento de escoamentos laminares, o número deReynolds constitui-se um valor relativamente pequeno, em comparação aos valores deescoamentos turbulentos. Por essa razão, empregam-se, no experimento do duto simples,𝑅𝑒 = 100 e, nos demais experimentos, 𝑅𝑒 = 10.

A Tabela 3 indica os parâmetros fundamentais para o modelo numérico da equaçãode transporte: o tempo final 𝑡𝑓 (em s), o número de partições na direção 𝜏 , o coeficientede decaimento 𝜅, os coeficientes de difusão 𝐾𝑥 e 𝐾𝑦 nas direções 𝑥 e 𝑦 (em m2/s), respec-tivamente, e a concentração de injeção prescrita 𝐶𝑖 (em kg/m3).

Tabela 3 – Parâmetros para a resolução numérica da equação de transporte

Duto 𝑡𝑓 Partições em 𝜏 𝜅 𝐾𝑥 𝐾𝑦 𝐶𝑖

Simples (caso 1) 100.0 100000 0.0 0.0 0.0 1.0Simples (caso 2) 100.0 100000 0.5 0.0 0.0 1.0Simples (caso 3) 100.0 100000 0.0 0.5 0.0 1.0Simples (caso 4) 100.0 100000 0.0 0.0 0.5 1.0Com um obstáculo 10.0 100000 0.1 0.3 0.1 1.0Com dois obstáculos 10.0 100000 0.1 0.3 0.1 1.0

É importante salientar, na Tabela 3, que os valores de 𝜅, 𝐾𝑥, 𝐾𝑦 e 𝐶𝑖, nos trêsexperimentos, constituem apenas valores para teste, como forma de comprovação teóricados processos dispersivos de advecção, difusão e reação, associados à equação de trans-porte. O coeficiente de decaimento necessita ser calculado quimicamente, pois indica ataxa de perda de uma determinada concentração, dada a reação química que o poluenteé capaz de encadear com outra substância presente no ar atmosférico. Os coeficientes dedifusão 𝐾𝑥 e 𝐾𝑦 também necessitam de uma análise química, uma vez que são dependen-tes do meio em que o poluente se propaga. Finalmente, utiliza-se um valor unitário paraa concentração 𝐶𝑖, visto que seu valor se deriva da quantidade de poluente que está sendolançada sobre a geometria.

Por fim, os valores de tempo final 𝑡𝑓 e o número de partições na direção 𝜏 presentesnas Tabelas 2 e 3 são determinados de acordo com dois critérios: o instante que se deseja2 O número de Reynolds permite classificar os escoamentos de fluidos em laminares ou turbulentos,

dependendo do seu valor [43].

75

finalizar a execução do algoritmo de resolução e o valor de Δ𝜏 , calculado pela divisão entreo tempo final e o número de partições. Muitas vezes, existe a necessidade de estipular umΔ𝜏 bastante pequeno, dependendo do problema em questão.

Nas Seções 5.1, 5.2 e 5.3, apresentam-se os três experimentos executados separada-mente, demonstrando as características da malha computacional, do campo de velocidadee das concentrações de um poluente qualquer lançado sobre a região.

5.1 Duto Simples

Como foi mencionado na introdução deste capítulo, para o experimento envolvendoum duto simples como geometria, foram tratados quatro casos, variando-se as configu-rações de 𝜅, 𝐾𝑥 e 𝐾𝑦. Para todos os casos, utilizam-se as mesmas definições da malhacomputacional e do campo de velocidade. Esse experimento foi realizado com a finali-dade de possibilitar a visualização da influência dos processos dispersivos no campo deconcentração.

A Figura 11 ilustra a malha computacional do duto simples, possuindo um formatoretangular, cujas células possuem dimensões iguais de comprimento e altura.

Figura 11 – Malha do duto simples

A Figura 12 representa o campo de velocidade, em m/s, do fluxo de vento para amalha do duto simples, na qual os maiores valores de velocidade se concentraram na regiãocentral do duto, como se pode identificar na Figura 12a. As Figuras 12b e 12c indicam adireção e o sentido do vetor velocidade, dado respectivamente sobre a geometria completado duto e sobre apenas um setor da geometria.

76

(a) Módulo do vetor velocidade (b) Direção e sentido do vetor velocidade

(c) Setor do campo de velocidade

Figura 12 – Campo de velocidade do duto simples

No primeiro caso, onde os valores 𝜅, 𝐾𝑥 e 𝐾𝑦 correspondem a 0.0, obteve-se ocampo de concentração, em kg/m3, do poluente no duto simples mostrado na Figura 13.Esse teste tem como objetivo demonstrar o processo advectivo da equação de transporte,uma vez que, ao se anularem os valores de 𝜅, 𝐾𝑥 e 𝐾𝑦, os processos reativo e difusivo sãoeliminados.

(a) 𝜏 = 0.0 (b) 𝜏 = 0.2

(c) 𝜏 = 0.4 (d) 𝜏 = 0.6

(e) 𝜏 = 0.8 (f) 𝜏 = 1.0

Figura 13 – Concentrações do poluente no duto simples (caso 1)

Tendo em vista que o poluente é injetado na fronteira esquerda da geometria,percebe-se, na Figura 13, que, com 800 iterações no tempo (𝜏 = 0.8), o poluente jáalcança a fronteira direita e que, com 1000 iterações (𝜏 = 1.0), ocupa uma parte bastante

77

significante do duto. Porém, o mais importante a se notar constitui-se no fato de queo poluente se desloca conforme o campo de velocidade. Nos locais onde a velocidade émaior, o poluente se desloca de maneira mais rápida; analogamente, nas localizações ondea velocidade é menor, o poluente se desloca vagarosamente.

No segundo caso, onde apenas o valor do coeficiente de decaimento 𝜅 é não nulo,obteve-se o campo de concentração, em kg/m3, do poluente no duto simples descritona Figura 14. Esse teste objetivou a demonstração do processo reativo da equação detransporte, quando se comparam os resultados encontrados com aqueles apresentados noprimeiro caso.

(a) 𝜏 = 0.0 (b) 𝜏 = 0.2

(c) 𝜏 = 0.4 (d) 𝜏 = 0.6

(e) 𝜏 = 0.8 (f) 𝜏 = 1.0

Figura 14 – Concentrações do poluente no duto simples (caso 2)

Comparativamente à Figura 13, pode-se verificar que, com 𝜅 = 0.5, houve grandeperda de concentração do poluente, visto que, para uma mesma partição 𝜏 , o campodeixou de ter uma coloração avermelhada e passou a ter uma coloração alaranjada ouamarelada, significando a diminuição dos seus valores. Como esperado, o coeficiente dedecaimento aplicado teve grande relevância nessa diminuição.

No terceiro caso, onde apenas o coeficiente de difusão 𝐾𝑥 é não nulo, obteve-se o campo de concentração, em kg/m3, do poluente no duto simples representado naFigura 15. Esse teste procurou demonstrar o processo difusivo da equação de transportena direção 𝑥, ao se compararem os resultados encontrados com aqueles ilustrados noprimeiro caso, dado na Figura 13.

78

(a) 𝜏 = 0.0 (b) 𝜏 = 0.2

(c) 𝜏 = 0.4 (d) 𝜏 = 0.6

(e) 𝜏 = 0.8 (f) 𝜏 = 1.0

Figura 15 – Concentrações do poluente no duto simples (caso 3)

Em relação à Figura 13, pode-se notar que, com 𝐾𝑥 = 0.5, houve uma influênciapouco perceptível nos campos de concentrações apresentados. O motivo para tal aconte-cimento foi devido ao fato de que o campo de velocidade está mais direcionado no sentidodo eixo 𝑥, assim como o coeficiente de difusão empregado, fazendo que o poluente seespalhe na mesma direção do fluxo de vento e impossibilitando uma melhor visualizaçãodo funcionamento do processo.

No entanto, no quarto caso, apenas o coeficiente de difusão 𝐾𝑦 é não nulo, obtendoum campo de concentração, em kg/m3, do poluente no duto simples apresentado na Figura16. Assim como o terceiro caso, objetivou-se demonstrar o processo difusivo da equaçãode transporte pela comparação com os resultados do primeiro caso, porém na direção 𝑦.

79

(a) 𝜏 = 0.0 (b) 𝜏 = 0.2

(c) 𝜏 = 0.4 (d) 𝜏 = 0.6

(e) 𝜏 = 0.8 (f) 𝜏 = 1.0

Figura 16 – Concentrações do poluente no duto simples (caso 4)

Em comparação às Figuras 13 e 15, pode-se observar que, com 𝐾𝑦 = 0.5, os camposde concentrações sofreram alterações relevantes no processo de dispersão. Pelo fato de queo poluente, desta vez, está se espalhando no sentido do eixo 𝑦, perpendicular ao campo develocidade, percebe-se que houve uma redução considerável na concentração decorrente doprocesso difusivo, uma vez que o campo anteriormente estava representado na coloraçãoavermelhada e se alterou para uma coloração amarelada ou esverdeada.

Na Seção 5.2, apresenta-se o segundo experimento realizado, no qual um duto comum obstáculo é empregado.

5.2 Duto com um Obstáculo

No experimento do duto com um obstáculo, foi inserida uma barreira na malhado duto simples, na qual o comprimento possui 0.15 metros, e a altura, 1.25 metros.Esse experimento foi aplicado com o intuito de proporcionar a visualização da dispersãode poluentes em geometrias idênticas ou semelhantes. Ainda, essa geometria auxiliarácomo base no estudo do processo dispersivo para outra geometria, na qual dois obstáculosbastante próximos um do outro formam a estrutura denominada como cânion urbano,comumente encontrada na literatura.

A Figura 17 ilustra a malha computacional do duto com um obstáculo. De acordocom a Tabela 1, percebe-se que se aumentou o número de partições na direção 𝜉, devidoà presença do obstáculo. Além disso, procurou-se manter um aglomerado de linhas nas

80

proximidades da barreira, com a finalidade de melhorar a captura das concentraçõesnessas localizações. Como consequência, nota-se, na Figura 17, que essa atração afetou asfronteiras esquerda e direita da geometria.

Figura 17 – Malha do duto com um obstáculo

A Figura 18 demonstra o campo de velocidade, em m/s, do fluxo de vento paraa malha do duto com um obstáculo, na qual os maiores valores de velocidade se concen-traram nas proximidades do obstáculo, após o fluido ultrapassá-lo, como se percebe naFigura 18a. Houve também maiores velocidades em parte da fronteira direita da região,resultante da configuração da malha computacional.

(a) Módulo do vetor velocidade (b) Direção e sentido do vetor velocidade

(c) Linhas de corrente do campo de velocidade

(d) Setor do campo de velocidade

Figura 18 – Campo de velocidade do duto com um obstáculo

Ainda, na Figura 18a, é importante levar em consideração a formação de umvórtice à direita do obstáculo, representada na coloração azul clara. Nas Figuras 18b e18d, nas quais se descrevem a direção e o sentido do vetor velocidade, e, na Figura 18c,na qual são mostradas as linhas de corrente do campo de velocidade, nota-se, com maiorclareza, a formação do vórtice. O número de Reynolds, calculado pelos valores de 𝜌, 𝜇 e𝑣𝑖 estipulados, além da geometria do problema, está diretamente relacionado à formaçãodesse efeito [43].

81

Nesse experimento, os valores de 𝜅, 𝐾𝑥 e 𝐾𝑦 são não nulos, significando a existênciade todos os processos envolvidos na dispersão do poluente. A Figura 19 representa o campode concentração, em kg/m3, de um poluente hipotético conforme o tempo, dada sua malhacomputacional e seu campo de velocidade.

(a) 𝜏 = 0.0 (b) 𝜏 = 0.25

(c) 𝜏 = 0.5 (d) 𝜏 = 0.75

(e) 𝜏 = 1.0 (f) 𝜏 = 1.5

Figura 19 – Concentrações do poluente no duto com um obstáculo

Inicialmente, é interessante verificar, na Figura 19, que o poluente alcança o obstá-culo pelas proximidades da fronteira inferior da geometria. A razão pela qual isso aconte-ceu, foi devido ao fato de que o campo de velocidade nessa região é levemente maior, com-parado à região superior. Além disso, após o poluente ultrapassar o obstáculo, percebe-seuma maior concentração do poluente no vórtice formado, também causada pelos valoresde velocidade da região. Comparado ao experimento do duto simples, esse experimentoexigiu uma quantidade maior de iterações para o poluente alcançar a fronteira direita,por conta de uma variação temporal Δ𝜏 menor e consistente ao problema, a fim de evitarinstabilidades numéricas.

Na Seção 5.3, realiza-se um último experimento e mostram-se os resultados encon-trados. No experimento, um duto com dois obstáculos de formatos diferentes são empre-gados.

5.3 Duto com dois Obstáculos

No experimento do duto com dois obstáculos, foram inseridas duas barreiras namalha do duto simples: uma, em um formato aproximado de uma chaminé, e outra, como

82

forma de ilustração de uma montanha. Assim como no duto com apenas um obstáculo,esse experimento foi executado para possibilitar a visualização dos processos dispersivos deum poluente em geometrias com características idênticas ou semelhantes. Essa geometriaservirá como base para um estudo de caso, no qual o primeiro obstáculo se constituirá afonte de emissão do poluente, e o segundo desempenhará o papel de uma barreira. Entrea fonte e a barreira, localizar-se-á uma cidade.

A Figura 20 ilustra a malha computacional do duto com dois obstáculos. O pri-meiro possui a altura de 0.75 metros, o comprimento da base de 0.2 metros e o compri-mento do topo de 0.15 metros. O segundo possui a altura de 1.2 metros e o comprimentoda base de 1.1 metros.

Figura 20 – Malha do duto com dois obstáculos

Nota-se, na Figura 20, que se aumentou o número de partições na direção 𝜉, coma finalidade de se manter um conjunto de linhas nas adjacências dos dois obstáculos parauma melhor aquisição dos resultados. Da mesma maneira que no duto com um obstáculo,a atração provocou mudanças nas linhas 𝜂 próximas às fronteiras direita e esquerda daregião.

A Figura 21 descreve o campo de velocidade, em m/s, do fluxo de vento para amalha do duto com dois obstáculos, onde as maiores velocidades se concentraram nasproximidades dos obstáculos, como se observa na Figura 21a. Ainda, foram produzidasmaiores velocidades em parte da vizinhança das fronteiras esquerda e direita.

83

(a) Módulo do vetor velocidade (b) Direção e sentido do vetor velocidade

(c) Linhas de corrente do campo de velocidade

(d) Setor do campo de velocidade

Figura 21 – Campo de velocidade do duto com dois obstáculos

Apesar de pouco perceptível, pode-se reparar, na Figura 21a, uma leve alteração nacoloração, de azul escuro para azul claro, entre os dois obstáculos e à direita do segundoobstáculo, significando a formação de dois vórtices na região. As Figuras 21b e 21d,descrevendo a direção e o sentido do vetor velocidade, e a Figura 21c, mostrando as linhasde corrente do campo de velocidade, possibilitam uma melhor visualização do fenômeno.

Utilizando os mesmos valores para os parâmetros 𝜅, 𝐾𝑥 e 𝐾𝑦 do experimentoanterior, esse experimento obteve como resultado o campo de concentração, em kg/m3,demonstrado na Figura 22, cuja solução também se assemelha bastante ao campo deconcentração para o duto com um obstáculo.

84

(a) 𝜏 = 0.0 (b) 𝜏 = 0.25

(c) 𝜏 = 0.5 (d) 𝜏 = 0.75

(e) 𝜏 = 1.0 (f) 𝜏 = 1.5

Figura 22 – Concentrações do poluente no duto com dois obstáculos

Tal como o experimento anterior, o poluente alcança inicialmente o primeiro obs-táculo pelas proximidades da fronteira inferior da geometria. Além disso, percebe-se tam-bém uma maior concentração do poluente nas duas zonas de recirculação de ar formadas,quando o poluente ultrapassa cada uma das barreiras. Dessa forma, pode-se constatarque, tanto neste experimento quanto no anterior, a presença de um vórtice no campo develocidade exerce forte influência sobre o campo de concentração, com a tendência de ele-var seus valores nessa região. No estudo de caso que se pretende analisar, a existência deum vórtice entre a fonte e a barreira poderia, portanto, tornar a cidade bastante poluída,ocasionando graves consequências à população.

Para uma melhor análise do algoritmo de resolução numérica da equação de trans-porte, efetua-se um estudo acerca dos tempos de processamento na Seção 5.4, de acordocom os experimentos realizados.

5.4 Análise dos Tempos de Processamento

Realizando uma pequena modificação no Algoritmo 4.1, calculam-se os temposde processamento para os três experimentos executados, lembrando que o primeiro ésubdividido em quatro casos. Esse cálculo é realizado pela diferença entre dois instantesobtidos, um antes da chamada (CALL) à subrotina GHOST_MESH_GENERATION e outro apósà chamada à subrotina SOLVER, abrangendo somente os módulos principais da aplicação.

A Tabela 4 apresenta os resultados encontrados para os tempos de processamento,

85

em segundos, com base no tipo de duto empregado.

Tabela 4 – Tempos de processamento dos experimentos realizados

Duto Tempo de processamentoSimples (caso 1) 13.6590004Simples (caso 2) 13.9740000Simples (caso 3) 57.3759995Simples (caso 4) 57.7270012Com um obstáculo 71.3229980Com dois obstáculos 77.7490005

As configurações do computador utilizado para a realização dos experimentos sãomostradas a seguir:

• Processador: Intel® Core™ i5-3317U 1.70GHz × 4;

• Memória RAM: 3.8GB;

• Sistema Operacional: Ubuntu 14.04 LTS 64-bit.

De acordo com as Tabelas 3 e 4, pode-se notar, em um primeiro momento, que oscasos 1 e 2 para o experimento do duto simples obtiveram um tempo de processamentorelativamente baixo, em comparação com os demais experimentos, por não possuir oprocesso difusivo da equação governante.

Quando o processo de difusão foi adicionado, através da determinação de valoresde 𝐾𝑥 ou 𝐾𝑦 não nulos, o tempo de processamento foi bastante elevado, devido à grandequantidade de condicionais e cálculos vinculados a esse processo, como se pode percebernos valores dos casos 3 e 4 do duto simples.

Além disso, quando a geometria se tornou mais complexa pela inclusão de obstá-culos, o aumento no número de partições na direção 𝜉 e 𝜂 foi responsável pela elevaçãodo tempo de processamento. Caso se aumentasse o número de partições na direção 𝜏 , aconsequência seria semelhante, aumentando o número de iterações para a resolução doproblema.

Discutidos os experimentos efetuados e os respectivos resultados, o Capítulo 6encerra este trabalho, apontando algumas considerações finais.

87

6 CONSIDERAÇÕES FINAIS

Este trabalho desenvolveu um novo modelo numérico para a simulação do problemade dispersão de poluentes na atmosfera. O modelo empregou malhas computacionais des-critas no sistema de coordenadas generalizadas e a discretização da equação diferencialgovernante através do método de diferenças finitas. Contudo, com a finalidade de evitarproblemas de instabilidade, o termo convectivo da equação foi aproximado por meio doesquema upwind de primeira ordem.

De acordo com o Capítulo 5, os resultados mostraram-se bastante satisfatóriosno sentido de comprovação teórica dos processos dispersivos. O experimento utilizandoum duto simples, subdividido em quatro casos, possibilitou a visualização de cada pro-cesso envolvido. Os experimentos mais complexos, com a presença de um obstáculo namalha computacional, proporcionou verificar a influência de um vórtice, formado pelascaracterísticas da geometria e do escoamento, sobre o campo de concentração do poluente.

A partir deste trabalho, pretende-se, como trabalhos futuros, realizar uma com-provação prática do modelo desenvolvido. Essa comprovação poderá ser efetuada de duasformas: a primeira, utilizando as características de um problema resolvido por um trabalhorelacionado e verificando, por meio de comparações, os resultados; a segunda, calculandoa solução analítica existente para a equação de transporte unidimensional, dadas as con-dições iniciais e um campo de velocidade uniforme, e comparando unidimensionalmenteos resultados, através do cálculo de erros. Evidentemente, a segunda maneira possibilitauma validação mais precisa do modelo.

Como mencionado no Capítulo 5, planeja-se construir dois experimentos: um câ-nion urbano e um estudo de caso. O primeiro experimento necessita a implementação deuma técnica denominada como multiblocos, devido à complexidade da malha computa-cional. O método de multiblocos consiste na separação de malha em pequenos blocos,sobre os quais se aplicam os modelos numéricos. Dessa forma, a utilização dessa técnicatambém exige uma reformulação dos algoritmos de resolução numérica das Equações deNavier-Stokes e de transporte.

O segundo experimento, no qual uma chaminé se contitui a fonte de emissão dopoluente, sendo transportado pelo fluxo de vento em direção a uma barreira, não foi apli-cado no momento, pela impossibilidade de alteração necessária no algoritmo de resoluçãonumérica da Equação de Navier-Stokes. Além disso, a dificuldade de obtenção dos pa-râmetros necessários em ambas as equações diferenciais governantes, por meio de umaanálise física e química, impede a execução de tal experimento neste instante.

Outros trabalhos futuros podem ser listados, tais como a utilização de técnicas

88

de paralelização dos algoritmos e o desenvolvimento de um modelo numérico em trêsdimensões. A execução de forma paralela dos algoritmos permite a obtenção de resultadosmelhores em um tempo de processamento reduzido. Claramente, um modelo numérico 3Dpossibilita resolver uma gama ainda maior de problemas, solucionando problemas maisrealísticos.

No que diz respeito às contribuições esperadas deste trabalho, deseja-se que oestudo realizado sirva como referência na avaliação dos modelos existentes, e o sistemadesenvolvido, atrelado às propostas de modelagem, em coordenadas generalizadas, paraa simulação numérica da equação de transporte bidimensional, auxilie como uma novatécnica para detecção e prevenção de impactos ambientais e na saúde humana. Dessamaneira, espera-se que incentive a criação de medidas de controle eficazes contra a emissãode poluentes na atmosfera e, por conseguinte, melhore a qualidade de vida da populaçãomundial e preserve os ecossistemas.

89

REFERÊNCIAS

[1] STOCKIE, J. M. The mathematics of atmospheric dispersion modeling. SIAMReview, v. 53, n. 2, p. 349–372, 2011. Disponível em: <http://dblp.uni-trier.de/db/journals/siamrev/siamrev53.html#Stockie11>.

[2] BARBA, A. N. D. Estudo e implementação de esquema upwind na resolução deum modelo de dinâmica dos fluidos computacional em coordenadas generalizadas.Dissertação (Mestrado) — Universidade Estadual de Londrina, 2015.

[3] PAOLI, F. D. Simulação em Túnel de Vento da Dispersão de uma Pluma Emitidapor uma Chaminé Isolada. Dissertação (Mestrado) — Universidade Federal do RioGrande do Sul, 2006.

[4] DALY, A.; ZANNETTI, P. An Introduction to Air Pollution – Definitions,Classifications, and History. [S.l.]: The Arab School for Science and Technology(ASST) and The EnviroComp Institute, 2007.

[5] BOÇON, F. T. Modelagem Matemática do Escoamento e da Dispersão de Poluentesna Microescala Atmosférica. Tese (Doutorado) — Universidade Federal de SantaCatarina (UFSC), Florianópolis, SC, jun 1998.

[6] STERN, A. Fundamentals of Air Pollution. Elsevier Science, 2014. ISBN9780323147675. Disponível em: <https://books.google.com.br/books?id=\_mkV\_0NNYQIC>.

[7] ZANNETTI, P. Air Pollution Modeling: Theories, Computational Methods andAvailable Software. Springer US, 2013. ISBN 9781475744651. Disponível em:<https://books.google.com.br/books?id=q6vwBwAAQBAJ>.

[8] SAYMA, A. Computational Fluid Dynamics. [S.l.]: ApS, 2009.

[9] CIRILO, E. R.; BORTOLI Álvaro L. Geração da malha da traquéia e dos tubosbronquiais por splines cúbico. Semina: Ciências Exatas e Tecnológicas, v. 27, n. 2,p. 147–155, 2006.

[10] BOYCE, W. E.; DIPRIMA, R. C. Equações Diferenciais Elementares e Problemasde Valores de Contorno. eighth. [S.l.]: LTC Editora, 2002.

[11] MIERSEMANN, E. Partial Differential Equations. CreateSpace IndependentPublishing Platform, 2014. ISBN 9781502910967. Disponível em: <https://books.google.com.br/books?id=ntb2oQEACAAJ>.

[12] MENDELSON, E.; AYRES, F. Teoria e problemas de cálculo. BOOKMANCOMPANHIA ED, 2007. ISBN 9788560031092. Disponível em: <https://books.google.com.br/books?id=4bsoo853eEgC>.

[13] CORRÊA, L.; LIMA, G. A. B. de; FERREIRA, V. G. Solução numérica deequações diferenciais parciais via o método das diferenças finitas. In: II Colóquio deMatemática do Centro Oeste. [S.l.: s.n.], 2011.

90

[14] CENEDESE, E. Solução das Equações de Burgers e de Navier-Stokes BidimensionaisUtilizando a Técnica da Transformada Integral Generalizada. Dissertação (Mestrado)— Universidade Estadual Paulista, Ilha Solteira, feb 2005.

[15] FORTUNA, A. de O. Técnicas Computacionais para Dinâmica dos Fluidos:Conceitos Básicos e Aplicações. São Paulo: Edusp - Editora da Universidade de SãoPaulo, 2000.

[16] COSTA, L. T.; RIBEIRO, M. C. C. Propriedades dinâmicas de fluidos por simulaçãocomputacional: Métodos híbridos atomístico-contínuo. Química Nova, Scielo, v. 33,p. 938 – 944, 00 2010. ISSN 0100-4042. Disponível em: <http://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-40422010000400032&nrm=iso>.

[17] NACHBIN, A.; TABAK, E.; JANEIRO, B. C. B. de Matemática (21. : 1997 :Rio de. 21° Colóquio Brasileiro de Matemática: Equações diferenciais em modelagemmatemática computacional. Conselho Nacional de Desenvolvimento Científico eTechnológico, Instituto de Matemática Pura e Aplicada, 1997. ISBN 9788524401275.Disponível em: <https://books.google.com.br/books?id=VIBCAQAACAAJ>.

[18] KUZMIN, D. A Guide to Numerical Methods for Transport Equations. [S.l.: s.n.],2010.

[19] STOCKER, T. Introduction to Climate Modelling. Springer Berlin Heidelberg,2011. (Advances in Geophysical and Environmental Mechanics and Mathematics).ISBN 9783642007736. Disponível em: <https://books.google.com.br/books?id=D4zulgFb5JwC>.

[20] BERN, M.; PLASSMANN, P. Mesh generation. In: Handbook of ComputationalGeometry. Elsevier Science. [S.l.: s.n.], 2000. p. 291–332.

[21] MARQUES, A. C. Desenvolvimento de Modelo Numérico Utilizando o Métododos Volumes Finitos em Malhas Não-estruturadas. Dissertação (Mestrado) —Universidade Federal do Espírito Santo (UFES), Vitória, ES, 2005.

[22] MALISKA, C. R. Transferência de Calor e Mecânica dos Fluidos Computacional.second. [S.l.]: LTC Editora, 1995.

[23] KOOMULLIL, R.; SONI, B.; SINGH, R. A comprehensive generalized meshsystem for cfd applications. Mathematics and Computers in Simulation, v. 78,n. 5–6, p. 605 – 617, 2008. ISSN 0378-4754. Applied Scientific Computing:Numerical Grid Generation and Field Simulation. Disponível em: <http://www.sciencedirect.com/science/article/pii/S0378475408001626>.

[24] THOMPSON, J. F.; THAMES, F. C.; MASTIN, C. W. Tomcat — a codefor numerical generation of boundary-fitted curvilinear coordinate systems onfields containing any number of arbitrary two-dimensional bodies. Journal ofComputational Physics, v. 24, n. 3, p. 274 – 302, 1977. ISSN 0021-9991. Disponívelem: <http://www.sciencedirect.com/science/article/pii/0021999177900389>.

[25] THOMPSON, J. F.; WARSI, Z. U. A.; MASTIN, C. W. Numerical Grid Generation:Foundations and Applications. [S.l.]: Elsevier Science Publishing Co., 1985.

91

[26] FARLOW, S. An Introduction to Differential Equations and Their Applications.Dover Publications, 2012. (Dover Books on Mathematics). ISBN 9780486135137.Disponível em: <https://books.google.com.br/books?id=r0QI5DMr6WAC>.

[27] FERZIGER, J.; PERIĆ, M. Computational methods for fluid dynamics.Berlin, New York: Springer, 2002. ISBN 978-3-540-42074-3. Disponível em:<http://opac.inria.fr/record=b1135325>.

[28] QUEIROZ, R. A. B. de; FERREIRA, V. G. Esquemas polinomiais upwind e suasaplicações em escoamentos incompressíveis 3d transientes. In: UNIVERSIDADEFEDERAL DE UBERLâNDIA (UFU). XIV Congresso Nacional de Estudantes deEngenharia Mecânica. UberlÂndia, MG, 2007.

[29] FANG, F. et al. Reduced order modelling of an unstructured mesh airpollution model and application in 2D/3D urban street canyons. AtmosphericEnvironment, v. 96, p. 96–106, out. 2014. ISSN 13522310. Disponível em:<http://www.sciencedirect.com/science/article/pii/S1352231014005421>.

[30] FERRAGUT, L. et al. An efficient algorithm for solving a multi-layerconvection–diffusion problem applied to air pollution problems. Advances inEngineering Software, v. 65, p. 191–199, nov. 2013. ISSN 09659978. Disponível em:<http://www.sciencedirect.com/science/article/pii/S0965997813001087>.

[31] MADALOZZO, D. et al. Numerical simulation of pollutant dispersion in streetcanyons: Geometric and thermal effects. Applied Mathematical Modelling,v. 38, n. 24, p. 5883–5909, dez. 2014. ISSN 0307904X. Disponível em:<http://www.sciencedirect.com/science/article/pii/S0307904X1400211X>.

[32] STEFFENS, J. T. et al. Modeling the effects of a solid barrier on pollutantdispersion under various atmospheric stability conditions. AtmosphericEnvironment, v. 69, p. 76–85, abr. 2013. ISSN 13522310. Disponível em:<http://www.sciencedirect.com/science/article/pii/S1352231012011193>.

[33] ZHONG, J.; CAI, X.-M.; BLOSS, W. J. Modelling the dispersion and transportof reactive pollutants in a deep urban street canyon: Using large-eddy simulation.Environmental Pollution, v. 200, p. 42–52, maio 2015. ISSN 02697491. Disponívelem: <http://www.sciencedirect.com/science/article/pii/S0269749115000779>.

[34] GOUSSEAU, P.; BLOCKEN, B.; HEIJST, G. van. Large-Eddy Simulationof pollutant dispersion around a cubical building: Analysis of the turbulentmass transport mechanism by unsteady concentration and velocity statistics.Environmental Pollution, v. 167, p. 47–57, ago. 2012. ISSN 02697491. Disponívelem: <http://www.sciencedirect.com/science/article/pii/S0269749112001339>.

[35] GOUSSEAU, P. et al. Near-field pollutant dispersion in an actual urbanarea: Analysis of the mass transport mechanism by high-resolution LargeEddy Simulations. Computers & Fluids, v. 114, p. 151–162, jul. 2015. ISSN00457930. Disponível em: <http://www.sciencedirect.com/science/article/pii/S0045793015000614>.

[36] KIKUMOTO, H.; OOKA, R. A numerical study of air pollutant dispersion withbimolecular chemical reactions in an urban street canyon using large-eddy simulation.

92

Atmospheric Environment, v. 54, p. 456–464, jul. 2012. ISSN 13522310. Disponívelem: <http://www.sciencedirect.com/science/article/pii/S1352231012001549>.

[37] LATEB, M. et al. Simulation of near-field dispersion of pollutants usingdetached-eddy simulation. Computers & Fluids, v. 100, p. 308–320, set. 2014.ISSN 00457930. Disponível em: <http://www.sciencedirect.com/science/article/pii/S0045793014002229>.

[38] FENG, X. et al. Artificial neural networks forecasting of PM2.5 pollution using airmass trajectory based geographic model and wavelet transformation. AtmosphericEnvironment, v. 107, p. 118–128, abr. 2015. ISSN 13522310. Disponível em:<http://linkinghub.elsevier.com/retrieve/pii/S1352231015001491>.

[39] SCHAAP, M. et al. Performance of European chemistry transport models asfunction of horizontal resolution. Atmospheric Environment, v. 112, p. 90–105, jul.2015. ISSN 13522310. Disponível em: <http://www.sciencedirect.com/science/article/pii/S1352231015300066>.

[40] OLIVEIRA, F. S. D. et al. Malhas móveis para solução numérica de equaçõesdiferenciais parciais. Revista de Sistemas de Informação da FSMA, 2013.

[41] FERREIRA, V. G. Análise e Implementação de Esquemas de Convecção e Modelosde Turbulência para Simulação de Escoamentos Incompressíveis EnvolvendoSuperfícies Livres. Tese (Doutorado) — Universidade de São Paulo, oct 2001.

[42] TAYLOR, J. Mecânica Clássica. Bookman Editora, 2013. ISBN 9788582600887.Disponível em: <https://books.google.com.br/books?id=A6Y5AgAAQBAJ>.

[43] REYNOLDS, O. An experimental investigation of the circumstances whichdetermine whether the motion of water shall be direct or sinuous, and of the law ofresistance in parallel channels. Philosophical Transactions of the Royal Society ofLondon, The Royal Society, v. 174, p. 935–982, 1883. ISSN 02610523. Disponívelem: <http://www.jstor.org/stable/109431>.