New Computational Geometry Methods Applied to Solve Complex Problems of Radiative Transfer

: Diverse problems of radiative transfer remain yet unsolved due to the difficulties of the calculations involved in them, especially if the intervening shapes are geometrically complex. The main goal of our investigation in this domain, is to convert the equations that were previously derived, into a graphical interface based on the projected solid-angle principle. Such procedure is now feasible by virtue of several widely diffused programs for Algorithms Aided Design (AAD). Accuracy and reliability of the process is controlled in the basic examples by means of subroutines from the analytical software DianaX, developed at an earlier stage by the authors, though mainly oriented for closed cuboidal or curved volumes. With this innovative approach the often cumbersome calculation procedure of lighting, thermal or even acoustic energy exchange can be simplified and made available for the neophyte, with the undeniable advantage of reduced computer time.


Form factors
To our knowledge, one of the main problems in Science as applied to aerospace, solar and industrial design technology has been to determine how the existing physical fields are transformed according to the physiognomy of fixtures or products and in which direction we should orient our developments to seek a more correct transmission of the environmental and heat transfer phenomena, or in other words, in which way the design of three-dimensional forms can be improved to obtain an optimal and coherent distribution of energies that effectively contributes to mitigate global warming and thereby climatic issues.
We must outline that radiation in a physical or spatial environment is made manifest through fields of a fundamentally vector nature. Therefore, our first objective would lie in the assessment of these types of fields in an unaltered state. On such issues, there exist relevant contributions of Yamauchi [1] and Moon [2,3] amongst others to seek radiative potential. However, successive modifications of the spatial features of the elements entailed, possess the capacity to substantially alter the layout of the field. Such geometric alterations contribute to changes in energy distribution that manufacturers and users are heavily demanding as it is already vital for them to acquire accurate and simple notions on the issue.
The problem is not recent at all. In fact, Lambert's statement of the famous theorem XVI of his treatise Photometry [4], wondered about the amount of rays (flux) that issues from any two equally luminous surfaces onto the adjacent and established that if the two surfaces A1 and A2 were equally luminous and faced each other in some way, the amount of incident rays from either of the two surfaces on the second is identical. The former is also known as reciprocity theorem and to develop it in some extent is crucial to assess the mathematical reach of the matter that we are discussing.
The energy that leaves surface A1 and reaches surface A2 will be: and, reciprocally, the energy that passes from surface A2 to A1 is: where E1 and E2 are the equivalent amounts of energy (in W/m2) emitted by surfaces A1 and A2; F12 or F21 are dimensionless entities called "form factors". If we assume, in principle, that there are no inter-reflections or re-emissions from one surface to the reciprocal, and theoretically no other significant sources of radiation have access, by any means, to the boundaries of the problem, all incident flux will be absorbed and the flow of energy (dΦ), according to Lambert's theorem should be zero, that is: To this aim, we would consider the surface elements dA1 and dA2. The angles θ1 and θ2 are measured between the line from dA1 to dA2 (direction of propagation) and the normal to each surface; r is the distance from the center of one differential area to the other ( Figure.  The cosine emission law, also attributed to J. H. Lambert, states that the flux per unit of solid angle in a given direction θ (or radiant intensity) is equal to the flux in the normal direction to the surface (maximum intensity) multiplied by the cosine of the angle between the normal and the direction considered.
Luminance or radiance is the radiant flux or power emitted per unit area and solid angle (steradian) in a given direction. To obtain the flow emitted by the surface element dA1 it is necessary to multiply the radiance or luminance by the projected surface in the fixed direction of the angle θ1, and thus: where I, is the radiant intensity, dΩ is the solid angle and L1 represents the radiance of surface A1.
It is easily deducted that L1 remains independent of the direction and therefore L1θ=L1.
Since the previously found flux refers to the solid angle unit, we must now evaluate the amount of flux that reaches a differential surface element dAn whose distance to A1 is precisely r.
It is often useful to establish a relationship between luminance or radiance and lighting or irradiance, to this aim, we could simply extend the flux to a hemisphere of radius r in which dA1, is inscribed. The surface dAn in spherical co-ordinates amounts to: dA n = r 2 sin θ 1 dθdφ, then from 5 and 6: Thus, E1=πL1 (E is emitted radiation per area unit, expressed in W/m 2 ).
Returning to the problem of energy exchange, if we assimilate the differential surface element dAn with dA2, we would find that: In this manner, the flux or radiant power that goes from dA1 to dA2, substituting L1 in eq. (5) will be: And respectively, dΦ 2−1 = E 2 cos θ 2 cos θ 1 dA 2 dA 1 πr 2 .
The energy exchange is eventually set at: The above integral equation responds to the fundamental or canonical formula of radiation which in differential terms is: In such situation, the integral found by virtue of the fundamental expression, adopts the value of A1F12 or A2F21, following the above definitions. In general, these are quadruple integrals and to solve them becomes so complex that extreme dexterity in geometry problems joined by haphazard techniques of integration is required. Therefore, we want to avoid this difficulty as much as possible by virtue of the procedure presented in this article.
The fundamental equation of the much-sought form factors is: and its asymmetric reciprocal in algebraic terms, ]. (14) If we assume for easier visualization and since it is a frequent case in luminous radiative transfer (but not so in thermal problems involving for instance gases), that the surface A2 does not emit energy and only receives it from A1, the total flow received in A2 will be: where M1= N2/F21 and N2 is the energy received by surface 2 in W/m2. radiation "map" and to be able to control the possible deficit or surplus which could, in turn, become critical.
Therefore, if we look for the point distribution of power that reaches surface A2 instead of the mere average, we could reduce the fourfold integrals in the expression of the form factor to a simpler surface integral as follows: will be called configuration factor and it presents several fundamental applications in the method that we will develop hereby.
The first elementary consequences are: and Then, Accordingly, Eq. (21) can be clearly identified with the "average" of f12"over surface" A2, especially since we had already obtained that N2 = F21M1. In quantum dynamics terms, the form factor expresses the probability of a surface to be stricken with photonic energy emitted by an adjacent source.
And similarly we would define N2 as the average of n2 on the surface A2; This is important finding represents the ignored relationship between form factors and configuration factors that will allow us to easily simplify and visualize many complex calculations in the following chapters by restricting them only to second order primitives and the ensuing "average" is ready accessible by numerical procedures. To our humble knowledge, it has solely been defined by Cabeza-Lainez [5]. United triangles (A1 and A2) and a point belonging to a parallel plane.

Rectangle
The full rectangle is just the sum of solutions for surfaces A1 (24) and A2 (see However, we need to remind that we have just obtained the value of radiative exchange at a single point located precisely under the vertex of the rectangle and in a particular direction (perpendicular) to the surface source.
Deft combinations of triangles, would give us the equilateral triangle, hexagon and as a limit case that we will discuss in Annex 1 the circle, which has been only obtained in such "polygonal" fashion by Cabeza-Lainez.

Calculations in a plane perpendicular to the figure
In the above discussion, we have solely found the component of radiation that was perpendicular to the emitting surface at an isolated point.
It is easy to imagine the severe complexity that arises when trying to find the solution for points moving freely in a plane, let alone space.
But even so, being radiance a vector, to complete the problem means finding values of other auxiliary components.
In this situation, the principle integral turns out to be (if we multiply by the corresponding By virtue of the former, we have exposed how complex it turns to find an analytical expression to solve the problem if only we alter the coordinate plane of calculation not to mention figures of geometry other than those limited by straight lines. That is why we concluded that there was heavy need for a graphical solution based on parametric design as this would imply a universal solution and an extraordinary simplification of the issue.

Calculation of integral equation for a 'rectangle' over a horizontal sqaure
Previously (Eq. 27) we had found the solution in terms of arctangent for the entire rectangle. The  For example, in the Y axis and at level -H1, the adapted solution would give:

Solid Angle Projection Law
To obtain the energy exchange between a given surface and a point P, we would simply trace a cone whose vertex is the point considered P, and its base is the surface in question (dΩ), and proceed to intersect the said cone with a sphere of unit radius (figure 5). The enclosure within this intersection, projected on any reference plane that we may require and divided by the horizontal area of the aforementioned unit sphere (that is, divided by π), will be equate the value of the configuration factor determined by virtue of the analytical methods of exact integration described in the previous chapters.
If, by means of what we have discussed in detail in the above sections, we were trying to find certain geometric factors or proportions, it seems equally reasonable to achieve these same values through graphic procedures such as the ones employed in Descriptive Geometry.
The advantages are obvious to the project engineer, because if, due to the difficulty of the mathematical problem, hesitation appears, we are allowed to somehow 'visualize' the entire process.
In summary, what has been proven is that the configuration factor is a mere proportion, a dimensionless number, consisting of a ratio of areas, one of them the area of a certain projection and the other the total area involved expressed in circular terms, usually π.
Such area divided by πr 2 (r = 1) equates the configuration factor.
The most important consequence of all the former is that the problem, regardless of its complexity, will always present a unique non-trivial solution, since the projection will have a univocal value and cannot be non-existent.
In addition, as the configuration factors are actually projections, they hold the additive property, which is useful to know when dealing with sundry and simultaneous emitting surfaces, because it is possible to add their effects. The average over the receiving surface of the configuration factors thus obtained, will constitute the form factor that had been so laboriously sought-for by means of mathematical formulation and solving of fourfold integral equations.
Alternatively, graphic methods and analytical methods or a combination thereof can be used.
This is the gist of our procedure.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 December 2020 doi:10.20944/preprints202012.0084.v1 It also seems gratifying to assume that we have solved the fundamental problem of radiation with geometric methods. Why is that? Because it means that 'form' is very important in this type of problem and that not all forms act in the same way from a radiative point of view, some geometries favor exchanges more than others and the mission of the heat transfer analyst is to adequate these forms predicting their impact. Such repercussion can also be verified by observers, that is, manufacturers, society, etc., using, as we shall see, a simple tool.
But the sole hindrance that we had identified at the time, was that numerical validation of the simple rectangle over a horizontal plane (Fig. 4) was not readily available, thus casting shadows over the viability of our approach. Fortunately, Cabeza-Lainez was able to work out an exact solution and such relevant finding is presented in Annex 2.
Then a similar problem arose with circular emitters in a free position, Cabeza-Lainez has completely solved it anew, as described in Annex 3 as it may help to validate other curved geometries.

Algorithms Aided Design
Although today, most Computer Aided Design (CAD) programs incorporate what is known as parametric design (which implies that it is not necessary to redraw the objects that we design if they change its shape, position or size), we should not consider all of them as Algorithm Aided Design programs (AAD).
An algorithm is a logical expression that leads to procedures to be performed through finite sets of closed basic instructions. Algorithms simulate human abilities to break down a problem into a series of simple steps that can be easily executed and, although they are closely related to computers, algorithms can be defined independently of programming languages and applications.
However, not just any procedure can be properly considered as an algorithm, since many procedures are far from being well defined or contain ambiguities, logical contradictions or infinite loops.
Some important properties that are required for a procedure to be considered an algorithm are: • An algorithm is a well-defined set of properly concise instructions.
• An algorithm demands a defined set of inputs that may or may not come from the output of a previous algorithm.
• An algorithm generates a precise output ( Figure 6): Figure 6. Algorithm programmed in Grasshopper® to obtain all the straight lines.
• And finally, an algorithm can produce warnings and error messages through the appropriate editor. If the inputs are not appropriate, e. g. if we enter text instead of numbers, the algorithm will return an error message instead of the expected output with the appropriate editor [10].
In addition, through parameterization, we have not merely obtained a straight line, (

Calculation by Algorithms Aided Design through Finite Element Method
The main problem of analytical calculation is that, although we may employ a convenient algebra of form factors [6] to find the radiance received by a surface, issuing from basic figures limited by straight lines, still the number of shapes of emitters whose direct integration is very difficult to perform or impossible, remains numberless.
To solve such an impasse, we have successfully developed a Grasshopper ® algorithm that allows us to calculate the configuration factor subtended by an arbitrary point of any surface (i. e. irrespective of its shape and orientation) from an emitter or 'three-dimensional figure', with manifold geometry or inclination (including those whose integration is deemed impossible).
To calculate the incident radiation n (W/m2), we will only have to add up all the configuration factors that we may require and multiply them by the luminance that departs from the emitter (W/sr·m 2 ).
The steps which we have followed to develop the algorithm can be summarized below: First we establish the plane of the irradiated surface defining, by means of sliders, the components of its normal vector and using the 'Vector XYZ' function, constructing the said plane with the help of the 'Plane Normal' function as shown in the upper part of the following figure (Fig.   7).
Subsequently, we divide the irradiated surface into N×(N-1) squares (finite elements) using the 'Square' function that generates a bi-dimensional grid with square cells and retains the points at the corners of the squares.
Note that (in this case) the first row of corner points has been cancelled, since in the plane of the 'rectangle', the radiance will necessarily be equal to zero and the size of the squares is divided by its number for a better performance of the algorithm (Figures 12, 13):

Exporting the data to Excel and representing them by color maps
To transfer the data to Excel ® , we first sort them using the 'Partition List' command and then the 'GhExcel' plugin that allows us to export an array of data to the said spreadsheet through the 'ExcelWrite' command, (Figure 11). The outcome for such an algorithm is shown below (Figure 13).

Numerical results
When exporting the values obtained by the Grasshopper® algorithm to an Excel® sheet, it is relatively easy to sort them in a matrix and the figures obtained are displayed in Table 1. In this manner, the form factor F12 is obtained as the average of all the configuration factors f12, which yields a value (for a discretization of 11×11=110 points) of F12 = 0.00556888.
If we consider the point with the highest configuration factor (point 5,3), where f12 = 0.01144 and we want to know the radiation of that point, assuming a radiant power in the rectangle of 1000 W/m2, we will obtain it by simply multiplying the value of f12 by 1000, that is, n = 11.44 W/m2.
The array presented in Table 1 coincides with the one that we can find by the analytical procedures ( Table 2) described in the references [6]. In this way the procedure is validated

Convergence
To test the convergence of the method we need to gradually decrease the size of the elements and show how by doing so, the form factors converge towards a certain value which must be increasingly closer to that calculated analytically.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 December 2020 doi:10.20944/preprints202012.0084.v1 One of the advantages of our algorithm is that simply modifying the value of the sliders to increase the partitions we obtain the new results rapidly and proof of convergence is readily available. Figure 14. The color map of the configuration factor obtained by discretization for each 10 centimeters (10201 values) over an area of 10×10 m 2 radiated by a 'rectangle' of 2×2 m 2 , 5 meters over the horizontal surface at its middle axis.
In the previous figure (Fig. 14), we found the results by raising the number of partitions or finite elements from 10×10 = 100 elements to 100×100 = 10000 elements, with a computing time of only 136 s (conventional computer).

Versatility of the method
As noted above, the difficulties in obtaining analytical solutions for emitters with forms other than the trivial ones are enormous, and even in many cases impossible. However, the method presented here allows to obtain solutions with high accuracy, limited only by the performance of the computer, for any type of emitters and irrespective of the orientation.
As a more complex example, we have calculated the configuration factors on a 10x10 m horizontal surface, radiated by an elliptical emitter whose major and minor axes are 2 m and 1 m respectively and whose lowest point is located at a height of 5m. As before, vertical axis of the ellipse coincides with the middle point of the horizontal surface in 5m ( Figure 16). Notice that, for the time being, the elliptical shape cannot be solved analytically due the appearance of elliptic integrals in the process [8,9].  Notice that in the case of the ellipse and the scalene triangle unlike the rectangle, it is not possible nowadays to calculate the configuration factors by any other procedure.

Conclusions and future aims
Given the accurate convergence of the values obtained for the particular case of rectangular emitters and receivers within a cuboid, we humbly consider that the graphical interface program that we have proposed is fully accomplished. However, for the time being, we have to warn that there is no feasible way to validate the other figures presented such as ellipse and scalene triangle.
Nevertheless, we believe that a complex physic-mathematical problem which has occupied the mind of scholars for decades if not centuries, is solved with perfect ease in a semi-automatic manner.
Ensuing researches will extend the procedure to other types of three-dimensional spaces whence analytical calculations were proposed but still remain underdeveloped, for instance curvilinear sources in non-cuboidal domains affected by interferences and obstacles.
A much wider repertoire of possible radiative exchanges would be achieved in this fashion, with relative simplicity and accuracy and without alternative. The repercussions of such findings would be wide-ranging in fields as diverse as aerospace technologies, Light Emitting Diodes (LED), radiotherapy medicine, risk prevention, acoustics and Architecture. On an industrial basis, we consider that the vast field of artificial lighting and design of luminaries will be immediately benefited by our methods.
Whereas the findings hereby exposed would represent a major achievement of Optics, precisely for the realms of Heat Transfer and Lighting, it is understood that our contribution has been timely supplied by the enormous potential of Algorithms Aided Design in the numerical resolution of physical-mathematical problems that, when approached directly, are extremely difficult to accomplish (if not altogether impossible).

A.1. Triangles and shapes composed of triangles
In the first case, which is also the most elementary, we would substitute the value of the cosine as expressed in Eq. 16 and find the fundamental double-integral equation, Note that we write again E (surface-to-point received power) instead of the previous magnitude n because this nomenclature is more usual. We have chosen an instant of time t in which the luminance or radiance L can be considered constant, keeping in mind that Lπ =M.
We will solve it as an example of the general ensuing procedure to attest to its difficulty.
The upper line of triangle A1 (Fig. A1)  Imagine if you can the severe complexity that arises when trying to find the solution for points moving freely in space, it is almost unfathomable.
But even so, as it is logical, radiance should be treated like a vector, thus, if we want to complete the problem we will have to find also the values of the two components that are missing and remain parallel to the area considered. The defining integral then turns out to be (if we multiply by the Adding and subtracting rectangles would allow us to obtain the result for a rectangular shape over an unlimited horizontal surface as described in figure 4 and Eq.28. Only a few months ago Cabeza-Lainez obtained a similar solution for a triangle with a horizontal side and also for circular sectors if the center of the vertical circle lies on the horizontal surface.
In this quite elaborate fashion, we have discovered how complex it turns to find an analytical expression to solve the problem if only we alter the coordinate plane of calculation not to mention figures of geometry other than the regular triangle or rectangle. That is why we concluded that a graphical solution based on parametric design, would imply a much more general solution and extraordinarily simplify the issue. This is the reason behind our method. Since f12 must be a dimensionless configuration factor, we ought to divide by the area of the projected unit sphere, that is πr 2 = π, from which we would arrive to: We can geometrically obtain that the lengths of the minor semi-axes of the projection of the two ellipses involved which are, respectively, b1 = 0.8941 and b0 = 0.7071. The major semi-axis is logically Strangely enough, the only way to solve the problem is changing the coordinates of the problem to (very rare) generalized polar coordinates, and thus we would have: x = a ρ cos θ; y = b ρ cos θ, thus, the differential element yields dA = a b ρ dρ dθ.
This is, in fact, the only possible procedure to obtain the area of an elliptic sector: For the larger ellipse y/x = 2/1.5 (by proportion); a = 1, b = 0.8944.
The area of the larger ellipse is 0.5984·0.8944 = 0.5284 and accordingly for the small ellipse the area would be 0.3447; therefore, the final area yields 0.5284 − 0.3447 = 0.1837.
It is important to stress that the values obtained for the configuration factors through the application of the principle of the projected solid angle are not approximations to those obtained by analytical calculations, but are identical and exactly the same. Such little demonstration, simple as it is, due to the extreme difficulties arisen with integration of elliptic geometries, is a finding solely attributed to Cabeza-Lainez. In this way he has dissipated the clouds that spread over the validity of such a convenient principle for our method and beyond.

ANNEX 3. The recurrent Problem of Circular emitters
Circular emitters can be treated in the same way as any type of emitter virtue of the law of the projected the solid angle ( Figure A3.1). Nonetheless, explicit solutions were very unassailable. Cabeza-Lainez proposed that through a so-called "pair of dejection" (Figure A3.2), we could determine the parameters of the cone that subtends the emitting circle of radius R, at point P. The resulting equation of the elliptical cone is: The intersection between the oblique cone of parameters α and β and a sphere of unit radius, defines itself by the expression: x 2 + y 2 + z 2 = 1 Although it is actually a fourth degree curve as shown ( Figure A3.5) in the ensuing representation and even it is not possible in most cases to find its enclosed area by numerical procedures.