Two-dimensional numerical simulation of channel flow with submerged obstacles using the lattice Boltzmann method

Subm: 25/04/2017	

Accep.: 27/06/2017

==Resumo==

Apresenta-se a simulação numérica bidimensional do escoamento de água em um canal retangular com chicanas submersas distribuídas alternadamente nas suas margens. As equações governantes do escoamento são as equações de águas rasas, que serão resolvidas por meio do método do reticulado de Boltzmann (LBM) com múltiplos tempos de relaxamento (MRT). Utilizou-se condições de contorno adequadas para escoamento em canal, vazão e altura da água constante na entrada e saída, respectivamente. Devido às características do problema que pretende-se simular, incorporou-se ao código computacional uma técnica de simulação em grandes escalas (LES - ''Large Eddy Simulation'') a qual permite obter resultados mais próximos do comportamento real do escoamento. Além disso, avalia-se a estabilidade da simulação em todos os pontos da malha para cada passo de tempo e, juntamente com a propriedade da consistência do LBM, obtém-se a convergência da solução. A simulação fornece a profundidade, velocidades nas direções <math display="inline">x</math> e <math display="inline">y</math>, e a magnitude da vorticidade da água.

''Palavras-chave'': Método do reticulado de Boltzmann; equações de águas rasas; escoamento em canal; obstáculos submersos.

==Summary==

A two-dimensional numerical simulation of the water flow in a rectangular channel with submerged obstacles distributed alternately along its banks is presented. The governing equations of flow are the shallow water equations, which will be solved by the Boltzmann lattice method (LBM) with multiple relaxation times (MRT). The non-slip bounce-back scheme was used on walls and obstacles, constant discharge at the inlet and fixed depth at the outlet of the channel. Due to the characteristics of the problem to be simulated, a large eddy simulation (LES) technique was incorporated into the computational code, which allows to obtain results that are closer to the actual behavior of the flow. In addition, the stability of the simulation at all points of the mesh is evaluated for each step of time and, together with the property of the consistency of the LBM, the convergence of the solution is obtained. The simulation provides the depth, velocities in the x and y directions, and the magnitude of water vorticity.

==1. Introdução==

O Método do reticulado do Boltzmann é eficiente e mostra-se competitivo nas aplicações de mecânica dos fluidos computacional, sendo especialmente útil em aplicações com geometrias complexas <span id='citeF-1'></span>[[#cite-1|[1]]]. Uma de suas principais características é a abordagem diferenciada quando comparado aos métodos numéricos tradicionais (volumes finitos <span id='citeF-2'></span>[[#cite-2|[2]]], diferenças finitas <span id='citeF-3'></span>[[#cite-3|[3]]], método das características <span id='citeF-4'></span>[[#cite-4|[4]]], elementos finitos <span id='citeF-5'></span>[[#cite-5|[5]]]), pois não utiliza discretizações nas equações macroscópicas governantes do escoamento <span id='citeF-6'></span>[[#cite-6|[6]]]. O LBM considera o comportamento de um conjunto de partículas caracterizando-o como um método em mesoescala, possui equação governante caracterizada pelas etapas de transmissão e colisão, e pode ser visto como uma forma discretizada da equação cinética de Boltzmann em que somente são mantidos os detalhes moleculares essenciais para recuperar o comportamento macroscópico do escoamento.

Derivadas do sistema de Navier-Stokes <span id='citeF-7'></span>[[#cite-7|[7]]], as equações de águas rasas modelam escoamentos muito comuns na natureza e importantes, como por exemplo, canais, rios e estuários. Canais de navegação, de irrigação, de drenagem e canal da piracema, são exemplos de canais artificiais que ajudam no desenvolvimentos das cidades e da agricultura. Objeto de estudo deste trabalho, o canal da piracema, construído no parque da usina hidrelétrica de Itaipu, possui grande importância para a região. Pode-se perceber os benefícios que esta obra trouxe, permitindo que diversas espécies de peixes existentes na Costa Oeste do Estado, possam transpor a barragem da usina e, assim, voltam a utilizar sua antiga rota migratória para a desova <span id='citeF-8'></span>[[#cite-8|[8]]].

Devido a natureza do problema de estudo deste trabalho e também ao fato que o operador BGK ser menos estável <span id='citeF-9'></span>[[#cite-9|[9]]], optou-se por utilizar o método do reticulado de Boltzmann com múltiplos tempos de relaxamento (LBM-MRT). A capacidade do LBM-MRT foi testada em simulações de escoamento em canais abertos governados pelas equações de águas rasas e apresentaram resultados satisfatórios <span id='citeF-10'></span><span id='citeF-6'></span><span id='citeF-1'></span><span id='citeF-11'></span>[[#cite-10|[10,6,1,11]]]. Além disso, incorporou-se o fenômeno da turbulência, por meio da abordagem LES e do modelo de Smagorinsky. De acordo com Zhou <span id='citeF-7'></span>[[#cite-7|[7]]], a incorporação da turbulência fornece soluções mais próximas da situação real. Utilizou-se o reticulado quadrado <math display="inline">D2Q9</math>, considerado uma boa opção para simulação em canais <span id='citeF-12'></span>[[#cite-12|[12]]]. As condições de contorno ''Bounce-back'' para as paredes do canal, vazão e altura constantes na entrada e saída do canal, respectivamente,  também são consideradas levando em conta as características do canal <span id='citeF-13'></span>[[#cite-13|[13]]].

O estudo de caso abordado neste trabalho e simulado por meio do LBM, se ratifica por dois principais motivos. O primeiro deles está no fato de que o LBM-MRT é uma alternativa potencial para simular escoamento de água em canal, com fundo descontínuo, governado pelas equações de águas rasas. Com código de simples implementação e a possibilidade de utilizar computação paralela, o LBM mostra-se um método promissor na dinâmica de fluidos computacional, respeitando as suas restrições. O outro motivo está relacionado com a importância do canal da piracema, permitindo que diversas espécies de peixes possam retomar a rota migratória que utilizavam antes da construção da usina hidrelétrica. Além disso, a contribuição principal deste trabalho está no fato de que não exite na literatura simulação numérica bidimensional do problema proposto.

O trabalho está organizado conforme segue. Na seção [[#3 Equações de Águas Rasas|3]] têm-se as equações governantes do escoamento. Em seguida, na seção [[#4 Método do Reticulado de Boltzmann|4]] apresenta-se o método numérico que será usado na simulação do estudo de caso que é apresentado na seção [[#5 Estudo de caso|5]]. As simulações numéricas e discussões foram detalhadas na seção [[#7 Resultados e Discussões|7]], e, em seguida, na seção [[#8 Conclusão|8]], apresentam-se as conclusões.

==2. Equações de águas rasas==

Neste trabalho, têm-se um canal em que a água está em contato com a atmosfera e a profundidade da água é bem menor do que a largura do canal. Este escoamento é governado pelas equações de águas rasas ([[#eq-1|1]]) e são apresentadas aqui com a notação tensorial <span id='citeF-7'></span>[[#cite-7|[7]]],

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle\frac{{\partial h}}{{\partial t}} + \frac{{\partial \left({h{u_j}} \right)}}{{\partial {x_j}}} = 0\\ \displaystyle\frac{{\partial \left({h{u_i}} \right)}}{{\partial t}} + \frac{{\partial \left({h{u_i}{u_j}} \right)}}{{\partial {x_j}}} =  - g\frac{\partial }{{\partial {x_i}}}\left({\frac{{{h^2}}}{2}} \right)+ \nu \frac{{{\partial ^2}\left({h{u_i}} \right)}}{{\partial {x_j}\partial {x_j}}} + {F_i}  \end{array},    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

onde <math display="inline">t</math> é o tempo, <math display="inline">h</math> é a profundidade da água, <math display="inline">x_j</math> são as coordenadas espaciais <math display="inline">x</math> e <math display="inline">y</math>, <math display="inline">u_i</math> são as componentes do vetor velocidade nas direções <math display="inline">x</math> e <math display="inline">y</math>, isto é, <math display="inline"> u</math> e <math display="inline">v</math>, <math display="inline">g</math> é a aceleração da gravidade, <math display="inline">\nu </math> é a viscosidade cinemática e <math display="inline">{F_i}</math> é o termo de força dado por,

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{F_i} =  - gh\frac{{\partial {z_b}}}{{\partial {x_i}}} - \frac{{{\tau _{bi}}}}{\rho },    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

em que <math display="inline">{\partial {z_b}}/{\partial {x_i}}</math> é a declividade do fundo, <math display="inline">\rho </math> é a densidade e <math display="inline">\tau _{bi}</math> é a tensão de cisalhamento com o fundo do canal, definida por,

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\tau _{bi}} = \rho {C_b}{u_i}\sqrt {{u_j}{u_j}}   ,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

onde <math display="inline">C_b</math> é o coeficiente de atrito do leito do canal, que pode ser um valor constante ou estimado por <math display="inline">{C_b} = g /C_z^2</math>, em que <math display="inline">C_z</math> é o coeficiente de Chezy.

==3. Método do reticulado de Boltzmann==

A equação do reticulado de Boltzmann (LBE) ([[#eq-4|4]]), com o operador de colisão BGK é,

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f_\alpha }\left({\displaystyle x + {{\displaystyle e}_\alpha }\Delta t,t + \Delta t} \right)- {f_\alpha }\left({\displaystyle x,t} \right)= - \frac{1}{\tau }\left({{f_\alpha } - f_\alpha ^{eq}} \right)+ \frac{{\Delta t}}{{{N_\alpha }{e^2}}}{e_{\alpha i}}{F_i}\left({\displaystyle x,t} \right),    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

onde <math display="inline">e = {{\Delta x} \mathord{\left/{{{\Delta x} {\Delta t}}} \right.} {\Delta t}}</math> é a velocidade no reticulado, <math display="inline">\Delta x = \Delta y</math> é o espaçamento da malha, <math display="inline">\Delta t</math> é o incremento no tempo, <math display="inline">\displaystyle e_{\alpha }</math> são as direções de movimento das partículas, <math display="inline">f_\alpha </math> é a função distribuição de partículas, <math display="inline">F_i</math> é a componente da força na direção <math display="inline">i</math>, <math display="inline">N_\alpha = \frac{1}{{{e^2}}}\sum \limits  {{e_{\alpha i}}{e_{\alpha i}}}</math> é uma constante estabelecida de acordo com o reticulado, <math display="inline">\tau </math> é o parâmetro de relaxamento e <math display="inline">f_\alpha ^{eq}</math> é a função distribuição de equilíbrio.

Utilizar-se-a o reticulado <math display="inline">D2Q9</math>, conforme figura [[#img-1|1]], cujas direções do movimento são <span id='citeF-12'></span><span id='citeF-7'></span>[[#cite-12|[12,7]]],

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\displaystyle e_\alpha } = \left\{{\begin{array}{l}   {\left({0,0} \right)}&{{{se }}\,\,\,\alpha  = 0}\\ \displaystyle {e \left({\cos \frac{{(\alpha  - 1)\pi }}{4},\, sen\frac{{(\alpha  - 1)\pi }}{4}} \right)}&{{{se }}\,\,\,\alpha  = 1,3,5,7}\\  \displaystyle{e\sqrt 2 \left({\cos \frac{{(\alpha  - 1)\pi }}{4},\, sen\frac{{(\alpha  - 1)\pi }}{4}} \right)}&{{{se }}\,\,\,\alpha  = 2,4,6,8}   \end{array}} \right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-D2Q9_new.png|120px|Reticulado D2Q9]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Reticulado <math>D2Q9</math>
|}

A função distribuição de equilíbrio ([[#eq-6|6]]) é determinada de forma que se recuperem as equações de águas rasas. Este processo é feito por meio da análise de Chapman-Enskog que conecta a LBE com as equações macroscópicas governantes <span id='citeF-7'></span><span id='citeF-13'></span>[[#cite-7|[7,13]]],

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_\alpha ^{eq} = \left\{{\begin{array}{l}   \displaystyle {h - \frac{{5g{h^2}}}{{6{e^2}}} - \frac{{2h}}{{3{e^2}}}{u_i}{u_i}, \,\,\,\,\, \alpha  = 0}\\  \begin{array}{l} \displaystyle   \frac{{g{h^2}}}{{6{e^2}}} + \frac{h}{{3{e^2}}}{e_{\alpha i}}{u_i} + \frac{h}{{2{e^4}}}{e_{\alpha i}}{e_{\alpha j}}{u_i}{u_j} + \\ \displaystyle  - \frac{h}{{6{e^2}}}{u_i}{u_i},\,\,\,\,\,\,\,\, \alpha  = 1,\, 3,\, 5,\, 7   \end{array}\\  \begin{array}{l}   \displaystyle \frac{{g{h^2}}}{{24{e^2}}} + \frac{h}{{12{e^2}}}{e_{\alpha i}}{u_i} + \frac{h}{{8{e^4}}}{e_{\alpha i}}{e_{\alpha j}}{u_i}{u_j} + \\ \displaystyle  - \frac{h}{{24{e^2}}}{u_i}{u_i},\,\,\,\,\,\,\, \alpha  = 2,\, 4,\, 6,\, 8   \end{array}   \end{array}} \right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

As quantidades físicas de interesse, a profundidade ([[#eq-7|7]]) e a velocidade da água ([[#eq-8|8]]), são determinadas de forma a recuperar as equações de águas rasas,

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>h\left({\displaystyle {\bf x},t} \right)= \sum \limits _\alpha  {{f_\alpha }\left({\displaystyle {\bf x},t} \right)}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

e

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{u_i}\left({\displaystyle {\bf x},t} \right)= \frac{1}{{h\left({\displaystyle {\bf x},t} \right)}}\sum \limits _\alpha  {{e_{\alpha i}}{f_\alpha }\left({\displaystyle {\bf x},t} \right)}.    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Como o operador BGK é menos estável <span id='citeF-9'></span>[[#cite-9|[9]]], optou-se por utilizar o MRT. Proposto por d'Humieres <span id='citeF-14'></span>[[#cite-14|[14]]], O MRT captura as informações físicas do escoamento de forma mais completa e alcança estabilidade numérica de maneira mais eficiente do que o LBM-BGK. Foram Lallemand e Luo <span id='citeF-15'></span>[[#cite-15|[15]]] que fizeram uma análise sobre a estabilidade do LBM-MRT, mostrando como melhorar o modelo de forma satisfatória. A principal mudança entre o LBM-BGK e o LBM-MRT está no operador de colisão. Migra-se do espaço de velocidades, utilizado do LBM-BGK, para o espaço dos momentos, no LBM-MRT <span id='citeF-15'></span>[[#cite-15|[15]]]. A equação do reticulado de Boltzmann com múltiplos tempos de relaxamento é <span id='citeF-1'></span>[[#cite-1|[1]]],

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{f_\alpha }\left({\displaystyle {\bf x} + {{\displaystyle e}_\alpha }\Delta t,t + \Delta t} \right)- {f_\alpha }\left({\displaystyle {\bf x},t} \right)= - {\left({{M^{ - 1}}\hat S} \right)_{\alpha i}}\left({{m_{\alpha }} - m_{\alpha }^{eq}} \right)+ \frac{{\Delta t}}{{{N_\alpha }{e^2}}}{e_{\alpha i}}{F_i}\left({\displaystyle {\bf x},t} \right) ,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}
em que <math display="inline">M</math> é a matriz que transforma um vetor do espaço <math display="inline">{V} = {R}^{i}</math> , gerado pelas velocidades discretas, em um vetor do espaço dos momentos <math display="inline">M = R}{^{i}</math>, <math display="inline">M^{-1}</math> é a matriz inversa de <math display="inline">M</math>, <math display="inline">{{m_i{\alpha } }}</math> são as funções de distribuição no espaço dos momentos, onde <math display="inline">\left|m \right\rangle  = M\left|f \right\rangle </math> e <math display="inline">\left|. \right\rangle</math> é uma notação para vetor coluna, <math display="inline">m_{\alpha }^{eq}</math> são as funções de distribuição de equilíbrio no espaço dos momentos, onde <math display="inline">\left|{{m^{eq}}} \right\rangle  = M\left|{{f^{eq}}} \right\rangle </math> e <math display="inline">\widehat S</math> é uma matriz diagonal <math display="inline">\widehat S = diag\left({{s_0},\,{s_1},\,...,\,{s_{{\alpha } - 1}}} \right)</math>, onde os seus valores são definidos de acordo com o problema estudado.

Pode-se separar a equação ([[#eq-9|9]]) nas etapas de colisão,

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline f_\alpha  \left( {\bf x}, t\right)= - {\left({{M^{ - 1}}\widehat S} \right)_{ \alpha i}}\left({{m_{\alpha } } - m_{\alpha } ^{eq}} \right)+ \frac{{\Delta t}}{{6{e^2}}}{e_{\alpha i}}{{\bf F}_i},      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

com <math display="inline">\alpha  = 1,...,\,M</math>, e a propagação,

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> f_{\alpha }\left({\bf x}  + {\bf e}_{\alpha } \Delta x,\,\,t + \Delta t \right)= \overline f_{\alpha } \left({\bf x} ,\,\,t \right),.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

O termo de força em ([[#eq-10|10]]), <math display="inline">\frac{\Delta t}{N_\alpha e^2} e_{\alpha i}{F_i}\left({\bf x},t \right)</math> está sendo somado ao operador de colisão, ou seja, ele não faz parte do operador de colisão. Deste modo, o termo de força pode ser posicionado tanto somando juntamente com o operador de colisão, como na etapa de transmissão.

Considerando-se o reticulado <math display="inline">D2Q9</math>, são escolhidos momentos para cada uma das 9 velocidades <span id='citeF-15'></span>[[#cite-15|[15]]]. Estes momentos formam a matriz <math display="inline">M_{9\times{9}}</math> em ([[#eq-12|12]]), que é a matriz transformação do espaço das velocidades para o espaço dos momentos,

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>M = \left({\begin{array}{c}   1&1&1&1&1&1&1&1&1\\  { - 4}&{ - 1}&2&{ - 1}&2&{ - 1}&2&{ - 1}&2\\  4&{ - 2}&1&{ - 2}&1&{ - 2}&1&{ - 2}&1\\  0&1&1&0&{ - 1}&{ - 1}&{ - 1}&0&1\\  0&{ - 2}&1&0&{ - 1}&2&{ - 1}&0&1\\  0&0&1&1&1&0&{ - 1}&{ - 1}&{ - 1}\\  0&0&1&{ - 2}&1&0&{ - 1}&2&{ - 1}\\  0&1&0&{ - 1}&0&1&0&{ - 1}&0\\  0&0&1&0&{ - 1}&0&1&0&{ - 1}   \end{array}} \right).    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

Assim, como <math display="inline">\left|m \right\rangle  = M\left|f \right\rangle </math>, pode-se caracterizar as funções distribuição de equilíbrio do espaço dos momentos (<math display="inline">m_{\alpha }^{eq}</math>) da seguinte forma,
<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>m_{(1...9)}^{eq} = \left({\begin{array}{c}   h\\  \displaystyle { - 4h + \frac{{3g{h^2}}}{{{e^2}}} + \frac{{3h\left({{u^2} + {v^2}} \right)}}{{{e^2}}}}\\ \displaystyle  {4h - \frac{{9g{h^2}}}{{2{e^2}}} - \frac{{3h\left({{u^2} + {v^2}} \right)}}{{{e^2}}}}\\ \displaystyle  {\frac{{hu}}{e}}\\  \displaystyle { - \frac{{hu}}{e}}\\  \displaystyle {\frac{{hv}}{e}}\\ \displaystyle  { - \frac{{hv}}{e}}\\  \displaystyle {\frac{{h\left({{u^2} - {v^2}} \right)}}{{{e^2}}}}\\ \displaystyle  {\frac{{huv}}{{{e^2}}}}   \end{array}} \right).    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

Com relação a estabilidade linear do LBM, não é possível fazer sua caracterização completa <span id='citeF-16'></span>[[#cite-16|[16]]], porém Zhou <span id='citeF-7'></span>[[#cite-7|[7]]] estabelece indica algumas condições necessárias para se obter estabilidade na simulação. A respeito das condições para promover a estabilidade do LBM, um primeiro olhar deve ser dirigido para a viscosidade cinemática <math display="inline">\nu </math>. Ela deve ser sempre positiva em problemas de escoamento envolvendo fenômenos de difusão. É uma propriedade do fluido e relaciona-se com o parâmetro de relaxamento, <math display="inline">\tau </math>, presente no termo de colisão da LBE,

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{{\Delta t}}{{\Delta {x^2}}} = \frac{{2\tau  - 1}}{{6\nu }} ,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

com <math display="inline">\tau > 1/2</math> <span id='citeF-16'></span>[[#cite-16|[16]]].

Deve-se garantir que a magnitude da velocidade física resultante do fluido é menor do que a velocidade no reticulado, conforme ([[#eq-15|15]]) <span id='citeF-7'></span>[[#cite-7|[7]]]. O lado esquerdo de ([[#eq-15|15]]) é o denominado número de Courant <span id='citeF-16'></span>[[#cite-16|[16]]],

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{{{{\left|{\displaystyle {\bf u} } \right|}^2}}}{{{e^2}}} < 1,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

e também registra-se a condição ([[#eq-16|16]]) para a celeridade,

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{{gh}}{{{e^2}}} < 1.    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

onde <math display="inline">{\bf u} </math> é o vetor de velocidade física do fluido, <math display="inline">g</math> é a aceleração da gravidade e <math display="inline">h</math> é a altura da água.

Além das condições necessárias citadas acima, o LBM está limitado a fluxos de baixa velocidade e isto sugere que o LBM é adequado para escoamentos subcríticos de águas pouco profundas <span id='citeF-7'></span>[[#cite-7|[7]]]. Desta forma impõe-se que,

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{{{{\left|{\bf u } \right|}^2}}}{{gh}} < 1,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

As condições ([[#eq-15|15]]), ([[#eq-16|16]]) e ([[#eq-17|17]]) serão incorporadas ao código computacional e o teste será feito em todos os pontos da malha para cada passo de tempo, garantindo-se, deste forma, que a simulação seja estável.

Atualmente, a limitação ([[#eq-17|17]]) está em estudo, como pode ser visto nos trabalhos de Chopard, Pham e Lefèvre <span id='citeF-17'></span>[[#cite-17|[17]]] e <span id='citeF-18'></span>[[#cite-18|[18]]], onde o escoamento é supercrítico, e no trabalho de Rocca et al. <span id='citeF-19'></span>[[#cite-19|[19]]], onde o escoamento é transcrítico. As pesquisas indicam o uso de um reticulado assimétrico para tornar possível a abordagem de problemas envolvendo escoamento supercrítico, pois são apropriados para descrever fluxos com grandes velocidades em uma direção, além disso, são necessários mais estudos para a aplicação em simulações bidimensionais <span id='citeF-17'></span>[[#cite-17|[17]]].

Utiliza-se, com frequência, no LBM uma viscosidade mais elevada do que a viscosidade cinemática do fluido e os resultados estão em concordância com as simulações disponíveis na literatura <span id='citeF-7'></span><span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span>[[#cite-7|[7,20,21,22,23]]]. A justificativa para esta abordagem, de acordo com Rocca, Adduce e Sciortino <span id='citeF-21'></span>[[#cite-21|[21]]], é que o LBM é um modelo intrinsecamente viscoso e a viscosidade numérica pode ser adequadamente ajustada em concordância com os parâmetros da discretização, possibilitando evitar instabilidades numéricas. Em particular, Janben, Grilli e Krafczyk <span id='citeF-22'></span>[[#cite-22|[22]]], simularam escoamento de água com viscosidade no LBM várias vezes mais elevada do que a viscosidade cinemática. Simulações foram feitas no mesmo problema, com o objetivo de manter a estabilidade e aproximar a viscosidade do LBM com a viscosidade cinemática. Observou-se que uma redução ainda maior na viscosidade no LBM não altera significativamente os resultados.

A incorporação do fenômeno da turbulência é realizada na equação governante do LBM-MRT ([[#eq-9|9]]), em que o parâmetro de relaxamento <math display="inline">\tau </math> é substituído pelo parâmetro de relaxamento total <math display="inline">\tau _t=\tau + \tau _e</math> ([[#eq-18|18]]), onde <math display="inline">\tau _e</math> é o parâmetro de relaxamento turbulento. Utilizou-se a simulação de grandes escalas (LES), por meio do modelo de Smagorinsky. É uma técnica que oferece resultados refinados e utiliza menos recursos computacionais. A incorporação da turbulência apresentou resultados satisfatórios em diversos trabalhos <span id='citeF-7'></span><span id='citeF-13'></span>[[#cite-7|[7,13]]],

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\tau _t} = \tau  + {\tau _e} = \frac{\tau }{2} + \frac{1}{2}\sqrt {{\tau ^2} + \frac{{18C_s^2}}{{{e^2}h}}\sqrt {{\Pi _{ij}}{\Pi _{ij}}} },    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

onde <math display="inline">C_s</math> é a constante de Smagorinsky escolhida de forma coerente com o problema e <math display="inline">{\Pi _{ij}}</math> faz parte do tensor taxa de deformação nas simulações em grandes escalas,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\Pi _{ij}} = \sum \limits _\alpha  {{e_{\alpha i}}{e_{\alpha j}}\left({{f_\alpha } - f_\alpha ^{eq}} \right)}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

==4. Estudo de caso==

O problema que pretende-se simular é baseado no primeiro trecho do canal de alimentação em trincheira, conforme figura [[#img-2|2]], com comprimento de <math display="inline">30\,m</math> e largura de <math display="inline">8\,m</math>, o qual faz parte do sistema de transposição de peixes (canal da piracema) construído na usina hidrelétrica de Itaipu em 2002. Itaipu está localizada no rio Paraná, em um trecho de fronteira entre o Brasil e o Paraguai <span id='citeF-24'></span>[[#cite-24|[24]]]. O canal de alimentação em trincheira está localizado entre o dique de regulagem e o lago das Grevilhas. O dique de regulagem atua junto com a estrutura de tomada da água e controlam a vazão que escoa pelo canal. O lago das Grevilhas é uma lagoa artificial construída como local de descanso para os peixes <span id='citeF-8'></span>[[#cite-8|[8]]].

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-Itaipu_Piracema.png|288px|Canal da Piracema incluindo o rio Paraná, rio Bela Vista(RIBE), riacho Brasília, Canal de Deságüe do rio Bela Vista (CABV), Canal de Iniciação (CAIN), Lago Inferior (LAIN), Lago Principal (LAPR), Canal de Alimentação em Aterro (CAAT), Lago das Grevilhas (LAGR), Canal de Alimentação em Trincheiras (CATR), Dique de Regulagem (DIRE) e Reservatório (RESE).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-top:15px;| '''Figure 2:''' Canal da Piracema incluindo o rio Paraná, rio Bela Vista(RIBE), riacho Brasília, Canal de Deságüe do rio Bela Vista (CABV), Canal de Iniciação (CAIN), Lago Inferior (LAIN), Lago Principal (LAPR), Canal de Alimentação em Aterro (CAAT), Lago das Grevilhas (LAGR), Canal de Alimentação em Trincheiras (CATR), Dique de Regulagem (DIRE) e Reservatório (RESE).
|}

As dimensões e parâmetros reais do trecho escolhido, podem ser vistos no relatório técnico elaborado por Muller e Gibertoni <span id='citeF-25'></span>[[#cite-25|[25]]] e encontram-se listadas na tabela [[#table-1|1]]. As chicanas possuem as mesmas características do problema simulado, porém a altura real é igual a <math display="inline">0,6\, m</math>.

<div class="center" style="font-size: 75%;"> '''Tabela 1'''. Características reais no primeiro trecho do canal de alimentação em trincheira que compõe o Canal da Piracema na Usina Hidrelétrica de Itaipu.</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%" 
|-
|  style="border: 1pt solid black;text-align: center;width:85px;padding-left:5px;"|<span style="text-align: left; ">Extensão (m)</span>
|  style="border: 1pt solid black;text-align: center;width:75px;"|<span style="text-align: center; ">Decliv. (%)</span>
|  style="border: 1pt solid black;text-align: center;width:95px;"|<span style="text-align: center; ">Larg. Do fundo (m)</span>
|  style="border: 1pt solid black;text-align: center;width:125px;"|<span style="text-align: center; ">Distância entre as
chicanas (m)</span>
|  style="border: 1pt solid black;text-align: center;width:100px;"|<span style="text-align: center; ">Núm. de chicanas</span>
|-
|  style="border: 1pt solid black;padding-left:5px;text-align: center;"|<span style="text-align: center; ">24</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">8</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">6</span>
|-
|  style="border: 1pt solid black;padding-left:5px;text-align: center;"|<span style="text-align: center; ">6</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">0</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">8</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">-</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">-</span>
|}


Devido a problemas de instabilidade numérica, ocasionada principalmente pelo número de Froude maior do que 1, escolheu-se a caracterização para o trecho de acordo com a Tabela [[#table-2|2]].


<div class="center" style="font-size: 75%;"> '''Tabela 2'''. Características baseadas no primeiro trecho
do canal de alimentação em trincheira que compõe o Canal da Piracema na Usina Hidrelétrica de Itaipu.</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%" 
|-
|  style="border: 1pt solid black;text-align: center;width:85px;padding-left:5px;"|<span style="text-align: left; ">Extensão (m)</span>
|  style="border: 1pt solid black;text-align: center;width:75px;"|<span style="text-align: center; ">Decliv. (%)</span>
|  style="border: 1pt solid black;text-align: center;width:95px;"|<span style="text-align: center; ">Larg. Do fundo (m)</span>
|  style="border: 1pt solid black;text-align: center;width:125px;"|<span style="text-align: center; ">Distância entre as
chicanas (m)</span>
|  style="border: 1pt solid black;text-align: center;width:100px;"|<span style="text-align: center; ">Núm. de chicanas</span>
|-
|  style="border: 1pt solid black;padding-left:5px;text-align: center;"|<span style="text-align: center; ">24</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">1</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">8</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">6</span>
|-
|  style="border: 1pt solid black;padding-left:5px;text-align: center;"|<span style="text-align: center; ">6</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">0</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">8</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">-</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; ">-</span>
|}


Cada chicana possui <math display="inline">7\,m</math> de comprimento, largura da base igual a <math display="inline">1\,m</math>, largura do topo igual a <math display="inline">0,5\,m</math> e altura igual a <math display="inline">0,4\,m</math>. Na figura [[#img-3|3]] ilustra-se ver o trecho do canal da piracema que será simulado e na figura [[#img-4|4]], tem-se a representação computacional para o fundo do canal.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-Itaipu_CATrincheira.png|288px|Trecho do canal da piracema que será simulado.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Trecho do canal da piracema que será simulado.
|}

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-Fundo_Canal.png|288px|Representação computacional para o fundo do trecho escolhido.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Representação computacional para o fundo do trecho escolhido.
|}

A ocorrência da piracema está condicionada às seguintes restrições: velocidades médias não superiores a <math display="inline">3\,m/s</math> nas seções transversais ao longo do complexo do parque de barragem, profundidade da água superior a <math display="inline">0,8\,m</math> e áreas molhadas nas seções transversais superiores a <math display="inline">6\,m^2</math> (pode ser de <math display="inline">4\,m^2</math> em trechos curtos). O não atendimento aos requisitos citados somente serão permitidos em trechos bastante curtos e em condições especiais, como por exemplo, sobre os obstáculos (chicanas) construídas ao longo do canal e nas aberturas das chicanas <span id='citeF-25'></span>[[#cite-25|[25]]]. Simulações computacionais do escoamento unidimensional realizadas por Muller e Gibertoni <span id='citeF-25'></span>[[#cite-25|[25]]], concluíram que nas proximidades das chicanas ocorre a diminuição da profundidade da água acompanhada de um aumento de velocidade e no trecho onde o canal não apresenta obstáculo, o fluxo tende a voltar a condição anterior acontecendo um aumento de profundidade e consequentemente diminuição de velocidade. Além disso, a vazão ideal no canal é de <math display="inline">Q = 11,5 \,m^3/s</math>.

==5. Algoritmo para o LBM-MRT==

O algoritmo para o programa desenvolvido com base no LBM-MRT possui todas as etapas do algoritmo para o LBM-BGK. Porém, são necessárias mais etapas no LBM-MRT, pois o operador de colisão MRT necessita de mais cálculos. Portanto, a principal diferença entre os dois algoritmos está na abordagem do operador de colisão.

Na equação governante do LBM-BGK ([[#eq-4|4]]), tanto no lado esquerdo como no direito, têm-se a mesma função distribuição de partículas <math display="inline">f_{\alpha }</math>, isto porque todo o processo de simulação do LBM-BGK ocorre no espaço das velocidades. Observando-se agora a equação governante do LBM-MRT ([[#eq-9|9]]), percebe-se que as funções (<math display="inline">{m_{\alpha }}</math> e <math display="inline">m_{\alpha } ^{eq}</math>) que compõem o termo de colisão não são as mesmas que estão do lado esquerdo. Isto ocorre pelo fato de que <math display="inline">{m_{\alpha }}</math> e <math display="inline">m_{\alpha } ^{eq}</math> são funções distribuição que estão no espaço dos momentos.

A transformação do espaço de velocidades para o espaço dos momentos é feita por meio da matriz <math display="inline">M</math> apresentada em ([[#eq-12|12]]). Além disso, pode-se proceder com a obtenção de funções que estão no espaço de velocidades a partir das funções que estão no espaço dos momentos e isto pode feito conforme apresentado na seção [[#4 Método do Reticulado de Boltzmann|4]] e que será descrito a seguir com mais detalhes,

* Escrever a função distribuição de partículas, que está no espaço de velocidades, no espaço dos momentos: <math display="inline">\left|m \right\rangle  = M\left|f \right\rangle </math>;

* Escrever a função distribuição de equilíbrio, que está no espaço de velocidades, no espaço dos momentos: <math display="inline">\left|{{m^{eq}}} \right\rangle  = M\left|{{f^{eq}}} \right\rangle </math>;

* Escrever a função distribuição de partículas, que está no espaço dos momentos, no espaço gerado pelas velocidades: <math display="inline">\left|f \right\rangle  = {M^{ - 1}}\left|m \right\rangle </math>;

* Escrever a função distribuição de equilíbrio, que está no espaço dos momentos, no espaço gerado pelas velocidades: <math display="inline">\left|f^{eq} \right\rangle  = {M^{ - 1}}\left|{m^{eq}} \right\rangle </math><br/>


Portanto, somente o operador de colisão do LBM-MRT opera no espaço dos momentos. Neste sentido, observa-se que na etapa do operador de colisão, o código recebe <math display="inline">f_{\alpha }</math> e <math display="inline">f_{\alpha }^{eq}</math>, mas precisa transformá-las em <math display="inline">{m_{\alpha }}</math> e <math display="inline">m_{\alpha } ^{eq}</math>. Procedimento análogo acorre na saída do operador de colisão, têm-se <math display="inline">{m_{\alpha }}</math> e <math display="inline">m_{\alpha } ^{eq}</math> e é necessário obter <math display="inline">f_{\alpha }</math> e <math display="inline">f_{\alpha }^{eq}</math>.

O algoritmo para o programa desenvolvido utilizando o LBM-MRT é fornecido a seguir:

<ol>

<li>Escolha o reticulado. Identifique o tamanho da malha, o passo de temporal e o espacial. Estabeleça o valor para a constante de Smagorinsky, valor inicial para o parâmetro de relaxamento, valores da matriz diagonal <math display="inline">\widehat S</math> e caracterize a geometria do canal.    </li>

<li>Inicialização (<math display="inline">t=0</math>):    </li>

* Inicializa-se o valor da velocidade macroscópica nas direções <math display="inline">x</math> e <math display="inline">y</math>, <math display="inline">u \left({x, y, 0} \right)</math> e <math display="inline">v \left({x, y ,0} \right)</math>, respectivamente;

* Inicializa-se o valor da profundidade da água <math display="inline"> h \left({x, y, 0} \right)</math>;

* Calcular <math display="inline">f_{\alpha }^{eq}</math>, utilizando <math display="inline"> h \left({x, y, 0} \right)</math>, <math display="inline">u \left({x, y, 0} \right)</math> e <math display="inline">v \left({x, y ,0} \right)</math>, por meio de ([[#eq-6|6]]);

* Inicializar <math display="inline">{f_{\alpha }}\left({\displaystyle x ,0} \right)</math> utilizando <math display="inline">f_{\alpha }^{eq}\left({h , u, v } \right)</math>, ou seja <math display="inline">{f_{\alpha }}\left({\displaystyle x ,0} \right)=f_{\alpha }^{eq}\left({h , u, v } \right)</math>.  

<li>Loop principal (<math display="inline">t=1</math> até <math display="inline">t_{max}</math>):    </li>

* Informar quais condições de contorno serão utilizadas e incorporar no código;

* Atualizar as variáveis macroscópicas, profundidade e velocidades nas direções <math display="inline">x</math> e <math display="inline">y</math>, utilizando ([[#eq-7|7]]) e ([[#eq-8|8]]), respectivamente. Além disso, atualiza-se o parâmetro de relaxamento total por meio de ([[#eq-18|18]]);

* Calcular <math display="inline">m_{\alpha } ^{eq}</math> utilizando ([[#eq-13|13]]);

* Cálculo do termo de força utilizando ([[#eq-2|2]]);

* Cálculo da colisão usando a expressão ([[#eq-10|10]]). Neste momento ocorrem todas as etapas de transformação descritas nesta seção;

* Cálculo da propagação usando ([[#eq-11|11]]);

* Verificar as condições de estabilidade descritas em ([[#eq-15|15]]), ([[#eq-16|16]]) e ([[#eq-17|17]]);

* calcular o erro relativo entre a simulação no tempo atual e no tempo anterior.

* Escrever os resultados.

<li>Plotar os resultados.    </li>

</ol>

==6. Resultados e Discussões==

===6.1 Etapa I: Validação do modelo===

Antes de avaliar o problema proposto, procede-se com a validação do modelo desenvolvido por meio da simulação numérica do escoamento de água ao redor de um cilindro, localizado perpendicularmente ao fundo e no centro do canal. Este problema clássico da Engenharia Hidráulica é bastante utilizado na validação de simulações numéricas. Em Yulistiyanto <span id='citeF-26'></span>[[#cite-26|[26]]] encontra-se a comparação entre a solução numérica, por meio de diferenças finitas, e os dados experimentais. Em 2004, Zhou <span id='citeF-7'></span>[[#cite-7|[7]]] refez esta simulação utilizando o LBM com um tempo de relaxamento. Os parâmetros utilizados são encontrados em Yulistiyanto <span id='citeF-26'></span>[[#cite-26|[26]]], possibilitando a comparação com os dados experimentais. O comprimento do canal é <math display="inline">L_x = 4 \,m</math>, a largura é <math display="inline">L_y = 2 \,m</math>, a vazão na entrada é <math display="inline">Q = 0,248 \,m^3/s</math>, a altura da água na saída é <math display="inline">h = 0,185 \,m</math>, a inclinação do fundo do canal é <math display="inline">\partial z_b /\partial x =  - 6,25 \times {10^{ - 4}}</math> na direção do escoamento. O coeficiente de Manning é <math display="inline">n_b=0,012\,m/s^{1/3}</math>. A constante de Smagorinsky escolhida é <math display="inline">C_s=0,3</math> <span id='citeF-7'></span>[[#cite-7|[7]]], <math display="inline">{s_0} = ... = {s_6} = 1,0</math> e <math display="inline">{s_7} = {s_8} = 1 / \tau _t </math> <span id='citeF-1'></span>[[#cite-1|[1]]]. A simulação foi realizada utilizando-se passo espacial <math display="inline">\Delta x = \Delta y = 0,035\,m</math> , o temporal <math display="inline">\Delta t = 0,007\,m</math> e o parâmetro de relaxamento, <math display="inline">\tau = 0,60</math>. Com relação as condições de contorno, nas paredes do cilindro, é conveniente utilizar a condição de contorno ''bounce-back'' para representar velocidade nula e nas paredes do canal utilizou-se condição de escorregamento. Na entrada do canal utilizou-se vazão constante. Na saída do canal têm-se altura fixa <span id='citeF-7'></span><span id='citeF-13'></span>[[#cite-7|[7,13]]]. O número de iterações realizadas foi de 40000, o que corresponde a um tempo de 280 segundos de medição. Utilizou-se o compilador GNU Fortran 5.3.0 em um computador com processador intel i7, CPU 2,2 GHz, 16 GB de memória RAM e sistema operacional 64-bit Windows 10. O tempo de processamento gasto foi de 298 segundos.

Na figura [[#img-5|5]], tem-se o perfil da profundidade da água no meio do canal e na direção do escoamento com o parâmetro de relaxamento <math display="inline">\tau = 0,60</math>, e os dados experimentais. O erro relativo na simulação é menor do que <math display="inline">1%</math> no decorrer de todo o canal. Na figura [[#img-6|6]], é possível observar que a profundidade da água aumenta ao se aproximar do cilindro, diminui nas laterais e atinge seu valor mínimo logo depois do cilindro. Tem-se na figura [[#img-7|7]] a velocidade da água na direção do escoamento, nota-se que ela aumenta nas laterias do cilidro e atinge o valor mínimo logo depois do cilindro.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-CilindroCanal.png|328px|<span style="text-align: center; font-size: 75%;">Perfil da profundidade da água no centro do canal e na direção do escoamento.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' <span style="text-align: center; font-size: 75%;">Perfil da profundidade da água no centro do canal e na direção do escoamento.</span>
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-Profundidade.png|328px|<span style="text-align: center; font-size: 75%;">Solução numérica para a profundidade da água.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' <span style="text-align: center; font-size: 75%;">Solução numérica para a profundidade da água.</span>
|}

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-Velocidade.png|328px|<span style="text-align: center; font-size: 75%;">Solução numérica para a velocidade da água na direção u.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' <span style="text-align: center; font-size: 75%;">Solução numérica para a velocidade da água na direção <math>u</math>.</span>
|}

Observa-se a ótima concordância entre a simulação e os dados experimentais, demonstrando a capacidade, eficiência e precisão do modelo desenvolvido. Procede-se então, na simulação do problema principal, descrito na seção [[#5 Estudo de caso|5]], onde será usado o mesmo modelo desenvolvido. Devido às instabilidades numéricas ocasionadas principalmente pelo número de Froude ser maior do que um, não foi possível utilizar todos os dados reais na simulação e deste modo,  os parâmetros foram escolhidos de forma mais próxima possível ao modelo real, de modo a obter estabilidade numérica. Conforme apresentado anteriormente, mais estudos são necessários para simular escoamentos bidimensionais supercríticos no LBM.

===6.2 Etapa II: Simulação do estudo de caso proposto===

Os parâmetros reais do problema proposto na seção [[#5 Estudo de caso|5]] são, o comprimento do canal é <math display="inline">L_x = 30 \,m</math>, a largura do canal é <math display="inline">L_y = 8 \,m</math>, a altura inicial da água é <math display="inline">h_0 = 1,17 \,m</math>, a velocidade inicial na direção <math display="inline">x</math> é <math display="inline">u_0 = 1,2 \,m/s</math> e a velocidade inicial na direção <math display="inline">y</math> é <math display="inline">v_0 = 0,0 \,m/s</math>.

Nesta simulação, os parâmetros utilizados são, o comprimento do canal é <math display="inline">L_x = 30 \,m</math>, a largura do canal é <math display="inline">L_y = 8 \,m</math>, a altura inicial da água é <math display="inline">h_0 = 1,40 \,m</math>, a velocidade inicial na direção <math display="inline">x</math> é <math display="inline">u_0 = 1,0 \,m/s</math>, a velocidade inicial na direção <math display="inline">y</math> é <math display="inline">v_0 = 0,0 \,m/s</math> e o parâmetro de relaxamento é <math display="inline">\tau = 1.1</math>. Os parâmetros foram escolhidos de forma mais próxima possível ao modelo real. O parâmetros de relaxamento foi escolhido o mais próximo possível do valor mínimo limite, que é <math display="inline">\tau > 0,5</math>, de forma que a simulação retorne resultado estáveis. A constante de Smagorinsky escolhida é <math display="inline">C_s=0,3</math> <span id='citeF-27'></span>[[#cite-27|[27]]], <math display="inline">{s_0} = ... = {s_6} = 1,0</math> e <math display="inline">{s_7} = {s_8} = {1 \mathord{\left/{{1 \tau }} \right. } \tau }</math> <span id='citeF-1'></span>[[#cite-1|[1]]].

Com relação as condições de contorno, nas paredes do canal e nas chicanas, é conveniente utilizar a condição de contorno ''bounce-back'' para representar velocidade nula <span id='citeF-7'></span><span id='citeF-13'></span>[[#cite-7|[7,13]]]. Na entrada do canal utilizou-se vazão constante. Na saída do canal tem-se altura constante. Além disso, é necessário um tratamento para cada um dos quatro pontos localizados nos cantos da malha. Estas configurações podem ser vistas em Liu <span id='citeF-13'></span>[[#cite-13|[13]]].

Foram feitas simulações conforme apresentadas na tabela [[#table-3|3]], porém os resultados apresentados nas próximas figuras são para o caso em que a velocidade no reticulado é de <math display="inline">e = 5,7 \,m/s</math>, o espaçamento na malha é <math display="inline">\Delta x = \Delta y = 0,04 \, m</math> e o passo de tempo <math display="inline">\Delta t = 0,007 \, s</math>. O número de iterações realizadas foi de <math display="inline">50000</math>, o que corresponde a um tempo de <math display="inline">350</math> segundos de medição. Utilizou-se o compilador GNU Fortran 5.3.0 em um computador com processador intel i7, CPU 2,2 GHz, 16 GB de memória RAM e sistema operacional 64-bit Windows 10. O tempo de processamento gasto no processamento foi de <math display="inline">8988</math> segundos.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Parâmetros de simulação.
|- style="border-top: 2px solid;"
| style="border-right: 2px solid;" |    
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Parâmetros de simulação.
|-
| Veloc. no 
|}

|-
|

reticulado 


|-
|

(m/s) 


|


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|-
| Passo 
|}


|-
|

espacial 


|-
|

(m) 


|


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|-
| Passo de
|}


|-
|

tempo 


|-
|

(s) 


|


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|-
| Estabil.    
|}


|


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|-
| Erro 
|}


|-
|

relativo (<math display="inline">%</math>) 


|-
|

Máx. e Min.    


|- style="border-top: 2px solid;"
|

5,0 


|

0.04 


|

0,00800 


|

Instável 


|

- 


|-
|

5,1 


|

0.04 


|

0,00780 


|

Instável 


|

- 


|-
|

5,2 


|

0.04 


|

0,00770 


|

Instável 


|

- 


|-
|

5,3 


|

0.04 


|

0.00752 


|

Instável 


|

- 


|-
|

5,4 


|

0.04 


|

0,00741 


|

Instável 


|

- 


|-
|

5,5 


|

0.04 


|

0,00727 


|

Estável 


| <math>0,13</math>

- <math display="inline">0,01</math> 


|-
|

5,6 


|

0.04 


|

0,00714 


|

Estável 


| <math>0,13</math>

- <math display="inline">0,01</math> 


|-
|

5,7 


|

0.04 


|

0,0070 


|

Estável 


| <math>0,13</math>

- <math display="inline">0,01</math> 


|-
|

6 


|

0.04 


|

0,00670 


|

Estável 


| <math>0,14</math>

- <math display="inline">0,01</math> 


|-
|

6,5 


|

0.04 


|

0,00615 


|

Estável 


| <math>0,14</math>

- <math display="inline">0,01</math> 


|- style="border-bottom: 2px solid;"
|

7,0 


|

0.04 


|

0,00570 


|

Estável 


| <math>0,15</math>

- <math display="inline">0,01</math>



|}

A figura [[#img-8|8]] apresenta a profundidade da água ao longo do canal. Nota-se que a restrição de altura mínima <math display="inline">0,80 \, m</math> foi satisfeita. Verifica-se também que, ao se aproximar da chicana, a profundidade da água aumenta. Ligeiramente depois da chicana, a profundidade diminui. Em cima da chicana a lamina de água é menor do que no restante do canal. Além disso, observa-se que na abertura da chicana, a profundidade diminui e torna a subir quando se aproxima da próxima chicana.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-AlturaAgua.png|288px|<span style="text-align: center; font-size: 75%;">Profundidade da água no canal com chicanas.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' <span style="text-align: center; font-size: 75%;">Profundidade da água no canal com chicanas.</span>
|}

Nota-se que na figura [[#img-9|9]], somente nos pontos em cima e na abertura das chicanas, a velocidade na direção <math display="inline">x</math> atinge seu ponto máximo. Esta previsão é feita no trabalho de Muller e Gibertoni <span id='citeF-25'></span>[[#cite-25|[25]]]. Na maioria dos pontos do canal a velocidade é menor do que <math display="inline">1,2 \, m/s</math>. Este resultado mostra que, de fato, as chicanas exercem seu papel de conter a velocidade da água. Verifica-se também que, a velocidade máxima é menor do que <math display="inline">3 \,m/s</math>, que é a velocidade máxima permitida para a ocorrência da piracema. Esta velocidade só pode ser excedida, na abertura ou em cima da chicana <span id='citeF-25'></span>[[#cite-25|[25]]]. A figura [[#img-10|10]] mostra que, utilizando os mesmos parâmetros de simulação, porém sem as chicanas, a velocidade na direção <math display="inline">x</math>, é máxima em uma grande região do canal. A comparação entre as figura [[#img-9|9]] e [[#img-10|10]], ratifica a importância das chicanas.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-VelocidadeDirecaoX.png|288px|<span style="text-align: center; font-size: 75%;">Velocidade da água na direção x no canal com chicana.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' <span style="text-align: center; font-size: 75%;">Velocidade da água na direção <math>x</math> no canal com chicana.</span>
|}

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-VelocidadeDirecaoX_SChic.png|288px|<span style="text-align: center; font-size: 75%;">Velocidade da água na direção x no canal sem chicana.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' <span style="text-align: center; font-size: 75%;">Velocidade da água na direção <math>x</math> no canal sem chicana.</span>
|}

Tem-se também, a velocidade da água na direção <math display="inline">y</math>, conforme figura [[#img-11|11]]. Observa-se a simetria da simulação e percebe-se que nas aberturas das chicanas, a velocidade é mais representativa. Na figura [[#img-12|12]], os vórtices acorrem com maior intensidade nas paredes do canal próximo da abertura das chicanas.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-VelocidadeDirecaoY.png|288px|<span style="text-align: center; font-size: 75%;">Velocidade da água na direção y no canal com chicana.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' <span style="text-align: center; font-size: 75%;">Velocidade da água na direção <math>y</math> no canal com chicana.</span>
|}

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_481518576-MagnitudeVorticidade.png|288px|<span style="text-align: center; font-size: 75%;">Magnitude dos vórtices no canal com chicana.</span>]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' <span style="text-align: center; font-size: 75%;">Magnitude dos vórtices no canal com chicana.</span>
|}

==7. Conclusão==

Considerando que o estudo de caso apresentado neste trabalho é um escoamento complexo, os resultados apresentados mostram a capacidade e a precisão do LBM na simulação de escoamento em canais com obstáculos submersos, governados pelas equações de águas rasas. Observa-se que a simulação cumpriu as condições para a ocorrência da piracema descritas na seção [[#5 Estudo de caso|5]], ou seja, as velocidades não ultrapassaram <math display="inline">3\,m/s</math> nas seções transversais ao longo do trecho analisado, a profundidade da água é superior a <math display="inline">0,8\,m</math> e áreas molhadas nas seções transversais superiores a <math display="inline">6\,m^2</math>. Além disso, conforme mencionado na seção [[#5 Estudo de caso|5]], ao se aproximar da chicana, a profundidade da água aumenta. Em cima da chicana e na abertura da chicana, a profundidade diminui e a velocidade aumenta. Nota-se também que a profundidade da água logo depois da chicana é menor que antes da chicana. Este processo se repete ao longo do trecho analisado. Outro ponto interessante é que o LBM-MRT foi capaz de confirmar que as chicanas exercem o papel de reduzir a velocidade da água.

O modelo de Smagorinsky e o MRT permitiram a diminuição do parâmetro de relaxamento e consequentemente, os resultados descrevem a situação física de forma mais realista. Sem esta implementação, utilizando o LBM-BGK, o valor mínimo alcançado para o parâmetro de relaxamento foi <math display="inline">1,5</math>, enquanto que no LBM-MRT <math display="inline">\tau=1,1</math>. Além disso, a presença no código da rotina de estabilidade trouxe, juntamente com a consistência do método, a garantia da convergência da solução numérica. Observou-se que o erro relativo é maior apenas nos dois cantos da entrada do canal. Mesmo com o tratamento aplicado nestes pontos, o erro relativo nesta região é <math display="inline">10</math> vezes maior do que o erro ao longo do canal. Pretende-se desenvolver mais estudos e melhor abordar a condição de contorno na entrada do canal.

Neste estudo, o escoamento é permanente. Como trabalhos futuros, pretende-se repetir as simulações para o caso de vazão variável na entrada do canal, onde tem-se um escoamento não-permanente. Além disso, a presença das chicanas com abertura lateral no canal da piracema, faz com a que água descreva um caminho que se assemelha a meadros. Pretende-se então incorporar a pressão não-hidrostática na simulação e como isto, representar de maneira mais real o escoamento. Outros pontos de interesse são paralelizar o código LBM-MRT e compará-lo com o LBM-BGK e, adaptá-lo para que seja possível simular escoamento supercrítico bidimensional.

==Bibliography==

<div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

<div id="cite-1"></div>
[[#citeF-1|[1]]]  Y. Peng, J. G. Zhou, J. M. Zhang, H. Liu, Lattice Boltzmann modeling of shallow water flows over discontinuous beds, International Journal for Numerical Methods in Fluids 75 (2014) 608-619.

<div id="cite-2"></div>
[[#citeF-2|[2]]]  E. Bladéa, L. Cea, G. Coresteina, E. Escolanoc, J. Puertas,  E. Vázquez-Cendónd, J. Dolz, A. Coll, Iber: herramienta de simulación numérica del flujo en ríos, Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería 30 (1) (2014) 1-10.

<div id="cite-3"></div>
[[#citeF-3|[3]]]  F. U. Prieto, J. J. B. Muñoz, L. G. Corvinos, E. S. Casino, A. C. Acevedo, Estudio de la estabilidad y dispersión del problema de propagación de ondas sísmicas en 2-d utilizando el método de diferencias finitas generalizadas, Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería 27 (4) (2011) 269-277.

<div id="cite-4"></div>
[[#citeF-4|[4]]]  A. M. Lobeiro, Solução das Equações de Saint-Venant em Uma e Duas Dimensões Usando o Método das Características, Ph.D. thesis, Universidade Federal do Paraná, Curitiba, Brasil (2012).

<div id="cite-5"></div>
[[#citeF-5|[5]]]  M. Dourado, J. Meireles, Obtenção de um modelo de elementos finitos simplificado para representação de juntas rebitadas em análise dinâmica de estruturas usando uma ferramenta de updating, Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, 32 (3) (2016) 173-180.

<div id="cite-6"></div>
[[#citeF-6|[6]]]  Y. Peng, Lattice Boltzmann Simulations of Environmental Flow Problems in Shallow Water Flows, Ph.D. thesis, University of Liverpool, 2012.

<div id="cite-7"></div>
[[#citeF-7|[7]]]  J. G. Zhou, Lattice Boltzmann Method for Shallow Water Flows, Springer, New York, 2004.

<div id="cite-8"></div>
[[#citeF-8|[8]]]  M. F. Friedrich, F. R. Terabe, J. J. Ota, N. do C. Saks, Levantamento de informações durante visitas a sistemas de transposição de peixes no Brasil, in: XVII Simpósio Brasileiro de Recursos Hídricos, São Paulo, Brasil, 2007.

<div id="cite-9"></div>
[[#citeF-9|[9]]]  Y. Peng, J. M. Zhang, J. G. Zhou, Lattice Boltzmann Model Using Two Relaxation Times for Shallow-Water Equations, Journal of Hydraulic Engineering, 142 (2) (2016) 1-6.

<div id="cite-10"></div>
[[#citeF-10|[10]]]  H. Liu, M. Li, A. Shu. Burrows, Large eddy simulation of turbulent shallow water flows using multi-relaxation-time lattice Boltzmann model, International journal for numerical methods in fluids 70 (2012) 1573-1589.

<div id="cite-11"></div>
[[#citeF-11|[11]]]  C. Z. Zhang, Y. G. Cheng, J. Y. Wu, Lattice Boltzmann simulation of the open channel flow connecting two cascaded hydropower stations, Journal of Hydrodynamics, 3 (28) (2016) 400-410.

<div id="cite-12"></div>
[[#citeF-12|[12]]]  S. Chen, G. D. Doolen. Lattice Boltzmann Method for fluid flows, Annual Review of Fluid Mechanics 30 (1) (1998) 329-364.

<div id="cite-13"></div>
[[#citeF-13|[13]]]  H. Liu, Lattice Boltzmann Simulations for Complex Shallow Water Flows. Ph.D. thesis, University of Liverpooly, Liverpooly, 2009.

<div id="cite-14"></div>
[[#citeF-14|[14]]]  D. d'Humieres, Generalized lattice Boltzmann equations. In Rarefied Gas Dynamics: Theory and Simulations, Shizgal BD, Weaver DP (eds), IEEE Latin America Transactions, 159 (1992) 450-458.

<div id="cite-15"></div>
[[#citeF-15|[15]]]  P. Lallemand, L-S. Luo, Theory of the Lattice Boltzmann Method: Dispersion, Dissipation, Isotropy, Galilean Invariance, and Stability, Physical Review E 61 (2000) 6546-6562.

<div id="cite-16"></div>
[[#citeF-16|[16]]]  J. Sterling, S. Chen, Stability Analysis of Lattice Boltzmann Methods, Journal of Computational Physics, 123 (1) (1996) 196-206.

<div id="cite-17"></div>
[[#citeF-17|[17]]]  B. Chopard, V. T. Pham, L. Lefèvre, Asymmetric lattice Boltzmann model for shallow water flows, Computers & Fluids 88 (2013) 225-231.

<div id="cite-18"></div>
[[#citeF-18|[18]]]  A. H. Hedjripour, D. P. Callaghan, T. E. Baldock, Generalized transformation of the lattice Boltzmann method for shallow water flows, Journal of Hydraulic Research 54 (4) (2016) 371-388.

<div id="cite-19"></div>
[[#citeF-19|[19]]]  M. La Rocca, A. Montessori, P. Prestininzi, S. Succi, A multispeed Discrete Boltzmann Model for Transcritical 2D Shallow Water Flows, Journal of Computational Physics 284 (2015) 117-132.

<div id="cite-20"></div>
[[#citeF-20|[20]]]  J. G. Zhou, H. Liu, S. Shafiai, Y. Peng, R. Burrows, Lattice Boltzmann Method for Open-channel Flows, Engineering and Computational Mechanics 163 (2010) 243-249.

<div id="cite-21"></div>
[[#citeF-21|[21]]]  M. La Rocca, C. Adduce, G. Sciortino, Development of a Lattice Boltzmann Method for Two-layered Shallow-water Flow, International Journal for Numerical Methods in Fluids 70 (2012) 1048-1072.

<div id="cite-22"></div>
[[#citeF-22|[22]]]  C. F. Janben, S. T. Grilli, M. Krafczyk, Efficient simulations of long wave propagation and runup using a LBM approach on GPGPU hardware, In: Proceedings of the Twenty-second International Offshore and Polar Engineering Conference, Rodes, Grécia, 2012.

<div id="cite-23"></div>
[[#citeF-23|[23]]]  J. G. Zhou, H. Liu, Determination of bed elevation in the enhanced lattice Boltzmann method for the shallow-water equations, Physical Review E 88 (2013) 1-6.

<div id="cite-24"></div>
[[#citeF-24|[24]]]  S. Makrakis, O canal da Piracema como sistema de transposição. Ph.D. thesis, UEM, 2007.

<div id="cite-25"></div>
[[#citeF-25|[25]]]  I. I. Muller, R. C. Gibertoni Projeto HG-169: Estudos computacionais do comportamento hidráulico dos canais do complexo do parque da barragem e do rio Bela Vista. Universidade Federal do Paraná. Relatória técnico, Centro de Hidráulica e Hidrologia Professor Parigot de Souza (CEHPAR), Curitiba, Brasil, 2002.

<div id="cite-26"></div>
[[#citeF-26|[26]]]  B. Yulistiyanto, Flow around a cylinder installed in a fixed-bed open channel. Ph.D. thesis, École Polytechnique Fédérale de Lausanne, Lausanne, 1997.

<div id="cite-27"></div>
[[#citeF-27|[27]]]  H. Liu, J. G. Zhou, R. Burrows, Lattice Boltzmann simulations of the transient shallow water flows, Advances in Water Resource 33 (2010) 387-396.
</div>