πŸ“œ ⬆️ ⬇️

OpenFOAM in terms of programmer physics


Intel has developed quite a lot of software for modeling various physical processes. In some of them we use the OpenFOAM package, and in this post I will try to give a brief description of its capabilities.
What is OpenFOAM? This, using the term Wikipedia, is an open (GPL) platform for numerical modeling, primarily for modeling related to solving partial differential equations by the finite-volume method, and first of all for solving problems of continuum mechanics.
EFC: the evolution of two immiscible liquids of different densities, initially separated by a thin partition (for example, β€œlockExchange” from the standard OpenFOAM delivery). Transition colors denote grid cells where there is a proportion of both liquids (more precisely: the liquid volume method is used in the simulation).

In fact, OpenFOAM, in fact, consists of two parts: firstly, it is a class library for many operations necessary for numerical simulation, and secondly, a library of programs (β€œsolvers”) using these classes and allowing to solve specific modeling problems . Accordingly, my post will consist of two sections: in the first, I will (briefly) describe the program library and talk about using OpenFOAM from the perspective of a person who does not want to go into the program code directly; in the second, I'll try a little more detail on programming under OpenFOAM.

Using OpenFOAM


I will not elaborate on step-by-step instructions on how to prepare data for OpenFOAM and directly on its launch. The simplest example is described in the official manual and gives a pretty good idea of ​​how to work with almost any OpenFOAM-solver. I can only say that to start OpenFOAM you need to prepare a whole directory with a certain file structure in which to place various parameters of the task (description of the grid, simulation duration, parameters for choosing a time step, etc.), physical properties of the system under study (viscosity / heat capacities of liquids or the like) and initial / boundary conditions. The specific structure of directories and files, of course, is determined by a specific solver program, but the general scheme is almost always the same. In particular, all files containing different data are recorded in a special hierarchical format of β€œdictionaries”, which have a fairly simple structure; examples are on the link above.
')
After preparing the file structure, you need to run a solver. All standard OpenFOAM solvers are console applications with minimal user interaction. In the end, if all goes well, several new subdirectories will be created in the working directory, in which the results of the calculations will be written at the indicated intermediate points in time. These files can be converted to various formats (for example, to the popular VTK), for this, there are many converters in OpenFOAM, and you can directly view, for example, in ParaView.

I will also list some of several standard OpenFOAM-solvers. OpenFOAM includes programs for:
β€’ solutions of β€œsimple” problems of motion of compressible and incompressible liquids with the support of different thermodynamic models and different models of turbulence;
β€’ solving problems of the movement of liquids / gases, taking into account chemical processes;
β€’ solving problems of magnetic hydrodynamics;
β€’ and much more, even to the Black – Scholes financial model.

A full list of standard solvers with a brief description is available on the OpenFOAM website . However, unfortunately, many of them are documented rather badly, and in order to understand which particular equations they solve and which parameters they require, one has to go into the code. In general, this applies to many aspects of OpenFOAM β€” many features are either not documented at all, or are documented quite poorly, and you can learn about them either from code or from examples.

Programming under OpenFOAM


Base classes

What does OpenFOAM look like from a programmer's point of view? Well, first, of course, there is a set of many elementary classes and support for elementary operations that can be expected from any physical code. Vectors, tensors; scalar, vector, tensor product; trace of a tensor, etc. True, unfortunately, everything is only three-dimensional, and with tensors above the second, if I am not mistaken, the rank is bad.

Then, of course, the mesh class (more precisely, the whole hierarchy of classes that store different mesh data is purely geometric coordinates; information necessary for working with sparse matrices; information necessary for the finite volume method; etc.). The class (or rather, the class hierarchy again) for the physical fields defined on the grid. Hierarchy of classes for boundary conditions - all boundary conditions are defined by classes with a standard interface, which allows you to select the desired condition at the stage of launching a solver, and not at all in the process of writing it.

Naturally, for all these classes, the addition / multiplication operators are redefined, etc., the addition of two fields or the calculation (pointwise) of the Bessel function from the field is done in one line.

Differentiation

Actually, already these classes make it much easier to write code that performs numerical simulation. But in fact, the pearl of OpenFOAM, an opportunity that none of its reviews can do without is the fvc and fvm . They define functions such as grad , div , curl , laplacian , etc. to calculate various spatial derivatives (respectively, gradient, divergence, rotor, Laplacian, etc.), as well as ddt and d2dt2 for time derivatives. As a result, if you have some psi field, then its gradient is calculated by simply calling fvc::grad(psi) , etc.

Actually fvc::grad , etc. - this is an β€œexternal interface” for a number of different functions that calculate the gradient according to different schemes; The choice of a specific scheme is performed at the stage of the program execution in accordance with what is written in the task parameter files. Accordingly, it is possible, without changing the code of the solver, to change the schemes used for the numerical differentiation.

So, if I need to calculate the field \ nabla \ vec U + \ alpha \ cdot \ Delta \ varphi I write in the code:
 fvc::div(U) + alpha * fvc::laplacian(phi) 

Please note that, firstly, no grid parameters appear anywhere in this code - functions like div will recognize them (due to the fact that each field stores a link to its own grid) and will correctly take into account. Secondly, likewise, there are no explicitly granular conditions here - but they will also be correctly taken into account.

OpenFOAM can work with grids of arbitrary structure (well, it’s actually a finite volume method), but even if the standard code had limitations on the grid structure, you could get rid of them by rewriting functions like grad without correcting the solver code.

fvm :: and solving equations

But, to solve equations, it is not enough just to be able to calculate expressions with spatial derivatives. Let's start with a fairly simple case: the equation (classical for OpenFOAM reviews)
{\ partial \ rho \ vec U \ over \ partial t} + \ nabla (\ vec \ varphi \ otimes \ vec U) - \ mu \ Delta \ vec U = - \ nabla p
We want to solve it by successive steps in time, calculating values \ vec U for each next time point according to the values ​​of variables at the previous time point. (Most likely, we will have more equations in the system, but with them everything is done similarly.)

Regardless of the methods of discretization of derivatives, it is possible to search for β€œnew” values, as is well known, by an β€œexplicit” and β€œimplicit” method. In the first case, we simply calculate, based on the "current" values \ vec U , all spatial derivatives, from them we obtain the derivative \ partial (\ rho \ vec U) / \ partial t - and in accordance with it we find the "new" value \ vec U . In the second case, we are looking for such "new" values. \ vec U so that when they are substituted into the equation in all places of entry \ vec U it was executed - i.e. we get (in this case, linear, see below) the system of equations for the β€œnew” value \ vec U which we decide.

In general, it is clear that both cases can be represented as an equation A \ vec U ^ {new} = \ vec B where A - some matrix, and \ vec B - some field. In the first case, the matrix A will be diagonal \ rho / \ Delta t ), and the field \ vec B will contain a member \ rho \ vec U ^ {old} / \ Delta t and fully all calculated values ​​of spatial derivatives. And in the second case, the matrix A It will also contain non-diagonal elements defined by spatial derivatives, and, accordingly, the field \ vec B will contain fewer components.

It is precisely for distinguishing these two approaches that there are fvc and fvm . Namely, the fvc namespace fvc return the corresponding field β€” i.e. just what you need for an explicit method; and the fvm functions return a special object storing a contribution from the corresponding member to the matrix. A and in the field \ vec B .

Accordingly, on OpenFOAM, the obvious way is as follows:
 solve ( fvm::ddt(rho,U) + fvc::div(phi,U) - fvc::laplacian(mu,U) == - fvc::grad(p) ); 

and implicit so:
 solve ( fvm::ddt(rho,U) + fvm::div(phi,U) - fvm::laplacian(mu,U) == - fvc::grad(p) ); 

And that's all. Very nice and easy to read. Well, it is easy to change the explicit representation of each member to an implicit one and vice versa.

(Both of these codes take only one time step, i.e. calculate new values \ vec U by current values ​​of all fields. To get a complete simulation code, it is necessary to solve this call β€” along with similar calls for other solved equations β€” to be enclosed in an external time loop.)

Three notes. First, about why differentiation functions sometimes take two arguments. On this account there are some subtleties associated with how the form of the corresponding member (for example, laplacian(mu, psi) corresponds to \ nabla (\ mu \ nabla \ psi) , and similarly for vector fields), and with the fact that the fields, in fact, are specified both in the centers of the grid cells and on the edges of the grid.

Secondly, of course, everything is so beautiful only while the equation is linear. But it is clear that non-linearity is easily taken into account in an explicit way, but if we want to take it into account somehow implicitly, then we have to invent something manually, and there is no general reasonable method here.

Third, of course, similar approaches can be applied to other equations and problems. For example, the simple equation \ Delta \ psi + k ^ 2 \ psi = f is solved by solving one call:
 solve( fvm::laplacian(psi) + k*k*fvm::Sp(1, psi) = f ); 

- moreover (since the equation is linear and non-evolutionary), you immediately get the fully sought field \ psi .

(Here, fvm::Sp is a function that introduces the desired function directly into the equation implicitly without any derivatives. If I wrote just psi , I would receive its explicit accounting, that is, simply substituting its β€œcurrent” values.)

Dimensions

Another very convenient feature of OpenFOAM is native dimension support. All fields, all numerical parameters of the task keep their dimension. (More precisely, it is possible to get to the depths of the code where the dimension disappears; but even there, as a rule, you can add a dimension check.) Therefore, in OpenFOAM you have to try very hard to add meters to kilograms.

At first glance, this seems at least not very necessary functionality - and sometimes even frankly interfering; for example, when I first started writing simple solvers on OpenFOAM for equations like the one above \ Delta \ psi + k ^ 2 \ psi = f I was very annoyed by the need to set the correct dimensions for k and f .

But as soon as expressions become a bit more complicated, dimensional control really helps, eliminating problems and errors like β€œthis E - is the total energy in the cell? Or local energy per unit mass of a substance? Or one mole? Or one molecule? And here I calculate the force from the electric field and substitute it into the equation - do I have to multiply it by density or not? ”

True, everything is not so rosy - in some places there are pieces of code where no dimension control is performed, and sometimes you have to write such code yourself (especially if you need to do such nontrivial operations on fields that you still have to write a loop on all the cells and not use overloaded with field operators). In one place, in general, in a large piece of code, the dimension "mol" is lost, i.e. all calculations are actually carried out in moles, but this does not figure anywhere. But even that which is is very convenient.

Models for calculations of continuum mechanics

In fact, everything that was discussed above is not directly related to continuum mechanics, but refers only to the finite volume method. This allows OpenFOAM to be considered not just an MSS calculation platform, but generally a finite volume calculation platform - and indeed, there are examples of using OpenFOAM for non-MSS modeling.

But nevertheless, in addition to the classes for the finite volume method, a huge part of the OpenFOAM code refers specifically to the problems of continuum mechanics. These are already briefly mentioned models of thermodynamic properties of substances, realizing, for example, a model of constant heat capacity or calculating heat capacities according to JANAF tables. These are different models of turbulence. This is a whole hierarchy of classes for chemical calculations. These are models of surface films, etc.

Moreover, all these models are usually implemented in the full traditions of object-oriented programming: all classes for different models of the same physical process have a common interface, and the code using these classes should not be aware of which particular model is used. Moreover, the choice of a specific model is determined at the execution stage - the user can select the desired model in the settings file.

These capabilities for dynamic model selection apply not only to MCC models, but also to the main code of the finite volume method. For example, as already mentioned, specific gradient calculation schemes, or specific types of boundary conditions, are also chosen at the execution stage.

Pattern chains

OpenFOAM code makes extensive use of templates. In this case, the following interesting structure is very actively used (if someone tells you what it is called, I will be grateful). There is, say, some base class. And there is a whole set of template classes that inherit their parameter:
 template <class T> class Foo: public T 

at the same time, each such class adds to its ancestor any new properties / features, and as a result, a whole chain of such templates of the type Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
 Foo<Bar<Baz > >. 

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
 Foo<Bar<Baz > >. 

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .
Foo<Bar<Baz > >.

, , , constTransport<specieThermo<hConstThermo > >. perfectGas
β€” , () ..; hConstThermo () ; specieThermo ; constTransport () .

«» , , . , constTransport sutherlandTransport , . , , . OpenFOAM.


, .

-, , , , , β€” main . , , icoFoam PISO. , PISO main , .

-, .H- , . , main :
int main(int argc, char *argv[]) { #include "setRootCase.H" ...
setRootCase.H ( ):
Foam::argList args(argc, argv); if (!args.checkRootCase()) { Foam::FatalError.exit(); }
( , OpenFOAM β€” .)

-, OpenFOAM . , , , : , , β€” , , Field. ! ( -), - , .

, -, OpenFOAM , , , , . , OpenFOAM, (, ) , (wmake) ..


, - .. β€” . , OpenFOAM , β€” .


OpenFOAM
(User Guide) Doxygen- Programmer's guide (- , ) OpenFOAM-extend β€” OpenFOAM, , , , . OpenFOAM cfd-online.com β€” OpenFOAM; , , . OpenFOAM, pdf- ; .

Source: https://habr.com/ru/post/170675/


All Articles