Resumo

No presente trabalho é realizado um estudo sobre a aplicação do método CBS (“Characteristic Based Split”) na análise de escoamentos incompressíveis, transientes e turbulentos com o objetivo de avaliar a eficiência computacional e precisão da formulação em problemas tipicamente encontrados na Engenharia do Vento Computacional (EVC). As equações fundamentais do escoamento são aqui discretizadas a partir da aplicação do método CBS no contexto do Método dos Elementos Finitos, onde elementos triangulares e tetraédricos lineares são empregados. Uma versão totalmente explícita é obtida a partir do método de compressibilidade artificial e uma versão semi-implícita é também testada, onde as equações de momentum são resolvidas explicitamente enquanto que a equação de massa é resolvida implicitamente através de uma equação de Poisson. Para a análise de escoamentos turbulentos utiliza-se a Simulação de Grandes Escalas (“Large Eddy Simulation” – LES) com o modelo clássico de Smagorinsky para as escalas inferiores à resolução da malha. São analisados primeiramente alguns problemas clássicos da Dinâmica do Fluido Computacional (DFC) para verificação e avaliação do código. Posteriormente, o esquema numérico é validado em casos típicos da EVC, como a ação do vento sobre uma seção de ponte de grande vão e escoamento de camada limite sobre um modelo topográfico de morro. Conclui-se que a ferramenta numérica proposta é capaz de obter resultados consistentes com os fenômenos envolvidos e de forma mais eficiente que outros métodos de estabilização disponíveis.

Palavras-chave: Dinâmica de Fluidos Computacional (DFC), Simulação de Grandes Escalas (LES), Método CBS, Método dos Elementos Finitos (MEF).

Abstract

In the present work, an investigation on Characteristic Based Split (CBS) method applied to the analysis of incompressible, transient and turbulent flows is proposed aiming at the evaluation of computational efficiency and accuracy of the present formulation considering problems typically found in Computational Wind Engineering (CWE). The flow fundamental equations are here discretized using the CBS method in the context of the Finite Element Method (FEM), where triangular and linear tetrahedral elements are employed. A fully explicit scheme is obtained applying the artificial compressibility method and a semi-implicit scheme is also tested, where the momentum equations are solved explicitly while the mass equation is solved implicitly through a Poisson equation. For the analysis of turbulent flows Large Eddy Simulation (LES) is adopted with the classic Smagorinsky’s model for sub-grid scales. Some classical Computational Fluid Dynamics (DFC) problems are analyzed first for code verification and evaluation. Subsequently, the numerical model is validated considering typical CWE applications, such as the wind action on a large-span bridge deck and the boundary-layer flow over a hill model. It is concluded that the proposed numerical tool is able to obtain results consistent with the phenomena involved and more efficiently than other stabilization methods available.

Keywords: Computational Fluid Dynamics (CFD), Large Eddy Simulation (LES), CBS Method, Finite Element Method (FEM)

1. Introdução

O estudo de problemas da Engenharia do Vento tem grande importância em várias áreas da Engenharia Civil, possibilitando a determinação das condições ambientais ao redor e no interior das edificações, além da avaliação de cargas aerodinâmicas e fenômenos de instabilidade induzidos pela ação do vento em estruturas [1]. Embora o emprego de técnicas experimentais em túneis de vento tem sido a prática usual neste campo de pesquisa, cresce cada vez mais o uso de simulações numéricas em razão do avanço da tecnologia dos computadores e dos métodos numéricos, levando ao desenvolvimento da chamada Engenharia do Vento Computacional (EVC).

As diferenças encontradas entre as aplicações usuais da Dinâmica de Fluidos Computacional (DFC) e as aplicações da EVC são significativas. Löhner et al. destaca que as exigências computacionais para os problemas de Engenharia Civil são muito maiores do que, por exemplo, nas engenharias aeroespacial e naval, principalmente pela física mais complexa do escoamento e das estruturas envolvidas, levando a níveis de refinamento superiores tanto em relação a malha como em relação a passos de tempo usados [2]. Dificuldades encontradas em aplicações da EVC em comparação com outras áreas da DFC já haviam sido apontadas por Murakami [3,4] e Stathopoulos [5], onde os motivos principais indicados para isso seriam a existência de obstáculos no campo de escoamento, produzindo separação da camada limite e recirculação, a presença de camada limite junto à superfície dos terrenos, a presença de obstáculos com arestas vivas, a turbulência, a incompressibilidade e a tridimensionalidade do escoamento, número de Reynolds elevados, além das fortes interações fluido-estrutura.

As simulações da DFC são extremamente sensíveis aos parâmetros computacionais definidos pelo próprio usuário. Para uma simulação típica o usuário deve selecionar as variáveis adequadas, as equações que caracterizam o escoamento, o esquema de discretização, o domínio computacional, o modelo de turbulência, as condições de contorno, critérios de convergência, além de diversas outras definições. Desta forma, as análises experimentais continuam indispensáveis e a EVC, portanto, deve ser pensada inicialmente como uma ferramenta complementar [6], sempre necessitando que os resultados numéricos obtidos das simulações sejam verificados e validados utilizando dados experimentais. Algumas instruções para melhor aplicação e validação de códigos da DFC podem ser encontradas em [7,8,9,10].

O emprego do método de Galerkin clássico nas equações fundamentais do escoamento produz soluções com oscilações espúrias sobre o espaço de análise devido à presença do termo convectivo, sobretudo quando o número de Peclet crítico é excedido. No contexto do MEF, esquemas de estabilização foram desenvolvidos a partir da modificação das funções de base utilizadas no método de Galerkin, obtendo-se o chamado método de Petrov-Galerki [11] e, posteriormente, o método SUPG – “Streamline Upwind/Petrov-Galerkin” [12]. Em problemas transientes, o esquema de Lax-Wendroff, empregado em diferenças finitas, tem seu equivalente no contexto do MEF através do método de Taylor-Galerkin [13].

Por outro lado, o método das características pode também ser utilizado para resolução dos problemas das oscilações espúrias em escoamentos, pois ele elimina o termo convectivo das equações do tipo advecção-difusão, deixando-a assim, auto-adjunta. Porém, a solução dessa nova equação apresenta algumas dificuldades, as quais podem ser solucionadas através da aplicação de séries de Taylor, possibilitando, então, a aplicação do método de Galerkin em sua forma clássica [14]. Para que seja realizada a correta compatibilização entre os campos de pressão e velocidade, Brezzi comprovou a necessidade de imposição de determinadas restrições às funções de interpolação utilizadas em elementos finitos [15]. No caso do acoplamento pressão-velocidade em escoamentos incompressíveis, é adotado geralmente um processo de decomposição/projeção das equações de Navier-Stokes, proposta originalmente por Chorin [16] (ver também [17,18]). A aplicação do procedimento de decomposição no contexto do método das características foi realizada inicialmente por Zienkiewicz e Codina [19], sendo chamado de método CBS (“Characteristic Based Split”).

Chorin propôs ainda uma alternativa para o tratamento numérico da equação da continuidade através da hipótese de pseudo-compressibilidade, onde se considera a presença de uma leve compressibilidade no escoamento a fim de obter o campo de pressão explicitamente [16]. Entretanto, muitas vezes o parâmetro de pseudo-compressibilidade, a velocidade do som, apresenta valores muito altos para escoamentos incompressíveis, levando a grandes restrições no passo de tempo em esquemas explícitos. Para contornar este inconveniente, Nithiarasu propõe a utilização de um parâmetro de compressibilidade artificial calculado localmente, baseado nas velocidades convectiva e difusiva do escoamento, deixando assim o esquema mais eficiente [20]. Para resolver problemas transientes, é necessário ainda a utilização de um esquema de passo de tempo dual. A utilização de um parâmetro de compressibilidade artificial em conjunto com o método CBS leva ao método CBS-AC (“Characteristic Based Split-Artificial Compressibility”).

Assim, problemas envolvendo escoamentos incompressíveis podem ser resolvidos usando-se uma formulação mista (variáveis de pressão e velocidade) através do método da projeção [17], o qual, sendo aplicado de forma consistente, permitirá livremente a escolha das funções de interpolação para pressão (Np) e velocidade (Nv) [14], sendo, neste caso, adotada geralmente a opção Np = Nv.

As variantes do método CBS incluem as formas explícita, semi-implícita e quase-implícita. Nos três casos, a pressão é removida da equação de balanço de momentum, estabelecendo um campo intermediário de velocidade em um primeiro passo. Após a obtenção do campo de pressão, a partir da equação de balanço de massa no segundo passo, o campo de velocidades inicial é corrigido em um terceiro passo. Esta decomposição acaba produzindo um efeito de aumento na eficiência computacional em relação a abordagens monolíticas. No entanto, sabe-se que métodos de projeção clássicos introduzem um erro temporal de primeira ordem na pressão quando ela é completamente removida da equação de momentum no processo de decomposição [21]. Caso a pressão seja parcialmente removida, instabilidades no campo de pressão são introduzidas, as quais podem ser corrigidas usando uma estabilização extra.

Para os esquemas semi-implícitos, as equações de momentum são resolvidas de uma maneira parcial ou totalmente explícita e o campo de pressão é obtido a partir da solução de uma equação de Poisson. No caso de um esquema totalmente explícito, adota-se um parâmetro de compressibilidade artificial [20] e um procedimento de integração temporal com passo de tempo dual para obter a solução transiente, onde os passos de tempo são obtidos localmente e um processo iterativo é considerado para cada passo de tempo global. Em princípio, o uso de uma abordagem de passo de tempo dual possibilitaria uma solução totalmente acoplada entre todos os passos de um esquema fracionado, aproximando-se de um esquema monolítico [22]. Embora o termo de estabilização possa ser essencial para métodos totalmente explícitos a fim de evitar oscilações nas respostas para números de Reynolds elevados, sua importância ainda não é clara para os demais esquemas [23]. Por outro lado, em um trabalho recente, Nithiarasu et al. apontam que uma estabilização adicional para a pressão não é necessária para a obtenção de uma precisão temporal de alta ordem quando emprega-se o conceito de compressibilidade artificial em esquemas explícitos [24].

O tratamento numérico de escoamentos turbulentos na EVC pressupõe o uso de modelos de turbulência. Os modelos mais utilizados na área são os modelos RANS (“Reynolds-Averaged Navier-Stokes”) e os modelos de simulação de grandes escalas (“Large-Eddy Simulation – LES). Como o uso dos modelos RANS apresenta comprovadamente uma série de limitações em escoamentos com separação [3], embora sendo mais eficientes computacionalmente, tem-se utilizado cada vez mais os modelos LES na EVC [25] em razão de sua precisão e de uma maior disponibilidade de recursos computacionais, já que esta metodologia requer malhas muito mais refinadas que simulações com modelos RANS. Os modelos LES tiveram início no trabalho de Smagorinsky [26] e consistem na separação das escalas do escoamento em função da resolução da malha (“Subgrid-Scales” – SGS). Nesses modelos, as grandes escalas do escoamento turbulento são simuladas diretamente, enquanto as pequenas escalas, associadas a escalas inferiores à resolução da malha, aparecem na formulação em termos de componentes de tensão sub-malha, as quais devem ser representadas através de modelos de turbulência [27]. Para a simulação das pequenas escalas pode-se empregar o modelo clássico de Smagorinsky [26]. Revisões sobre a aplicação de modelos LES na EVC podem ser encontradas, por exemplo, em [3,4,6].

Neste trabalho, apresenta-se uma ferramenta numérica baseada no método CBS (“Characteristic Based Split”) e no conceito de compressibilidade artificial a fim de obter-se um esquema numérico mais eficiente para futuras aplicações na EVC, onde características de escoamento tais como incompressibilidade e turbulência tem um papel preponderante. O Método dos Elementos Finitos é empregado na discretização espacial das equações fundamentais do escoamento, onde triângulos e tetraedros lineares são utilizados. Neste contexto, propõe-se o uso de um esquema explícito de integração no tempo, permitindo captar um espectro mais amplo de frequências do escoamento e facilitando a programação do código numérico para processamento em paralelo. Além disso, o emprego de elementos finitos triangulares e tetraédricos favorece uma futura implementação de esquemas de adaptação de malha a fim de se obter uma distribuição mais racional e econômica de elementos sobre o domínio computacional. Uma versão semi-implícita do método CBS é também testada, onde as equações de momentum são resolvidas explicitamente enquanto a equação de massa é resolvida implicitamente através de uma equação de Poisson. Exemplos clássicos da DFC são analisados para verificação do código. Posteriormente, alguns exemplos da EVC são simulados para validação da formulação proposta.

2. Equações fundamentais e métodos numéricos

As equações de conservação para um escoamento isotérmico e incompressível em uma descrição Euleriana são: a equação de conservação de momentum e a equação de conservação de massa [28]:

246px
(1)
(2)


onde é massa específica, é a componente da velocidade no fluido na direção , é o delta de Kronecker, são as componentes do vetor de forças de corpo na direção do eixo , sendo as direções dos eixos coordenados em um espaço cartesiano ortogonal e a variável tempo. As componentes do tensor de tensões totais são representadas por:

(3)


sendo a pressão termodinâmica e as componentes do tensor de tensões viscosas, descritas de acordo com o modelo de fluido newtoniano por:

(4)


sendo a viscosidade dinâmica e a viscosidade volumétrica.

Para escoamentos turbulentos, emprega-se a metodologia LES, onde é adotado um procedimento de decomposição entre variáveis de grandes e pequenas escalas do escoamento, além da filtragem espacial sobre o sistema de equações fundamentais, obtendo-se:

(5)


onde a barra superior indica uma variável de grandes escalas, é o termo correspondente às componentes do tensor de tensões sub-malha, o qual consiste nos efeitos das escalas inferiores à resolução da malha, requerendo um modelo de fechamento para sua solução. No contexto da abordagem LES, esse termo pode ser calculado através de:

(6)


sendo as componentes do tensor de taxa de deformação de grandes escalas e a viscosidade turbulenta, que deve ser determinada através de um modelo de turbulência para as escalas inferiores à malha. Empregando-se o modelo de viscosidade turbulenta proposto por Smagorinsky [26], é obtido através da seguinte expressão:

(7)


onde é a constante de Smagorinsky, que varia entre 0,1 e 0,25 [3], dependendo do tipo de escoamento, é a magnitude do tensor de tensões de deformação de grandes escalas e representa a largura do filtro em nível da malha. No caso tridimensional, emprega-se usualmente e no caso bidimensional, , onde e representam, respectivamente, o volume e a área do elemento finito .

2.1 Método CBS ("Chacteristic – Based Split")

O método das características consiste na aplicação de uma mudança de coordenadas ao longo das linhas características, a fim de remover os termos advectivos das equações fundamentais do escoamento. Desta forma, ao aplicar o método CBS sobre a forma não conservativa das equações de Navier-Stokes e adotando uma decomposição , chega-se a:

(8)


Obtendo em x os termos acima definidos em x - Δx a partir de aproximações em séries de Taylor no espaço, considerando-se que e desprezando os termos de ordem 3 ou superiores, obtém-se:

(9)


Desta forma, o processo de decomposição pode ser definido pelas seguintes expressões:

(10)


(11)


No presente trabalho são abordadas duas soluções para equação de conservação de massa, a primeira consiste em um esquema explícito, com o uso de um parâmetro e compressibilidade artificial e processo de avanço dual no tempo, quando analisados escoamentos transientes [20]. A outra adota uma solução para a equação de conservação de massa através de um esquema implícito [23]. Nas expressões acima, considera-se que as equações são escritas em termos das variáveis de grandes escalas e quando o escoamento é turbulento.

2.1.1 Método explícito

O método explícito parte da equação de balanço de massa escrita de acordo com o método da pseudo-compressibilidade de Chorin [16], que considera a presença de uma compressibilidade artificial, através da substituição do termo transiente da densidade por um termo equivalente de pressão relacionado a densidade, conforme apresentado a seguir:

(12)


Na presente formulação, a velocidade do som c constante na equação de conservação de massa é substituída pelo parâmetro de compressibilidade , calculado localmente em função das condições locais e atuais do escoamento. Consequentemente, será necessária a utilização de passos de tempo locais no procedimento de integração no tempo [20].

Assim, a equação de balanço de massa, assumindo e desprezando os termos 3ª ordem ou superiores, pode ser definida de forma explícita conforme apresentada abaixo, onde varia de 0,5 a 1,0:

(13)


No presente esquema, é calculado levando-se em conta a condição apresentada abaixo:

(14)


sendo uma constante que deve garantir que seu valor não se aproxime de zero em nenhum caso, a velocidade de convecção e a velocidade de difusão, definidas como:

(15)


(16)


onde é a dimensão característica do elemento finito, sendo calculada para o caso de elementos triangulares e tetraédricos da seguinte forma:

(17)
(18)


onde é o número de elementos conectados ao nó , é o volume do elemento tetraédrico, é a área da face do tetraedro, ou, no caso do triangulo, a área do elemento, e é o comprimento da aresta do triangulo oposta ao nó .

O passo de tempo local é calculado conforme a condição apresentada abaixo:

(19)


Os valores de passo de tempo de convecção e passo de tempo de difusão são definidos por:

(20)


(21)


sendo Re o número de Reynolds do escoamento. O valor de a ser utilizado é calculado multiplicando-se um fator de segurança sobre a Equação (19), o qual varia entre 0,5 e 2,0 dependendo do problema e da malha empregada.

A aplicação do método do CBS com compressibilidade artificial em problemas envolvendo escoamentos estacionários pressupõe a solução dos 3 passos indicados pelo esquema de decomposição de forma sequencial, até que o estado estacionário seja atingido. No entanto, para a aplicação do método em escoamentos transientes, um processo de avanço dual no tempo deve ser empregado, onde o tempo físico é divido em um determinado número de passos de tempo, denominados passos de tempo real, os quais são subdivididos em passos de tempo fictícios necessários para a convergência das variáveis do escoamento. Além disso, um termo adicional associado ao passo de tempo real deve ser considerado na equação de momentum a fim de estabelecer a solução transiente verdadeira [29,30].

Aplicando o método de Galerkin convencional (Bubnov-Galerkin) sobre as equações discretizadas pelo método CBS-AC no contexto do MEF, obtêm-se os seguintes passos de solução para o escoamento para cada passo de tempo:

Passo 1:

(22)


Passo 2:

(23)


Passo 3:

(24)


onde e são matrizes de massa discreta para as equações de momentum e massa, respectivamente, A é a matriz de convecção, D ij são as matrizes de difusão, G i são matriz gradientes, f é o vetor que inclui forças de corpo e termos de contorno, S ν e S p são matrizes de estabilização para as equações de momentum, H é a matriz laplaciana integrada por partes, os vetores V i e ΔV i contêm os valores nodais de velocidade e incremento de velocidade, p e Δp são vetores com valores nodais de pressão, é o passo de tempo fictício e é o passo de tempo. O incremento de velocidade referente ao processo transiente real , é aproximado neste trabalho por:

(25)


onde é o vetor velocidade tomado no passo fictício atual e e referem-se a vetores velocidade tomados no passo real atual e anterior, respectivamente.

2.1.2 Método semi-implícito

O esquema semi-implícito consiste na solução das equações de conservação de momentum de forma explícita, sendo somente o termo de pressão solucionado de forma implícita. Logo, a decomposição das equações de conservação de momentum sofre alteração somente na parcela da correção da velocidade, no terceiro passo, na qual o termo de pressão é definido no tempo n+ e não há necessidade da adição do termo transiente real, já que nesse método é utilizado um passo de tempo único para toda a malha. O processo de decomposição no caso semi-implícito pode ser definido pelas seguintes expressões:

(26)
(27)


As diferenças na formulação deste esquema em relação ao esquema explícito apresentado anteriormente estão somente no segundo e terceiro passos do método de projeção. Assim, ao substituir por na equação de balanço de massa, além de desprezar os termos de 3ª ordem ou superiores, obtém-se:

(28)


onde o valor adotado para varia de 0,5 a 1,0 e o valor de varia entre 0 e 1.

Aplicando o método convencional de Galerkin no contexto do MEF sobre as equações acima e integrando por partes os termos de ordem 2, obtêm-se os 3 passos do método de projeção para o esquema semi-implícito:

Passo 1:

(29)


Passo 2:

(30)


Passo 3:

(31)


O passo de tempo é calculado conforme a condição apresentada abaixo:

(32)


onde os valores de passo de tempo de convecção e passo de tempo de difusão são definidos por:

(33)
(34)


O valor de a ser utilizado é calculado multiplicando-se um fator de segurança sobre a Equação (32), o qual varia entre 0,5 e 2,0 dependendo do problema e da malha empregada.

3. Exemplos clássicos

Foram analisados aqui dois exemplos clássicos da Dinâmica de Fluidos Computacional, a fim de verificar a funcionalidade do algoritmo elaborado, baseado no método CBS. Foram utilizadas as duas formulações apresentadas: o esquema explícito, para o qual os exemplos referem-se a escoamentos estacionários e transientes, a fim de verificar a efetividade do esquema de decomposição do método de projeção em escoamentos estacionários e o processo de avanço dual no tempo para escoamentos transientes, onde se considera um termo adicional associado ao passo de tempo real; e o método semi-implícito, para o qual os exemplos envolvem a aplicação do esquema de decomposição do método de projeção em escoamentos transientes, sem a adoção de um parâmetro de compressibilidade artificial. Exemplos envolvendo escoamentos turbulentos também são analisados, verificando a estabilidade e precisão do presente modelo quando se utiliza a abordagem LES juntamente com o modelo sub-malha clássico de Smagorinsky.

3.1 Escoamento bidimensional em uma cavidade

Foram realizadas três análises com diferentes números de Reynolds (RE100, RE1000, RE10000), empregando-se o método CBS explicito com duas malhas, uma composta por elementos triangulares e a outra com tetraedros. A primeira malha é formada por 8.836 nós e 17.298 elementos triangulares, com menor dimensão do elemento próxima a 5x10-3 m,a segunda malha é composta por 50.000 elementos tetraédricos e 20.402 nós com menor dimensão de elemento próxima a 5x10-3 m.As características geométricas e as condições de contorno empregadas nestes exemplos estão apresentadas na Figura 1, onde as coordenadas são representadas em metros. A espessura da cavidade com a malha tridimensional corresponde a um elemento e a terceira componente da velocidade é nula em todo domínio; portanto os resultados obtidos neste caso devem ser idênticos aos resultados obtidos nos casos com malha bidimensional.

Review 391264738980-image85.png Review 391264738980-image86.png
(a) (b)
Figura 1. Condições geométricas e condições de contorno para escoamento em uma cavidade: (a) Malha composta por elementos triangulares. (b) Malha composta por elementos tetraédricos com um elemento de espessura


As constantes físicas e geométricas utilizadas nestas análises encontram-se relacionadas na Tabela 1, sendo que a constante de Smagorinsky é utilizada apenas para Reynolds igual a 10.000. Além disso, o número de Reynolds desejado foi obtido variando o valor de viscosidade dinâmica, mantendo constantes a massa específica, a velocidade característica de referência e a dimensão característica de referência.

Tabela 1. Constantes físicas e geométricas utilizadas nas análises do escoamento bidimensional em uma cavidade.
Massa específica (ρ) 1,0 kg/m³
Viscosidade dinâmica (µ) 10,0/Re kg/m.s
Viscosidade volumétrica (λ) 0,0 kg/m.s
Velocidade da placa superior (Vo) 10,0 m/s
Dimensão característica (L) 1,0 m
Constante de Smagorinsky (Cs) 0,10 ou 0,20


O primeiro caso analisado foi a cavidade com RE100, para o qual a viscosidade dinâmica vale 0,10 kg/m.s e o escoamento é estacionário. Na Figura 2 estão apresentados os resultados de linhas de corrente e campo de pressão obtidos para este caso. Verifica-se que as linhas de corrente do escoamento e as linhas isobáricas apresentadas estão de acordo com as referências [31,32]. Quanto às linhas de corrente, observa-se a formação de um grande vórtice principal, com centro bem definido, e outros dois vórtices, conhecidos como secundários, nos cantos inferiores da cavidade.

Review 391264738980-image87.png Review 391264738980-image88-c.png
(a) (b)
Figura 2. Resultados no interior da cavidade com RE100: (a) Linhas de corrente. (b) Linhas isobáricas


Nas linhas médias vertical e horizontal da cavidade foram tomados os perfis de velocidade para RE100 e plotados em conjunto com os valores obtidos por Ghia et al. [31], onde pode-se observar uma ótima concordância entre os valores de referência e as predições aqui apresentadas com elementos triangulares e tetraédricos, conforme apresentado na Figura 3.

Review 391264738980-image89-c.png Review 391264738980-image90-c.png
(a) (b)
Figura 3. Perfis de velocidade nas linhas médias no interior da cavidade RE100: (a) Perfil de velocidade em X = 0,5 m. (b) Perfil de velocidade em Y = 0,5 m


A segunda análise consiste no escoamento transiente na cavidade com RE3200 e viscosidade dinâmica igual a 3,125.10-3 kg/m.s, sem a utilização de um modelo de turbulência. Na Figura 4 podem ser observadas as linhas de corrente no interior da cavidade com a formação mais centralizada do vórtice primário, assim como o aumento dos vórtices secundários nos cantos inferiores e superior. Os perfis de velocidade nas linhas médias aqui obidos para malha de triângulos e tetraedros são plotados em conjunto com os resultados obtidos por Ghia et al. [31], onde verifica-se uma boa concordância.

Review 391264738980-image91.png Review 391264738980-image92-c.png
(a) (b)
Figura 4. Resultados no interior da cavidade com RE3200: (a) Linhas de corrente. (b) Perfis de velocidade nas linhas médias no interior da cavidade em X = 0,5 m e em Y = 0,5 m


Por fim, a última análise consiste em um escoamento turbulento com Reynolds igual a 10.000 e viscosidade dinâmica igual a 10-3 kg/m.s. Para a simulação deste escoamento foi necessária a utilização de um modelo de turbulência, no qual foram testados dois valores diferentes de constante de Smagorinsky (CS), para cada uma das malhas utilizadas. Os perfis de velocidade média obtidos a partir das linhas médias vertical e horizontal da cavidade estão apresentados na Figura 5, em conjunto com os valores de referência obtidos por Ghia et al. [31]. Observa-se que os valores apresentam uma certa diferença nas regiões próximo às paredes laterais da cavidade em função dos diferentes valores adotados para a constante de Smagorinsky (CS). Para a malha de triângulos, um melhor ajuste foi obtida usando-se CS = 0,2, enquanto que para a malha de tetraedros o melhor resultado foi obtido com CS = 0,1.

Review 391264738980-image93-c.png Review 391264738980-image94-c.png
(a) (b)
Figura 5. Perfis de velocidade nas linhas médias no interior da cavidade RE10000: (a) Perfil de velocidade em X = 0,5 m. (b) Perfil de velocidade em Y = 0,5 m


Review 391264738980-image95.png Review 391264738980-image96.png
(a) (b)
Figura 6. Resultados no interior da cavidade com RE10000: (a) Linhas de corrente. (b) Linhas isobáricas


Na Figura 6 (a) são apresentados o campo de pressão e as linhas de corrente referentes ao campo médio do escoamento. Pode-se verificar que as linhas de corrente obtidas no presente trabalho estão de acordo com os resultados apresentados pela referência [31]. Além disso, percebe-se que os vórtices secundários nos cantos inferiores e superior da cavidade tornam-se mais evidentes à medida que o número de Reynolds aumenta. Na Figura 6 (b) são apresentadas as linhas isobáricas no interior da cavidade, onde percebe-se uma pequena instabilidade no campo de pressões junto ao canto superior esquerdo da cavidade.

Portanto, para as análises realizadas envolvendo escoamentos bidimensionais, verifica-se que as formulações com elementos triangulares e tetraédricos conseguiram apresentar resultados similares aos obtidos pela referência. A próxima etapa da verificação é a análise de um escoamento no interior de uma cavidade tridimensional.

3.2 Escoamento tridimensional em uma cavidade

A simulação deste escoamento tridimensional considera um número de Reynolds igual a 3.200, com modelo de turbulência clássico de Smagorinsky. Resultados numéricos para cavidade tridimensional foram obtidos anteriormente por Zang et al. [33]. As características geométricas e as condições de contorno adotadas estão apresentadas na Figura 7, onde as coordenadas são dadas em metros. As dimensões da cavidade foram iguais a 1 metro na direção X, 1 metro na direção Z e 1 metro na direção Y. Porém, para redução do número de elementos necessários para representação do domínio, foi considerado um plano de simetria no escoamento em Y = 0,5 m, simulando apenas metade do domínio. Desta forma, a malha é composta por 312.500 elementos tetraédricos e 67.676 nós, com menor dimensão de elemento igual a 3x10-2 m.

Review 391264738980-image97.png
Figura 7. Condições geométricas e condições de contorno para escoamento em uma cavidade tridimensional


Os resultados aqui obtidos sem a utilização de um modelo de turbulência não foram satisfatórios para o caso tridimensional e, desta forma, foram utilizados dois valores de constante de Smagorinsky, correspondente a 0,1 e 0,2 para analisar sua influência sobre os resultados. As constantes físicas e geométricas utilizadas nesta análise encontram-se relacionadas na Tabela 2.

Tabela 2. Constantes físicas e geométricas utilizadas nas análises do escoamento em uma cavidade tridimensional.
Massa específica (ρ) 1,0 kg/m³
Viscosidade dinâmica (µ) 0,03125 kg/m.s
Viscosidade volumétrica (λ) 0,0 kg/m.s
Velocidade da placa superior (Vo) 100,0 m/s
Constante de Smagorinsky (Cs) 0,10 e 0,20
Dimensão característica (L) 1,0 m


A seguir, na Figura 8 estão apresentados os resultados obtidos para os perfis de velocidade média nas linhas do centro geométrico vertical e horizontal no plano de simetria da cavidade, com os diferentes valores da constante de Smagorinsky, os quais se referem ao campo médio do escoamento. Em conjunto com os resultados obtidos no presente trabalho, estão apresentados os resultados numéricos obtidos por Zang et al. [33]. Nota-se que os resultados referentes a CS = 0,1 apresentam uma muito boa concordância com os resultados de referência. A variabilidade da solução em função da constante de Smagorinsky evidencia a necessidade de um estudo prévio sobre as condições de escoamento a fim de obter o valor da constante que melhor reproduz as suas características. Para evitar isso, pode-se utilizar um modelo dinâmico para avaliação automática da constante de Smagorinsky em função das condições instantâneas do escoamento (ver, por exemplo, Braun [25]).

Review 391264738980-image98-c.png Review 391264738980-image99-c.png
(a) (b)
Figura 8. Perfis de velocidade nas linhas médias do plano de simetria (Y = 0,5 m) em uma cavidade tridimensional – RE3200: (a) Perfil de velocidade em X = 0,5 m. (b) Perfil de velocidade em Z = 0,5 m


Foram obtidas também as linhas de corrente e linhas isobáricas no plano de simetria da cavidade, conforme pode-se observar na Figura 9. Neste plano é possível verificar a existência de um vórtice principal na cavidade, localizado mais próximo à parede lateral direita, e de dois vórtices secundário nos cantos inferiores.

Review 391264738980-image100-c.png Review 391264738980-image101-c.png
(a) (b)
Figura 9. Resultados obtidos no presente trabalho no plano de simetria (Y = 0,5 m) de uma cavidade tridimensional – RE3200: (a) Linhas de corrente. (b) Linhas isobáricas

3.3 Escoamento bidimensional sobre um cilindro fixo

Neste item aborda-se o problema do escoamento ao redor de um cilindro fixo. Ao analisar esse escoamento, os resultados de maior interesse são os coeficientes de sustentação (CD) e de arrasto (CL), que são, respectivamente, as resultantes das forças paralelas e perpendiculares sobre o corpo em relação à direção do escoamento, além dos coeficientes de pressão (CP) e o número de Strouhal (St), que avalia o fenômeno de formação e desprendimento de vórtices.

O escoamento ao redor de um cilindro apresenta uma dependência entre o número de Reynolds e a configuração do escoamento a jusante do corpo [34]. Em casos onde o número de Reynolds é menor que 5, observa-se que não há separação do escoamento sobre a superfície do cilindro. Em casos onde o número de Reynolds é inferior a 40, aproximadamente, espera-se a formação de dois vórtices simétricos em uma região de circulação atrás do cilindro. Para Reynolds superiores a 40, aproximadamente, há o fenômeno de desprendimento regular e alternado de vórtices, formando a esteira de vórtices de Von Kármán.

No presente trabalho foram analisados o problema de escoamento sobre um cilindro fixo com números de Reynolds 40, 100 e 1000. É importante lembrar que escoamentos sobre cilindros passam a apresentar um comportamento turbulento e tridimensional mais evidente a partir de Re = 300, aproximadamente. Assim, a presente simulação para Re = 1000 tem um caráter puramente acadêmico, procurando reproduzir resultados numéricos obtidos por outros autores nestas mesmas condições. As simulações foram realizadas com duas malhas diferentes, uma composta por 41.400 elementos triangulares e 20.940 nós, e a outra composta por 103.500 elementos tetraédricos e 41.880 nós. Ambas apresentam a mesma distribuição de elementos ao redor do corpo, com menor espaçamento entro os nós próximo de 10-2 m. O escoamento é considerado bidimensional em todos os casos, sendo que no caso da malha com tetraedro, o domínio é composto por apenas um elemento de espessura. Além das diferentes malhas, foram utilizadas duas formulações diferentes (explícita e semi-implícita) para o método CBS para o caso com elementos triangulares. Na Figura 10, estão apresentadas as características geométricas e condições de contorno do problema. O centro do cilindro, com diâmetro igual a 1 m, está localizado nas coordenadas X = Y = 8 m. As coordenadas estão em metros e os valores de velocidade em metros por segundo.

Review 391264738980-image102-c.png
(a)
Review 391264738980-image103-c.png
(b)
Figura 10. Condições geométricas e condições de contorno em um escoamento ao redor de um cilindro: (a) malha composta por elementos triangulares. (b) Malha composta por elementos tetraédricos com um elemento de espessura


As constantes físicas e geométricas são consideradas idênticas nas duas configurações da malha, sendo apresentadas na Tabela 3. Cabe salientar que para simulação de diferentes números de Reynolds os parâmetros massa específica, velocidade característica de referência e dimensão característica de referência são mantidos constantes, alterando-se apenas o valor da viscosidade dinâmica. A única diferença adotada nas constantes físicas das simulações entre as formulações explicita e semi-implícita é que, no caso semi-implícito, tem-se uma definição prévia de um passo de tempo único para toda malha, igual a 2,3x10-3 s.

Tabela 3. Constantes físicas e geométricas adotadas no escoamento ao redor de um cilindro.
Massa específica (ρ) 1,0 kg/m³
Viscosidade dinâmica (µ) 1,0/Re kg/m.s
Viscosidade volumétrica (λ) 0,0 kg/m.s
Dimensão característica do cilindro (D) 1,0 m
Velocidade de entrada (Vo) 1,0 m/s


Para simulação do problema com Reynolds 40 foi utilizada uma viscosidade dinâmica igual a 0,025 kg/m.s. Verifica-se que os resultados de campo de pressão e linhas de correntes obtidos após a convergência ao estado estacionário foram idênticos com as duas configurações de malha aqui empregadas. Na Figura 11 estão apresentadas as linhas de correntes ao redor do cilindro, apresentadas de forma mais detalhada, em conjunto com as dimensões características dos vórtices formados à jusante do cilindro. Na Tabela 4 estão apresentadas as dimensões características aqui obtidas para os vórtices fixos à jusante do cilindro, bem como os respectivos resultados obtidos pelas referências [25,35,36], onde pode-se verificar uma boa concordância.

Review 391264738980-image104.png Review 391264738980-image105.png
(a) (b)
Figura 11. Escoamento ao redor de um cilindro – RE40: (a) Detalhe das linhas de corrente ao redor do cilindro. (b) Dimensões características dos vórtices a jusante do cilindro


Tabela 4. Características dos vórtices e coeficientes de arrasto no escoamento ao redor de um cilindro – RE40.
Grandeza
Presente trabalho 1,62 2,22 0,72 0,60
Braun [25] 1,77 2,10 0,71 0,58
Wanderley et al. [35] 1,56 2,29 0,73 0,60
Tonon [36] 1,58 2,21 0,71 0,59


Na Figura 12 tem-se a distribuição do coeficiente de pressão (Cp) sobre a superfície do cilindro obtida neste trabalho, referente às duas configurações de malha e às duas formulações CBS com elementos triangulares aqui utilizados, bem como os resultados de referência [25,35]. Pode-se observar que os resultados obtidos com o modelo implementado estão de acordo com as referências e, além disso, não foram observadas grandes diferenças entre a malha bidimensional e a malha tridimensional com tetraedros e um elemento de espessura.


Review 391264738980-image110-c.png Review 391264738980-image111-c.png
(a) (b)
Figura 12. Distribuição do coeficiente de pressão (Cp) sobre a superfície do cilindro – RE40: (a) Elementos triangulares. (b) Elementos tetraédricos


A outra análise aqui realizada foi com número de Reynolds igual a 100. Para essa simulação, a viscosidade dinâmica adotada foi igual a 0,01 kg/m.s. Os valores de coeficiente de arrasto e sustentação, obtidos com as diferentes malhas foram idênticos. Na Figura 13 os coeficientes de arrasto e sustentação com as diferentes formulações utilizadas com o elemento triangular são apresentados, onde é possível perceber uma pequena diferença nos valores de arrasto obtidos e uma pequena defasagem no histórico de sustentação, embora o comportamento geral seja bastante próximo. Na Figura 14 é apresentado o campo de pressão já desenvolvido em t = 200 s, onde a formação e desprendimento de vórtices regulares e alternados é observada, além da esteira de von Karman na região à jusante do cilindro.

Review 391264738980-image112-c.png Review 391264738980-image113-c.png
(a) (b)
Figura 13. Resultados dos coeficientes aerodinâmicos sobre um cilindro – RE100: (a) Histórico de coeficiente de arrasto (CD). (b) Histórico de coeficiente de sustentação (CL)


Review 391264738980-image114.png
Figura 14. Linhas de corrente e campos de pressão para escoamento ao redor de um cilindro fixo, t=200s – RE100


Na Figura 15 tem-se a distribuição do coeficiente de pressão (Cp) médio sobre a superfície do cilindro obtida com as duas malhas e formulações CBS aqui empregadas para o presente escoamento, além dos resultados obtidos por Tonon [36] usando Análise Isogeométrica. Pode-se observar que os resultados obtidos com o modelo implementado para ambos os elementos (triângulos e tetraedros) estão de acordo, assim como com as diferentes formulações utilizadas para malha com elementos triangulares.

Review 391264738980-image115-c.png Review 391264738980-image116-c.png
(a) (b)
Figura 15. Distribuição do coeficiente de pressão (Cp) sobre a superfície do cilindro – RE100: (a) Elementos triangulares. (b) Elementos tetraédricos


Por fim, foi analisado o escoamento com número de Reynolds igual a 1.000. Na Figura 16 são apresentados os campos de pressão e linhas de corrente após o início do fenômeno de desprendimento de vórtices, em t = 60 s, onde pode-se observar a formação de uma esteira bem definida e o fenômeno de desprendimento alternado de vórtices na superfície do cilindro.

Review 391264738980-image117.png
Figura 16. Linhas de corrente e campos de pressão em um escoamento ao redor de um cilindro fixo em t = 60 s – RE1000


Outros resultados obtidos através desta análise foram os coeficientes de arrasto e de sustentação, conforme apresentado na Figura 17. O valor médio do coeficiente de arrasto obtido no presente trabalho foi de 1,50. A partir do histórico dos coeficientes de sustentação, pode-se obter a frequência de desprendimento de vórtices, igual a 0,236 Hz. Logo, o número de Strouhal (St = fD/V) obtido é igual a 0,236, considerando-se diâmetro e velocidade de referência com valores unitários. Verifica-se que o comportamento geral dos históricos obtidos com as diferentes formulações CBS aqui usadas é muito semelhante.

Review 391264738980-image118-c.png Review 391264738980-image119-c.png
(a) (b)
Figura 17. Resultados dos coeficientes aerodinâmicos sobre um cilindro – RE1000: (a) Histórico de coeficiente de arrasto (CD). (b) Histórico de coeficiente de sustentação (CL)


Na Figura 18 são apresentados os resultados de coeficiente de pressão (Cp) médio sobre a superfície do cilindro, onde pode-se observar a coerência entre os resultados obtidos no presente trabalho com as diferentes configurações. Os resultados são comparados com predições obtidas por Braun [25], que utilizou uma formulação em elementos finitos com hexaedros lineares e integração reduzida.

Review 391264738980-image120-c.png Review 391264738980-image121-c.png
(a) (b)
Figura 18. Distribuição do coeficiente de pressão (Cp) sobre a superfície do cilindro – RE1000: (a) Elementos triangulares. (b) Elementos tetraédricos


Tabela 5. Resultados de coeficientes de arrasto e número de Strouhal no escoamento ao redor de um cilindro.
Presente trabalho Herjford [37] Braun [25] Wanderley et al. [35] Tonon [36]
40 CD 1,62 1,60 1,77 1,56 1,58
100 CD 1,36 1,36 - 1,30 1,37
St 0,167 0,168 - 0,158 0,165
1000 CD 1,50 1,47 1,42 0,96 1,50
St 0,236 0,234 0,212 0,193 0,23


Na Tabela 5 estão apresentados de forma resumida os coeficientes aerodinâmicos e o número de Strouhal para os diferentes números de Reynolds analisados no presente estudo, assim como os resultados obtidos por outros autores [25,35,36,37]. A partir dos resultados pode-se comprovar que o modelo implementado com ambas as formulações do método CBS é capaz de reproduzir valores referentes a coeficientes aerodinâmicos, além de descrever adequadamente os processos de formação e desprendimento de vórtices para escoamentos incidindo sobre um corpo imerso.

4. Comparação entre os métodos

Durante as análises realizadas para verificação do algoritmo, alguns parâmetros tiveram bastante influência no tempo de convergência dos problemas no método CBS-AC. Para uma melhor avaliação foram realizadas algumas análises adicionais tomando por base o escoamento incidindo sobre cilindro com RE1000 e variando alguns desses parâmetros.

A primeira análise foi quanto ao passo de tempo real. Foram adotados quatro passos reais diferentes, variando entre 0,05 s e 0,5 s analisando, primeiramente, o número de passos fictícios necessários para a convergência de um passo real, conforme apresentado na Figura 19 (a). Observa-se que, quanto menor o passo de tempo real, maior o número de passos fictícios necessários para convergência do passo real. Além disso é necessário um maior número de passos fictícios para a convergência dos primeiros passos reais. O que pode ser explicado pelo fato do termo transiente da equação de momentum, apresentado na Eq. (25), sempre requerer campos de velocidade em dois instantes tempo anteriores ao atual. Portanto, ao iniciar a simulação, sempre há essa dificuldade em relações às condições iniciais. Foram analisados também os coeficientes de arrasto ao longo do tempo para os diferentes passos reais (Figura 19b), onde observa-se que para um passo real mais alto, o escoamento demora mais para desenvolver e o valor do arrasto não está de acordo com as referências. Além disso, o valor obtido para o número de Strouhal torna-se impreciso, indicando que o fenômeno de desprendimento de vórtices não é bem representado. A razão para isso é a perda de precisão na discretização temporal para o escoamento analisado.

Review 391264738980-image123-c.png Review 391264738980-image124-c.png
(a) (b)
Figura 19. Análises do método CBS-AC com diferentes passos reais: (a) número de passos fictícios versus o tempo; (b) Coeficiente de arrasto (CD) versus o tempo


A partir dessas análises foi observado que o passo real de 0,1 s apresentou bons resultados quando comparado com passos maiores, não havendo necessidade de um maior gasto computacional, como no caso do passo igual a 0,05s, para obtenção de uma precisão adequada para o caso analisado. Porém, constatou-se durante as análises de verificação que, dependendo do refinamento da malha e das condições do escoamento, o passo real de 0,1s apresenta problemas de convergência, principalmente para números de Reynolds mais altos. Logo, em alguns casos, como na cavidade com elementos tetraédricos, foi adotado um passo real inferior no valor de 0,05s. Não foram observados problemas de convergência quando adotado um passo de tempo real mais baixo, no entanto o número de passos necessários para a convergência é maior, não havendo necessidade da utilização de um passo tão baixo quando bons resultados podem ser obtidos com um passo real mais alto.

A aplicação do método CBS-AC para solucionar escoamentos transientes necessita a adição de um termo transiente no terceiro passo, conforme apresentado na Eq. (24). A adição deste termo na equação torna a formulação transiente, onde diversos estados estacionários devem ser obtidos ao longo do tempo, com a solução convergindo para resíduos de pressão e velocidade inferiores a um determinado valor prescrito. Ou seja, o tempo total é dividido em diversos passos de tempo chamados reais, sobre os quais não incide qualquer limitação quanto ao tamanho, que por sua vez são divididos em passos fictícios locais para cada elemento e onde a convergência a um estado estacionário deve ser obtida. Ao propor o método, Nithiarasu [20] não impôs limitações quanto ao passo de tempo real, porém pode-se observar que, ao adotar uma malha mais refinada, o passo de tempo real deve ser menor, caso contrário surgem problemas na convergência. A natureza explícita do termo transiente acarreta em limitações no passo de tempo local, devendo ser inferior a um quarto do passo de tempo real. Desta forma, uma malha mais refinada apresenta um pequeno passo de tempo local, deixando o esquema instável quando utilizado com passos de tempo reais muito altos [38].

O método CBS-AC utiliza um passo de tempo fictício para cada elemento, no entanto para melhor caracterização do modelo foi realizado um teste com o mesmo passo de tempo para todos elementos da malha, sendo este correspondente ao menor passo de tempo calculado entre os elementos. Verifica-se que o desenvolvimento do escoamento neste caso é mais lento do que utilizando um passo de tempo local, observando-se algumas perturbações no campo de pressão, embora nenhuma perda de precisão seja observada quanto a resultados referentes a coeficientes de arrasto e sustentação.

Outra consideração a ser realizada é quanto às diferentes formulações aplicadas na resolução do campo de pressão. Nos resultados acima foram utilizadas tanto a formulação explícita (CBS-AC) com um parâmetro de compressibilidade artificial e passo de tempo dual quanto a formulação semi-implícita, onde somente a equação de pressão é resolvido de forma implícita. Observa-se que não houve diferenças significativas entre os resultados, no entanto a formulação semi-implícita apresentou um desenvolvimento do escoamento mais rápido que a explícita. Nota-se também que o método sofre com regiões de menor refinamento da malha, sendo necessário utilizar um processo de suavização sobre o campo de pressão para eliminar as oscilações existentes. De forma geral, a formulação semi-implícita apresentou um melhor desempenho computacional.

5. Aplicações da EVC

Primeiramente é estudado um escoamento de camada limite incidindo sobre um modelo topográfico de um morro, tendo como principal objetivo a obtenção das linhas de corrente e os perfis de velocidade média sobre sua superfície. Na sequência, é realizada a simulação da ação do vento sobre uma seção de ponte de grande vão, onde o principal objetivo é a obtenção dos coeficientes aerodinâmicos, campos de pressão e linhas de corrente ao redor do corpo.

5.1 Morro isolado hipotético proposto por Ferreira et al. [39]

O morro simétrico estudado tem relação entre altura (H) e base (b) igual a 2, conforme apresentado na Figura 20, sendo a altura (H) igual a 4,0 metros e a base (b) igual a 2,0 metros. O domínio computacional, juntamente com as condições de contorno, é apresentado na Figura 21. Foi elaborada uma malha não estruturada, composta por 232.135 elementos triangulares e 116.877 nós, com a menor dimensão do elemento na ordem de 0,05 m. Inicialmente, a pressão e a velocidade são assumidas nulas em todo espaço de análise e o perfil vertical de velocidades foi obtido respeitando a lei potencial proposta pela referência [40].

Review 391264738980-image125.png
Figura 20. Geometria para modelo de morro isolado


Review 391264738980-image126.png
Figura 21. Condições geométricas e condições de contorno para morro isolado com elementos triangulares


O número de Reynolds utilizado na análise é de 105 e as demais constantes físicas e geométricas estão apresentadas na Tabela 6. Em todas as análises foi adotado o modelo de turbulência clássico de Smagorinsky e o valor da constante foi variado entre 0,1 e 0,2 para obtenção de um resultado mais próximo da referência [40].

Tabela 6. Constantes físicas e geométricas adotadas no escoamento em um morro isolado.
Massa específica (ρ) 1,0 kg/m³
Viscosidade dinâmica (µ) 0,000385 kg/m.s
Viscosidade volumétrica (λ) 0,0 kg/m.s
Dimensão característica (H) 4,0 m
Velocidade de entrada (Vo) 10,0 m/s


As linhas de corrente e o campo médio de pressão normalizado obtidos em t = 30 s são apresentadas na Figura 22 juntamente com um detalhe próximo da região do morro. Conforme esperado, as linhas de corrente tornam-se próximas na região junto ao topo do morro devido à aceleração do escoamento nesta região, observando-se o ponto de separação no cume e a formação de uma zona de recirculação à sotavento com desprendimento de vórtices de grande extensão.

Review 391264738980-image127-c.png Review 391264738980-image128-c.png
(a) (b)
Figura 22. Linhas de corrente e campo de pressão normalizado na simulação de um morro isolado: (a) em t = 30 s. (b) Detalhe das linhas de corrente na simulação próximo da superfície do morro


Na Figura 23 são observados os perfis de velocidades normalizados obtidos neste trabalho com os diferentes valores para a constante de Smagorinsky, sendo comparados em conjunto com o perfil obtido experimentalmente por Matuella [40]. Os resultados obtidos através da simulação numérica apresentaram comportamento similar a referência [40], porém pode-se notar que as predições aqui obtidas apresentam uma estimativa maior que os dados experimentais para alturas inferiores a 6,20 metros. A partir dessa cota, a simulação apresentou resultados de acordo com os dados experimentais. Apesar desta diferença, nota-se que o modelo numérico consegue reproduzir a aceleração do escoamento junto à superfície do morro através do aumento de velocidade observado em 4 < Y < 4,5 m.

Review 391264738980-image129-c.png
Figura 23. Perfis de velocidade sobre o morro


Pode-se observar que a análise que mais se aproximou dos resultados apresentados na referência [40] foi aquela realizada com a constante de Smagorinsky igual a 0,1. Na Figura 24 estão apresentados os perfis de velocidade ao longo da superfície do morro, com Cs igual a 0,1. Na região de recirculação é observada a redução e inversão no sentido das velocidades ocasionadas pelo gradiente adverso de pressão após o ponto de separação do escoamento próximo ao cume do morro.

Review 391264738980-image130.png
Figura 24. Perfis de velocidade média ao longo do morro

5.2 Ponte Great Belt East

Um caso de ponte suspensa que teve seu comportamento aerodinâmico e aeroelástico bastante analisado por diversos autores [41,42,43] é a Great Belt East Bridge, situada na Dinamarca, sendo composta por dois vãos de aproximação de 535 metros cada e um vão central com 1.624 metros de extensão, com seção transversal do tipo caixão único trapezoidal com 31 metros de largura e 4,4 metros de altura, conforme apresentado na Figura 25.

Review 391264738980-image131.png
Figura 25. Seção transversal da ponte Great Belt East


Na simulação deste problema foi utilizado o método CBS com formulação semi-implícita. O escoamento foi considerado bidimensional, com uma velocidade constante na entrada. A malha estruturada é composta por 110.400 elementos triangulares e 55.780 nós, sendo que as menores dimensões de elementos estão na ordem de 5x10-2 m. As dimensões do domínio computacional e a configuração da malha, bem como as condições de contorno, são apresentadas na Figura 26. O número de Reynolds utilizado na simulação é de 3x105. Na Tabela 7 são apresentadas as demais constantes empregadas nesta análise. O incremento de tempo usado é de 3,5x10-4 s.

Review 391264738980-image132.png
Figura 26. Condições geométricas e condições de contorno da seção de Great Belt East com elementos triangulares


Tabela 7. Constantes físicas e geométricas adotadas no escoamento na Great Belt East para RE 3x105.
Massa específica (ρ) 1,0 kg/m³
Viscosidade dinâmica (µ) 0,00058667 kg/m.s
Viscosidade volumétrica (λ) 0,0 kg/m.s
Constante de Smagorinsky (Cs) 0,1
Dimensão característica (H) 4,40 m
Velocidade de referência (Vo) 40,0 m/s


Na Figura 27 são apresentadas as linhas de corrente e campo de pressão instantâneos em t = 40 s, sendo que os valores de pressão são normalizados pela pressão dinâmica na região não perturbada do escoamento. Para comparação, são apresentadas as linhas de corrente obtidas numericamente por Kuroda [42]. Pode-se observar que o campo obtido no presente trabalho apresenta uma boa semelhança com o resultado de referência conseguindo, inclusive, reproduzir vórtices de menor magnitude não captados pela referência.

Review 391264738980-image133-c.png Review 391264738980-image134-c.png
(a) (b)
Figura 27. Campo de pressão e linhas de corrente no entorno da Great Belt East em t=40s: (a) Presente trabalho. (b) Kuroda [42]


A partir do histórico de sustentação foi obtida a frequência de desprendimento de vórtices e com o valor da velocidade na entrada igual a 40 m/s e a dimensão característica igual a 31 metros, o valor do número de Strouhal (St = fD/V) fica em 0,13. O valor médio dos coeficientes aerodinâmicos e o resultado para o número de Strouhal são relacionados com os valores de referência na Tabela 8. Cabe salientar que o valor de coeficiente de arrasto obtido por Larsen e Walther [41] foi calculado considerando a dimensão característica igual a 31 metros.

Tabela 8. Coeficientes aerodinâmicos e número de Strouhal para a ponte Great Belt East.
Presente trabalho Kuroda [42] Larsen e Walther [41] Braun [43]
CD 0,63 0,49 0,061 0,62
CL -0,20 -0,18 - 0,06
St 0,13 - 0,100 e 0,168 0,18

6. Conclusões

Neste trabalho foi desenvolvida uma ferramenta numérica baseada no esquema CBS (“Characteristics-based split”) para análise de escoamentos incompressíveis e turbulentos no contexto da Engenharia do Vento Computacional, onde o Método dos Elementos Finitos foi empregado na discretização espacial das equações fundamentais do escoamento através de elementos triangulares e tetraédricos. A solução da equação de conservação de massa foi obtida através de um esquema explícito, no qual foi utilizado o método de compressibilidade artificial, e através de um esquema semi-implícito no qual o campo de pressão é obtido de forma implícita. A abordagem LES permitiu a simulação de problemas com escoamentos turbulentos, onde o modelo sub-malha clássico de Smagorinsky foi usado. Os exemplos clássicos apresentados neste trabalho possibilitaram a verificação através de comparações com resultados obtidos por diversas referências, podendo-se afirmar que o esquema numérico implementado mostrou ser capaz de simular escoamentos complexos encontrados na EVC, inclusive com corpos imersos.

Quanto ao esquema CBS-AC, pôde-se observar a grande importância do passo de tempo real, pois ao se adotar valores muito altos o problema não consegue obter a convergência no processo fictício. Quanto mais refinada a malha, menor o passo de tempo local, fazendo com que a utilização de passos reais muito altos leve o esquema a instabilização. Além disso, o esquema demostrou um desenvolvimento lento do campo de pressão devido ao termo transiente fictício (Eq. 25), porém os valores obtidos não são perturbados quanto a precisão, não sendo necessária a utilização da pressão suavizada como no caso de algumas formulações.

Quanto ao esquema semi-implícito, pode-se concluir que a necessidade de um passo de tempo global para todos os elementos acarreta problemas na convergência quando considerado um passo de valor muito alto devido a restrições de estabilidade, embora estas restrições sejam menores do que as observadas em esquemas totalmente explícitos. Por outro lado, o esquema CBS semi-implícito apresentou um desenvolvimento muito mais rápido do escoamento e um tempo de processamento significativamente menor, tornando-o uma alternativa a ser considerada em futuras aplicações na EVC. Em futuros trabalhos, deverá ser implementado um algoritmo de adaptação de malha para um melhor aproveitamento e distribuição de elementos sobre o domínio computacional. Outras variantes do método CBS devem também ser testadas em problemas da EVC, como o método quase-implícito. Para a versão explícita do método, condições iniciais mais adequadas devem ser buscadas para um desenvolvimento mais rápido do escoamento, já que o termo transiente inserido no terceiro passo dificulta a convergência do processo fictício nos primeiros passos reais da simulação. Para um melhor tratamento da turbulência, principalmente junto a superfícies, o modelo dinâmico de Smagorinsky (Lilly [44]) deve ser implementado na presente formulação.

7. Referências

[1] SIMIU, E.; SCANLAN, R. H. Wind Effects on Structures. 3. ed. New York: John Wiley & Sons Inc, 1996.

[2] LÖHNER, R.; HAUG, E.; MICHALSKI, A.; MUHAMMAD, B.; DREGO, A.; NANJUNDAIAH, R.; ZARFAM, R. Recent advances in computational wind engineering and fluid–structure interaction. Journal of Wind Engineering and Industrial Aerodynamics, v. 144, p. 14-23, 2015.

[3] MURAKAMI, S. Current status and future trends in computational wind engineering. Journal of Wind Engineering and Industrial Aerodynamics, vol. 67-68, pp. 3-34, 1997.

[4] MURAKAMI, S. Overview of turbulence models applied in CWE–1997. Journal of Wind Engineering and Industrial Aerodynamics, vol. 74, pp. 1-24, 1998.

[5] STATHOPOULOS, T. Computational wind engineering: Past achievements and future challenges. Journal of Wind Engineering and Industrial Aerodynamics, vol. 67&68, pp. 509-532, 1997.

[6] BLOCKEN, B. 50 years of Computational Wind Engineering: Past, present and future. Journal of Wind Engineering and Industrial Aerodynamics, vol. 129, pp. 69-102, 2014.

[7] MENTER, F. HEMSTROM, B.; HENRIKKSON, M.; KARLSSON, R.; LATROBE, A.; MARTIN, A.; MUHLBAUER, P.; SCHEUERER, M.; SMITH, B.; TAKACS, T.; WILLEMSEN, S. CFD best practice guidelines for CFD code validation for reactor safety applications. Deliverable D01 of the ECORA project, 2002.

[8] MENTER, F. R. Best practice: scale-resolving simulations in ANSYS CFD. Technical Report: ANSYS, v. 1, 2012.

[9] FRANKE, J.; HIRSCH, C.; JENSEN, A. G.; KRÜS, H. W.; SCHATZMANN, M.; WESTBURY, P. S.; MILES, S. D.; WISSE, J. A.; WRIGHT, N. G. Recommendations on the use of CFD in wind engineering. In: Proceedings of the International Conference on Urban Wind Engineering and Building Aerodynamics. COST Action 14. 2004.

[10] FRANKE, J.; HELLSTEN, A.; SCHLÜNZEN, H.; CARISSIMO, B. The COST 732 Best Practice Guideline for CFD simulation of flows in the urban environment: a summary. International Journal of Environment and Pollution, v. 44, n. 1-4, p. 419-427, 2011.

[11] HEINRICH, J. C.; HUYAKORN, P. S.; ZIENKIEWICZ, O. C.; MITCHELL, A. R. An ‘upwind’ finite element scheme for two‐dimensional convective transport equation. International Journal for Numerical Methods in Engineering, vol. 11, pp. 131-143, 1977.

[12] BROOKS, A. N.; HUGHES, T. J. R. Streamline upwind/Petrov-Galerkin formulations for convection-dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, vol 32, pp. 199-259, 1982.

[13] DONEA, J. A Taylor-Galerkin method for convective transport problems. International Journal for Numerical Methods in Engineering, vol. 20, pp. 101-119, 1984.

[14] ZIENKIEWICZ, O. C.; TAYLOR, R. L.; NITHIARASU, P. The Finite Element Method for Fluid Dynamics. 7. ed. Oxford: Elsevier Butterworth-Heinemann, 2014.

[15] BREZZI, F. On the existence, uniqueness and approximation of saddle-point problems arising from Lagrandian multiplies. RAIRO, vol. 8, pp. 129-151, 1974.

[16] CHORIN, A. J. A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics, vol. 2, pp. 12-26, 1967.

[17] CHORIN, A. J. Numerical solution of the Navier-Stokes equations. Mathematics of Computation, vol. 22, pp. 745-762, 1968.

[18] CHORIN, A. J. On the convergence of discrete approximation to the Navier-Stokes equations, Mathematics of Computation, vol. 23, pp. 341-353, 1969.

[19] ZIENKIEWICZ, O. C.; CODINA, R. A general algorithm for compressible and incompressible flow. Part I: The split characteristic-based scheme. International Journal for Numerical Methods in Fluid, vol. 20, pp. 869-885, 1995.

[20] NITHIARASU, P. An efficient artificial compressibility (AC) scheme based on the characteristic based split (CBS) method for incompressible flows. International Journal for Numerical Methods in Engineering, vol. 56, pp. 1815-1845, 2003.

[21] CODINA, R. Pressure stability in fractional step finite element methods for incompressible flows. Journal of Computational Physics, vol. 170, pp. 112-140, 2001.

[22] BEVAN, R. L. T.; NITHIARASU, P. A dual time stepping approach to eliminate first order error in fractional step methods for incompressible flows. International Journal of Numerical Methods for Heat and Fluid Flow, vol. 26, pp. 556-570, 2016.

[23] BEVAN, R. L. T.; BOILEAU, E.; LOON, R. V.; LEWIS, R. W.; NITHIARASU, P. A comparative study of fractional step method in its quasi-implicit, semi-implicit and fully-explicit forms for incompressible flows. International Journal of Numerical Methods for Heat and Fluid Flow, vol. 26, pp. 595-623, 2016.

[24] NITHIARASU, P.; BEVAN, R.; MURALI, K. An artificial compressibility based fractional step method for solving time dependent incompressible flow equations. Temporal accuracy and similarity with a monolithic method. Computational Mechanics, vol. 51, pp. 255-260, 2013.

[25] BRAUN, A. L. Simulação Numérica na Engenharia do Vento Incluindo Efeitos de Interação Fluido-Estrutura. Porto Alegre: UFRGS, 2007. Tese (Doutorado em Engenharia Civil) – Programa de Pós-Graduação em Engenharia Civil, Universidade Federal Rio Grandes Sul, Porto Alegre, 2002.Tese de Doutorado, PPGEC/UFRGS, Porto Alegre, 2007.

[26] SMAGORINSKY, J. General circulation experiments with the primitive equations, I, the basic experiment. Monthly Weather Review, vol. 91, pp. 99-135, 1963.

[27] KNIGHT, D.; ZHOU, G.; OKONG’O; SHUKLA, V. Compressible large eddy simulation using unstructured grids. AIAA JOURNAL, n° 980535, 1998.

[28] SCHLICHTING, H. Boundary Layer Theory. 7. ed. New York: McGraw-Hill Inc, 1979.

[29] MALAN, A. G.; LEWIS, R. W.; NITHIARASU, P. An improved unsteady, unstructured, artificial compressibility, finite volume scheme for viscous incompressible flows: Part I. Theory and implementation. International Journal for Numerical Methods in Engineering, vol 54, pp. 695-714, 2002a.

[30] MALAN, A. G.; LEWIS, R. W.; NITHIARASU, P. An improved unsteady, unstructured, artificial compressibility, finite volume scheme for viscous incompressible flows: Part II. Application. International Journal for Numerical Methods in Engineering, vol 54, pp. 695-714, 2002b.

[31] GHIA, U.; GHIA, K. N.; SHIN, C., T. High-Re solutions for incompressible flow using Navier-Stokes equations and multigrid method. Journal of Computational Physics, vol. 48, pp 387-411, 1982.

[32] RAMASWAMY, B. Theory and implementation of a semi-implicit finite element method for viscous incompressible flow. Computers and Fluids, vol. 22, pp. 725-747, 1993.

[33] ZANG, Y.; STREET, R. L.; KOSEFF, J. R. A Dynamic Mixed Subgrid-Scale Model and its Application to Turbulent Recirculating Flows. Physics of Fluids A, vol. 5, pp. 3186-3196, 1993.

[34] BLEVINS, R. D. Flow-Induced Vibration. 2nd ed. Florida: Krieger Publishing Company, 1990.

[35] WANDERLEY, J. B. V.; SOUZA, G. H. B.; SPHAIER, S. H.; LEVI, C. Vortex-induced vibration of an elastically mounted circular cylinder using an upwind TVD two-dimensional numerical scheme. Ocean Engineering, vol. 35, pp. 1533-1544, 2008.

[36] TONON, P. Simulação Numérica de Escoamentos Incompressíveis Usando a Análise Isogeométrica. Porto Alegre: UFRGS, 2016. Dissertação (Mestrado em Engenharia Civil) – Programa de Pós-Graduação em Engenharia Civil, Universidade Federal do Rio Grande do Sul, Porto Alegre, 2016.

[37] HERJFORD, K. A study of two-dimensional separated flow by a combination of the finite element method and Navier-Stokes equations. PhD thesis, The Norwegian Institute of Technology, Trondheim, Norway, 1995.

[38] ARPINO, F.; CORTELLESSA, G.; DELL’ISOLA, M.; MASSAROTTI, N.; MAURO, A. A robust artificial compressibility (AC) algorithm for transient incompressible flows. In: Second International Conference on Computational Methods for Thermal Problems, 2011.

[39] FERREIRA, A. D.; SILVA, M. C. G.; VIEGAS, D. X.; LOPES, A. G. Wind tunnel simulation of the flow around two-dimensional hills. Journal of Wind Engineering and Industrial Aerodynamics, vol. 38, pp. 109-122, 1991.

[40] MATUELLA, J. L. M. Avaliação em Túnel de Vento do Comportamento da Camada Limite Atmosférica em Terrenos Complexos. Porto Alegre: UFRGS, 2012. Tese (Doutorado em Engenharia Civil) – Programa de Pós-Graduação em Engenharia Civil, Universidade Federal do Rio Grande do Sul, Porto Alegre, 2012.

[41] LARSEN, A.; WALTHER, J. H. Aeroelastic analysis of bridge girder sections based on discrete vortex simulations. Journal of Wind Engineering and Industrial Aerodynamics, vol. 67&68, pp. 253-265, 1997.

[42] KURODA, S. Numerical simulation of flow around a box girder of a long span suspension bridge. Journal of Wind Engineering and Industrial Aerodynamics, vol. 67&68, pp. 239-252, 1997.

[43] BRAUN, A. L. Um modelo para simulação numérica da ação do vento sobre seções de ponte. Porto Alegre: UFRGS, 2002. Dissertação (Mestrado em Engenharia Civil) – Programa de Pós-Graduação em Engenharia Civil, Universidade Federal Rio Grandes Sul, Porto Alegre, 2002.

[44] LILLY, D. K. A proposed modification of the germano subgrid-scale closure method, Phys. Fluids A, vol. 4, pp. 633–635, 1992.

Back to Top