Simulação de cilindros em tandem no OpenFOAM utilizando a biblioteca libAcoustics

Simulação de cilindros em tandem no OpenFOAM utilizando a biblioteca libAcoustics

Observações:

  • Este tutorial tem função meramente didática, visando auxiliar usuários iniciantes, e não representa garantias de qualidade de resultados da simulação. A verificação e validação dos resultados, assim como estudos de refino de malha e demais análises, são de responsabilidade do usuário final.
  • This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software via www.openfoam.com, and owner of the OPENFOAM® and OpenCFD® trade marks.
  • Este tutorial foi criado dentro do contexto de um projeto de extensão da UFSC e não é homologado pelos desenvolvedores dos códigos GMSH, OpenFOAM® e LibAcoustics.
  • Foi utilizada a versão 1912 do OpenFOAM.
  • Este tutorial foi elaborado pelo aluno de graduação Julio Victor Vieira e Widmark Kauê Silva Cardoso sob supervisão de Filipe Dutra da Silva.
  • Ao utilizar este material, solicita-se que sejam dadas as devidas referências à esta página. Dúvidas e sugestões podem ser enviadas aos endereços de contato da página.

Neste tutorial serão apresentados os passos para a simulação de um caso de cilindros em tandem no OpenFOAM juntamente com a biblioteca libAcoustics criada por Epikhin et al., utilizando um caso de dipolo 3D. Esta biblioteca contém a implementação de analogias acústicas que permitem o cálculo do ruído em campo distante a partir de dados aerodinâmicos coletados próximos à fonte. Primeiramente é calculada a pressão acústica no domínio do tempo. Em seguida, esses dados são convertidos para o domínio da frequência utilizando o algoritmo Fast Fourier Transform (FFT). Aqui faremos o uso da analogia de Ffowcs Williams-Hawkings (FWH) e da analogia Curle.

Obs: Antes de instalar a biblioteca libAcoustics, é necessário ter instalado a biblioteca FFTW.

Vamos utilizar como base um caso do próprio OpenFOAM e também um caso da biblioteca de acústica. Vá até o local onde o OpenFOAM está instalado e abra o seguinte diretório: tutorials/incompressible/pisoFoam/RAS. Copie a pasta “cavity” e cole onde você deseja salvar o tutorial com o nome “tutorialCilindros”.

Dentro de “tutorialCilindros” apague a pasta “system”, pois vamos utilizar um modelo de um caso da biblioteca de acústica. Abra, então, a pasta “Tutorials” que está junto com os arquivos de instalação da biblioteca libAcoustics. Abra o caso “dipole3D”, copie a pasta “system” e cole dentro da pasta “tutorialCilindros”.

O primeiro passo é importar a malha gerada previamente no Gmsh (Ver tutorial). Para isso, copie o arquivo “TandemCylinders.msh” e cole dentro da pasta “tutorialCilindros”. Assim, os seguintes arquivos devem ser vistos dentro da pasta:

Apague a pasta “system”, pois vamos utilizar um modelo de um caso da biblioteca de acústica. Abra, então, a pasta “Tutorials” que está junto com os arquivos de instalação da biblioteca libAcoustics (não confundir com a pasta “tutorials” do OpenFOAM). Abra o caso “dipole3D”, copie a pasta “system” e cole dentro da pasta “tutorialCilindros”.

Abra a pasta “constant” e dentro dela crie uma pasta com o nome “triSurface”. Agora copie as superfícies geradas no Gmsh (“superficie1.stl” e “superficie2.stl”) e cole dentro da pasta “triSurface”.

Em seguida, abra o terminal e vá até a pasta do caso. Digite “gmshToFoam TandemCylinders.msh” e aperte “enter”. Esse comando converte o formato da malha do Gmsh para o formato utilizado no OpenFOAM.

Agora digite “checkMesh” no terminal e aperte “enter”. Se tudo estiver certo com a malha, a mensagem “Mesh OK” deve aparecer.

O próximo passo é editarmos os arquivos do modelo copiado para as configurações do caso dos cilindros. Dentro da pasta “tutorialCilindros” abra a pasta “system”. Os arquivos que serão utilizados estão mostrados logo abaixo. Os demais arquivos podem ser excluídos.

Abra, então, o arquivo “controlDict”. Vamos utilizar o solver “pisoFoam”. Faça todas as modificações de acordo com as configurações mostradas abaixo:


 application pisoFoam;
 startFrom latestTime;
 startTime 0;
 stopAt endTime;
 endTime 0.2;
 deltaT 0.25e-4;
 writeControl timeStep;
 writeInterval 1000;
 purgeWrite 0;
 writeFormat ascii;
 writePrecision 6;
 writeCompression on;
 timeFormat general;
 timePrecision 6;
 runTimeModifiable true;
 adjustTimeStep false;
 maxCo 2;
 maxDeltaT 1e-3;

 functions
 {
    #include "curleControl"
    #include "fwhControl"
    #include "soundPressureSampling"
 }

Obs: em “functions” são inseridas as funções da biblioteca de aeroacústica.

Feitas as devidas mudanças, salve o arquivo “controlDict” e feche-o.

Ainda dentro da pasta “system”, abra o arquivo “fvSchemes”. Os seguintes esquemas numéricos serão utilizados (os demais esquemas que aparecem no modelo devem ser apagados):


ddtSchemes
 {
      default backward;
 }

d2dt2Schemes
 {
      default Euler;
 }

gradSchemes
 {
      default Gauss linear;
      grad(p) Gauss linear;
      grad(U) Gauss linear;
 }
 divSchemes
 {
      default Gauss limitedLinear 1;
      div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 laplacianSchemes
 {
      default Gauss linear corrected;
 }
 interpolationSchemes
 {
      default linear;
      interpolate(U) linear;
 }
 snGradSchemes
 {
      default corrected;
 }
 fluxRequired
 {
      default no;
      p ;
      pcorr ;

 }
 wallDist
 {
       method meshWave;
 }

Salve o arquivo “fvSchemes” e feche-o.
Agora, abra o arquivo “fvSolution”. Vamos usar as seguintes configurações:


solvers
{
p
{
 solver PCG;
 tolerance 1e-07;
 relTol 1e-3;
 preconditioner
 {
 preconditioner GAMG;
 tolerance 1e-5;
 relTol 0.01;
 smoother GaussSeidel;
 nPreSweeps 0;
 nPostSweeps 2;
 cacheAgglomeration no;
 nCellsInCoarsestLevel 100;
 agglomerator faceAreaPair;
 mergeLevels 1;
 }

}

pFinal
{
    $p;
    relTol 1e-6;
}
"(U|k|omega)"
{
    solver          smoothSolver;
    smoother        DILUGaussSeidel;
    tolerance       1e-08;
    relTol          0;
}
 } PISO
  {
  nOuterCorrectors 1;
  nCorrectors 3;
  nNonOrthogonalCorrectors 1;
  pRefCell 0;
  pRefValue 0;
  momentumPredictor true;
  }

Salve as modificações feitas e feche o arquivo

Abra o arquivo “decomposeParDict”. Em “numberOfSubdomains” mude o valor para o número de processadores que você deseja utilizar em sua máquina para rodar a simulação em paralelo. Salve a mudança e feche o arquivo.

Agora abra o arquivo “curleControl”. Nele serão definidas as configurações da analogia Curle. Ele deve ser configurado da seguinte maneira:


CurleAnalogy
 {
 libs (“libAcoustics.so”);
 type Curle;
 log true;
 probeFrequency 1;
 patches (cyl-1 cyl-2);
 timeStart 0.16;
 timeEnd 0.2;
 c0 343;
 dRef -1;
 pName p;
 pInf 0;
 rho rhoInf;
 rhoInf 1.225;
 CofR (0 0 0);
 writeFft true;
 U0 (44 0 0);

observers
{
    microphone-A
    {
    position    (-0.4761 1.5896 0.4572);
    pRef        2.0e-5;
    fftFreq        1024;
    }
    microphone-B
    {
    position    (0.5206 1.8288 0.457);
    pRef        2.0e-5;
    fftFreq        1024;
    }
    microphone-C
    {
    position    (1.5173 1.5896 0.457);
    pRef        2.0e-5;
    fftFreq        1024;
    } 
}
 }

Em “patches” são especificadas as superfícies nas quais serão armazenadas os dados de entrada para a analogia de Curle , que nesse caso são as superfícies dos dois cilindros. Para pular o regime transiente, iniciamos a analogia apenas em 0.16 segundos. O “c0” indica a velocidade de propagação do som no ar. O “dRef” se refere à espessura do domínio para simulações 2D, mas como esse caso é 3D, ignoramos essa função alterando o seu valor para um número negativo. Por fim, temos em “observers” os pontos no campo distante em que serão medidos os níveis de ruído. O “fftFreq” indica a quantidade de dados presente em uma janela em que é aplicada a FFT.

Após fazer as modificações, salve e feche o arquivo “curleControl”.

Abra, então, o arquivo “fwhCommonSettings” para definir as configurações da analogia FWH. Vamos utilizar definições semelhantes a da analogia Curle. Para isso, substitua o conteúdo do arquivo pelo seguinte:


 libs ("libAcoustics.so");

log             true;
writeFft        true;
probeFrequency  1;
timeStart       0.16;
timeEnd         0.2;
c0              343;
dRef            -1;
pName        p;
pInf            0;
rho             rhoInf;
rhoInf          1.225;
CofR            (0 0 0);
U0            (44 0 0);

observers
{
    microphone-A
    {
    position    (-0.4761 1.5896 0.4572);
    pRef        2.0e-5;
    fftFreq        1024;
    }
    microphone-B
    {
    position    (0.5206 1.8288 0.457);
    pRef        2.0e-5;
    fftFreq        1024;
    }
    microphone-C
    {
    position    (1.5173 1.5896 0.457);
    pRef        2.0e-5;
    fftFreq        1024;
    }   
}

Salve o arquivo “fwhCommonSettings” e feche-o.

Agora abra o arquivo “fwhControl”. Para a analogia Curle as superfícies utilizadas para a aquisição dos dados são as próprias superfícies dos cilindros (definidas em “patches” no arquivo “curleControl”). Já para a analogia FWH as superfícies a serem utilizadas são aquelas geradas no Gmsh no formato .stl e inseridas na pasta /constant/triSurface. As configurações dessas superfícies são definidas no arquivo “fwhControl”. Substitua seu conteúdo pelo seguinte:


superficie1
 {
 type FfowcsWilliamsHawkings;

#include            "fwhCommonSettings";
patches             (cyl-1 cyl-2);
formulationType     GTFormulation;
cleanFreq           100;

interpolationScheme insideCells;

surfaces
(
    superficie1
    {
        type            sampledTriSurfaceMesh;
        surface         superficie1.stl;
        source          insideCells;
        interpolate     false;
    }

);

     nonUniformSurfaceMotion false;
     Ufwh                    (0 0.0 0.0);

     fixedResponseDelay  true;
     responseDelay       1e-6;
}  
superficie2   
{
  type FfowcsWilliamsHawkings;
  #include            "fwhCommonSettings";
  patches             (cyl-1 cyl-2);
  formulationType     GTFormulation;
  cleanFreq           100;

  interpolationScheme insideCells;

  surfaces
  (
    superficie2
    {
        type            sampledTriSurfaceMesh;
        surface         superficie2.stl;
        source          insideCells;
        interpolate     false;
    }

  );

  nonUniformSurfaceMotion false;
  Ufwh                    (0.0 0.0 0.0);

  fixedResponseDelay  true;
  responseDelay       1e-6;
 }

As superfícies de aquisição em “surfaces” são aquelas dentro da pasta “triSurface” e não as superfícies dos cilindros. O “Ufwh” é utilizado quando se deseja uma superfície móvel, porém nesse caso vamos utilizar superfícies estacionárias.

Feitas as alterações, salve o arquivo “fwhControl” e feche-o.

O último arquivo pertencente à biblioteca é o “soundPressureSampling”. Abra-o e substitua seu conteúdo pelo mostrado abaixo. Nele são configuradas os dados que as superfícies devem coletar e também a forma em que devem ser salvos.


surfacePressure  
{  
type soundPressureSampler;  
libs (“libAcoustics.so”);  
writeControl timeStep;  
writeInterval 1; outputGeometryFormat gmsh;

fields
(
    p
);

pName p;

log   true;

interpolationScheme cellPoint;

surfaces
(
    superficie1
    {
        type        sampledTriSurfaceMesh;
        surface     superficie1.stl;
        source      cells; 
        interpolate true;
    }

    superficie2
    {
        type        sampledTriSurfaceMesh;
        surface     superficie2.stl;
        source      cells; 
        interpolate true;
    }
);
 }

Após aplicar as mudanças, salve e feche o arquivo.

As configurações dos arquivos na pasta “system” estão concluídas. Volte para o diretório do “tutorialCilindros” e abra a pasta “constant” e depois abra o arquivo “transportProperties”. Nele, insira a viscosidade e a densidade do ar, como mostrado abaixo.


transportModel Newtonian;
nu nu [ 0 2 -1 0 0 0 0 ] 1.47e-05;  
rho rho [ 1 -3 0 0 0 0 0 ] 1.225;  
rho [ 1 -3 0 0 0 0 0 ] 1.225;

Ainda na pasta “constant”, abra o arquivo “turbulenceProperties”. O modelo de turbulência a ser utilizado é o “kOmegaSST”.


simulationType RAS;  
RAS  
{  
RASModel kOmegaSST;  
turbulence on;  
printCoeffs on;  
}

Salve as mudanças e feche o arquivo.

Vá até a pasta “polyMesh” e abra o arquivo “boundary”. Vamos alterar o tipo da face dos dois cilindros para parede. Para fazer isso, troque o “type” e “physicalType” de “patch” para “wall” em ambos os cilindros. As demais configurações permanecem inalteradas, como mostradas abaixo.


6  
(  
frontAndBack  
{  
type patch;  
physicalType patch;  
nFaces 39536;  
startFace 569932;  
}  
free  
{  
type patch;  
physicalType patch;  
nFaces 3320;  
startFace 609468;  
}  
inlet  
{  
type patch;  
physicalType patch;  
nFaces 1120;  
startFace 612788;  
}  
outlet  
{  
type patch;  
physicalType patch;  
nFaces 1120;  
startFace 613908;  
}  
cyl-1  
{  
type wall;  
physicalType wall;  
nFaces 560;  
startFace 615028;  
}  
cyl-2  
{  
type wall;  
physicalType wall;  
nFaces 560;  
startFace 615588;  
}  
)

Salve o arquivo “boundary” e feche-o. Retorne para o diretório do “tutorialCilindros”.

Por fim, vamos definir as condições de contorno. Abra a pasta “0”. Os arquivos “epsilon” e “nuTilda” devem ser excluídos, restando apenas os arquivos mostrados logo abaixo.

Substitua o conteúdo de cada arquivo conforme mostrado abaixo.

k” :


dimensions [0 2 -2 0 0 0 0];
internalField uniform 0.00325;

boundaryField

{

inlet

{

type freestream;

freestreamValue $internalField;

}

outlet

{

type zeroGradient;

}

cyl-1

{

type kqRWallFunction;

value uniform 0.00325;

}

cyl-2

{

type kqRWallFunction;

value uniform 0.00325;

}

free

{

type freestream;

freestreamValue $internalField;

}

frontAndBack

{

type zeroGradient;

}

}


nut” :

dimensions [0 2 -1 0 0 0 0]; 
internalField uniform 0; 
boundaryField 
{ 
inlet 
{ 
type calculated; 
value uniform 0; 
} 
outlet 
{ 
type calculated; 
value uniform 0; 
} 
cyl-1 
{ 
type nutUSpaldingWallFunction; 
value uniform 0; 
} 
cyl-2 
{ 
type nutUSpaldingWallFunction; 
value uniform 0; 
} 
free 
{ 
type calculated; 
value uniform 0; 
} 
frontAndBack 
{ 
type zeroGradient; 
} 
}

omega” :


dimensions [0 0 -1 0 0 0 0]; 
internalField uniform 43; 
boundaryField 
{ 
inlet 
{ 
type freestream; 
freestreamValue $internalField; 
} 
outlet 
{ 
type zeroGradient; 
} 
cyl-1 
{ 
type omegaWallFunction; 
value $internalField; 
} 
cyl-2 
{ 
type omegaWallFunction; 
value $internalField; 
} 
free 
{ 
type freestream; 
freestreamValue $internalField; 
} 
frontAndBack 
{ 
type zeroGradient; 
} 
}

Obs.: utilizamos função parede nos cilindros para cada um dos arquivos acima.

p” :


dimensions [0 2 -2 0 0 0 0]; 
internalField uniform 0; 
boundaryField 
{ 
inlet 
{ 
type freestreamPressure; 
freestreamValue uniform 0; 
} 
outlet 
{ 
type fixedValue; 
value uniform 0; 
} 
cyl-1 
{ 
type zeroGradient; 
} 
cyl-2 
{ 
type zeroGradient; 
} 
free 
{ 
type freestreamPressure; 
freestreamValue uniform 0; 
} 
frontAndBack 
{ 
type zeroGradient; 
} 
}

U” :


dimensions [0 1 -1 0 0 0 0];
 internalField uniform (44 0 0); 
boundaryField 
{ 
inlet 
{ 
type freestreamVelocity; 
freestreamValue uniform (44 0 0); 
} 
outlet 
{ 
type zeroGradient; 
} 
cyl-1 
{ 
type noSlip; 
} 
cyl-2 
{ 
type noSlip; 
} 
free 
{ 
type freestreamVelocity; 
freestreamValue uniform (44 0 0); 
} 
frontAndBack 
{ 
type slip; 
} 
}

Obs: Para a pressão na entrada utilizamos a pressão de corrente livre. Para a velocidade, utilizamos a condição de não escorregamento em ambos os cilindros.

Já podemos iniciar a simulação. Abra o terminal e vá até o diretório do caso “tutorialCilindros”. Primeiro vamos decompor o domínio para rodar em paralelo. Digite “decomposePar” no terminal.

Em seguida digite “mpirun -np 5 pisoFoam -parallel > log &” para rodar a simulação (substitua o número 5 pela quantidade de processadores configurados no arquivo “decomposeParDict”).

Obs: A simulação pode demorar algumas horas. Os detalhes do processo vão sendo registrados no arquivo “log” dentro da pasta “tutorialCilindros”.

Pós-processamento

Terminada a simulação, digite “reconstructPar -latestTime” no terminal para obter o último passo da simulação.

Para verificar o valor de y+, digite “pisoFoam -postProcess -func yPlus”.

Antes de abrir os resultados no Paraview, é necessário criar um arquivo com a extensão .foam. Para isso, digite “nano cilindros.foam” e aperte “enter”. Em seguida aperte “ctrl + o”, “enter” e por fim “ctrl + x”.

Para abrir o Paraview digite “paraview cilindros.foam

No Paraview, clique em Apply. Para visualizar o campo de pressão, primeiro selecione o último passo de tempo e também clique na opção “Rescale to data range” como indicado na figura abaixo.

Já o campo de velocidade tem o seguinte resultado:

Os dados acústicos encontram-se na pasta “acousticData” dentro da pasta “tutorialCilindros”. Para plotar os gráficos relacionados aos dados acústicos será utilizado aqui um código em Matlab.

Primeiramente é necessário importar os dados. Com o Matlab aberto vá até a aba “Home” e clique em “ImportData”. Vá até o diretório do “tutorialCilindros” e abra a pasta “acousticData”. Vamos plotar os resultados para o microfone A. Selecione então o arquivo “fftCurleAnalogymicrophoneA” e clique em “Open”. Na janela que abrir clique na opção ”ImportSelection”. Repita o mesmo procedimento para abrir os arquivos “fft-superficie1-microphone-A” e “fft-superficie2-microphone-A”.

Para plotar o gráfico comparativo entre os três resultados, isto é, utilizando a analogia de FWH com superfícies 1 e 2, e também utilizando a analogia Curle, é executado o seguinte código:


% Frequências
fCurle = fftCurleAnalogymicrophoneA{:,1};
fS1 = fftsuperficie1microphoneA{:,1};
fS2 = fftsuperficie2microphoneA{:,1};

% Nível de Pressão Sonora
splCurle = fftCurleAnalogymicrophoneA{:,3};
splS1 = fftsuperficie1microphoneA{:,3};
splS2 = fftsuperficie2microphoneA{:,3};

% Espectro de frequência
semilogx(fCurle,splCurle,’-r’)
hold on
semilogx(fS1,splS1,’-b’)
hold on
semilogx(fS2,splS2,’-k’)
hold on
xlabel(“Frequência (Hz)”)
ylabel(“NPS (dB)”)
xlim([50 500])
legend(“Curle”, “FWH-S1”, “FWH-S2”)


Dessa forma, o resultado obtido é:

Obs: Para obter melhores resultados pode-se aumentar o número de pontos por bloco da FFT (nfft), tanto no arquivo “fwhCommonSettings” quanto no arquivo “curleControl”. Porém ao aumentar o nfft também pode ser necessário aumentar o tempo de simulação. Vale ressaltar, ainda, que a malha utilizada tem fins meramente didáticos e que se o objetivo for obter resultados próximos de valores experimentais, refinamentos no modelo de simulação devem ser aplicados.