Prevention of Hazards Induced by a Radiation Fireball through Computational Geometry and Parametric Design

: Radiation ﬁreballs are singular phenomena which involve severe thermal radiation and, consequently, they need to be duly assessed and prevented. Although the radiative heat transfer produced by a sphere is relatively well known, the shadowing measures implemented to control the ﬁreball’s devastating effects have frequently posed a difﬁcult analytical instance, mainly due to its speciﬁc conﬁguration. The objective of this article is to develop a parametric algorithm that provides the exact radiative conﬁguration factors for the most general case in which the ﬁreball is located at any distance and height above the ground, partially hidden by a protective wall over an affected area at different positions with respect to the said ﬁreball. To this aim we use methods based on Computational Geometry and Algorithm-Aided Design; tools that, departing from the projected solid-angle principle, provide exact conﬁguration factors, in all cases, even if they do not present a deﬁnite analytical solution. This implies dealing with spatially curved radiative sources which had not been addressed formerly in the literature due to their mathematical difﬁculties. Adequate application of this method may improve the safety of a signiﬁcant number of facilities and reduce the number casualties among persons exposed to such risks. As a similar radiative problem appears in volcanic explosions; we hope that further extensions of the method can be adapted to the issue with advantage.


Introduction
A Boiling Liquid Expanding Vapor Explosion (BLEVE) is a type of explosion that occurs in tanks that store liquefied and overheated gases under pressure, in which due to rupture or leakage of the tank, the liquid inside boils and is massively incorporated into the expanding vapor. The most frequent cause of this type of explosion is an external combustion that engulfs the pressurized tank, weakens it mechanically, raises the temperature of the liquid contained and increases the pressure inside the tank. There comes a point at which the pressure reaches values that the container cannot withstand, producing a fissure or rupture of the same. This causes a sudden drop in pressure, the spontaneous nucleation process begins, and all the liquid contained changes its state to gas almost instantaneously, increasing its volume hundreds or thousands of times. BLEVE occurs even though the liquid contained is not a flammable product. The combustion of the content will take place as long as the product is combustible and flammable, but this is a second explosion that constitutes a different phenomenon named as Unconfined Vapor Cloud Explosion (UVCE), also known as fireball, and is a consequence of a BLEVE and not part of it.
The formation of a fireball, isolated or as part of a BLEVE, is characterized by the emission of intense thermal radiation, capable of causing lethal and irreversible damage even to people located at significant distances. Although knowledge about this type of accident has improved substantially, it continues to occur, so it is necessary to develop countermeasures aimed at reducing the effects of its obnoxious consequences.
The process of formation and dynamics of fireballs, including their size, duration, elevation, temperature variation and their concentration range and velocities have been studied by numerous authors [1][2][3][4][5][6][7]. Their initial phase is characterized by the typical silhouette of a mushroom, the most widely accepted model for calculating the ensuing radiative configuration factors and their effects consist of a sphere as a representative figure of the so-called fireball, sometimes at ground level and sometimes elevated.
With several caveats, a similar problem appears in some volcanic explosions that require protection for visitors and researchers.
In recent years, ways for modeling heat transfer processes have been registered in the literature [8][9][10]. However, procedures for finding radiative heat transfer through the space are scarce and inaccurate, mainly due to the difficulties posed by configuration factors.
It is true that there are some exact analytical calculations of configuration factors for the particular geometry of certain fragments of spheres [11][12][13][14][15][16], collected mainly in the catalog of Howell and Pinar Mengüç [17]. However, the problem of the effects and, therefore, of the calculation of the configuration factors, from a fireball partially hidden by a protection wall with an arbitrary elevation above the ground has been (to our best knowledge) approached only by J. M. Bonilla [18]. Said sphere affects an area located at any height with respect to its base-level although as the author himself recognizes (p. 97 of his doctoral thesis).
In the case under study, a fireball with an obstacle constituted by a flat wall, the resulting integral equation for the configuration factor F 12 (Equation (1)), a concept developed below in Section 2 (Equations (2)- (6)), that measures the non-dimensional ratio of energy transferred from surface 1 to surface 2 (see Figure 1) remains analytically unsolvable [11,16]. Therefore, it is necessary to resort to numerical methods to obtain an adequate solution.
F 12 = π 0 π 0 R 0 R 0 (r 2 1 * r 2 * sinθ 1 * r 2 2 * r 1 * sinθ 2 )dr 1 dr 2 dθ 1 dθ 2 π * [(r 2 1 + r 2 2 + 2 * r 1 * r 2 * cosθ 1 * cosθ 2 )] 2 (1) Mathematics 2022, 9, x FOR PEER REVIEW 2 of 21 The formation of a fireball, isolated or as part of a BLEVE, is characterized by the emission of intense thermal radiation, capable of causing lethal and irreversible damage even to people located at significant distances. Although knowledge about this type of accident has improved substantially, it continues to occur, so it is necessary to develop countermeasures aimed at reducing the effects of its obnoxious consequences.
The process of formation and dynamics of fireballs, including their size, duration, elevation, temperature variation and their concentration range and velocities have been studied by numerous authors [1][2][3][4][5][6][7]. Their initial phase is characterized by the typical silhouette of a mushroom, the most widely accepted model for calculating the ensuing radiative configuration factors and their effects consist of a sphere as a representative figure of the so-called fireball, sometimes at ground level and sometimes elevated.
With several caveats, a similar problem appears in some volcanic explosions that require protection for visitors and researchers.
In recent years, ways for modeling heat transfer processes have been registered in the literature [8][9][10]. However, procedures for finding radiative heat transfer through the space are scarce and inaccurate, mainly due to the difficulties posed by configuration factors.
It is true that there are some exact analytical calculations of configuration factors for the particular geometry of certain fragments of spheres [11][12][13][14][15][16], collected mainly in the catalog of Howell and Pinar Mengüç [17]. However, the problem of the effects and, therefore, of the calculation of the configuration factors, from a fireball partially hidden by a protection wall with an arbitrary elevation above the ground has been (to our best knowledge) approached only by J. M. Bonilla [18]. Said sphere affects an area located at any height with respect to its base-level although as the author himself recognizes (p. 97 of his doctoral thesis).
In the case under study, a fireball with an obstacle constituted by a flat wall, the resulting integral equation for the configuration factor F12 (Equation (1)), a concept developed below in Section 2 (Equations (2)-(6)), that measures the non-dimensional ratio of energy transferred from surface 1 to surface 2 (see Figure 1) remains analytically unsolvable [11,16]. Therefore, it is necessary to resort to numerical methods to obtain an adequate solution.
Equation (1), describes the complexity of the factor F12 (concept defined below in Equations (2)- (6), to evaluate the radiative exchanges that occur in Figure 1, as such exchanges typically appear in the problems of the fireball that we are bound to study.
Precisely, the object of this article is to develop a parametric algorithm that provides the exact radiative configuration factors for the most general case in which the fireball can Equation (1), describes the complexity of the factor F 12 (concept defined below in Equations (2)-(6), to evaluate the radiative exchanges that occur in Figure 1, as such exchanges typically appear in the problems of the fireball that we are bound to study.
Precisely, the object of this article is to develop a parametric algorithm that provides the exact radiative configuration factors for the most general case in which the fireball can be located at any position over the ground, partially hidden or not, by a protective structure in an affected area surrounding said fireball. The protection measures, possessing a length of several meters and a considerable height, become costly and for that reason they are usually built as vertical walls, However, our procedure could be applied to other shapes with which we are currently experimenting.
In order to solve the complex issues of obstructions in the way to radiation, we would use methods based on Computational Geometry and Algorithm-Aided Design; tools that, originating from the principle of projection of the solid angle, provide exact configuration factors, in all cases, even when they do not present a known analytical solution. Under such procedure, the cases with a definite solution become the results of our work.
In 1928, W. Nußelt, in an article [19] with 54 lines, 5 equations and only one figure drawn by hand, suggested that the configuration factors could be calculated graphically by (1): placing a sphere of unit radius centered on the surface element dA 1 ; (2): obtaining the intersection between the cone with its vertex in the center of said sphere and any planar trapezoidal surface, in any orientation and position that radiates energy (A 2 ) as directrix, and the unit sphere; (3): projecting orthogonally the previously calculated intersection onto the plane containing dA 1 and, finally; (4): dividing by π the area enclosed by this last orthogonal projection. In this way we would obtain the configuration factor f dA 1 −A 2 ( Figure 2). Despite the short length of Nußelt's article, his intuition was taken at face value without proper demonstration, and consequently introduced the possibility of a geometric approach to a difficult analytical problem. It is convenient to stress that Nußelt did not mention non-planar or curved emitting sources. Unlike Nußelt, we needed to investigate precisely this possibility and, starting from his conjecture, we have developed our article by means of algorithm aided design to provide a solution for such complex problems.
Mathematics 2022, 9, x FOR PEER REVIEW 3 of 21 be located at any position over the ground, partially hidden or not, by a protective structure in an affected area surrounding said fireball. The protection measures, possessing a length of several meters and a considerable height, become costly and for that reason they are usually built as vertical walls, However, our procedure could be applied to other shapes with which we are currently experimenting. In order to solve the complex issues of obstructions in the way to radiation, we would use methods based on Computational Geometry and Algorithm-Aided Design; tools that, originating from the principle of projection of the solid angle, provide exact configuration factors, in all cases, even when they do not present a known analytical solution. Under such procedure, the cases with a definite solution become the results of our work.
In 1928, W. Nußelt, in an article [19] with 54 lines, 5 equations and only one figure drawn by hand, suggested that the configuration factors could be calculated graphically by (1): placing a sphere of unit radius centered on the surface element dA1; (2): obtaining the intersection between the cone with its vertex in the center of said sphere and any planar trapezoidal surface, in any orientation and position that radiates energy (A2) as directrix, and the unit sphere; (3): projecting orthogonally the previously calculated intersection onto the plane containing dA1 and, finally; (4): dividing by π the area enclosed by this last orthogonal projection. In this way we would obtain the configuration factor fdA1−A2 ( Figure  2). Despite the short length of Nußelt's article, his intuition was taken at face value without proper demonstration, and consequently introduced the possibility of a geometric approach to a difficult analytical problem. It is convenient to stress that Nußelt did not mention non-planar or curved emitting sources. Unlike Nußelt, we needed to investigate precisely this possibility and, starting from his conjecture, we have developed our article by means of algorithm aided design to provide a solution for such complex problems.

Fundamental Equations: Form Factors and Configuration Factors
Since the fundamental equations governing radiant energy transfer are well established and systematically discussed in most articles dealing with radiant transfer, we will present here in short only those relationships strictly necessary for a complete understanding of our work. For more in-depth references on the subject, which do not address nonplanar radiative emitters, see for example [20][21][22][23]. Nußelt's conjecture about radiative emission of a trapezoid (A,B,C,D) and its respective projection onto a hemisphere of center (X 0 ,Y 0 ,0), radius 1 and Equation (X − X 0 ) 2 + (Y − Y 0 ) 2 + Z 2 = 1 for Z > 0. Points A, B, C and D which define surface A 2 , are coplanar in the plane z = k. Such points must be joined with said center (X 0 ,Y 0 ,0) to obtain the projection.

Fundamental Equations: Form Factors and Configuration Factors
Since the fundamental equations governing radiant energy transfer are well established and systematically discussed in most articles dealing with radiant transfer, we will present here in short only those relationships strictly necessary for a complete understanding of our work. For more in-depth references on the subject, which do not address non-planar radiative emitters, see for example [20][21][22][23].
Suppose two surfaces A i and A j that emit diffuse radiation as shown in Figure 3. We wish to obtain the energy exchange between both surfaces. If, initially, we ignore the details of the transmission, the problem essentially becomes one of determining the amount of energy that leaves one of the surfaces and reaches the other.
Suppose two surfaces Ai and Aj that emit diffuse radiation as shown in Figure 3. We wish to obtain the energy exchange between both surfaces. If, initially, we ignore the details of the transmission, the problem essentially becomes one of determining the amount of energy that leaves one of the surfaces and reaches the other. where θi and θj represent, respectively, the angles formed by the normals to the differential elements of area dAi and dAj and rij is the vector that joins the differential surfaces dAi and dAj.
If we assume for convenience that there are no inter-reflections from one surface to the other, all the incident radiation will be absorbed and the differential flux of energy gives:

cos cos
(2) Being Ei and Ej the amounts of energy (in W/m 2 ) emitted by surfaces i and j.
Next, let us define the view factors (or form factors), Fij as the fraction of energy that radiates from surface i intercepted by surface j. Then the energy that leaves the surface i and reaches the surface j will be EiAiFij and, reciprocally EjAjFji represents the energy that passes from surface j to surface i.
Then Equation (2) could be rewritten in the form: .
(3) Figure 3. Arrangement of surface-source elements. where θ i and θ j represent, respectively, the angles formed by the normals to the differential elements of area dA i and dA j and r ij is the vector that joins the differential surfaces dA i and dA j .
If we assume for convenience that there are no inter-reflections from one surface to the other, all the incident radiation will be absorbed and the differential flux of energy gives: Being E i and E j the amounts of energy (in W/m 2 ) emitted by surfaces i and j.
Next, let us define the view factors (or form factors), F ij as the fraction of energy that radiates from surface i intercepted by surface j. Then the energy that leaves the surface i and reaches the surface j will be E i A i F ij and, reciprocally E j A j F ji represents the energy that passes from surface j to surface i.
Then Equation (2) could be rewritten in the form: Let us now consider a stationary phase of the UCVE with constant dimensions of the radiant sphere, under which there is no increase or decrease in flux and that all possible losses due to the mode of transmission are negligible. Thus dφ = 0, which implies that E i A i F ij = E j A j F ji , an expression known as Lambert's reciprocity theorem [24]. If we also assume (since it is a frequent case of heat transfer for UVCEs, the negligible emission of affected bodies) that surface A j does not emit energy and only receives heat from A i , the total flux will be: Thereupon, and finally we get: Which reduces the problem to the calculation of purely geometric expressions. A particular case of this equation for the problem of the fireball was presented in Equation1 where we expressed our concerns about the difficulties of solution that it posed.
However, the form factor that Equation (6) gives is the net radiant heat exchange from surface to surface. A surface that emits a certain amount of energy E i successively induces a certain "average" energy E 2 on an adjoining surface, which tells us very little about the detailed distribution of said energy.
When analyzing the effects of a UVCE, it is necessary to describe the flux function point by point, i.e., the radiant field, and thus be able to evaluate the possible damage to people and buildings. Therefore, if we look for the point to pint distribution of energy that reaches the surface j, instead of the final value averaged over the receiving surface, we can effectively reduce the integral of Equation (6) as A i is condensed in a point. In other words, we only need to perform the inner part of Equation (6). Such integral will consequently be double, not quadruple, resulting in a new expression, also strictly geometric, known as pointwise view factor or configuration factor:

Nußelt's Analogy or Solid Angle Proyection Law
As explained in the introduction, Nußelt described in 1928 a hypothesis to calculate form factors from planar figures, that acquired wide predicament but had not been properly demonstrated since its publication. This proposal is also known as the Solid Angle Projection Law.
His ideas follow the geometric logic explained in the segueing lines. Let us consider the planar surface element dA j and connect all the points of its contour with the center of a unit sphere. In this way we will obtain a conical surface, which in turn delimits an area dω over the sphere (Figure 4). The solid angle under which a surface is seen from a point is defined as the area of the conical projection of said surface over a sphere of unit radius centered at that same point.
The vector e r is the unit vector in the direction r ij and the scalar product (dot product) e r ·dA j = cos θ j dA j represents the projection of the vector dA j in the radial direction e r , that is, cos θ i dA j gives the projection of the element of area dA j over a perpendicular plane to the direction e r . Now, a simple similarity relation between dA j at a distance r ij and the subtended area in the unit sphere, allows us to calculate the solid angle element as:  The vector er is the unit vector in the direction rij and the scalar product (dot product) er•dAj = cos θjdAj represents the projection of the vector dAj in the radial direction er, that is, cos θidAj gives the projection of the element of area dAj over a perpendicular plane to the direction er. Now, a simple similarity relation between dAj at a distance rij and the subtended area in the unit sphere, allows us to calculate the solid angle element as: Then the solid angle under a finite surface ω viewed from the center of the unit sphere will be, cos , and its projection onto a plane that forms an angle θi with er will have a surface (with respect to the area π, of the unit circle): cos , Thus, substituting ω for its value we obtain: Since θi is independent from Aj, the orthogonal projection of ω onto said plane (divided by the area of the unit circle π) yields the dimensionless configuration factor previously obtained (Equation (7)).
Furthermore, taking into account Lambert's reciprocity theorem, any diffusive surface (Ak) adscribed to the cone under which Aj is seen, will have the same configuration factor since, as we have assumed, dϕ = 0 and therefore fdAi−Aj = fdAi−Ak ( Figure 5). An important consequence of the above reasoning is that configuration factor problems, can be reduced in most cases to the calculation of a projected area. Then the solid angle under a finite surface ω viewed from the center of the unit sphere will be, and its projection onto a plane that forms an angle θ i with e r will have a surface (with respect to the area π, of the unit circle): Thus, substituting ω for its value we obtain: Since θ i is independent from A j , the orthogonal projection of ω onto said plane (divided by the area of the unit circle π) yields the dimensionless configuration factor previously obtained (Equation (7)).
Furthermore, taking into account Lambert's reciprocity theorem, any diffusive surface (A k ) adscribed to the cone under which A j is seen, will have the same configuration factor since, as we have assumed, dφ = 0 and therefore f dAi−Aj = f dAi−Ak ( Figure 5). An important consequence of the above reasoning is that configuration factor problems, can be reduced in most cases to the calculation of a projected area.
We would like to outline that the values obtained for the configuration factors through the application of the projected solid-angle principle are not a mere graphic approximation but should be identical to the results of analytical methods That said, demonstrating the former principle by integral calculus becomes, in most cases, extremely difficult [18]. Cabeza-Lainez [20,25] has exposed this hindrance of Nußelt's work for the apparently simple case of the configuration factor of a rectangular emitter as defined in Figure 6b. The particulars of such complex integral follow. We would like to outline that the values obtained for the configuration factors through the application of the projected solid-angle principle are not a mere graphic approximation but should be identical to the results of analytical methods That said, demonstrating the former principle by integral calculus becomes, in most cases, extremely difficult [18]. Cabeza-Lainez [20,25] has exposed this hindrance of Nußelt's work for the apparently simple case of the configuration factor of a rectangular emitter as defined in Figure 6b. The particulars of such complex integral follow. The exact analytical formula obtained by direct integration of the configuration factor [23] is ( Figure 6a): Equation (11) substituting the data consigned in Figure 6b gives a result of:  We would like to outline that the values obtained for the configuration factors through the application of the projected solid-angle principle are not a mere graphic approximation but should be identical to the results of analytical methods That said, demonstrating the former principle by integral calculus becomes, in most cases, extremely difficult [18]. Cabeza-Lainez [20,25] has exposed this hindrance of Nußelt's work for the apparently simple case of the configuration factor of a rectangular emitter as defined in Figure 6b. The particulars of such complex integral follow. The exact analytical formula obtained by direct integration of the configuration factor [23] is ( Figure 6a): Equation (11) substituting the data consigned in Figure 6b gives a result of: The exact analytical formula obtained by direct integration of the configuration factor [23] is ( Figure 6a): Equation (11) substituting the data consigned in Figure 6b gives a result of: On the other hand, and surprisingly enough, if we want to verify these values using strictly the projected solid angle principle, the only viable approach to solve the problem is to use the "generalized polar coordinates" in which: x = aρcosθ; y = bρsinθ, other more usual coordinates produced no positive results for reasons that might be related with the special nature of elliptic integrals. Those reasons will be analyzed in the discussion section.
For the said coordinates, the equation of the ellipse takes the form ρ = 1 and the unit surface element is: dA = abρdρdθ.
Note the singular fact that the θ angles in these coordinates are different for each ellipse as they are colligated to the respective sizes of their major and minor axes. As a consequence, the elliptical sector of the minor ellipse area will be, a = 1 m, b = 0.7071 m and θ = tan −1 (2/(1.5 ×·0.7071)) = 1.0832. Its complementary with respect to π/2 is (π/2) − θ = θ m = 0.4876 therefore, A m = abθ m = 1·0.7071 × 0.4876 = 0.3448 m 2 .
Unfortunately, we spent a considerable amount research time to attain this finding which is only valid for the rectangle. Still, other simple shapes of emitters such as the circle (Figure 7) or the ellipse remain unsolved, by direct planar integration or any other method [25]. Moreover, calculations for three-dimensionally curved emitting sources, such as for instance the fragment of sphere that precisely appears in the fireball, become even more complicated.
is to use the "generalized polar coordinates" in which: x = a ρ cos θ; y = b ρ sin θ, other more usual coordinates produced no positive results for reasons that might be related with the special nature of elliptic integrals. Those reasons will be analyzed in the discussion section.
For the said coordinates, the equation of the ellipse takes the form ρ = 1 and the unit surface element is: dA = abρdρdθ.
Unfortunately, we spent a considerable amount research time to attain this finding which is only valid for the rectangle. Still, other simple shapes of emitters such as the circle (Figure 7) or the ellipse remain unsolved, by direct planar integration or any other method [25]. Moreover, calculations for three-dimensionally curved emitting sources, such as for instance the fragment of sphere that precisely appears in the fireball, become even more complicated. In the case of a circular emitting source displaced from the center of the unit sphere, the intersection with the said sphere has been defined by Cabeza-Lainez as a Tomomi curve [25]. This curve (Figure 8), presents the following parametric coordinates, sin cos  (15)) of an oblique cone with circular basis and a unit radius sphere.
In the case of a circular emitting source displaced from the center of the unit sphere, the intersection with the said sphere has been defined by Cabeza-Lainez as a Tomomi curve [25]. This curve (Figure 8), presents the following parametric coordinates, Its projection on a plane is actually a fourth-degree curve as shown in Figure 9 and can be deducted from Equation (15); therefore, it is not possible to find its enclosed area by direct integration in an exact manner. These findings convinced us of the necessity of establishing other approaches for the issue that will be explained in detail in the ensuing section. The above results for the rectangle and the circle have been thoroughly proven by the authors with the proposed new method and they can be found in reference [25].

Computational Geometry and Algorithm-Aided Design
Amongst the different classifications of widespread methods for the determination of configuration factors, a clear line can be traced between analytical and numerical methods, with an almost general omission of experimental methods, arising in the literature.
Although it is true that at the time that Nußelt proposed the hypothesis that bears his name, the available techniques lacked sufficient accuracy to develop graphic methods (so Its projection on a plane is actually a fourth-degree curve as shown in Figure 9 and can be deducted from Equation (15); therefore, it is not possible to find its enclosed area by direct integration in an exact manner. sin sin (15) cos cos cos sin Its projection on a plane is actually a fourth-degree curve as shown in Figure 9 and can be deducted from Equation (15); therefore, it is not possible to find its enclosed area by direct integration in an exact manner. These findings convinced us of the necessity of establishing other approaches for the issue that will be explained in detail in the ensuing section. The above results for the rectangle and the circle have been thoroughly proven by the authors with the proposed new method and they can be found in reference [25].

Computational Geometry and Algorithm-Aided Design
Amongst the different classifications of widespread methods for the determination of configuration factors, a clear line can be traced between analytical and numerical methods, with an almost general omission of experimental methods, arising in the literature.
Although it is true that at the time that Nußelt proposed the hypothesis that bears his name, the available techniques lacked sufficient accuracy to develop graphic methods (so These findings convinced us of the necessity of establishing other approaches for the issue that will be explained in detail in the ensuing section. The above results for the rectangle and the circle have been thoroughly proven by the authors with the proposed new method and they can be found in reference [25].

Computational Geometry and Algorithm-Aided Design
Amongst the different classifications of widespread methods for the determination of configuration factors, a clear line can be traced between analytical and numerical methods, with an almost general omission of experimental methods, arising in the literature.
Although it is true that at the time that Nußelt proposed the hypothesis that bears his name, the available techniques lacked sufficient accuracy to develop graphic methods (so that these methods fell into disuse). Currently, with the development of Computational Geometry and Algorithm-Aided Design, under programs such as Grasshopper ® from Rhinoceros a recovery of the origins of graphical methods and their preferential use in a myriad of problems has been enforced [26]. We firmly believe that the solid angle principle after Cabeza Lainez' generalization effected in Section 2.2, can be applied to advantage to all kinds of non-planar sources and, consequently, to the problems posed by the fireball, assuming diffuse radiation.

Calculation Algorithm
The ensuing steps necessary to obtain the code that gives the configuration factors of arbitrary points on differential surfaces, arbitrarily oriented and positioned, for a fireball of variable height are described below.

Fireball Position and Radius
First, we would write a simple algorithm that sets and constructs a spherical fireball of radius R (which basically depends, as we will see, of the amount and type of combustible substance) and (for convenience's sake), center at point (0, 0, H f ) ( Figure 10). Note that the R and H f parameters are provided by two sliders of Grasshopper, so the radius of the fireball and its height above the ground are modifiable at will. noceros a recovery of the origins of graphical methods and their preferential use in a my iad of problems has been enforced [26]. We firmly believe that the solid angle princip after Cabeza Lainez' generalization effected in Section 2.2, can be applied to advantage t all kinds of non-planar sources and, consequently, to the problems posed by the fireba assuming diffuse radiation.

Calculation Algorithm
The ensuing steps necessary to obtain the code that gives the configuration factors o arbitrary points on differential surfaces, arbitrarily oriented and positioned, for a fireba of variable height are described below.

Fireball Position and Radius
First, we would write a simple algorithm that sets and constructs a spherical fireba of radius R (which basically depends, as we will see, of the amount and type of combu tible substance) and (for convenience's sake), center at point (0, 0, Hf) ( Figure 10). No that the R and Hf parameters are provided by two sliders of Grasshopper, so the radius o the fireball and its height above the ground are modifiable at will. The result of conducting the previous algorithm is shown in Figure 11. The result of conducting the previous algorithm is shown in Figure 11. The second step (Figure 12), needed for the complete definition of the problem, is to set the protective wall at a certain distance from the fireball (yw) and define its height (Hw) and its extremes. Finally, to give more generality to these types of problems, we have conceived the possibility that the fuel tanks or the vehicles that transport the combustible Figure 11. Fireball of radius R and height above the soil H f .

Position, Height, Extremes of the Protective Wall and Height of the Affected Area
The second step (Figure 12), needed for the complete definition of the problem, is to set the protective wall at a certain distance from the fireball (y w ) and define its height (H w ) and its extremes. Finally, to give more generality to these types of problems, we have conceived the possibility that the fuel tanks or the vehicles that transport the combustible could be located on a hill that overlooks the affected area by means of the H g parameter. Again, all these values are modifiable by the user depending on the specific problem to be addressed. An overview of the parameters involved in this type of situation is shown in Figure 13. Figure 11. Fireball of radius R and height above the soil Hf.

Position, Height, Extremes of the Protective Wall and Height of the Affected Area
The second step (Figure 12), needed for the complete definition of the problem, is to set the protective wall at a certain distance from the fireball (yw) and define its height (Hw) and its extremes. Finally, to give more generality to these types of problems, we have conceived the possibility that the fuel tanks or the vehicles that transport the combustible could be located on a hill that overlooks the affected area by means of the Hg parameter. Again, all these values are modifiable by the user depending on the specific problem to be addressed. An overview of the parameters involved in this type of situation is shown in Figure 13.

Calculation of the Shadow Zone Provided by the Protective Wall
To calculate the shadow area provided by the protective wall (Figure 14), we first calculate the intersection of the fireball with a vertical plane that passes through the center of the sphere and contains the midpoint of the highest line of the wall by means of the brep|plane command, obtaining a circumference perpendicular to the wall enclosed in the sphere.

Calculation of the Shadow Zone Provided by the Protective Wall
To calculate the shadow area provided by the protective wall (Figure 14), we first calculate the intersection of the fireball with a vertical plane that passes through the center of the sphere and contains the midpoint of the highest line of the wall by means of the brep|plane command, obtaining a circumference perpendicular to the wall enclosed in the sphere.
calculate the intersection of the fireball with a vertical plane that passes through the center of the sphere and contains the midpoint of the highest line of the wall by means of the brep|plane command, obtaining a circumference perpendicular to the wall enclosed in the sphere.
Next, we draw the tangent lines to the said circumference from the midpoint of the upper line of the wall and select the highest of the two tangents found (tangent lines function). The intersection of the plane that contains the highest line of the wall and the tangent previously chosen with the ground plane (line|plane command), gives us the y coordinate of the shadow line on the exposed area, which separates the real affected zone from the area of shadow. A graphic representation of the result of the previous process is shown in Figure 15.  Next, we draw the tangent lines to the said circumference from the midpoint of the upper line of the wall and select the highest of the two tangents found (tangent lines function). The intersection of the plane that contains the highest line of the wall and the tangent previously chosen with the ground plane (line|plane command), gives us the y coordinate of the shadow line on the exposed area, which separates the real affected zone from the area of shadow. A graphic representation of the result of the previous process is shown in Figure 15.

Creation of a Grid of Control Points to Calculate Configuration Factors
When systematizing and automating the process, it is convenient to create a mesh of control points on which to calculate the different configuration factors between the spherical cap, resulting from the section of the fireball produced by the shadow plane and that given by the plane perpendicular to the line connecting each point with the center of the

Creation of a Grid of Control Points to Calculate Configuration Factors
When systematizing and automating the process, it is convenient to create a mesh of control points on which to calculate the different configuration factors between the spherical cap, resulting from the section of the fireball produced by the shadow plane and that given by the plane perpendicular to the line connecting each point with the center of the fireball. As can be seen in the code depicted in Figure 16, the affected surface to be evaluated has the same width as the protection wall and its length can be freely modified by means of a slider. We also want to note that the first row of control points have been omitted, since in that entire row the configuration factors are null by definition (Figure 17).

Creation of a Grid of Control Points to Calculate Configuration Factors
When systematizing and automating the process, it is convenient to create a mesh control points on which to calculate the different configuration factors between the sph ical cap, resulting from the section of the fireball produced by the shadow plane and t given by the plane perpendicular to the line connecting each point with the center of fireball. As can be seen in the code depicted in Figure 16, the affected surface to be eva ated has the same width as the protection wall and its length can be freely modified means of a slider. We also want to note that the first row of control points have been om ted, since in that entire row the configuration factors are null by definition (Figure 17)

Determination of the Radiant Heat Cones with Shadow
We will now proceed to calculate the radiant heat cones with vertex at the control point tangent to the fireball, subtracting from them the part that lies below the protective wall. To perform this, given that the radiant heat cones and the fireball are two quadrics with two planes of common symmetry-one of them vertical that passes through the center of the fireball and through the control point-we would cut the fireball by the said plane to obtain a circumference contained in the fireball and in the aforementioned vertical plane. Tracing the tangents to above-defined circumference from the control points and forming the surface of revolution that generates by rotating a 2π angle, one of the previous tangents (generatrix), around the line that joins the control point and the center of the fireball (axis), (by means of the revolution function) we obtain the complete radiant

Determination of the Radiant Heat Cones with Shadow
We will now proceed to calculate the radiant heat cones with vertex at the control point tangent to the fireball, subtracting from them the part that lies below the protective wall. To perform this, given that the radiant heat cones and the fireball are two quadrics with two planes of common symmetry-one of them vertical that passes through the center of the fireball and through the control point-we would cut the fireball by the said plane to obtain a circumference contained in the fireball and in the aforementioned vertical plane. Tracing the tangents to above-defined circumference from the control points and forming the surface of revolution that generates by rotating a 2π angle, one of the previous tangents (generatrix), around the line that joins the control point and the center of the fireball (axis), (by means of the revolution function) we obtain the complete radiant heat cone.
Next, we must cut the newly calculated cone and the fireball along the plane defined by the control point and the upper line of the protection wall (Brep|Plane command) by means of the Surface Split function, keeping the upper part of the three surfaces found by this function (Figure 18). The result is shown in Figure 19a. point tangent to the fireball, subtracting from them the part that lies below the protect wall. To perform this, given that the radiant heat cones and the fireball are two quadr with two planes of common symmetry-one of them vertical that passes through the ce ter of the fireball and through the control point-we would cut the fireball by the sa plane to obtain a circumference contained in the fireball and in the aforementioned ve cal plane. Tracing the tangents to above-defined circumference from the control poi and forming the surface of revolution that generates by rotating a 2π angle, one of previous tangents (generatrix), around the line that joins the control point and the cen of the fireball (axis), (by means of the revolution function) we obtain the complete radia heat cone.
Next, we must cut the newly calculated cone and the fireball along the plane defin by the control point and the upper line of the protection wall (Brep|Plane command) means of the Surface Split function, keeping the upper part of the three surfaces found this function (Figure 18). The result is shown in Figure 19a.  In order to achieve this, we would place a sphere of unit radius at the chosen point and obtain the intersection between the sector of the cone calculated in the previous section with the said sphere, obtaining an arc of circumference over the unit sphere ( Figure  19b). To close the solid angle, we need to compute the intersection circumference between the unit sphere and the plane that yielded the radiating cone sector above the wall.
Of this circumference we are only interested in the smallest portion that is included In order to achieve this, we would place a sphere of unit radius at the chosen point and obtain the intersection between the sector of the cone calculated in the previous section with the said sphere, obtaining an arc of circumference over the unit sphere (Figure 19b).
To close the solid angle, we need to compute the intersection circumference between the unit sphere and the plane that yielded the radiating cone sector above the wall.
Of this circumference we are only interested in the smallest portion that is included between the initial and final points (End Points) of the intersection formed by the radiating cone sector and the unit sphere. To find the value of such piece, we take the midpoint of the segment that links the initial and final points and determine the location of the closest point of the unit sphere (Curve Closest Point). Thus, we will arrive to the arc that passes through the ends of the intersection circumference segment, between the unit sphere and the radiating cone sector and its nearest point.
Finally, it remains to connect the arcs previously calculated (Join Curves) and draw their orthogonal projection on any desired plane (vertical, horizontal or maximum exposure situation). In this case we have selected the orthogonal projection on the affected ground plane.
To complete the problem, Grasshopper calculates the area of said projection on the ground plane and divides it by π, yielding the configuration factor of the arbitrary point chosen ( Figure 20). By automating the process through a loop for all the points of the control gri obtain a matrix of n × n elements that constitute each configuration factor of the aff zone considered, in this case on a meter by meter basis. The parametric nature of the developed, allows us to obtain them for any distance and plane (horizontal, verti maximum).
In the following example, the parameters that have been considered are: Returning to the initial problem of radiative heat transfer, we must reconside physical parameters on which the heat exchange depends in order to control the By automating the process through a loop for all the points of the control grid, we obtain a matrix of n × n elements that constitute each configuration factor of the affected zone considered, in this case on a meter by meter basis. The parametric nature of the code developed, allows us to obtain them for any distance and plane (horizontal, vertical or maximum).
In the following example, the parameters that have been considered are: Returning to the initial problem of radiative heat transfer, we must reconsider the physical parameters on which the heat exchange depends in order to control the likelihood of a fireball inducing damages on the personnel present at the facilities.
To this aim we will use an empirical formula that according to [27,28] gives Thermal Radiation Dose (DTU) that produces a small probability of lethality for an average population, with second degree burns (1% lethality), estimated at 1000 DTU: where I w is the intensity of the radiation in kW/m 2 and, is the time required for the fireball to reach its maximum size and M stands for the mass of fuel. I w must be defined as a function of the configuration factors and the emissive power: wherein f max is the maximum configuration factor at the point considered (given directly by the algorithm proposed in the methods section), and we need to complete the expression with, Which gives the atmospheric transmissivity, P w is the partial pressure of vapor (typically 1.15 Pa for a 50% of relative humidity of the air) and d f is the distance from the center of the fireball to the control point considered (also automatically provided by our algorithm).
On the other hand, E p is the emissive power and is better estimated by the equation: where, as expressed in references [27,28], η rad is the fraction of the fuel mass involved in the explosion (typically equal to 0.25). M, as we have mentioned, is the fuel mass and H c is the calorific power of said fuel and D is provided by the expression: A calculation example with the relevant quantities is detailed in the results chapter.

Results
Following what has been explained in the previous section about our methods, we have been able to obtain an innovative representation ( Figure 21) of the horizontal configuration factors (i.e., the factors projected onto the horizontal plane) for a shadowed fireball. All relevant parameters of the figure are mentioned in the caption.
In this way we can determine the radiative intensities due to a shadowed fireball of variable dimensions and positions and, at the same time, provide the disaster prevention department with an effective design tool. In Figure 22, (vide infra) we have compared these results with the case of the unobstructed fireball [25].

Applied Mathematics. Safety Checks
Although the most complex issue when calculating the damage that a fireball can produce is to obtain the configuration factors, we must not forget that it is still an issue of radiative heat transfer. Thus, we must take into account the physical magnitudes on which said transfer depends to estimate the hazards that a fireball may cause on personnel present at the site.
Using equations 16 to 21, as previously detailed, we have proceeded with the following case study. We would suppose a tank located in an airport at ground level loaded with 40,000 L of the most widely used kerosene in commercial aviation, Jet A-1, with a density of 0.8 kg/L and a calorific power H c = 42,800 kJ/kg. We will thus have a mass of fuel M = 32,800 kg that would cause a fireball (assumed at ground level) of an estimated diameter D of = 6.14 × 32,800 (0.325) ∼ = 179 m.

Results
Following what has been explained in the previous section about our methods, we have been able to obtain an innovative representation (Figure 21) of the horizontal configuration factors (i.e., the factors projected onto the horizontal plane) for a shadowed fireball. All relevant parameters of the figure are mentioned in the caption.
In this way we can determine the radiative intensities due to a shadowed fireball of variable dimensions and positions and, at the same time, provide the disaster prevention department with an effective design tool. In Figure 22, (vide infra) we have compared these results with the case of the unobstructed fireball [25].

Applied Mathematics. Safety Checks
Although the most complex issue when calculating the damage that a fireball can produce is to obtain the configuration factors, we must not forget that it is still an issue of radiative heat transfer. Thus, we must take into account the physical magnitudes on which said transfer depends to estimate the hazards that a fireball may cause on personnel present at the site.
Using equations 16 to 21, as previously detailed, we have proceeded with the following case study. We would suppose a tank located in an airport at ground level loaded with 40,000 L of the most widely used kerosene in commercial aviation, Jet A-1, with a density of 0.8 kg/L and a calorific power Hc = 42,800 kJ/kg. We will thus have a mass of fuel M = 32,800 kg that would cause a fireball (assumed at ground level) of an estimated diameter D of = 6.14 × 32,800 (0.325) ≅ 179 m.
Imagine also that a protective wall is located at 100 m from the kerosene tank and presents a height of 6 m above the ground. Substituting in Equation (18) t and Iw by the values provided by Equations (18)-(21), through processing and reworking the expressions found on Ref. [15], we would obtain a DTU of:  Imagine also that a protective wall is located at 100 m from the kerosene tank and presents a height of 6 m above the ground. Substituting in Equation (18) t and I w by the values provided by Equations (18)-(21), through processing and reworking the expressions found on Ref. [15], we would obtain a DTU of: from where we will find a safety distance with a lethality of less than 1% (D tw = 996.8252 (kW/m 2 ) 4/3 s) of 105.75 m measured from the protection wall, with a maximum configuration factor of f max = 0.1528 (as provided by our method) at d f = 224.17021 m.

Discussion
The radiant field caused by a sphere which constitutes the fireball in its most usual appearance cannot be determined by analytical methods in the presence of defensive frames that obstruct the radiation. The fragments of sphere that induce the emitting surface are of an irregular nature not suited to conventional integral calculus.
The main reason for this difficulty is, in the first place, uncertainty on the geometrical boundary of the source that can be no longer taken as a sort of disk but as complex circular or spherical fragment. In the second place, the variability of the said emitter depends on the receiving point considered. In other words, the geometry of the source is not fixed in space.
To overcome such major drawbacks, we have resorted to the principle of the projected solid angle, which was suggested by Nuβelt but was not properly demonstrated by means of area integrals. In fact, in reference [21] it is stated literally that "the integration (of solid angle projection) is generally considerably less involved than the one with general equations". Unfortunately, we believe that such redundant adverbs suggest more wishful thinking than a fact responding to reality for most cases of configuration factors.
Thus, we have extended the validity of the principle, first by providing an analytical proof for the case of a simple rectangular source.
In this common case for instance, polar coordinates, such as x = ρcosθ; y = ρsinθ, are not well suited to the singularity of elliptic fragments. Initially they were employed for the problem since they offer the undoubted advantage that the parameter angle θ is the same for all the figures that we can represent with such coordinates. On the contrary, with generalized polar coordinates, the angle θ, must be defined for each separate figure or ellipse in this case, as it is not constant or depends on external values. This seemed a strange fact and enough reason to discard it in the initial steps, although it proved essential as we would demonstrate below. We believe that such a situation might be related to the singularities of elliptic geometry which led, for instance, to open questions about elliptic integrals or the volume of a non-revolution ellipsoid.
Subsequently, we have applied the initial steps of our proposed algorithm-aided method to a variety of configuration factors for planar shapes as described in reference [25]. In a similar fashion, and to complete the study at the demand of the editors, the unobstructed radiant field due to a fireball has been compared with our method through Figure 22, showing a total coincidence.
The main novelty of this article is that in it we have been able to extend for the first time our method to 3D spherical irregular sources and we have added the shadow systems and obstructing frames to our calculations. We believe that such findings will mitigate the casualties experienced in this kind of accident.
In this manner, as the proposed computational algorithm is parametric, we could evaluate the distance and height of the protective wall or structure to be built, for any given fuel mass that causes a previously established DTU. In addition, the effect of several fireballs acting simultaneously can be evaluated with a similar procedure. In the future we hope to be able to evaluate the dynamic performance of the fireball phenomenon and to develop some experimental validation methods for these findings.

Conclusions
This article proposes and develops a new graphic algorithm that gives way to new codes useful for calculating the configuration factors produced by a fireball of any diameter and height obstructed by a protective wall in an innovative and accurate manner.
The viability of the method resides firstly in the novel extension of the project solidangle principle to non-planar emitters, as shown in Section 2.2, and secondly in identifying the hidden geometric features of the surfaces involved.
The radiant shapes in these cases are evolving in a rather unpredictable way and could not be addressed analytically.
Our procedure solves the configuration factors exactly, as was previously shown, even in the presence of dynamically evolving fireballs. The simpler cases of the phe-nomenon have been evaluated with our method showing a total correspondence that we are convinced that will bring us the momentum for further attempts of extension.
In a manner of conclusion our manuscript presents an innovative calculation for an exacting problem that induces many casualties and causes much economic damage each year at a worldwide scale. Besides, as the thermal radiation dose received by each point of the spatial zone defined is found and deftly determined, critical aspects of the design of protections of various types can be implemented, enhancing the safety and well-being of the personnel and industries concerned in the near future by virtue of computer aided numerical-geometric methods.