PDF

Radiative Heat Transfer in Finite Cylindrical Media
This application uses the discrete ordinates method (DOM) to solve a 3D radiative transfer problem in an emitting, absorbing, and linearly anisotropic-scattering finite cylindrical medium. Using the S6 quadrature of DOM leads to accurate results, which are needed in combined modes of heat transfer. The calculated incident radiation and heat fluxes agree well with published results obtained by transformed integral methods (see Ref. 1). In addition, The DOM formulation of the Heat Transfer Module easily handles the effects of boundary emission and reflection.
Introduction
There are numerous engineering applications of radiative transfer in absorbing, emitting, and anisotropically scattering media with variable radiation properties. Examples include, among others, coal-fired combustion systems, light-weight fibrous insulations, and heat transfer systems containing small scattering particles. Furthermore, the efficiency of radiative transfer depends on the boundary conditions, for example, the temperature and the emissivity of the surrounding walls, and the target where heat transfer is desired. Studies have shown that radiative transfer is highly sensitive to the wall emissivity. In this study you build a validation model representing a cylinder with homogeneous walls. Then you go on to consider walls with space-dependent emissivity and investigate the effects of albedo and scattering.
Model Definition
In this tutorial you validate the DOM formulation in COMSOL Multiphysics by examining three benchmark cases. You investigate the method’s efficiency through parametric analyses by changing the single-scattering albedo, wall emissivity, and linear function. In particular, in the final case, you approximate the scattering phase function by a linear function reflecting highly backward, isotropic, and highly forward scattering.
The model geometry, shown in Figure 1, is a cylinder of radius R = 0.5 m and height L = 1 m.
Figure 1: Schematic diagram of the physical model.
These examples use the S6 discrete ordinate method for predicting the heat flux on the enclosure side walls and the incident radiation distribution inside the domain.
Verification case
For the initial study, assume cold boundaries, that is, Ii = 0 where i = 1, 2, 3 refers to the surface index in Figure 1. Furthermore, assume that the walls diffusively reflect radiation, that is, εi = 0.5 for i = 1, 2, 3. The medium is at a uniform temperature T such that the blackbody radiation intensity in an arbitrary direction per unit area and solid angle  Ib(T) = σT4 ⁄ π equals 1 W/(m2·sr).
Walls With variable Emissivity
Two cases are computed for comparison purposes. These 3D cases represent opaque partial side wall diffuse emission and/or reflection. In both cases, the radiosity on the side walls (Surface 3) varies with the angular coordinate along the full height of the finite cylinder according to ε3 = (1 − y ⁄ R) ⁄ 2. Both cases also have cold walls at the cylinder top (Surface 1) and bottom (Surface 2).
Case 1 has isotropic scattering function and compares results for different albedos. Case 2 has constant albedo and compares results with highly forward, isotropic, and highly backward scattering function parameterized by the Legendre coefficient a1.
albedo = 0.1, 0.5, 0.9
a1 = 0
albedo = 0.5
a1 = −0.99, 0, 0.99
Thermal Analysis
The discrete ordinates method relies on the discrete representation of the directional dependence of the radiation intensity. It involves solving the radiative equation for a set of directions that span the full solid angle range of 4π around a point in space.
The radiation transport equation (RTE) for this type of configuration can be written as
where
I(Ω, s) is the radiation intensity at a given position s in the direction Ω
T is the temperature
κ, β, and σs are absorption, extinction, and scattering coefficients, respectively
Ib(T) is the blackbody radiation intensity
φ(Ω, Ω′) = 1 + a1μ0 where μ0 = Ω ⋅ Ω′ is the cosine of the scattering angle.
The boundary intensities at the cylinder walls are given by the sum of the effective emitted intensity and the reflected incident intensities in the given direction:
where
εw is the surface emissivity, which is in the range [0, 1]
ρd = 1 − εw is the diffusive reflectivity
n is the outward normal vector
qout is the heat flux striking the wall:
The above equations can be discretized in Cartesian coordinates for monochromatic or gray radiation as
The Sn approximation of the RTE in the direction i can be expressed as
For a discrete direction, Ωi, the values of Ωi,1, Ωi,2, and Ωi,3 define the direction cosines of Ωi obeying the condition Ωi,12 + Ωi,22 + Ωi,32 = 1. The index j in the above equation denotes the direction of incoming radiation contributing to the direction Ωi.
For a diffuse reflecting surface on a wall boundary, the boundary condition equation is transformed as
Results and Discussion
The results demonstrate that the DOM procedure for the prediction of radiation is an elegant and accurate method for modeling multidimensional radiative heat transfer in cylindrical geometries.
Verification case
This case treats the effects of the scattering albedo on the incident radiation and heat fluxes. Figure 2 shows the distribution of the net heat flux qrnet(R0z) versus axial optical thickness. There is good overall agreement of the present work with published literature results (Ref. 2, Ref. 3, and Ref. 4).
Figure 2: The effects of the scattering albedo on the radial heat flux qrnet(R0z); for a hot cylindrical medium enclosed by cold walls, ε1 = ε2 = ε3 =0.5.
The effects of albedo on the distribution of centerline incident radiation in axial direction are shown in Figure 3. The incident radiation is symmetric with respect to z = L ⁄ 2 plane and decrease with increasing scattering albedo. Furthermore, results become more uniform with larger scattering albedo.
Figure 3: The effects of the scattering albedo on the distribution of centerline incident radiation G(00z) for a hot cylindrical medium enclosed by cold walls, ε1 = ε2 = ε3 = 0.5.
Figure 4 displays the distributions of the incident radiation across the midplane radius G(x0L ⁄ 2) with respect to normalized optical thickness x ⁄ R.
Figure 4: The effects of the scattering albedo on the distributions of the incident radiation G(x, 0, L ⁄ 2) with respect to normalized optical thickness x ⁄ R for a hot cylindrical medium enclosed by cold walls, ε1 = ε2 = ε3 = 0.5.
Walls With variable Emissivity
The incident radiation at the midplane z = L ⁄ 2 at the radial position R ⁄ 2 is shown in Figure 5. At the side surface R ⁄ 2 distance from the cylinder axis, G(R0L ⁄ 2) changes with the azimuthal angle. The changes in scattering albedo are also illustrated. The smallest albedo makes the biggest change around the azimuthal angle.
Figure 5: The effects of scattering albedo on the distribution of incident radiation at midplane z = L ⁄ 2 at R ⁄ 2 distance from the cylinder axis with respect to an azimuthal coordinate for Case 1.
Figure 6 shows the effect of a nonzero linear anisotropic scattering coefficient a1. As expected, the differences between isotropic scattering, forward scattering, and backward scattering are most accentuated far from the boundary.
Figure 6: The effects of linear anisotropic scattering coefficient a1 on the distribution of incident radiation G(0yL ⁄ 2) with respect to normalized optical thickness y ⁄ R for Case 2.
References
1. X.L. Chen and W.H. Sutton, “Radiative Transfer in Finite Cylindrical Media Using Transformed Integral Equations,” J. Quantitative Spectroscopy and Radiative Transfer, vol. 77, pp. 233–271, 2003.
2. X. Chen, Transformed Integral Equations of Radiative Transfer and Combined Convection-radiation Heat Transfer Enhancement with Porous Insert, Doctoral Thesis, University of Oklahoma, 2003.
3. S.T. Thynell and M.N. Ozisik, “Radiation Transfer in Absorbing, Emitting, Isotropically Scattering, Homogeneous Cylindrical Media,” J. Quantitative Spectroscopy and Radiative Transfer, vol. 38, no. 6, pp. 413–426, 1987.
4. H.Y. Li, M.N. Ozisik, and J.R. Tsai, “Two-dimensional radiation in a cylinder with spatially varying albedo,” J. Thermophysics and Heat Transfer, vol. 6, no. 1, pp. 180–182, 1992.
Application Library path: Heat_Transfer_Module/Verification_Examples/cylinder_participating_media
Modeling Instructions — Validation Case
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Heat Transfer>Radiation>Radiation in Participating Media (rpm).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Stationary.
6
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Cylinder 1 (cyl1)
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, locate the Size and Shape section.
3
In the Radius text field, type R.
4
In the Height text field, type L.
5
Click  Build All Objects.
Materials
Add a material to specify the absorption and scattering coefficients inside the cylinder.
Domain
1
In the Materials toolbar, click  Blank Material.
2
In the Settings window for Material, type Domain in the Label text field.
3
Locate the Material Contents section. In the table, enter the following settings:
Analogously, specify the emissivity of the walls using a material.
Walls
1
In the Materials toolbar, click  Blank Material.
2
In the Settings window for Material, type Walls in the Label text field.
3
Locate the Geometric Entity Selection section. From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose All boundaries.
5
Locate the Material Contents section. In the table, enter the following settings:
Radiation in Participating Media (rpm)
1
In the Model Builder window, under Component 1 (comp1) click Radiation in Participating Media (rpm).
2
In the Settings window for Radiation in Participating Media, locate the Participating Media Settings section.
3
Find the Radiation settings subsection. From the Pindex list, choose 0.3.
With this performance index value, solving the model requires roughly 2 GB of RAM. If your computer has less available memory than that, try a value between 0.5 and 1.
4
From the Discrete ordinates method list, choose S6 (48 directions).
Participating Medium 1
1
In the Model Builder window, expand the Radiation in Participating Media (rpm) node, then click Participating Medium 1.
2
In the Settings window for Participating Medium, locate the Model Input section.
3
In the T text field, type T0.
Opaque Surface 1
1
In the Model Builder window, click Opaque Surface 1.
2
In the Settings window for Opaque Surface, locate the Model Input section.
3
In the T text field, type Tw.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Coarse.
Free Triangular 1
1
In the Mesh toolbar, click  Boundary and choose Free Triangular.
2
Swept 1
1
In the Mesh toolbar, click  Swept.
2
In the Settings window for Swept, click  Build All.
Study 1
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
If your computer has 6GB of RAM or more, it is possible to greatly decrease the solution time by using a fully coupled solver.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
The tolerance is decreased to reach better accuracy.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Stationary Solver 1.
3
In the Settings window for Stationary Solver, locate the General section.
4
In the Relative tolerance text field, type 1e-4.
5
Right-click Study 1>Solver Configurations>Solution 1 (sol1)>Stationary Solver 1 and choose Fully Coupled.
6
Right-click Study 1>Solver Configurations>Solution 1 (sol1)>Stationary Solver 1>AMG, radiation variables and choose Enable.
7
In the Study toolbar, click  Compute.
Results
Incident Radiation (rpm)
The first default plot shows the incident radiation distribution in 3D and the second default plot represents the net radiative heat flux distribution in 3D, see figure below.
Add a new 1D Plot to represent the net radiative heat flux along the z-coordinate and compare with Figure 2.
Net Radiative Heat Flux vs. z, 1D
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Net Radiative Heat Flux vs. z, 1D in the Label text field.
Line Graph 1
1
In the Net Radiative Heat Flux vs. z, 1D toolbar, click  Line Graph.
2
Click the  Transparency button in the Graphics toolbar.
3
4
Click the  Transparency button in the Graphics toolbar again to return to the original state.
5
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1)>Radiation in Participating Media>Boundary fluxes>rpm.qr_net - Net radiative heat flux - W/m².
6
Click Replace Expression in the upper-right corner of the x-Axis Data section. From the menu, choose Component 1 (comp1)>Geometry>Coordinate>z - z-coordinate.
7
Click to expand the Legends section. Select the Show legends check box.
Finish the plot by adjusting the title and axis labels.
Net Radiative Heat Flux vs. z, 1D
1
In the Model Builder window, click Net Radiative Heat Flux vs. z, 1D.
2
In the Settings window for 1D Plot Group, click to expand the Title section.
3
From the Title type list, choose None.
4
Locate the Plot Settings section. Select the x-axis label check box.
5
Select the y-axis label check box.
6
In the Net Radiative Heat Flux vs. z, 1D toolbar, click  Plot.
Cut Line 3D 1
1
In the Results toolbar, click  Cut Line 3D.
2
In the Settings window for Cut Line 3D, locate the Line Data section.
3
In row Point 2, set X to 0, and z to L.
4
The Graphics window shows the location of the line in the model geometry.
Add a new 1D Plot to represent the incident radiation along the z-coordinate and compare with Figure 3.
Incident Radiation vs. z, 1D
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Incident Radiation vs. z, 1D in the Label text field.
3
Locate the Data section. From the Dataset list, choose Cut Line 3D 1.
Line Graph 1
1
In the Incident Radiation vs. z, 1D toolbar, click  Line Graph.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the x-Axis Data section. From the menu, choose Component 1 (comp1)>Geometry>Coordinate>z - z-coordinate.
3
Locate the Legends section. Select the Show legends check box.
Incident Radiation vs. z, 1D
1
In the Model Builder window, click Incident Radiation vs. z, 1D.
2
In the Settings window for 1D Plot Group, locate the Title section.
3
From the Title type list, choose None.
4
Locate the Plot Settings section. Select the x-axis label check box.
5
Select the y-axis label check box.
6
In the Incident Radiation vs. z, 1D toolbar, click  Plot.
Cut Line 3D 2
1
In the Results toolbar, click  Cut Line 3D.
2
In the Settings window for Cut Line 3D, locate the Line Data section.
3
In row Point 1, set Z to L/2.
4
In row Point 2, set X to R, and z to L/2.
5
Add a new 1D Plot to represent the incident radiation along the x-coordinate and compare with Figure 4.
Incident Radiation vs. x, 1D
1
Right-click Incident Radiation vs. z, 1D and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Incident Radiation vs. x, 1D in the Label text field.
3
Locate the Data section. From the Dataset list, choose Cut Line 3D 2.
4
Locate the Plot Settings section. In the x-axis label text field, type x-coordinate (m).
Line Graph 1
1
In the Model Builder window, expand the Incident Radiation vs. x, 1D node, then click Line Graph 1.
2
In the Settings window for Line Graph, click Replace Expression in the upper-right corner of the x-Axis Data section. From the menu, choose Component 1 (comp1)>Geometry>Coordinate>x - x-coordinate.
3
In the Incident Radiation vs. x, 1D toolbar, click  Plot.
Modeling Instructions — Case 1
Radiation in Participating Media (rpm)
Opaque Surface 2
1
In the Physics toolbar, click  Boundaries and choose Opaque Surface.
2
In the Settings window for Opaque Surface, locate the Model Input section.
3
In the T text field, type Tw.
4
These are the vertical wall boundaries. To reach all of them, you can rotate the geometry or click either the Transparency button or the Wireframe Rendering button in the Graphics toolbar.
5
Locate the Surface Radiative Properties section. From the ε list, choose User defined. In the associated text field, type ew*(1-y/R).
Now, disable Opaque Surface 2 in Study 1 to be able to re run the same Study 1 configuration.
Study 1
Step 1: Stationary
1
In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step check box.
4
In the tree, select Component 1 (Comp1)>Radiation in Participating Media (Rpm)>Opaque Surface 2.
5
Click  Disable.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Stationary.
4
Click Add Study in the window toolbar.
Study 2
Step 1: Stationary
In the Home toolbar, click  Add Study to close the Add Study window.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
Solution 2 (sol2)
In the Study toolbar, click  Show Default Solver.
Solver Configurations
In the Model Builder window, expand the Study 2>Solver Configurations node.
Solution 2 (sol2)
1
In the Model Builder window, expand the Study 2>Solver Configurations>Solution 2 (sol2) node, then click Dependent Variables 1.
2
In the Settings window for Dependent Variables, locate the General section.
3
From the Defined by study step list, choose User defined.
4
Locate the Initial Values of Variables Solved For section. From the Method list, choose Solution.
5
From the Solution list, choose Solution 1 (sol1).
6
In the Model Builder window, expand the Solution 2 (sol2) node.
7
Right-click Stationary Solver 1 and choose Fully Coupled.
8
Right-click AMG, radiation variables and choose Enable.
9
In the Study toolbar, click  Compute.
Results
Incident Radiation (rpm) 1
The first default plot shows the incident radiation distribution in 3D and the second default plot represents the net radiative heat flux distribution in 3D, see figure below.
Parameterized Curve 3D 1
1
In the Results toolbar, click  More Datasets and choose Parameterized Curve 3D.
2
In the Settings window for Parameterized Curve 3D, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
4
Locate the Parameter section. In the Name text field, type phi.
5
In the Maximum text field, type 2*pi.
6
Locate the Expressions section. In the x text field, type R*cos(phi)/2.
7
In the y text field, type R*sin(phi)/2.
8
In the z text field, type L/2.
9
Add a new 1D Plot to represent the incident radiation in function of the azimuthal angle and compare with Figure 5.
Incident Radiation vs. Azimuthal Angle, 1D
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Incident Radiation vs. Azimuthal Angle, 1D in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameterized Curve 3D 1.
Line Graph 1
1
In the Incident Radiation vs. Azimuthal Angle, 1D toolbar, click  Line Graph.
2
In the Settings window for Line Graph, locate the Legends section.
3
Select the Show legends check box.
Incident Radiation vs. Azimuthal Angle, 1D
1
In the Model Builder window, click Incident Radiation vs. Azimuthal Angle, 1D.
2
In the Settings window for 1D Plot Group, locate the Title section.
3
From the Title type list, choose None.
4
Locate the Plot Settings section. Select the x-axis label check box.
5
6
Select the y-axis label check box.
Line Graph 1
1
In the Model Builder window, click Line Graph 1.
2
In the Settings window for Line Graph, locate the x-Axis Data section.
3
From the Parameter list, choose Expression.
4
In the Expression text field, type phi.
5
In the Incident Radiation vs. Azimuthal Angle, 1D toolbar, click  Plot.
Modeling Instructions — Case 2
Global Definitions
Parameters 1
For case 2, you need to modify the value of the single-scattering albedo that you defined previously and add a parameter for the Legendre coefficient a1 in the scattering phase function.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Radiation in Participating Media (rpm)
Participating Medium 1
1
In the Model Builder window, under Component 1 (comp1)>Radiation in Participating Media (rpm) click Participating Medium 1.
2
In the Settings window for Participating Medium, locate the Scattering section.
3
From the Scattering type list, choose Linear anisotropic.
4
In the a1 text field, type a1.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Stationary.
4
Click Add Study in the window toolbar.
Study 3
Step 1: Stationary
In the Home toolbar, click  Add Study to close the Add Study window.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
Note that the value a1 = 0 gives the same solution as for ω = 0.5 in case 1.
Solution 3 (sol3)
In the Study toolbar, click  Show Default Solver.
Solver Configurations
In the Model Builder window, expand the Study 3>Solver Configurations node.
Solution 3 (sol3)
1
In the Model Builder window, expand the Study 3>Solver Configurations>Solution 3 (sol3) node, then click Dependent Variables 1.
2
In the Settings window for Dependent Variables, locate the General section.
3
From the Defined by study step list, choose User defined.
4
Locate the Initial Values of Variables Solved For section. From the Method list, choose Solution.
5
From the Solution list, choose Solution 2 (sol2).
6
In the Model Builder window, expand the Study 3>Solver Configurations>Solution 3 (sol3)>Stationary Solver 1 node.
7
In the Model Builder window, expand the Solution 3 (sol3) node.
8
Right-click Stationary Solver 1 and choose Fully Coupled.
9
Right-click AMG, radiation variables and choose Enable.
10
In the Study toolbar, click  Compute.
Results
Incident Radiation (rpm) 2
The first default plot shows the incident radiation distribution in 3D and the second default plot represents the net radiative heat flux distribution in 3D, see figure below.
Parameterized Curve 3D 2
1
In the Results toolbar, click  More Datasets and choose Parameterized Curve 3D.
2
In the Settings window for Parameterized Curve 3D, locate the Data section.
3
From the Dataset list, choose Study 3/Solution 3 (sol3).
4
Locate the Expressions section. In the y text field, type -s*R.
5
In the z text field, type L/2.
With the above definition, s = y/R equals the optical thickness along the negative y-axis for the given constant values of x and z.
6
Add a new 1D Plot to represent the incident radiation in function of the normalized optical thickness and compare with Figure 6.
Incident Radiation vs. Normalized Optical Thickness, 1D
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Incident Radiation vs. Normalized Optical Thickness, 1D in the Label text field.
3
Locate the Data section. From the Dataset list, choose Parameterized Curve 3D 2.
Line Graph 1
1
In the Incident Radiation vs. Normalized Optical Thickness, 1D toolbar, click  Line Graph.
2
In the Settings window for Line Graph, locate the x-Axis Data section.
3
From the Parameter list, choose Expression.
4
In the Expression text field, type s.
5
Locate the Legends section. Select the Show legends check box.
Incident Radiation vs. Normalized Optical Thickness, 1D
1
In the Model Builder window, click Incident Radiation vs. Normalized Optical Thickness, 1D.
2
In the Settings window for 1D Plot Group, locate the Title section.
3
From the Title type list, choose None.
4
Locate the Plot Settings section. Select the x-axis label check box.
5
In the associated text field, type Normalized optical thickness (y/R).
6
Select the y-axis label check box.
7
In the Incident Radiation vs. Normalized Optical Thickness, 1D toolbar, click  Plot.