Este tutorial apresenta o passo a passo para construir um caso do OpenFOAM para o escoamento conhecido como Backward Facing Step.
Esse problema consiste em um problema bidimensional de um escoamento entre placas planas onde uma entrada (inlet) com um perfil parabólico encontra um degrau, conforme esquema mostrado na figura abaixo.
Este é um problema amplamente estudado na literatura tanto do ponto de vista numérico quanto experimental. Nosso intuito será reproduzir um conjunto de resultados experimentais publicados no artigo
TIHON, Jaroslav; PĚNKAVOVÁ, V.; PANTZALI, M. The effect of inlet pulsations on the backward-facing step flow. European Journal of Mechanics-B/Fluids, v. 29, n. 3, p. 224-235, 2010.
No openFOAM, um novo caso consiste em criar uma pasta com um conjunto de arquivos de configuração específicos organizados em uma estrutura de subpastas, mostradas abaixo.
Os arquivos e as opções contidas nestes arquivos variam de acordo com o solver utilizado. Apesar de ser possível criar estas pastas e arquivos do zero (from scracth), o modo mais prático é copiar um tutorial para o solver em questão, alterando as opções desejadas.
Como faremos uso do icoFoam faremos uma cópia do tutorial cavity através do comando foamCloneCase
foamCloneCase $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity $FOAM_RUN/bfs
O icoFoam é um solver para escoamentos incompressíveis sem modelo de turbulência, ou seja, utiliza as equações de Navier-Stokes para escoamentos incompressíveis. Mais informações sobre o solver podem ser encontradas no endereço:
OpenFoam.com User Guide: icoFoam
Possuindo a estrutura de arquivos passaremos a configuração do nosso caso que englobará:
Assim o primeiro passo é acessar a pasta criada:
cd $FOAM_RUN/bfs
Os arquivos de configuração do OpenFOAM permitem a criação de variáveis, assim como a realização de cálculos no interior dos arquivos o que auxilia na parametrização dos arquivos e no fácil estudo da variação de parâmetros.
Para declarar uma variável basta especificar uma entrada simples. Por exemplo para criar uma variável Re (número de Reynolds) com o valor 1500:
Re 1500.0;
Para utilizar a variável em algum campo da configuração do openfoam utiliza-se então o nome da variável precedido do símbolo de crifão ($). Por exemplo, para utilizar a variável definida acima utiliza-se $Re.
Adicionamente, o OpenFOAM permite a realização mos operações algébricas básicas assim como utilizar funções matemática básicas. Para realizar um cálculo em alguma parte do arquivo utiliza-se a sintaxe #eval “expressão”.
Por exemplo para calcular a viscosidade cinemática a partir do número de Reynolds, da velocidade e de um comprimento
Re 1500.0;
L 20.0;
U 1.0;
nu #eval "$U * $L / $Re";
Para maiores informações sobre a sintaxe e operações que a diretriz #eval aceita acessar a documentação no link abaixo.
https://www.openfoam.com/documentation/guides/v2112/doc/openfoam-guide-expression-syntax.html
O openfoam também possui uma diretriz chamada #calc que aceita o uso de funções importadas da linguagem C. Adicionalmente, é possível inserir código de várias linhas na linguagem C dentro dos arquivos (#codestream) , o que permite, por exemplo, a criação de funções e a utlização de estruturas for, while e if.
O blockMesh é uma ferramenta para geração de malhas estruturadas onde a geometria é especificada dividindo o domínio em hexaedros, sendo que as arestas destes podem ser curvas.
A especificação da geometria é realizada no arquivo system/blockMeshDict. Tal arquivo pode ser editado usando o nano
nano system/blockMeshDict
Editaremos tal arquivo realizando na seguinte sequência a:
Antes de iniciarmos a especificação da geometria, vamos alterar o fator de escala para que a unidade utilizada na geração da malha seja o metro alterando a variável scale para 1
scale 1;
Este parâmetro multiplicará todas as distância contidas neste arquivo durante o processo de geração da malha. Ou seja, se scale for 0.1, todos os comprimentos definidos no arquivo blockMeshDict serão dividos por 10. O objetivo deste parâmetro é facilitar a conversão de unidades da geometria.
Obs: Historicamente este parâmetro se chamava ConvertToMeters, mas foi renomeado pois seu nome levava a interpretações errôneas sobre seu significado.
A geometria do Backward Facing Step (BFS) pode ser descrita por meio de três blocos, conforme ilustrado abaixo:
Para facilitar a especificação dos vértices podemos criar variáveis conforme a figura. Logo adicionamos após o cabeçalho do arquivo
h 1;
l1 4;
l2 40;
L #eval "$l1+$l2";
H #eval "2*$h";
dz 0.1;
Onde dz é aprofundidade da geometria no eixo z, o que é um valor arbitrário, haja vista que a simulação é na bidimensional na prática.
Na lista deminada vertices acrescentamos os vertices conforme numeração e descrição da figura:
vertices
(
(0 $h 0) // vertex number 0
($l1 $h 0) // vertex number 1
($l1 0 0) // vertex number 2
($L 0 0) // vertex number 3
($L $h 0) // vertex number 4
($L $H 0) // vertex number 5
($l1 $H 0) // vertex number 6
(0 $H 0) // vertex number 7
(0 $h $dz) // vertex number 8
($l1 $h $dz) // vertex number 9
($l1 0 $dz) // vertex number 10
($L 0 $dz) // vertex number 11
($L $h $dz) // vertex number 12
($L $H $dz) // vertex number 13
($l1 $H $dz) // vertex number 14
(0 $H $dz) // vertex number 15
);
Note que os arquivos do OpenFOAM seguem os padrões da linhagem C, logo o símbolo // representa comentários no arquivo.
O próximo passo é a definição dos blocos a partir dos pontos dados. Nossa geometria será dividida em três blocos hexaédricos. Isto se faz necessário pois a face de um bloco deve coincidir com a fase de um outro bloco (ou faze parte do contorno do problema).
Um bloco é especificado a partir de uma lista de vértices, com uma ordem espefíca. Por ser um hexaedro, cada bloco será constituido de 8 vértices. Especifica-se primeiro os 4 vértices que formam uma das face do bloco, depois específica-se os 4 vértices da face diametralmente oposta, seguindo a mesma sequência e começando com o vértice da segunda face que está conectado ao primeiro vertice especificado da primeira face.
Na figura acima escolhemos especificar os dois planos destacados começando pelo plano azul e pelo vértice a, seguindo o sentido a b c d. O plano verde (diametralmente oposto ao azul) deve ser especificado no mesmo sentido e começando do ponto e pois o mesmo está conectado com o ponto inicial a do plano azul. Assim o segundo plano fica e f g h. Assim esse bloco ficaria especificado na forma:
a b c d e f g h
Note que poderíamos escolher outros dois plano e começar de pontos diferentes, logo não há um modo único de realizar tal descrição. Adicionalmente, a ordem de especificação definirá os eixos. A aresta que liga o primeiro ao segundo vértice (ab no exemplo) será o eixo 1, a aresta que liga o segundo ao terceiro vértice (bc no exemplo) será o eixo 2 e a aresta que liga o primeiro a quinto vértice ( ae no exemplo) será o eixo 3.
Após definir o bloco, e consequentemente os seus eixos, definimos em quantos intervalos cada eixo será divido para formar as células da malha.
Também podemos definir qual a razão entre o tamanho dos invervalos. Por exemplo, se a razão for 0.5 então o segundo intervalo será metade do primeiro e assim por diante, seguindo uma progressão geométrica de razão 0.5.
A especificação do bloco é feita no interior da lista blocks por uma entrada simples no formato
hex (a b c d e f g h) (N1 N2 N3) simpleGrading (r1 r2 r3)
A palavra hex especifica o formato do bloco (único aceito até o momento pelo blockMesh), seguido pela lista dos vértices entre parênteses. Os valores N1, N2, N3 são respectivamente o número de divisões nos eixos 1, 2 e 3 do bloco. A palavra simpleGrading especifica o método de divisão do bloco (único aceito até o momento pelo blockMesh) seguido pela razões entre os tamanhos dos intervalos em cada eixo entre parênteses.
É importante notar que deve haver consistência na divisão dos blocos, isto é, uma aresta que pertence a dois blocos diferentes deve ser particionada da mesma maneira nos dois blocos.
Para especificar a divisão nos blocos de forma alterá-las facilmente em futuras simulações definimos as seguintes variáveis
dx 0.1;
dy 0.1;
Nh #eval "10 ";
N1 #eval "$l1 / $dx ";
N2 #eval "$l2 / $dx ";
Desse modo, para nossa geometria temos a seguinte especificação dos blocos
blocks
(
hex (0 1 6 7 8 9 14 15) ($N1 $Nh 1) simpleGrading (1 1 1)
hex (1 4 5 6 9 12 13 14) ($N2 $Nh 1) simpleGrading (1 1 1)
hex (2 3 4 1 10 11 12 9) ($N2 $Nh 1) simpleGrading (1 1 1)
);
Os blocos podem ter arestas curvas, o que permite criar geometria mais complexas. Isso é realizado por meio da lista edges.
O padrão do blockMesh é considerar as arestas como retas. Logo para o nosso problema, onde as arestas são linhas retas, este o bloco ficará em branco.
Caso quiséssemos que a aresta que ligas os vértices a e b for curva basta especificar uma entrada no seguinte forma no interior da lista edges:
a b curva (parâmetros)
Os tipos de curvas e parâmetros são apresentados na tabela abaixo:
Curva | Description | Parâmetros |
---|---|---|
arc | Arco de um Círculo | Um ponto do círculo |
spline | Spline cúbica | Lista de Pontos de interpolação |
polyLine | Conjunto de linhas | Lista de Pontos que definem as linhas |
BSpline | Curva B-Spline | Pontos de Controle da B-Spline |
line | Linha reta |
Antes de criarmos a malha podemos verificar se os blocos estão corretamente definidos. Isto pode ser realizado utilizando o paraFoam (paraview)
blockMesh -write-vtk
Este comando gerará dois arquivos blockTopology.vtu e blockFaces.vtp que podem ser abertos usando o paraview e mostram a estrutura de blocos e a faces definidas. Para visualizar os blocos abra o arquivo blockTopology.vtu
Para que permitir a aplicação de condições de contorno na malha gerada é preciso separar as superfícies do blocos em conjuntos. Depois da malha gerada configuraremos as condições de contorno por meio deste conjuntos.
Isto é realizado no interior da lista boundary por meio de dicionários com o seguinte formato:
nome_do_conjunto
{
type tipo_da_superficie;
faces
(
(a b c d)
);
}
As superfícies podem ser dos seguintes tipos:
Outros tipos superfícies podem ser especificados, maiores informações podem ser encontradas no link abaixo:
Para o problema BFS, temos condição de entrada (inlet) e saída (outlet) que são do tipo patch. A parte superior (top) e inferior (bottom) que são do tipo wall. Já a parte da frente e de trás (frontAndBack) são do tipo empty
Cada face de cada bloco que não é uma face interna (entre blocos) deve ser especificada. A face pode ser especificada partindo que qualquer vértice da face mas deve ser uma ordem tal que os vérticies são percorridos no sentido horário se estamos olhando para a face estando posicionados no interior do bloco.
Assim ficamos com os seguintes contornos:
inlet
{
type patch;
faces
(
(0 8 15 7)
);
}
outlet
{
type patch;
faces
(
(12 4 5 13)
(11 3 4 12)
);
}
top
{
type wall;
faces
(
(14 6 7 15)
(13 5 6 14)
);
}
bottom
{
type wall;
faces
(
(2 3 11 10)
(9 1 2 10)
(1 9 8 0)
);
}
frontAndBack
{
type empty;
faces
(
(0 7 6 1)
(6 5 4 1)
(1 4 3 2)
(13 14 9 12)
(14 15 8 9)
(12 9 10 11)
);
}
Executando novamente o comando
blockMesh -write-vtk
Agora abrimos o arquivo blockFaces.vtp para visualizar as faces.
O OpenFOAM fornece um aplicativo denominado checkMesh para analisar a qualidade da malha gerada.
Executando
checkMesh
Tem-se informações como razão de aspecto, não ortogonalizada, skewness entre outros.
Podemos fazer uma inspeção visual na malha gerada utilizando o paraview (ou o paraFoam).
touch bfs.foam
paraview bfs.foam
Porém é necessário observar que ainda não especificamos as condições de contorno e isso pode acarretar um erro. Assim, antes de pressionar apply é necessário desabilitar a visualização dos campos p e U na seção volume fields e também colocar a representação (representation) como Wireframe, para uma melhor visualização da malha.
É importante notar que o blockMesh é um gerador de malhas e que diversos geradores de Malhas (seja do OpenFOAM ou externo) podem ser utilizados em conjunto com o OpenFOAM.
Quando o blockMesh é executado com sucesso, ele cria uma série de arquivos na pasta constant/polyMesh que descrevem todos as células, vertices e faces da malha. Uma breve visualização do conteúdo destes arquivos lhe dará uma idéia de como a malha é especificada em seu interior.
boundary faces neighbour owner points
Arquivo | Descrição |
---|---|
constant/polyMesh/boundary | Lista das faces que pertencem a cada tipo de contorno. |
constant/polyMesh/faces | Lista com todas as faces das células/volumes. |
constant/polyMesh/neighbour | Descreve a relação de vizinhança entre as células. |
constant/polyMesh/owner | Relaciona as faces e as células |
constant/polyMesh/points | Lista dos os nós/vértices da malha |
Tendo a malha gerada e as superfícies de contorno devidamente nomeados passamos a definir as condições de contorno.
Estas definições nos arquivos 0/<nome do campo>, onde os campos, no caso do solver icoFOAM, são p e U.
Nestes arquivos também são impostas as condições iniciais (*internalField *), que para o nosso caso escolheremos velocidade nula e pressão nula em todo o domínio.
As condições de contorno são então impostas no dicionário boundaryField
desses arquivos. No formato,
nome_do_conjunto
{
type tipo_da_condição;
opções_do_tipo
}
Uma descrição detalhada dos diferentes tipos de condição de contorno podem ser encontradas na documentação do OpenFOAM.
Basic Boundary Conditions (OpenFOAM Foundation User Guide)
Boundary Conditions (OpenCFD Ltd User Guide)
As condições de contorno para a pressão são definidas no arquivo (0/p), assim
nano 0/p
Na superfície de entrada (inlet) queremos impor fluxo constante de entrada de fluido, com uma velocidade prescrita, usaremos neste caso então um gradiente nulo de pressão. O mesmo ocorrerá nas superfícies da parte superior (top) e inferior (bottom), onde usaremos a condição de não escorregamento, ou seja, de velocidade nula na superfície.
inlet
{
type zeroGradient;
}
top
{
type zeroGradient;
}
bottom
{
type zeroGradient;
}
Na superfície de saída (outlet) iremos impor uma pressão fixa, tentando emular uma condição de escoamento plenamente desenvolvido.
outlet
{
type fixedValue;
value uniform 0;
}
Para as superfícies na frente e atrás (frontAndBack) só precisamos especificar que elas são fazias (empty)
frontAndBack
{
type empty;
}
As condições de contorno para a velocidade são definidas no arquivo (0/U), assim
nano 0/U
Na entrada (inlet) queremos impor um fluxo volumético prescríto. Realizaremos impondo um campo de velocidade uniforme na direção x. Antes definiremos uma variável para tornar fácil o controle dessa condição
Re 124.0;
Ui $Re;
Como , podemos simplificar o problema fazendo m e m/s, assim a velocidade será numericamente igual ao número de Reynolds porém com a unidade de m/s. Assim a condição de contorno para a entrada (inlet) pode ser inserida no dicionário boundaryField na forma:
inlet
{
type fixedValue;
value uniform ($Ui 0 0);
}
Na saída queremos impor uma condição de escoamento plenamente desenvolvido, ou seja, que não há mais variação do campo de velocidade na direção . Para isso impomos uma condição de gradiente nulo:
outlet
{
type zeroGradient;
}
Como mencionado anteriormente, nas superfície superior e inferior queremos impor uma condição de não escorregamento, ou seja, velocidade nula. Isto poderia ser alcançado com a condição do tipo fixedValue, porém devido a seu uso freqüente, os desenvolvedores do OpenFOAM criaram o tipo noSlip.
top
{
type noSlip;
}
bottom
{
type noSlip;
}
Para as superfícies na frente e atrás (frontAndBack) só precisamos especificar que elas são fazias (empty)
frontAndBack
{
type empty;
}
O único parâmetro físico do modelo matemático utilizado pelo icoFoam é a viscosidade cinemática (). Este parâmetro é alterado no arquivo constant/transportProperties. Como mencionado anteriormente faremos a viscosidade igual a 1 m/s.
nu [0 2 -1 0 0 0 0] 1.0;
O arquivo system/fvSchemes contém a configuração sobre qual esquema de discretização é utilizado para transformar as diferentes operações diferenciais em operações numéricas.
nano system/fvSchemes
Para o caso em questão não faremos alterações neste arquivo mas é interessante fazer algumas observações.
O dicionário ddtSchemes se refere ao esquema de discretização no tempo, sendo o método de Euler dada por,
Informação sobre outros esquemas de discretização da derivada temporal podem ser encontrados em
Time Schemes - OpenCFD Ltd User Guide
Quanto aos esquemas de discretização das derivadas espacias como gradientes (gradSchemes), divergentes (divSchemes) e laplacianos (laplacianSchemes) o método principal é a utilização do teorema de Gauss,
Como os valores das variáveis só são conhecidos no centro das células é necessário utilizar um método de interpolação. Para este caso em questão utilizaremo uma interpolação linear. Maiores informações sobre outras possibilidades de esquemas de discretização espacial podem ser encontrados em
Grad Schemes - OpenCFD Ltd User Guide
Divergence Schemes - OpenCFD Ltd User Guide
Laplacian Schemes - OpenCFD Ltd User Guide
Como discutido em nossa Aula 0, o método dos volumes finitos aproximará o conjunto de Equações Diferenciais Parciais do modelo matemático por um sistema linear, cuja solução fornece o valor das variáveis no centro das células.
Assim o ponto de maior demanda computacional é de fato a resolução destes sistemas lineares. A definição dos métodos utilizados para obtenção desta solução é realizada no arquivo (system/fvSolution)
Neste caso também não faremos nenhuma alteração neste arquivo. Porém nota-se que podemos definir desses arquivos o método numérico para solução do sistema linear (solver)
um algorítimo de pré-condicionamento do sistema (preconditioner)
uma tolerância (tolerance), uma tolerância relativa (relTol).
Note que opta-se por um método para cada uma das variáveis (no caso p e U), ou seja, o algoritmo assume cada equação gera um sistema linear. Como as duas variáveis são acopladas, é necessário definir um esquema de acoplamento pressão-velocidade. Em geral cada aplicativo tem um esquema pré-definido, podendo o usuário alterar somente os seus parâmetros .
As configurações sobre a saída da simulação (sobre a escrita de arquivos com os campos U e p) assim como o controle geral do tempo de simulação e do passo de tempo será realizada pelo arquivo (system/controlDict).
Os parâmetros são em sua maioria autoexplicativos. Alteraremos aqui o passo de tempo e o tempo final. O passo de tempo deve ser calculado para que o número de Courant seja menor que 1.
Portanto,
Onde é o tamanho do intervalo utilizado na discretização espacial.
Como a nossa discretização espacial é uniforme podemos utilizar como = =
Assim
Para definir um tempo final podemos usar como referência o tempo que o fluido leva para atravessar a geometira usando a velocidade de entrada:
Para finalizar podemos editar o arquivo system/controlDict
nano system/controlDict
Definindo então os parâmetros tempos
dtCo 8E-4;
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime #eval "44/124";
deltaT #eval "0.5*$dtCo";
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
Podemos executar a simulação com o comando
icoFoam
A saída no console mostrará algumas informações como o número de Courant Máximo, resíduo, numero de iterações para convergência, para cada passo da simulação. É interessante salvar essa informação em um arquivo de log, isso pode ser realizado executando o icoFoam na: forma:
icoFoam > bfs.log
Podemos visualizar o resultado com o paraFoam
paraFoam
Ou diretamente com o paraview
touch bfs.foam
paraview bfs.foam
Notamos através da simulação, que o regime permanente ainda não atingido, sendo necessário calcular mais passos na simulação.
Se a simulação ainda não atingiu um regime permanente, é possível alterar o parâmetro controlDict para ele utilizar o último tempo disponível como tempo de inicio e alterar o tempo final da simulação, a fim de continuar a simulação. Ou seja, o OpenFOAM procura as pastas de tempo e utilizará como condição de inicio da simulação a com a pasta do último tempo já simulado.
Isso pode ser realizado mudando o parâmetro startFrom no arquivo controlDict.
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime #eval "2 * 44 / 124";
O openFoam permite a utilização de funções para cálculos de pós-processamento como a determinação do coeficiente de arrasto, cálculo da tensão, entre outros. Tais funções também pode ser utilizadas em tempo de execução, ou seja, os cálculos são realizados e os resultados são salvos em arquivos a medida que a simulação vai sendo executada.
A forma recomendada é configurar tais funções em arquivos textos individuais, colocados na pasta system. Caso os cálculos sejam realizados após a execução da simulação, com base nos resultados salvos, utiliza-se este arquivo de configuração com o aplicativo postProcess. Quando deseja-se que os cálculos seja realizado durante a simulação este arquivo de configuração deve ser importado para o dicionário functions do arquivo controlDict.
Estas funções são chamadas de functions objects e uma ajuda (bastante limitada) sobre sua configuração pode ser encontrada no user guide do OpenFOAM
Function Objects - User Guide OpenCFD Ltd
No caso precisamos calcular a taxa de cisalhamento nas paredes superior e inferior. Isso pode ser aproximado numericamente fazendo
Pegando a velocidade na primeira camada de células próxima a parede,
e sabendo que a velocidade na parede é nula:
Assim precisamos extrair a velocidade na primeira camada de sítios na parte inferior e na parte superior. Isso pode ser realizado usando a função do tipo sets.
Para realizar tal operação criamos um arquivo sample na pasta system
nano system/sample
Nele configuramos a functionObject desejada que neste caso é a sets que faz parte da biblioteca sampling:
sample
{
type sets;
libs ("libsampling.so");
writeControl timeStep;
writeInterval 20;
interpolationScheme cellPoint;
setFormat csv;
sets
(
bottom
{
type uniform;
axis x;
start (4.05 0.05 0.05);
end (43.95 0.05 0.05);
nPoints 400;
}
);
fields ( U );
}
Para que esta função seja executada durante a simulação ela deve ser incluída no arquivo controlDict, em um dicionário denominado functions que pode ser inserido no final do arquivo.
functions
{
#includeFunc sample;
}
Existem duas maneiras de utilizar a função sample recém criada via o aplicativo postProcess ou durante a própria execução do programa.
Durante a simulação o solver irá executar por padrão todas as functions incluídas no controlDict. Para instruir o solver a não executar as functions objects deve-se adicionar a flag noFunctionObjects
icoFoam -noFunctionObjects
As functions objects também pode ser executadas após o termino da simulação com os dados salvos, através do comando
postProcess -func name
Para comparação dos resultados usuaremos os dados da tensão de cisalhamento podem ser encontradas no artigo,
Tihon, Jaroslav, V. Pěnkavová, and M. Pantzali. "The effect of inlet pulsations on the backward-facing step flow."European Journal of Mechanics-B/Fluids 29.3 (2010): 224-235.
Para facilitar a comparação os dados foram já extraídos e podem ser baixados com o comando
wget https://filedn.com/lB6c5SdGE1nuKljnQmnY2YS/Bottom_less.csv
Neste arquivo são colocadas a posição adimensional versus a , onde é a posição do degrau.
Para a comparação dos resultados usaremos aqui o programa gnuplot que é bastante simples e pode ser encontrando em quase toda distrubuição linux.
gnuplot
O gnuplot funciona através de comandos. Para fazer o gráfico comparando dados de simulação e
set term png
set out "resultado.png"
set datafile separator ','
plot "postProcessing/sample/<Last Time Step>/bottom_U.csv" using ($1-4):($2/(0.05*124)) , "Bottom_less.csv" using 1:2
set out
exit