2.1. 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
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.
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 Ei and Ej the amounts of energy (in W/m2) 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:
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
EiAiFij =
EjAjFji, 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
Aj does not emit energy and only receives heat from
Ai, 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 Ei successively induces a certain “average” energy E2 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
Ai 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:
2.2. 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
dAj 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
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,
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):
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.
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:
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θ.
Thus, to calculate the area of the elliptical sector we will have to solve:
since, by symmetry,
θ1 =
θ2.
For the larger ellipse, a = 1 m, b = 0.8944 m and θ = tan−1(2/(1.5·0.8944)) = 0.9799 and its complement with respect to π/2 is (π/2) − θ = θL = 0.5908, then the larger ellipse elliptical sector area is: Al = abθL = 1 × 0.8944 × 0.5908 = 0.5284 m2.
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, Am = abθm = 1·0.7071 × 0.4876 = 0.3448 m2.
Finally, subtracting the areas of the previously calculated elliptical sectors and dividing by π, we obtain exactly the same result as in Equation (13), i.e., ((0.5284 − 0.3488)/π) = 0.0584.
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,
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].
2.4. 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.
2.4.1. 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,
Hf) (
Figure 10). Note that the
R and
Hf parameters are provided by two sliders of Grasshopper, so the radius of the fireball and its height above the ground are modifiable at will.
The result of conducting the previous algorithm is shown in
Figure 11.
2.4.2. 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.
2.4.3. 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.
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.
2.4.4. 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).
2.4.5. 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.
2.4.6. Application of the Projected Solid Angle Principle and Determination of the Configuration Factor for an Arbitrary Point on the Surface Partially Affected by a Fireball
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 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:
Fireball radius (R): 5 m.
Fireball height (Hf): 0 m.
Wall height: (Hw): 2 m.
Wall distance (Yw): 10 m.
Ground level (Hg): 0 m.
Horizontal affected area: 40 × 40 m: 1681 points.
The output is presented in the following results section.
2.4.7. Application of the Proposed Method to Find the Excess of Thermal Radiation Due to a Fireball
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
Iw 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.
Iw must be defined as a function of the configuration factors and the emissive power:
wherein
fmax 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, Pw is the partial pressure of vapor (typically 1.15 Pa for a 50% of relative humidity of the air) and df is the distance from the center of the fireball to the control point considered (also automatically provided by our algorithm).
On the other hand,
Ep 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
Hc 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.