Para explicar a estrutura e o algoritmo de solução de um código no OpenFOAM, vou usar como base o solver laplacianFoam. Para acompanhar melhor esse post, é interessante que o leitor tenha alguns conhecimentos básicos da sintaxe de C++. Porém, vou apresentar alguns detalhes referentes aos comandos e funções que são membros das classes e templates, facilitando a leitura do código para os leitores sem experiência em linguagens orientadas a objetos.
O solver laplacianFoam é usado para resolver o problema da difusão pura de um campo escalar T, sem considerar nenhum termo fonte. Esta equação está colocada abaixo, sendo D o coeficiente de difusão.
Os arquivos referentes aos
solvers do OpenFOAM ficam no diretório
OpenFOAM-version/applications/solvers., onde
version se refere a versão do OpenFOAM. O código do
laplacianFoam fica no diretório
basic e está colocado abaixo.
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "fvCFD.H"
00034
00035
00036
00037
00038 int main(int argc, char *argv[])
00039 {
00040
00041 # include "setRootCase.H"
00042
00043 # include "createTime.H"
00044 # include "createMesh.H"
00045 # include "createFields.H"
00046
00047
00048
00049 Info<< "\nCalculating temperature distribution\n" << endl;
00050
00051 for (runTime++; !runTime.end(); runTime++)
00052 {
00053 Info<< "Time = " << runTime.timeName() << nl << endl;
00055 # include "readSIMPLEControls.H"
00056
00057 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
00058 {
00059 solve
00060 (
00061 fvm::ddt(T) - fvm::laplacian(DT, T)
00062 );
00063 }
00064
00065 # include "write.H"
00066
00067 Info<< "ExecutionTime = " << runTime.elapsedCpuTime() <<" s"
00068 << " ClockTime = " << runTime.elapsedClockTime() <<" s"
00069 << nl << endl;
00070 }
00071
00072 Info<< "End\n" << endl;
00074 return(0);
00075 }
00076
00077
00078
Para acompanhar a leitura, sugiro que você abra diretamente o arquivo
laplacianFoam.C em um editor de textos de sua preferência. De forma mais direta, você também acessar o código do
laplacianFoam abrindo
este link em uma outra aba do seu navegador.
A primeira linha a ser executada no código (linha 33) declara a biblioteca
fvCFD.H, que fornece ao
solver acesso a todas as classes e propriedades descritas em
meu post anterior. Deve-se frisar que é vital declarar esta biblioteca. Na linha 38, a função
main engloba todo o código fonte principal e possui dois argumentos de entrada: o inteiro
argc e a cadeia de caracteres
argv. Estes parâmetros contém informações sobre a simulação, como o diretório e o nome do caso a ser simulado. Os argumentos são lidos pelo programa diretamente na linha de comando para execução do
solver.
A biblioteca
setRootCase.H é usada para testar a validade dos argumentos
argc e
argv da simulação. O conteúdo desta biblioteca está colocado abaixo.
00001 Foam::argList args(argc, argv);
00002
00003 if (!args.checkRootCase())
00004 {
00005 Foam::FatalError.exit();
00006 }
O primeiro comando do código acima declara a variável
args construída com os argumentos
argc e
argv, a partir da classe
argList. Em seguida, o comando
checkRootCase() verifica a validade e a existência do diretório e do nome do caso simulado. Caso o retorno de
checkRootCase() seja
False, a execução do
solver é interrompida pelo comando padrão de erro do OpenFOAM
FatalError.
As duas próximas bibliotecas declaradas,
createTime.H e
createMesh.H, são responsáveis pela criação de bancos de dados para armazenar dados sobre o caso simulado e a estrutura da malha utilizada. O código colocado abaixo refere-se à biblioteca
createTime.H.
00001 Foam::Info<< "Create time\n" << Foam::endl;
00002
00003 Foam::Time runTime
00004 (
00005 Foam::Time::controlDictName,
00006 args.rootPath(),
00007 args.caseName()
00008 );
Sendo construída a partir de informações sobre o nome e o diretório do caso simulado (provindas de
args), a variável
runTime do template
Time é declarada no código acima. Desta forma,
runTime obtém a localização do arquivo de configuração
controlDict do caso e, utilizando as informações contidas neste último, monta um banco de dados para controle da simulação. Por exemplo, pode-se citar alguns dos dados contidos em
runTime: (i) instante inicial e final; (ii) controle do passo de tempo (fixo, adaptativo, etc.); (iii) diretórios que contém os arquivos com os campos iniciais das propriedades transportadas; (iv) controle de escrita em arquivo (formato de saída, compressão de dados, etc.); entre outros.
A biblioteca
createMesh.H usa o template
fvMesh para declarar a variável
mesh, construída a partir de outro template chamado
IOobject, como mostra o código a seguir.
00001 Foam::Info<< "Create mesh for time = "
00002 << runTime.timeName()<< Foam::nl << Foam::endl;
00003
00004 Foam::fvMesh mesh
00005 (
00006 Foam::IOobject
00007 (
00008 Foam::fvMesh::defaultRegion,
00009 runTime.timeName(),
00010 runTime,
00011 Foam::IOobject::MUST_READ
00012 )
00013 );
O template IOobject define os atributos de um objeto de modo a fornecer meios para entrada e/ou saída de dados (usualmente em arquivo). Com as informações de runTime, o template fvMesh é capaz de localizar os arquivos cells, faces, points e boundary para construção da malha. Note que o último parâmetro na construção de IOobject refere-se à regras de leitura e escrita de arquivos. Neste caso, a leitura dos dados deve ser obrigatória (MUST_READ).
Todas as bibliotecas supracitadas são gerais e podem ser usadas em qualquer código do OpenFOAM (salvo pequenas modificações já comentadas). Contudo, a biblioteca
createFields.H é usada para criação e leitura dos campos iniciais das incógnitas do problema e leitura de propriedades físicas aplicadas a cada caso (propriedades de transporte, termodinâmicas, termofísicas, etc.). Desta forma, este
header é específico para cada
solver e deve ser desenvolvido com cuidado pelo programador, pois todas as incógnitas e todas as propriedades físicas do problema devem ser definidas neste arquivo. Portanto, deve-se ter domínio do modelo fluidodinâmico e qual a melhor forma de armazenar suas variáveis de modo a otimizar o código. A biblioteca
createFields.H específica para o
laplacianFoam está colocada abaixo.
00001 Info<< "Reading field T\n" << endl;
00002
00003 volScalarField T
00004 (
00005 IOobject
00006 (
00007 "T",
00008 runTime.timeName(),
00009 mesh,
00010 IOobject::MUST_READ,
00011 IOobject::AUTO_WRITE
00012 ),
00013 mesh
00014 );
00015
00016
00017 Info<< "Reading transportProperties\n" << endl;
00018
00019 IOdictionary transportProperties
00020 (
00021 IOobject
00022 (
00023 "transportProperties",
00024 runTime.constant(),
00025 mesh,
00026 IOobject::MUST_READ,
00027 IOobject::NO_WRITE
00028 )
00029 );
00030
00031
00032 Info<< "Reading diffusivity DT\n" << name="l00033">00033
00034 dimensionedScalar DT
00035 (
00036 transportProperties.lookup("DT")
00037 );
Para resolver o problema, a temperatura deve ser definida em um campo geométrico (template geometricField). Fisicamente, a temperatura é uma variável escalar e, portanto, deve-se criar um campo geométrico de um escalar. O template volScalarField constrói um campo escalar da variável T, definida no centro das células da malha mesh. A criação deste campo depende dos templates IOobject e fvMesh como argumentos de entrada. O primeiro constrói o objeto, definindo o nome da variável e do arquivo que contém os valores do campo inicial de T (O nome do arquivo que contém o campo inicial da variável é idêntico ao nome da própria variável), o registro das informações contidas em runTime (por exemplo, o diretório e a periodicidade da saída de resultados em arquivo) e as opções de entrada e saída de dados (MUST_READ indica que o campo inicial de T deve ser necessariamente lido e AUTO_WRITE configura a saída automática de T em arquivo ao longo do tempo). O segundo parâmetro para a criação de T é o template fvMesh definido como mesh, que indica onde o campo escalar T será alocado e inserido.
A equação da difusão possui apenas uma propriedade física de transporte, a condutividade térmica DT (representada pela letra D na equação da difusão). A leitura das propriedades de transporte é realizada através do template IOdictionary, construída a partir do template IOobject como argumento. O template IOdictionary, por sua vez, é derivada de outros dois templates, dictionary e IOobject, proporcionando funcionalidade na entrada e saída de dados automática a partir de um banco de dados. O template dictionary define uma lista de palavras chave, onde cada uma destas é associada a um número arbitrário de valores. A construção de transportProperties declara o nome do arquivo que contém as propriedades de transporte e seu diretório (runTime.constant()) e as regras de leitura (MUST_READ) e saída (NO_WRITE) de dados em arquivo. Por fim, a criação da variável escalar dimensional DT é construída pelo comando transportProperties.lookup("DT"), que procura no arquivo transportProperties a palavra-chave "DT" e associa um valor dimensional a essa variável.
As
etapas descritas acima apenas inicializam os dados para
runTime, a malha
mesh, o campo inicial de
T e a propriedade de transporte
DT do código principal. O próximo passo é programar o algoritmo de solução do problema específico.
O laço
for, iniciado na linha 51 do código principal, tem o intuito de repetir as instruções no interior do laço para cada passo de tempo (incrementado por
runTime++). O laço é realizado até que a condição de
runTime.end() seja satisfeita (retorne
true).
O primeiro comando no interior do laço é a declaração do
header readSIMPLEControls.H para leitura dos parâmetros do método SIMPLE de acoplamento pressão-velocidade e das condições de ortogonalidade da malha. Estes parâmetros são lidos no arquivo
fvSolution do caso analisado. Apesar de não ser necessário aplicar o método SIMPLE para resolver a equação da difusão pura, esta biblioteca lê o número de iterações para correção dos fluxos devido à não-ortogonalidade da malha. Isso é necessário pois o OpenFOAM divide o cálculo do fluxo através das faces em duas parcelas chamadas contribuição ortogonal e a correção não ortogonal. Esta correção é realizada um laço para ajustar o fluxo das propriedades nas faces dos volumes, semelhante ao esquema
defferred correction. Maiores detalhes podem ser encontrados na
tese do Prof. Jasak.
A discretização por volumes finitos é realizada pelo template
fvm, armazenando as equações discretizadas em sua forma matricial com o template
fvMatrix e montando um sistema de equações lineares resolvido pelo comando
solve. Este último comando retorna na tela para o usuário dados estatísticos da solução, como a convergência do sistema, número de iterações, etc. A definição das formulações de discretização usadas na simulação estão colocadas no arquivo
fvShemes. A solução do sistema algébrico, cujo método de solução está selecionado no arquivo
fvSolution, retorna a temperatura em cada célula da malha. A biblioteca
write.H, colocada abaixo, escreve os arquivos de resultados quando o comando de classe
runTime.output() for válido.
00001 if (runTime.outputTime())
00002 {
00003 volVectorField gradT = fvc::grad(T);
00004
00005 volScalarField gradTx
00006 (
00007 IOobject
00008 (
00009 "gradTx",
00010 runTime.timeName(),
00011 mesh,
00012 IOobject::NO_READ,
00013 IOobject::AUTO_WRITE
00014 ),
00015 gradT.component(vector::X)
00016 );
00017
00018 volScalarField gradTy
00019 (
00020 IOobject
00021 (
00022 "gradTy",
00023 runTime.timeName(),
00024 mesh,
00025 IOobject::NO_READ,
00026 IOobject::AUTO_WRITE
00027 ),
00028 gradT.component(vector::Y)
00029 );
00030
00031 volScalarField gradTz
00032 (
00033 IOobject
00034 (
00035 "gradTz",
00036 runTime.timeName(),
00037 mesh,
00038 IOobject::NO_READ,
00039 IOobject::AUTO_WRITE
00040 ),
00041 gradT.component(vector::Z)
00042 );
00043
00044
00045 runTime.write();
00046 }
Com o intuito de também escrever em arquivo os componentes do gradiente da temperatura, calcula-se
gradT pela classe
fvc que realiza operações tensoriais explícitas com os dados da malha. Note que, como a temperatura é uma variável escalar, seu gradiente será uma entidade vetorial, e portanto o tipo da variável
gradT deve ser um
volVectorField. As linhas 5, 18 e 31 do código acima decompõem
gradT nas direções X, Y e Z do domínio e definem esta variável como
AUTO_WRITE para que a saída em arquivo seja automática. Cada decomposição é, portanto, um
volScalarField e será formado pela componente do vetor de
gradT da respectiva direção. O comando
runTime.write() escreve em arquivos os valores de temperatura nos centros das células, seus gradientes e os componentes do gradiente.
O algoritmo repete-se até o fim do laço no tempo e, por consequência, o final da simulação. Na minha opinião, o
solver laplacianFoam possui o código mais simples de todos. A complexidade dos códigos aumenta junto com o detalhamento dos fenômenos físicos considerados, sendo necessário o desenvolvimento de algoritmos mais elaborados, o uso de mais templates, classes e comandos.
Espero que este post seja de ajuda aos que estão começando a estudar e desbravar a programação no OpenFOAM. Se quiserem discutir mais sobre o assunto ou mesmo trocar idéias, deixe um comentário aqui.
Um abraço!
