EP2941726A1

EP2941726A1 - Method for simulating thermal radiation between surfaces

Info

Publication number: EP2941726A1
Application number: EP14701297.5A
Authority: EP
Inventors: Jakob FAINBERG
Original assignee: MAGMA Giessereitechnologie GmbH
Current assignee: MAGMA Giessereitechnologie GmbH
Priority date: 2013-01-07
Filing date: 2014-01-07
Publication date: 2015-11-11
Also published as: KR102177554B1; CN104956370B; KR20150104577A; CN104956370A; US10970429B2; WO2014106670A1; US20160048614A1

Abstract

The invention relates to a method for simulating the influence of thermally coupled surface radiation on a solid body, which solid body has at least one surface that can be irradiated, by calculating the radiative transfer between grey, diffusive surfaces, characterised in that the surface or surfaces to be irradiated is/are divided adaptively and hierarchically into radiation tiles of equal or virtually equal radiation intensity, and the surface temperature brought about by irradiation is produced by means of a hierarchical view factor method. Said view factor method comprises the evaluation of a solid angle integral using a primary solid angle division. Said primary solid angle division comprises a homogeneous view factor discretisation, wherein that solid angle division is discretised adaptively and hierarchically into the part-regions thereof by spherical projection and wherein the sum of all the partial amounts of that solid angle integral can be determined by means of ray tracing.

Description

METHOD FOR SIMULATING WARMING RADIATION BETWEEN SURFACES

DESCRIPTION

The invention relates to a method and an associated algorithm for the calculation of the thermal surface radiation, and their application for simulating the influence of the thermally coupled surface radiation, with particular reference to

Casting processes.

Further, the invention will be described in the context of other helpful methods and algorithms, which are described herein as parts of the method and algorithm for calculating surface thermal radiation, and their application for simulating the influence of thermally coupled surface radiation

However, applications are not limited to the local method and algorithm.

The invention therefore also relates to a method for discretizing a solid angle for use in a simulation or computation process for savings in computer time and consumed

To achieve computer memory.

Further, the invention therefore relates to a method of ray tracing for use in a simulation or computation process to achieve acceleration and savings in computer time and used computer memory. Especially is an unconventional method using Parallel Computing is used to accelerate a Ray Tracing calculation alone and in its combination with Anisotropic Chebychev Distance calculations and / or additional Acceleration proposed by Tile Association.

TABLE OF CONTENTS Background

Description of the pictures

principle

Energy balance at the radiating surface

Calculation of the incident radiation flux

radiation tile

Discretization of the solid angle

Hierarchical system

Adaptation to the geometry and the temperature distribution

Ray tracing

acceleration

Ray tracing

Tile Association

Anisotropic Chebychev distance

Parallel ray tracing

Thermal calculation

applications

BACKGROUND

Many important properties, such as the viscosity, the electrical conductivity, the specific volume or the chemical reactivity of a material or Media are determined by the temperature. For example, in the field of melting, refining and / or the shaping processing of solid materials, such as metals, glass, ceramics, compliance with certain temperatures, often time-dependent, is of great importance.

Casting processes. The invention further relates to a method for simulating the influence of the thermally coupled surface radiation by calculating the radiative exchange between gray, diffuse surfaces, characterized in that the surface or surfaces to be irradiated is adaptively, hierarchically subdivided into radiation tiles of the same or nearly the same radiation intensity and by irradiation to calculate the resulting surface temperature by means of ray tracing as the sum of all partial amounts of the radiation tiles.

The invention further relates to a method as mentioned above wherein the solid angle is adaptively and hierarchically discretized into its partial areas by spherical projection. The invention further relates to a method as mentioned above wherein the ray tracing procedure is accelerated. Further, the invention will be described in the context of other helpful methods and algorithms, which are described herein as parts of the method and algorithm for calculating surface thermal radiation, and their application for simulating the influence of thermally coupled surface radiation

However, applications are not limited to the local method and algorithm. Furthermore, the invention relates to a method for simulating the influence of thermally coupled surface radiation on a solid, which solid has at least one irradiable surface; by calculating the radiative exchange between gray, diffuse surfaces; characterized in that the surface or surfaces to be irradiated is adaptively subdivided hierarchically into radiation tiles of the same or nearly identical radiation intensity, and the surface temperature resulting from irradiation is effected by means of a hierarchical view factor method; which visual factors method the evaluation of a

Solid angle integral using a primary solid angle subdivision; which primary solid angle subdivision comprises homogeneous visual factor discretization, where those

Solid angle subdivision is subdivided into its subareas adaptively and hierarchically by spherical projection and wherein the sum of all partial amounts of that solid angle integral is to be determined by means of ray tracing. The invention further relates to a method as stated above, after which ray tracing is accelerated, in particular after which ray tracing is accelerated by radiation tiling, after which ray tracing is accelerated by an anisotropic Chebychev distance method, and / or after that ray tracing by means of Parallel computing is accelerated.

Furthermore, the invention further relates to a computer software product on a computer readable

Medium, which is a software code to perform a procedure above.

Bindick, Ahrenholz and Krafczyk have printed in "Heat Transfer - Mathematical Modeling, Numerical Methods and Information Technology", ISBN 978-953-307-550-1, Aziz Belmiloudi Editor, 2011-02-14, Chapter 7 - "Efficient Simulation of Transient Heat Transfer Problems in Civil Engineering "summarizes the state of the art within the technical field of thermally coupled radiation in three dimensions (see in particular Section 7.3). Their considerations are hereby incorporated into their entirety by reference in this work.

EP 1 667 069 A1 describes a method for determining the distribution of the local radiation intensity in a semitransparent medium by using ray tracing, wherein a significantly faster method with simultaneously less memory requirement for determining the distribution of the local radiation intensity in a semitransparent, at least one interface having medium was brought here. DESCRIPTION OF THE FIGURES

Illustration 1.

Illustrate the solid angles by a half unit sphere.

Figure 2.

Example of the firmly impressed discretization of the solid angle.

Figure 3.

Illustration of the chosen principle of the fixed subdivision of the unit sphere, the view "from above" along the normal direction of the radiant tile The subdivision resembles a target and consists of the azimuthally arranged segments, which in turn form radial rings.

Figure 4.

An example of the system of direction vectors for parameterization resulting from equation (7) (r = 4, n = 10). The total number of directions according to (7a) is equal to N = 86. Figure 5.

A: The second refinement plane created from the first plane with subdivision (r, n) = (4, 10), see Fig. 4. B: the third plane, created from the second plane, see Fig. 4. C: Directional vectors from all three levels. Vectors from different levels are labeled with different symbols at the origin of the vectors, cf. with Fig. 4, 5.A, 5.B. The rays for all directions shown in Fig. 5.C are sent by the ray tracing method. Figure 6.

The flow chart for the selection of the radiation sources from different refinement levels of the angle discretization.

Figure 7.

An example of the adapted discretization of the solid angle. The spherical quadrilaterals of different planes are indicated in the figure in the direction of the refinement.

Figure 8.

The scheme of ray tracing on the Cartesian grid. It consists of the detection of a sequence of intersections with the orthogonal grid.

Figure 9.

Supported variants of the tile construction, starting from single radiating sides of the orthogonal grid cells (crosses).

Figure 10.

Illustration of the Chebychev distance based on 2D examples. 10A: Isotropic Chebychev distance. 10B Anisotropic Chebychev distance.

Figure 11.

An example of ray tracing in 2D, supported by anisotropic Chebychev distance. Figure 12.

Example of a thermal result calculated with the radiation model for a precision casting project. Example of a real application.

Figure 13.

Elephant family: test of model efficiency, effects of occlusions and distance to heat source. The heat comes from the ball in the upper left corner of the picture.

Figure 14.

Modeling the simultaneous heat treatment of 4 gear units in a semi-closed cavity furnace. A: the whole model. B: tempered aggregates. The heat sources are attached to the left and right walls.

Figure 15.

Calculated scene with three cars. The heat source is up front, not shown in the picture.

Figure 16.

Heat treatment of 5 cast pump housings in an oven.

PRINCIPLES

Energy balance at the radiating surface

The consideration of the surface radiation in a thermal model is made by adapting the equation for the condition of the energy balance at the radiating surface. Without radiation, this condition is expressed solely by the conduction fluxes, (1). Ä _L T-Ä _R T) N = 0

The quantities with the indices L and R in (1) characterize the thermal conductivities and gradients of the temperature correspondingly evaluated to the left and right of the surface. The vector N is the normal vector of the surface.

An additional term, the net radiation flux density, comes into the balance equation by radiation (1):

(2). (A _LL T - λ _Κ Τ _κ Τ) iV + q _net = 0

Equation (2) must be solved in a thermally coupled model with surface radiation.

The net radiant flux consists of the difference between the absorbed and emitted radiant heat

(3). q _ne t = e (Rin ^~ ^"" Γ ⁴⁾ e in (3) is the emissivity of the surface. It is averaged hemispheric in the case of so-called gray and diffuse radiating surface over the entire electromagnetic spectrum and over the solid angle. It will also qi _n in (3) is the incident heat flux density numerical model for surface radiation is to determine the values of incident radiation at each piece of the radiating surface. Calculation of the incident radiation flux

The incident radiation flux is coupled by an integral to the outgoing radiation flux: (4). q _in = f _{i> 0} (ÜN) q _out d [l

It is integrated via the solid angle 2π. The result applies to the center of the piece of the radiating surface.

The method of view factors replaces the integral in (4) with a sum over the contributions of the individual subdivisions of the radiating surface, for which the direct visibility is given by linear optics. As a single element of the surface, in a normal case, the areas of the numerical grating apply, and the combination of several adjacent areas of the grating in such radiating elements is also known according to the hierarchical view factor method.

The solid angle is often illustrated by a half unit sphere as shown in Figure 1. A unit sphere is placed around the center of the radiating surface. The surrounding visible radiating surface lattice is centrally projected onto the unit sphere. According to the Nusselt principle, the component of the projection onto the unit sphere parallel to the equatorial plane is equal to the respective visual factor. The calculation of the visual factors is illustrated geometrically in Figure 1. To find out the value of a visual factor between surface j outside and surface i (shaded in gray in the picture), the visible part of surface j should first be projected onto the unit sphere around the center of surface i. The calculation of the total integral (4) means the projection of the entire radiating grating onto the unit sphere, as illustrated in Fig. (1).

The complexity of the quadrature (4) is generally proportional to the square of the subdivisions of the radiating surface through the numerical grid. Therefore, the calculation in (4) can lead to excessive memory and computation, especially for large, complex geometries. The geometric scheme of view factor method of Fig. (1) can be reversed, which in principle can lead to a smaller memory and computational effort.

In the case of visual factors, the subdivision (= discretization) of the solid angle is fixed by the previously generated numerical grid. The fineness of the subdivision is given completely by the numerical grid.

Conversely, it is possible to impart an arbitrary subdivision of the solid angle first independently of the numerical grating at the edge of the transparent hollow space in the first step. In the second step, this subdivision of the unitary sphere then becomes on the surrounding grids projected out. The principle is explained in Fig. 2.

This scheme leads to a different logic in the integration of the radiation flux according to equation (4). While the task in view factor methods consists of determining the individual subdivisions of the solid angle, where the subdivisions are an image of the visible subdivisions of the surface, the individual subdivisions of the solid angle are fixed in the alternative methods and therefore known in advance. The actual task now is to search for a representative element of the radiating grating that takes on the role of the radiation source for the given piece of subdivision of the solid angle. It is thus assumed that the radiation emanating from the radiation source found ^¬ radiation intensity is present homogeneously throughout the subdivision of the solid angle (See Fig. 2).

In the firmly impressed discretization of the solid angle is proceeded in the following. In the first step, a certain subdivision of the unit sphere around the center of the selected radiating surface i is performed independently of the numerical grid. In the second step, a central projection of this subdivision is performed on the surrounding grid. The direction of the projection is indicated by arrows pointing upwards, compare with the reverse projection direction in the method of the visual factors in Fig. 1. The projection of the center of each part of the solid angle is assigned to a radiating grating element, so that a mapping of the entire subdivision of the Solid angle on the radiating surface elements comes about.

An error in the discretization of the distribution of the radiation flux over the solid angle arises if the considered part of the solid angle, projected onto the radiating grating, has surface elements with very different outgoing radiation flux densities, which leads to a large variation of the radiation intensity within the solid angle.

The advantage of this method is that it eliminates the complicated geometric analysis of the occlusions. For each subdivision of the solid angle starting from the center of the surface, a single sample beam is sent. The beam direction corresponds to the center of the given part of the solid angle.

On the sending of the rays, or ray tracing, discussed in more detail in the section "Ray Tracing".

radiation tile

The method of subdividing the radiating surfaces is shown below by way of example by means of the currently used numerical grating in the program MAGMAsoft (MAGMA Foundry Technology GmbH). However, the method is in principle applicable to any type of lattice without limitation.

MAGMAsoft uses a tensor product grid. The 3D grid consists of three systems of grid lines in 3 Cartesian spatial directions X, Y, and Z, which run through the entire model and put it into a cuboid embeds. The grid thus consists of the orthogonal cuboid cells.

The numerical grid therefore consists of an orthogonal grid and is therefore completely predetermined solely by the material distribution over the grid cells in this grid and three rows of coordinates along the three Cartesian directions. The individual radiating surfaces in the radiation model are found during initialization by material neighborhoods in the lattice. A radiating surface is defined as a rectangular facet of a lattice cell with the normal vector in one of the 6 spatial directions + X, -X, + Y, -Y, + Z, -Z if one of the two conditions is met:

1. The facet separates two grid cells, one of which is covered with an opaque material and another with a transparent material. The normal vector points towards the transparent one

Cell. The transparent material in MAGMAsoft is normally air, in which case only the heat conduction is calculated.

2. The facet separates two grid cells, one with an opaque material and another with the

Material ID "boundary." The cell with the material "boundary" lies outside the computing area of MAGMAsoft. Such a facet is defined as a radiant surface if it is not at the boundary of the bounding box of the

Grid is because of the fact that there is no other facet seen from the surface. In this case, the determination of net radiant flux is trivial. The radiating surfaces defined in this way contribute to the energy balance according to equation (2). They are referred to below as radiation tile or tile.

DISCRETE THE ROOM ANGLE

The complete solid angle is subdivided in this method so that each subdivision corresponds to the same visual factor VF, namely VFj = 1 / N = const, where N is the total number of subdivisions. Such a homogenous subdivision of the unit sphere in terms of visual factors is not a homogeneous subdivision of its area because of the term (HN) in the corresponding integral:

(5). / _Ω (HV) an = / / _θ sin (G) cos (6) dcpdQ = ^ = const

The advantage of such a subdivision arises in the calculation of the incident radiation flux according to equation (4). The always same visual factor can be excluded after the determination of the incident flow is converted to an averaging of the outgoing flows.

A homogeneous subdivision of the unit sphere according to the visual factors is not unambiguous and can take place in innumerable ways. The method chosen here additionally asks for a certain symmetry in the subdivision relative to the normal vector of the radiation tile and is easy to handle. The unit sphere is first axially symmetric, starting from a circle at the north pole, divided into a series of radially successive rings. Subsequently, each ring is divided into a different number of ring segments in the azimuthal direction. Each segment is a spherical quadrilateral bounded by 2 azimuthal and 2 meridional coordinate lines (circular arcs) of the spherical coordinate system. The number of subdivisions in individual rings form an arithmetic progression (see Fig. 3).

The subdivision can be completely parameterized by the number of rings n in the meridional direction and the number of azimuthal segments of the first ring at the north pole r. For every possible parameterization (n, r), there is a single solution for the meridional coordinates of the rings, so that each spherical segment can always be assigned the same visual factor:

, _^ , 1, _Λ (n- £ + l) (2r + n + £ -2 \

7. Θ; = - ₂ arccosfvl - - -)

N in (7) is the total number of divisions after the arithmetic progression,

₍ 7a). _N = ² ^ l _{n + 1}

2

The number of subdivisions N and thus also the number of rays per one tile grows quadratically according to (7a) with the number of predetermined meridional rings n (see Fig. 4). The segments in a ring of Fig. 3 can be rotated in the azimuthal direction by a free angle, so that a larger angular distance between the segments of the adjacent rings is formed.

Hierarchical system

Coarse subdivision of the unitary sphere leads to the numerical errors in the integration of the radiant flux, as explained in the section "Energy balance at the radiating surface." In order to achieve a higher angular resolution with respect to the geometry of the radiating grating, this process becomes This method is similar to the different refinement levels of the numerical grid in an explicit multigrid method.

The first plane represents the system of spatial directions generated according to equation (7), see Fig. 4. The next plane is created by dividing each segment of the first plane by the double halving in the azimuthal and meridional directions, ie by a quarter ,

The exception is the first refinement of the circular area at the North Pole, which is divided into four azimuthal divisions in four spherical triangles. For further subdivisions, there are no exceptions, all spherical areas are refined as spherical squares, as described above. The four spherical triangles around the North Pole are treated as degenerate spherical quadrilaterals from the third level of refinement. Thereafter, the next finer level is further subdivided recursively according to the same rules. This results in a geometric progression of the number of rays on each next level. For k levels of refinement arise overall

, ^ "N _{» rtot} = N _»r -4 ^k -l = ^, 2r + nl +," .- ^k -l

(7b).

Radiance at all refinement levels.

The directional vectors after the refinement are shown in Figure 5. Three refinement levels are used in the current implementation.

In the hierarchical method of ray tracing, at first rays from all refinement levels are sent one after the other and all results obtained thereby are stored. The results are then evaluated.

Adaptation to the geometry and the temperature distribution ^¬

It would be quite possible to make virtually the hierarchical subdivision of the solid angle described in the previous section, and to send exclusively the rays for the vectors of the finest level reached. The consequence would be that all radiation sources found after ray tracing have to be processed in the thermal model, which would have led to an unnecessarily high computational effort for a large number of beams. If, for example, the set angle resolution is too high due to the number of rays in one spatial direction, the same radiation source defined by a radiation tile would be hit several times by different rays. At this point, therefore, the angular resolution from the coarser refinement levels would be sufficient.

However, if the rays are sent from all defined refinement levels, there is the possibility for local adjustment of the angular resolution. The Ray Tracing result for thermal calculation without fitting would be a complete list of radiation sources, all resulting from either the primary subdivision (Figure 4) or one of the refinement planes (Figure 5).

The local adaptation of the angle discretization replaces such a list with sublists from all discretization levels. The procedure is the following.

The subdivisions of all refinement levels are pre-marked as "unrefined." First, it is assumed that the discretization contains only the finest-level vectors, then iterates through the subdivisions of the next coarser level until it reaches the coarsest primary plane.

Each time, the radiation sources of the previous finer plane corresponding to the four quarters of the subdivided spherical quadrangle of the finer plane are examined. The previously defined active sources of radiation from the 4 quarters become either replaced by a radiation source from the current spherical segment, or maintained.

1. If at least one of the 4 subdivisions is already marked as "refined," any already defined radiation sources will be retained at the finer levels within the current spherical segment, the source from the current plane will be discarded.

2. When all 4 rays of the finer plane have hit the same opaque material, or all 4 rays have "left", ie left the semi-open transparent cavity, the already defined radiation sources of the finer plane are discarded and by the radiation source of the current plane replaced.

3. When 4 beams both leave the cavity and hit an opaque material, the already defined radiation sources from the finer planes are retained.

4. When 4 beams hit different opaque materials, the temperatures of the struck radiation tiles are compared. When adjusting the hierarchical system of directions for the first time prior to thermal simulation, sources of temperature come from the initial initialization of temperature for different materials. Otherwise, if the adjustment takes place dynamically during the thermal calculation, these are the currently calculated temperatures. If the absolute difference between the maximum and the minimum temperature with respect to the maximum temperature is less than the fixed limit, the already defined Radiation sources of the finer level also discarded and replaced by a source from the current level. Otherwise they will be retained.

5. After editing a subdivision in the current layer, it will either be marked as "refined" or "not refined" depending on the result. If it comes to replacing the sources of the finer levels by a source of the current level, this is included in the list of sources considered the current

Refinement level entered. If one of the sources is marked as "not refined," then Source is included in the list of considered sources for the finer level, and entry is made only if the new source belongs to a surface of the opaque material, otherwise the visibility factor of the outer one will not hidden space by a constant value of the current space

1

Level of refinement ι same-- increased.

The procedure described above is shown in Figure 6. The algorithm provides k lists of radiation sources for k levels of homogeneous

Discretization. The adapted angular discretization is chosen such that a compression of the representative radiation sources arises in the angular ranges, where a discontinuous dependence of the radiation intensity on the spatial direction is potentially to be expected.

The jump in radiation intensity occurs in the half-open cavities at a boundary between the hot opaque material and the open space. One such jump is also caused by the thermal contrast of two opaque materials adjacent in angular space. The materials labeled with different ID ^x s may either be in direct thermal contact along their boundary in the angular space, or partially obscure another.

The adapted discretization of the solid angle for the observed radiation tile is done in this way by the selection of the radiation sources. The subdivided surface of the unit sphere can be used for the representation of angular discretization. It is covered completely and without overlapping by spherical segments from different refinement levels after adaptation.

An example of the adjusted discretization of the solid angle is shown in Figure 7. A visual factor of an H-shaped surface, shown in the figure on the top left, is to be approximated by ray tracing. For this, 3 discretization levels with the first level parameterized with (n = 15, r = 4) according to equation (7) were used. The number of beams is equal to 166, 664, and 2656 for the corresponding 1-e, 2-e, and 3-e discretization levels, so that a total of 3486 beams are sent for hemispheric scanning. After performing the geometric fit, 218 beams have been selected, which is smaller by a factor of 16 than at the finest level. The spherical quadrilaterals of different planes are marked in the picture sorted in the direction of the refinement by their size. The outgoing radiation fluxes, given by sources from different planes, must be weighted differently in the calculation of the incident radiation flux. The equation (6) for the incident radiation flux is supplemented by prefactors of different levels: For three levels, k = 3. Three sums with the weight factors 1, H and 1/16 give (8) the total incident radiation flux. The terms q _o ^] _ut in the inner sum are the contributions from the list of radiation sources of the respective discretization levels "level".

RAY TRACING A voxel-assisted ray tracing method with backward tracking of the rays is used.

A voxel supports the test of whether the ray hits an object in the given voxel. A voxel is a cube-shaped volume where the information about geometric objects is contained in it. When used with MAGMAsoft, a tensor product lattice offers itself as a raster consisting of the individual voxels. Each grid cell becomes a voxel. However, the method according to the invention is not bound to the program MAGMAsoft but can be used for any voxel-based ray tracing method. Backward tracking means that all rays from the tile receiving radiation are first sent out. The radiation sources hit by radiation are determined. However, the radiation sources physically send energy to the receiver. The energy reaches the tile in the "straight" path, which is the same as in Ray Tracing's ray path, so this is referred to as ray tracing with backward tracking.

The geometric objects that absorb and reflect the rays are the defined radiation tiles that are defined by the material neighborhood in the grid cells.

Due to the structured nature of the numerical grid, it is sufficient to store 3 ID ^x s of the tile per grid cell. When referring to 6 sides of a grid cell as front, back, west, east, north, south, the 3 IDs refer to ^x s on three possible tiles on the back, east, and north sides of the cell. The information about the tile on three remaining sides of the cell Top, West and South can be taken from the neighboring cell. The positive values are assigned to the real existing tiles on one side of the cell. The assignment takes place during the definition of the tile.

Before the start of the ray tracing, the entire once-defined bundle of the vectors, which specify directions of individual rays, is centered around the normal vector of the respective tile. This is achieved by multiplying the vectors with a rotating matrix that is the central vector of the bundle in Cartesian direction + Z converted into the normal vector of the tile.

Thereafter, the rays from all refinement levels are sent in order. The ray tracing through the numerical grid consists of finding intersections between the continuation of the ray and the individual sides of the grid cell in which the ray is currently located (see Fig. 8).

For this, three possible sides of the grid cell are checked for the next overlap. The possible sides are characterized by the sign of the three components of the

Beam direction Ω given. The intersection with the minimum distance from the current position of the beam is considered the next point on the beam. The side of the cell with the minimum length to the point of intersection is searched, and the ray to the found

Intersection point in the direction of Ω continued around the found length X ^ray :

(9a). X _MIN = M / N _{L = (ij} - _fc) - ^ -

(9b). l ^ X ™ * = x _MIN n _t

After the next intersection is found, ID of the corresponding tile is queried on the side of the grid cell. · If ID equals a real tile, d. H. the

Side of the cell separates an opaque and a transparent material, the pursuit of the Beam stopped and the found global ID of the radiation source returned.

• If ID corresponds to a plane of symmetry with a Cartesian normal vector, the normal direction component to the plane of symmetry is inverted and the reflected beam is traced.

Otherwise, the ray tracing procedure in the next grid cell is repeated until the ray hits a tile or leaves the boundaries of the numerical grid. In the latter case, a fixed ID of the outer space is returned.

ACCELERATION Ray Tracing Tile Association

The planes of finely networked radiating surfaces (plates) often occur in the modeled geometries. They produce many radiant elements. The calculation effort in the radiation modeling is linearly proportional to the number of radiation tiles with the presented model. Computing time can be saved if the fineness of the cross-linking provides the necessary resolution in the distribution of the integrated radiation flux

(= irradiance) exceeds the surface. In this case, it would be sufficient to have less coarse tile than imprinted by the numerical grid.

In published procedures such as the DTRM model of FLUENT is based on clustering of Radiation tiles pointed. There, the union of several adjacent radiation tiles is made, if they deviate from a plane only by a minimum angle (planarity). The ray paths and the radiation fluxes are calculated not for each component of the cluster, but for the entire union at the location of its geometric center. A similar technique is implemented for "cheap" acceleration in current practice.

All radiation tiles have a Cartesian orientation in the current model. The union of the neighbors is allowed if they are both in one plane and their grid cells have the same opaque material.

It is allowed up to 3 united tiles in each of the two lateral directions along the surface. Thus, from 1 to 9 sides of the grid cells can be combined into one tile. The approved variants are shown in Fig. 9.

Clustering of the radiation tiles usually only occurs massively on flat surfaces. For curved geometries, the surface lattice is staged, which reduces the number of coplanar sides of the lattice cells.

Acceleration by a factor of 2 or slightly above is achieved on complex geometries on average by the implemented method of tiling. Program-technically, the tiling in clusters consists of the following steps: • The list of grid cells belonging to a cluster radiation tile is stored in the data structure of a radiation tile.

• The ray tracing starts from its geometric center for each cluster.

• Information about the found and matched radiation sources is also stored for the entire cluster rather than for every smallest component.

· The emitted radiation flux of the cluster is with

Consideration of temperature distribution calculated:

(10). q _em = e ^ - with S area of the side i and T absolute temperature at the surface of the grid cell.

The calculated incident radiant flux is the same for each member of the cluster resulting in localized smearing of the smog

Radiation flux leads. But that's the price for the acceleration.

· The average incident radiation flux is included in the same way for each side of the grid cell from the cluster into the energy balance (2).

Anisotropic Chebychev distance

The ray tracing according to the voxel-based ray tracing method requires the visit of each transparent grid cell lying on the beam path, see Fig. 8. The processing of many empty (ie transparent) cells can take a substantial part of the computation time, especially if the model is large Air-filled spaces or the fine networking in the spaces has. A procedure would be the Ray Accelerate tracing that instantly "tunnel" a ray over a larger block of empty cells instead of visiting each grid cell hitting the road, and the information about the size of such empty cell blocks should be present.

Such empty blocks can be explicitly defined in the grid. A better method is the so-called Chebychev distance, also called chessboard distance. It is about a distance measure for discrete objects such as grid cells, which is measured in whole numbers.

An illustration of the Chebychev distance using 2D examples is shown in Figure 10. The procedure is as follows.

The distance from the painted grid cells (null cells) is measured. Figure 10A shows the isotropic classical Chebychev distance. The cells of equal distance around the zero distance zero-lined cells are arranged in a layer on the edge of a square. Figure 10B shows the anisotropic Chebychev distance shown for the upper left quadrant. Accordingly, the anisotropic distance is defined exclusively only for the cells lying in this quadrant with respect to one of the null cells. Otherwise, it is measured the same way as the isotropic distance in Fig. 10A. The considered directions (quadrant) are marked in both cases with red arrows in the zero cells. It is noticeable that the anisotropic distance reaches greater values than the isotropic distance.

Using the Chebychev distance makes steps in ray tracing adaptive, depending on the current one Make distance to a surface. If the current position of the tracked beam is far from the nearest surface, the jump in the grid will also be correspondingly large. As the beam approaches a surface, the Chebychev distance decreases, and so does the jump magnitude. The closer a beam comes to the surface, the slower its movement becomes. The isotropic classical Chebychev distance is measured in relation to the labeled grid cells with the distance zero. Such cells are referred to below as null cells. The cells of equal Chebychev distance from a selected null cell are arranged in square layers around the null cell.

The distance is obtained by searching for the not yet marked grid cells. If, at the end of the first step of the algorithm, all next cells from the first layer with a distance equal to 1 are found and marked, the next layer is occupied by distance 2. For this, all unmarked cells with direct neighbors or neighbors are searched for "over edge" with distance 1. The same applies to the other layers until no cells without a defined distance are found.

The effort for a direct calculation of the Chebychev distance behaves correspondingly as where N is the number of grid cells.

However, this dependence is not compulsory and can be optimized if, when occupying each next layer of the cells with distance, only the previously found and temporarily stored grid cells with a distance i be processed. Only their neighboring cells are examined. This propagates the front of the last found "active" cells, but it assumes that write access to the outer cells of the stencil is possible, so it enforces restrictions on parallelization through area decomposition.

The procedure for calculating the anisotropic Chebychev distance is very similar. However, a distance for a cell in the isotropic case is now replaced by 2 ² = 4 distances in a 2D case and 2 ³ = 8 distances in a 3D case. Each of the 4 distances is calculated in the 2D case only for grid cells located in the corresponding quadrant, seen from zero cells, see Figure 10B. The assignment of the distances there is calculated for the quadrant between the Cartesian directions -X and + Y. In 3D, a distinction is made between the 8 octants. In both cases, in contrast to the isotropic Chebychev distance, there is no assignment of the distances in the cells, which in the given quadrant (2D) or octant (3D) have no zero cells from which the distance is measured.

The determination of an octant is conditional on the use of the anisotropic distance in the ray tracing. It is first checked in which

Octant the beam direction Ω is. Then one of the 8 anisotropic distances in the 3D case is accessed, which corresponds to an opposing octant. If z. If, for example, a ray comes from the direction (0.5, 0.2, -0.1), it lies in the octant (+ X, + Y, -Z). The corresponding distance, starting from the zero cells, was measured in an octant opposite to the beam direction. The necessary distance is thus stored in the octant (-X, -Y, + Z). The advantage of an anisotropic distance over an isotropic is a more accurate default for a possible jump across a block of empty cells in the ray tracing. The value of the anisotropic distance in an octant is often larger than in the isotropic one, see Figs. 10A and 10B. The probability that a null cell lies in an octant cube of the grid is smaller than in the isotropic case where the corresponding cube is composed of 8 such octants. Therefore, ray tracing can make larger cracks in the grid.

A disadvantage of the anisotropic distance is the additional memory required for this

Ray tracing, which is an extra 96 bytes per grid cell for the type "int" of the distance.To save memory, distances are taken with the type "unsigned char" in this method. This type allows a maximum size of the distance of 128 grid cells. Theoretically, larger values can occur in one application. Because of this, the assignment of the distances in their calculation is interrupted after reaching the layer with value 128, and the remaining unoccupied transparent cells are assigned the same value of 128. This limits the maximum possible length of a jump in the grid during ray tracing. In the method, the null cells are defined as the opaque cells with defined radiation tiles and the cells at the outer edge of the numerical grid in contact with a transparent cell. The last one ensures that all transparent cells are seamless are enclosed by a layer of zero cells, and for each transparent cell at least one reference zero cell is also present in the anisotropic case. An example of 2D ray tracing, supported by anisotropic Chebychev distance, is shown in Figure 11. The grid and the calculated distance are taken from Fig. 10 right. The ray in the octant (+ X, - Y) is processed with distance in the octant (-X, + Y). The starting position of the beam at the intersection with one cell side at the top left is in the cell with distance 3. Therefore, the next intersection with the side of the 3x3 cube is searched. After shifting the position to the intersection, the Chebychev distance is taken from the next cell in the direction of the beam path. She is equal to 3 again, and the next jump leads the beam directly to the target on the lower right. In this example, 2 intersections were searched; without applying the Chebychev distance would be the whole 9.

The additional cost of Chebychev distance ray tracing is to determine the discrete position of the beam after the jump over a block of grid cells. The coordinates of the point of intersection can be similarly calculated according to Equations (9). Instead of the Cartesian coordinates on the edge of the grid cell X _L , however, the coordinate at the edge of the octant cube, which is the coordinate of a more distant grid cell is the same.

The indices of the cell hit after the jump in 3 Cartesian directions are unknown. They are determined according to the binary search algorithm. The index in the grid at the intersection is calculated from an equidistant Grid in the cube estimated in each direction. Then the coordinate comparison checks whether the found intersection point lies in the estimated index cell. If not, the corresponding half of the Chebychev distance interval is further halved. The procedure is repeated recursively in this way until the found index of the cell belongs to the intersection. The effort of the binary search in a fixed interval depends logarithmically on the interval length, ie for a maximum jump length of 128 cells, ~ ln (128) = 81n (2) ~ 5 comparison operations per direction would be necessary.

Despite the overhead of gridding, the Chebychev distance provides a considerable acceleration of the process, especially for models where a lot of space is occupied by transparent material.

Parallel ray tracing

Load distribution with ray tracing:

Good scalability in parallel operation is the important prerequisite for initialization of the model (Chebychev Distance, Ray Tracing) and its application in thermal computation. The decisive factor is the computing time for ray tracing.

The common practice for parallelization with distributed memory like MPI in the solution of partial differential equations is the decomposition of the computational domain. Their application for ray tracing, however, has shown poor scalability in this development. The reason for this are high costs of Communication. Data of the millions of rays must be transmitted again and again between the individual partitions. The specific weight of the communication grows as the number of CPUs increases, eventually leading to poor scalability. Because of this, another method of parallelization has been developed in this model.

Ray tracing and Chebychev distance calculation are performed simultaneously on each CPU for the entire grid. However, the radiation tiles remain locally defined on each CPU.

The global grid model is created in the first section. The grid data from the entire grid is communicated to all CPUs. On the one hand, three rows of coordinates are communicated in three Cartesian directions. On the other hand, a transparency flag of each grid cell is determined by root and communicated to all CPUs in a numeric field of type "unsigned char." This data is sufficient to perform ray tracing on the tensor product grid independently on each CPU The symmetry limits in the thermal model are also considered here.

However, the data locally determined by Ray Tracing about the radiation sources of the radiation tiles should finally reach the CPU, where the corresponding tile is locally defined. At this point costs for communication arise. To minimize them, the procedure is as follows. 1. The number of existing on a CPU locally defined radiation tiles N is determined.

2. Then the arithmetic mean N _{av is} calculated via the CPUs.

3. From a CPU donor with the excess of tiles

> N _av , a portion of AN of the tiles is virtually subtracted from the donor and attributed to the next found CPU acceptor with N ₂ <N _av , so that after this operation either the condition N ₁ = N _av or N ₂ = N _{av is} met becomes. The transferred

Portion is equal to AN = M1N (N ₁ - N _av , N ₂ - N _av ). The indices of the transferred tiles and the ID ^x s of CPUs acceptors are saved each time. 4. Step 3 is repeated until it is no longer possible to balance the tiles between the CPUs.

After the load balancing parameters become known to each CPU as described above, the local data about the location of the radiation tiles in the grid are communicated from the CPUs donors to the CPUs acceptors. In this way, the CPUs acceptors receive a list of the imported tiles to be processed in addition to the list of their own tiles. The lists of custom tiles to be edited on the CPUs donors are trimmed accordingly.

Parallel ray tracing:

1. Before ray tracing, the anisotropic Chebychev distance is calculated independently on each CPU. No communication is required here, calculations and results on each CPU are identical at this point. The effort for This calculation takes a few seconds, so that the organization of a work distribution between the individual CPUs and subsequent communication would hardly be worthwhile here.

During ray tracing, each beam independently sends each beam from its own original tiles to be processed, finds the radiation sources, and stores the selected sources according to their geometric and thermal adaptation. The storage of the radiation sources takes place in this step directly into the data structure of the own radiation tiles.

After all of your own tiles have been edited, it is the imported tiles on the CPUs acceptors' turn. The same procedure of ray tracing and radiation source detection is used in the processing of these tiles. The found and compressed by the procedure of adaptation radiation sources are buffered in a buffer.

When all the tiles have been edited, the results are communicated. The CPUs acceptors return the found radiation sources of the imported tiles to the CPUs donors. The MPI communication is asynchronous, and with an irregular distribution of tiles, each donor may have multiple acceptors and each acceptor may have multiple donors. The radiation sources received from the donor are written into the data structure of their own tiles in this step just after receiving from a CPU acceptor. 5. At the end of ray tracing, data that is no longer needed, such as chebychev distances and communication buffer, is deleted., THERMAL CALCULATION

In order to solve the summation of the outgoing radiation fluxes according to the multi-level method according to equation (8), a suitable algorithm for the approximate solution of linear equation systems is used. The summation is exemplified below on the basis of the Gauss-Seidel iteration, but can also be solved by other algorithms, for example the Jacobi method or the SOR method.

In one embodiment, the Gauss-Seidel iteration is used for the calculation of the absorbed

Radiation flow used on each radiation tile. This flux represents the end result for the thermal model from the radiant side. The knowledge of the outgoing radiant fluxes at this point is required across CPUs since the considered beam path was calculated for each beam regardless of the area decomposition.

The advantage is that the Gauss-Seidel iteration on each CPU is done only for the local radiant tiles. The summation of the outgoing radiation fluxes according to the multi-level method according to equation (8) includes the contributions of all tiles determined by ray tracing, also from other CPUs.

This allows the outgoing radiation fluxes to be made globally visible in the algorithm. Everyone Radiation tile is assigned a continuous global index, and advantageously the same indices can be stored instead of the radiation sources for each tile during ray tracing.

The global outgoing radiation flows can be indexed in an array after the global index, where before each pass of the Gauss-Seidel loop, the outgoing flows from all the local tiles are communicated to the global array.

POTENTIAL APPLICATIONS, EXAMPLES

The above-presented model of thermally coupled surface radiation for diffuse gray radiators may find application in many engineering processes where voids or semi-transparent materials and high temperatures are present. The attractive side of the process is the ease and very short computation time, in which the most demanding tasks in the modeling of the long-range effect of the thermal radiation are processed.

The most interesting applications in the simulation of the casting processes are investment casting, heat treatment, ingot casting in the aligned arrangements of the plants.

Further applications are related to the production and operation of industrial furnaces, for example in glass production or ceramics production, in chemistry, crystal growing, metallurgy or even in the food industry, eg ovens. Various examples of application of the method are shown in Figures 12 to 16. As shown, the present method is not limited to just one heat source but may be used for a variety, eg, two (Figure 14) and a higher number of heat sources.

The method given here for discretizing a solid angle for use in a simulation process is to be considered in its application not only in connection with thermally coupled surface radiation, but can also be used in other calculation and simulation processes where acceleration and savings of computer time and used computer memory, such as For example, rendering computer software (rendering software) in general and computer games in particular or scientific imaging tools is of importance. The method of ray tracing given here for use in a simulation process is equally applicable in its application not only in connection with thermally coupled surface radiation, but can also be used in other calculation and simulation processes where acceleration and savings of computer time and used computer memory, such as For example, rendering computer software (rendering software) in general and computer games in particular or scientific imaging tools is of importance. In particular, the unconventional method of parallel computing proposed here for accelerating a ray tracing calculation alone and in combination with anisotropic chebychev distance calculations and / or additional acceleration by tile union is not known to the inventor of the prior art.

Claims

Method for simulating the influence of thermally coupled surface radiation on a

Solid, which solid has at least one irradiable surface; by calculating the radiative exchange between gray, diffuse surfaces; characterized in that the surface or surfaces to be irradiated is adaptively, hierarchically subdivided into radiation tiles of the same or almost identical radiation intensity, and the surface temperature resulting from irradiation is effected by means of a hierarchical view factor method; Which

Visual factor method comprises the evaluation of a solid angle integral using a primary solid angle division; which primary solid angle subdivision comprises homogeneous visual factor discretization, where those

Solid angle subdivision is subdivided into its subareas adaptively and hierarchically by spherical projection and wherein the sum of all partial amounts of that solid angle integral is to be determined by means of ray tracing.

The method of claim 1, wherein said ray tracing is accelerated.

The method of claim 1 or claim 2, wherein said ray tracing is accelerated by radiation tile merging. A method according to any one of claims 1 to 3, wherein said ray tracing is accelerated by an anisotropic Chebychev distance method.

Method according to one of claims 1 to 4, wherein that ray tracing is accelerated by means of parallel computing.

A computer software product on a computer readable medium comprising software code for carrying out a method according to any one of claims 1 to 5.