Remarks on generation of the orthogonal structured grids Текст научной статьи по специальности «Математика»

UDC 519.632.4

DOI: 10.18698/1812-3368-2019-1-16-26

REMARKS ON GENERATION OF THE ORTHOGONAL STRUCTURED GRIDS

S.I. Martynenko [email protected]

[email protected]

Central Institute of Aviation Motors, Moscow, Russian Federation

Institute of Problems of Chemical Physics, Russian Academy of Sciences, Chernogolovka, Moscow Region, Russian Federation

Abstract

Grid generation techniques have contributed signifi-

cantly toward the application of mathematical modeling in large-scale engineering problems. The structured grids have the advantage that very robust and parallel computational algorithms have been proposed for solving (initial-)boundary value problems. Orthogonal grids make it possible to simplify an approximation of the differential equations and to increase computation accuracy. Opportunity of the orthogonal structured grid generation for solving two- and three-dimensional (initial-)boundary value problems is analyzed in the article in assumption that isolines or isosurfaces of d (=2,3) functions form this grid. Condition of the isolines/isosurfaces orthogonality is used for formulation of the boundary value problems, the solutions of which will be form the orthogonal grid. A differential substitution is proposed to formulate the boundary value problems directly from the orthogonality condition of the grid. The substitution leads to the general partial differrential equations with undetermined coefficients. In the two-dimensional case, it is shown that the orthogonal grid generation is equivalent to the solution of partial differential equations of either elliptic or hyperbolic type. In three-dimensional domains, an orthogonal grid can be generated only in special cases. The obtained results are useful for mathematical modeling of the complex physicochemical processes in the technical devices

Keywords

Orthogonal grids, mathematical modeling, boundary value problems

Received 04.04.2018 © Author(s), 2019

The work was supported by the RNF under the agreement no. 15-11-30012 of 07/08/2015 on the topic: "Supercomputer simulation of physical and chemical processes in the high-speed direct-flow propulsion jet engine of the hypersonic aircraft on solid fuels''

Introduction. The rapid development of computers and mathematical modeling methods had a strong influence on mechanical engineering: a significant amount of experimental research in the new technology development can be replaced by mathematical modeling results. This significantly reduces the time and cost of research and development work. However, in the early 1990s, it became clear that the greatest effect from the use of mathematical modeling methods in mechanical engineering was achieved using specialized software built on the ''black box'' principle. The user of such software is busy only with setting the task (i.e., setting the geometry region, initial and boundary conditions, choosing materials and equations describing processes, etc.) analyzing the results, and computational experiment details are not available to him. Software packages, arranged according to the ''black box'' principle, are a powerful tool in the hands of an engineer that will allow him to concentrate on optimizing work processes occurring in developed technical devices, rather than wasting time writing and debugging computer programs.

All the problems associated with the creation of modern mathematical modeling software can be divided into three groups: 1) ''physical''; 2) "mathematical"; 3) "computer". ''Physical'' problems are associated with the mathematical description difficulty of complex physicochemical processes, such as hydrodynamics and heat and mass transfer processes in multiphase reactive media, turbulent transfer, etc. ''Computer'' problems arise due to compatibility difficulties of various and quickly updated software and hardware, and ''mathematical'' problems are associated with the formalizing complexity of computational experiment main stages: building a computational grid, approximating the fundamental (non-)linear (integro-)differential equations and the effective solution of ill-conditioned systems of high order linear algebraic equations on a sequential or parallel computer. The complexity of ''physical'' problems is due to the diversity of the processes being modeled and their mathematical description depth, the ''mathematical'' problems are the result of insufficient knowledge of the underlying (non-)linear (integro-)differential equations and insufficiently developed methods for the numerical solution of mathematical physics equations.

Currently, several software packages, such as STAR-CD, CFX, FLUENT, PHOENICS, SIGMAFLOW and mesh generators, such as ICEM CFD, Gambit, NETGEN, etc., are developed and widely used for engineering calculations. The listed programs are still very far from perfect, although they can significantly reduce the mathematical modeling disadvantages.

Computational grid characteristics strongly influence the computation time and the obtained solution accuracy. On the one hand, in a region with a comp-

lex geometry, it is easier to construct an unstructured grid, but it is more difficult to solve grid equations (in modern packages algebraic multigrid methods are most often used [1], which do not use information about the grid). On the other hand, it is more difficult to construct a regular grid in the same area, but it is easier to solve grid equations (for example, to apply geometric multigrid methods [1] that use grid information). In modern software packages, it is often the grid type that determines the choice of an iterative method for solving (linearized) grid equations obtained as a result of approximation of the original (non)linear (integro-)differential equations of a mathematical model.

Previously it was assumed that the Seidel method is completely unsuitable for numerical solution of boundary value problems, since it requires performing O(N2) arithmetic operations in the two-dimensional case, where N is the number of unknowns vector components of the resulting linear algebraic equations system (for a single equation, N is the number of grid nodes) [2]. The use of upper relaxation method makes it possible to reduce arithmetic operations cost to O(N3/2) arithmetic operations [2], however, the problem of finding the relaxation parameter optimal value immediately arises [3]. The complexity of classical multigrid methods is even lower. These methods cannot be represented as a unified algorithm, but with optimal adaptation of their problem-dependent components, have a minimal (unimprovable) complexity O(N)of arithmetic operations. In other words, the lower the computational algorithm complexity, the lower the computation formalization level and the more difficult it is to use such an algorithm in software built on the ''black box'' principle. It is obvious that the fundamental idea of the first multigrid algorithm author, R.P. Fedorenko, described in [4], has a powerful potential, but the classical multigrid methods in their present form are difficult to use in ''black box'' software, due to problems associated with optimizing their problem-dependent components. In addition, the classical multigrid methods do not allow to effectively parallelizing smoothing iterations, especially on coarse grids [1].

In the early 1990s, it was done the attempt to realize a different form of R.P. Fedorenko fundamental idea, for this purpose, the basic multigrid principle was used in the single-mesh algorithm [5-7]. The developed iterative method for the numerical solution of boundary value problems is called robust multigrid technique (RMT). It was shown that if the computational grid is structured (that is, regular, but generating a special sequence of subgrids — a multigrid structure), the complexity of the Seidel method can be reduced to O(N lg N) arithmetic operations without involving problem-dependent components. In addition, it was shown in [5, 6, 8], that the maximum parallelization efficiency

(as compared to an unimproved sequential algorithm, which has the complexity of O(N) arithmetic operations) will be independent of the number of processors used. Thus, RMT has a minimum of problem-dependent components, such as the Seidel single-grid method (criterion for stopping iterations, the lower relaxation parameter for nonlinear problems and the ordering of unknowns for aniso-tropic problems), and close to the optimal O(N lg N) arithmetic operations comparable to the classical multigrid methods. The most serious limitation for RMT is the structured computational grid requirement.

Currently, there are many ways to build grids, but all known grid generation algorithms have a low level of formalization and parallelism, and, as a rule, require expert evaluation of the built mesh quality [9, 10]. Often the grid should be changed in the process of solving mathematical model equations, adapting to the features of the solution and (or) changing the region. The development of a robust algorithm for constructing adaptive grids in multiply connected areas is one of the key problems of the mathematical apparatus for modern software tools designed to solve engineering problems.

In terms of practical application, the key problem is the development of a robust structured grids generator in areas with complex geometry, which can later be used to numerically solve boundary and initial-boundary problems with the help of RMT. Of particular theoretical and practical interest are orthogonal grids, which make it possible to significantly simplify the approximation of the original (non)linear (integro-)differential equations of a mathematical model.

The purpose of this work is to study the possibility of constructing orthogonal computational grids in two- and three-dimensional areas.

Governing equations for orthogonal computational grids. Consider the three-dimensional case and assume that the computational grid is formed by isosurfaces U(x, y, z) = const, V(x, y, z) = const and W(x, y, z) = const, where U, V and W are some functions. Let (x0, y0, z0) be some intersection point of three isosurfaces. We construct tangent planes to isosurfaces U (x, y, z) = const, V(x, y, z) = const and W(x, y, z) = const, at a point (x0, y0, z0) and require that these planes be mutually perpendicular. Due to the fact that the point (x0, y0, z0) is chosen arbitrarily, the condition of mutual perpendicularity of the tangent planes will take the form

U'x Vx + U'y V'y + U'z V'z = 0; Vx W'x+ V'y W'y + V'z W'z = 0;

w'x u'x+w; u ;+w'z u z = 0.

(1a) (1b) (1c)

A computational grid, formed by isosurfaces U(x, y, z) = const, V(x, y, z) = = const and W(x, y, z) = const, is called orthogonal if the functions U, V and W satisfy the orthogonality condition (1a)-(1c). The task is to build such a system of differential equations, the solution of which are functions U, V and W satisfying the orthogonality condition (1a)-(1c). Then the generation of an orthogonal grid in a certain region will be reduced to the solution of this system and the construction of isosurfaces U (x, y, z ) = const, V ( x, y, z ) = const and W (x, y, z ) = const.

We represent functions U, V and W in the form

U = an^'x + a^'y + a^'z + bn®'x + bu®'y +

+ b1sOz + cnAx + C12Ay + C13Az ; (2a)

V = a21^'x + a22^'y + a23^'z + b21®'x + b22®'y +

+ b23&z + C21^x + C22^y + 023^'z ; (2b)

W = a31x¥'x + a32^'y + a33x?'z + b31®'x + b32®'y +

+ b33®'z + C31K + C32^'y + C33K, (2c)

where aj, bj, Cj, i, j = 1, 2, 3, are some coefficients; Y, O and A are some functions. For convenience of transformations, we introduce the notation

x = Ç1, y = Ç2, z = Ç3,

U = 01, V = 02, W = 03. Then equations (2a)-(2c) take the form

®i = X (ain^ + binO'qn + CinA'qn ), i = 1,2,3, (3)

n = 1

i.e., values i = 1, 2, 3 correspond equations (2a), (2b) and (2c). Differentiation of(3)leads to

(®i)'çt = i (an+ hn^'L+ CinKn?k ), uk = 1,2,3

n=1

(®j )'çt =i (ajtfh çk + bji®h çk + Cji^h ^ j, k = 1,2,3

Taking into account the accepted notation, the condition of the grid orthogonality (1) takes the form

+ ainCji ^ n qk A^ ; çk + ajibin^l „ ^ Yç i + knb^ n ^ O^, ^ + + bincjl®ln Çk Açl Çk + ajlCinAçn Çk Çk +

+ bjlcin^qn Çk ^Çl Çk + Cinc l Açn Çk Açl Çk ) , (4)

where i = 1, j = 2 correspond to (la), i = 2, j = 3 correspond to (lb) and i = 3, j = 1 correspond to (1c). In view of the identity

3 3

X ainbjl^"n Çk Çk + X ajlbin®l Çk ^Çl Çk = n, l = 1 n, l = 1

3

= E (<ainbjl + ajnbil )^ç'n Çk Çk n, l = 1

individual terms in (4) can be combined, then

3

0 = S ( Çk ^ Çk +Ы)1 + üjnhl Çk Ф"

Çl Çk +

k,n,l = 1

+ (ainCjl + ajnCil Çk Açl Çk + binbjl®qn Çk °çl Çk +

+ (binCjl + bjncil )Ф"n çk Açl Çk +CinC jj Aç n Çk Aç l Çk ). (5)

Now we transform the "diagonal" terms in (5):

3 3 2

X ainajl^çn Çk ^Çl qk = X aimajm( ^çm qk ) + k ,n,l = 1 k ,m=1

3 2 3

+ 1 I X (ainajl + auaj n ) n Çk ^Çl Çk ;

(6a)

k = 1 n = 1 l = n +1 3 3

X binbjl®ln Çkl Çk = X bimbjm (Oçm Çk )2 + k,n,l = 1 k,m = 1

+ X X X (binbji + bilbjn ÇkФ& Çk ; (6b)

k = 1 n = 1 l = n +1

3 „ „ 3 ( » \2

X cincjl Aç„ Çk Açl Çk = X cimcjm (Açm Çk ) + k, n, l=1 k,m=1

3 2 3

cincjl + cilcjn ) Açn Çk Açl Çk . (6c)

k=1 n=11=n+1

We note that the first terms in the right-hand side of (6a)-(6c) can be converted as follows:

3 3

X aimajm m ) "X aikajk (^qк Чк ) + m = 1 к = 1

2 3 2

X X (ainajn + ailajl)( ^i « ); (7a)

n = 1 l = n + 1

ч2 3

X bimbjm (Фдm ?k ) = X bikbjk (Фдк ?k )

к, m = 1 к = 1

+ X X (binbjn + bilbjl) ( 1 st ); (7b)

n = 1 l = n +1

3 2 3 2

X cimcjm (A?m qk ) = X cikcjk (A?k qk ) + k ,m = 1 k=1

2 3 ,2

+ X X (CinCjn + CilCjl) (A^ I qk ) . (7c)

n=1 l=n+1

Substituting (7a)-(7c) into (6a)-(6c), and (6a)-(6c) into (5), we obtain the final equation

3 2 2 3 2

0 =X aikajk (Yqk?k) +X X (ainajn + auajl)(Vg^) +

k = 1 n = 11 = n +1

3 2 3

+ X X X {ainajl + ailajn ?k ^ll ?k +

k = 1 n = 11 = n +1

+ X bkbjk K ?k )2 + X X (binbjn + bilbjl) (o^l ,n )2 + k = 1 n = 11 = n +1

3 2 3 3

+ X X X (binbjl + bilbjn ) O"n gt l ?k + X CikCjk (Ar,r, f + k = 1 n = 11 = n +1 k = 1

XXX (binbjl + bilbjn ) Ф" n ?к Ф? l Як + X С1кС,к к ?к ) к = 1 n = 11 = n +1 к = 1

3 2 3 2 3

+ Х X (cinc)n + wjl)(л"l?n) + XX X ( cincjl + cilcjn ) n дк A(;l дк + n = 11 = n +1 к = 1 n = 11 = n +1

3

X ((ainbjl + ajnbil n ?к Ф"l ?к + (aincjl + ajncil n ?к l ?к + к ,n,l = 1

(bincjl + bjncil )Ф"п « l Ск ) • (8)

Recall that expression (8) contains 27 indefinite coefficients by choosing which properly this equation can be significantly simplified. Next, we consider the two- and three-dimensional cases.

2D orthogonal grids. The two-dimensional case matches the set i = 1, j = 2, O = A = 0 and (-)Z = 0. Then (8) is reduced to

2 2 2 aiia2i (Yxx ) + ai2a22 y ) + (ana2i + ai2a22) (^y ) +

+ (aiia22 + ai2a2i)^Xx ^"xy + (aiia22 + a^i^Xy^y = 0.

We assume aii = ai2 = a2i = - a22 = i, we arrive at the following statement. Statement 1. Let the function Y(x, y) satisfies equality

№ )2 =(Vyy )2. (9)

Then the isolines of thefunctionsU and V, defined as

U = V'x ; (i0)

V = x¥'x-W'y, (ii)

form an orthogonal grid.

As 0 = (x¥"xx )2 -(W'yy )2 = (WZx + ^"yy ) (Y"xx - ^"yy ), twospecialcasesare possible.

Statement 2. Let the function x, y) satisfy the Laplace equation

vjcx +v;y = 0. (i2)

Then the isolines of the functions U and V defined as (10) and (11) form an orthogonal grid, and these functions satisfy the Laplace equation and the Cauchy — Riemann conditions.

< From (i0) and (ii) follows

U + V = 2 W'x; (i3)

U - V = 2Y'y. (i4)

Differentiating (i3) by x and differentiating (i4) byy then summing them up considering (i2), we get

U'x+ Vx + U'y - V'y= 2 (Wxx +V"yy ) = 0. (i5)

Likewise differentiating (i3) by y and differentiating (i4) by x and taking into account the theorem on independence of a mixed derivative of the differentiation order, we obtain

U y + V'= U'x- vx.

(16)

Summation and subtraction (i5) and (i6) leads to the Cauchy — Riemann conditions U'x - Vy = 0 and U'y + V'x = 0, which lead to the Laplace equations U"xx + Uyy = 0 and V£x + V"y = 0. In this case, a general view of the functions U and V is U(x, y) = fi(x + iy) + f2(x - iy) and V(x,y) = gi(x + iy) + + g2(x - iy). ►

Statement 3. Let a function Y(x, y) satisfy an equation of hyperbolic type

Then the isolines of the functions U and V, defined as (10) and (11) form an orthogonal grid, and the functions satisfy an equation of hyperbolic type.

< Analogically the differentiating (i3) and (i4) based on (i7) leads to Ux - Uy = 0 and VX + Vy = 0, which implies that U"xx - U"y = 0 and V'Xx - Vyy = 0. In this case, a general view of the functions U and V is U(x, y) = = fi(x + y) + f2(x - y) and V (x, y) = gi(x + y) + g2(x - y). ►

In the two-dimensional case, orthogonal grids are most often constructed using the Laplace equations, and the grid is thickened by the appropriate construction of contour lines.

3D orthogonal grids. Let us return to (8) and analyze the possibility of orthogonal grids constructing in the three-dimensional case. At first, we require that aikajk ^ 0 at ainajn + auajl = 0. We have previously considered two-dimensional grids as a special case of three-dimensional grids. Now three-dimensional grids will be considered as a generalization of two-dimensional grids, i.e., we assume that some of the coefficients are already known from the analysis of the two-dimensional case, in particular aii = ai2 = a2i = -a22 = i. Taking into account the symmetry of the transformations, expressions (2a)-(2c) take the form

U = V'x +% + a^Z + biiO'x + bi20y + bisOZ + Ax + ci2Ay + AZ;

V = n + a23x¥'z + b2iOx + Oy + OZ + c^Ax + c^My + c2sAZ;

W = a3i^'x + asi^'y + ass^Z + bsi^'x +®'y + c32A'y - AZ.

Condition ainajn + auaji = 0 for i = i, j = 2 (i.e., equation (ia)) takes the form

(17)

аца21 + ai2a22 = 0, n = 1, l = 2; ^11^21 + ^13^23 = 0, n = 1, l = 3; Ü12Ü22 + 013^23 = 0, n = 2, l = 3.

If an = «12 = «21 = ~a22 = 1, then the second and third equations are incompatible. It follows that the grid will be orthogonal only when Wz = 0. Further, it follows from (1b) that a31 = a32 = 0, i.e., the equations describing a three-dimensional orthogonal grid are

U (x, y) =n +%; y (x, y) = n -^y; W (z) = 9(z).

Thus, in the three-dimensional case, an orthogonal grid can be constructed only in special cases:

1) if a three-dimensional region is formed by the movement of a certain surface along a guideline 9 (z) with preservation of the condition Wz = 0. We can give the following example: orthogonal grids in the simplest bodies — parallelepiped, cylinder, etc.;

2) if the three-dimensional region is axisymmetric, that is, it is formed by rotating the two-dimensional region.

Conclusion. It is shown that there are two methods for constructing orthogonal grids (based on partial differential equations of elliptic and hyperbolic type or in an equivalent formulation related to minimizing the functional) in two-dimensional domains and two particular cases (displacement and rotation of an orthogonal two-dimensional mesh) in three-dimensional domains. Equation (8) can also be used to construct structured grids that are close to orthogonal.

Acknowledgments

The author is sincerely grateful to M.P. Galanin, Dr. Sc. (Phys.-Math.), Professor, for the interest in the work and critical discussion of the results obtained.

Translated by V. Shumaev

REFERENCES

[1] Trottenberg U., Oosterlee C.W., Schüller A. Multigrid. Academic Press, 2001.

[2] Marchuk G.I. Metody vychislitel'noy matematiki [Computational mathematics methods]. Moscow, Nauka Publ., 1989.

[3] Hageman L.A., Young D.M. Applied iterative methods. Academic Press, 1981.

[4] Fedorenko R.P. A relaxation method for solving elliptic difference equations. USSR Comput. Math. & Math. Phys., 1962, vol. 1, iss. 4, pp. 1092-1096.

DOI: 10.1016/0041-5553(62)90031-9

[5] Martynenko S.I. Mnogosetochnaya tekhnologiya: teoriya i prilozheniya [Multigrid technology: theory and applications]. Moscow, Fizmatlit Publ., 2015.

[6] Martynenko S.I. Robust multigrid technique for Black-Box software. DeGruyter, Berlin, 2017.

[7] Martynenko S.I. Robust multigrid technique for solving partial differential equations on structured grids. Vychislitel'nye metody i programmirovanie [Numerical Methods and Programming], 2000, vol. 1, pp. 82-103 (in Russ.).

[8] Martynenko S.I., Volokhov V.M., Yanovskiy L.S. Parallel geometric multigrid. Int. J. Comp. Sci. Math., 2016, vol. 7, no. 4, pp. 293-300. DOI: 10.1504/IJCSM.2016.078741

[9] Thompson J.F., Soni B.K., Weatherill N.P. Handbook of grid generation. CRC Press, 1998.

[10] Liseikin V.D. Grid generation methods. Springer, 1999.

[11] Liseykin V.D., Likhanova Yu.V., Shokin Yu.I. Raznostnye setki i koordinatnye preobrazovaniya dlya chislennogo resheniya singulyarno vozmushchennykh zadach [Difference grids and coordinate transforms for numerical solution of singularly perturbed problems]. Novosibirsk, Nauka Publ., 2007.

Martynenko S.I. — Dr. Sc. (Phys.-Math.), Researcher, Department of Special Engines and Chemmotology, Central Institute of Aviation Motors (Aviamotornaya ul. 2, Moscow, 111116 Russian Federation), Senior Researcher, Department of Combustion and Explosion, Institute of Problems of Chemical Physics, Russian Academy of Sciences (Academician Semenov prospekt 1, Chernogolovka, Moscow Region, 142432 Russian Federation).

Please cite this article as:

Martynenko S.I. Remarks on Generation of the Orthogonal Structured Grids. Herald of the Bauman Moscow State Technical University, Series Natural Sciences, 2019, no. 1, pp. 16-26. DOI: 10.18698/1812-3368-2019-1-16-26