Finite volume method. Finite volume method Properties of discrete circuits

Some time ago I was looking for a description of operations and processes occurring in the OpenFOAM numerical modeling library. I found many abstract descriptions of the operation of the finite volume method, classical difference schemes, and various physical equations. I wanted to know in more detail - where did these values come from in such and such an output file at such and such an iteration, what expressions are behind certain parameters in the fvSchemes, fvSolution settings files?
For those who are also interested in this - this article. Those who are well acquainted with OpenFOAM or the methods implemented in it - write about the errors and inaccuracies found in a personal message.

There were already a couple of articles about OpenFOAM on Habré:

Therefore, I will not dwell on the fact that it is “an open (GPL) platform for numerical simulation, designed for simulations associated with solving partial differential equations using the finite volume method, and is widely used to solve problems in continuum mechanics.”

Today I will use a simple example to describe the operations that occur during calculations in OpenFOAM.

So, given the geometry - a cube with a side of 1 meter:

We are faced with the task of modeling the flow-propagation of a certain scalar field (temperature, amount of matter), which is given by the following transport equation (1) inside the volume of the body.

(1)
,

Where a scalar quantity, for example, expresses temperature [K] or the concentration of a certain substance, and expresses the transfer of a substance, mass flow [kg/s].

This equation is, for example, used to model heat propagation
,
where k is thermal conductivity, and is temperature [K].

The divergence operator is actually

operator .
Let me remind you that there is a nabla operator (Hamilton operator), which is written as follows:
,

Where i, j, k are unit vectors.
If we scalarly multiply the nabla operator by a vector quantity, we obtain the divergence of this vector:

“From the point of view of physics, the divergence of a vector field is an indicator of the extent to which a given point in space is a source or sink of this field”

If you multiply the nabla operator by a scalar, you get the gradient of that scalar:

A gradient shows an increase or decrease in some direction in the magnitude of a scalar.

The boundary conditions of the problem are as follows: there is an input face, an exit face, and the remaining faces are smooth walls.

Dividing the volume of a cube into finite volumes

Our grid will be very simple - we divide the cube into 5 equal cells along the Z axis.

Lots of formulas

The finite volume method provides that (1) in integral form (2) will be satisfied for each finite volume.

(2)
,

Where is the geometric center of the final volume.

Center of final volume

Let us simplify and transform the first term of expression (2) as follows:

(2.1) (HJ-3.12)*

As you can see, we assumed that the scalar quantity changes linearly inside the finite volume and the value of the quantity at some point inside the finite volume can be calculated as:

To simplify the second term of expression (2), we use the generalized Gauss-Ostrogradsky theorem: the integral of the divergence of the vector field over the volume is equal to the vector flux through the surface bounding the given volume. In human language, “the sum of all flows into/from a finite volume is equal to the sum of flows through the faces of this finite volume”:

(2.3)
,

Where is the closed surface limiting the volume,
- vector directed along the normal from the volume.

Vector S

Considering that the finite volume is limited by a set of flat faces, expression (2.3) can be transformed to the sum of integrals over the surface:

(2.4) (HJ-3.13)
,

Where expresses the value of the variable at the center of the face,
- area vector, coming out from the center of the face, directed away from the cell (locally), away from the cell with a lower index to the cell with a higher index (global).

A little more about vector S

In order not to store the same vector parameters twice, because It is obvious that for two neighboring cells the normal vector to the edge between the cells, directed towards the center of the cell, will differ only in direction-sign. Therefore, an owner-neighbor relationship was created between the edge and the cell. If the area vector (global, positive direction from a cell with a lower index to a cell with a larger index) indicates FROM the center of the cell, such a relationship between the cell and the vector, or more precisely between the cell and the edge, is denoted owner). If this vector points inside the cell in question, then the neighbor. The direction affects the sign of the value (+ for owner and - for neighbor) and this is important when summing, see below.

About difference schemes

The value at the center of the face is calculated through the values at the centers of adjacent cells - this method of expression is called a difference scheme. In OpenFOAM, the type of difference scheme is specified in the file /system/fvSchemes:

DivSchemes ( default none; div(phi,psi) Gauss linear; )

Gauss- means that the central difference scheme is selected;
linear- means that interpolation from the centers of cells to the centers of faces will occur linearly.

Let us assume that our scalar quantity changes linearly inside the finite volume from the center to the edges. Then the value approximated at the center of the face will be calculated according to the formula:

Where are the weights and are calculated as

Where are the cell volumes.
For cases of skewed cells, there are more complex formulas for calculating approximation weights.

Thus, the phi_f values at the cell edge centers are calculated based on the values at the cell centers. Gradient values grad(phi) are calculated based on phi_f values.
And this entire algorithm can be represented in the form of the following pseudocode.
1. We declare an array of gradients of finite volumes, initialize it with zeros 2. We go through all the internal faces (which are not boundary) > We calculate flux_f = phi_f*S_f. Calculate phi_f values based on phi values in cell cents > Add flux_f to the gradient of the owner element and -flux_f to the gradient of the neighbor element 3. Iterate over all the boundary faces > Calculate flux_f = phi_f*S_f > Add flux_f to the gradient of the owner element (neighbor -the boundary faces have no elements) 4. Let’s go through all the elements > Divide the resulting gradient sum by the volume of the element

Time sampling

Taking into account (2.1) and (2.4), expression (2) takes the form:

(3)

According to the finite volume method, time discretization is carried out and expression (3) is written as:

(4)

Let's integrate (4):

(4.1)

Let's divide the left and right sides into:

(5)

Data for sampling matrix

Now we can obtain a system of linear equations for each finite volume.

Below is the numbering of the grid nodes that we will use.

Node coordinates are stored in /constant/polyMesh/points

24 ((0 0 0) (1 0 0) (0 1 0) (1 1 0) (0 0 0.2) (1 0 0.2) (0 1 0.2) (1 1 0.2) (0 0 0.4) (1 0 0.4) (0 1 0.4) (1 1 0.4) (0 0 0.6) (1 0 0.6) (0 1 0.6) (1 1 0.6) (0 0 0.8) (1 0 0.8) (0 1 0.8) (1 1 0.8) (0 0 1) (1 0 1) (0 1 1) (1 1 1))

Numbering of nodes-centers of cells (50, 51 - centers of boundary faces):

Numbering of face center nodes:

Element volumes:

Interpolation coefficients needed to calculate values on cell faces. The subscript "e" denotes the "right edge of the cell". Right relative to the view, as in the figure “Numbering of nodes-centers of cells”:

Formation of the sampling matrix

For P = 0.
Expression (5) describing the behavior of the quantity

Will be transformed into a system of linear algebraic equations, each of the form:

Or, according to the indices of points on the faces

And all flows to/from a cell can be expressed as a sum

Where, for example, is the flow linearization coefficient at the center point of cell E,
- flow linearization coefficient at the center point of the face,
- nonlinear part (for example, constant).

According to the numbering of the faces, the expression will take the form:

Taking into account the boundary conditions for the element P_0, the linear algebraic equation can be represented as

...substitute the previously obtained coefficients...

The flux from inlet"a is directed into the cell and therefore has a negative sign.

Since in our control expression we also have, in addition to the diffusion term, a time term, but the final equation looks like

For P = 1.

For P = 4.

A system of linear algebraic equations (SLAE) can be represented in matrix form as

A(i,j) === 40.5 0.5 0 0 0 -0.5 40 0.5 0 0 0 -0.5 40 0.5 0 0 0 -0.5 40 0.5 0 0 0 -0.5 40.5

Psi = dimensions; internalField nonuniform List 5(0.0246875 0.000308546 3.85622e-06 4.81954e-08 5.95005e-10);

On the basis of which the values for the vector are obtained

Then the vector is substituted into the SLAE and a new iteration of the vector calculation occurs.

And so on until the discrepancy reaches the required limits.

Links

* Some equations in this article are taken from the dissertation of Jasak Hrvoje (HJ is the equation number) and if anyone wants to read more about them (

Previously, the subdomain method was mentioned, which served as the starting point for a number of numerical methods. One such method is the finite volume method. This same method is a representative of another widespread class - integral methods. From the classical form of notation of the subdomain method, the division of the computational domain into subdomains and the integration of the residual over the subdomain are taken. The difference is the absence of an explicit recording of the approximating (test) function. But, as before, we are trying to “exactly” solve the equation in each subdomain. Therefore, the original equation is integrated over the subdomain. Integral methods are characterized by the fact that first the integral of the differential equation is taken, and an integral form of writing the equation is obtained. The equation in this form is then applied to individual grid cells. In this case, cells and subareas are one and the same.

In fact, the integral form of writing equations has (from the point of view of physics) an even wider range of application than the differential one. The fact is that in the presence of function discontinuities, differential equations are not applicable, and their integral analogues continue to work, work and work…. Unfortunately, when they are implemented numerically, this advantage is sometimes lost.

As a rule, integrals from equations have a simple and understandable physical meaning. For example, consider the continuity equation. The original differential equation is written

Let's integrate it over the volume V, which has a surface S, and over time in the interval from t 0 to t 1. When integrating derivatives, we use the Stokes formula (its special cases are called the Green and Ostrogradsky-Gauss formulas). As a result we get

In this notation, the difference between the first two integrals means the change in mass in a given volume over the time interval under consideration. And the double integral shows the mass flowing into a given volume through the surface bounding it over the same period of time. Naturally, since we are talking about numerical methods, these integrals are calculated approximately. And here the questions of approximation begin, similar to those considered in the finite difference method.

Let's consider one of the simplest cases - a two-dimensional rectangular uniform grid. In the finite volume method, the values of functions are usually determined not at grid nodes, but at the centers of cells. Accordingly, it is also not the grid lines in each direction that are indexed, but the layers of cells (see figure).

j-1

j+1

k-1

k+1

For this case, the integral form of the equation will be written as follows

As you can see, in this case we received an ordinary equation, which we could also write using the finite difference method. This means that the same methods of studying stability can be applied to it. (A quick question: is this scheme stable?)

But if we got the same thing, then was it worth building this whole garden? In the simplest cases, we really don’t get any benefits. But in more complex situations, the benefits emerge. First, as noted above, such methods (even in such a simple implementation) describe discontinuities and areas with high gradients much better. At the same time, the fulfillment of the laws of conservation of mass, momentum and energy is guaranteed, since they are observed in each cell. Secondly, these methods can withstand a wide variety of abuses on the grid. Even curvilinear, uneven and irregular grids do not throw these methods off track. These benefits are especially often felt when boundary conditions are specified.

j-1

j+1

k-1

k+1

For example, for the case shown in the figure, the integral form of the equation will have the form

that is, simply where we took the integral over the area of the full cell, now we take it over the “trimmed” area, where we took the integral over the full edge, now we take it over the remaining part of it. An integral over the boundary section was added. But it is easily found from the boundary conditions. In particular, if no mass flow is supplied through the wall (and also no mass is carried away from the surface and/or we neglect the mass flow of ions losing charge on the wall), then such an integral is simply equal to zero. In a similar form of energy equation, the flow through the wall, as a rule, has to be taken into account. But it is also not difficult to find from the boundary conditions (if they are set correctly).

To reinforce this, let us describe what the application of the finite volume method to one of the momentum conservation equations will look like. Let us take the flat stationary case for singly charged ions. We neglect viscosity and elastic collisions. We get the equation

For a rectangular mesh (see figure above) we get

The simplest approximation of such an equation can be written as follows:

after reductions we get the formula

algorithm modeling program

The starting point of the finite volume method (FVM) is the integral formulation of the laws of conservation of mass, momentum, energy, etc. Balance relationships are written for a small control volume; their discrete analogue is obtained by summing over all faces of the selected volume of flows of mass, momentum, etc., calculated using some quadrature formulas. Since the integral formulation of conservation laws does not impose restrictions on the shape of the control volume, MCM is suitable for discretizing fluid dynamics equations on both structured and unstructured grids with different cell shapes, which, in principle, completely solves the problem of the complex geometry of the computational domain.

It should be noted, however, that the use of unstructured meshes is quite complex in algorithmic terms, labor-intensive to implement and resource-intensive to carry out calculations, especially when solving three-dimensional problems. This is due both to the variety of possible shapes of the cells of the computational grid, and to the need to use more complex methods for solving a system of algebraic equations that does not have a specific structure. The practice of recent years shows that the advanced development of computing tools based on the use of unstructured grids is only possible for fairly large companies with the appropriate human and financial resources. It is much more economical to use block-structured grids, which involves dividing the flow region into several subregions (blocks) of a relatively simple form, in each of which its own computational grid is constructed. In general, such a composite mesh is not structured, but within each block the usual index numbering of nodes is retained, which allows the use of efficient algorithms developed for structured meshes. In fact, to move from a single-block grid to a multi-block one, you only need to organize the joining of blocks, i.e. exchange of data between adjacent subareas to take into account their mutual influence. Note also that dividing a task into separate relatively independent blocks naturally fits into the concept of parallel computing on cluster systems with processing of individual blocks on different processors (computers). All this makes the use of block-structured meshes in combination with MCM a relatively simple but extremely effective means of expanding the geometry of the problems being solved, which is extremely important for small university groups developing their own programs in the field of fluid dynamics.

The above-mentioned advantages of MKO served as the basis for the fact that in the early 1990s. It is this approach, focused on the use of block-structured grids, that was chosen by the authors as the basis for developing their own wide-profile software package for problems of fluid dynamics and convective heat transfer.

Description

Informal

A certain closed region of liquid or gas flow is selected, for which a search is made for fields of macroscopic quantities (for example, velocity, pressure) that describe the state of the medium in time and satisfy certain laws formulated mathematically. The most commonly used are conservation laws in Euler variables.

For any value, at every point in space, surrounded by some closed finite volume, at the moment of time there is the following relationship: the total amount of a quantity in the volume can change due to the following factors:

In other words, when formulating the MKO, the physical interpretation of the quantity being studied is used. For example, when solving problems of heat transfer, the law of conservation of heat in each control volume is used.

Mathematical

Modifications

Literature

Patankar S.V. Numerical solution of problems of thermal conduction and convective heat transfer during flow in channels = Computation of conduction and Duct Flow Heat Transfer: Transl. from English - M.: MPEI Publishing House, 2003. - 312 p.

5. Properties of discrete circuits

5.1 Accuracy

Accuracy characterizes the acceptability of the numerical scheme for its practical use. Assessing the accuracy of a discrete circuit seems to be a very difficult task, since it turns out to be almost impossible to separate errors that arise as a result of the properties of the circuit from errors that arise as a result of other factors (such as rounding errors, inaccuracy in specifying boundary and initial conditions, etc.).

When talking about the accuracy of a discrete scheme, they usually mean the error in approximation of derivatives 27 . In particular, if the approximation error is comparable to the second power of the computational grid step, then the discrete scheme is said to have second order accuracy. This issue was discussed in more detail in § 3.

5.2 Consistency

The discrete circuit is called agreed upon with the original differential equation, if, when the computational mesh is refined, the approximation error (see § 3) tends to zero,

There are known calculation schemes in which additional conditions must be met to achieve consistency [Anderson and K]. Since checking the consistency of calculation schemes is the task of software developers (and not users) of the software, this issue will not be discussed in more detail here.