# Understanding Classical Gray Body Radiation Theory

November 1, 2018

Even though the first man-made light source used thermal radiation, the effect wasn’t fully understood until the discovery of quantum mechanics. Nowadays, it’s a well-known physics concept. In this blog post, we discuss surface-to-surface radiation theory for the so-called gray body, how to implement it in the COMSOL Multiphysics® software, and an interesting use of this theory.

Candlelight was the first man-made light source. Much later, gas lights emerged. Then, in the 19th century, the first electric incandescent lamp was invented. 140 years later, the incandescent lamp is being replaced by the LED.

Candles, gas lights, and incandescent lamps all use thermal radiation, the only heat transfer phenomenon that works without the presence of heat-conductive media. Without this special property, we wouldn’t be able to live on Earth: Solar energy would pour in by this mechanism through a vacuum.

An incandescent light bulb.

One of the most interesting stories in physics is when Max Planck discovered the formula of the spectral distribution, now called the Planck distribution, for black body radiation. The black body is an idealized physical body that absorbs all incident electromagnetic radiation; equivalently, a body that emits maximized thermal radiation.

The discovery of the black body theory opened the door to quantum mechanics. With Planck’s work, we now know how much power we get from a black body object at temperature T in an equilibrium state, which is called the Stefan-Boltzmann law:

I = n^2\sigma T^4 \hspace{1cm} (W/m^2),

where \sigma \simeq 5.67 \times 10^{-8} \; \; \; (W/m^2 \cdot K^4) is the Stefan-Boltzmann constant and n is the refractive index of the media.

We now consider more practical objects, called gray bodies. A gray body is an imperfect black body; i.e., a physical object that partially absorbs incident electromagnetic radiation. The ratio of a gray body’s thermal radiation to a black body’s thermal radiation at the same temperature is called the emissivity of the gray body.

The emissivity of a black body is unity, while that of a gray body is larger than 0 and smaller than 1. The emissivity is a function of the geometry of the radiative surface, its physical properties, and the wavelength. So, how is emissivity determined from certain structures?

For the sake of simplicity, we will consider the diffuse gray surface, which doesn’t account for the spectral dependency of the radiation, throughout this section. The following theory is based on Gouffé’s paper (Ref. 1), which outlines the classical gray body radiation theory. The math is something that you might have learned in middle or high school, but the physical concept is a little more complicated. Gouffé’s paper was cited by many researchers at the time it was published (Ref. 7). In this blog post, we will try to include what is missing from the original paper in order to provide a complete explanation.

First of all, let’s define our terms and conventions. When we mention reflectivity and absorptivity in this blog post, we must differentiate between the “material” and “apparent” quantities. This is very important in order to avoid confusion when learning about thermal radiation. Because gray bodies typically have certain (tiny) structures on their surfaces, the “apparent” quantity viewed from a far field is different from the raw material quantity. For example, the apparent reflectivity of a rough surface is always lower than its material reflectivity.

Our goal is to calculate the apparent emissivity from a given material reflectivity and the geometry of the structure. Here, we will differentiate “material” quantities from “apparent” quantities by appending “_0” in the notation. Otherwise, it is understood that we are talking about a quantity in general. The nomenclature of the terms is as follows:

• Material reflectivity, \rho_0
• Material absorptivity, \alpha_0
• Material emissivity, \varepsilon_0
• Apparent reflectivity, \rho
• Apparent absorptivity, \alpha
• Apparent emissivity, \varepsilon
• Cavity’s internal area, S
• Opening area, A
• View angle, \theta
• View factor (normalized solid angle), F
• View factor (area ratio), G

From the conservation law of energy, we have the following relation for opaque materials between the reflectivity \rho and the absorptivity \alpha:

1= \rho+\alpha

Kirchhoff’s law states that the material emissivity \varepsilon at thermodynamic equilibrium is equal to the absorptivity \alpha; i.e.,

\varepsilon = \alpha.

From the two relations above, we get the emissivity from the reflectivity; i.e.,

\varepsilon = 1-\rho.

We consider a single surface structure, as shown in the following figure. The structure can be anything, but in this blog post, we use a spherical cavity with an opening on top with area A (circle of radius R) at distance L from the bottom of the cavity. We assume that the material reflectivity \rho_0 is uniform over the internal surface of the cavity and the reflection takes place according to Lambert’s law; i.e., the intensity of the reflection is I_r = \rho_0 \cos \theta, where \theta is the viewing angle, as depicted in the figure. We want to calculate how much apparent reflection we get out of an incident light of the energy of unity.

Surface structure for calculating the apparent reflectivity.

As the first-order approximation, the reflection from the bottom of the cavity is:

\rho_0F,

where F is called the view factor.

Note that in Gouffé’s paper, the view factor is assumed to be uniform over the cavity surface. In the COMSOL® software, it is instead a function of the position, which is always true. (So, we need integration to compute the emissivity.)

In this case, the view factor is the normalized solid angle to the opening from the first reflection point. The solid angle \Omega is calculated as:

\Omega = \int_0^\theta \int_0^{2\pi} \sin \theta ^\prime \cos \theta ^\prime d\theta ^\prime d\varphi ^\prime=\pi \sin^2 \theta

Note that the total solid angle for the hemisphere is \pi, not 2\pi. This is due to the Lambertian factor, \cos \theta.

As a result, the view factor F is

F=\sin^2 \theta.

To this approximation, the apparent reflectivity is

\rho^{(1)}=\rho_0 F,

from which the apparent emissivity is derived as

\varepsilon^{(1)} = 1-\rho_0 F.

Roughly speaking, the smaller the opening, the more the cavity becomes a black body due to the view factor.

Next, let’s improve the approximation. After the first reflection, which we already calculated, the rest is absorbed by the cavity material or contributes to further reflections. The material absorption \alpha_0 is \alpha_0=1-\rho_0 from the conservation of energy. The energy left for the subsequent reflections is:

1-\rho_0F-\alpha_0= \rho_0(1-F).

Now, assuming again that further reflections take place in a uniform way, the reflection that gets out of the cavity at the second reflection should be the above quantity multiplied by the material reflectivity one more time and another view factor defined by the area ratio, G= A/S, where S is the area of the cavity, including the opening area A. The apparent reflectivity for the second reflection is

\rho_0^2(1-F)G

Similarly to the first-order approximation, the apparent emissivity for the second-order approximation is:

\varepsilon^{(2)}=1-\rho_0F-\rho_0^2(1-F)G

To this approximation, the apparent emissivity should become more accurate with the additional term.

Finally, we can take all of the reflections into account by calculating the following converging infinite series:

\varepsilon^{(\infty)}=1-\rho_0F-\rho_0^2(1-F)G-\rho_0^3(1-F)G(1-G)-\rho_0^4(1-F)G(1-G)^2-\cdots =\frac{(1-\rho_0)( 1+\rho_0(G-F) ) }{1-\rho_0(1-G)}.

In the case of a sphere, it can be shown that F=G, which reduces this result to

\varepsilon^{(\infty)}=\frac{1-\rho_0 }{1-\rho_0(1-G)}.

Now, let’s rewrite F and G explicitly in terms of the geometry parameters R and L. From the geometry, it’s easy to prove that \sin \theta = R/\sqrt{L^2+R^2}, from which F can be rewritten as

F=\sin^2 \theta = \frac{R^2}{R^2+L^2}= \frac{1}{1+\left (\frac LR \right )^2}.

Next, the opening area A can be calculated by

A=\int_\varphi^{\pi/2} 2\pi r^2 \cos \varphi ^\prime \ d\varphi ^\prime = 2\pi r^2 (1-\sin \varphi),

where r is the radius of the sphere and \varphi is the angle between the plane intersecting the sphere center, which is parallel to the opening perimeter.

From the geometry, there are relations that R^2+(L-r)^2= r^2, which yields that r=(R^2+L^2)/(2L) and \cos \varphi = R/r. By using these relations, the opening area A can be rewritten as

A=2\pi r(2r-L).

Since the total sphere area is S=4\pi r^2, the view factor G is given as

G=\frac A S = \frac{2\pi r(2r-L)}{4\pi r^2} = \frac{1}{1+\left( \frac L R \right )^2},

which proves that F=G.

Note that in the COMSOL® software, we don’t have the view factor F, only the view factor G, which is the purely geometrical quantity.

### COMSOL Multiphysics® Simulation Versus Theory

So far, we’ve learned the classical gray body theories. Now, let’s compute our final goal, the apparent emissivity for a gray body consisting of a spherical shell with an opening, by using the Heat Transfer with Surface-to-Surface Radiation interface. Then, we can compare the computed dependency of the opening size on the emissivity with the approximate formulas.

Before moving to the main results, let’s check one thing. In COMSOL Multiphysics, we can calculate the view factor values (corresponding to G in the above theory, but as a function of position) using the radopu and  radopd operators. The quantity that is obtained by using these operands is purely geometric, so there is no need to run a study to use this operator. If you change the geometry, you can just update the solution before using the operator.

In this model, there is only one radiating surface because the model is an open cavity. So, COMSOL Multiphysics can calculate only the view factor for the cavity surface itself. We’ll call this the self view factor. The view factor G, which we have discussed so far as the factor that views outside of the cavity from a point on the cavity surface, can be called the ambient view factor. We can calculate the view factor G by subtracting the self view factor from unity; i.e., 1-intop1(comp1.ht.radopu(1,0))/intop1(1).

As mentioned earlier, the view factor is position dependent in general, and so is the view factor in the COMSOL® software. To calculate the contribution from all points on the cavity surface, we need to integrate the view factor as a function of position. The operator intop1 is an integration operator defined on the cavity surface. The following figure compares the calculated results with the classical view factor theory, showing very good agreement.

Comparison between the classical theory and the radopu operator calculation in COMSOL Multiphysics for the ambient view factor.

Now, let’s move on to the settings. To calculate surface-to-surface radiation, a Surface-to-Surface Radiation physics interface and a Diffuse Surface node for the internal surface of the sphere need to be added, and a Heat Transfer with Surface-to-Surface Radiation multiphysics coupling node needs to be added (see below). COMSOL Multiphysics will compute the view factor by performing surface integration for each mesh node except when the geometry is axisymmetric where a closed-form expression is used. For the sake of simplicity, we don’t include ambient radiation; i.e., T_{\rm amb}=0. A Temperature boundary condition is set to the outer boundary of the spherical shell with T_0 = 2500 deg K. The material emissivity is set to 0.5 in order to see the difference more easily.

Settings for the Diffuse Surface node.

The following plots show the computed ambient view factor, radiosity, and temperature for various opening radii. Following Gouffé’s paper, the geometry aspect ratio defined by L/R, where R is the opening’s radius and L is the cavity height, is used as a sweep parameter.

Computational results for the ambient view factor, radiosity, and temperature versus L/R.

These results qualitatively suggest that the smaller the size of the opening, the more the cavity emits thermal radiation and the higher the surface temperature.

To calculate the apparent emissivity for the computed radiosity result, we can now use the Stefan-Boltzmann law. The radiosity of the black body at temperature T_0 is \sigma T_0^4. Dividing the computed radiosity by this number, we get the apparent emissivity of the gray body spherical cavity, which is intop1(ht.J)/intop1(1)/(sigma_const*T0^4). Let’s now compare the results with the various theories we have discussed so far.

Comparison between COMSOL Multiphysics results and various theoretical approximations for the apparent emissivity for a spherical cavity.

The cyan line shows the first-order approximation, which underestimates the amount of reflection causing the higher emissivity, which is inaccurate. The second order obviously improves the accuracy, but it’s not enough. The plot for the infinite series approximation (orange line) gives a much better result, but it’s still not very close to the COMSOL Multiphysics results.

The reason for this discrepancy is mentioned in the last part of Gouffé’s paper. As some readers might have noticed in the previous figure, the surface temperature is different for each geometry, regardless of the same outer boundary temperature of 2500 deg K. This is called radiation cooling. We have yet to include this effect, which takes place for any real material with a finite thermal conductivity. The effect changes the temperature of the inner surface of the shell temperature, depending on the opening area. Therefore, the temperature needs to be compensated in order to calculate the emissivity at the same temperature for all opening areas.

The correction factor is, owing to the Stefan-Boltzmann law, the fourth power of the temperature ratio; i.e., (maxop1(T))^4/T0^4, where  maxop1 is an operator defined on the cavity surface that finds the maximum value on the surface. Finally, with this correction, the red curve is the most accurate theoretical prediction, which agrees very well with the COMSOL Multiphysics results (blue line).

### Designing an Efficient Incandescent Lamp with the COMSOL® Software

The light source in an incandescent lamp is created by a twisted tungsten filament. The material emissivity of tungsten is 0.462 at 2500 deg K for a 0.467-um wavelength (Ref. 2). In the past, researchers proposed that we can design more efficient incandescent lamps if we fabricate some microstructures on the surface of the tungsten filament. This is true. As we just learned, the apparent emissivity can be close to 1 if we have a cavity with a very tiny opening, which we can call a black body filament. In addition, research (Ref. 3–4, 8–9) also proposed that if the maximum cavity size is about a half of 0.78 um, then any infrared light with a wavelength larger than 0.78 um may be suppressed due to the waveguide’s cutoff effect. Then, the efficiency for the visible light may be significantly improved.

This effect may be more than the emissivity enhancement, because thermal radiation below tungsten’s melting temperature consists of mostly infrared light that is wasted to heat. It would be fantastic if we could really cut the infrared off.

Our “dream” incandescent lamp with an infrared suppressing black body filament. Background image in the public domain, via Wikimedia Commons.

Unfortunately, this “dream lamp” can’t be a reality for various reasons.

First, it isn’t possible for the surface to be made up of “all” holes. To make holes with a low view factor, meaning a smaller opening area than the hole size, the surface can’t be densely filled with the openings; i.e., the surface needs a flatter area between the openings, which radiate infrared. Alternatively, we can make deep holes to decrease the ambient view factor and then make each hole as close to the adjacent holes as possible. However, this makes it more difficult to fabricate deeper holes.

Also, the surface energy of the surface with holes seems to be higher than that of a flat surface. So, the surface with holes tends to melt down at a lower temperature than the bulk melting temperature.

Third, it is very difficult to make holes on a twisted wire filament because it’s a 3D structure. It is easier to make holes on a flat ribbon filament because it’s 2D, but flat filaments are not electrically convenient because they need more current to achieve the same temperature than the twisted wire filament (higher voltage and lower current is convenient for our current infrastructure).

### Concluding Remarks

Although gray body radiation theory was developed a long time ago, the theory is still well organized. There are many sources for closed-form view factor formulas (Ref. 5; note that it’s also called “configuration factor”). This theory has been proven by experiments and is still applied to many real-world applications, but it can seem too difficult to understand. In this blog post, we learned the theory for the simplest example and saw a very successful benchmark result with COMSOL Multiphysics. The Heat Transfer with Surface-to-Surface Radiation interface is a reliable tool that makes complicated radiation computations easier.

We can come up with a more interesting incandescent lamp by making a special view factor. If we make tiny periodic grooves on a flat ribbon filament, we will get a “directional” incandescent lamp. Usually, a filament emits light to all directions, but this special filament illuminates only a certain direction (Ref. 6). We can explore many more ideas and many more applications of gray body radiation theory with simulation.

### Next Step

Try the Gray Body model yourself by clicking the following button, which will take you to the Application Gallery. There, you can download the related MPH-file.

### References

1. A. Gouffé, “Corrections d’ouverture des corps-noirs artificiels compte tenu des diffusions multiples internes (Corrections of emissivity for the artificial black-body considering multiple internal diffusions)” (in French), Revue d’Optique, t. 24, no. 1–3, 1945.
2. C.J. Smithells, Tungsten, Chapman & Hall, Ltd., 1926.
3. J.F. Waymouth, Proceedings of LS-5, 1989.
4. J.F. Waymouth, J. Illum. Engng. Inst. Jpn., 74, 12, pp. 800–805, 1990.
5. J.R. Howell, “A Catalog of Radiation Heat Transfer Configuration Factors“, University of Texas at Austin.
6. W.Z. Black, “Radiative heat transfer characteristics of specially prepared V-groove cavities”, PhD Thesis to Purdue University, 1968.
7. Precision Measurement and Calibration, Radiometry and Photometry (Vol. 7), U.S. Department of Commerce, National Bureau of Standards.
8. M. Sugimoto, T. Fujioka, T. Inoue, H. Fukushima, Y. Mizuyama, S. Ukegawa, T. Matsushima, and M. Toho, “The Infrared Suppression in the Incandescent Light from a Surface with Submicron Holes”, Journal of Light & Visual Environment, Vol. 18, No. 2 (1994).
9. T. Kondo, S. Hasegawa, T. Yanagishita, N. Kimura, T. Toyonaga, and H. Masuda, “Control of thermal radiation in metal hole array structures formed by anisotropic anodic etching of Al”, Optics Express Vol. 26 No. 21 (2018).

#### Categories

##### Marcel vanWijk
February 6, 2019

Would this mean that’t is possible to have much more efficient incandescent lights in the future?

##### Yosuke Mizuyama
February 6, 2019

Dear Marcel,
Thank you for reading my blog.
It’s obviously possible to make such efficient incandescent lights but not with the same quality, that is, the life can be very short for example.

##### Robert Bedard
January 1, 2020

Wow, that is a great read! How would you suggest to proceed regarding the amount of energy radiating from a flame vs convection in an heat exchanger (water) ? The combustible being methane and the flame being supplied with the appropriate amount of air (oxygen) for an ideal burner. I guess Radiation in the UV and visible light wavelength regions would not satisfy a black nor a grey body definition, though would the post flame region qualify as a grey body, if the reaction is in thermodynamic equilibrium? How could I simulate the heat exchange in the simplest way possible with Comsol MP, taking in account radiation, convection and thin layer conduction ?

##### Yosuke Mizuyama
January 2, 2020 COMSOL Employee

Dear Robert,
Thank you for reading my blog. I’m not familiar with radiation from a flame and I assume that some or all parts of the spectrum of a methane gas flame wouldn’t follow the Stefan-Boltzmann law. But if you think you can approximate it as a gray body, you can define it by the user defined emissivity function. As I wrote in the blog post, the spectral shape of a radiation can be arbitrarily changed by changing the emissivity (as far as it satisfies the Stefan-Boltzmann law), you can model any radiation. Then the tool you need to calculate the radiation propagating back and forth in your heat exchanger is COMSOL’s Surface-to-Surface Radiation interface.

##### Robert Bedard
January 8, 2020

Mr Mizuyama,
Thank you for the tips, they are really helpful. I think your assumption is correct, some parts might not, though a 2011 study from L.Rossy seems to indicate there is a way for outdoor flame, so I figure there should be one too for heat exchanger. I am going to work on this.

##### Yosuke Mizuyama
January 8, 2020 COMSOL Employee

Dear Robert,
Great! I hope your research goes well.
Best regards,
Yosuke

##### Tianle Yuan
December 11, 2021

Mr. Mizuyama,

Could you explain where the solid angle expression comes from? Honestly, for me when I try to calculate the solid angle for a sphere it should be sin(theta)d(theta)d(fi). Where is your cos(theta) come from? By the way, seems like the sequence of integration should swap.

##### Yosuke Mizuyama
December 12, 2021 COMSOL Employee

Dear Tianle,
Thank you for reading my blog. The factor cos(theta) is due to the Lambert cosine law. I had written about it in the note in the blog. You can’t literally view a surface at 90 degree. Thank you for your note the integration notation out. I agree.
Best regards,
Yosuke

##### Tianle Yuan
December 12, 2021

Mr. Mizuyama,

I have one more question here for the final diagram. Do you think the apparent value that is independent of spatial position is reasonable? Whether it is possible that Comsol can show apparent reflection for a different position at the hole? Right now I think the “theta” used for the integration is a constant which will cause the result with only a uniform value. So if we watch the cavity from the outside, the hole will only have a uniform dark value.

Best,
Tianle

##### Yosuke Mizuyama
December 17, 2021 COMSOL Employee

Hi Tianle,
Sorry for the delay. I think I understand your question. I mentioned about it at the last part of the blog. You can make your view factor angle dependent. So you can have angle– or position dependent apparent reflectivity.
Best regards,
Yosuke

##### Tianle Yuan
December 21, 2021

Mr. Mizuyama,

I got what you mean from the last reply. For double-checking, do you mean the F configuration factor should be integrated according to the position of the point inside the cavity?

Best,
Tianle

##### Yosuke Mizuyama
December 24, 2021 COMSOL Employee

Dear Tianle,
I think you made a good point. In this Gouffe’s theory, two things are assumed. 1. Only the bottom point is considered to calculate the view factor F, 2. All the internal reflections are uniform so there is no position dependence inside the cavity. So, you don’t need to integrate F. Only one point represents F. This is Gouffe’s model. This may not be a good approximation for larger opening cases or some applications with flat structures or cases where there is certain position dependences for the reflection/absorption. I encourage you to explore such cases using COMSOL.
Thank you for the good question!
Best regards,
Yosuke

##### Tianle Yuan
December 26, 2021

Mr. Mizuyama,

Best,
Tianle Yuan

##### Rana
October 13, 2023

Is it possible to model grey body radiative transfer if the materials have a wavelength and temperature dependent absorption coefficient?

##### Yosuke Mizuyama
October 13, 2023 COMSOL Employee

Dear Rana,
Thank you for reading my blog. Typically the wavelength and temperature dependence of the absorption coefficient comes from the geometry of the material surface as described in this blog. So, you can model it as you read in this blog. But I assume that you are talking about the “bulk” material properties. If so, yes, you can define a wavelength and temperature dependent material properties under Definition by using an interpolation function. So you would include both (geometrical) surface and bulk material properties in your model.
Best regards,
Yosuke

EXPLORE COMSOL BLOG