Modelling Techniques for ultra-porous heat-protective materials’ spectral properties

A model of the interaction of electromagnetic radiation with representative orthogonal elements of highly porous thermal protection fibrous materials and reticulated materials is developed. This model, which is based on the Mie theory and its consequences, makes it possible to calculate such important optical characteristics as spectral absorption and scattering coefficients, spectral and spectral transport diffusion coefficients of radiation, phase scattering function, radiation thermal conductivity, and others.


List of symbols
Greek symbols α Angle of radiation incidence α λ , β λ , γ λ Spectral absorption, scattering and extinction volumetric coefficients Plane δ Dirac function Solid angle element in the grid of directions ε Dielectric constant ζ Angle of the scattering plane orientation θ Polar angle, angle of scattering λ Wavelength μ Magnetic permeability π n , τ n Auxiliary functions of the scattering angle τ xy Angle between the x-and y-cylinders apertures, etc. ϕ Azimuth, angle between the incidence and scattering planes χ Angle between the polarization plane and incidence plane

Introduction
Ultra-porous materials, such as light fibrous materials (FM, porosity is up to 90 %) and reticulate materials (RM, porosity is up to 97 %) are attractive for many applications. These materials consist of chaotically oriented fibres or have rigid spatial skeleton formed by junctures and struts. Such fragments may be made of the same or different materials. The space between the fibers or elements of the reticulated skeleton of the material may be filled with gas.
In Cherepanov (2009), Cherepanov (2011), discrete statistical models of the structure, thermophysical, and electrooptical properties of highly porous fibrous and reticulated thermal protection materials were described based on the Mie theory and direct simulation Monte Carlo method. In these models, the authors were successful in combining a fairly accurate description of radiation processes, inclusion of anisotropy, and microstructure features of a material with the consistent application of the concept of representative element. The application of models made it possible to considerably improve the informativeness of experimental studies of materials, as well as determination and prediction of heir most important thermophysical properties (Alifanov et al. 2011). However, the models were restricted in what concerns the direction of illumination of representative elements. Similar restrictions are removed in more advanced model of spectral optical properties of ultra-porous heat-protective materials. The description of this model and of some its possibilities makes a subject of the given paper.
The practical experience of statistical modelling is showing that the features of geometry of separate representative element, especially of the its fragments, are not so significant for the optical and thermophysical properties of ultra-porous materials. At modelling, when the sample volume includes thousand or even tens of thousands of a representative element, the result is essentially defined by the statistics. Therefore, at the modelling of FM and RM it is quite possible to use the universal enough orthogonal representative element, as in works (Alifanov and Cherepanov 2009;Cherepanov 2011;Alifanov et al. 2011).
When the analytical description of the problem of interaction with electromagnetic field, cylindrical and spherical body allows you advance the most far ahead, maintaining the accuracy of the description at a level sufficient to account for such subtle but significant effects, as resonance scattering and absorption of radiation (Wait 1955;Lind and Greenberg 1966;Bohren and Huffman 1983). The optimized forms of such an analytical solutions presentation that can significantly reduce the computational time are known (Bohren and Huffman 1983). Nevertheless, prisms often form the representative elements of some materials (Zhao et al. 2004;Placido et al. 2005;Loretz et al. 2008;Öchsner et al. 2008).
However often the RM skeleton, for example reticulate vitreous carbon (RVC) skeleton, is used as a solid and heat-resistant basis for various substances, that besieged from gas phase. Such a process considerably is changing the struts cross sections, doing their more rounded (Colombo 2006). Therefore, when we modelled the spectral optical properties of RM representative elements we used struts of circular cross-section, instead of struts with the form of cross sections close to a triangular form.
The structure analysis of RM RVC shows a significant number of accrete junctures. Owing to it, in the model of RM, the quantity of the struts converging in one juncture, varies from three to seven. The representative element of FM (a), variant of RM representative element (b), and distribution of the junctures versus number of their struts for RM RVC (c) are shown on Fig.1. We make the following assumptions, which are conventional to a certain degree in discrete modelling of the optical properties of such materials (see Wait 1955;Lind and Greenberg 1966;Bohren and Huffman 1983): 1. On the spatial scales from fractions of micron and above, the components of the incident field and the field scattered by particles of the material can be related using the rigorous electromagnetic theory or, sometimes, scalar diffraction theory, which is its consequence. 2. The electromagnetic wave incident on a representative element of the material (see Alifanov and Cherepanov 2009;Cherepanov 2011) is scattered by individual fragments elastically and independently according to the Mie theory. The radiation scattered by a fragment and the secondary radiation emitted by it are neglected in the radiation that is incident on other fragments. 3. The sizes of the scattering fragments of the material is taken into account wherever possible. 4. The direction of illumination of the representative elements must be consistent with the local properties of the radiation. 5. The optical characteristics of the volume element of the material are typically obtained by averaging the characteristics of the individual fragments with weights proportional to their cross sections in the absorption and scattering processes. However, various hypotheses, for example, the assumption that the radiation is in local equilibrium, can affect the averaging procedure (Moiseev et al. 1991). 6. The medium around the material fragments does not absorb the electromagnetic radiation.
At the solving of radiation transfer problems, the material, as a rule, is replaced with system of independent scattering bodies (Öchsner et al. 2008;Zeghondy et al. 2006;Petrasch et al. 2007). However, for environments with spherical inclusions known examples of correcting received independent solutions by introducing some amendments to the cooperative effects (the near-field factors) (Drolen and Tien 1994). Such an approach has not some worthy alternative today, because the exact solution of a scattering problem for a structures type of representative elements, within the limits of the rigorous electromagnetic theory does not exist. In offered above models, these cooperative amendments are carried out by means of the corrective multipliers for the spectral efficiencies of absorption and scattering, which were calculated for fragments of the representative element based on the Mie theory. These multipliers possess the simple enough physical sense and are defined by comparison of the modelling results with the experimental data for spectral absorption coefficient and transport diffusion coefficient of radiation (Moiseev et al. 1991).
In the Mie theory, the spectral optical characteristics of particles are represented in terms of coefficients of the expansion of the incident and scattered fields in the complete Maxwell system of vector eigenfunctions. The specific form of the eigenfunctions is determined by the geometry of the scattering body. The ideas behind the Mie theory are described in the classical works (Wait 1955;Lind and Greenberg 1966;Bohren and Huffman 1983), and some details that are important for the present paper can be found in Cherepanov (2011). In the statistical models FM and RM (see Alifanov and Cherepanov 2009;Cherepanov 2011), orthogonal representative elements formed of such simple bodies as a ball and a right circular cylinder ( Fig. 1) are used. That is why we will show the most important results of the Mie theory for these objects.

Scattering of radiation by a ball and a cylinder
Scattering by a ball is independent of the polarization of the incident field (Bohren and Huffman 1983). In this case, the following expressions for the cross section C, efficiency Q, asymmetry parameter cosθ , and the scattering phase function p for the monochromatic radiation hold (see Wait 1955;Bohren and Huffman 1983): [n(n + 2)Re(a n a * n+1 + b n b * n+1 ) + 2n + 1 n Re(a n b * n ), [(a n a * m + b n b * m )(τ n τ m + π n π m ) +(a n b * m + a * m b n )(τ n π m + π n τ m )]. The quantities determined by the absorption, scattering, and extinction processes are marked by the indices a, s, and e, respectively. The angle brackets · · · denote the averaging operation. The coefficients a n and b n of the scattering series appearing in these expressions are as follows (see Wait 1955;Bohren and Huffman 1983): . (1) Here, j n (x) and h (1) n are, respectively, the spherical Bessel and Hankel functions of the first kind and order n; m = N 1 /N is the complex refractive index of the scattering body relative to the environment; x= ka = 2π N a/λ is the parameter of diffraction for the environment; λ is the wavelength in vacuum; N , and N 1 are the relative complex refractive indexes of the environment and the ball relative to the vacuum; and ε and μ are the dielectric constant and magnetic permeability of the ball relative to the environment. Since the polarization is irrelevant in the case of scattering, the strengths of electric E and magnetic H fields are transformed almost identically due to the symmetry of the Maxwell equations, with the natural replacement of ε with μ. Relations (1) emphasize this symmetry of the fields (Lind and Greenberg 1966;Bohren and Huffman 1983).
The functions π n and τ n , which depend on the polar angle θ (in the problem of scattering by the ball, the polar axis Oz is oriented in the direction of its illumination) can be conveniently calculated using the recurrences (Bohren and Huffman 1983): The Mie theory does not give scattering formulas for cylinders of a finite length. Such a solution is known for a straight infinite circular cylinder of radius a (Lind and Greenberg 1966;Bohren and Huffman 1983); however, its form considerably depends on the polarization of the incident field. Two independent versions of polarization of the incident field are distinguished.
In version E, the electric vector lies in the incidence plane xOz (the axis Oz is parallel to the cylinder axis); in version H, the magnetic vector lies in this plane.
If the incident wave is polarized linearly and its electric vector E i is at an angle χ to the plane of incidence xOz, then the complete picture of the relation between the amplitudes of the incident and scattered fields can be obtained (Bohren and Huffman 1983): This picture is needed, for example, to obtain the scattering phase function in the case of the oblique (with respect to the cylinder) incidence wave. In (2), α is the is the incidence angle of the wave with respect to the cylinder axis; for the components of the incident (i) and scattered (s) field, different spherical coordinate systems are conventionally used (Wait 1955;Lind and Greenberg 1966;Bohren and Huffman 1983). The front of wave (2) (see Fig. 2) is a conical surface with the axis coinciding with Oz. Angular half-width of the cone constant phase is equal to the wave incidence angle α on the cylinder axis. It is distributed in the direction of wave normals e s , which also form a cone, coaxial to cylinder. Its generatrix forms the angle of π/2-α with the cylinder axis. The elements T n of the amplitude scattering matrices in (2), the Poynting vector of the scattered field S sχ and the direction vector of its propagation e s are determined by the following relations (Lind and Greenberg 1966;Bohren and Huffman 1983;Cherepanov 2011): a H n cos(nφ), (3) ], e s = e r cos α + e z sin α. (4) Fig. 2 Scattering of a plane wave by a cylinder with the axis Oz: k i and χ are the wave vector and the polarization angle of the incident wave, e i , e s , e s0 are the unit vectors of incidence, scattering, and forward scattering directions, α is the incidence angle, θ is the scattering angle, ϕ is the angle between the incidence and scattering planes (4), E b H are the complexvalued strengths and * denotes the complex conjugation), one can introduce the efficiencies of radiation extinction and scattering per unit length of the cylinder as They satisfy the following equalities: in which the superscripts mark the characteristics of radiation in independent polarization states, and the efficiencies without superscripts relate to the unpolarized (natural) radiation. Expressions for the coefficients a n and b n of the scattering series can be found (as well as formulas (5)) in Lind and Greenberg (1966), Bohren and Huffman (1983). The spherical angles θ and ϕ, which determine the direction of scattering e s , and the incidence angle α are related by the equation cos θ = e s · e i = sin 2 α − cos ϕ cos 2 α.
Due to the dependence θ = θ 1 (α, ϕ) determined by (6), scattering by an infinite cylinder is singular; therefore, independent integration over the polar angle for the evaluation of scattered energy fluxes for different directions is impossible. However, we can integrate over the cone formed by vectors e s and determine the scattering phase function p of monochromatic radiation (Bohren and Huffman 1983) where δ is the Dirac function and x-the cylinder diffraction parameter. The asymmetry parameter for scattering this radiation is Cherepanov (2011) cos If the cylinder radius is so small that the conditions x, |m| x << 1 are fulfilled, then it is sufficient to take into account only the lowest order terms (with respect to the parameter x) in the expansions of the coefficients a n , b n and in the matrix elements T n . Following Bohren and Huffman (1983), using the asymptotic representations of the Bessel and Hankel functions, and retaining the terms of order O(x 2 ), we obtain In the case of normal incidence of the radiation on axis of the cylinder and μ =1, these relation become identical to those given in Bohren and Huffman (1983). Within the scalar diffraction theory, cylinders of finite length can be considered. The amplitude scattering function for the cylinder can be found in Bohren and Huffman (1983). The scalar theory describes the position of the first several maximums of the phase Mie function, which depends on the scattering angle θ , fairly well. However, as the scattering angle increases, the accuracy of the scalar theory steadily decreases. It is also noted in Bohren and Huffman (1983) that the infinite length approximation for the cylinder of length L may be used if K = Lcosα/2a > 10.
However, this criterion is incomplete because scattering by the cylinder intricately depends on the incidence angle and the diffraction parameter. The paper (Lind and Greenberg 1966), which contains the results of both theoretical and experimental study of the optical properties of cylinders, is a key work in this respect. In particular, it is noted in that paper that at α= 90 • the efficiency of extinction Q e for the cylinder of length L decreases with increasing length as L −1 ; thus, the cylinder illuminated along its generatrix does not attenuate the radiation flux. Experiments showed that, for example, at m= 1.26, 2 ≤ x ≤ 5, and α ≤ 50 • , the deviation of the experimental results from the Mie theory did not exceed 15 % already for L= 7λ. For L > 20λ the efficiency of extinction for the longitudinally illuminated cylinder became less than unity.
As a rule, the radiation incident on the cylinder at right angles is scattered most intensively. However, an interesting effect of resonance that occurs in scattering at large incidence angles close to 90 • was mentioned in Lind and Greenberg (1966). As the relative length of the cylinder L/λ decreases, the resonance builds up; therewith, the resonance domain moves to greater values of the diffraction parameter x. For infinitely long cylinders, such a resonance is not observed for large values of the diffraction parameter. These facts must be taken into account when the infinitely long cylinder model is used.

Scattering of radiation by representative elements of highly porous materials
Suppose that unpolarized monochromatic radiation (Fig. 3) is incident on a representative element (Fig. 1) of a highly porous material in an arbitrary direction determined by the angles θ i and ϕ i of the global spherical coordinate system. The polar axis Oz of such a system is typically oriented in the direction in which the external heat flux is incident on the material. Note that the cylindrical fragments of representative elements in models (Alifanov and Cherepanov 2009;Cherepanov 2011) are oriented along the axes of the global Cartesian system of coordinates. Each such fragment of the representative element has an individual incidence plane i . For the ball, it contains the axis Oz and the e i -the vector of the radiation incidence direction . In this case, the scattering results are independent of the incident azimuth. Therefore, the juncture of the representative element of foam material can be described using the relations of the ball presented above in which θ must be replaced with arccos((e s , e i )).
For cylindrical fragments, the incidence plane iν (ν = x, y, z) contains the vector e i and the axis direction vector e ν of the scattering cylinder. We consider the vector e ⊥iν , which is coaxial with e i × e ν (Fig. 3), as normal to the incidence plane iν .
For cylinder with the axis Oz (z-cylinder), the relations written in the preceding section can be used with ϕ replaced by ϕ − ϕ i and the algebraic incidence angle defined as α z =90 • -θ i . Indeed, the system of coordinates of the representative element is chosen and fixed so that the positive directions along the axes of the cylinders cannot be chosen arbitrarily as in Lind and Greenberg (1966), Bohren and Huffman (1983). The use of algebraic incidence angles α ν ∈ [−90 • , 90 • ] is caused by the necessity of differentiating between the cases of radiation slip in the positive and negative directions relative to the constantly oriented axes of the cylinders.
Let the normal to the scattering plane e ⊥sν = e sν × e ν /|e sν × e ν | form an angle ζ ν ∈ [0, 2π] with the axis Oz (Fig. 3a). As this angle changes, the rays of the scattering directions form a cone with the angular half-width 90 • − |α ν | and the director circle that is shown by a dotted line in Fig. 3b for ν = x. The scattering direction vectors can be represented in terms of the angles ζ ν and α ν it is easy to calculate both the scattering angle and the angle between the incidence and scattering planes, which are the arguments of the phase function in the scattering problem by individual fragments e sx · e i = sin θ i (cos ϕ i sin α x − cos ζ x sin ϕ i cos α x ) − cos θ i sin ζ x cos α x , e sy · e i = sin θ i (cos ζ y cos ϕ i cos α y + sin ϕ i sin α y ) − cos θ i sin ζ y cos α y , e ⊥ix · e ⊥sx = − sin ζ x cos θ i + cos ζ x sin θ i sin ϕ i sin 2 θ i sin 2 ϕ i + cos 2 θ i , e ⊥iy · e ⊥sy = − sin ζ y cos θ i + cos ζ y sin θ i cos ϕ i sin 2 θ i cos 2 ϕ i + cos 2 θ i , θ s,ν = arccos(e sν · e i ), ϕ si,ν = ±ϕ 0ν , ϕ 0ν = arccos(e ⊥iν · e ⊥sν ).
The angle θ s,ν (α ν , ϕ si,ν ) can also be determined from (6). The last equalities in (10) ensure the symmetry of the scattering pattern by the cylinders about the incidence plane. The actual radiation scattering goes in the direction of the vectors e sν . The spherical angles θ sν , ϕ sν corresponding to these vectors can be easily calculated given the Cartesian coordinates of vectors (8). For example, for the x-cylinder we obtain Taking into account (10), we can write for the ν-cylinder (ν =x, y) Here and hereinafter, we use the notation R ν for the radius of the ν-cylinder. Note that (11) yields the probability density depending on the angle ζ ν between the normal e ⊥sν to the scattering plane and the axis Oz. Therefore, (11), as well as the similar relation for the z-cylinder, determines the probability of scattering in the direction θ , ϕ under the condition that scattering is in the plane ζ ν = const; that is, it determines the conditional probability. However, the scattered wave is a conical one (or cylindrical if the radiation is incident perpendicularly to the fiber (strut) (Lind and Greenberg 1966;Bohren and Huffman 1983). It includes the scattering direction vectors corresponding to all admissible values of ζ ν . For that reason, the phase scattering function should establish a relationship between the direction θ, ϕ and all admissible scattering directions e sν , and it must include the components corresponding to the density of probability distribution over the admissible scattering direction vector e sν . It is clear that formula (11) does not satisfy this condition. The desired relation can be obtained using the Bayes formula taking into account the fact that all the values of ζ ν are realized and are, therefore, equally probable. Hence, the complete phase function for ν-cylinder has the following form: The incident radiation is not changed by the part of the representative element that is free of the fragments of the material skeleton (empty part). Therefore, the phase function and the scattering efficiency of this part have the following form: Using the phase functions and optical coefficients of the fragments and the phase function of the empty part of the representative element, we can find the characteristics of the representative element as a whole: Here, the sums must include both the components corresponding to the fragments of the material skeleton and the components corresponding to the empty part of the representative element. The quantities S k⊥ are equal to the area of the normal projections of the fragments on the plane that is orthogonal to the incidence direction of radiation. They must satisfy the equality in which the summation is over all the fragments and the subscripts V b em indicate the representative element as a whole and its empty part, respectively. Since the external boundary of the representative element is a rectangular parallelepiped, not more than three faces are illuminated at any direction of radiation incidence. Taking into account (9), we obtain where the subscript b refers to the juncture of a representative element (if any). Now, the quantity S em⊥ can now be determined from Eqs. (15), (16). In formulas (14) and (16), the following notation is used: a 1,2 are the characteristics of the material anisotropy (Alifanov and Cherepanov 2009;Cherepanov 2011); L 3 is the size of the representative element in the direction of the external heat flux; R ν and l ν are the radius and the length of ν-cylinder, respectively, and α ν is the angle of incidence on this element. In the scalar theory, the problems described above can be solved in a simpler way because one operates with apertures of the scattering bodies rather than with the bodies themselves, and the apertures lie in the plane that is orthogonal to the direction of illumination of the representative element. Almost all the necessary results can be found in Bohren and Huffman (1983). However, in using the corresponding relations one must take into account the fact that they are written in the coordinate system Oxỹz, in which the axis Ox coincides with the aperture axis of the illuminated cylinder because the azimuth is measured from this axis and the radiation is incident on the aperture in the positive direction of the axis Oz. It is clear that this coordinate system is different from the global system Oxyz. Therefore, to use the relations presented in Bohren and Huffman (1983), one needs the following: 1. To know the transformation from the coordinate system Oxyz of the representative element to the system Oxỹz whose aperture axis of the x-cylinder coincides with the axis Ox and the axis Oz is oriented in the direction of radiation e i (9). 2. To know the angles τ xy and τ xz between the apertures of the x, y and z-cylinders.
The angles between the fiber apertures and the unit direction vectors e ν,⊥i can be obtained using the following obvious equalities: Since e i is a unit vector, we obtain which form a right-handed triplet. Then, the transition to this system is determined by the orthogonal (which can be easily verified) transformatioñ The determinant of its matrix is equal to unity. Writing this transformation for the director vectors specified in the spherical coordinates related to the systems Oxỹz and Oxyz, we obtain which makes it possible to determine for each pair of spherical angles θ , ϕ related to the coordinate system Oxyz their analogsθ,φ, related to the system Oxỹz. The anglesθ,φ are arguments of the amplitude scattering function in the scalar theory (Bohren and Huffman 1983) if the scattering by the x-cylinder is considered. In the case of y and z-cylinders, the angleφ in these relations must be replaced with the anglesφ + τ xν , ν = y, z. Unfortunately, these are not all the difficulties. Phase scattering functions (12)-(14) are singular, and they cannot be used in computations. We now show how this difficulty can be overcome in practical computations.

Nonsingular picture of scattering by the representative element
The main idea behind the proposed method is based on the fact that the probability of scattering by the representative element integrally depends on the phase function, and it must continuously relate the directions of illumination and scattering. Therefore, a computational algorithm can be designed that can form such probabilities for a certain set of discrete directions and then, after the renormalization of these probabilities, the corresponding "continuous" values of the phase scattering function can be obtained.
Consider the discrete directions used in S 2m , which is an approximation of the discrete ordinate method (DOM) (see Placido et al. 2005), these directions cover the sphere of directions fairly uniformly because the parts of identical area = π/[m(m + 1)] formed by "parallels and meridians" contain one discrete direction each (hereinafter, the symbol . . . denotes the exact division) Let us arbitrarily choose one scattering direction vector e sν (ν = x, y, z, b) for each fragment of the representative volume element. They are determined by the arbitrarily chosen angles ϕ, ζ ν ∈ [0, 2π], ν = x, y, θ ∈ [0, π] and fixed angles θ i , ϕ i . These vectors are functions of the following angles: for cylindrical fragments, e sν (θ i , ϕ i , ζ ν ), ν = x, y, e sz (θ i , ϕ − ϕ i ); for a juncture (Cherepanov 2011) of the foam material, e sb (θ − θ i , ϕ). When choosing scattering direction vectors, we use discrete directions (17) for the juncture. Denote by [..] the operation of taking the integral part of a number. Then, the discrete directions on the scattering cones of cylinders ζ ν,n = h ζ,ν (n − 0.5), h ζ,ν = 2π n ζ,ν , n = 1 . . . n ζ,ν , have the same density as the directions of grid (17). For example, for the z-cylinder, the discrete directions on the cone for the layer θ sz of grid (17) are spaced with the step h ϕ,z . The scattering direction vector of the empty part of the representative element coincides with the incidence direction vector (4), i. e., e se = e i . In scattering by the representative element, each chosen scattering direction vector of fragments is realized randomly with a probability proportional to the fragment cross section Scattering in each chosen direction occurs with the probability that can be written on grid (17) without using δ-functions, namely |T nb (arccos(e sb , e i ))| 2 , P se = 1.
In the scalar theory, the probabilities of scattering by the fragments are calculated by the formulas from Bohren and Huffman (1983) multiplied by with regard to the relations obtained above.
By averaging the direction vectors of fragment with probabilities (18), we obtain the effective scattering direction vector e V s by the representative element. Due to the assumed independence of scattering by fragments, the vector e V s is associated with the probability equal to the product of probabilities (19) The probability P V s is equivalent to the statistical weight of the direction sample {θ k , ϕ n,k } of grid (17) to whose neighborhood the vector e V s is pointing. Going over all possible values of the discrete elements of the set ζ x , ζ y , ϕ sz , θ sb , where ϕ, ζ ∈ [0, 2π], θ ∈ [0, π], and summing the statistical weights related to the same direction sample, we can assign to each direction sample (17) the accumulated statistical weight P V En,k . Upon its renormalization, it is easy to obtain the probability of scattering P for the discrete directions {θ k , ϕ n,k } and the nonsingular phase function p P(θ k , ϕ n,k |θ i , ϕ i ) = P V n,k / n,k P V n,k , p(θ k , ϕ n,k |θ i , ϕ i ) = P(θ k , ϕ n,k |θ i , ϕ i )/ . (20) 5 Computational experiment

Testing key programs
For the practical implementation of the optical mathematical model described above, software was developed, and its key components were thoroughly tested. Where possible, the available well-proven algorithms and programs were used. For example, in regard to the Mie theory for the cylinder and ball, the Bessel and Hankel functions, as well as their logarithmic derivatives, were calculated using iterative sequences; the reasons for using these sequences and principles of their construction are thoroughly described in Bohren and Huffman (1983). The values of the Bessel and Hankel functions thus obtained were compared with the results provided by the standard programs, which compute those values by summing the appropriate series up to the accuracy much higher than that of the values given in handbooks on special functions (e.g., Loretz et al. 2008). The double precision results produced by different programs coincided up to the decimal digits that could be affected by processor rounding errors.
A standard generator for pseudorandom numbers uniformly distributed on the interval [0, 1] was used. To generate pseudorandom numbers with nonuniform distributions, the wellknown von Neumann method (Öchsner et al. 2008) was used. In all the cases, the generated sequences had the proper first moments when sufficiently ample statistics was gathered.
For the optical characteristics of the ball within the Mie theory, a version of the program presented in Bohren and Huffman (1983) was used. The program was tested by comparing the resultant efficiencies with those published in Bohren and Huffman (1983). When rounded, the results coincided in all the decimal digits with the published results.
For the computations within the Mie theory of the optical characteristics of infinitely long normally illuminated cylinders, the version of the program published in Bohren and Huffman (1983) was also used. However, for the case of arbitrary illumination of the cylinder, new programs were developed, which automatically used the simpler programs from Bohren and Huffman (1983) in the case of near-normal illumination direction. To verify the results produced by the new programs, the following comparisons were made: 1. In the case of the arbitrary illumination direction of the cylinder, the coefficients of the scattering series were compared with those provided in Table 2 in Lind and Greenberg (1966). 2. In the same case, the calculated values of the efficiencies were compared with the values obtained for x = ka = 0.7 and x =3.1 taken from Figs. 10 and 11 in Lind and Greenberg (1966). 3. In the case of the normally illuminated cylinder, the series coefficients, the efficiencies, and the matrix elements were also compared with those presented in Lind and Greenberg (1966). In addition, the efficiencies calculated for x =5.213 were compared with the results presented in (Bohren and Huffman, 1983, p.497). 4. For the near-normally illuminated cylinder, the results obtained by newly developed programs were compared with those obtained using an analog of the program published in Bohren and Huffman (1983). The results produced by the newly developed programs always created a picture of continuous transformation into the results produced by the program from Bohren and Huffman (1983) as the incidence angle varied from an arbitrary one to zero.
Comparison of the computed numerical values of the coefficients of the scattering series with the results presented in Table 2 in Lind and Greenberg (1966) showed that they coincide in all the digits up to rounding. However, the comparison of the computed efficiencies with the results published in Bohren and Huffman (1983) showed that they differ in the fifth decimal digit. On the one hand, this difference is quite understandable because the computations were performed by a different program than that used in Bohren and Huffman (1983) for the normally illuminated cylinder. On the other hand, the efficiencies computed by the program identical to that published in Bohren and Huffman (1983) also differed from the results presented in Bohren and Huffman (1983) in the fifth digit; however, they coincided with the results obtained by the newly created program up to the tenth digit inclusive. This difference can be due to different rounding errors on different processors and to the computational instability of logarithmic derivatives of cylindrical Bessel and Hankel functions (as a matter of fact, the iterative algorithm was proposed in Bohren and Huffman (1983) to suppress this instability). For example, the values of the expression H ( ) computed using the seemingly equivalent relations (Cherepanov 2011;Bohren and Huffman 1983).
can be significantly different. Since the values of the cylindrical Bessel and Hankel functions rapidly decrease at constant arguments as their order increases, only the last of these relations is correct; indeed, in the other relations, the rounding error increases manifold at the division by rapidly decreasing values as the order of the cylindrical functions increases. Figure 4 shows spectral dependences of the absorption coefficients α λ and scattering coefficients β λ computed for a representative element (Cherepanov 2011) of glassy carbon foam ETTI-CF-ERG at the temperature T =500 R illuminated in the direction of the axis Oz. The values of the wave length on the horizontal axis are indicated in meters. The domain in Fig. 4a was scanned with the step λ = 2 · 10 −6 m and in Fig. 4b, with the step λ = 2 · 10 −8 m.

Some modelling results for local spectral properties of highly porous materials
In all the figures, the lower pulsating line, which has narrow and relatively high resonance peaks, corresponds to the spectral absorption coefficient α λ . The pulsating line above shows the spectral scattering coefficient β λ divided by 10 3 . Computations show that the highest absorption resonance peaks are located in the spectral domain that contains the values of the diameters of the juncture and fibers. In the long wave part of the spectrum in Fig. 4a, the width of the resonance peaks increases, while their height (especially in absorption resonances) decreases. Thus, certain coherence of resonance phenomena in absorption and scattering processes occurs, which was observed in spectra of normally illuminated slabs and which is well known (see Bohren and Huffman 1983). In this domain, the representative element (and the material as a whole) behaves as homogeneous medium because the scale of its spatial inhomogeneity becomes insignificant compared with the wave length. As the wave length increases further, the resonance effects become weaker. The absorption and scattering become negligibly weak and the material becomes optically transparent. In the domain of relatively short waves that contains the support of weighting function (21) at the temperature T = 500 K, the resonance effects at absorption are not very significant. In this domain, the variations of spectral coefficients are low amplitude fluctuations about mean values, absorption is insignificant, and the material behaves as conservative medium with invariable properties with respect to the spectrum. The values of such characteristics of the material can be generally obtained by averaging the corresponding values of the representative elements. However, the scattering and absorption spectra of the representative elements of glassy carbon foam in the range of temperatures from 500 to 1000 • C are practically invariable (Cherepanov 2011). However, this is not the case for function (21) whose support shifts to the domain of shorter waves as the temperature increases. For that reason, the accuracy of description based on averaged radiation characteristics of the material under examination improves with increasing temperature. Figure 5a shows the probability of scattering (20) by a representative element of the fibrous thermal protection material TZMK-10. The corresponding phase function integrated over the azimuth is shown in Fig. 5b in polar coordinates. All the results were obtained using the grid of directions S 16 (17). The phase function is very sensitive to changes in the direction of illumination. The relatively high impulse in Fig. 5a corresponds to the part of the incident radiation that was not diffracted, and it smoothly shifts as the direction of illumination changes.
The figure clearly shows the incident and reflected components and the peculiar "diffraction tails" of fibers. As θ i changes, the polar phase function changes drastically. It turns out to be much more complex than the commonly used model phase functions (Colombo 2006). This fact must be taken into account when the radiation transfer in such materials is modeled.
We emphasize that the results presented here are only illustrations of the some capabilities of the developed methods and reflect only a small part of the study of local optical properties of thermal protection materials.

Conclusions
In the framework of the Mie theory and its consequences, a mathematical description of the interaction of electromagnetic radiation with orthogonal representative elements of structure models of lightweight highly porous materials is obtained. Some mathematical difficulties that arise due to the relative geometrical complexity of representative elements and the singularity of one of their key characteristic-the phase scattering function-are overcome. This significantly extends the capabilities of direct Monte Carlo simulation of materials.