Tutorial Backward Facing Step

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.

enter image description here
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.

Criando um novo caso

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

Variáveis e cálculos nos arquivos OpenFOAM

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.

Definindo a geometria no blockMesh

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.

Listando os vértices

A geometria do Backward Facing Step (BFS) pode ser descrita por meio de três blocos, conforme ilustrado abaixo:

Variáveis e vertíces da geometria

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.

Especificando os blocos e o particionamento (blocks)

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).

Divisão do domínio em blocos hexaédricos

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.

Exemplo de um bloco
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)
);

Especificando a curvatura das arestas

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

Visualizando os blocos

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

Blocos definidos no arquivo blockMeshDict

Categorizando as superfícies de contorno.

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:

OpenFOAM.org User Guide:

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:

Entrada (inlet)

Superfície de entrada

inlet
{
    type patch;
    faces
    (
            (0 8 15 7)
    );
}

Saída (outlet)

Superfícies da Saída (outlet)

outlet
{
    type patch;
    faces
    (
        (12 4 5 13)
        (11 3 4 12)
    );
}

Superior (top)

Fronteira Superior

top
 {
     type wall;
     faces
     (
         (14 6 7 15)
         (13 5 6 14)
     );
 }

Inferior (bottom)

Superfície inferior

bottom
{
    type wall;
    faces
    (
		(2 3 11 10)
		(9 1 2 10)
		(1 9 8 0)
    );
}

Frente e trás (frontAndBack)

Frentre e trás

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)
   );
}

Visualizando as faces criadas

Executando novamente o comando

blockMesh -write-vtk

Agora abrimos o arquivo blockFaces.vtp para visualizar as faces.

Faces definidas no arquivo blockMeshDict

Checando a qualidade da Malha

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.

Visualizando a Malha

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.

Desabilita Volume Fields

Entendendo os arquivos de malha do OpenFOAM

É 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

Configuração das condições de contorno

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)

Pressão (0/p)

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;
}

Velocidade

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 Re=Uh/νRe = U h / \nu, podemos simplificar o problema fazendo h=1h = 1m e ν=1\nu = 1m/s2^2, 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 xx. 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;
}

Configuração dos parâmetros físicos

O único parâmetro físico do modelo matemático utilizado pelo icoFoam é a viscosidade cinemática (ν\nu). Este parâmetro é alterado no arquivo constant/transportProperties. Como mencionado anteriormente faremos a viscosidade igual a 1 m2^2/s.

nu [0 2 -1 0 0 0 0] 1.0;

Esquemas de discretização/interpolação

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,

∂ϕ∂t≈ϕ−ϕoΔt \frac{\partial \phi }{\partial t} \approx \frac{\phi - \phi_o}{\Delta t}

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,

∫V(∇⋅u)dV=∮S(n^⋅u)dS \int_V (\bm \nabla \cdot \bm u) dV = \oint_S (\hat{\bm n} \cdot \bm u) dS

∫V(∇2ϕ)dV=∫V(∇⋅∇ϕ)dV=∮S(n^⋅∇ϕ)dS \int_V (\nabla^2\phi ) dV = \int_V (\bm \nabla \cdot \bm \nabla \phi ) dV = \oint_S (\hat{\bm n} \cdot \nabla \phi) dS

∫V(∇ϕ)dV=∮S(n^ϕ)dS \int_V (\bm \nabla\phi ) dV = \oint_S (\hat{\bm n} \phi) dS

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

Esquemas de solução dos sistemas lineares

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 .

Configuração dos arquivos de saída e controle geral da simulação

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.

C=U∗ΔtΔs<1 C = \frac{U* \Delta t}{ \Delta s } < 1

Portanto,

Δt<ΔSU \Delta t < \frac{ \Delta S}{U}

Onde ΔS\Delta S é o tamanho do intervalo utilizado na discretização espacial.

Como a nossa discretização espacial é uniforme podemos utilizar como ΔS\Delta S = Δx\Delta x = Δy\Delta y

Assim

Δt<0,1124=8×10−4s \Delta t < \frac{0,1}{124} = 8 \times 10^{-4} s

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:

tb=LU=44124 t_b = \frac{L}{U} = \frac{44}{124}

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;

Executando a simulação

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

Visualizando o resultado

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.

Rodando mais passos

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";

Extraindo dados com o Paraview

Extraindo dados com o OpenFOAM

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

∂ux∂y≈u(y+δ)−u(y)δ \frac{\partial u_x}{\partial y} \approx \frac{ u(y + \delta ) - u (y) }{\delta}

Pegando a velocidade na primeira camada de células próxima a parede,

Camada de volumes próximas a parede.
e sabendo que a velocidade na parede é nula:

∂ux∂y≈u(y+Δy/2)Δy/2 \frac{\partial u_x}{\partial y} \approx \frac{ u(y + \Delta y/2 )}{\Delta y / 2}

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 (x−xD)/h(x-x_D)/h versus a (h/U)(∂ux/∂y)(h/U)(\partial u_x/\partial y), onde xDx_D é a posição xx 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