Simulação Numérica Direta de Turbulência com Corte e Estratificação Escalar em Malhas Não-Ortogonais

Simulação Numérica Direta de Turbulência com Corte e Estratificação Escalar em Malhas Não-Ortogonais D. F. Jorge 1, D. F. Jorge 2 1Faculdade de Engenharia da Universidade do Porto Universidade do Porto Rua Dr. Roberto Frias, s/n 4200-465 Porto e-mail: [email protected] 2Escola Superior de Ciência e Tecnologia Instituto Superior Politécnico Gaya Av. dos Descobrimentos 333, 4400-103 Santa Marinha, Vila Nova de Gaia e-mail: [email protected] Resumo Neste trabalho propomos a discretização dos termos de corte e da estratificação escalar, a partir das equações de Navier-Stokes e transporte escalar, utilizando o método dos volumes finitos em coordenadas não-ortogonais no contexto da Simulação Numérica Directa (DNS) na simulação de escoamentos com corte e estratificação escalar. As técnicas numéricas utilizadas foram validadas através da simulação de um escoamento turbulento com corte e estratificação escalar num domínio cúbico periódico. Nas simulações foram impostas condições de fronteira de corte periódico na direcção vertical z e condições periódicas nas restantes direcções. As simulações iniciam-se com campos isotrópicos, com distribuições Gaussianas e variância igual para a velocidade e o escalar, utilizando um número de Reynolds Reλ igual a 40, número de Prandlt P r igual a 0 . 7, e malhas com 64 3 nós. A boa concordância entre os resultados obtidos com malhas ortogonais e não ortogonais, mesmo sob diferentes graus de distorção, permite concluir que as técnicas numéricas aplicadas são adequadas. Keywords: Simulação Numérica Directa, Métodos Numéricos, Condições de Entrada, Corte, Turbulência, Dinâmica de Fluidos. 1 1 INTRODUÇÃO A análise de escoamentos turbulentos tem sido realizada com base em simulações computacionais. Esta abordagem tem ganho relevância em relação a outras técnicas, nomeadamente as experimentais, por ser mais rápida, económica e por permitir uma análise mais detalhada dos escoamentos. Os programas utilizados nestas simulações baseiam-se nas equações de Navier-Stokes, na equação da continuidade e na equação de transporte escalar, que são resolvidas numericamente. A maioria dos escoamentos de interesse para a engenharia ocorrem em geometrias complexas e são turbulentos, o que torna o estudo destes escoamentos desafiante e exige métodos de simulação adequados. Neste trabalho, utilizámos o método dos volumes finitos em coordenadas não ortogonais para discretizar os termos de corte e de estratificação escalar, com base nas equações de transporte de quantidade de movimento e escalar no contexto da Simulação Numérica Directa (DNS). O uso de condições de fronteira periódicas é habitual em DNS. Contudo, na presença de corte, um campo inicialmente periódico perde periodicidade ao longo do tempo. Por exemplo, duas partículas de fluido colocadas em diferentes alturas z deslocam-se a velocidades médias horizontais U distintas e acabem por afastar-se gradualmente. Isso implica desafios adicionais na aplicação de condições de fronteira. Em Rogallo [ ? ], o problema das condições de fronteira é resolvido através de uma transformação de coordenadas dependente do tempo (Lagrange), que impõe periodicidade na direção transformada. Esse procedimento foi aplicado em diversos estudos do grupo de Stanford (p.ex., Rogers and Moin [ ? ]). No entanto, a necessidade de remeshing com frequência 1 2dU/dz introduz erros de interpolação. Utilizámos uma aproximação alternativa proposta por Baron [ ? ] para fluxo com corte puro, e também empregada por Gerz et al. [ ? ]. Nessa abordagem, as equações são discretizadas na formulação Euleriana com condições de fronteira de corte periódico, cuja periodicidade evolui no tempo conforme a equação (30). As técnicas numéricas de discretização foram validadas por meio de simulações de escoamentos turbulentos com corte e estratificação escalar num domínio cúbico periódico submetido a diferentes distorções de malha. A boa concordância observada em todas as malhas reforça a eficácia dos métodos aplicados aos termos de corte e estratificação escalar. 2 DESCRIÇÃO MATEMÁTICA 2.1 Determinação dos Termos de Corte Vamos analisar o termo não linear das equações de Navier-Stokes,                          ∂(vxvx) ∂x +∂(vxvy) ∂y +∂(vxvz) ∂z ∂(vyvx) ∂x +∂(vyvy) ∂y +∂(vyvz) ∂z ∂(vzvx) ∂x +∂(vzvy) ∂y +∂(vzvz) ∂z .(1) 2 Admitimos que a velocidade instantânea vi se decompõe nas componentes média Ui e flutuação ui, ou seja,        vx=Ux+ux, vy=Uy+uy, vz=Uz+uz. (2) Substituindo a equação 2 na equação 1, obtemos um sistema de três componentes para o termo convectivo. Temos a componente x,                      ∂ ∂x(U+u)(U+u) = ∂ ∂x(U·U |{z} 0 +Uu +uU +uu) = ∂ ∂x2uU +∂ ∂xuu ∂ ∂y(U+u)(V+v) = ∂ ∂y(U·V |{z} 0 +Uv +uV |{z} 0 +uv) = ∂ ∂yvU +∂ ∂yuv ∂ ∂z (U+u)(W+w) = ∂ ∂z (U·W | {z } 0 +Uw +uW |{z} 0 +uw) = ∂ ∂z wU +∂ ∂z uw ,(3) a componente y,                      ∂ ∂x(V+v)(U+u) = ∂ ∂x(V·U |{z} 0 +V u |{z} 0 +vU +vu) = ∂ ∂xvU +∂ ∂xvu ∂ ∂y(V+v)(V+v) = ∂ ∂y(V·V |{z} 0 +V v |{z} 0 +vV |{z} 0 +vv) = ∂ ∂yvv ∂ ∂z (V+v)(W+w) = ∂ ∂z (V·W | {z } 0 +V w |{z} 0 +vW |{z} 0 +vw) = ∂ ∂z vw ,(4) e a componente z,                      ∂ ∂x(W+w)(U+u) = ∂ ∂x(W·U | {z } 0 +Wu +wU +wu) = ∂ ∂xwU +∂ ∂xwu ∂ ∂y(W+w)(V+v) = ∂ ∂y(W·V | {z } 0 +Wv |{z} 0 +wV |{z} 0 +wv) = ∂ ∂ywv ∂ ∂z (W+w)(W+w) = ∂ ∂z (W·W | {z } 0 +Ww |{z} 0 +wW |{z} 0 +ww) = ∂ ∂z ww .(5) Nos sistemas de equações 3, 4 e 5 os termos sem componente média constituem o termo convectivo das equações de Navier-Stokes já conhecido. Para tratar os restantes termos vamos decompor o segundo termo da primeira equação do sistema de equações 3 da seguinte forma, ∂ ∂x2uU =∂ ∂xuU | {z } A +∂ ∂xuU | {z } B .(6) O termo (A) da equação 6 e aos termos dos sistemas de equações 4 e 5 constituem o primeiro termo de corte, U∂ui ∂x1 .(7) Adicionando o termo (B) da equação 6 aos restantes termos do sistema de equações 3 resulta, u∂U ∂x +U∂u ∂x +v∂U ∂y +U∂v ∂y +w∂U ∂z +U∂w ∂z .(8) 3 Considerando o gradiente de velocidade nulo nas direcções x e y e agrupando os termos que constituem a equação da continuidade, podemos obter o segundo termo de corte, u∂U ∂x |{z} 0 +v∂U ∂y |{z} 0 +U ∂u ∂x +∂v ∂y +∂w ∂z ! | {z } 0 +w∂U ∂z ,(9) De notar que este termo só existe na componente xdas equações de Navier-Stokes. Os termos das equações 7 e 9 têm de ser adimensionalizados. Para isso começamos por definir o número de corte adimensional, S=L ∆U dU dz⇔dU dz=S∆U L.(10) Integrando a equação 10, obtemos uma expressão para o cálculo da velocidade. Zx3 0 dU dx3 dx3=Zx3 0 S∆U Ldx3⇔U=S∆U Lx3(11) Substituindo a equação 11 na equação 7 e adimensionalizando a derivada da velocidade ∂ui/∂x1resulta, S∆U Lx3 | {z } U ∆U L ∂u∗ i ∂x∗ 1 =S∆U2 L x3 L |{z} x∗ 3 ∂u∗ i ∂x∗ 1 ,(12) Dividindo por ∆ U2/L , à semelhança de todos os termos da equação de quantidade de movimento, a equação 12 reduz-se à equação 13, que constitui o primeiro termo de corte adimensional. Sx∗ 3 ∂u∗ i ∂x∗ 1 (13) Procedendo de modo idêntico com o segundo termo de corte, equação 9, w∂U ∂z =w∗∆UdU dz,(14) e dividindo por ∆U2/L, chegamos ao segundo termo de corte adimensional. w∗∆UL ∆U2 dU dz=w∗L ∆U dU dz | {z } S =Sw∗=Su∗ 3(15) 2.2 Modelo Matemático As equações dimensionais que regem o transporte de escalares passivos na presença de corte e estratificação são a equação da continuidade, ∂ρui ∂xi = 0 ,(16) equação de Navier-Stokes, ∂ρui ∂t +∂(ρujui) ∂xj +Sx3 ∂ρui ∂x1 +Sρu3δi1=∂τij ∂xj −∂p ∂xi , i = 1,2,3,(17) 4 Figura 1: Perfis do campo médio de velocidade e escalar. e equação de transporte do escalar, ∂ρϕ ∂t +∂(ρujϕ) ∂xj +Sx3 ∂ρϕ ∂x1 +sρu3=∂ ∂xj Γ∂ϕ ∂xj!.(18) A componente média horizontal de velocidade e do escalar varia linearmente com a altura z , conforme ilustrado na figura 1. O fluido apresenta viscosidade ν e difusividade Γ constantes. Denotamos por ui , ϕ e p as flutuações em relação aos campos médios de velocidade, escalar e pressão, respectivamente. Introduzimos ainda o número de corte adimensional S=L ∆U dU dz , e o parâmetro de estratificação s=L ∆Φ dΦ dz , cuja escolha de valores S = 0 , 1e s = − 1 , 0 , 1permite distinguir os casos sem corte, com corte puro, e estratificação instável, neutra e estável, respectivamente. 2.3 Discretização dos Termos Convectivos de Corte O termo convectivo não linear das equações de Navier-Stokes, na presença de um gradiente médio de velocidade, gera dois novos termos na equação de transporte de quantidade de movimento: Sv,1 i=Z V ρSx3 ∂ui ∂x1 dV= (ρSx3β11)·(ui,e −ui,w),(19) Sv,2 i=Z V ρSu3δi1dV≈X j ρSu3,P δi1∆Vj.(20) 5 Os indices n, s, e, w, t, b denotam, respectivamente, as faces norte, sul, este, oeste, topo e base do volume de controle. Cada factor da equação 19 foi avaliado separadamente. Em particular, ρSx3β11 =ρS "(k−0.5Nk −0.5) L (Nk −2)#L2 (Nj −2)(Nk −2),2≤k≤Nk −1.(21) As componentes da velocidade na face Este e Oeste ( ui,e e ui,w ), são obtidas por interpolação linear nos nós P E eWde acordo com o sistema de equações (22), (ui,e =λeui,P +(1−λe)ui,E ui,w =λwui,P +(1−λw)ui,W ,(22) em que os coeficientes de interpolação linear λeeλwsão determinados por, λe=xE−xe xE−xP ;λw=xW−xw xW−xP ,(23) para malhas cartesianas e por, λe=∥xE−xe∥ ∥xE−xe∥+∥xe−xP∥;λw=∥xW−xw∥ ∥xW−xw∥+∥xw−xP∥,(24) para manhas não-ortogonais. Usando malhas regularmente espaçadas, λe=λw= 0.5e, (ui,e −ui,w) = 1 2(ui,E −ui,W ).(25) Conjugando as equações 21 e 25 resulta uma expressão final para o cálculo do primeiro termo de corte, Sv,1 i≈X j ρS 1 2(ui,E −ui,W )"(k−0.5Nk −0.5) L3 (Nj −2)(Nk −2)2#.(26) Tal como no transporte de quantidade de movimento, o uso de gradiente médio de velocidade e escalar vão dar origem a dois novos termos na equação de transporte do escalar. Um termo devido ao gradiente médio de velocidade, definido por, Sv=Z V ρSx3 ∂ϕ ∂x1 dV≈X j ρS 1 2(ϕE−ϕW)"(k−0.5Nk −0.5) L3 (Nj −2)(Nk −2)2#(27) e outro devido ao gradiente médio do escalar, sϕ=Z V ρsu3dV≈X j ρsu3,P ∆Vj.(28) 2.4 Domínio de Cálculo O domínio físico é um cubo de aresta π , discretizado por uma malha uniforme com 64 pontos cartesianamente espaçados. Para avaliar os efeitos da não-ortogonalidade da malha, realizamos simulações num domínio distorcido, em que as direcções computacionais ξ e η coincidem com as direcções 6 Figura 2: Domínio físico, inclinado a 30°entre as direcções zand ζ. físicas z e ζ , ver figura 2. Definindo α como o ângulo entre z e ξ , o tensor de transformação de coordenadas B/J é dado por: 1 ∆Bij, Bii = 1, B31 =−tan α, Bij = 0 (i=j), onde Jé o Jacobiano e ∆é o fator de escala da malha. A viscosidade cinemática foi ν = 0 , 01189 cm2s−1 de modo que o número de Reynolds baseado na micro-escala de Taylor Rλ foi aproximadamente 40. O número de Prandtl foi Pr = 0,7e o passo temporal adotado para a discretização foi ∆t= 0,0004s. 2.5 Campos Iniciais e Condições de Fronteira Os campos iniciais foram gerados segundo o método de Rogallo, garantindo isotropia e distribuição gaussiana por meio do espectro E(k)=Ak4exp(−Bk2),(29) em que A e B são determinados de modo que RE ( k ) dk = 3 / 2 (cm−1) e o máximo de E ( k ) ocorra para kp= 29/4(cm−1). Aplicaram-se condições de fronteira de corte periódico. Para qualquer variável genérica f∈ {ui, ϕ, p}, f(t, x1+m1, x2+m2, x3+m3)=f(t, x1−Sm3∆t, x2, x3),(30) sendo mium inteiro arbitrário. Numa situação de condições de fronteira periódicas com corte, a variável f no nó 1é transportada uma distância igual a Sx1 3 ∆ t . Atendendo a que o novo valor de f no nó 1deriva do valor de f no nó Nk − 1e devido ao corte a variável f no nó Nk − 1é transportada SxNk−1 3 ∆ t , significa que o valor corrigido de f no nó 1é S ( x1 3−xNk−1 3 )∆ t . Isto corresponde à aplicação de dois shift um para anular o efeito do corte no nó Nk − 1 e outro para aplicar a condição de fronteira devido ao corte no nó 1. Para a condição 7 de fronteira no nó Nk procedeu-se do mesmo modo. Assim, as condições de fronteira aplicadas no espaço de Fourier foram, f(t, I, J, 1)=f(t, I, J, Nk −1)e−Shift, f(t, I, J, Nk) = f(t, I, J, 2)eShift ,(31) onde, Shift = 2πκS (x1 3−xNk−1 3) |{z } L ∆t L= 2π(I−1)S∆t . (32) 3 RESULTADOS E DISCUSSÃO Nas próximas subseções, analisamos a evolução temporal de três indicadores principais: energia cinética da turbulência E , variância escalar ϕ2 e fluxos escalares ( uϕ, wϕ ). Realizamos simulações em malhas ortogonais, para turbulência em decaimento, e em malhas não-ortogonais inclinadas a 15 ◦ ,30 ◦ e45 ◦ , para escoamentos com corte e gradiente vertical do escalar. Para interpretar esses resultados, utilizamos as seguintes equações de transporte: dE dt=−Suw −ϵ(33) dϕ2 dt=−2swϕ −ϵϕϕ (34) duϕ dt=−suw −Swϕ +P1ϕ−ϵ1ϕ(35) dwϕ dt=−sww +P3ϕ−ϵ3ϕ(36) Definimos ainda ϵij =2 Re ∂ui ∂xk ∂uj ∂xk , ϵiϕ =1+Pr RePr ∂ϕ ∂xk ∂ui ∂xk , ϵϕϕ =2 RePr ∂ϕ ∂xk ∂ϕ ∂xk eϵ=ϵii/2,(37) que correspondem às taxas de dissipação, respectivamente, das tensões de Reynolds, dos fluxos escalares e da variância escalar. Além disso, as correlações pressão-escalar são dadas por Piϕ =p∂ϕ ∂xi .(38) Aqui, a barra denota média de Reynolds. 3.1 Energia Cinética da Turbulência e Variância Escalar A figura 3 mostra que, em presença de gradientes médios de velocidade e escalar, tanto a energia cinética da turbulência E quanto a variância escalar ϕ2 atingem valores superiores aos observados no caso sem corte ( S = 0). Esse comportamento decorre diretamente dos termos de corte das equações (17) e (18) e (33) e (34) , que, ao dependerem dos parâmetros de corte S e de estratificação s , amplificam as flutuações de velocidade e do escalar. Do ponto de vista numérico, verifica-se que a inclinação da malha (0°, 15°, 30° e 45°) praticamente não afeta a evolução temporal dessas quantidades, evidenciando a robustez tanto do modelo matemático quanto do método de discretização adotado. 8 Figura 3: Energia cinética da turbulência (a) e variância escalar (b), perante o corte. 3.2 Campo de Velocidade no Plano de Corte A Figura 4 ilustra como o campo de velocidade é distorcido pela presença do gradiente médio de velocidade. Em instantes sucessivos no plano de corte médio xz , observa-se que as estruturas de grande escala são transportadas e alongadas na direção x pela componente média U , especialmente em cotas z mais elevadas, onde U é maior. Nas regiões de maior velocidade média, as estruturas maiores se fragmentam em filamentos de menor área superficial, refletindo a cascata de energia típica da turbulência. Quanto aos aspectos numéricos, note-se que a inclinação do domínio de integração não afeta a variação temporal dessas quantidades, o que constitui uma primeira evidência de que o modelo matemático e as técnicas numéricas adotadas são adequados para a resolução deste problema. Além disso, não foram encontradas diferenças significativas entre domínios ortogonais e não-ortogonais. 3.3 Campo Escalar no Plano de Corte A figura 4 apresenta o comportamento do campo escalar no plano de corte médio xz . Assim como no caso da velocidade, as estruturas escalares são transportadas pela componente média do fluxo, porém logo se observam alongamentos de aproximadamente 45°resultantes da combinação dos gradientes de velocidade e escalar. Nas regiões de maior velocidade média, identifica-se um acúmulo local dessas estruturas, sugerindo a formação de padrões pontuais que podem refletir tanto efeitos físicos de estratificação quanto artefatos numéricos. Ao comparar as figuras 4 e 5, fica claro que, enquanto o corte estende as estruturas de velocidade ao longo de x , o campo escalar inclina-se em 45°e concentra-se nos cantos do domínio. Para aprofundar essa observação, analisamos a evolução temporal dos fluxos escalares uϕ e wϕ (figura 6). Constatamos que uϕ cresce continuamente ao longo do tempo, indicando transporte escalar predominante na direção x , enquanto wϕ atinge rapidamente um valor assintótico, reforçando a ideia de acúmulo escalar nas regiões de maior velocidade. Quanto aos aspectos numéricos, a inclinação da malha (0°, 15°, 30°e 45°) não afeta significativamente a evolução temporal do campo escalar, nem emergem diferenças perceptíveis entre domínios ortogonais e não-ortogonais. Esse comportamento confirma a 9