1 / 45

Summary of Analytical and Computational Work Performed in Development of Thesis

Summary of Analytical and Computational Work Performed in Development of Thesis. Johnathan R. Williams. General Overview. Objective of Thesis was to demonstrate the application of Monte Carlo methods for Radiant Heat Exchange. Three general problems were worked on to demonstrate this.

zahir-melendez
Download Presentation

Summary of Analytical and Computational Work Performed in Development of Thesis

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Summary of Analytical and Computational Work Performed in Development of Thesis Johnathan R. Williams

  2. General Overview • Objective of Thesis was to demonstrate the application of Monte Carlo methods for Radiant Heat Exchange. • Three general problems were worked on to demonstrate this. • The first problem featured use of the Monte Carlo method to determine the configuration factor for Radiation Heat transfer from an elemental area to a circular area. • The second problem featured use of the Monte Carlo method to determine the configuration factor for Radiation Heat Transfer between two areas where shielding was present. • The third problem involved application of the Monte Carlo method for Radiative Heat Exchange in a problem where conduction and convection were present at the boundaries.

  3. Problem #1 Configuration Factor for Radiation Heat Transfer from an Elemental Area to a Circular Area • Published configuration factors already exist for a number of simple geometries. • This simple problem serves as initial validation for the application of the Monte Carlo method. • Configuration Factor was analytically determined using published configuration factor and compared to Monte Carlo solution. • Monte Carlo script was prepared using Visual Basic within Excel.

  4. Problem Geometry & Analytical Solution According to Thermal Radiation Heat Transfer, 4th Edition, by Robert Siegel and John Howell, the Configuration Factor is: It was arbitrarily assumed that r = 4 and h = 10 (units are unimportant as they will cancel, resulting in:

  5. Monte Carlo Geometry

  6. Development of Monte Carlo Code • For an elemental area, initial position is not a factor, thus we have two parameters of concern (q and j). • q is the angle at which the photon bundle is emitted relative to the elemental areas normal. It can vary from 0 to p/2. • j is the cone angle. It can vary from 0 to 2p. • For this problem j is unimportant. As long as q results in a hit, there will be a hit anywhere within the cone angle. • The configuration factor is independent of surface properties. It only depends on geometry. This allows for assumption of blackbody properties, which is computationally the simplest assumption. • Random numbers between 0 and 1 are assigned to the cumulative density function for theta, Rq. • For a blackbody (or a graybody) Rq = sin2q. • q is then calculated using the above relationship and compared to an acceptance standard for a hit.

  7. Acceptance Standard for q 4 4 q 10 10 Acceptance standard, qs is the solution of the following:

  8. Flowchart for Code Execution User inputs number of trials (N), h and r Number of hits & counter initially set to 0 Determine acceptance standard based on r and h Is counter = N? F = Hits/N counter = counter + 1 No No Randomly assign Rq and calculate q Yes Hits = Hits + 1 Is q less than standard?

  9. Obtained Mean and Variability Shown above are means, standard deviations, and variances obtained using 1000, 10000, 100000, and 1000000 trials. The solution for 1000000 trials was obtained using the visual basic macro and also using crystal ball. Recall that the predicted configuration factor was .1379. Note that the standard deviation/variance tends to get smaller with an increasing number of trials.

  10. Obtained Distribution of q for 10,000,000 Trials Note that 13.76907% of the data has q less than an angle of .3805. In other words, the simulated configuration factor is .1376907 for this case.

  11. Fit of Data Trend Line The statistics package in crystal ball reveals that the distribution of q for 1,000,000 trials can be represented by a beta distribution with a min of -.01, a max of 1.58, an alpha of 2.21323 and a beta of 2.21116.

  12. Summary of Data for Problem #1 • Data obtained from Monte Carlo simulation correlated well to predicted configuration factor. • For 1,000,000 trials, the mean configuration factor for 10 observations was .137818 (a percent error of 0.0595%). • The variance predictably drops with increasing number of trials. For 1,000,000 trials, the variance for 10 determinations of the configuration factor is only 9.4192*10-8.

  13. Problem #2 Configuration Factor for Radiation Heat Transfer from One Area to Another With Shielding Present • This is another problem type that can be easily solved analytically. • This serves as additional validation for the application of the Monte Carlo method. • Analytical Solution is via Hottel’s Crossed String Method. • Monte Carlo script was modified to account for additional parameters and hit acceptance criteria.

  14. Problem Geometry Looking for the fraction of thermal radiation leaving A1 that is incident on A2

  15. Analytical Solution – Left Enclosure .7906c Sum of crossed lines = c*4*.7906 Sum of crossed lines = 3.1624c Sum of uncrossed lines = c*(1+.5+2*.7906) Sum of uncrossed lines = 3.0812c .7906c c c/2 .7906c .7906c

  16. Analytical Solution – Right Enclosure .3536c Sum of crossed lines = c*2*1.4142 Sum of crossed lines = 2.8284c Sum of uncrossed lines = c*(1+.5+2*.3536) Sum of uncrossed lines = 2.2072c 1.4142c c/2 c 1.4142c .3536c

  17. Analytical Solution – Combined A1*F1-2 = .5*(sum of crossed lines – sum of uncrossed lines) c*F1-2 = .5*(3.1624c + 2.8284c – 3.0812c – 2.2072c) c*F1-2 = .3512c F1-2 = .3512

  18. Monte Carlo Geometry Initial Position Randomly Chosen r2 = r*cos(j1) As a result of cone angle, photon now collides with the shield

  19. Development of Monte Carlo Code • Now that there are two areas, initial position has become a factor. We now have three parameters of concern (Initial Position, q and j). • q is the angle at which the photon bundle is emitted relative to the elemental areas normal. It can vary from 0 to p/2. • j is the cone angle. It can vary from 0 to 2p. • For this problem j is now important. A non-zero j will cause the distance traveled parallel to A1/A2 for a photon bundle per distance traveled perpendicular to A1/A2 to be less (Reference the proceeding slide). • The configuration factor is independent of surface properties. It only depends on geometry. This allows for assumption of blackbody properties, which is computationally the simplest assumption. • Random numbers between 0 and 1 are assigned to the cumulative density function for initial position, Rpos. • Initial position = c*Rpos. • The acceptance standard for q is calculated based on initial position (the photon bundle can’t hit the shield or leave the system boundary if a hit is to occur). • Random numbers between 0 and 1 are assigned to the cumulative density function for theta, Rq. • For a blackbody (or a graybody) Rq = sin2q. • Random numbers between 0 and 1 are assigned to the cumulative density function for phi, Rj. • For a diffuse material, Rj = 2pj.

  20. Flowchart for Code Execution User inputs number of trials (N) Number of hits & counter initially set to 0 Randomly assign initial position, Rq1 and Rj Determine acceptance standard based initial position Is counter = N? No counter = counter + 1 No F = Hits/N Hits = Hits + 1 Yes Is q2 within the standard range? Calculate q1, j1 and q2 Yes

  21. Results for 1,000,000 trials

  22. Summary of Data for Problem #2 • Data obtained from Monte Carlo simulation correlated well to predicted configuration factor. • For 1,000,000 trials, the mean configuration factor for 10 observations was .3510947 (a percent error of 0.0300%).

  23. Problem #3 Solution of a Problem with Conduction and Convection at the boundaries • Most problems feature combinations of all three modes of heat transfer to varying degrees. • A basic cylinder shape (similar to a combustor or many other flow problems) was analyzed. • A finite difference model was constructed in cylindrical coordinates with conduction and convection. • The finite difference model was validated via comparison to a similar model created in ANSYS. • Next, Radiation was introduced to the finite difference model. • For the first attempt, a macroscopic configuration factor was determined. It was assumed for the first attempt that radiation from each node would uniformly affect the other nodes. • After validating this model, a Monte Carlo solution for the configuration factor from each node to all other nodes was applied. • This is a computationally effective solution for opaque radiation heat transfer.

  24. Problem #3 Parameters Cylinder Geometry Length (L) = 3 ft = 36 in = 91.44 cm Diameter (D) = 1.5 ft = 18 in = 45.72 cm Thickness (t) = 1 in = 2.54 cm Material Properties Assume that the material is a high temp Ni-Alloy Density (r) = 8900 kg/m3 Specific Heat (cp) = 446 J/kg-K Conductivity (k) = 91 W/m-K Emissivity (e) = 0.41 Fluid Properties Inlet Temperature (Tin) = 2500 F = 1644 K Mass Flow Rate ( ) = 80 lbm/s = 36.288 kg/s Density (r) = .2141 kg/m3 Dynamic Viscosity (m) = 57.4x10-6 Pa-s Specific Heat (cp) = 1257 J/kg-K Conductivity (k) = 10.25x10-2 W/m-K

  25. Determination of Film Coefficient Flow is Turbulent

  26. Determination of Film Coefficient (Cont.)

  27. Governing Equations Without Radiation, the governing equations are: q” = 0 at the left boundary q” = -500,000 W/m2 at the right boundary q” = 0 at the lower boundary At the upper boundary: Through the interior: The right boundary condition is a little bit odd, but was intentionally chosen for future model validation, a relatively easy convective solution, and it still allows for a suitable demonstration of the application of the Monte Carlo method.

  28. Solution for Temperatures After discretizing the governing equations, the following system of equations is obtained: [A][T] = [C] Where: [A] is a coefficient matrix for nodal temperatures [T] is the matrix containing the nodal temperatures [C] is the results matrix (mainly populated with zeros or heat flux where appropriate) This system of equations is linear and can readily be solved either explicitly or implicitly.

  29. Results of Finite Difference Model

  30. Results of Finite Difference Model

  31. Validation of Finite Difference Model

  32. Results of Finite Difference Model (English Units)

  33. Results of Finite Difference Model (English Units)

  34. Extension of Radiation to Model • For initial extension of Radiation Heat Transfer to finite difference model, it is assumed that Radiative Heat Transfer only involves opaque surface to surface exchange. • A macroscopic configuration factor for exchange along the entire cylindrical area to the same cylindrical area is determined analytically. • The configuration factor is applied to each node and it is assumed that each node will uniformly distribute thermal radiation to all other nodes (including itself). • This is used to construct a coefficient matrix for Radiation Heat Transfer. • Because Radiation Heat Transfer is a function of the difference of temperature to the fourth power, the governing equation for the finite difference model is now non-linear and requires an iterative solution. A successive under-relaxation technique is employed.

  35. Determination of Configuration Factor Definitions: R = r/a; X = (2R2+1)/R2 Everything else leaving surface 1 impacts the cylinder wall! Calling the cylinder wall surface 3 F1-3 = 1-F1-2 = 1-0.055728 = 0.944272

  36. Determination of Configuration Factor (Cont.) A1F1-3 = A3F3-1 F3-1 = (A1F1-3)/A3 A1 = pD2/4 = p*(18in)2/4 = 254.47in2 A3 = pDL = p*(18in)*(36in) = 2035.8in2 F3-1 = (254.47in2*0.944272)/2035.8in2 = .11803 F3-2 = F3-1 = .11803 F3-3 = 1-F3-1-F3-2 = 1-2*.11803 = .76394 For this model, there are 144 axial divisions (145 nodes) Assuming that each node radiates uniformly to each other node, the configuration factor for node x to node y is: Fx-y = .76394/145 = .005269

  37. Governing Equations • The same boundary conditions and governing equations that were applied for the No Radiation case are applied. The only difference is that along the upper boundary, we also have a Radiation term. • The governing equation for Radiative Exchange between node x and node y is: • q”r = seFx-y(Tx4 – Ty4) • Where: • = Stefan-Boltzman constant (5.67x10-8 W/m2K4) • = Emissivity • Fx-y = The configuration factor for node x to y

  38. Solution for Temperatures After discretizing, the system of equations now becomes: [A][T] + [B][T4] = [C] The difference being a T4 term and its associated coefficient matrix To solve this, the system can be rewritten as: [A*][T] = [C] Where: [A*] = [A + BT3] This system of equations is no longer linear and must be solved implicitly. A SUR subroutine is employed.

  39. Results of Finite Difference Model (English Units)

  40. Results of Finite Difference Model (English Units)

  41. Development of Monte Carlo Code • We again have three parameters of concern (Initial Position, q and j). • q is the angle at which the photon bundle is emitted relative to the elemental areas normal. It can vary from 0 to p/2. • j is the cone angle. It can vary from 0 to 2p. • The configuration factor is independent of surface properties. It only depends on geometry. This allows for assumption of blackbody properties, which is computationally the simplest assumption. • Random numbers between 0 and 1 are assigned to the cumulative density function for initial position, Rpos for an interval around each node. • Initial position = c*Rpos + the Initial Position of the Node Interval. • Random numbers between 0 and 1 are assigned to the cumulative density function for theta, Rq. • For a blackbody (or a graybody) Rq = sin2q. • Random numbers between 0 and 1 are assigned to the cumulative density function for phi, Rj. • For a diffuse material, Rj = 2pj. • What makes this problem tricky is that the cylinder wall has curvature (it’s circular). • We know that if q is 90 degrees and j is 0 or 180 degrees, a photon bundle will pass without hitting the cylinder. Otherwise, if q is 90 degrees, it will immediately impact the cylinder wall (at the node location of emission). If q is 0 degrees, it will impact the cylinder wall (directly across from the initial position at the node location of emission). If j is 90 degrees or 270 degrees, it will impact the node location of emission.

  42. Determination of Impact Location • Unique results for special angles of q and j were discussed on the proceeding slide. • For all other angles, the photon impact location must be identified. • The location of photon emission can be treated as bottom dead center. • q1 determines the “height” of photon travel per unit length traveled. • j1 determines the distance the photon bundle travels radially per unit length traveled. • The photon impacts the cylinder wall when the radial location of the photon meets the radial location of the cylinder wall at that height.

  43. Determination of Impact Location Emission Location & Origin Equation of circle is: We now have two simultaneous equations to solve: From l, we know the node that is hit!! From the fraction of photon bundles leaving node x and impacting node y, we have the configuration factor for x to y. This permits integration into the finite difference model.

  44. Results of Finite Difference Model (English Units)

  45. Results of Finite Difference Model (English Units)

More Related