Latest revision as of 16:22, 23 February 2021

Resumo

Este trabalho apresenta a simulação hidráulica do escoamento em um reservatório, utilizando o método dos volumes finitos para resolver o sistema de equações que modelam o fluxo bidimensional em águas rasas, negligenciando as tensões tangenciais. Para resolver as equações, utilizou-se um esquema híbrido de volumes finitos. A solução do sistema de equações linearizadas foi obtida por implementação computacional do método iterativo de Gauss-Seidel. Para ilustrar a aplicação do esquema numérico em engenharia hidráulica, apresenta-se o estudo do escoamento em um reservatório subterrâneo com pilares internos, cujo fluxo é gerado pela abertura de dois portões de saída. O código desenvolvido permite a análise temporal do volume e do fluxo no reservatório, a interpretação gráfica do fenômeno físico investigado bem como o cálculo da precisão do modelo. Os resultados obtidos mostram a interpretação física esperada, boa concordância com os dados da literatura, boa precisão, estabilidade e baixo custo computacional.

Palavras chave: Equações de águas rasas, volumes finitos, interpolação híbrida, estabilidade numérica, simulação de leito seco, reservatório de contenção

Abstract

This work presents the numeric simulation of flow in a reservoir, using the finite volume method for solver the system of equations that model the two-dimensional flow in shallow water, neglecting the tangential tensions. The solution of the system of linearized equations was obtained by computational implementation of the Gauss-Seidel iterative method. To illustrate the application of the numerical scheme in hydraulic engineering, it was considered the study of flow in an underground reservoir with internal pillars, whose flow is generated by the opening of two outlet gates. The employed code allowed the temporal analysis of the volume and the flow in the reservoir, the graphical interpretation of the investigated physical phenomenon as well as, the calculation of the precision of the model. The results obtained show the expected physical interpretation, good agreement with literature data, good precision,stability and low computational cost.

Keywords: Shallow water equations, finite volumes, hybrid interpolation, numerical stability, dry bed simulation, containment reservoir

1. Introdução

A ocorrência ou agravamento de inundações nas grandes cidades, bem como a quebra de barragens, tornaram-se problemas cada vez mais recorrentes, causando mortes e fortes impactos ambientais. Nesse cenário, fica evidente a necessidade de conhecimento prévio do comportamento dos fluxos de superfície livre e simulações de abertura de comportas de barragens para manejo de recursos hídricos, avaliação de riscos e usos para obtenção de melhor disponibilidade visando a preservação ambiental. Para tanto, é necessário utilizar metodologias que possam descrever as variáveis hidrodinâmicas envolvidas, como campo de velocidade do fluxo, turbulência e mudanças temporais na altura da superfície da água [1,2]. Em [3] sugere-se o uso de modelos matemáticos para entender os padrões de fluxo em corpos de água. Modelos matemáticos adequados são ferramentas úteis no gerenciamento de recursos hídricos, pois reduzem o tempo de análise, custos e riscos envolvidos em levar dados experimentais em fluxos reais auxiliando na tomada de decisões e na definição de ações corretivas diante de vazamentos iminentes. Matematicamente, os fluxos de superfície livre podem ser descritos de várias maneiras. Uma abordagem mais simples e muito aplicada é dada pelas Equações de Águas Rasas, que representam uma simplificação das Equações de Navier-Stokes e é aplicável aos fluxos nos quais as tensões e acelerações verticais podem ser negligenciadas sem prejuízo do resultado. Essas equações são derivadas por meio da integração das equações governantes (Navier-Stokes e Continuidade) na direção vertical. Sua derivação está fora do escopo deste trabalho, mas pode ser vista em detalhes em [4] . Esse sistema de equações, embora definido em um domínio bidimensional, fornece a altura da coluna de água e as componentes de velocidade em um fluxo instável. No entanto, a obtenção de uma solução analítica para o sistema de equações de águas rasas é, na maioria dos casos, uma tarefa difícil, pois trata-se de um conjunto de equações diferenciais parciais hiperbólicas não lineares com termos fontes. O que requer um tratamento numérico robusto para evitar instabilidades numéricas. Ao longo das últimas três décadas, muitos avanços foram realizados na solução numérica das equações de águas rasas para simulação de escoamentos unidimensionais aplicados a hidráulica fluvial em rios, canais e tubulações [5,6,7,8,9], e modelos bidimensionais aplicados à análise de escoamentos causados pela abertura instantânea de comportas e ruptura de barragens [10,11,12,13,14,15]. Em [10] utilizaram esquemas explícitos de diferenças finitas, na discretização das equações de águas rasas bidimensionais, adicionando condições de estabilidade e viscosidade artificial para suavizar as oscilações de alta frequência. Para validar, os esquemas discretos foram aplicados em um problema típico de engenharia hidráulica: abertura parcial de uma barragem e passagem de uma onda de inundação através de uma contração do canal, adotando a razão entre o nível de água a jusante ${\displaystyle (h_{r})}$ e à montante ${\displaystyle (h_{l})}$ como sendo, ${\displaystyle h_{r}/h_{l}=0.5}$. Os resultados dos diferentes esquemas foram comparados e a ocorrência de oscilações numéricas foi observada, evidenciando a necessidade de se utilizar esquemas implícitos. O estudo de caso apontado em [10] foi analisado por [15] usando um esquema de elementos finitos de alta ordem e uma melhoria na relação ${\displaystyle h_{r}/h_{l}}$ para ${\displaystyle 0,05}$ foi obtida, buscando simular a condição de leito seco em uma determinada região do domínio. Em [12] usaram dados experimentais em um evento de quebra de barragem para verificar a precisão, estabilidade e confiabilidade de uma solução numérica das equações de águas rasas obtidas por um esquema de volumes finitos bidimensional. O modelo numérico permitiu uma precisão de segunda ordem, tanto no espaço como no tempo, acarretando porém um alto custo computacional, apresentando algumas discrepâncias entre os resultados reais dos medidores e os resultados numéricos nos pontos pesquisados. Os autores destacam a necessidade de simular a propagação de uma onda de inundação para um fundo seco e áspero e otimizar o algoritmo para reduzir o tempo de computação. Um esquema de volumes finitos não estruturado para malhas triangulares para simular fluxos bidimensionais instáveis em águas rasas com frentes molhadas/secas sobre topografia complexa é apresentado em [11]. Utiliza-se em [11] para o esquema temporal o método explícito de Runge-Kutta e para o espacial a malha não estruturada. Comparações satisfatórias foram obtidas com dados mensurados e numéricos. Todavia, os autores sugerem a incorporação de um procedimento de escalonamento num esquema implícito temporal, e uma implementação paralela do algoritmo, a fim de reduzir o tempo computacional. A eficiência de um esquema implícito de volumes finitos é explorada em [13] para simulação de fluidos em águas rasas bidimensionais, com frentes secas/molhadas. Adotou-se em [13] um tratamento complexo para a discretização espacial, fazendo uso de uma malha não estruturada flexível para obter um melhor ajuste à geometrias irregulares. No entanto, é feito ressalvas quanto ao tempo computacional para obtenção da solução por esquemas implícitos, pois embora um método implícito normalmente exija menos etapas da solução, cada uma exige mais tempo computacional do que as soluções obtidas por um esquema explícito. Os autores apontam a necessidade de se utilizar um solucionador de matrizes que com um número reduzido de iterações alcance a convergência. Neste contexto, o objetivo principal deste trabalho é formular um esquema híbrido e totalmente implícito de volumes finitos para discretização das equações de águas rasas, a fim de obter uma solução numérica estável com um algoritmo iterativo otimizado, com convergência rápida e baixo custo computacional, para aplicação na simulação de fluxos bidimensionais transitórios em topografias complexas (com declives e sujeitos a resistência ao fluxo). Condições iniciais descontínuas representando leitos úmidos/secos foram adotadas. Para tanto, o trabalho está organizado da seguinte forma: A Seção 2 apresenta as equações hiperbólicas, que regem a natureza das águas rasas, considerando as perdas de atrito e declives da topografia do domínio computacional e detalhando a configuração numérica destas equações; na Seção 3, o modelo é validado e uma comparação é feita com os resultados obtidos com o trabalho apresentado em [10]. A robustez da solução numérica é verificada com uma extrapolação na relação ${\displaystyle h_{r}/h_{l}}$ para ${\displaystyle 0.00002}$, na Seção 4; o esquema é aplicado na análise de um problema hidráulico real na Seção 5 e; as conclusões são apresentadas na Seção 6.

2. Modelo matemático e solução numérica

Na análise do escoamento bidimensional transiente em um canal, cuja topografia do fundo é descrita por uma superfície ${\displaystyle z_{0}(x,y)}$, busca-se determinar o nível do fluido de superfície livre ${\displaystyle h(x,y,t)}$ e o campo de velocidade do fluxo ${\displaystyle U(u,v)}$ para cada instante de tempo ${\displaystyle t}$, conforme Figura 1.

Figura 1. Notação adotada na análise de fluxo em águas rasas

2.1 Sistema de Equações Governantes

O sistema de equações diferenciais que descreve o escoamento hidrostático de um fluido incompressível em águas rasas com superfície superior livre é obtido pela integração ao longo da profundidade das equações tridimensionais que modelam a quantidade de movimento e conservação de massa num escoamento, assumindo distribuição de velocidade uniforme na direção vertical. Como em Mahmood and Yevjevich [3], adota-se as seguintes simplificações: (i) a distribuição de pressão na vertical é puramente hidrostática; (ii) não há contribuições laterais ao escoamento no reservatório analisado; (iii) a massa específica e a viscosidade do fluido no reservatório não sofrem variações significativas ao longo do tempo e do espaço; (iv) as tensões tangenciais são desprezíveis. Deste modo, as equações governantes do modelo hidrodinâmico, também ditas equações de águas rasas ou equações bidimensionais de Saint Venant, podem ser escritas na forma conservativa (com os fluxos dentro do sinal da derivada):

onde ${\displaystyle x}$ e ${\displaystyle y}$ representam as componentes cartesianas nas direções longitudinal e transversal, respectivamente; ${\displaystyle h}$ é profundidade da água no reservatório; ${\displaystyle u}$ e ${\displaystyle v}$ são as componentes cartesianas da velocidade; ${\displaystyle g}$ é a aceleração da gravidade, ${\displaystyle t}$ é o tempo, e ${\displaystyle S_{y}}$ e ${\displaystyle S_{x}}$ são os declives do fundo do reservatório definidos como:

enquanto os declives de resistência ao escoamento (perda de energia), nas direções ${\displaystyle x}$ e ${\displaystyle y}$, respectivamente, são dadas por:

onde ${\displaystyle c_{M}(m/s^{\frac {1}{3}})}$ é o coeficiente de rugosidade de Manning, que quantifica o atrito do fluxo com o fundo do canal. As equações de águas rasas (1) - (3) não possuem solução analítica fechada, portanto, uma abordagem numérica é empregada para resolvê-las de forma discreta e linearizada, onde o domínio é dividido em um número finito de células resultando em um sistema linear algébrico que é resolvido por um algoritmo computacional apropriado codificado em FORTRAN 90. Neste estudo, diferentes condições de contorno e topografia de fundo foram impostas para verificar a robustez e aplicabilidade do código. Na etapa inicial, o nível da água é sempre conhecido e todas as velocidades são iguais a zero para todo o domínio. Essas condições estão detalhadas nas Subseções 3.1 e 3.2.

2.2 Discretização do domínio e linearização das equações governantes

A solução numérica do modelo hidrodinâmico descrito pelas equações (1)-(3) é obtida por meio do método dos volumes finitos (MVF). Os fundamentos matemáticos e físicos do MVF são bem discutidos por [16,17,18,19,20]. De acordo com [20], o MVF consiste na integração espacial e temporal das equações diferenciais na forma conservativa (com todos os fluxos envolvidos pelo operador derivada) em um volume genérico de controle. O sistema de equações (1)-(3) é integrado espacialmente (sobre um volume de controle retangular) e no tempo (de ${\displaystyle t}$ a ${\displaystyle t+\Delta t}$) e é empregado para a integração temporal a seguinte aproximação numérica:

Adota-se aqui uma notação mais compacta:

E assim, a Eq. (6) pode ser reescrita como:

Neste artigo, adota-se a formulação totalmente implícita, ${\displaystyle (\theta =1)}$, ou seja, todas as derivadas são avaliadas no nível de tempo ${\displaystyle t+\Delta \,t}$. Com uma implementação mais complexa, a formulação totalmente implícita é limitada apenas pela precisão e não sofre instabilidades numéricas devido à mudança de coeficientes de sinal pela variação do intervalo de tempo. E assim, permite o uso de etapas de tempo maiores com menor tempo de computação. Esse é um método incondicionalmente estável no tempo. Obviamente, esta estabilidade pode ser alterada devido às não linearidades do acoplamento e ao comportamento dos termos de fontes, como a perda de energia devido à declividade, o que pode gerar oscilações numéricas. Portanto, faz-se necessário um tratamento cuidadoso na obtenção de uma solução numérica estável e também precisa. A derivada temporal é avaliada com uma aproximação de Euler regressiva de ${\displaystyle 1^{a}}$ ordem:

onde ${\displaystyle \phi }$ representa uma variável. Na integração espacial é utilizado o Teorema Fundamental do Cálculo para reescrever as integrais das equações (1)-(3) como um sistema de equações constituído pelos valores das variáveis avaliadas nos pontos de integração ${\displaystyle w,e,s,n}$, que por sua vez, são avaliados por funções de interpolação. Para o cálculo destes valores, adota-se uma combinação dos esquemas de diferenciação central de ${\displaystyle 2^{a}}$ ordem (CDS - Central Differencing Scheme) e um esquema upwind (UDS - Upwind Differencing Scheme) de ${\displaystyle 1^{a}}$ ordem. Na interpolação central, os valores das variáveis nas faces do volume de controle são dados pelo valor médio dos centróides vizinhos, isto é,

No esquema upwind, os valores das variáveis nas faces do volume de controle são dados considerando primeiramente a direção do fluxo que transporta alguma propriedade do centróide vizinho para a face onde o fluxo é avaliado, por exemplo,

O sistema (1) - (3) apresenta produtos de variáveis envolvendo ${\displaystyle u}$, ${\displaystyle v}$ e ${\displaystyle h}$. Esta característica é a natureza não-linear das equações de Saint-Venant e uma das razões pelas quais não existe uma solução analítica para problemas complexos de águas rasas, e assim, o sistema deve ser linearizado para permitir que o método numérico gere um sistema linear resolvido por qualquer solucionador linear. A linearização é feita considerando no produto de variáveis apenas uma delas como a variável ativa, e as outras são constantes ou variáveis secundárias calculadas a partir da condição inicial da iteração de looping do tempo anterior e/ou do último passo de tempo, por exemplo,

onde neste exemplo, o índice ${\displaystyle 0}$ representa os valores disponíveis do cálculo anterior e esses valores são usados para calcular o coeficiente da matriz para a variável ativa ${\displaystyle h}$. Naturalmente, outras equações estão disponíveis para resolver ${\displaystyle u}$ e ${\displaystyle v}$ como a variável ativa (ou variável principal), e estas são empregadas para calcular cada coeficiente, assim, quando uma nova solução é alcançada, um novo cálculo é iniciado e novos coeficientes são definidos, e uma nova solução é executada até a convergência. Na combinação CDS/UDS, o esquema UDS é escolhido para a variável principal e a aproximação CDS para a variável secundária conhecida, por exemplo na integral da Eq. (1) o UDS é usado para a variável principal ${\displaystyle h}$, e para ${\displaystyle u}$ e ${\displaystyle v}$ é aplicada a interpolação CDS usando os valores nodais vizinhos conhecidos ${\displaystyle u^{0},\,v^{0}}$ do nó ${\displaystyle P}$. Todos os valores nas faces são funções dos valores nodais ${\displaystyle P,W,E,S}$, usando coeficientes ${\displaystyle \alpha ^{,}s}$ que representam o sinal das componentes ${\displaystyle u}$ e ${\displaystyle v}$ nas faces do volume de controle, isto é, a direção do vetor de velocidade:

onde,

Portanto, a partir da formulação totalmente implícita e da interpolação híbrida CDS/UDS, a equação algébrica discreta para a Eq. (1) é dada por:

onde,

Para integrar a Eq. (2) aplica-se o Teorema Fundamental do Cálculo considerando a formulação totalmente implícita. Deste modo, o lado esquerdo desta equação é reescrita em função das variáveis calculadas em ${\displaystyle P}$ e nas faces ${\displaystyle w,e,s,n}$, conforme a expressão ${\displaystyle (}$26${\displaystyle )}$:

Na expressão ${\displaystyle (}$26${\displaystyle )}$ o termo ${\displaystyle u^{2}}$ é escrito como ${\displaystyle (u^{0}\,u)}$ onde ${\displaystyle u^{0}}$ é conhecido e ${\displaystyle u}$ é a variável ativa a ser calculada. A interpolação UDS é usada para avaliar ${\displaystyle u}$, enquanto a interpolação CDS para avaliar ${\displaystyle u^{0}}$ nas faces. Por exemplo,

com ${\displaystyle \alpha _{e}}$ definido na equação (17). O termo quadrático da expressão (26) da variável secundária ${\displaystyle h}$ é calculado usando a interpolação CDS com os valores disponíveis da Eq. (19) e assim,

A mesma abordagem é usada nos termos advectivos da expressão (26). O produto ${\displaystyle (uvh)}$ é dividido em ${\displaystyle (hv^{0})u}$ com UDS para variável ativa ${\displaystyle u}$ e CDS para as variáveis secundárias, ${\displaystyle h}$ atualizada e ${\displaystyle v^{0}}$. Para a integração do lado direito da Eq. (2), que representa a resistência à declividade, o valor de ${\displaystyle h}$ é calculada no centróide ${\displaystyle P}$, pois fisicamente a fonte representa a média integral no volume discreto centrado em ${\displaystyle P}$, isto é, o fonte não flui, é um termo volumétrico. Portanto, ${\displaystyle h=h_{P}}$ e a magnitude da velocidade é constante. Dessa forma, a integral do lado direito da Eq. (2) é reescrita da seguinte forma:

onde os valores da variável principal ${\displaystyle u}$ nas faces são obtidos por UDS. Assim, a forma discreta linearizada da Eq. (2) é dada por:

onde,

com

A mesma abordagem é usada para o cálculo de ${\displaystyle v}$, integrando a Eq. (3). Assim, são obtidos três sistemas lineares a serem resolvidos:

onde ${\displaystyle A^{\phi }}$ representa a matriz de coeficientes da variável ${\displaystyle \phi }$ (${\displaystyle h,u}$ ou ${\displaystyle v}$), enquanto ${\displaystyle B^{\phi }}$ é o vetor que contém as condições de contorno e termos fontes.

3. O algoritmo

O código de simulação computacional foi programado em Linguagem Fortran, compilado na versão Fortran Power Station 4.0 licenciada pela Microsoft Corporation. A escolha por esse programa deve-se à sua velocidade de execução e alta capacidade de armazenamento de dados. Neste trabalho, o domínio físico retangular ${\displaystyle D=[x_{a},x_{b}]\times [y_{a},y_{b}]}$ é dividido em ${\displaystyle N_{x}}$ e ${\displaystyle N_{y}}$ volumes de controle com comprimentos dados por ${\displaystyle \Delta x,\Delta y}$ nas direções ${\displaystyle x}$ e ${\displaystyle y}$, respectivamente. Da mesma forma, o espaço temporal ${\displaystyle [t_{0},t_{f}]}$ é dividido em ${\displaystyle N_{t}}$ passos de tempo uniformes iguais a ${\displaystyle \Delta t}$. A malha cartesiana ordenada e estruturada é varrida de acordo com o par ${\displaystyle (i,j)}$ da esquerda para a direita aumentando o comprimento ${\displaystyle (i-1)\Delta x}$ na direção ${\displaystyle x}$, e de baixo para cima aumentando o comprimento ${\displaystyle (j-1)\Delta y}$ na direção ${\displaystyle y}$, começando no centróide ${\displaystyle (1,1)}$ que representa o ponto ${\displaystyle (x_{a}+\Delta x/2,y_{a}+\Delta y/2)}$. Desta forma, faz-se a seguinte correspondência:,

Após a discretização do domínio, são lidas a topografia da parte inferior, representada pela superfície ${\displaystyle z_{0}(x,y)}$, e as condições iniciais ${\displaystyle h_{0},u_{0}}$ e ${\displaystyle v_{0}}$. O sistema linear é resolvido por um método de Gauss-Seidel de acordo com a direção varrida anteriormente descrita. Ele itera até que um dos dois critérios seja atingido: o número máximo residual ${\displaystyle tol}$ ou máximo de iterações ${\displaystyle itmax}$. Se for atingido resíduo ${\displaystyle R<tol}$ em um número ${\displaystyle k}$ de iterações e ${\displaystyle 2k<itmax}$, então continua-se o processo iterativo até ${\displaystyle 2k}$ iterações. Ou seja, ao atingir a tolerância do resíduo máximo, faz-se mais ${\displaystyle k}$ iterações, isso torno o resíduo bem inferior ao residual máximo permitido. Para este trabalho, o número máximo de iterações para cada intervalo de tempo é fixado em ${\displaystyle itmax=1000}$ iterações, enquanto a tolerância para o resíduo máximo global permitido é igual a ${\displaystyle tol=10^{-11}}$, sendo o resíduo global ${\displaystyle R}$ definido como:

onde ${\displaystyle N_{x}N_{y}}$ é o número total de volume de controles, e os resíduos ${\displaystyle r_{h},r_{u},r_{v}}$ são a soma dos resíduos individuais de cada volume de controle,calculados incrementalmente de acordo com a Eq. (41):

onde ${\displaystyle d_{\phi }}$ representa a diferença absoluta entre os dois lados da equação linearizada para a variável ${\displaystyle \phi }$. O algoritmo completo é apresentado abaixo. O código também permite usar a última solução para inicializar uma nova execução com novas condições de contorno, se necessário.

4. Validação da solução numérica

Com o objetivo de validar a abordagem numérica apresentada na seção anterior, esta seção apresenta uma comparação com os resultados obtidos para um modelo de abertura instantânea de comportas apresentado em [10]. Este problema consiste em um domínio quadrado com ${\displaystyle 200\,m}$ de comprimento com uma barragem interna com ${\displaystyle 10\,m}$ de largura, disposta na faixa de ${\displaystyle 90-100\,m}$, que divide o domínio como é mostrado na Figura 2. Para simulação da abertura de comportas, considera-se uma abertura não simétrica e instantânea com ${\displaystyle 75\,m}$ de largura entre ${\displaystyle 30}$ e ${\displaystyle 105\,m}$. Essa abertura causa uma liberação sutil da coluna de água a jusante, resultando em um fluxo com gradientes espaciais e temporais acentuados.

Figura 2. Domínio do modelo

Todos os parâmetros numéricos e físicos são mostrados na Tabela (1). De acordo com [10], neste problema a resistência ao escoamento é desprezada e sua parte inferior é plana.

4.1 Condições iniciais

Para o tempo inicial, o fluido é quiescente, ou seja, ${\displaystyle u=v=0}$ para todo o domínio, enquanto a coluna d'água inicial é dada por:

com o nível da água à montante ${\displaystyle h_{l}=10\,m,\,}$ e a jusante ${\displaystyle h_{r}=5\,m}$. Consequentemente, a razão entre os níveis de água é ${\displaystyle h_{r}/h_{l}=0.5\,m.}$

4.2 Condições de Contorno

Por uma questão de clareza, esta seção é apresentada aqui para mostrar de forma prática como as condições de contorno são tratadas neste trabalho usando para isso, o caso de validação como exemplo. As abordagens apresentadas aqui também são aplicadas a outros casos de maneira similar, respeitando, é claro, a configuração geométrica de cada um. Como apontado em [21], existem pelo menos três maneiras de lidar com as condições de contorno: realizar um balanço para cada volume de controle na fronteira, chamado simplesmente de equilíbrio de contorno; considerar um meio volume na fronteira com seu centróide pertencente ao próprio limite; e, por último, considerar volumes fictícios em torno dos limites, isso significa que o valor de limite prescrito deve ser igual à interpolação (qualquer que seja) entre os centróides do volume real de controle e o novo volume de controle fictício. Esta última abordagem é preferida, pois a primeira implica um tratamento matemático distinto para cada fronteira, o que torna o código mais complexo de acordo com a geometria. O caminho de ${\displaystyle 2^{a}}$ ordem, o meio volume, não é conservador porque o valor é imposto, em vez disso, como resultado de um balanço de fluxo. A abordagem de volume fictício permite que toda a malha seja constituída por volumes inteiros, o que é fácil de codificar e tem a mesma natureza conservadora para todos os volumes. Embora este artifício aumente o número de incógnitas do problema, acarretando custo computacional, é um procedimento útil com codificação fácil. Na Figura 3 é apresentado um esboço deste conceito.

Figura 3. Volumes fictícios ao redor da malha real

Utiliza-se a Condição de Neumann nula para o nível da água nas paredes externas e nos contornos da represa. Desse modo,

onde ${\displaystyle w}$ e ${\displaystyle e}$ são os pontos dos contornos esquerdo e direito, respectivamente; e ${\displaystyle s}$ e ${\displaystyle n}$ são os pontos dos contornos inferior e superior, respectivamente, como descrito na Figura 4.

Figura 4. Esquema da condição de contorno para o nível ${\displaystyle h}$ do fluido

Todas as derivadas nas faces fronteiras (${\displaystyle w,e,n,}$ e ${\displaystyle s}$) podem ser interpoladas como funções dos centróides vizinhos destacados na Figura 4: um externo fictício (${\displaystyle W,E,S}$ ou ${\displaystyle N}$) e um interno ${\displaystyle (P)}$ pertencente ao domínio. Consequentemente,

o que resulta em,

Para as variáveis de fluxo (${\displaystyle u}$ e ${\displaystyle v}$), as paredes são consideradas impermeáveis e todas as componentes de velocidade normal aos contornos são nulas. Esta condição de Dirichlet é imposta considerando o esquema de diferenciação central, ou seja,