On the classification of discontinuities by the pattern recognition methods

Computers & Fluids Vol. 18, No. I, pp. 35-74, 1990

0045-7930/90 $3.00+ 0.00 Copyright© 1990PergamonPressplc

Printed in Great Britain. All rights reserved

ON

THE THE

CLASSIFICATION PATTERN

OF

DISCONTINUITIES

RECOGNITION

BY

METHODS

E. V. VOROZHTSOV Institute of Theoretical and Applied Mechanics, U.S.S.R. Academy of Sciences, Novosibirsk 630090, U.S.S,R. (Received 27 June 1988; in revised form 30 January 1989) A~tract--It is proposed to use the deterministic discriminant methods of the pattern recognition theory

for a solution of the problem of classifying the discontinuities into several types (shock waves, contact discontinuities, etc.) in the numerical solutions of two-dimensional gas dynamic problems obtained by finite-difference shock-capturing schemes. The formulae for the computation of the features are derived. Three different schemes for the classification of singularities in gas flows are presented. The efficiencyof these schemes is illustrated by the example of a number of applied two-dimensional problems, including transonic potential flow around an airfoil, the diffraction of a strong shock wave on a wedge and supersonic flow in a wind tunnel with a lower-wall step.

1. INTRODUCTION At present shock-capturing finite-difference schemes are widely used in the numerical studies of problems on the dynamics of inviscid compressible non-heat-conducting gas. In the numerical solutions obtained the strong discontinuities are approximated by some transition regions, the size of which in the direction of a normal to the discontinuity surface, is typically equal to several intervals of a spatial computing mesh. This gives rise to a problem of interpretation of the results obtained and to a problem of increasing the accuracy of numerical solutions in the vicinity of discontinuities. The information on the flow structure is in many cases of primary importance. The data on the discontinuity surfaces: their disposition, shape, propagation speed, etc, may be the elements of this information. The methods for localisation of discontinuity fronts on the basis of shock-capturing numerical computations that have been developed until now may be subdivided into two groups: (1) localisation methods the implementation of which is closely related to the use of a priori information on the orientation of the discontinuity surface with respect to the spatial coordinate axes [1]; (2) the methods which do not require the above a priori information for their implementation [2-7]. The role of the methods for localisation of singularities, which belong to the second group, increases substantially in the cases when the phenomena or the processes under study prove to be inaccessible for other investigation methods (e.g. for experimental methods). Consequently this is when the mathematical modeling becomes the only possible means for the investigation of these phenomena or processes [8]. When one discusses the problem of the interpretation of the results of gas-dynamical computations, one usually means not only the determination of the location and shape of the discontinuity surfaces, but also determination of the types of discontinuities. Shock waves and contact discontinuities are such types of discontinuities in the case of the flows of an inviscid compressible gas. It is the work [9] in which it was formulated for the first time that the higher levels of processing and the analysis of the gas-dynamical computation results are to include pattern recognition algorithms and an automatic classification of new objects, that is they should possess specific capacities of the artificial intelligence. In Ref. [10] the pattern recognition and mathematical modeling methods are considered as the two basic methods which were introduced by the informatics into the solution of the problem of prognostication (diagnosis) of the phenomena being investigated by descriptive sciences. Using the terminology of Ref. [11], we shall consider the systems of processing the results of gasand hydrodynamical-calculations, which enable one not only to localize the singularities in the multidimensional flows, but also to automatically classify them into several types; as the intellectual 35

36

E.V. VOROZHTSOV

systems of information extraction. It is possible to indicate at least three ways of using the information obtained at the output of such systems for processing the digital data of gas-dynamical computations. Firstly the implementation of computer systems of an automated analysis of the flow structure enables one to achieve a better understanding of the physics of a flow under study [12]. It was pointed out in Ref. [13] that "flow phenomena such as flow-reversals, shocks, shear-layers and vortices can often be difficult to identify and visualize, especially if the flow is also unsteady. As a consequence, much of the computational aerodynamicist's time is now being devoted to extracting, displaying, and analyzing various features of the solution." The same thought has been expressed in Ref. [14] where it was emphasized that a rapid interpretation of CFD results is of crucial importance for the efficiency of the work of a researcher concerned with numerical modeling of 2-D and, especially, 3-D fluid flows, and that "additional methods are needed to extract the full measure of understanding for CFD data". The author of Ref. [15] reaches a similar conclusion that an active or "smart" software display package is needed "that searches the data base for interesting flow phenomena and displays them. This would require, for example, the program to look for different combinations of the flow variables and their gradients, or derivatives of their gradients, to uncover interesting regions of the flow that might be lurking in the small-scale portions of the mesh, such as secondary, separated flow regions, and call the viewers attention to it and display it". The second way of using the intellectual system of information extraction from CFD data is related to a direct use of the information obtained at the output of such systems, for the control of a process of the numerical solution of the basic aerohydrodynamic problem, with the purpose of increasing the accuracy of computation. A number of approaches to the control of different stages of a computational process have been enumerated in Ref. [16]. Firstly we would like to mention a few works [17-19] devoted to the application of artificial intelligence (AI) concepts to the grid adaptation. In particular, Dannenhoffer and Baron [17 18] report on an expert system embedded in a 2-D grid adaptation scheme. This system in particular contains a knowledge on when and how to refine or coarsen the computational grid (see a description of this system also in Ref. [19]). The expert system EZGrid presented in Ref. [19] is an Expert Zonal Grid Generation system that partitions a 2-D flow field into four-side, well-shaped zones that are then individually discretized. The AI methods, ideas and concepts can also be used in CFD for the choice of a "rational" mathematical model for simulating a specific aerohydrodynamics problem [8, 16, 20, 21]. Some techniques of constructing "rational" numerical models for solving nonlinear and the manydimensional problems of fluid mechanics have been discussed in Ref. [8]. A specific feature of the system described in Refs [20, 21] is the use of a vast data base containing the information on simple engineering mathematical models which fit the flow parameters in different subdomains of viscous gas jets. Finally, there exists the third way of using the AI systems in combination with CFD: the application of AI concepts and expert systems in combination with the numerical solution of aerodynamic problems in the design of new aerodynamic shapes [19, 22-24]. In particular, earlier work [19, 22-23] present the expert system EXFAN (Expert Cooling Fan Design System). This system performs aerodynamic design of turbomachinery components by starting with an initial design and then iterating through analysis and redesign until the design goals are met within the specified constraints [19]. An automated design system MAVR has been presented in Ref. [24], which is now adapted to the design of heating and refrigerating engineering systems. A positive experience of practical use of this system is described in Ref. [24]. A number of requirements for the structure of future expert systems coupling AI and CFD have been formulated in Ref. [25]. It is stressed therein that the main goal of such expert systems is to aid the organization of aerodynamic design. According to Ref. [26] the pattern recognition techniques may be subdivided into three directions: heuristic, discriminant and syntactic. The basis of the heuristic approach to pattern recognition is constituted by the experience and intuition of a research worker. The systems built by such methods incorporate specific procedures developed for concrete recognition problems [27]. Shortcomings of heuristic methods have been enumerated in Ref. [26], in particular, the absence of strict validation, although the efficiency of such a method can be sufficiently high in practice.

Pattern recognitionmethods

37

The basis of discriminant methods (sometimes they are called mathematical methods [28, 29] or vector or geometric methods [30]) is constituted by the use of certain mathematical formalisms for the description of original data on the pattern and by the application of well-known mathematical methods for obtaining a solution. This approach may be subdivided into two classes: deterministic class and statistical class. The methods belonging to the first of these classes do not utilize in an explicit form the statistical properties of the pattern classes under study. The statistical algorithms of pattern recognition are based on the application of the Bayesian classification rule or of its varieties [26, 28, 31, 32]. The syntactic (linguistic, structural) methods [11, 27, 29, 33] are applied in the cases when not only the class to which a pattern belongs is of interest, but also a description of each pattern with the aid of the elements (subpatterns) and their relationships. The pattern can be described with the aid of a hierarchical structure of subpatterns being similar to some extent to the syntactic structure of a language. This enables one to employ in the solution of pattern recognition problems the theory of formal languages developed specifically for the communication between the man and the computer. As was pointed out in Ref. [26], this approach is especially useful in the work with patterns which cannot be described by numerical data, or they are so complicated that their local features are difficult to identify and one has to apply to global properties of the patterns. An earlier work [34] appears to represent the first attempt of using the ideas of syntactic pattern recognition to the flow analysis. The authors of Ref. [34] use features such as critical points and dividing streamlines as a basis which enables them to obtain a representation of the global topology of the flow by a graph with the various structures represented by the nodes and their relationships in the flow by the connecting lines of the graph. Once the flow field has been placed in this form, it can be studied and compared with other data sets using techniques of syntactic pattern recognition. As in Ref. [4], for the solution of a problem on the classification of discontinuities in the 2-D gas flows modeled on the basis of finite-difference shock-capturing schemes, we use the methods of the pattern recognition theory that belong to the class of discriminant methods. Here a question arises: which of the discriminant methods--deterministic or statistical ones--are preferable? We shall consider, as in Refs [2-7], the results of shock-capturing calculation of 2-D problems as digital images (they will be called in the following finite-difference images for brevity). Consider these images from the point of view of the presence of some statistical or probabilistic properties in such images. Finite-difference images differ substantially from the digital images obtained, for example in aerial photography, in medicine and biology, in that the properties of the noise and blur in finite-difference images have under deterministic initial and boundary conditions a deterministic character. Thus they are determined by truncation errors of a difference scheme and by a length of the machine word of a computer employed in computations [35]. This determinism of finite-difference images leads in practice to the fact that the results of two repeated computations of the same gas-dynamical problem obtained on the same computer by the same finite-difference method coincide in all digits of the mantissas of numerical results. Although, as already stated, the finite-difference images obtained in the numerical gas-dynamical simulations have a deterministic character, it is of interest to elucidate a question on the fundamental applicability of statistical pattern recognition methods in the case of finite-difference image recognition taking into account a comparatively wide acceptance of statistical methods in the pattern recognition systems. Investigations of the properties of finite-difference schemes carried out by Vatolin [35-36] using mathematical statistics, the Shannon's theory of information and the probability theory, enable us to answer this question positively. In particular, it was found in Ref. [35] from the requirement of the growth of the correlation coefficients and the criterion for the uncertainty (entropy) reduction in the solution, u n÷~ obtained by an explicit difference scheme for the heat equation u, = auxx the well-known Neumann stability condition az/h 2 <~ 1/2. Where r is the time step and h is a step of the mesh on the x-axis. Applying the method of Vatolin to an explicit difference "corner" scheme with one-sided differences, which approximates the equation u, + aux = 0, it is easy to obtain the well-known Courant-Friedrichs-Lewy condition a'c/h <~ 1. Below we apply for the classification of gas-dynamical discontinuities the deterministic discriminant methods; some of them were presented previously by us in Ref. [4]. A more detailed validation than in Ref. [4] of the chosen system of the features is carried out (see Section 4). The methods for localisation and classification of singularities presented in Refs [2-4] in the case of a uniform

38

E. V, VOROZHTSOV

rectangular spatial computing mesh are extended below for the case of a curvilinear grid. Such an extension is carried out at the example of a problem on the transonic potential flow around an airfoil. We also present a number of examples of automatic analysis of the structure of some 2-D flows in which the intensities of shock waves and contact discontinuities strongly vary in space and in time. 2. I M A G E

FORMATION

The input data for the algorithms, which realize the pattern recognition methods presented in Section 5 are formed in a lower-level image processing system. The basic stages of this system are: (1) the input of a digital image, that is the input of 2-D numerical arrays which present the distributions of the density p, the velocity components u, v and the gas pressure p obtained as a result of the numerical solution of a gas-dynamical problem; (2) image preprocessing which may include one or several procedures from the set: image smoothing, noise filtration, image restoration; (3) image segmentation. An automatic analysis of the data obtained at the output of a lower-level processing system is performed in the higher-level processing system the basic elements of which are: (1) calculation of the features; (2) calculation of discriminant functions; (3) classification of the flow singularities. This system of processing and of automatic analysis of the finite-difference images is conceptually close to the machine vision system of a manipulation robot presented in Ref. [37]. The input of an image into the system of image processing is usually performed from the basic program of the numerical calculation of a gas-dynamic flow. Depending on the mode of using the information obtained at the output of an intellectual system of digital image processing there can be implemented various input regimes. For example, in the cases when this information is to be used in the process of solving a nonstationary problem for the computation control, the processsing of the digital images of the flow field may be performed at each time step. In cases when a stationary problem is solved and the postprocessing of the numerical solution obtained is needed only for obtaining an information on the flow structure, then the input of an image into the system of information processing may be performed only once. The simplest algorithms of image preprocessing are reduced to the smoothing [38]. Note that the smoothing is applied for a long time in gas-dynamical computations. It is usually performed with two purposes: either stabilization of computations by a difference scheme, or the reduction of the amplitude of parasitic post-shock oscillations. Write the Euler equation system governing the flow of an inviscid compressible non-heat-conducting gas in the case of two space variables x, y and the time t as follows:

Ou/~t + c~F(u)/Ox + OG(u)/Oy = 0

(1)

where P

I

pu

U =

pv

pE

I

I P +pupu 2 F(u) =

puv

pu + puE

pv puv

G(u) =

p +pv 2 pv + pvE

(2)

In eqns (1)-(2) E = ~ + 0 . 5 ( u 2 + v2), c is the specific internal energy. In this paper we restrict ourselves to the ideal gas equation of state: p = (~ - 1)p~,

(3)

where y is the ratio of specific heats. In some of the computational examples presented below we have used the MacCormack scheme [39] augmented by the smoothing process from Ref. [40]. The corresponding formula for the smoothing may be written down by using the notions of an image window and of a mask that are widely used in the theory of digital image processing [2-5, 31, 38]. Let f b e any of the quantities p, pu, pv, pE, and let F s ( f ) be the 3 x 3 image window centered in the pixel fj where the subscript i is a number of a cell of a uniform rectangular computing mesh in the (x, y) plane along the x-axis, and j is the number of this cell along the y-axis. Then the formula for the smoothing operator from Ref. [40] may be represented in the form of a discrete

Pattern recognition methods

39

fi~+' = Hs x F3(u~ +')

(4)

convolution:

where fi~+~ denotes a smoothed value, n is the number of time step, and H3 is the smoothing (defocusing) mask of the form H3 = 1/(q~ + 4)

IZ' ] ~p

.

(5)

1

In eqn (5) tp is a constant positive parameter chosen empirically. Some effects of the image restoration have been considered by us in Refs [2, 5] at the example of restoration of piece-wise constant functions. The results of corresponding numerical experiments presented in Refs [2, 5] are in our opinion very encouraging. The problem of the restoration of arbitrary gas-dynamical discontinuous profiles is much more complicated, since in the general case one has to take into account a number of constraints: (1) the laws of conservation of the mass, momentum and energy; (2) the non-negativity of the density, pressure and specific internal energy; (3) the increase of the entropy across the shock waves. 3. IMAGE SEGMENTATION Let f ( x , y) be any of the functions sought to enter the system eqns (1-3) and are computed for fixed t. It was proposed in Refs [2-5] to consider the f u n c t i o n f f r o m the point of view of the digital image processing theory, that is as the function of intensity (or brightness) of the image; further, it was proposed to consider an edge in the image brightness as an analog of a gas-dynamical strong discontinuity smeared by a shock-capturing computation. The localisation of the edges in brightness constitutes as a rule the main goal of the algorithms for the segmentation of the digital images. At first let us briefly present the segmentation of the digital images given on a rectangular uniform mesh in the (x,y) plane. For this purpose, we have used in Refs [2-6] the Sobel edge detector, which is one of the most popular gradient edge detectors. The magnitude g,j of the function f gradient at a point (i, j ) is computed with the aid of this detector by the formulae: S, = H, x F3 ( f ) ,

$2 = (h2/h,)H2 × r 3 ( f ) ,

g,j= IS, I + 1521

(6)

where H~,//2 are the Sobel masks: [ Ht=

1 0

-1

0 -2

,

H:=

-

-2 -1

0 0

,

(7)

hi, h2 are the mesh steps along the x-axis and the y-axis, respectively. The orientation of the vector of the gradient at the centre of the (i, j ) cell is computed by the formula ~t~= arctan(Sj/$2).

(8)

At the first stage of a procedure for the edge detection the magnitudes of the gradient gij and its orientation % are computed at all the points (i, j). At the second stage the points (i, j ) are determined which satisfy the threshold limitation g¢> T ~'). The value T tv~ was determined iteratively. Let N~ be the overall number of the pixels in an image. The zeroth approximation T t°~ was calculated as an arithmetic mean: T I°) = (l/N/) ~ go,

(9)

i,j

where the summation was carried out over all the pixels of a digital image. Let ~v~ be a set of the pixels f j of the image, for which the inequality gu > T~v)is satisfied, and let N[ v) be the number of pixels entering in y~"), v = 0 . . . . . Ng (0 <<.N g ~ 2). Then the next approximation T tv+l) was computed as: Tt"+')= (1/N~ ~')

~

(i, j) ~ ~l(V)

gij,

v = 0 . . . . . Ng.

(10)

40

E.V. VOROZHTSOV

According to Ref. [4], the use of more than one iteration (that is N e > 0) by the formula (10) may be recommended (with the purpose of reducing the number of artefacts) in the cases when it is known a priori that the intensity of all strong discontinuities arising in the flow is approximately the same. Even in these cases the value Ng = 2 was found in practice to be the maximum allowable value of Ng, because the number of pixels satisfying the threshold limitation go > TIN'~ 1 rapidly decreases with further increase of Ne. At the third stage the thinning of the edge contours is performed by using the nonmaximum suppression method of Rosenfeld [5, 41], which represents a rather natural and straightforward procedure. In this procedure at first the one of the eight principal directions which is the closest one to the direction determined by the angle ~ , see formula (8), is determined. Since the angle ~i, determines the direction normal to the edge orientation, the found principal direction approximates the true direction of a normal to the edge and, in addition, it enables one to determine two neighbouring pixels (i~,j~) and (it, Jr) closest to (i, j ) and approximately lying on the normal to the edge and on different sides of it. If the inequalities: g i / - gi~.j~ > 0,

g,1 - gi~• I~ > 0,

( I 0a)

are simultaneously satisfied, the node (i, j ) is included in the set of edge points. This simple approach eliminates the competition of neighbouring pixels located on the edge, i.e. the node (i, j ) at which the inequalities eqn (10a) are satisfied will not be eliminated from the set of edge points in the case when the neighbouring pixel lying on the edge possesses a higher value of g than the pixel (i, j). At the fourth stage the isolated artefacts are removed with the aid of the following simple procedure [4]: the point (x~j, y!j) is considered to be an isolated artefact, if they are within a circle of the radius R = ~ .(h~ + h~) °5 no other edge points. Here 0 is a constant multiplier; we have used, as in Ref. [4], the value 0 = 1.5. A more detailed presentation and discussion of the above multistage procedure for the segmentation of digital images on a rectangular uniform grid may be found in [2, 5]. Curvilinear grids are presently widely used in the solution of problems in the domains with complex geometric form [1, 19, 42]. In this connection it appears that the problem of development of the methods for automatic analysis of difference soituions computed on the arbitrary curvilinear grids (both adaptive and non-adaptive ones) is of present interest. Let us present an extension of the above image segmentation procedure for the case when the coordinates of the points of a curvilinear grid with the aid of the functions: x = x(~,,7),

Y =Y(~,'7),

(11)

can be uniquely mapped onto the points of the lines of a uniform rectangular grid in the plane of curvilinear coordinates (~,r/). Let J be the Jacobian of the transformation eqn (11). Then: (12)

J = l/(x¢y,~ - x,7)'~)

where x~ = ~x(~, r/)/0~, y, =Oy(~,r/)/Oq, etc. Let f ( x , y ) arguments. Then:

be a differentiable function of both

f , = J . (f~y, - f , y + ), L = J . ( f , x ¢ - f ~ x , ~ ) .

(13)

Employing the image window F 3 ( f ) and the masks of the Sobel edge detector eqn (7), we can obtain with regard for an analysis of Refs [2, 5], the following difference approximations of formula (13): L = (1/8)2 • (y~'S2 - y ~ . S ~ ) ; L = (1/8)J • (x~. St - x . . S z)

(14)

where S ] = H t x F 3 ( f ) ~ - 8 o f / o r l; $2 = 112 x F:,,(f ) ~- 88f/8~.

For the computation of the derivatives x¢, y , in eqn (14) one can use the formulae: x~ = (1/8)H2 x F3(x),

y, = (I/8)H~ x F.,(y),

(]5)

Pattern recognition methods

41

etc. Taking into account a relative complexity of these formulas we have used in practical computations in the interior nodes (i, j ) more simple formulae [42]: x~ ~--- 0.5(Xi+ I.j -- Xi-- l,j);

y, = 0.5(yi.j+, - yl.j_ t),

(16)

etc. Let the nodes (xq, y0) be numbered as follows: 1 ~< i ~ i max, 1 ~
i = 1. . . . . i max;

i=2,...,/max-1; i = 1. . . . . /max,

(17)

etc. in the boundary nodes. The basic ideas of the segmentation algorithm for an image given on a curvilinear grid will be illustrated at the example of a problem on the computation of a 2-D transonic potential flow around an airfoil NACA 0012. The merits and shortcomings of the mathematical model of potential flow have been discussed in a number of surveys, see for example [43-44]. Here we note only the following limitation of this model: it is generally applicable only for the description of those transonic flows in which the maximal value of the Mach number M = IUI/c (IUI is the magnitude of the gas velocity, c is the local speed of sound) does not exceed the value 1.2 [45-46]. This means that only shock waves of a relatively weak intensity are allowed in a potential transonic flow around an airfoil. In Ref. [47] interesting considerations are presented in favour of using the numerical simulation on the basis of the full equation for the velocity potential go(x, y) in the design of efficient shockless transonic profiles. This equation may be written in conservative form as follows: (Pgox)~ + (P%).,. = 0,

(18)

where x, y are rectangular Cartesian coordinates, and the u, v components of the velocity vector U are expressed by the potential go as follows: u = ~go/Ox,

v = ~go/~y.

(19)

The density p is expressed with the aid of the Bernoulli integral in terms of go~, % by the formula: p = [ l _ , _ 7 + 11 (go.2+ go2)]-,;I;.- ,),

(20)

where the density p and the velocity components go.y,go,. are nondimensionalized with respect to the stagnation density Ps and the critical speed of sound c , , respectively. The full-potential equation written in ¢, q coordinates has the form [47]: (pU/J)~ + ( p V / J ) , = O,

(21)

where: p

L[-1- ~ + - 11 (Ugo~+ Vgo.)j /T/(~'-')

U=Atgo¢+A2go~, AI = ~x~ +. , . , 22 ~x=ly~,

V=A2go~+A3go,,

A2 = ~.JL + ~y/7,.,. A3 =/7.~t /+y , 22

~,.=-1%,

rlx=-lye,

q.,.=lx~,

(22)

J is the Jacobian of the transformation eqn (11) given by the formula (12). For the generation of a curvilinear computing mesh around an airfoil we have used the method of Ref. [48]. This method enables one to generate a C-type grid with the aid of a transformation x = B + A cosh(q)cos ~, y = A sinh(~/)sin ~ where A, B are constants. Applying the algorithm of Ref. [48], one can create with the aid of a number of user-specified parameters the concentrations of the curvilinear grid lines in the domains near the profile, near its leading and trailing edges, see Fig. 1.

42

E.W. VOROZHTSOV

(a)

(b) Fig. 1. (a) The partial view of the mesh; (b) the C-type 95 × 15 grid around the NACA 0012 airfoil.

We have used a modified approximate factorization scheme AF2 of Hoist [47, 49] for the numerical solution of eqn (21). In the process of computer implementation of this scheme we have taken into account some ideas on the improvement of this scheme that were proposed in Refs [50-53], namely: consistent approximations of the derivatives of x, y and ~o with respect to ~ and of the derivatives of x, y and q~ with respect to r/have been used which ensure the satisfaction of the conditions for the exact reproduction of a homogeneous flow when it is calculated by the AF2 scheme [50-51]. We have also implemented the improved intermediate boundary conditions on the airfoil surface which are derived in Ref. [52] (see also Ref. [53]). These conditions were derived in Ref. [52] as the result of a systematic treatment of the problem on the effect of the intermediate boundary condition on the convergence of the AF type schemes. Our choice of the AF2 scheme for the numerical integration of eqn (21) is related to the fact that this scheme possesses a higher computational efficiency (in terms of a number of iterations needed for convergence and of computer expenses for the calculation of the solution q~ at one iteration) in comparison with other numerical methods of solving the full-potential equation which were developed during the period of 1978-1985 [51, 54]. In Fig. 2 we present the examples of comparison of the pressure coefficient cp graphs with the results obtained by Lock [55], Fig. 2(a), and by Jameson [56], Fig. 2(b). Note that in the Hoist method there are three user specified parameters. The first is an empirical fix used for the intermediate boundary condition, the second, co, is a relaxation factor to improve the convergence and the third is a parameter controlling the amount of the artificial viscosity.

Pattern recognitionmethods

(a)

43

(b) ~---h - 0.8'

-0.6- ~

~ - 0,6'

-0.4 -

-0.4.

-0.2"

-0.2.

-0.0

,

-O.O

0.2.

0.2.

0.4-

0.4.

,

,

V__ - - ~ . . ~.,,~

0.6. 0,60.8" 0.8" 1.0. 1.0-

Fig. 2. Pressure coefficient cp as a function of x. (a) The freestream Mach number M~ = 0.72, the angle of attack ~ = 0: ( ) the present method; (. . . . . ) Lock, 1970; (b) M~ = 0.8, e = 0. ( ) The present

method, 95 x 15 mesh; (. . . . ) Jameson, 1981; 256 x 64 mesh.

Recently, a number of other methods for solving eqn (21) have been proposed, see for example earlier work [54, 57-64]. In Ref. [54] a parameter free conjugate gradient type method was introduced to avoid the user specified relaxation factor in Hoist original work. The third parameter in the Hoist Code, controlling the amount of artificial density, can be replaced by a flux biasing scheme, satisfying an entropy inequality [57]. The authors of Refs [58-60] also avoid using the artificial compressibility method in AF type schemes, because they have shown that this method gives a significant overestimation of the gas flow rate while calculating the nozzle flows. The AF schemes proposed in Ref. [58] contain only one user specified parameter ~o = 4/3. The stability of these schemes was studied in Ref. [58] by considering their first differential approximations (or modified equations, in another terminology). The authors of Ref. [61] solve eqn (21) with the aid of finite-element method and the multi-grid method of Fedorenko. In addition, the equations of shock polar are used in Ref. [61] to increase the domain of agreement between the potential equation solutions and the solutions of the Euler equations in cases of transonic flows with shock waves. The AF2 scheme of Hoist was also used in Refs [62-64] for the computations of transonic airfoil and cascade flows on the spatial curvilinear grids of the H-, C-, and O-types. For example, in a computation of a cascade flow on a 64 x 32 H-type mesh the residual dropped by 5 orders after 90 iterations [63]. To reduce the level of classification errors while detecting the shock waves in transonic potential flows with the aid of pattern recognition method whose different stages are presented in this section and in Section 5.2, it is desirable to choose a numerical method for solving eqn (21), which captures shocks without wiggles. The AF2 scheme meets this requirement, as this may be seen in Fig. 2(b). This requirement can be weakened, if one makes use of an image preprocessing stage (see at the beginning of Section 2). However, it should be kept in mind that the implementation of an image preprocessing, especially image restoration, involves additional computer time expenses [2, 4, 5]. At the first stage of a procedure for segmentation of an image given in the nodes of a curvilinear grid the gradient g[x(~, rl), y(¢, r/)] of the image was computed by the formulae of the Sobel edge detector eqns (7), (14), (15) where one should set f = p, p is the gas density. In Fig. 3 we present the gradient surface g[x(~, q),y(~, ~/)] mapped onto a rectangular computational domain in the (~, r/) plane. The shaded parallelogram is placed in Fig. 3 under the nodes ( x , , y , ) lying on the airfoil surface. We can see from Fig. 3 that the use of one-sided and central differences of the form eqn (17) in the boundary nodes, leads to a smooth transition in the image from the boundary points to interior points which justifies the use of the formulae eqn (17).

44

E.V. VOROZHTSOV 12

i q9

) 95

Fig. 3. Dimetrics of the gradient surface g[x(~, r/), y(¢, ~1)].

Since the shock waves which may be present in a transonic flow around an airfoil are characterized by a relatively weak intensity, it is reasonable to perform only one iteration eqn (9) while calculating the threshold T. For the application of the Rosenfeld's nonmaximum suppression method it is first necessary to number the principal directions. The scheme of numbering the principal directions in the case of a rectangular uniform mesh was presented in Refs [2, 5]. In the case of a curvilinear grid we associated the kth principal direction with a certain node of a curvilinear grid independently of the orientation of the grid lines in the neighbourhood of a node (x~, Yij) under consideration with respect to the axes of a Cartesian coordinate system (x, y), see Fig. 4. For example, by number 3 we encoded a direction from the node (i,j), to the node ( i - I , j + 1), etc. Denote the corresponding angles by a~k~ k = 0.1, 7. For example: gl3) = arctan[(yi_ ,.j+l - Yu)/(xi

l.j+l

-

-

Xij)];

if a~) < 0, then we replace the value a~) by -u ~7~kl- a-~~kl+ 2n. As a consequence of this, the quantity ~ ) undergoes a jump of the magnitude 2n when going from the negative values ~I~I to the positive ones. To determine more exactly the number of the principal direction in this situation, we

i+ i+

1,J-

I

1,J

i + 1,J + 11

i,J

+ 1

i-

-

i-

1,j - I

1,j

1,j + 1

Fig. 4. Principal directions numbering in the case of a curvilinear grid.

Pattern recognition methods

45

introduced an auxiliary ninth direction determined by the angle:

~

I°) + 2~,

ctts)u= J 0,

L~01

if - (0) ~< A~/2 if 2rr - ct~) ~< Ae/2

(23)

otherwise.

The increment Act in eqn (23) was set as an arithmetic mean of the differences ~(k)_ ~,(k-I) k = l , . . . , 8. It is obvious that Act = r~/4. Figure 4 explains how the introduced principal directions are used in the edge thinning process. The dashed line in Fig. 4 shows the orientation of a normal to the shock wave front, passing through the point (i, j). The orientation angle ~to of this normal is computed via the formula: ~u = arctan[(c~p

It?>,)/(~p/Ox)],

(24)

where the derivatives Op/t?y, Op/Ox are approximated by the eqns (14), (15) in which one should set f - p. Similarly to the case of rectangular grid, the angles e(ik), k = 0 . . . . . 7, are used for the determination of two neighbouring nodes with indices (i~,jt) and (it,jr) closest to ( i , j ) and approximately lying on the normal to the edge and on different sides of it. For example, in the case shown in Fig. 4 we should set in eqn (10a) i~= i + l,j~ = j + 1, ir = i - l,jr = j . Figure 5 shows the points obtained after execution of the first (a), the second (b) and the third (c) stages of the segmentation of a digital image given on a curvilinear 95 x 15 mesh presented in Fig. l(b). The fourth stage was not realized by us in the case of a curvilinear grid. As we shall see in the following, the classification stage enables one to exclude from the set of points, obtained as a result of executing the three segmentation stages (see Fig. 5c), a large number of points and as a result of this a set of shock wave fronts is obtained in which there are no isolated artefacts. This is related to the fact that the AF2 scheme used by us enables one to obtain numerical solutions which do not contain parasitic oscillations in the neighbourhood of a shock wave, see Fig. 2(b). 4. FEATURE SPACE The points (xm, y,,), m = l . . . . . N4 (N4 > l) which are obtained at the segmentation stage and which may potentially belong to the lines of strong discontinuities, are along with the 2-D arrays p, u, v, p the input data for the pattern recognition procedure. As was pointed out in [28], the choice and extraction of features play a central role in the pattern recognition. Let us enumerate a number of the requirements for the features which in the literature on pattern recognition are taken into account while determining a system of features ensuring a possibility of an efficient objects classification by the algorithms of automatic recognition. While formulating these requirements we shall immediately try to "attach" them to our specific problem of recognition and automatic classification of singularities on the basis of the results of a shock-capturing calculation of 2-D gas dynamics problems [4]. 1. Feature invariance [26]. This property means that all the objects belonging to the same class (for example, to the shock waves class) should possess this feature. 2. High information content of the features. There are various quantitative and qualitative definitions of the information content of the features. It was proposed in Ref. [65] to use the probabilistic measures of the information quantity for the quantitative characteristic of the information content of the features (the measurements of Hartley, Fano, Kotelnikov, Shannon, Kullback, etc.). In Ref. [26] the feature of an object is considered to have a high information content. The removal of this feature from the classification algorithm increases substantially the measure of uncertainty in the classification. This is similar to the information criteria underlying the definitions of information content in Ref. [65]. Some techniques for the construction of the features having a high information content were discussed in Ref. [31], in particular, the means of the factored analysis. An algorithm for an automatic estimation of the quality of the features has been presented in Ref. [66], where the methods of the theory of fuzzy sets have been used. 3. Independence of the form of the functional dependencies to which the features obey on the specifics of the gas-dynamical problem being solved. This property determines in the end the universality of a system of pattern classification.

E. V. VOROZHTSOV

46

.•

(a)

.::

×

~××

j

~ :

oao

.

•

•

°

°

,oo

° ° ° °

.

j

k x ×

×

x~xx \~xxxx

~

Fig. 5. The results of the execution of the first three stages of the procedure for the segmentation of the p(x, y) image. ( ) Sonic lines and the airfoil contour.

Pattern recognition methods

47

4. The requirement of a minimal overall number of features [31, 67, 68]. As already mentioned, the features may be subdivided into features with rich information content and features with little information content. The inclusion into a recognition system of a too large number of features, firstly increases the computer time expenses and the requirements for the computer memory resources; and secondly it may lead in some cases to the deterioration of classification [31, 68]. 5. Feature invariance with respect to the location and orientation of the objects with respect to the spatial coordinate axes [37]. 6. The features should possess a sufficiently high noise immunity in order to provide a correct classification of singularities also in the presence of the numerical solution oscillations in the neighbourhood of discontinuities. The above 6 requirements are directed in the end to the goal of diminution of the level of classification errors under simultaneous satisfaction of the requirement of the universality and computational efficiency of a recognition system. We have tried in Ref. [4] to take into account the above requirements while deriving the formulae for the features applicable for the recognition of discontinuities in the 2-D gas flows modeled on the basis of finite-difference shock-capturing schemes. We used for this derivation available information on the mathematical properties of solutions of the equation system (1)-(3) that take place in the subdomains of a continuous flow, and an algebraic system representing the dynamic compatibility conditions at the discontinuities [69, 70]. It follows from the consideration of these properties, in particular, that it is reasonable to use the gas density p as an image intensity function f in the Sobel edge detector eqns (6), (7). Suppose that we have obtained, with the aid of a lower-level system of processing of the gas-dynamical computation results, a set of objects--the points (Xm, Ym), m = 1. . . . . N4. Among these points one can also find such points which are indeed located in the subdomains of a continuous flow, for example, in the domain of a rarefaction wave. From this it follows that one should take into account in a system of recognition and automatic classification of the objects (Xm, y,,,) also the object classes corresponding to such subdomains. We shall consider, as in Ref. [4], the following classes: o9~, the class of normal shock waves; e92, the class of oblique shock waves; 0)3, the class of tangential discontinuities; 0)4, the class of purely contact discontinuities; 0)5, the class of compression waves; ~o6, the class of rarefaction waves; 0)7, other subdomains of continuous flows. Here by a tangential discontinuity and by a purely contact discontinuity we mean, following Ref. [70], two possible types of contact discontinuity in an ideal gas. Let us denote by p~, p~, E~, u,,, u,~ and P2, P2, £2, u,2, ur2 the values of the pressure, density, specific internal energy, normal velocity component and tangential velocity component, respectively, on different sides of a contact discontinuity. Then in the case of a purely contact discontinuity the following relations are valid:

unl=u,2,

ttrl=Ur2,

pl=P2,

(pl--p2)2q-(EI--E2)2#0.

In the case of a tangential discontinuity the relations:

u,,j=u,,2,

u.~lvau,:2, pt=p2,

(pl-p2)2+(¢l-E2)2~:0,

are valid. The purely contact discontinuities occur relatively seldom in applied problems of gas dynamics. Such discontinuities may be present for example in a flow in which a jump in the density or in the specific internal energy of an isobaric fluid is advected by a constant velocity field, see Ref. [5] for a more detailed discussion and computational examples. Since the normal and oblique shock waves are the kinds of shock waves, it is generally possible to unite the classes co~ and 0)2 into the one "big" class f~, of shock waves. Similarly, the classes o93 and 0.)4 may be united into one class f12 of contact discontinuities and it is possible to interprete o93 and 0)4 as the subclasses of the f~2 class. And, finally, the classes 0)5, 0)6, 0)7 can be united into a class f~3 of the points belonging to the subdomains of continuous flow. In its turn, the classes f~ and f~2 have a common property that they include only the points belonging to the discontinuity surfaces. Therefore, it is worthwhile considering a class C~ of all the points of discontinuity surfaces. Then f~ and ~2 prove to be subclasses of the class C~ and 0)~, 0)2, e93,0)4 are subsubclasses of C~. Such a hierarchical grouping of classes, as we shall see in the following, may be put into the basis of hierarchical systems of discontinuities classification. C A F 18 I - - D

48

E.V. VOROZHTSOV

Let the coordinates of the (xm, Ym) points coincide with the coordinates of a geometric centre of some cell (i, j ) of a rectangular computing mesh in the (x, y) plane. Consider along with the (i, j ) cell the cells with indices (it, Jr) and (it, Jr) which lie on a normal to the discontinuity line passing through the (x,,, Ym) point and on different sides of the discontinuity line. The above indices are determined at the thinning stage of the procedure of the segmentation of a digital image p (x, y, t), see the details in Ref. [2, 5]. Taking into account the property of an increase in the entropy S across a shock wave [69, 70], let us ascribe to the indices (it, Jr) and (it, Jr) the following meaning: if S~,.j~> St,jr, then we leave the above indices unchanged; otherwise the pairs of indices (it, Jr) and (i, Jr) should be interchanged. Thus, the pair of indices (it, Jr) always corresponds to a larger entropy. In practical computations this procedure in the case of using the ideal gas equation of state eqn (3) was implemented as follows. Consider a well-known relationship: exp[(S - So)/Cr] = p /( p~),

(25)

where So is an arbitrary constant and cv is the gas specific heat under constant volume. Let S, p, p in eqn (25) be differentiable functions of some variable ¢ and let A = p/(p:'). Then it follows from eqn (25) that s i g n ( d S / O ~ ) = s i g n ( O A / # ¢ ) . Thus, from the inequality Sg,.j, >S~r.j, it follows the inequality Aig.h > A,,.j, and conversely. Strictly speaking, the inequality S~,j~ > Sg,/, is ensured on a shock wave, if Sg~,j~ = S~, S~,~ = $2 where the subscripts "1" and "2" refer to the medium states behind and before the front of an original nonsmeared shock wave, respectively. Let us show that the inequality A~,,j, > A~,j~ can be violated, if the values of the quantities p, c at the points with the indices (it, Jr) and (it, Jr) deviate substantially from the values p~, c~ and P2, e2, respectively. In the following we shall write for brevity At and Ar instead of A~:.~, and Ae,.j. It is not difficult to show that: A t - Ar = (7 - I ) E r P ) ;[~t/~r

(Pt/Pr)" - t].

--

(26)

From the properties of the Hugoniot adiabat it follows that at the shock wave p~ > P2, q > q [69]. Suppose that the components p, E of a difference solution have a monotone behaviour in the zone of smearing of a shock wave. Then in this zone also the inequalities: ~l > ~r,

Pl

> Pr

(27)

take place. The relationships eqn (27) are insufficient to provide the satisfaction of the inequality At > At. In fact, it follows from eqn (26) that in the case when: 1 < Et/Er < (pl/p,.) ~' I, we have that At < A~. Such a phenomenon was observed in the computations of some 2-D shocked flows. In connection with the foregoing we have realized a search procedure for determining such values Pt, Pt,Pr, Pr which provide the satisfaction of the entropy inequality At > At. A search for the values Pt, Pt and Pr, P; is subdivided into two stages. At the first stage the quantities Pt and Pt are sought; if such values Pt, Pt that At > Ar are not found, one proceeds to the second stage, the stage of search for Pr and p~. Let us introduce the image windows F ~ ( f , i, j ) and F+(f, i, j ) of the form: I f i , j + v~(m) f i + v~(m), j + v~(m)~

Fm(f,i,j)=

f.i

fi+v~{m),i

l-I

f i + #~(,,),j+ u,(,n)1 /

F ~ ( f , i, j ) = ~ f.j+~,,I,,,l

L ~.j

f t + ~v(m),)

-I

where v,~(m)=m.(it--i),

vv(m ) = m . ( j j - j )

#.,:(m) = rn .(i,. -- i),

,u, (rn) = m "(j,. - - j ) .

Let us also make use of a smoothing mask H4 of the form: H4 = (1/4)

[' :l 1

'

(28)

Pattern recognition methods

49

Then we compute sequentially the quantities:

p~l)=pil.jl,

p~1)= pil.j~;

p~:)= Ha x F=z(p, i,j),

p}2)= Ha x F - i ( p , i,j);

p~3) = 1t4 x F?(p, i, j),

p/3) = n 4 x F~-I(p, i, j);

pl 4) = H4 × F{(p, i, j),

p~4) = H4 × F { ( p , i, j);

pl 5) = 0.5(pi,.g, +pij),

p~5)= 0.5(p~,.j~ + p~j);

pl 6) = 0.5(p~.j~ +Pi+,.,:(l).y+Vy(O, p~6)= 0.5(p~,,j~ + P~+vx(l).j+~,,.O))"

(29)

As soon as the value p[k), pl~) is found, such that A~k) > Ar the further search is stopped. If such a value is not found, we setpt = p~k0),Pt = p~k0)where k0 is such a number, 1 ~ A~k) at k ~: k0, 1 ~
p~O = Pit.j,,

p~l) = Pi,.j, ;

p~2)= H4 × F++_,(p,i,j),

p~2)= H4 x F+_~(p,i,j);

p~3)=//4 x F~-(p, i,j),

p~4) = H4

X

p~3)=H 4 × F~-(p, i,j);

F~(p, i, j),

p~5) = 0.5(p~,j r +p~j),

p~4)= 1-14 × F~(p, i, j); p~S)= 0.5(p,,.j, + p,j);

p~6)= 0.5(p~.jr + P~+,xo).J+,,.~)), p(6) , = 0.5(pi,.j, + Pi+~(I).j+~,~l)).

(30)

As soon as the value p~k),p~ _(k) is found, such that At > ---(k) ~ , the further search is stopped. If we connect by straight lines segments the cell centres of a spatial grid, involved in the computations by eqns (29) and (30), we obtain an array of rectangular contours shown in Fig. 6. In cases when the edge passing through the point (i, j ) is parallel with the x- or y-axis, the rectangular contours depicted in Fig. 6 degenerate into straight line segments. Of course, the use of contours shown in Fig. 6 is justified only in the cases when all the points of a contour involved in the computation of a quantity, for example, pl 4) or p~3), lie on one side of a generally curvilinear discontinuity line. The above search procedure was applied only in the cases when et > E,, p, > Pt, but At < A,. In the remaining cases we have used simpler formulae for the computations of At and A,:

At= pt/(p~), .4r= Pr/(p~).

(31)

If At> A,, the pairs of indices (it, jr) and (i, jr) were left unchanged, otherwise they were interchanged. Below we introduce, for a quantitative characteristic of different types of the singularities in inviscid gas flows, seven quantitative features a~. . . . . . a7,n of the objects ( x , . , y , . ) , m = 1 . . . . . N 4 .

/ .,,,-q

q

-/

f f f

J ,,O' 1

/\

L.-"

Fig. 6. Contours for the search interpolation procedure.

50

E.V. VOROZHTSOV

Here we basically follow our work [4], but with a more detailed validation of the formulae for the computation of the features proposed in Ref. [4]. Thus, we consider the derivation of the first feature a~,,. We shall make use of a well-known relationship (see for example Ref. [70]):

un2- u,,1 = 2 ( s i g n j ) c 2 ( M 2 - 1/M2)/(7 + 1)

(32)

where M2 = l u . 2 - D I/c2 > 1 . j = p~(u.~ - D ) = P2(U.2- D), D is the magnitude of the shock wave velocity, c is the sound speed, the subscript " n " by u shows that u. is a gas velocity component being normal to the discontinuity surface. Consider at first the c a s e j < 0 in eqn (32). Then it follows from eqn (32) that:

u.2 - u.t < 0.

(133)

Suppose that n is a coordinate measured along a normal to the shock wave at a point (x... y~). such that n increases while moving from the state "1" to the state "2". Suppose also that the velocity component u. changes monotonously in the zone of smearing of a shock wave. Then with regard for eqn (33) we can write the inequality:

Ou./On < 0.

(34)

Taking eqn (34) into account we shall compute the feature arm by the formula:

a~., = 1 - sign(Ou./On)o.

(35)

The quantity (Ou./On) o was computed in Ref. [4] by the formula

(Ou./On )ij = (Ou,,/Ox)ij cos :t,j + (Ou./Oy )iJ sin s 0.

(36)

where c%. is the angle between the normal to the discontinuity surface at the point (xij, Yo) and the positive direction of the x-axis. For the calculation of 8u./Ox, 8u./Oy we have used two versions of computational formulas. In the first version being simple in terms of computational efforts the derivatives Ou./Ox and 8u./Oy were approximated by the simplest one-sided differences:

(Ou./Ox) o = (1/h,)[(u.)i+ , . j - (u.)~], (Ou./Oy ),j = (1/h2)[(u.)~.j + , - (u.)0],

(37)

where (U.)k.t = UktCOS Ct0 + Vkt sin ct~j. In the other version the edge detector of M6r6 and Vassy [2, 4-5. 71] was used for the computation of Ou./Ox, Ou./Oy:

(Ou./Ox) o = [1/(2hl)]H5 x Fz(u., i, j); (Ou./Oy)o = [1/(2h2)]H 6 × F2(u., i, j),

(38)

where F2(u.. i, j ) is the 2 x 2 image window of the form:

F2(f,i.j)=I

f'f"+~.J f+"J+'].~+,.J'

and /-/5,/46 are the masks of M6r6 and Vassy: Hs=

E~l ll]' H6=[ l 1] 1

-i

-1

The practical use of the formulas eqns (36)-(38) at the smeared shock waves in the 2-D flows showed that the inequality 8u./On < 0 was always satisfied independently of the sign of the quantity f in eqn (32). To explain this phenomenon let us turn to the formula similar to eqn (8): tg ~to.= (Op ~By )/(Op /Ox ).

(39)

Consider two cases of the propagation of a shock wave (oblique or normal wave) depicted in Fig. 7. Let the point (x, y) in the case of Fig. 7(a) lie on the shock front line. Suppose that at this point there is no intersection of the shock wave under study with other shock wave or with a contact discontinuity. Note that pt > P2 since the density p increases behind the shock wave front. Let us

Pattern recognition methods

51

set two small positive quantities Ax and Ay and consider the differences:

Axp = p(x + Ax, y) -- p(x -- Ax, y) Ayp = p(x, y + Ay) - p(x, y - Ay). If within a small neighbourhood of the point ( x , y ) there are no discontinuities of the flow, the function p(x, y) is continuous near the different sides of the shock wave under study, therefore: lim A,p = lim Ayp = Pl - P2 > 0.

Axe0

A)'~0

F r o m this it follows that at sufficiently small values of the steps h I and h2 of the spatial grid in eqn (39) we always will have dp/t3y > O, 8p/t3x > 0 if the difference solution for p changes monotonously in the zone of a smeared shock wave. Thus, the unit normal vector n = {cos ~0, sin ~tu} will have a direction being opposite to the direction of the velocity vector D of the shock wave, see Fig. 7(a). In the case under study f > 0 in eqn (32), so that un2 - un~ > 0. In the image segmentation procedure described in Section 3 the coordinate n always increases in the direction indicated by the vector n. Let the value n = 0 correspond to the (x, y) point at a shock front. Let us set a small increment An > 0. Then: sign(t~u~/t3n) = sign {[u~ (An) - un ( - An)]/(2An )} = sign(u,i - u,2) = - l, that is c~u,/t3n < 0 in the case shown in Fig. 7(a). Considering similarly the case f < 0 presented in Fig. 7(b), we can write with regard for eqn (33) that:

sign(,:3u,/c~n) = sign(u,2 - u~l ) = - 1. Thus, we have shown that in the zones of smeared shock waves aIm = 2 under the condition of a monotone behaviour of grid functions p, u, v in these zones and in the absence of other discontinuities in some neighbourhood of the point (x, y) under consideration. Let us now proceed to the construction of the second feature a2m. Let the function F(p, E) in the equation of state p = F(p, E),

(40)

be continuously differentiable at p >I 0, E >i 0 and let:

dF/dp>O,

dF/&>O

atp>0,

E>0.

(41)

As is known, the condition of the pressure continuity p~ = P2 is satisfied at the original "ideal", that is nonsmeared, contact discontinuity. Employing the Taylor formula we obtain:

F(p~, El) = F(p2, E2) + ((3F(~, ¢)/(~p)'(p, - P2) + (dF(~, E)/t~E)'(E, -- E2),

(42)

where ~ = p t + O ( p 2 - - p ~ ) , ~ = E ~ + O ( E 2 - - E ~ ) , O < O < l . Since at the contact discontinuity F(p,, c , ) = F(p2,E2), we obtain from eqn (42) the relationship: (OF(D, E)/3p)'(p, - P2) + (t~F(fi, E)/OE). (Ej - E2) = 0.

~(1) S ~" 5"

5'

(a)

(b) Fig. 7. (a) f > 0; (b) ~ < 0.

(43)

52

E.V. VOROZHTSOV

In the case of satisfaction of the inequalities eqn (41) it follows from eqn (43) that at the contact discontinuity: sign(pl -- P2) = -- sign(el -- ~2), so that: sign(p~ - p~) + sign(¢~ - ¢~) = O.

(44)

sign(p1 -- P2) + sign(El -- E2) = 2,

(45)

In the case of a shock wave:

since p~ > P2, Ej > E2 at a shock wave. This consideration leads to a conclusion that it is reasonable to calculate the second feature by the formula: a2m

=

sign(flip,J, - fi,r, Jr ) + sign(Ei~,J, - E-i,.Jr )

(46)

where: ~,,,j, =P,,,j,/[~,,,j,(7

-

1)], C,r,jr

=P',,Jr/[&,Jr(~

--

1)],

in the case when the ideal gas equation of state eqn (3) is used, and the quantities fi, fi at the points with indices (i, j~) and (i~, Jr) are calculated in a general case by eqns (29) and (30). Note that the situation fi~,,j~=fii,.j;(,,,j~ =(i,,j, is excluded automatically, because at the (xm,ym) points a threshold limitation gq > T ~v) is satisfied where T(")> 0, which means that there are among the (xm, y,,) points no points located in a subdomain of constant flow. In this sense one can assert that the classification of pixels of the original digital image starts already at a stage of image segmentation with the aid of which the domains with nonzero gradients of the gas density are separated from the domains of constant flow. Let us proceed to the construction of the features a3m , a4m. By virtue of the continuity o f p and u, across a contact discontinuity the quantities (1/p)l~p/c~nl and IOu,/Onl should remain finite. At a shock wave of finite intensity of the order O(l) the quantity (h~/p)lOp/Onl computed on the basis of the difference solution will also be a quantity of the order O(l). At the same time, at a contact discontinuity lap~an I = 0(1), consequently (h~/p)lOp/c~n I = O(h~). Similarly, at a shock wave of finite intensity v lOu./Onl = O(z/h, ), and at a contact discontinuity v l,~u./Onl = O ( v ) . Therefore, let us introduce the features: a3m =

(hi/fi,j)'[ 6~p/~n [ij,

(47)

a4,, = r IOu,/~n 10,

(48)

keeping in mind that we are going to use in the following a classification technique based on the minimum distance principle. In eqn (47):

fi~=(ll4)(pt+t./+~ +pi ~.j+t +P,+Lj ~+P,-I.j-~)

(49)

and the quantity z in eqn (48) is the time step of a difference scheme used for the numerical integration of the Euler equations (1)-(2). The averaging eqn (49) was chosen rather arbitrarily, that is one can in principle also use other formulae for P0' The quantity h~lOp/Onl~j was computed with the aid of the formulae of the Sobel edge detector:

hll ap/On[ 0 = (1/8) {[Hi x F 3(p)]. (hi/h2)sin ~0 + [n2

X

F 3(p)]cos ctO},

(50)

where F3(p) is a 3 x 3 image window centred in the pixel p,~, the masks H~, H2 are determined by eqn (7). Finally, the feature eqn (48) was computed with the aid of eqns (36)-(37) or (36) and (38). Note that the quantity eqn (47) may also be interpreted as the approximation of the quantity

h~13 inp/On[~. In order to distinguish between the points of discontinuity surfaces and the points belonging to the subdomains of continuous flow (for example, to compression waves), let us introduce the binary feature:

IAt--h~l/Am,x>6 ~ O, ]A~-A~i/Am~,~6~. 1,

a~,,=

(51)

Pattern recognition methods

53

The quantities At and Ar in eqn (51) are calculated in accordance with eqns (29)-(31), Amax = max(At, At), 6~ is a user-specified positive constant. This constant is set with regarding the following considerations. First, IA t - A,I = O(1) at the shock waves and contact discontinuities of finite intensity; in the domains of continuous flow IA t - Ar I = O(An) where An is a distance in the ( x , y ) plane between the points with indices (i~, j~) and (ir,J~); in the case of a square mesh An ~< 2x/~h I . In order to distinguish between the points of purely contact discontinuities and the points of tangential discontinuities, we make use of the binary feature: 1, ( r / h l ) l u , t - - u ~ l > 6 2 O, (v/hl)lu~l-u~rl<<.62"

a6m=

(52)

In eqn (52) 62 is a user-specified positive constant: u¢t = vi~,j~ cos ~u - ui,.h sin O~ij

uCr = v,,j, cos ~ij - ui,,jr sin ~0'

(53)

In order to distinguish between the points of normal shock waves and the points of oblique shock waves, let us introduce the binary feature:

a7~=

1,

( z / h ~ ) ( l u , t l + lU~rl) >63,

0,

('c/h,)(lu=zl+lu~rl)<~63,

(54)

where 63 is a user-specified positive constant. To substantiate the use of this feature in the case of a normal shock wave it is necessary to prove that the magnitude of the tangential component of the gas velocity will be small in the zone of smearing of such a discontinuity. In the practical computation of the quantities u,t and u~ using eqn (53) we have used the grid values of the velocity components u, v in the nodes lying on a straight line corresponding to one of the eight principal directions (these directions for the case of a curvilinear grid are shown in Fig. 4). Let us number these grid values of the quantities u and v as u~ . . . . . UN; V~ . . . . . VN where the subscript "1" refers to the one end of the zone of a smeared normal shock wave, and the subscript " N " refers to the other end of the above zone. Assume that the sequences u~ . . . . . UN and v~ . . . . . VN are monotone: sign(uj-uj+l)=Sl;

sign(vj-vj+l)=s2,

j=l

.....

N-I

(55)

where sk = 1 or - 1 , k = 1, 2. The relationship eqn (55) is valid while using monotone difference schemes. In the case of using many well-known schemes of second and higher orders of approximation, one is usually able to indicate such a neighbourhood of a shock wave front in which the relationship (55) is satisfied. Outside this neighbourhood (usually behind the front of a smeared shock wave) the second-order schemes generate parasitic oscillations, see for example Ref. [1]. Thus, let us assume that the relationship (55) is satisfied. Assume also that the equations: u~l = 0,

u~u = 0,

(56)

take place at the ends of a zone of smearing of a normal shock wave, that is at j = 1 and j = N where: u~j = vj cos ~t - uj sin ~, j = 1. . . . . N,

(57)

~t is the angle between a normal to the wave front and the positive direction of the x-axis. Let us show that the only condition for u,j which is noncontradictory with the requirements eqns (55) and (56) is vanishing of the quantity u,j throughout the zone of smearing of a normal shock wave. Since uj, vj in eqn (55) may have arbitrary signs, it is convenient for the following to employ the formula: i u,i = ~ (u,j -- u~j_ l), j=2

i = 2 .....

(58)

N.

This formula is obtained with regard for the first relationship in eqn (56). Substituting the right-hand side of the formula eqn (57) instead of uj in eqn (58), we obtain: u~ = cos e"

( v j - vj_ j) - sin e . J

(uj - uj_ i) , j=2

i = 2 . . . . . N.

(59)

54

E.V. VOROZHTSOV

Consider now the case:

t!l-t ) ~>0,

uj-uj

~>0

Vj,

cosc¢>0,

sinct<0.

Then we obtain from eqn (59): U~N> 0, which contradicts the second of the relationships eqn (56). Assume now that: ej-vj

~>0,

uj-uj_~>O

Vj,

cos~>0,

sin~>0.

The orientation of the front line of a considered shock wave with respect to the axes Ox', Oy" of some other coordinate system x ' 0 y ' whose axes 0x', 0y' are rotated by a certain angle with respect to the axes of the original coordinate system xOy will change in the numerical computation of the same gas-dynamical problem in the x ' 0 y ' system. Let us rotate the axes 0x and 0y by such an angle that in the new coordinate system x ' 0 y ' the inequalities cos ~ > 0, sin ~ < 0 are satisfied. Then we arrive at a case already considered. All other possible combinations of the signs of the quantities s~, s2 and cos ~, sin ~ are considered in a similar way. In all of these cases the only requirement for u~j which does not contradict the conditions eqns (55) and (56) is the condition u~j = 0 Vj. Thus, we have introduced for the characteristic of the objects (Xm, y,,) in the feature space 7 features at . . . . . . . aT,,,. A nondimensional character of the chosen functional dependencies for aim gives reason to expect that the results of singularities classification with the aid of these features will comparatively weakly depend on the specifics of a gas-dynamical problem under study. In the following we shall use, as in Ref. [4], the features bl . . . . . . . bTm obtained from a~. . . . . . a7,, with the aid of a linear transformation:

bkm=ak,,,

k=1,2,5,6,7;

bkm=wkkak,~,

k=3,4

(60)

where Wkk = l/maxm ak,,, k = 3, 4. The transformation eqn (60) enables one to reduce the intraclass distances. The form of the transformation eqn (60) takes into account the smallness of the dispersions of the features a~m, azm, asm, a6m, a7m [4]. Note that since the shock wave speed D does not enter explicitly in the expressions for the features akin, k = 1. . . . . 7, these features may be used for the classification of singularities both in stationary and in nonstationary gas flows. 5. A L G O R I T H M S

OF

PATTERN

CLASSIFICATION

5.1. Minimum-distance classifier Pattern classification on the basis of the minimum-distance criterion is one of the first ideas of the automatic pattern recognition. As was pointed out in Ref. [28], this simple classification method proves to be a very efficient tool in the solution of the problems in which the classes are characterized by a degree of variability being within reasonable limits. For a realization of this method it is at first necessary to construct a training set of patterns. In our case we shall assume that there do not exist any other classes of objects besides the classes ~o~. . . . . ~o7 which were discussed in the preceding section. Thus, we assume that we know in advance a set of classes to which the points (Xm,Ym) under study may be assigned. This enables us to realize the idea of supervised classification, which is associated with the use of a notion of sample patterns [28, 31]. Taking into account the properties of the objects belonging to each of the classes ~ol . . . . . ~7, as well as the properties of the features b~. . . . . . b7m constructed in Section 4 it is possible to find the samples in each of the classes eo~. . . . . ~o7, see Appendix. The knowledge of the samples enables one to realize the parallel classification scheme [28] in which the assignment of the object (xm, y,,) to some class is carried out independently of the fact to which classes the other objects were assigned, therefore, the classification of all the objects (xm, y,,) can be carried out in this case with the aid of the minimum-distance principle simultaneously, which can be especially efficiently realized on a computer with parallel processors. From the point of view of the optimality of the classification systems [72] the scheme of parallel classification is not optimal, since it requires the classification to be a calculation of all the coordinates b~ . . . . . b 7 of the pattern vectors, which is not economical in terms of computer expenses. In this connection let us turn to the idea of hierarchical classification in which the data subject for classification is taken from the only class which may be subdivided into subclasses and subsubclasses [73]. In our case such a subdivision was discussed in Section 4 in the analysis of the

Pattern recognition methods

[Minimum /

55

-distance

classifier

Q

if/'4 class Shock waves and compression waves

J

~ s class Contact discontinuities rarefaction waves and other subdomains of continuous flow

Q

Fig. 8. The scheme of hierarchical classification of singularities in gas flows.

interrelationships between the classes ~o~. . . . . co7. In Fig. 8 we present one of the possible schemes of hierarchical classification. While implementing this scheme one feature blm is not sufficient for the assignment of the object (xm, y,,) to one of the two classes Q4 or Qs, because the equality bt,. = 2 can be satisfied both at the shock waves and at the contact discontinuities. Similarly, the equality b2,, = 2 can be satisfied both at the shock waves and in some subdomains of the continuous flow. In this connection we have used at the initial stage of the scheme of Fig. 8 a minimum-distance classifier in which four features blm ~ b2m, b3m, b4m are involved. Let = (bl, b2, b3, b4), ~m :

(bin, b2,n, b3m, b4m), ¢1 ~" (ZII, ZI2, ZI3, "714), ~(2k) = (Z[kl), ~22'-(k)"23"r(k)"24"~(k)~1'k = 1,2.

(61)

In eqn (61) (~ is a sample vector for the (xm, Ym) points belonging to the class Q4. We place in the class Qs, as in Ref. [4], two samples ~zr°l,~2r(z)having the coordinates: Z~II) : 0,

"~22"(I)= z'23"(1): Z24(1)_- 0," zt2) = 2,

z'22"(2): Zt 2) = ~24"¢2)_- v.n

(62)

Further: zlt = 2,

zl2 = 2,

z13 :

Z I 4 "~"

1.

(63)

A linear discriminant function d~ (/~) may then be written in the form [4, 28]: d, (/J) = prE, - (1/2)(r~, where the superscript T denotes the transposition operation. Since two samples enter in Qs, we compute the discriminant function d2(~) by the formula: d2(~) = max { p r ~ _ (1/2)((~0)r~o}, t

1 = 1, 2

(64)

in accordance with the formulas for the case of multiple samples presented in Ref. [28]. The pattern p is assigned to the class Q4, if d~(p)> d2(~). After the original objects ( x , , , y , , ) , m = 1 . . . . . N 4 , have been subdivided into two classes Q4 and Q5 with the aid of the above minimum-distance classifier, the set Q4 is considered. If it is not empty, the shock waves are separated from the

56

E.V. VOROZHTSOV


yes

iqO

®

®

yes

no

yes /

I

no

® rio

"~ I yes Q ~_ no

®

yes

}0 no

\ yes

Fig. 9. Decision tree in the sequential classification scheme.

compression waves with the aid of the feature b 5. After that the class f~5 is considered. If it is not empty, a subset of potential points of contact discontinuities is extracted with the aid of the feature bs. Then a subsubclass of contact discontinuity points is extracted with the aid of the feature b2.

5.2. Sequential classification Figure 9 presents a decision tree for the solution of the same classification problem as in the case of the classification scheme of Fig. 8, but using a different sequence of features. Here the idea of sequential classification [73] has been realized. We have denoted by 098 in Fig. 9 a class of rarefaction shock waves which in principle should not arise in the case of using a difference scheme ensuring the increase in the entropy at the shock waves. Nevertheless, in the case of substantial oscillations in the numerical solution the class 098 can prove to be non-empty upon completing the process of classification of the set (Xm, Ym). We shall call for convenience the classification schemes presented in Figs 8 and 9 "scheme I" and "scheme II", respectively. Following Ref. [72], denoted by cj is the cost of the measurement of the j t h feature bjm, j = 1. . . . . 7. In our case one can mean cj the computing time needed for the calculation of the j t h feature. In the decision tree presented in Fig. 9 are both interior and terminal nodes. The less the number of interior nodes in a tree, the shorter the way in the graph at the recognition of ruth object (xm,y,~). In this connection in Ref. [72] the questions of optimizations of the decision trees were considered both from the point of view of minimization of the cost of objects classification and from the point of view of minimization of the number of interior nodes. In Figs 8 and 9 the interior nodes are numbered by encircled figures. In scheme I there are five interior nodes and in scheme II eight such nodes. However, it is needed in scheme II at each stage to compute one feature. In Ref. [74] there are formal proofs of the optimality of such classification procedures, from the point of view of measurement costs (in our case of computer time needed for the computation) of the features. Thus, from the point of view of the "cost" the scheme I is worse than the scheme II. However, one should keep in mind, along with the factor of cost, the reliability of classification schemes. In the case of scheme I a decision on the assignment of the object (xm, ym) to one of the two classes f~4 or ~ is made on the basis of ensemble of the four features bl, b2, b3, b4, so that the random

Pattern recognition methods

57

errors in one or two features may not affect the classification result. In the case of scheme II the object (Xm,y,,) will be classified incorrectly in the case of an error in one feature. Thus, the classification reliability makes increased demands for the features used in the systems of sequential classification. In particular, the features should be less sensitive to the noise than in the case of a classification system based on the minimum-distance principle. Aiming at the satisfaction of this requirement we have applied the smoothing operators while calculating the majority of the features, see the eqns (28)-(30), (46), (47), (49)-(51). In the algorithm of sequential classification the scheme of which is presented in Fig. 9 we have used the formulae for computation of the features/76,.,/77,, which are different from the features a6m, a7m proposed in Section 4:

/76m=

1, O,

1, /77"= 0,

I(u~i- u~,)/q;jl >66 I(u~Z--U~r)/q;jl ~<66; (1/2) (I U,[ + I U~rI)/q~i > 67

(1/2)(lU~;l+lu~rl)/qa<~67.

(65) (66)

In eqns (65) and (66) q,j = (u~ + v~)°5, 06, 07 are user-specified positive constants, 0 < 06, 07 < 1. The quantities u~;/q;j, u,,/qo in the case of a purely contact discontinuity have a clear physical interpretation: they indicate a relative contribution of the tangential velocity component to the total magnitude of the fluid velocity. The quantities 06 and 07 may be set in the interval ]0, 1[ rather arbitrarily. The user should keep in mind that the computer will classify the contact discontinuities as purely contact or tangential discontinuities, and the shock waves as normal or oblique shocks in the sense of the definitions eqns (65)-(66). Note that at the desire of the user it is easy to realize in the scheme I a subdivision of a class of contact discontinuities into the subclasses of purely contact discontinuities and of tangential discontinuities and a subdivision of a class of shock waves into the subclasses of normal and oblique shock waves. For this purpose it is sufficient to implement in the scheme of Fig. 8 the features a6" and a7m [eqns (52), (54)] or the features/76m,/77m [eqns (65), (66)] in the same way as the features /76mand/77" have been implemented in the scheme of Fig. 9. The above algorithm of sequential classification requires also the specification of positive constants 03 and 04, see Fig. 9. The choice of these constants should be carried out in accordance with a general procedure for determining the values of features for the sample patterns, see for example Ref. [28]. In accordance with this procedure a relatively well studied problem was at first chosen while considering some class of gas-dynamical problems, such that for this problem there is an information on the types of singularities which may be present in a given problem. After that such sample values of the constants 03 and 64 were determined by test computations at which the level of classification errors was minimal. As a rule, for each of the constants 03 and 64 it is possible to indicate certain interval within which the classification result depends very weakly on the specific values 03 and 64 from the corresponding intervals. It appears that owing to this fact it proves to be rather easy to determine the values fi3 and 64 when going over to an analysis of other classes of aerohydromechanics problems. In computations of problems where the intensity of shock waves varies quickly with time it proves to be useful to set 03 and 64 as some functions of time (for example, as linear functions). Let us now describe a scheme of sequential classification for the extraction of shock waves in the computations of stationary potential transonic flows. The conditions at the shock waves in such flows are obtained as a particular case from the Rankine-Hugoniot conditions [69, 70] at D = 0: /91Unl ~---p2Un2;

(67)

u~1 = u~2.

(68)

In addition, it follows from the Zemplen's theorem that at D = 0:

l u,21/c2> 1, l u~,l/c, < I.

(69)

The components u, v of the velocity vector in a 2-D potential flow were computed by eqn (19). The orientation angle ~tU of the vector of a normal to the edge at the node (x~j,y,j) was computed by eqn (24). After that the normal velocity component u, and the tangential velocity component u~

58

E.V. VOROZXTSOV

at same node (xkz, Ykl) being a neighbour of the node (x~j,Yo) were determined by the formulae: (u,,)k~ = Ukt COS C~j+ Vktsin c%,

(u~)kt = v~t cos ~,j -- ukt sin ~ . Let (it, Jr) and (it, L) be the indices of the grid nodes (x o, yg) lying on the different sides of the edge and closest to the line of a normal to the edge (see Section 4). The subscripts "1" and "2" in eqn (69) refer to the quantities behind and before the shock wave front, respectively. Therefore, to check the inequalities eqn (69) it is at first necessary to determine which of the points (xt,,j~, Yi,.j,) and (x,,.j,,yt,.j,) is located behind the front of a shock wave in the potential flow. In the case of nonisentropic shock waves we have checked for this purpose the inequality At > At, see eqns (29)-(31). This criterion is inapplicable in the case of isentropic shock waves. In the literature several analogues of the entropy condition at a compression shock wave in a potential flow are known. In particular, it was shown in Ref. [75] that if the pressure function p ( V ) , V = l / p , satisfies the Weyl condition, namely, the function p ( V ) is convex, then at a discontinuity across which the pressure and the density increase (shock wave) the relationship: Pl + Pl U2nl --

(P2 "~ p2U2n2)~ 0

(70)

is satisfied where the subscripts "1" and "2" refer to the gas state behind and before the shock wave front, respectively. For an approximate determination of the quantities p~, p~, u,~, P2, Pz, u,2 entering eqn (70) we have realized a search interpolation procedure similar to the one presented in Section 4, see the formulae (28)-(30). This procedure enabled us to find the values Pt, Pz, u,~ and p,, p,, u,,. Setting then in eqn (70)p, = Pt, P, = P~, u,,~ = u,,t, P2 = P,, P2 = Pr, Un2 Unr, we have checked the satisfaction of the inequality eqn (70). If it is not satisfied, the pairs of indices (it, j~) and (it, Jr) were interchanged. To increase the reliability of classification we have used along with the condition (70) also a simpler criterion enabling us to exclude from further consideration the rarefaction shock waves and having the form [76]: =

u,,.[IUI] < 0,

u,2.[IUI] < 0,

(71)

where the symbol [I U I] denotes the jump of the magnitude of the velocity IUI across a shock wave. The conditions (71) were implemented in the classification computer code as follows. At first the values Pt, P~, u.~,pr, p~, u.~ were found with the aid of a search interpolation procedure (the corresponding points with the indices (il, Jr) and (i~, j~) were interchanged, if it was needed), such that Pt + PlU2nl - P, - P, U2,~ >t O. After that the inequalities: Utlr'[(U2nt'Ji - U2rl) -- (U2nr "3I- U~r)] < O;

u,,. [(uZ,,+ u2~t)- (u]~ + U~r)] < 0

(72)

were checked in accordance with eqn (71). If one considers a flow around an airfoil the chord of which is oriented along the x-axis and the flow moves in the direction of the x-axis, it can be assumed that u,, > 0, u,~ > 0 near the shock wave front, and then the formulae (72) are reduced to a trivial inequality IUtl z < l U r Iz which corresponds to a well-known fact that the flow slows down behind the shock wave front. Figure 10 shows the scheme of sequential classification that we have realized for the extraction of shock waves in 2-D stationary potential transonic flows. Note that in virtue of the fact that the mathematical model of potential flow is simpler than the model of a nonisentropic rotational flow the scheme of Fig. 10 is simpler than the scheme of Fig. 9. The features al,,, a2,, were computed with regard for eqns (67) and (68) by the formulae: aim = I[Pi.,y, "(u.)i..j. -- Pi,.j,' (u.)#.j, ]/p,j"

(u.)ql,

az,, = I [( u~ ),. j~ - (u.),,, j, ]/( u,)tjt .

6~,62 are user-specified positive constants. In order to determine the sample values of these constants we have carried out a number of computations of the transonic flows around an airfoil by the AF2 scheme for different values of the free stream Mach number M= for which the fact

Pattern recognition methods

59

/

IUnll /C1< I and lunzl /c2
IN°

Continuous ftow

lunll/cl> I and lu.2J /c2>I Y~[

INo

Interchange the pairs of indices

Pt + PL u2nL~> Pr + /Or Unr

(iL , j~), (ir, Jr)

U°r(l~l z- IZTrl2)
I Yes lun21/c2>t

t

N0

I Yes al m

azm

< 81

No

I Yes < ~2

No

I Yes Output of t h e coordinates of shock waves points

J

Fig. 10. The scheme of sequential classification of singularities in the stationary transonic potential flows.

of the presence of shock waves is established in the literature. As a result we have found a range of the variation of the quantities a~m, a~,, at the points (xij,yu) belonging to shock waves: 0.0014 ~< a~m ~< 0.013; 0.0012 ~ a2m ~ 0.0049. Taking this into account we have set 6~ = 0.02, 62 = 0.01. Note that the choice of the specific values 6~, 62 generally depends on a difference scheme employed, because the shock waves are smeared over different number of mesh intervals in the numerical solutions obtained by different difference schemes. After the completion of the work of the classification system presented in Fig. 10 one can use the data obtained on the presence or on the absence of shock waves in a problem under study in an expert system for the automation of the design of aerodynamic bodies which meet some set of requirements, in particular, a requirement of the absence of shock waves in a transonic flow around an airfoil. 6. E X A M P L E S O F P A T T E R N S C L A S S I F I C A T O N IN T H E N U M E R I C A L T W O - D I M E N S I O N A L GAS D Y N A M I C S P R O B L E M S

SOLUTIONS

OF

6. 1. Tests of the pattern classification algorithms Assume that Ns points (O<~Ns<~N4) from the set of N4 points (xm,ym) introduced into the system of singularities classification along with the p, u, v, p images should belong to shock waves, and N~ points should belong to contact discontinuities (0 ~< N~ ~< N4), so that N~ + N, <~N4. Let

60

E.V. VOROZHTSOV

N~ Points be assigned to the shock waves in the result of classification, where 0 ~< N~ ~< N4 and it can generally occur that N~ > N s. Similarly, let N~ be a number of points, assigned to the contact discontinuities by the classification system, 0 ~< N~ ~< N4, N~ + N~ ~< N4. Then the levels of the errors of classification of shock waves 3N~ and of contact discontinuities 6N~ may be determined quantitatively as:

6N~ = (IN~ - N,I/N~.). 100%, fiN, = (IN~ - N,I/N,.)" 100%.

(73)

It is of interest to obtain analytic estimates of the quantities eqn (73) while using a specific classification scheme. Having such estimates, one would be able to compare the reliability of different classification schemes. However, it was pointed out in Ref. [31] that the exact calculation of the level of classification error is very difficult even in cases when the probabilistic structure of the classification problem is known completely. In this connection it was proposed in Ref. [31] to apply an empirical approach which enables one to avoid the above difficulties and which consists in experimental tests of a classifier. It was shown in Ref. [31] at the examples that the reliability of the estimation of the true level of classification errors by means of eqn (73) depends substantially on the numbers Ns and N,. of the objects introduced in the classifier: the more are the values N~ and N,., the closer are the values (fN~ and fiN,. to their true values. In our case the quantities N~ and N,. depend substantially on the size of the steps of a spatial computing mesh. Indeed, the finer is this mesh, the more number of cells of the mesh will be traversed by a given discontinuity line. For example, if the shock wave front is oriented along the x-axis, then it is obvious that

N~ = O(l/h, ). For the estimation of the level of classification errors by the eqn (73) it is necessary to have a control set of the objects (Xm, Ym) for which the true values of N~ and N¢ are known. Here one can use two approaches. The first one reduces to obtaining of the above set (Xm, Ym) in the result of segmentation of a digital image p(x, y, t) obtained by the solution of a specific gas-dynamical problem by some shock-capturing difference scheme. Another approach consists in the generation of some artificial numerical solution of some hypothetical 2-D gas-dynamical problem. Note that such method of testing is widespread in the theory of digital pattern recognition, see for example, [31, 38, 41, 72, 73]. The first of the above approaches to testing of classifiers appears to be more natural, since the classification schemes presented in Section 5 are designed just for the classification of singularities on the basis of shock-capturing numerical solutions. However, it should be noted that each specific difference schemes has (i) its own intensity of discontinuities smearing (or of the edge blur, in terms of the theory of the digital image processing); (ii) its level of parasitic oscillations, which depends in addition for the same scheme on the value of the Courant number K, on the magnitude of the solution gradients in the domains adhering the discontinuity fronts, on the relation h2/h~ and on some other factors [1, 77]. In this connection it is of interest to carry out the tests of classifiers within the framework of the second of the above approaches. This subsection is just devoted to a brief presentation of the results of such testing carried out in Ref. [4]. The next subsections of this section may be regarded as the presentation of the results of testing of the classifiers of Section 5 within the framework of the first approach, that is on the basis of actual finite-difference shock-capturing solutions of a number of applied gas-dynamical problems. When constructing the test in Ref. [4] we aimed at the satisfaction of the following requirements: (1) both shock waves and contact discontinuities should be present in the "gas flow"; (2) the noise and the blur with controlled characteristics should be present in the numerical solution; (3) the shock waves intensity should be variable, changing from some finite magnitude to small values corresponding to a continuous flow; (4) the angles of orientation of the discontinuity lines with respect to the x-axis should vary within a wide range, that is the discontinuities fronts should be curvilinear. When generating the p, u, v, p digital images which meet the above four requirements we introduced an additive noise with controlled signal-to-noise ratio (we shall denote it by S / N in the following) in each of these four images by a formula taken from [38]:

S / N = d2/~2~,

(74)

where a,, is a mean quadratic deviation of an independent Gaussian noise, d is the height of a step

Pattern recognition methods

61

edge. Since in our case the profiles of the quantities p, u, v, p differ substantially from the step edge profile, there arise certain difficulties in determining the "height" of the edges of these quantities. In practical computations the height d ( f ) where f i s any of the functions p, u, v, p was computed by the formula [4]:

12dxdy) °5

75,

c_?I

where ~, is a domain in the (x, y) plane in whichf(x, y) ~f2, f2 is the value of the function f ( x , y) in an undisturbed medium before the shock wave front, sl is the area of the domain ~1. It is easy to see that in a particular case of a step edge the formula (75) yields the exact value of edge height d. The square spatial domain in the (x, y) plane was discretized in Ref. [4] by a 40 × 40 mesh, that is: this was a comparatively rough mesh. The quantities 6N~, 6N~ in eqn (73) were calculated on this mesh in Ref. [4] at different S/N ratios. In the result of application of the classification schemes I and II from Section 5 the following rough estimates were obtained in Ref. [4] at S/N = oo (that is there is no noise): &Ns = fiN, = 0 for scheme I; fiNs = 20%, fiN, = 0 for scheme II. When the S/N ratio diminishes, the characteristics fiNs and &No for both classification schemes become closer. For example, at S/N = 10 the following data may be obtained from the results of Ref. [4]: fiNs = 50%, fiN, = 20% for scheme I and fN~ = 45%, 6N, = 20% for scheme II. It should be noted that the case S/N = 10 corresponds to a highly oscillatory numerical solution; the oscillations were superposed in Ref. [4] not only in the domains behind the fronts of discontinuities, but also in the remaining spatial subdomains, that is also within the limits of the zones of discontinuities smearing and in the domains of an undisturbed flow. The properties of clustering of the patterns of points corresponding to (xm, Ym) points in a 4-D Euclidean space of the coordinates b~, b2, b3, b4 were also studied in Ref. [4]. At S/N = ~ there are in this space three distinct clusters of the (btrn, b2m , b3m , b4m ) points. As the S/N ratio decreases, the cluster diameters increase; in addition, the boundary between the clusters of points corresponding to shock waves and contact discontinuities becomes smeared, fuzzy. Nevertheless, even at S/N = 10 the presence of clusters is obvious. This fact shows that the formulae of Section 4 for the computation of the features are good from the point of view of the requirement of good clustering in the feature space. The main conclusion which may be drawn from the results of Ref. [4] is that in the case of using the oscillatory difference schemes for which the quantity eqns (74) has finite value of the order 10 ~< S/N <<.100, the both of the schemes I and II from Section 5 yield comparable levels of the classification errors. The second conclusion is that the most favourable situation for the classification of singularities with the aid of the schemes I and II is the situation when a monotone difference scheme is used in the computations of a gas-dynamical problem and in the numerical solutions obtained by this scheme all the discontinuities are smeared over not more than 3 - 4 cells of a spatial mesh.

6.2. Recognition of shock waves in transonic flow around airfoil In Fig. 1 l(a) we present the results of the application of the scheme of sequential classification of Fig. 10 to the set of the points (xm, Ym) obtained as a result of the image segmentation and shown in Fig. 5(c). We have used at a stage of debugging of the classification program instead of checking the inequality eqn (70) a simpler criterion taking into account the flow direction over the airfoil the chord of which is parallel with the x-axis. Since this flow moves in the direction of the x-axis, the inequality x ~2)< x t') should be satisfied for the points lying before and behind the shock wave front, respectively. Taking this into account we assumed that x c2)= xi~.j, if the inequality xi,,j, < x~,.j, was satisfied. Otherwise the pairs of indices (i~,jj) and (ir,jr) were interchanged. After that we checked the analogues eqn (72) of the entropy inequality. The result of the work of such simplified classification program (it does not need a search interpolation procedure described in Section 4) is presented in Fig. l l(b). Comparing Figs l l(a) and (b), we can see that the weakening of the classification criteria leads to an increase in the number of points included by the classification program in the class of the shock waves points. One can get some increase in the number of shock waves points while using the criterion eqn (70) by performing a more accurate computation of the

62

E.V. VOROZHTSOV (a)

(b)

Fig. 11. The application of the sequential classification scheme of Fig. 10 in the problem on potential flow around the NACA 0012 airfoil. ( ) Sonic lines and the airfoil contour; (***) the points of shock waves.

potential flow. For example, to obtain a stationary solution of the AF2 type scheme we made by analogy with [53] a number of iterations providing a three-order-of-magnitude drop in the residual as compared to the residual value for the initial iteration. It appears that the error in the quantities p, p, u, describing the states on both sides of a weak shock wave in a potential flow may be slightly reduced by performing a larger number of iterations. In Ref. [78] the coordinate x~ of the shock position on the airfoil surface, measured from the leading edge, was defined as the middle point of the j u m p in the pressure coefficient curve. The arrows in Fig. 11 indicate the shock positions on the airfoil surface, that were found by this method for the same problem in Ref. [78] on the basis of computational results of Ref. [56] (see Fig. 2b). If only the very fact of shock waves presence near the airfoil surface is of practical interest and there is no need in an accurate information on the extent of the shock front lines, then it is sufficient (for example, for the expert systems of automated design) to have an information on the shock waves of the form presented in Fig. 1 l(a). There arises a question whether it is possible to use in the gaps between the points (Xm,y,,) detected in the result of classification the points of sonic lines to fill the above gaps. A study carried out in Ref. [79] with the aid of matched asymptotic expansions for the case of viscous transonic flows shows that the sonic line position does not depend in the first approximation on small dissipative effects and corresponds to the shock front.

6.3. Cylindrical shock problem Let us now present three more examples of gas-dynamical flows for the automatic analysis of which we have used the pattern recognition algorithm described in the present chapter. A common feature of these three examples is that the solution gradients in the neighbourhood of discontinuities are nonzero; in addition, there are in the second and third examples presented below the lines of strong discontinuities along which the jumps of gas-dynamical parameters vary within rather wide limits, so that there may be present in the same problem, for example, both weak and intensive shock waves.

Pattern recognition methods

63

In order to estimate quantitatively the accuracy of localisation of a shock wave and a contact discontinuity by the image segmentation procedure presented in Section 3 it was necessary to find exact values of the abscissas of the points of strong discontinuities for certain chosen moments of time. For this purpose we have considered the problem within the framework of a mathematical model of a 1-D axisymmetric flow the equations of which may be written in the form [80, 81]: (76)

w, + [F(w)], = 7t (w, r)

where M F=

W ~

R-

_ (7 -- 1)

E

M2

I°

r

-2-R

7E+

2

=

,

\R.]J

(77)

0

R = rp,

M = rpu,

E = rp(E + u2/2) -- re,

p = [(~ -- 1)/r] [E - M 2 / ( 2 R ) ] = (~ - 1)(e - m2/(2p)),

(78)

m = pu.

In eqns (76)-(78) r is a spatial coordinate, the symmetry axis corresponds to the value r = 0, ~ is a constant in the equation of state, eqn (3). We shall determine the "exact" solution of the system (76)-(78) by analogy with [81] with the aid of a numerical integration of this system with very small steps h and r along the r- and t-axes, respectively. At the points of a spatial grid for which r / > h let us approximate the system (76) by the M a c C o r m a c k scheme [39]: r~7 +l = w ~ - (r/h)(F~+l - F ~ ) + r ~ 7 ; W n + l = ( 1 / 2 ) [ ( W n "71- I~ n + I) - - ( r / h ) ( F

n

+ l _ p n 4-11) + r l p n + 1].

(79)

Let us also use by analogy with eqns (4-5) the smoothing: n+ I w- ni+ I = [1/(c.p -+-2)](C,0wn+l "k ~"i+, -f- WT-+,'),

(80)

where if7 +~ is a vector of the smoothed solution, ~o is an empirically chosen parameter. In the example of computation by scheme (79), (80) presented below we have used the value cp = 200. In order to avoid spurious fluxes of the material across the boundaries of a computational domain on the r-axis while using the smoothing eqn (80), we have applied the formula (80) in a slightly modified form: n+l n+l if,"+' = w7 +' + [I/(~p + 2)](Q~+ ,:2 - Qi-,/2),

(81)

/'in + I

where ~:i+~/2r~"+t= ,+, ,j", ,+~~+-~w ~ . We have assumed the quantity ~+~,,2 to be equal to zero at the boundaries of the computational domain, that is Q " ~ .'- = ~O"+~m~x~+'/2= 0 where we have taken rm,x = im,x'h = 1.5 (the subscript " i " in eqns (79)-(81) refers to the cell boundary, so that ri = ih). It follows from eqn (78) that on the symmetry axis where r = 0, i = 0 it is necessary for the computation o f p and p to use the formula different from eqn (79). At the same time, the boundary conditions for the components of the vector w in eqn (77) have a simple form [81]: ug=O,

Rg=[rpg],=o=O,

M g = Rgug = O, Eg = [rpg/(7 - 1)]'=° + M g u g / 2 = 0. C.A F 18 I--E

(82)

64

E . V . VOROZHTSOV

Stable extra boundary formulae in the case of using the M a c C o r m a c k scheme in the interior nodes are simple extrapolation formulae [82]

p~= 2RT/rl - R~/r2;

eg = 2ET/rl - E~/r2;

p~ = (7 - l)e~.

(83)

The computational domain 0 ~< r ~< 1.5 was covered by a uniform mesh of 300 nodes. For the time advancement of the solution from t = 0 to t = 0.3 458 time steps were needed at such a mesh. Denote by r~,.(t) and r,(t) the abscissa of a shock wave front and of a contact discontinuity, respectively. The values of r,,,, r, were found on the basis of a numerical solution by scheme (78), (81)-(83) as the abscissas of the nodes at which the gradients Ic~p/c~r[ reached their local maxima in the zones of smearing of a shock wave and a contact discontinuity. The initial data for the system eqns (76)-(78) were set as in [83]:

p{,={4,

r<~r,,(O).

l,

r > r,,,(0)'

p0=

{5, 1,

r<~r~.,(0) r > r,.~,(0)

(84)

where r~,.(0) = 0.4; 7 = 2 in the equation of state eqn (3). As a result of integration of the system (76)-(78) by the method (79), (81)-(83) at an initial distribution eqn (84) the following values of the quantities r~,,, r, were obtained: r~,(0.2) = 0.770;

r,(0.2) = 0.510;

r,,(0.3) = 0.940;

r,(0.3) = 0.570.

(85)

The 2-D computations of the problem on the axisymmetrical breakdown of discontinuity were carried out as in Ref. [83] by the modified particle-in-cell method of Harlow (for the calculation of the pressure with the aid of the equation of state (3) the density p was used that was obtained by a recalculation by the F L I C type scheme approximating the continuity equation). A uniform spatial 26 x 26 mesh with the steps h~ = h2 = 0.05 was used in these computations. Sixteen particles were placed in each cell lying to the left of a discontinuity at t = 0, and 4 particles were placed in each cell lying to the right of a discontinuity at t = 0 taking into account the initial condition eqn (84). In Fig. 12 we present the results of using a hierarchical classification technique the scheme of which has been presented above in Fig. 8. It is to be noted that the scheme of Fig. 8 has a shortcoming which manifests itself at t---0.3, see Fig. 12(b): this scheme detected spurious rectilinear contact discontinuities near the x- and y-axes which intersect at a right angle. We explain this phenomenon by the effect of using the reflection boundary condition at the boundaries x --- 0 and y = 0 in the computations by Harlow's method, which in this case leads to large density gradients on the parts of the above lines adhering the coordinate origin (0, 0), see Fig. 13. As may be seen in Fig. 12, the classification scheme 1 not only correctly classifies the cylindrical discontinuity surfaces, but also gives with a rather good accuracy (with an error within the limits of the size of the steps h~ = h2 = 0.05) the position of these discontinuities.

(a)

(b)

X2

X2

1.0.

0,5,

"L0--

0.5

~

0.0 0.0

X~

0.5

.

I

1.0

Fig. 12. T h e c y l i n d r i c a l s h o c k p r o b l e m .

X1

0.0 0.0

0.5

1.0

(a) t = 4 0 r = 0.2; (b) t = 6 0 r = 0.3.

Pattern recognition m e t h o d s

65

2.8-

2.1

"1.4

0.7

X2 12 Fig. 13. Isometrics of the surface p = p(x, y, t) at t = 0.3.

6.4. Double Mach reflection of strong shock wave The flow under study can be set up experimentally, for example, in an impingement of a shock wave with a planar front on the wedge. In the result of an interaction of this wave with the wedge surface there occurs a diffraction of a shock wave and a system of secondary shock waves and contact discontinuities is formed. Depending on the values of some basic parameters (the initial angle of the incident shock wave, the shock wave Mach number, the adiabat constant 7, the atomic composition of a gas in which the shock propagates), various types of shock-wave configurations are realized. A survey of some results of experimental studies of this problem may be found in Refs [84-86]. The numerical studies were carried out by a number of the authors, see for example [87, 88]. We have used in the computations the results of which are presented below the same input parameters as the authors of Ref. [88]: 7 = 1.4 in the equation of state eqn (3), the shock Mach number in air M0 = 10, the initial angle of the shock wave incidence on the wedge surface is 60°; the density and the pressure of the undisturbed air ahead of the shock are 1.4 and 1, respectively. In the interior nodes of a uniform 100 x 30 mesh the computation was carried out by the MacCormack scheme [39] augmented by the smoothing process (4), (5). The value ~0 = 50 found empirically was used in eqn (5). The initial and boundary conditions in the computations by the scheme of Ref. [39] were set as in Ref. [88]. The input parameters in a problem under study correspond to a double Mach reflection of the shock wave from a wall. In the process of such a reflection two Mach shocks are formed together with two contact discontinuities. The second contact discontinuity is extremely weak, therefore, it is impossible to detect it in the computations on a relatively crude mesh of 100 x 30 cells. As was pointed out in Ref. [88], it is generally possible to observe in the 480 × 120 mesh computations by the P P M D E scheme developed by the authors of Ref. [88] the second contact discontinuity by examining the picture of the lines of constant values of the velocity components u and v. The second Mach shock is rather weak, and its amplitude is reduced to zero by the time it reaches the contact discontinuity emanating from the first triple point. The results of a computation of an irregular reflection presented in Fig. 14 have been obtained at Courant number K = 0.65. One can observe in Figs 14(b) and (d) the effect of the roll-up of a tangential discontinuity emanating from the main triple point. One can also see in these figures a fragment of a compression wave which in accordance with known information on the structure of the flow under consideration should be indeed a weak Mach wave emanating from the second triple point. The dashed lines in Figs 14(b), (d) show the positions of discontinuities, which we have found from the pictures of density isolines presented in Ref. [88] and obtained by the P P M D E method [88, 89] on a mesh with the steps h~ = h 2 = 1/120.

66

E.V. VOROZHTSOV

(o) 1.0 X2

0.5"

i

Xl

0.0

.......

O0

J

0.5

.

.

.

.

I

1.0

-,

- it ~

1.5

"

~

I

2.0

.

.

.

.

I

2.5

.

.

.

.

3.0

(b) 'I.0 X2

I

0.5-

0.0, 0,0

1.0

1.5

2.0

fl.o

1.5

2.o

2.5

3.0

(c) 1.o

0.5

o.o

I o.o

05

. . . .

I

. . . .

2.5

3.0

(d) I.o

+

0.5

0.0 0.0

0.5

1.0

-1.5

2.0

2.5

3.0

Fig. 14. The results of singularities classification by different schemes: (a), (b) the scheme ! (***--shock waves: F-lEE]--contact discontinuities; ~ O - - - - - c o m p r e s s i o n waves); (c), (d) the scheme II (× x x - normal shock waves; ***--oblique shock waves; F1F1Fl--tangential discontinuities; + + + - - p u r e l y contact discontinuities; O O O---compression waves). (a), (c) t = 0.1, the number of time level n = 75: (b), (d) t =0.2, n = 160.

Pattern recognition methods

67

The P P M D E method, as well as the P P M L R method which will be mentioned in Section 6.5, represent a higher-order extension of Godunov's method of a type first introduced by van Leer in his M U S C L algorithm [90]. The sequential classification scheme presented in Fig. 9 enables one to obtain more complete information on the properties of the flow singularities in comparison with the scheme of Fig. 8, since the scheme of Fig. 9 includes a subdivision of the shock waves class into the subclasses of the normal and oblique shock waves and of the contact discontinuities class into the subclasses of tangential and purely contact discontinuities. The results of sequential classfication presented in Figs 14(c) and (d) have been obtained at the following values of the constants 61, 63, 64, 6 6, 67 for the features eqns (47), (48), (51) and (60), eqns (65) and (66), respectively: 6~ = 0.01, 6 3 = 0. l , 6 4 = 0. l , 66 = 0.05, 67 = 0.05. The threshold value T ~') for the image segmentation was determined by eqn (9). The set of the objects (x,,, Ym) which are subject to classification and are obtained as a result of the work of the image segmentation procedure is here the same as in the case of Figs 14(a), (b). It is interesting to note that the tangential discontinuity emanating from the main triple point goes over into a purely contact discontinuity when it rotates towards the flow, but near the reflecting wall it again becomes a tangential discontinuity. According to Ref. [88], in a region limited by the second Mach wave, the curved reflected wave and the reflecting wall there is a very little vertical motion. As may be seen in Fig. 14(d), the main Mach stem proved to be a normal shock wave. From this a conclusion may be drawn that the vertical gas motion is very little also behind the front of this shock wave. Under the same conditions the incident shock wave which initially was a normal shock is at t = 0.2 also classified as a normal shock wave. 6.5. Supersonic flow in a wind tunnel with a lower-wall step This 2-D problem is often used as a test for finite-difference schemes [88, 90, 91]. In the computational examples presented below we have used the same input parameters as the authors of Ref. [88]: V = 1.4 in the equation of state (3), the Mach number of a uniform flow at the channel entrance M~ = 3; the density and the pressure of the undisturbed air are 1.4 and 1, respectively. In the interior nodes of a uniform 90 x 30 mesh the computation was carried out by the MacCormack scheme [39] augmented by the smoothing process eqns (4) and (5). The value ~0 = 30 was used in eqn (5). The Courant number K = 0.6. The initial and boundary conditions in the computations by the scheme of Ref. [39] were realized as in Ref. [88]. As was shown in Refs. [90], a steady flow in the problem under study develops by the time 12. We shall be interested in the flow for earlier times when the flow structure rapidly changes with time. The classification results presented in Fig. 15 are obtained with the aid of the hierarchical classification method the scheme of which has been presented above in Fig. 8. The dashed lines in Figs 15 and 17 show the positions of discontinuities, which we have found from the pictures of density isolines presented in Ref. [88] and obtained by the P P M L R method [88, 89] on a 240 × 80 mesh. In Fig. 16 the dimetrics of the surface p = p(x, y, z) are shown. It may be seen in this figure that the numerical solution gradients are nonzero in the domains behind the shock waves fronts. In addition, it is seen in Fig. 16 that the intensity of each shock wave changes significantly along the line of the wave front. In Fig. 17 we present some results of the application of the sequential classification algorithm the scheme of which has been given in Fig. 9. The most complicated flow structure develops by the time t = 3, when in the computational domain there are along with the detached shock wave also a Mach stem, the shock wave emanating from the triple point; a secondary shock wave resulting from the reflection of the foregoing shock from the lower wall; and, finally, a shock wave emanating from the upper wall and being the result of a reflection from this wall of a secondary shock just mentioned. In addition, one can see in Figs 15(c) and 17 a contact discontinuity beginning near the triple point. A numerical boundary layer along the upper surface of the step causes a spurious Mach reflection from the step, but the corresponding Mach stem is only two zones high, as this may be seen in Figs 15(c) and 17. This numerical effect arising in the case of using the MacCormack scheme was noted previously in Ref. [88]. One can also observe in Figs 15(c) and 17 the presence of a spurious contact discontinuity near the upper wall of the wind tunnel. This spurious contact discontinuity may be caused by two reasons: first, by the use of simple reflection conditions on the rigid wall in the MacCormack scheme computations. Second, the same conditions are used in the computations in the Sobel edge detector

68

E.V.

VOROZHTSOV

(a)

1.0 ,~

..... .... ::: :::x

x:xL

:

X2

0.5

X1

0.0

,

.

0.0

.

.

0.5

.

.

.

.

.

.

1.0

I

1.5

. . . .

I

2.0

I

. . . .

3.0

2.5

(b) 1.o X2

0.5

X1 0.0 O0

0.5

1.0

1.5

2.0

2.5

3.0

(c)

I0 X2

~

,~

0.5

I

X~

0.0. 0,0

0.5

'1.0

1.5

2,0

F i g . 15. (a) t = 0.5, n = 112; ( b ) t = 1.5, n = 3 2 5 ; (c) t = 3.0, n = 6 3 2 . *** S h o c k

2.5 waves;

3.0 D[5][SI c o n t a c t

discontinuities,

at the stage of segmentation of the digital image p (x, y, t). One can also see in Figs 15(c) and 17 the gaps in some lines of secondary shock waves. We explain these gaps by the fact that the contact strip containing a contact discontinuity which emanates from the triple point intersects with a shock wave emanating from the triple point and with the next secondary shock wave. In the region of intersection of this contact strip with the zone of smearing of an oblique shock wave the errors in the entropy may become so large that the sign of the entropy production may prove to be incorrect. Such points are assigned to the continuous flow in the classification algorithms used. And, finally, one can see in Fig. 15(c) a spurious contact discontinuity starting near the corner of the step and reaching the triple point. Judging by Fig. 16, this contact discontinuity is located in the domain of the rarefaction fan, near the bottom of a "valley" in the p(x, y, t) surface. The change of the entropy in the numerical solution along the detected spurious contact discontinuity obviously proved to be sufficiently large to assign the line under consideration to the class of discontinuity lines. Comparing Figs 15(c) and 17, one can note that the line of the above spurious contact discontinuity is substantially shorter while using the sequential classification scheme so that the level of classification errors in the case of using this scheme proves to be less than in the case

Pattern recognition m e t h o d s

69

6-

5432IO-

o

I

x

~,,I0.5

2

,.0

3 Fig. 16. D i m e t r i c s of the surface p = p(x, y, t) at t = 3.0.

of the hierarchical scheme presented in Fig. 8. It is to be noted that the reliable classification of contact discontinuities on the basis of shock-capturing numerical computations proves to be a more difficult matter than the classification of shock waves. We relate this to two reasons. First, the contact discontinuities are smeared in a numerical computation more intensively than the shock waves [1, 77]. Second, a part of algebraic equations representing the jump conditions in inviscid gas flows degenerate at the contact discontinuities in the sense that these equations are satisfied identically. As a result of this the number of basic relationships which could be used as the features of contact discontinuities is diminished. In conclusion of this section let us present some data on the computer codes that implement the algorithms presented in this chapter for automatic singularities classification in the 2-D inviscid gas flows. If one takes the length (in the number of machine words) of the program for solving the 2-D problem by the MacCormack scheme as one unit, then the length of the programs realizing the classification schemes of Figs 8-9 is 2.57 units. About 40% of this volume are constituted by the subroutines of graphical information output (plotting isolines, surfaces, the points of discontinuity lines). The length of the program for classification of singularities in the 2-D transonic potential flows (together with graphical subroutines) is approximately equal to the length of the program for computing the transonic potential flow around an airfoil by the AF2 type scheme on a given curvilinear mesh. To economize the core memory of a computer, the magnitudes eqn (6) of the gradients gu and of the orientation angles of the gradient vector eqn (8) needed at the stages of the image segmentation and of the singularities classification were not stored in the form of 2-D numerical arrays, as this was done for example in Refs [2-5], but were computed each time by the ~.o [

|

""

lllllllllllll ll,|l|llll,II

',/

0.5

o

. 0.0

0.5

0 10

~ 1.5

2.0

2.5

3.0

Fig. 17. The result of classification by scheme 1I at t = 3.0; × x x n o r m a l shock waves; *** oblique shock waves: []l-ql~ t a n g e n t i a l discontinuities; + + + purely contact discontinuities.

70

E.V. VOROZHTSOV

grid values p,~ of the density with the aid of corresponding subroutines. A comparison with the variant when gij and ~0 are computed once and stored in arrays shows an increase by a factor of about 2 in the computer time needed for segmentation and classification. The operations of the computation of the quantities go and ~,j by the formulae eqns (6) and (8) are carried out independently for each pixel (i, j), therefore, the computer time needed for the segmentation can be substantially reduced by using the cellular logic devices [92, 93] as well as (to a lesser extent) by using the parallel general-purpose computers. Let us now present some data on the computer time which was needed for the classification of singularities in various gas-dynamical flows and for the graphical display of the data in percents with respect to the time needed for the solution of a gas-dynamical problem (the computations were carried out on a sequential general-purpose computer): transonic potential flow--12.5%; double Mach reflection of a strong shock wave--10.3%; supersonic flow in a wind tunnel with a lower-wall step--2.7%. Although there are here significant percent differences, the absolute values (in the seconds of computer time) of the times needed for the classification of singularities in the above three problems were close to each other. This is related to the fact that the computer time needed for the digital processing of 2-D arrays is determined mainly by the size (that is by the number of pixels) of the images. In all the three problems listed above we have used the number of spatial computing mesh nodes having the same order of magnitude. 7. C O N C L U S I O N S

.Methods for localisation and automatic classification of singularities presented above are more laborious and complicated than any of the methods for shock localisation on fixed and dynamic grids in our book [1]. A doubltess advantage of the presented methods (see also Ref. [5]) is their universality and applicability for an automatic analysis of the numerical solutions of the problems for which there is no a priori information on the singularities. It should be noted that the above presented methods for localisation and classification of singularities in 2-D flows can be extended for the case of 3-D flows. Indeed, the algorithms of segmentation of 3-D digital images have already been considered in the literature, see for example Ref. [94]. As regards the computation of the features needed for the algorithms of singularities classification, we cannot imagine here any substantial difficulties of the extension of corresponding formulae for the case of three space variables, because all the formulas for the features are written down in Section 4 in the coordinates along a normal and along a tangential direction to the discontinuity surface. The variety of existing methods for the digital image processing and pattern recognition creates a huge field for the research of the applicability of these methods in the systems of localisation and automatic classification of the singularities in the numerical solutions of multidimensional gas dynamics problems. We shall not enumerate these methods here (a review of the most widespread of these methods has already been made above, see also Refs [2, 5]), but we indicate a number of problems which in our opinion can be successfully solved with their aid. First, we would like to point out the application of the methods for the restoration of digital images to improve the accuracy of shock-capturing numerical solutions. Such an improvement of the p, u, v,p images should increase the reliability of localisation and classification of singularities (see also Ref. [5]). The use of the methods of automatic analysis of flow structure in the expert systems for automation of the design of aerodynamic bodies satisfying some optimality requirements appears also to be very promising. One can expect here a further increase in the computational efficiency of such expert systems by an integration of the general-purpose computers and the computers based on the cellular logic principle.

REFERENCES 1. E. V. Vorozhtsov and N. N. Yanenko, Methods for Localisation oJ" Singularities in the Numerical Solutions of Gas Dynamics Problems. Nauka, Novosibirsk (1985). 2. E. V. Vorozhtsov, On shock localisation in difference solutions by the Sobel edge detector. Preprint Inst. Theor. Appl. Mech. Sib. Div. USSR Acad. Sci., Novosibirsk, N 12 (1985). 3. E. V. Vorozhtsov, Extraction of discontinuity lines in difference solutions by pattern recognition methods. Proc. Semin.

Pattern recognition methods

71

of Socialist Countries on the Computational Aerodynamics, Complex Problem "Cybernetics", p. 30. Scientific Council of the U.S.S.R. Academy of Sciences, Moscow (1985). 4. E. V. Vorozhtsov, Classification of discontinuities in gas flows as the pattern recognition problem. Preprint Inst. Theor. Appl. Mech. Sib. Div. USSR Acad. Sci., Novosibirsk, N 23 (1986). 5. E. V. Vorozhtsov, On shock localisation by digital image processing techniques. Comput. Fluids 15, 13 (1987). 6. E. V. Vorozhtsov, On automatic classification of singularities in the numerical simulations of two-dimensional gas flows. Image Processing and Remote Sensing. Abstracts of the Reports of the Regional Conf. Novosibirsk, p. 63. Comput. Centre Sib. Div. USSR Acad. Sci., Novosibirsk (1987). 7. E. V. Vorozhtsov and S. I. Mazurik, Application of artificial intelligence techniques for the classfication of discontinuities in gases and for the stability analysis of difference schemes. Soviet Union-Japan Symposium on Computational Fluid Dynamics, U.S.S.R., Khabarovsk, Book of abstracts (Edited by P. I. Chushkin), p. 124. Far Eastern Div. U.S.S.R. Acad. Sci., Vladivostok (1988). 8. O. M. Belotserkovskii, Mathematical modeling is the branch of informatics. Cybernetics. Formation of Informatics, p. 45. Nauka, Moscow (1986). 9. V. V. Rusanov, Processing and analysis of computation results for multidimensional problems of aerohydrodynamics. Lecture Notes Phys. 18, 154 (1973). 10. A. A. Dorodnitsyn, lnformatics: The Subject and Problems, N2, p. 85. Vestnik AN SSSR (1985). 11. N. J. Nilson, Principles of Artificial Intelligence. Tioga Publishing Company, Palo Alto, CA (1980). 12. P. G. Buning and J. Steger, Graphics and flow visualization in computational fluid dynamics. A1AA 7th Comput. Fluid Dyn. Conf., Cincinnati, Ohio (1985); Collect. Techn. Pap., p. 162, New York (1985). 13. P. Kutler, J. L. Steger and F. R. Bailey, Status of computational fluid dynamics in the United States. AIAA Paper 87-1135 (1987). 14. R. G. Belie, Some advances in digital flow visualization. AIAA Paper 87-1179 (1987). 15. P. Kutler, A perspective of computational fluid dynamics. Proc. Fourth Int. Conf. Boundary and Interior Layers-Computational and Asymptotic Methods, Novosibirsk, U.S.S.R. (Edited by S. K. Godunov, J. J. H. Miller and V. A. Novikov), p. 332. Boole Press, Dublin (1986). 16. P. Kutler, U. B. Mehta and A. Andrews, Potential application of artificial intelligence concepts to numerical aerodynamic simulation. Lecture Notes Phys. 218, 340 (1985). 17. J. F. Dannenhoffer III and J. R. Baron, Robust grid adaptation for complex transonic flows. AIAA Paper 86-0495 (1986). 18. J. F. Dannenhoffer III and J. R. Baron, A hybrid expert system for complex CFD problems. AIAA Paper 87-111 I, p. 99. Proc. A1AA 8th Comput. Fluid. Dyn. Conf., Honolulu (1987). 19. A. E. Andrews, Progress and challenges in the application of artificial intelligence to computational fluid dynamics. AIAA J. 26, 40 (1988). 20. I. L. Dobroserdov, Base models and methods for the construction of rapid algorithms and programs for computing non-isobaric viscous jet flows with reacting components. Manuscript dep. VINITI 4.05. 1987. N. 4143-V87. VINITI AN SSSR, Moscow (1987). 21. S. V. Bobyshev and I. L. Dobroserdov, Identification modeling of processes in the non-isobaric region of a turbulent jet. Modeling in Mechanics (in Russian) 1(18), No. 6, p. 3. Novosibirsk (1987), 22. S. S. Tong, Design of aerodynamic bodies using artificial intelligence/expert system technique. AIAA Paper 85-0112 (1985). 23. S. S. Tong, Coupling artificial intelligence and numerical computation for engineering design. AIAA Paper 86-0242 (1986). 24. G. S. Pospelov, Artificial intelligence. New information technology. Cybernetics. Formation of Informatics, p. 106. Nauka, Moscow (1986). 25. V. M. Kovenya, Some problems of computational fluid dynamics. Construction of Numerical Algorithms and the Solution of Mathematical Physics Problems (Edited by K. I. Babenko), p. 5. Institute of Applied Mathematics of the U.S.S.R, Academy of Sciences, Moscow (1987). 26. Z. B. Golembo, V. P. Zinkevich, Mathematical problems of the image analysis and recognition. Itogi Nauki i Tekhniki. Tekhni~. Kibernetica, Vol. 18, p. 123. VINITI AN SSSR, Moscow (1985). 27. V. A. Kovalevskii, Optimal Solution Methods in Image Recognition. Nauka, Moscow (1976). 28. J. T. Tou and R. C. Gonzalez, Pattern Recognition Principles. Addison-Wesley, Reading (1974). 29. V. I. Timohin, Application of the Computers to the Solution of Pattern Recognition Problems. Leningrad State University, Leningrad (1983). 30. L. Goldfarb, A unified approach to pattern recognition. Pattern Recognition 17, 575 (1984). 31. R. O. Duda and P. E. Hart, Pattern Classi~eation and Scene Analysis. Wiley, New York (1973). 32. Ya. A. Fomin and G. R. Tarlovskii, Statistical Theory of Pattern Recognition. Radio i Sviaz, Moscow (1986). 33. K. S. Fu, Syntactic Pattern Recognition and Application. Prentice-Hall, Englewood Cliffs, NJ (1982). 34. L. Hesselink, J. Helman, Evaluations of flow topology from numerical data. AIAA Paper 87-1181 (1987). 35. Yu. N. Vatolin, Information criterion of stability. ~isl. Metody Mekh. Splo~noi Sredy, Novosibirsk 2, 22 (1971). 36. Yu. N. Vatolin, On the application of the entropy estimates of the stability. Cisl. Metody Mekh. Splo.~noi Sredy, Novosibirsk 5, 5 (1974). 37. A. A. Petrov, Algorithmical software for information and control systems of adaptive robots (algorithms of technical robot vision). Itogi Nauki i Tehniki. Technig. Kibernetika, Vol. 17, p. 251. VINITI AN SSSR, Moscow (1984). 38. W. K. Pratt, Digital Image Processing. Wiley-lnterscience, New York (1978). 39. R. W. MacCormack, The effect of viscosity in hypervelocity impact cratering. AIAA Paper 69-354 (1969). 40. W. G. Habashi, Numerical methods of turbomachinery. Recent Advances in Numerical Methods in Fluids, Vol. 1, p. 245. (Edited by C. Taylor and K. Morgan). Pineridge, Swansea (1980). 41. A. Rosenfeld and A. C. Kak, Digital Picture Processing. Academic Press, New York (1976). 42. J. F. Thompson, Z. U. A. Warsi and C. W. Mastin, Numerical Grid Generation. Foundations and Applications. North-Holland, New York (1985). 43. L. R. Miranda, Application of computational aerodynamics to airplane design. J. Aircraft 21, 355 (1984).

72

E.V. VOROZHTSOV

44. J. Flores, J. Barton, T. Hoist and T. Pulliam, Comparison of the full-potential and Euler formulations for computing transonic airfoil flows. Lect. Notes Phys. 218, 213 (1985). 45. J. M. Longo, W. Schmidt and A. Jameson, Viscous transonic airfoil flow simulation. Z. Flugwiss. Weltraumforsch. 7, 47 (1983). 46. D. L. Whitfield and J. L. Thomas, Transonic viscous-inviscid interaction using Euler and inverse boundary-layer equations. Computer Methods Viscous Flows, p. 451. Pineridge Press, Swansea (1984). 47. T. L. Hoist, Solution of the transonic full potential equation in conservative form using an implicit algorithm. Numerical Methods for the Computation of lnviscid Transonic Flows with Shock Waves. A G A M M Workshop (Edited by A. Rizzi and H. Viviand), p. 28. Vieweg, Braunschweig (1981). 48. A. Rizzi, Computational mesh for transonic airfoils. Numerical Methods for the Computation of lnviscid Transonic Flows with Shock Waves. A G A M M Workshop (Edited by A. Rizzi and H. Viviand), p. 222. Vieweg, Braunschweig (1981). 49. T. L. Hoist, Implicit algorithm for the conservative transonic full-potential equation using an arbitrary mesh. AIAA J. 17, 1038 (1979). 50. J. Flores, T. L. Hoist, D. Kwak and D. M. Batiste, A new consistent spatial differencing scheme for the transonic full-potential equation. AIAA J. 22, 1027 (1984). 51. T. L. Hoist, Approximate-factorization schemes for solving the transonic full-potential equation. Recent Advances in Numerical Methods in Fluids Series, Advances in Computational Transonics (Edited by W. G. Habashi), Vol. 4, p. 59. Pineridge, Swansea (1985). 52. J. C. South Jr and M. M. Hafez, Stability analysis of intermediate boundary conditions in approximate factorization schemes. AIAA Paper 83-1898-CP (1983). 53. J. Flores, T. L. Holst and R. L. Sorenson, Transonic solutions for a multielement airfoil using the full-potential equations. J. Aircraft 22, 50 (1985). 54. Y. S. Wong and M. M. Hafez, Some iterative schemes for transonic potential flows. AIAA J. 23, 808 (1985). 55. R. C. Lock, Test cases for numerical methods in two dimensional transonic flow. AGARD Report R-575-70 (1970). 56. A. Jameson, D. A. Caughey, W. H. Jou and J. Steinhoff, Accelerated finite-volume calculations of transonic potential flows. Numerical Methods .[or the Computation of lnviscid Transonic Flows with Shock Waves. A G A M M Workshop (Edited by A. Rizzi and H. Viviand), p. 11. Vieweg, Braunschweig (1981). 57. M. Hafez, W. Whitlow Jr and S. Osher, Improved finite-difference schemes for transonic potential flow calculations. AIAA J. 25, 1456 (1987). 58. M. Ya. Ivanov and V. V. Koretskii, Computation of flows in two-and three-dimensional nozzles by the method of approximate factorization. Zh. vychisl. Mat. mat. Fiz. 25, 1365 (1985). 59. M. Ya. lvanov, V. V. Koretskii, A. S. Lieberson and D. B. Solovyov, A high-accuracy conservative approximate factorization scheme for the integration of the full equation for the velocity potential in a transonic range. Theoretical fundamentals and the construction of numerical algorithms for solving mathematical physics problems. Vlth All-Union School, Gor'ky, Abstracts of the Reports, p. 66. Gor'ky State University, Gor'ky (1986). 60. P. M. Byval'tsev and M. Ya. Ivanov, A rapid method for the computation of potential transonic flows in two- and three-dimensional cascades of blades. Theoretical fundamentals and the construction of numerical algorithms for solving mathematical physics problems. Vllth All-Union Seminar, Kemerovo, Abstracts of the Reports, p. 20. Kemerovo State University, Kemerovo (1988). 61. A. P. Aralov, Yu. B. Lifschitz and A. A. Shagaev, Solution of multi-dimensional problems of potential transonic flows. Theoretical fundamentals and the construction of numerical algorithms for solving mathematical physics problems. Vllth All-Union Seminar, Kemerovo, Abstracts of the Reports, p. 7. Kemerovo State University, Kemerovo (1988). 62. V. E. Makarov, The use of the approximate factorization method for the computation of plane and spatial sub- and transonic potential flows. Theoretical fundamentals and construction of numerical algorithms for solving mathematical physics problems. Vlth All-Union School, Gor'ky, Abstracts of the Reports, p. 102. Gor'ky State University, Gor'ky (1986). 63. A. N. Kraiko, Some questions of the construction of the numerical algorithms for computation of ideal gas flows. Construction of Numerical Algorithms and the Solution of Mathematical Physics Problems (Edited by K. I. Babenko), p. 33. Institute of Applied Mathematics of the U.S.S.R. Academy of Sciences, Moscow (1987). 64. O. N. Gus'kov and V. E. Makarov. Transonic potential cascade flow and cascade design, and simulation of a separated incompressible flow in a cascade of blades. Theoretical fundamentals and the construction of numerical algorithms for solving mathematical physics problems. Vllth All-Union Seminar, Kemerovo, Abstracts of the Reports, p. 38. Kemerovo State University, Kemerovo (1988). 65. B. N. Petrov, G. M. Ulanov and S. V. Ulyanov, Information content of the features and the compression of information process of the control. Itogi Nauki i Tehniki. Techni(eskaya Kibernetika, Vol. 13, p. 3. V1NITI AN SSSR, Moscow (1980). 66. S. K. Pal and B. Chakraborty, Intraclass and interclass ambiguities (fuzziness) in feature evaluation. Pattern Recognition Lett. 2, 275 (1984). 67. A. L. Gorelik, I. B. Gurevich and V. A. Skripkin, State-of-the Art of the Recognition Problems, Some Aspects. Radio i Sviaz, Moscow (1985). 68. G. 1. Razorenov and G. A. Poddubskii, Automated selection of the features for objects classification. Zav. Lab. 51, 48 (1985). 69. B. L. Ro£destvenskii and N. N. Janenko, Systems of quasilinear equations and their applications to gas dynamics. Translations of Mathematical Monographs, Vol. 55. Am. Mathematic. Society, Providence, Rhode Island (1983). 70. B. L. Ro~destvenskii, Mathematical theory of shock waves. Mathematical Encyclopedia, Vol. 5, p. 468. Soviet Encyclopedia, Moscow (1985). 71. L. MSr6 and Z. Vassy, A simplified and fast version of the Hueckel operator for finding optimal edges in pictures. Proc. 4th Int. Conf. Artific. Intelligence, Tbilisi, U.S.S.R.p. 650 (1975). 72. G. R. Dattatreya and L. N. Kanal, Decision trees in pattern recognition. Progress in Pattern Recognition, Vol. 2 (Edited by L. N. Kanal and A. Rosenfeld), p. ',89. Elsevier, Amsterdam (1985). 73. E. B. Hunt, Artificial Intelligence. Academic Press, New York (1975). 74. G. R. Dattatreya and V. V. S. Sarma, Bayesian and decision tree approaches for pattern recognition including feature measurement costs. IEEE Trans. Pattern Analysis and Machine Intelligence PAMI-3, 293 (1981).

Pattern recognition methods

73

75. J.-J. Chattot, Condition d'unicite pour les 6coulements irrotationnels stationnaires de fluid parfait compressible. C.r. Acad. Sci. Paris 2,86, A l l l (1978). 76. P. Morice and H. Viviand, Equations de conservation et condition d'irreversibilit~ pour les 6coulements transsoniques potentiels. C.r. Acad. Sci. Paris 291, B235 (1980). 77. P. J. Roache, Computational Fluid Dynamics. Hermosa, Albuquerque, New Mexico (1976). 78. A. Rizzi, H. Viviand, Collective comparison of the solutions to the workshop problems. Numerical Methods for the Computation of lnviscid Transonic Flows with Shock Waves. A G A M M Workshop (Edited by A. Rizzi and H. Viviand), p. 167. Vieweg, Braunschweig (1981). 79. J. Mace and T. C. Adamson, Shock waves in transonic channel flows at moderate Reynolds number. A I A A Paper 0369 (1985). 80. R. D. Richtmyer and K. W. Morton, Difference Methods for Initial-Value Problems, 2nd edn. Wiley Interscience, New York (1967). 81. S. Arbarbanel and M. Goldberg, Numerical solution of quasiconservative hyperbolic systems--the cylindrical shock problem. J. Comput. Phys. 10, 1 (1972). 82. Yu. I. Shokin and N. N. Yanenko. Method of Differential Approximation, Application to Gas Dynamics. Nauka, Novosibirsk (1985). 83. V. M. Fomin, E. V. Vorozhtsov and N. N. Yanenko, On the properties of curvilinear shock waves "smearing" in calculations by the particle-in-cell method. Comput. Fluids 7, 109 (1979). 84. G. V. Ba~enova, L. G. Gvozdeva, Yu. P. Lagutov, V. N. Liahov, Yu. M. Faresov and V. P. Fokeev, Nonstationary Interactions o f Shock Waves and Detonation Waves in Gases. Nauka, Moscow (1986). 85. M. Shirouzu and I. I. Glass, Evaluation of assumptions and criteria in pseudostationary oblique shock-wave reflections. Proc. R. Soc. London A406, 75 (1986). 86. I. I. Glass, Some aspects of shock-wave research. AIAA J. 25, 214 (1987). 87. Yu. M. Lipnitskii, V. N. Liahov, Numerical solution of a problem on the diffraction of a shock wave on the wedge. Izv. A N SSSR, Mehanika ~,idkosti i Gaza Vol. 6, p. 88 (1974). 88. P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks. J. Comput. Phys. 54, 115 (1984). 89. P. Colella, P. R. Woodward, The piecewise parabolic method (PPM) for gas-dynamical simulations. J. Comput. Phys. 54, 174 (1984). 90. B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to G o d u n o v ' s method. J. Comput. Phys. 32, 101 (1979). 91. J, P. Vila, Simplified G o d u n o v schemes for 2 × 2 systems of conservation laws. S l A M J. numer. Anal. 23, 1173 (1986). 92. K. Preston, Jr., M. J. Duff, S. Levialdi, P. E. Norgren and J.-I. Toriwaki, Basics of cellular logic with some applications in medical image processing. Proc. IEEE 67, 826 (1979). 93. A. Rosenfeld, Algorithms for cellular image processing. ICCD'83. Proc. IEEE Int. Conf. Comp. Des.: VLSI Comput., Port Chester, New York, p. 719. Silver Spring, M D (1983). 94. D. C. Morgenthaler and A. Rosenfeld, Multidimensional edge detection by hypersurface fitting. IEEE Trans. Pattern Anal. Mach. Intell. PAMI-3, 482 (1981).

APPENDIX

D e r i v a t i o n o f the F o r m u l a s f o r C o o r d i n a t e s o f S a m p l e s Let us denote the coordinates of sample patterns for each of the classes 09~, k = 1. . . . . 7, by sk~ . . . . . SkT. The question of determining the coordinates of samples in the feature space plays an important role for a successful application of the minimum-distance principle from the theory of pattern recognition [28]. Consider at first the question on determining the values ofskt for the classes 09~, 092, 093 and 094. As was shown in Section 4, at,, = 2 at the shock waves, thus one can take for the classes 09~ and 09., the sample value sk~ = 2, k = 1, 2. As is known, three possibilities can realize at the contact discontinuities: Ou,/On < O, Ou,/On >10. The situation Ou,/On = O, that is when the normal velocity component is constant in the neighbourhood of a contact discontinuity, is encountered very seldom in the actual applied problems of gas dynamics. Nevertheless, such a situation can occur in the computations of model problems on the convective transport of a fluid mass by a constant velocity field. Suppose that this velocity field is user-specified by the formulae: u=U0cosfl,

v=U0sinfl,

(AI)

where U0 = const > 0 and fl is real. The use of the machine floating-point arithmetic leads in practice to the fact that even in the case when the velocity field is given by eqn (AI) the equality Ou,/On = 0, where the Ou,/On derivative is approximated by finite differences, does not take place [4]. In this connection we have taken into account only two sample values of the coordinate s I for the classes 093 and 094: skt = 2 and skt= 0, k = 3, 4. Consider now the question on determining the values of ski for the classes co~, 092,093 and 094,J = 2, 3, 4. As was already discussed in Section 4, sk: = 2 for k = 1, 2 and s~2 = 0 for k = 3, 4. In a geneal case the sample values sk3 and Sk4 can be set as functions of time. In this paper we have used a simpler technique for determination of s~3 and sk4 which was based on a passage to a limit in the formulae (47)-(48) as h L--*, h2--* 0, z --*0. Indeed, let us consider the formulae (47) and (48) in the case of a contact discontinuity. From the properties of the functions p ( x , y , t ) and u,(x,y,t) at the contact discontinuities it follows that: lim as, , = 0, hi ~ 0

lira a4,, = 0.

(A2)

T~ 0

We obtain from eqn (A2) that at sufficiently small steps h~, h., and ~ it is possible to use the sample values sk3 = 0, sk4 = 0, k = 3,4, for the classes of contact discontinuities 093 and 094- At a shock wave of finite intensity

74

E.V. VOROZHTSOV

a~,,, = O(1). a4,,,= O(r/h~ ). With regard for eqn (60) we have that at the shock waves b30, = aa.,/max a3., = O( 1)/O( 1) ~ 1, m

b4m = a4m/max a4m = O ( z i h i ) / O ( z / h l ) = O(I) ~< I. m

Thus, if the intensities of shock waves in a problem under study vary insignificantly, then the values of b3., and b4m are close to 1, as this follows from eqn (A3). In the case of shock waves with an intensity which varies strongly in space and time the scatter of the values of b3., and b4m is more significant. Nevertheless it can be assumed also in this case that there will take place in the (b 3, b4) plane a clustering of the (b3, ., b4,.) points around the point with coordinates sk3 = 1, sk4 = 1, k = 1. 2. at sufficiently small h t, h2 and r. This clustering property of the (b3. ,, b4.,) points was corroborated in Ref. [4] by the tests on problems in which the shock wave intensity changed from zero to a finite value of the order 0(l). From the above procedure for obtaining the values of sk3, sA4,k = 1.2.3, 4 one can draw a conclusion that the minimum-distance classifier which makes use of these constant values ofsk3 and s~4 (1 ~