Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UF00090205/00001
## Material Information- Title:
- Tomographic two-phase flow measurement using Compton scattering of gamma rays
- Creator:
- Bodette, David E. (
*Dissertant*) Jacobs, Alan M. (*Thesis advisor*) Carroll, Edward E. (*Reviewer*) Dugan, Edward T. (*Reviewer*) Hanrahan, Robert J. (*Reviewer*) Roan, Vernon P. (*Reviewer*) - Place of Publication:
- Gainesville, Fla.
- Publisher:
- University of Florida
- Publication Date:
- 1986
- Copyright Date:
- 1986
- Language:
- English
## Subjects- Subjects / Keywords:
- Annular flow ( jstor )
Density ( jstor ) Gamma ray spectrum ( jstor ) Gamma rays ( jstor ) Inverted spectra ( jstor ) Noise spectra ( jstor ) Photons ( jstor ) Pixels ( jstor ) Spectral energy distribution ( jstor ) Spectral reconnaissance ( jstor ) Compton effect ( lcsh ) Dissertations, Academic -- UF -- Nuclear Engineering Sciences Gamma rays ( lcsh ) Nuclear Engineering Sciences thesis Ph.D Tomography ( lcsh ) Two phase flow ( lcsh ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Abstract:
- Using Compton scattered gamma-rays to measure local void fraction was first suggested by Kondic and Hahn in 1970. The Compton scattered gamma-ray densitometers they suggested employ a single, narrow beam source and either a well collimated detector or an uncollimated detector. The collimated detector configuration only gives the local void fraction measurement in the small volume where the source and detector collimators intersect. The uncollimated detector configuration is a more efficient design since the local void fraction along the source beamâ€™s path is measured in a single reading. A logical extension of the technique is to use wide beam illumination and uncollimated detectors to sample an even greater portion of the flow cross section. This report investigates and demonstrates several methods of inferring two-phase flow parameters using wide beam illumination coupled with two detectors placed symmetrically about he source and pipe. The spatial distribution of the fluid in a slice of the pipe is encoded with respect to energy in the singly scatter photon flux. Two basic techniques are detailed for decoding the spatial information: the method of spectral moments and the method of computed tomography (CT). Examination of the low-order moments of the singly scattered photon spectra from the two detectors provides sufficient information for flow regime identification. Flow asymmetries are revealed by comparison of the spectral moments of the two measured spectra. Further classification of the flow pattern is made on the basis of the first and second moments of the spectra. The real focus of this report is the adaptation and successful demonstration of CT techniques with Compton scattering. Three series expansion techniques are considered: the algebraic reconstruction technique (ART), the simultaneous iterative reconstruction technique (SIRT), and the iterative least squares (ILS) technique. Application of the modified ART, SIRT, and ILS algorithms to a hot-spot and a cold-spot reconstruction problem indicates that SIRT and ILS are more accurate than ART. A more extensive testing of the ILS algorithm for a variety of model flow regimes demonstrates the potential of Compton scatter tomography as a quantitative measurement technique.
- General Note:
- Typescript.
- General Note:
- Vita.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 1986.
- Bibliography:
- Bibliographies : leaves 191-195.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 030321815 ( alephbibnum )
16656185 ( oclc )
## UFDC Membership |

Downloads |

## This item has the following downloads:
Binder1 ( .pdf )
UF00090205_00001.pdf 00006.txt 00199.txt 00206.txt 00026.txt 00047.txt 00080.txt 00058.txt 00105.txt 00060.txt 00054.txt 00092.txt 00051.txt 00177.txt 00055.txt 00061.txt 00153.txt 00162.txt 00137.txt 00205.txt 00183.txt 00067.txt 00142.txt 00181.txt 00037.txt 00033.txt 00100.txt 00096.txt 00145.txt Binder1_pdf.txt 00108.txt 00174.txt 00062.txt 00002.txt 00112.txt 00146.txt 00076.txt 00057.txt 00148.txt 00182.txt 00158.txt 00087.txt 00066.txt 00186.txt 00073.txt 00075.txt 00194.txt 00007.txt 00127.txt 00027.txt 00063.txt 00114.txt 00091.txt 00071.txt 00120.txt 00059.txt 00136.txt 00150.txt 00042.txt 00012.txt 00201.txt 00156.txt 00125.txt 00023.txt 00167.txt 00039.txt 00122.txt 00163.txt 00133.txt 00210.txt 00072.txt 00081.txt 00020.txt 00038.txt 00188.txt 00179.txt 00193.txt 00151.txt 00101.txt 00011.txt 00190.txt 00160.txt 00034.txt 00010.txt 00083.txt 00157.txt 00143.txt 00024.txt 00110.txt 00093.txt 00117.txt 00152.txt 00184.txt 00022.txt 00204.txt 00119.txt 00189.txt 00168.txt 00111.txt 00154.txt 00207.txt 00019.txt 00203.txt 00126.txt 00135.txt 00172.txt 00191.txt 00170.txt 00169.txt 00070.txt 00032.txt 00138.txt 00068.txt 00107.txt 00128.txt 00140.txt 00064.txt 00008.txt 00035.txt 00095.txt 00200.txt 00090.txt 00196.txt 00016.txt 00116.txt 00118.txt 00005.txt 00103.txt 00208.txt 00166.txt 00197.txt 00017.txt 00139.txt 00178.txt 00097.txt 00050.txt 00121.txt 00085.txt 00195.txt 00018.txt 00098.txt Copyright.txt 00209.txt 00113.txt 00052.txt 00144.txt 00084.txt 00069.txt 00134.txt 00004.txt 00088.txt 00187.txt 00029.txt 00175.txt 00074.txt 00132.txt 00077.txt 00041.txt 00053.txt 00164.txt 00198.txt 00104.txt 00185.txt 00115.txt 00078.txt 00149.txt 00141.txt 00131.txt 00021.txt 00028.txt 00031.txt 00009.txt 00046.txt 00147.txt 00044.txt 00013.txt 00001.txt 00109.txt 00099.txt 00102.txt 00180.txt 00040.txt 00129.txt 00094.txt 00159.txt 00014.txt 00086.txt 00130.txt 00049.txt 00079.txt 00048.txt 00165.txt 00123.txt 00065.txt 00106.txt 00015.txt 00056.txt 00192.txt 00045.txt 00161.txt 00171.txt 00176.txt 00173.txt 00202.txt 00030.txt UF00090205_00001_pdf.txt 00089.txt 00082.txt 00155.txt 00036.txt 00124.txt 00043.txt 00025.txt 00003.txt |

Full Text |

TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT USING COMPTON SCATTERING OF GAMMA RAYS By DAVID E. BODETTE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1986 Copyright 1986 by David E. Bodette ACKNOWLEDGMENTS I thank Dr. Alan M. Jacobs, my adviser, for his patient guidance and careful attention to detail. This work was partially supported by a Liquid Metal Fast Breeder Reactor Technology Research Award from the Argonne National Laboratory under contract with the United States Department of Energy (DOE). The author was also supported by the DOE Energy Graduate Traineeship Program. iii TABLE OF CONTENTS ACKNOWLEDGMENTS . . . . . . . . . . LIST OF TABLES . . . . . . . . . . LIST OF FIGURES . . . . . . . . . ABSTRACT . . . . . . . . . . . . CHAPTERS 1 INTRODUCTION . . . . . . . . 2 THE FORWARD PROBLEM . . . . . . . Collimated Point Source/Collimated Detector . Collimated Point Source/Uncollimated Detector Uncollimated Point Source/Uncollimated Detector Parallel Beam Source/Uncollimated Detector. . Discrete Forward Equations 3 METHOD OF SPECTRAL MOMENTS. . . . Description of Method . . . . . Procedure and Results . . . . . Interpreting the Moments. . . . A Detailed Example . . . . Using the Moments . . . . . 4 THE RECONSTRUCTION PROBLEM. . . . General Solution Methods . . . . A Closer Look . . . . . . . Backprojection . . . . . . Monte-Carlo Backprojection . . . Matrix Methods . . . . Algebraic Reconstruction Technique Simultaneous Iterative Reconstruction Technique Iterative Least-Squares Technique . Comparison of Reconstruction Techniques 4 COMPTON SCATTER TOMOGRAPHY . . . An Example . . . . . . . Accuracy . . . . . . . Resolution . . . . . . . Image Contrast . . . . . . . . 35 . . 35 . . 42 . . 44 . . 46 . . 52 . . 54 . . 56 . . 58 . . 59 . . 64 . . 69 . . 72 . . 76 . . 77 . . 79 . . 90 . . 90 . . 92 . . 106 . . 109 Page iii vi vii xi 1 10 11 15 17 28 . . . . . 29 6 SUMMARY AND CONCLUSIONS . . . . . Summary . . . . . . . . . A Practical Application . . . . . Conclusions and Future Research . . . APPENDICES I MONTE-CARLO TRANSPORT PROGRAM . . . . II INSIGHTS AND FINER POINTS . . . . . Parallel Versus Diverging Source Beams. . Multiple Scatter Versus Detector Position . Source Placement . . . . . . Convergence Criterion . . . . . . III RECONSTRUCTION COMPUTER PROGRAMS .. Running the Programs . . . . . . Flow of Programs . . . . . . . IV MAXIMUM ENTROPY METHOD . . . . . V DATA . . . . . . . . . . REFERENCES . . . . . . . . . . BIBLIOGRAPHY . . . . . . . . . . BIOGRAPHICAL SKETCH . . . . . . . . Page 118 119 121 126 130 133 133 137 142 143 S. 144 S. 144 . 145 S. 146 S. 149 S. 191 S. 195 S. 196 LIST OF TABLES Page 2.1 Multiply scattered photon spectra . . .. 22 3.1 Spectral moment results . . .. . . 45 4.1 Data for hot-spot and cold-spot tests ... . 83 5.1 Test parameters and values . . .. . . 98 5.2 Summary of reconstruction results. . . . 105 5.3 Parameters affecting image resolution . . 107 LIST OF FIGURES Page 2.1 Collimated source/collimated detector configuration 12 2.2 Collimated source/uncollimated detector configuration 16 2.3 Uncollimated source/uncollimated detector configuration 19 2.4 Illustration of scattering angle blur 25 2.5 Illustration of an isogonic segment 27 2.6 Parallel beam source/uncollimated detector configuration 30 2.7 Overlaying of pipe cross section with an n x n pixel grid 30 2.8 Feathering of isogonic segments 34 3.1 Compton scatter densitometry configuration 36 3.2 Spectral moments of a single-scattered spectrum 39 3.3 Source-pipe-detector configuration 44 3.4 Illumination of stratified flow regime 44 3.5 Pipe illumination for various beam widths 48 3.6 Spectrum of upper detector versus beam width 49 3.7 Spectrum of lower detector versus beam width 50 4.1 Conventional CT scanner 60 4.2 Monte-Carlo backprojection 67 4.3 Spectra from Monte-Carlo backprojection 68 4.4 Flow chart of reconstruction algorithms 71 4.5 Example of an ART solution for two pixel image 75 4.6 Hot-spot and cold-spot test configuration 80 4.7 Hot-spot spectra 81 4.8 Cold-spot spectra 82 vii Page 4.9 Hot-spot reconstruction 84 4.10 Cold-spot reconstruction 85 4.11 Convergence of reconstruction algorithms 87 5.1 Calculated spectrum using incorrect source strength value 91 5.2 Reconstruction using incorrect source strength value 93 5.3 Experimental setup 95 5.4 Flow regime models 99 5.5 Singly scattered photon spectrum 100 5.6 Reconstruction with and without noise 101 5.7 Spectra with noise added 103 5.8 Detected spectrum with pipe wall effects 117 6.1 Spectrum from a 4 inch pipe 122 6.2 Detail of spectrum from the pipe 123 6.3 Effect of signal-to-noise ratio 125 6.4 Detector response versus bubble size 127 I.1 Example geometry of photon trajectory 132 II.1 Detector response for diverging source beam 135 11.2 Detector response for parallel source beam 136 11.3 Detector configuration for multiple scatter versus detector position 138 11.4 Spectra for multiple scatter versus detector position 139 V.1 Spectra from annular flow 151 V.2 Spectra from annular flow 152 V.3 Spectra from annular flow 153 V.4 Spectra from inverted annular flow 154 viii Page V.5 Spectra from inverted annular flow 155 V.6 Spectra from inverted annular flow 156 V.7 Spectra from bubbly flow 157 V.8 Spectra from bubbly flow 158 V.9 Spectra from bubbly flow 159 V.10 Spectra from stratified symmetric flow 160 V.11 Spectra from stratified symmetric flow 161 V.12 Spectra from stratified symmetric flow 162 V.13 Spectra from stratified asymmetric flow 163 V.14 Spectra from stratified asymmetric flow 164 V.15 Spectra from stratified asymmetric flow 165 V.16 Annular flow pattern 166 V.17 Reconstruction of annular flow 167 V.18 Reconstruction of annular flow 168 V.19 Reconstruction of annular flow 169 V.20 Inverted annular flow pattern 170 V.21 Reconstruction of inverted annular flow 171 V.22 Bubbly flow pattern 172 V.23 Reconstruction of bubbly flow 173 V.24 Reconstruction of bubbly flow 174 V.25 Reconstruction of bubbly flow 175 V.26 Stratified symmetric flow pattern 176 V.27 Reconstruction of stratified symmetric flow 177 V.28 Stratified asymmetric flow pattern 178 V.29 Reconstruction stratified asymmetric flow 179 V.30 Example of smooth convergence 180 Page V.31 Example of divergence 181 V.32 Example of divergence 182 V.33 Example of oscillating divergence 183 V.34 Example of non-convergence 184 V.35 Setup for laboratory experiments 185 V.36 Experimental spectrum of annular flow 186 V.37 Experimental and Monte-Carlo spectra 188 V.38 Experimental spectrum of stratified flow 189 V.39 Experimental spectrum of bubbly flow 190 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT USING COMPTON SCATTERING OF GAMMA RAYS By DAVID E. BODETTE December, 1986 Chairman: Dr. Alan M. Jacobs Major Department: Nuclear Engineering Sciences Using Compton scattered gamma-rays to measure local void fraction was first suggested by Kondic and Hahn in 1970. The Compton scattered gamma-ray densitometers they suggested employ a single, narrow beam source and either a well collimated detector or an uncollimated detector. The collimated detector configuration only gives the local void fraction measurement in the small volume where the source and detector collimators intersect. The uncollimated detector configuration is a more efficient design since the local void fraction along the source beam's path is measured in a single reading. A logical extension of the technique is to use wide beam illumination and uncollimated detectors to sample an even greater portion of the flow cross section. This report investigates and demonstrates several methods of inferring two-phase flow parameters using wide beam illumination coupled with two detectors placed symmetrically about the source and pipe. The spatial distribution of the fluid in a slice of the pipe is encoded with respect to energy in the singly scattered photon flux. Two basic techniques are detailed for decoding the spatial information: the method of spectral moments and the method of computed tomography (CT). Examination of the low-order moments of the singly scattered photon spectra from the two detectors provides sufficient information for flow regime identification. Flow asymmetries are revealed by comparison of the spectral moments of the two measured spectra. Further classification of the flow pattern is made on the basis of the first and second moments of the spectra. The real focus of this report is the adaptation and successful demonstration of CT techniques with Compton scattering. Three series expansion techniques are considered: the algebraic reconstruction technique (ART), the simultaneous iterative reconstruction technique (SIRT), and the iterative least squares (ILS) technique. Application of the modified ART, SIRT, and ILS algorithms to a hot-spot and a cold-spot reconstruction problem indicates that SIRT and ILS are more accurate than ART. A more exten- sive testing of the ILS algorithm for a variety of model flow regimes demonstrates the potential of Compton scatter tomography as a quantitative measurement technique. xii CHAPTER 1 INTRODUCTION The case for developing improved non-intrusive measur- ing technologies for sensing material/phase distribution in two-phase flow has been well argued by many authors. The problem is one which spans a broad number of fields inclu- ding manufacturing, engineering, security, and medicine. Each potential application imposes its own peculiar demands, which are often at odds with constraints in other applica- tions. This situation has resulted in many diverse techniques being pursued and often optimized for each specific problem. Examples of the technologies employed include ultrasonics, microwaves, lasers, optics, and irra- diation with particles, X rays, or gamma rays. One specific problem encountered in the nuclear power industry and its associated research industry is to measure the void fraction and/or liquid distribution (flow regime) of two-phase flows inside steel walled pipes. Any intrusive instrument which disturbs the flow pattern is certainly not desirable, but an even more strict requirement often encoun- tered in these circumstances is that penetrations through the pipe wall are not permitted. Another requirement of any measuring instrument is that its output should be unambi- guous. In the case of two-phase flow void fraction measurement this means the entire flow cross section must be interrogated. There may also be physical limitations in terms of size and orientation, as well as access to the pipe. A perusal of the recent literature on measurements in two-phase flow [1, 2, 3, 4] reveals that the choice of techniques for making non-intrusive local void fraction and/or phase distribution measurements is limited to irra- diation with particles, gamma rays, and laser light. When intrusion is not an issue, a variety of probe techniques are available for local property measurements. Reference [4] lists four basic probe methods which appear to be reliable enough for general use: electrical probes, optical probes, thermal anemometers, and micro-thermocouples. Laser scattering and interferometry is suggested as a very promising method when the flow boundary is not opaque. There are techniques for making average void fraction measurements which might be adaptable to the local void fraction problem. Ultrasonics is capable of measuring the average void fraction in very thick walled pipes, but only for low void fractions (i.e. < 0.20) [3]. An array of ultrasonic sensors around a pipe might produce an image of the flow provided the voids are not too big or too many. Microwave resonant cavities are also capable of measuring average void fraction. Using microwaves to sense local void fraction would probably run into difficulties from the metal pipes [1]. Nuclear magnetic resonance (NMR) is currently receiving a lot of attention in the medical field and certainly has potential for two-phase flow. The medical NMR machines involve a radio-frequency field and a large magnetic field which cause problems in nearby electronic equipment. The metal pipes seem to be another source of potential problems in implementing NMR for local measurements; however, plastic pipes can be used to avoid such problems [5]. Gamma-ray attenuation is the predominant method of choice for making average void fraction measurements and is often used as the calibration standard for other methods. Gamma-ray attenuation instrumentation and techniques are recognized as having reached the highest level of development. The reason for this dominance is that the method places a minimum of constraints on the user. The method is non-intrusive, fairly insensitive to flow pattern, and functions over a broad range of void fractions, flow regimes, and temperatures. There are many ways of utilizing gamma rays or neutrons to measure local void fraction and phase distribution. The basic methods of implementing tomographic two-phase flow measurement with gamma rays or neutrons can be grouped into the three following categories: 1) Measure the emission of gamma rays or neutrons from- activation products in the fluid, induced excitation, or injected radiotracers. 2) Measure the attenuation of narrow or wide source beams in traversing the pipe. 3) Measure the singly scattered gamma rays or neutrons from a primary source beam. Methods from category 1 could be roughly classified as emission computed tomography (ECT) techniques. An example of a category 1 method is a single photon ECT system which examines the ~6 MeV gamma rays from the activation-decay scheme pair 16(n,p)N16 and N16 O16 + y + 8. Although recent advances in the analytical techniques of ECT and in instrumentation suggest this as a viable method, it appears that only cursory attention has been given to the technique [6]. An ECT system could also be employed using injected radiotracers; however, the quantity of radio- nuclides required poses a serious drawback [7]. The most familiar tomographic instruments comprise category 2 and are generally referred to as computed axial tomography (CAT) scanners. The basic hardware setup of CAT scanners involves placing a source on one side of the pipe and a detector on the opposite side of the pipe. The whole source-detector assembly rotates in a controlled manner around the pipe while the attenuation of the source beam at different orientations is recorded. The rate at which CAT scanners can rotate about the pipe is a limiting factor. CAT scanners of the newest generation do not employ any moving parts, but instead use a magnetically focused flying spot electron beam and a fixed circular detector array [8]. These fifth generation CAT scanners can acquire an image in approximately 1/10 second and are very promising where cine tomography is necessary. The capital investment in these new CAT scanners is significant and combined with the fact that they are not compact or portable means they are not a viable alternative in many applications. Category 3 techniques are closely related to ECT except that the emitted radiation does not come from radioactive material within the fluid, but from radiation scattered out of an external illuminating source beam. Gamma-ray and neutron sources are possible candidates, and methods based on both are actively being developed [9]. The earliest reported efforts with scattered radiation focused on gamma rays scattering and involved measuring the local density in a small volume [10]. A well collimated source illuminated an object and a well collimated detector measured the Compton scattered gamma rays coming from the small volume where the source and detector collimators' fields-of-view intersected. The output of the detector was proportional to the electron density in the small volume which in turn is proportional to the material density for a constant material composition. The Compton scattering method of local density measurement developed, at first, independently in the medical and industrial imaging fields. Lale published two articles, one in 1959 [11] and another in 1968 [12], which demonstrated how internal organs could be imaged by exploiting the Compton scattering of gamma rays. In 1970, Kondic and Hahn [13] published the first suggested engineering application of Compton scattering, that of making local density measurements in two-phase flow. They proposed both a collimated source/collimated detector (as Lale) and a collimated source/uncollimated detector arrange- ment. The uncollimated detector offered the improvement that local density values were given for points along the illuminated path of the source beam from a single measurement. Acting independently, Farmer and Collins [14, 15] proposed and demonstrated the collimated source/uncolli- mated detector arrangement in determining anatomical cross sections. Their images were only qualitative since the calibration procedures used failed to completely remove the effects of multiple scattering or variations in attenuation. The analytic method for making quantitative measurements using this source-detector arrangement by deconvolving the substantial multiple scatter blur was developed by Anghaie [16, 17]. The use of Compton scattering to make local density measurements in two-phase flows was pursued earnestly in the late 70's and in to the 80's. Almost exclusively the focus was on the collimated source and either collimated or uncol- limated detector arrangement. Tomographic images produced with either of these arrangements required scanning the source and/or detector to cover all the points in the plane of interest. Kondic [18] in 1978 first proposed using an uncollimated or wide beam source coupled with an uncolli- mated detector to image an entire flow cross section. Such an arrangement made more efficient use of the source photons than the previous source-detector arrangements and required no moving parts. Both of these attributes made the arrange- ment more adapted to the high framing rates necessary for transient analysis. However, no critical study or viable implementation of the wide source beam/uncollimated detector arrangement has appeared in the literature. Imaging with scattered neutrons was developed in the late 70's and early 80's independently of and somewhat concur- rently of the present work with scattered gamma rays. The use of neutrons to make local void fraction measurements was suggested by Banerjee and Lahey [19] and separately by other researchers. Neutron scattering tomography was developed by Hussein [20] as his Ph.D. thesis. He used a wide parallel beam of 14 MeV neutrons to illuminate a simulated two-phase flow and then measured the fluence and energy of the scattered neutrons using small organic scinti- lator detectors at different positions around the flow volume. The focus of the present dissertation is the develop- ment and demonstration of the algorithms for reconstructing the local density distribution of a two-phase flow using Compton scattered gamma rays. The specific experimental setup considered uses a wide-beam collimated gamma-ray source to illuminate the flow and energy sensitive detectors to examine the emerging singly Compton scattered photon flux. To date, the reconstruction techniques proposed for this setup have been based on deterministic algebraic solu- tions and have not addressed the problems of multiple scat- tering, statistical error, image blur, and other errors. The contribution of this author is to apply the powerful reconstruction algorithms developed in the medical and optical imaging fields to this problem and to critically assay the effects of__ multiple scattering and statistical errors on the reconstruction. Tomographic imaging with Compton scattered gamma rays has two distinct features which distinguish it from previously encountered tomographic reconstruction problems. First, the imaging processes dealt with in the medical and optic fields have been linear or could be linearized by an analytic technique. Imaging the phase distribution in a two-phase flow with Compton scattered gamma rays is a non- linear problem. Second, the combination of using a point source and the Compton scattering process introduces a unique curved geometry to the flow reconstruction process. The forward problem of predicting the emerging singly Compton scattered photon flux given the local density dis- tribution is discussed in Chapter 2. Chapter 3 details a method of inferring tomographic information about a two- phase flow by examining the spectral moments of the Compton scattered photons [21]. Chapter 4 develops several more sophisticated tomographic reconstruction techniques and compares their performance under ideal conditions using data generated by Monte-Carlo calculations. Chapter 5 evaluates Compton scatter tomography under more realistic conditions and discusses some of the considerations in designing a system. Chapter 6 summarizes the significant results of this study and offers directions for future work. CHAPTER 2 THE FORWARD PROBLEM Fundamental to the problem of reconstructing the phase or density distribution using Compton scattered gamma rays is being able to model the gamma ray interaction with the density field. The aim of this chapter is to describe the basic equations predicting the emerging scattered photon flux from a known density distribution illuminated by a point gamma-ray source. The basic concepts of using gamma- ray scattering to make local density measurements are introduced first with the collimated point source/collimated detector arrangement. Once the equations for analyzing this setup are developed, the analysis is extended to collimated point source/uncollimated detector arrangements by introduc- ing the equivalent approach of energy collimation. The analysis is extended further to the wide-beam collimated point source/uncollimated detector arrangement. Lastly, parallel source beam illumination is considered and the differences between it and the diverging wide-beam illumination are discussed. Throughout, a detailed analysis is only given to the singly scattered photon flux. The presence of the multiply scattered photon flux is included in all the equations and discussed in general terms. However, for the case of the wide-beam collimated source /uncollimated detector configuration a specific method is given for removing the multiple scatter component. Collimated Point Source/Collimated Detector Figure 2.1 illustrates the arrangement of a collimated source and collimated detector for making local density measurements. Monoenergetic gamma rays leave the source and are collimated to a thin pencil beam which illuminates a chord of the pipe. A similar collimator on the detector is aligned so that its field-of-view intersects the source beam at the small volume denoted by V. The detector senses the source gamma rays scattered in the volume V. The detector response, D(E), is D(E) = ESOPwAiAoAZ + M(E) (2.1) where S p pw Az Ai and Ao M(E) is the source activity; is the density of the fluid in the volume V; is the incoherent differential scattering cross section; is the path length of source photons in V; are the total attenuation functions, including interaction attenuation and the solid angles of view, along the incoming and outgoing paths respectively; is the detector efficiency; and represents the contribution to the spectrum of multiply scattered photons. Pipe Collimator Source Detector Figure 2.1: A collimated source/collimated detector gamma- ray scattering configuration for measuring the density in the volume V. The functional form of both attenuation functions is Ax = Cx*exp(-f pxp(z)dz) (2.1a) zx where Cx contains the geometric or field-of-view factors, P, is the mass attenuation coefficient, p(z) is the local density at z, and x is a dummy subscript representing i and o for the incoming and outgoing paths, respectively. The line integral in equation (2.la) is carried out either along the incoming path, ki, from the source to the volume V or along the outgoing path, o from V to the detector. The value of the integral is in units of optical path length. The predominant interaction of gamma rays with energies between 0.1 and 1.0 MeV in water is Compton scattering. The energy of the singly Compton scattered photons arriving at the detector is related to the polar scattering angle, e, by the kinematic relation E = Ei/(l + (Ei/mec2)(1 cos 8)) (2.2) where Eo is the energy of the scattered photon, Ei is the energy of the source photon, and mec2 is the rest mass energy of an electron. The singly scattered photons, represented by the first expression on the right hand side (RHS) of equation (2.1), all arrive at the detector after scattering through approximately the same angle and thus have approximately the same energy. Assuming the attenuation along the incoming and outgoing paths remains constant in time, the detector response to the singly scattered photons is proportional to the density in the small volume V. The multiply scattered photons arriving at the detector, unlike the singly scattered photons, have a broad range of energies since they can suffer any number of collisions (with concomitant changes in energy) prior to reaching the detector. The functional form of the multiply scattered photon component of the detector response is M(E)= ffff ('(r ,E', ')K(r',E',Q-;rd,E,Q)dr'dE'dQ'drdd ( (2.1b) where M(..)is the angular photon flux and K(..) is the transport kernel. The kernel can be viewed as the probabi- lity that a photon at r going in direction n' with an energy E' suffers a scattering interaction at r'; the scat- tered photon goes in direction Q toward the detector with an energy E; the scattered photon reaches point rd in the detector's active volume, interacts with the detector and is counted. The integral with respect to r is carried out only over the detector collimator's field-of-view. The detector collimator serves not only to pick out the volume of interest, V, but to reduce the magnitude of the multiply scattered photon spectral component by restricting the region of space where the photons can suffer their last collision prior to reaching the detector. Collimated Point Source/Uncollimated Detector The function of the detector collimator in picking out the small volume V can also be achieved by electronic colli- mation. An example is to use a high energy resolution detector coupled to a multi-channel analyzer (MCA) instead of the detector collimator. The point of scattering of a singly scattered photon is related to its energy by equa- tion (2.2) and the polar scattering angle as determined from the measurement geometry. The singly scattered photons arriving at the detector now have a range of energies. The detector response becomes D(E) = Sop(x)pwAi(x)Ao(x)e(E) + M(E) (2.3) where S, p, 1,, Ai, Ao, and M are as defined for equa- tion (2.1) and x is the point of scattering. The x-axis is colinear with the source beam illumination path (Figure 2.2). Note that x and e can be expressed in terms of E using equation (2.2) as 2 Yd x = xd - 1 cos2e (2.4) cos e = 1 + (mec2/Ei)(l Ei/E) Pipe Collimator Source- - x- d Detector Figure 2.2: Configuration of a collimated source/uncollimated detector gamma-ray scattering densitometer. where xd and yd are the detector coordinates illustrated in Figure 2.2. The first term on the RHS of equation (2.3) is strictly a function of a single variable for a fixed geometry. The multiple scattering term in equation (2.3) has the same functional form as equation (2.1) except that the inte- gration with respect to r is over all space since there is no detector collimator. Therefore, the magnitude of the multiply scattered photon component of the spectrum is greater with an uncollimated detector than with a collimated detector. Offsetting this increased measurement structured noise is the fact that the singly scattered spectrum now contains the local density information for the entire source beam illumination path. The information which would have required many measurements with a collimated detector is available from a single measurement with an uncollimated detector. Uncollimated Point Source/Uncollimated Detector Just as eliminating the detector collimator increases the information contained in a single measurement, so does eliminating the narrow source collimator. An uncollimated or a wide angle-of-view collimated source can illuminate the entire flow cross section. The detector response then contains the local density information from every point in the flow cross section. The manner in which the local density information is encoded in the detector response, however, complicates the data processing necessary to reconstruct the local density distribution. Figure 2.3 illustrates a wide beam source Compton scat- tering densitometer setup. The curves drawn across the pipe represent lines of constant scattering angle, called isogonic lines by Kondic [18]. Primary source photons scattering at points along an isogonic line all scatter through the same angle in order to reach the detector. Using elementary geometry it is easy to show that the isogonic lines are all segments of circles which pass through both the source and the detector. The axis labeled Y in Figure 2.3 is the loci of all the centers of these circles. The radius, r, and the location of the center- point, yc, of each isogonic circle as a function of the scattering angle, 8, are r = SD sin 8 2 (2.5) y SD cot e C- 2 where SD is the distance between the source and the detector. The detector response is now represented by an integral equation D(E)= f D(x,y)dk + M(E) (2.6) Z(E) Pipe Source A TOP VIEW ELEVATION VIEW Figure 2.3: Configuration of an uncollimated source/ uncollimated detector gamma-ray scattering system for tomography. where M(E) is the multiply scattered photon spectrum, D(x,y) is the contribution to the detected photon spectrum of singly Compton scattered photons scattered at point (x,y), and dZ is the differential line element along the isogonic contour (E). The first term on the RHS of equation (2.6) is essentially the integral of the first term on the RHS of equation (2.3). The second term, as before, represents the multiply scattered photon contribution to the spectrum and has the same functional form as equation (2.1b). The multiple scatter component can comprise approxi- mately one-third of the total signal [22, 23] and must be removed. It may be possible to remove the multiple scatter- ing component analytically since the transport kernel is a known function and the angular photon flux can be generated from an order-of-scattering calculation once a good estimate of the density distribution is obtained. However, this would be a cumbersome process and does not appear to be necessary as there appears to be an acceptable alternative approximation. A multiply scattered photon arriving at the detector with energy E could have come from any point within a large volume of the pipe. The structure of the fluid distribution imposed on the multiply scattered photon spectrum would tend to be averaged out because of the many possible paths that the photon could have taken prior to arriving at the detector with energy E. The result is that for wide beam source/uncollimated detector arrangements the multiply scattered photon spectrum should be smooth and a simple background subtraction procedure can be employed to remove it. The validity of this intuitive argument has been tested by a series of Monte-Carlo calculations (see Appendix I). The multiply scattered photon spectrums have been com- puted for five flow regimes at both a high and low void fraction. The calculations assume that the flow is contained in a pipe whose thickness and composition are such that it can be neglected. The results, summarized pictorially in Table 2.1, show that to a first order approx- imation equation (2.6) can be rewritten as D(E)= J D(x,y)dk + B(E). (2.6*) R(E) where B(E) is some function of energy which is fitted to the smooth multiply scattered photon spectrum, and D(x,y) is the contribution to the detected photon spectrum of singly Compton scattered photons scattered at point (x,y). The function D(x,y) can be written more explicitly as D(x,y)= SO (E)pw(Ei,Eo)p(x,y)Ai(x,y)Ao(x,y) (2.7) where e(E) is the detector efficiency, P, is the incoherent differential scattering cross section for water, p(x,y) is the local fluid density at (x,y), TABLE 2.1 Monte-Carlo multiple scatter calculation results. Stratified Bubbly Annular Inverted Annular Table 2.1 continued. - - multiple scatter single scatter ," Ai and Ao are as defined for equation (2.1), and S is the source strength. There are several implicit assumptions expressed in equations (2.6 ) and (2.7). They are that a) the divergence and width of the source beam in the direction perpendicular to the imaging plane is small, b) the source is a point body, c) the detector is a point body, and d) the effects of the detector electronics and MCA are neglected. Careful construction of the source collimator can insure that a) is true. Including the effects of b), c), and d), poses no serious analytic hurdles as far as the forward problem of predicting the detector output given the density distribution. However, these details have serious implications on the resolution obtainable when trying to reconstruct the density distribution from the detector output. ,Equation (2.4) provides a one-to-one correspondence between the energy of a detected singly scattered photon and the point of scattering. The relationship is true only if the source and detector are point bodies. If the so-uce ana detector have finite extent, then there is some uncertainty associated with the scattering angle as depicted in Figure 2.4. This means that photons scattering at a point in the density field do not arrive at the detector with a single energy but with an energy over a small range, AE = (@E/@e)Ae so e d Figure 2.4: Illustration of the scattering angle blur problem arising from the use of non-point source and detector. This "blurring" effect is analogous to taking a picture with an out-of-focus camera. The amount of uncertainty in the scattering angle is not constant over the entire density field nor even along an isogonic line. Consequently, the blurring cannot be expressed as the convolution of some spread function with the spectrum and cannot De removed by a straignt forward deconvolutional procedure. There is also blurring of the detected photon spectrum by the detector and electronics. Unlike the blurring discussed above, this blurring is generally energy invariant and can be removed by a deconvolutional procedure. Methods for removing the detector system blur and correcting for the energy dependent efficiency are well developed and it is assumed that the measured data are automatically corrected for these effects. Another aspect of the data collection process which affects image quality is the use of a multi-channel analyzer (MCA) to record the photon spectrum. The impact of collec- ting the spectrum data in an MCA is that the spectrum is not available as a continuous function but as a discrete set of values. That is, D(E) becomes the set {Dc;c=l...I} where Dc = D(E)dE. AEc It is now convenient to introduce the isogonic segment as shown in Figure 2.5. An isogonic segment is the region bounded by the isogonic curves which correspond to the Figure 2.5: Illustration of an isogonic segment. (Reproduced by permission. Copyright ISA 1985 [21].) energy limits of a single channel in the MCA. Any features of the spectrum smaller than the width of a single MCA channel are smoothed over and cannot be recovered. Generally this does not pose a problem since modern equip- ment can collect data with very narrow channel widths. Parallel Beam Source/Uncollimated Detector The analysis so far has dealt exclusively with a small source and a diverging beam. It is relevant to extend the analysis to parallel source beams since work has been done using parallel beam illumination with neutron scattering [20]. Figure 2.6 illustrates one possible way of constructing a parallel beam source from an array of narrow beam sources. Alternatively, a parallel beam source can be seen as the result of moving a wide beam point source infinitely far away from the pipe. As the source is moved away, the source-to-detector distance goes to infinity and by equation (2.5) the radius of an isogonic curve also goes to infinity. The curvature (i.e. the reciprocal of the radius) goes to zero and the isogonic curves become lines or rays which originate at the detector. The exit paths of the scattered photons (or neutrons) are, therefore, colinear with the isogonic lines (see Appendix II). The same equations describing the detector response for the wide beam source/uncollimated detector case are valid provided the appropriate limits on the integrals are used. Discrete Forward Equation This section deals with the practical aspects of expressing equation (2.7) in a form which enables computing the integral of equation (2.6*) using a digital computer. The details of how the computations are done are presented in Appendix III. The purpose here is to show the relation- ship between the discrete equations, which are used in reconstructing an image of a flow cross section, and the "exact" integral form. The measured data are available as a discrete set from which the singly scattered flux spectrum must be extracted. The extraction process of removing the multiply scattered flux component and correcting for the detector system efficiency can be expressed as S1 = (D B )/c (2.8) where S1 is the number of photons that arrived at the detector having suffered only a single Compton scattering collision. The notation Sl denotes a vector with components Sc. The object here is to be able to calculate S1 given the density distribution in the illuminated cross section of the pipe. An expression for calculating Sc is derived using equations (2.7) and (2.7*) S = sOf dEf wPAiAodZ (2.9) AEc M(E) Beam Collimators Pipe Wall Monoenergetic Source \ Array Isogonic Line Detector Figure 2.6: A parallel beam illumination/uncollimated detector gamma-ray scattering system. Isogonic Segment ..i ....... ... iiiiiii nI K S.iiiiiiii.. .......... ....... 2 8 99 . 1 2 3 n Figure 2.7: An n x n pixel grid is overlayed the flow cross section in order to facilitate calculation or the singly scattered photon flux. Individual pixels are designated 2by a single number k ranging from 1 to K(=n ). The first step in numerically evaluating equation (2.9) is to overlay the illuminated flow cross section with an n x n pixel grid as shown in Figure 2.7. Equation (2.9) is rewritten as a sum of integrals over the individual pixels S = so I dEf wPkAiAodk (2.10) kc AEc Z(E) where pk(x,y) is equal to p(x,y) for points inside pixel k and is zero elsewhere. The notation kc indicates that the sum is to be carried out only over those pixels in the isogonic segment c. The integrals in equation (2.10) can be expressed equivalently by S = Sol fdxfdyPwPkAiAoPSFc(x,y) (2.11) k The term PSFc(x,y) is the characteristic function of the isogonic segment c defined as PSF (x,y) ={ 1 if (x,y) is inside isogonic segment c c 0 if (x,y) is outside isogonic segment c. Any real imaging system imposes a lower limit on the spatial resolution below which variations in the density can not be seen. If the size of the pixels are about the same as the spatial resolution limit, then it is reasonable (even necessary) to assume a constant density within each pixel. The remaining integral in equation (2.11) is approximated by S = S I pkAckPSFk (2.12) k The factor Ack is an average quantity defined by ffdxdypw(x,y)Ai(x,y)A (x,y) Ack = -k (2.13) -k ffdxdy k where the integrals are carried out only over the area of pixel k. The function PSFck is the fractional area of pixel k which is inside isogonic segment c or PSFck = ffdxdy PSFc(x,y). (2.14) k An implicit assumption in equation (2.12) is that the attenuation factors do not vary significantly across a pixel. This assumption is valid except for pixels near an edge in the density distribution and simultaneously on the edge of an isogonic segment. The assumption also breaks down when the size of the pixels are on the order of a photon mean-free-path. In practice there are no sharp edges to the isogonic segments due to the scattering angle blur. The scattering angle blur causes the isogonic segments to "feather" off at the edges (and overlap) as shown in Figure 2.8. The effect is that the assumption works in real applications where non-point source and detectors are used. The definitions for Ack and PSFck must be modified to account for the scattering angle blur. The new definition for Ack specifies that the attenuation factors must not only be averaged over the area of the pixel but also over the range of possible energies with which the photons can reach the detector. The new definition of PSFck is the fraction of photons reaching the detector from pixel k with energies in the range of channel c multiplied by the pixel area. Equation (2.12) gives an expression for calculating the singly scattered photon spectrum for a detector looking at the scattered photon flux of a known density distribution illuminated by a wide beam mono-energetic gamma-ray source. Each channel of the singly scattered photon spectrum is seen to represent an integral measure of the density distribution within the corresponding isogonic segment. The goal of the remainder of this report is to do the inverse operation of extracting information about an unknown density distribution given measurements of the singly scattered photon flux. Chapter 3 tackles this inverse problem by correlating the size and shape of spectra with the different possible flow regimes. Chapter 4 goes one step further in solving the inverse problem and looks at ways of inverting equation (2.12). PSFch(x,y) Figure 2.8: The scattering angle blur causes the isogonic segments to feather off and overlap. PSFch(x,y) is the fraction of those photons which scatter at point (x,y) toward the detector that have an energy in the spectrum channel ch. CHAPTER 3 METHOD OF SPECTRAL MOMENTS The basic premise of Compton scatter tomography is that the singly scattered photon spectrum contains information about the density distribution. If this is true then each flow pattern should impart some unique characteristics to the spectrum. Examination of the first few moments of the spectrum is a simple pattern recognition technique employed in this chapter to extract these unique characteristics. Description of the Method Consider a Compton scattering setup using a wide beam source and uncollimated detector as illustrated in Figure 3.1. Mono-energetic gamma rays leave the source and undergo Compton scattering interactions inside the pipe. Some of the scattered gammas are received by the detector and the energy spectrum of the detected gamma rays is accu- mulated in an MCA. The total count in each MCA channel is proportional to the density of the material in the isogonic segment for the channel. Note that the source beam is stopped down slightly in order to reduce the amount of the pipe wall which is illuminated. The typical two-phase fluid problem consists of a fluid and its vapor. The vapor is, in general, considerably less Figure 3.1: Compton scatter densitometry configuration. Scattered radiation from the illumination beam is received by the detector. (Reproduced by permission. Copyright o Instrument Society of America (ISA) 1985 [21].) dense than the liquid. Therefore, an MCA channel represents an integral measure of the liquid in the isogonic segment. The singly scattered photon spectrum is a projection of the liquid distribution onto the energy axis. The energy axis can, by suitable transformation, be related to any spatial axis intersecting the isogonic lines. A convenient reference axis perpendicular to the isogonic lines is the Y-axis. The transformation relating the energy axis to the Y-axis is available from equation (2.5) Y(E) = r(9(E)) + yc(e (E)) (3.1) Note that each detector has a characteristic Y-axis. The fluid distribution in the pipe can be inferred from the one dimensional image (i.e. projection) by examining the moments of the spectrum with respect to the Y-axis. The zeroth moment is a measure of the amount of liquid phase present. The first moment is a measure of where the center- of-mass of the liquid phase is. The standard deviation is a measure of how the liquid phase region is distributed about the center-of-mass. The moments are not direct indicators of their respective quantities because of the attenuation factors. The zeroth moment is calculated in a straight forward manner as c where Sc is a channel of the singly scattered photon spectrum [see equation (2.8)]. The first moment and standard deviation are calculated according to the formulae Sc Y YdY 1 C AYc c+l 1C c and S C Y c (Y 2 c AY Yc+l = (3.3b) c sc c Figure 3.2 graphically depicts the zeroth moment, first moment, and standard deviation for annular and inverted annular flow at 50% void fraction. The two peaks in the spectrum for the annular flow correspond to the liquid in the illuminated portions of the annulus separated by the low density vapor core region. The peak heights are not the same because the source beam is attenuated before it reaches the far side of the annulus which results in fewer scat- tering collisions taking place there. Geometric attenuation 1.00 z M rel. units 0 0.75 0 0.50- U T 0.25- 0 10 15 20 Y (cm) (a) Figure 3.2: Spectral moments of the singly scattered photon spectrum. a) Spectrum for an annular flow; b) Spectrum for an inverted annular flow. (Reproduced by permission. Copyright ISA 1985 [21].) Y (cm) (b) Figure 3.2--continued of the photon fluence (i.e. photons/area) for points further from the source also occurs. Similarly, these same attenua- tion terms are applicable with respect to the detector. Only a single peak for the liquid core region is seen in the inverted annular spectrum. The peak is not symmetric due to attenuation as the gamma rays pass through the liquid and to geometric attenuation. The zeroth moments are not identical even though the void fraction is the same. All of the liquid core region is illuminated in the inverted annular flow. Not all of the liquid annulus in the annular flow is illuminated. In addition, the different liquid distribution leads to differ- ences in attenuation, solid angle, and differential scattering cross section. Consequently, the zeroth moment of the inverted annular flow is larger than that of the annular flow. Both the annular and inverted annular flows are symmetrically distributed within the pipe and so, the first moment would be expected to lie at the center of the pipe. However, due to the attenuation effects the first moment is skewed closer to the source and detector. The standard deviation of the annular spectrum is larger than that of the inverted annular. This is expected since the annular spectrum consists of two widely separated peaks while the inverted annular spectrum is a more compact single peak. Procedure and Results Monte Carlo transport calculations are carried out to test whether flow regime identification is possible by examining the first few moments of the singly scattered gamma ray spectrum. Several different flow regimes are analyzed; in each case the energy spectrum from single scattering interactions with liquid or vapor inside a pipe is calculated for two detectors. The detectors are assumed to be placed symmetrically with respect to the pipe and source as shown in Figure 3.3. The source-pipe-detector geometry is the same for all cases. The source and detector are treated as points for purposes of calculating the scattering angle. However, for purposes of calculating the solid angle view of each interaction, the detector has a radius of 0.5 cm. Detected scattering interactions within the pipe wall are not considered since the primary focus is to evaluate the gamma rays coming from the fluid. However, attenuation due to the pipe is taken into account. The flow regimes considered are annular, inverted annular, stratified, bubbly, and mist. The annular and inverted annular flows are modeled as a liquid/vapor annulus and a vapor/liquid central core concentric with the center of the pipe, respectively. The stratified flow is modeled as a pipe half-filled with liquid as illustrated in Figure 3.4. The bubbly and mist flows are simulated by choosing fifteen randomly sized (0.85 bubbles of vapor/liquid in a matrix of liquid/vapor, respec- tively. New bubbles are periodically chosen during a cal- culation to simulate stochastic flow conditions. The void fraction is fixed at 50% for the annular, inverted annular, and stratified flow. Due to the random nature of choosing the bubbles, the bubbly flow has a void fraction of ~45% and the mist flow has a void fraction of ~58%. Four calcula- tions are made for each flow regime where the width of the source beam in the plane of the pipe cross section varies from 4.5 to 28 degrees. The results of the moment calculations are tabulated in Table 3.1. Interpreting the Moments The moments are functions of the magnitude, shape, and width of the spectra. The spectra change as the source beam width varies since the amount of the flow cross section illuminated changes. Stating general trends of the spectra based on the changing amount of illuminated fluid is diffi- cult due to the effects of attenuation. However, the following two observations can be made. First, the magni- tude of the spectrum from a narrow source beam is smaller than that from a wider source beam. Assuming the source strength remains constant, the number of photons illumina- ting the pipe is less for the narrow source beam. The result is that the zeroth moment becomes smaller as the source beam narrows. Second, the width of the spectrum decreases as the source width is narrowed. As illustrated in Figure 3.5, the Detector 1 Source '_, ipe Detector 2 Figure 3.3: One source, two symmetrically positioned detector configuration used for the Monte Carlo studies presented in this chapter. (Reproduced by permission. Copyright ISA 1985 [21].) Detector 1 Source Vapor Liquid Figure 3.4: 4 Detector 2 Source-detector configuration relative to the stratified flow regime case. (Reproduced by permission. Copyright ISA 1985 [21].) TABLE 3.1 Zeroth moment, first moment, and standard deviation of the singly scattered photon spectra for the five flow regimes listed. Source Beam Width Flow Regime ANNULAR Detector 1 Detector 2 INV. ANNULAR Detector 1 Detector 2 STRATIFIED Detector 1 Detector 2 BUBBLY Detector 1 Detector 2 MIST Detector 1 Detector 2 [relative] 280 [cm] 0.311 0.313 0.337 0.340 0.621 1.00 0.434 0.436 0.431 0.423 <02> [cm2] 140 0.0935 15 0.0935 15 3.7 0.194 16 3.7 0.193 16 6.1 0.360 17 5.4 0.353 15 8.4 0.204 16 8.4 0.205 16 6.6 0.183 16 6.9 0.184 16 80 15 15 3.7 3.6 0.0485 15 0.0484 15 0.108 16 0.108 16 4.50 3.5 3.6 6.6 0.207 16 6.9 6.3 0.158 16 6.6 7.6 0.108 16 8.1 7.5 0.107 16 8.0 6.5 0.108 16 5.9 6.7 0.112 16 5.7 0.0267 15 0.0267 15 0.0601 16 0.0600 16 0.116 16 0.0789 16 0.0597 16 0.0602 16 3.5 3.5 6.7 6.4 8.2 8.2 0.0614 16 5.3 0.0600 16 5.2 Note: The spectra are from a Cs-137 source illuminating a 12 cm diameter pipe with a 0.75 cm thick steel wall. The beam divergence normal to the flow cross section is 5 degrees in all cases. The void fraction is nominally 50% for all cases. region within the isogonic segments associated with the extreme channels of the spectrum is no longer illuminated as the beam narrows. The effect is more pronounced for the isogonic segments at low values of the Y-coordinate. If liquid is present in the affected regions, the result is that the spectrum narrows and the narrowing can be more pronounced at the low end of the Y-axis. Such an asymmetric change in the spectrum width means the first moment shifts to larger values of Y as the beam narrows. A narrower spectrum would also be expected to yield a smaller standard deviation. The magnitude of change in the first moment or standard deviation, or whether any change occurs at all, depends on any accompanying change in the shape of the spectrum. If there is no liquid in the affected regions then, of course, there is no change in the width of the spectrum. A Detailed Example The following is a discussion in physical terms of how the moments vary with changes in the source beam width. So that sufficient attention can be given to the changes in each moment, the discussion is confined to only the strati- fied flow regime. Enough insight should be gained to interpret the data in Table 3.1 for the remaining flow regimes. Figures 3.6 and 3.7 are the spectra for the upper and lower detector, respectively, viewing a stratified flow. These figures show the progressive changes in the magnitude, shape, and width of the singly scattered photon spectrum as the source beam width narrows. The magnitudes of the spectra decrease as the beam narrows because there are fewer photons illuminating the flow. The magnitudes of the two figures should not be compared since the vertical scales are diffe- rent. Instead the figures should be compared on the basis of the zeroth moments. The zeroth moment of the lower detector's spectrum is larger when the source beam is wide, while the zeroth moment of the upper detector's spectrum is larger when the source beam is narrow. The behavior of the zeroth moments can be explained by considering the effect of attenuation on the scattered gamma rays and the differences in solid angle each detector views of the scattering interactions. When the source beam is widest the average attenuation of the scattered gamma rays reaching each detector is about the same. The average position of scattering occurs in the middle of the liquid and the scattered gamma rays must pass through almost the same amount of liquid and steel to reach either detector. When the source beam is narrowed the average scattering interaction takes place toward the center of the pipe and the scattered gamma rays must pass through more liquid to reach the lower detector. Thus, attenuation of the scattered gamma rays explains why the zeroth moment of the upper detector is larger when the source beam is narrow. When the source beam is wide the zeroth moments should be the same since the average attenuation is the Figure 3.5: Darkened regions are those isogonic segments which are eclipsed when the more restricted beam illumination is used. (Reproduced by permission. Copyright O ISA 1985 [21].) Z beam line zeroth 3 width pattern moment 0 (degrees) (relative) .75 28 0.621 0 22 ---- 0.532 '- 14 --- 0.360 S8 -- 0.207 \ 1- 5 4.5 -- 0.116 LU S.25- h- U.1 -- - scattered gamma rays for upper detector viewing a stratified flow as in Figure 3.4. (Reproduced by permission. Copyright ISA 1985 [21].) 1.0 I-- beam line Izroth 2 width pattern moment (degrees) (relative) 0 0 .7 28 1.000 S22 ---- 0.708 o 14 --- 0.353 A -.. 0.158 0 4.5 ---- 0.079 F- .5- \ . LU \ p.25- 0 5 10 15 20 Y (cm) Figure 3.7: Detected photon energy spectra of once- scattered gamma rays for the lower detector viewing a stratified flow as in Figure 3.4. (Reproduced by permission. Copyright ISA 1985 [21].) same. The reason they are not is because the upper detector is further away and presents a smaller solid angle to each scattering interaction. The behavior of the first moments can be explained in the following manner. Observe that the spectrum for the widest source beam is peaked at about 13 cm for the lower detector and at about 17 cm for the upper detector. Due to the source-pipe-detector geometry the isogonic segments do not all intersect the liquid region. The projection of the liquid region onto the corresponding Y-axis for the two detectors is not symmetric. As the source beam narrows both spectra become more uniform since approximately the same amount of liquid is illuminated in each isogonic segment. The remaining slight peak is due to the non-uniform amount of attenuation and to solid angle considerations. The more uniform a spectrum the more the center-of-mass coincides with the center of the spectrum. The result is that the first moment shifts to higher Y values for the lower detector and shifts to lower Y values for the upper detector. Note also the preferential narrowing of the spectra which tends to enhance the changes in the first moments. The standard deviation for both detectors' spectra increases uniformly as the source beam narrows from 28 to 8 degrees and then decrease slightly as the beam narrows further to 4.5 degrees. When the source beam is widest both spectra are peaked near the center-of-mass which gives a relatively small standard deviation. The standard deviation increases as the source beam narrows because the spectra become more uniform. The slight narrowing of the spectra leads to a smaller standard deviation. The two competing effects result in the observed behavior. Using the Moments The purpose of calculating the moments is to identify the flow regime automatically and without any a priori knowledge. Inter-comparing the moments of all the flow regimes at each particular source beam width suggests that using the widest source beam offers the greatest discrimina- tion between flow regimes. A close agreement between the moments of the upper and lower detectors is seen in all the data for the annular, inverted annular, bubbly and mist flows where the fluid distribution is symmetric. The stratified flow is asymmetric, but only for the wider source beams (i.e. > 14 degrees) is the asymmetry discernible from comparing the moments of the upper and lower detector. The annular flow is discernible from the other symmetric flows with the wide source illumination because of its large standard deviation. The inverted annular flow conversely is distinguishable because of its small standard deviation. Discriminating between the bubbly and mist flows is not possible on the basis of the first moment or standard deviation due to their similar nature. Bubbly flow is expected when the void fraction is low (i.e. < 50%) while mist is expected when the void fraction is high. The zeroth moment is sensitive to the void fraction and could be used to determine which flow exists. The results presented so far show that flow regime identification can be made using the singly scattered photon spectra from two detectors sensing the emerging scattered photon flux from a flow cross section illuminated by a wide beam source. A tomographic image of the flow can be produced using a simple model of the flow cross section for each flow regime. The consistency between the actual flow regime and the flow regime as identified by the spectral moments can be checked using the solution to the forward problem [i.e. equation (2.12)]. The calculated singly scattered photon spectra for the two detectors based on the identified flow regime should match the measured spectra. If the calculated and measured spectra are not close enough, then the flow cross section model is adjusted accordingly. For example, if the spectral moments indicate that the flow is annular, then the size of the annulus is adjusted until the measured and calculated spectra match. The next chapter develops other methods of producing such tomographic images of the flow. These other methods are more systematic in adjusting the flow cross section image so that the measured and calculated singly scattered photon spectra match. CHAPTER 4 THE RECONSTRUCTION PROBLEM A general statement of the problem at hand is Given a set of singly scattered photon flux measure- ments, {S }, made using a wide-beam source and uncollimated detectors, reconstruct the liquid-vapor phase distribution {p} that produced {S }. Ideally the set of flux measurements is acquired with a single exposure. That is, the source and detector locations are fixed and there are no moving parts. Insight into this problem can be gained by considering the following formulation of the forward problem from equations (2.6*) and (2.7): SI=[T]p (4.1) where [T] is the matrix which mathematically describes the scattering transport process, S1 is the vector of singly scattered photon flux measurements, and p is the vector of pixel densities. The obvious solution of inverting [T], that is, p=[T]-1 S1 (4.2) is not strictly possible since the elements of [T] are functions of p. The problem at hand is nonlinear. Until now the general approach has been to avoid dealing directly with the nonlinearity of [T]. One often-proposed method (e.g. [24, 25]) is to use multiple source-detector locations and/or source energies in order to cancel out the nonlinear aspects of [T]. There are numerous practical limitations in implementing this approach for dynamic measurements. Another proposed method is to assume a density distribution based on prior experience and use this in computing [T]. Since [T] is a function of the density, a bias in p is likely for the high contrast reconstructions expected from two-phase flows. A charac- teristic of these methods is the idea of calculating the density distribution in an algebraic, once-through fashion. It is proposed here that the problem be viewed as a tomographic reconstruction problem where the measured data, S, represent projections of the phase distribution, p, onto non-physical energy axes. The projection data are incomplete since only a few detectors are employed. Simila- rly the elements of [T] are at best known only approximately since they depend on the unknown phase distribution p. When the available information is incomplete and/or imperfect, an algebraic solution may not exist or several solutions may exist. These are familiar situations in CT and the analytic techniques are already developed to handle them. The analo- gous projection operator, [T], of current CT applications is mathematically well-behaved (i.e. linear, single-valued, etc.). The CT reconstruction methods must, therefore, be modified to handle the nonlinearities of Compton scatter tomography (CST). Several CT reconstruction techniques are adapted to CST in this chapter. General Solution Methods The basic methods of CT for solving the reconstruction problem can be grouped into five categories: 1. Fourier transform: exponential Radon transform [26], filtered backprojection [27]. 2. Direct scanning [12, 14, 15]. 3. Matrix inversion: successive approximation [20], maximum entropy (see Appendix IV) [28], generalized inverse, and pseudo inverse [29]. 4. Series expansion [30]: algebraic reconstruction technique (ART), simultaneous iterative reconstruc- tion technique (SIRT), and iterative least squares (ILS). 5. Statistical: Monte-Carlo backprojection. The transform techniques are recognized as being the more elegant, mathematical reconstruction methods. Trans- form methods in general require a large number of projec- tions and do not show good noise immunity. The cost of detectors and the physical limitations in mounting and cooling the detectors imply a limit on the number of projec- tions available. Practical limits on source size, limits on detector size, and desired framing rate are all going to result in statistical noise in the measurements. These features of using Compton scatter to image dynamic two-phase systems make transform methods unsuitable as a recon- struction technique. Direct scanning techniques inherently require a narrow collimated source beam and a moving source detector arrange- ment. A narrow collimated source is less photon efficient than a wide collimated source. This inefficiency coupled with the time necessary to translate the source and detec- tors means that a direct scanning system is slower than a comparable wide-beam source, uncollimated detector system. Since one goal is to measure dynamic two-phase flow, the faster measurement system is the obvious choice. There is another reason for dismissing direct scanning techniques that involves accounting for the multiply scattered photon flux. The structure of the flow regime is more strongly imprinted on the multiply scattered photon flux when a narrow-beam source is used than when a wide-beam source is used. Any structure in the multiply scattered photon spectrum means that approximating the multiply scattered photon spectrum by a general function B(E) [see equation (2.6*)] is not possible. Removing the multiple scatter component in direct scanning systems, therefore, requires more computation and is more prone to error. Matrix inversion and series expansion reconstruction techniques are equivalent in that a solution is sought by iteratively improving an initial guess of p (and hence [T]) until some integral transformation of the reconstructed phase distribution and the true phase distribution are judged close enough with respect to the measurement error. The difference between matrix methods and series expansion methods is that at each iterative step matrix methods update all the pixel values at once, whereas series expansion methods update the pixel values sequentially. Matrix methods usually exhibit better convergence rates than series expansion methods since they can take into account higher order effects in a more straightforward manner. Techniques from either of these categories appear to be well suited for use in a CST system. A Closer Look This section details how specific CT techniques are adapted for use in CST. The first technique to be considered is backprojection. Generally the initial step in a reconstruction algorithm is a backprojection operation which yields a very smooth representation of the true image. Much of the fundamental differences between traditional CT and CST are evident in the way the projection data are handled in order to perform a basic backprojection opera- tion. The second technique considered is Monte-Carlo backprojection (an original concept by this author and Dr. E. E. Carroll at University of Florida). The result of Monte-Carlo backprojection is a higher contrast image than obtained with traditional backprojection. The third part of this section is a general discussion of matrix inversion techniques. The purpose of this discussion is to introduce the basic iterative nature of any practical solution strategy imposed by the nonlinearity of CST. The fourth part of this section adapts three series expansion techniques (ART, SIRT, and ILS) for use in CST. Backprojection Consider a projection formed by the conventional CT scanner shown in Figure 4.1. The pencil beam source irra- diates the object and the transmitted radiation is measured at points in the image plane. The source and detector are moved in the direction indicated to scan the entire object and form a single projection. Other projections are formed by rotating the source and detector around the object and taking another scan. The number of photons reaching the detector at each point in a projection is a function of the average optical path through the object. The functional relationship is Ijc = Ioexp{-f p(s)pt(s)ds} (4.3) kjc where Ijc is the output of the detector at position c in projection j, I is the number of photons leaving the source, p(s) is the density at point s, pt(s) is the total mass attenuation coefficient of the material at s, and Pjc indicates that the integral is to be carried out Pencil Beam Source -- ----Image Plane Detector Figure 4.1: Example of a conventional transmission CT scanner. The source and detector move parallel in the direction indicated by the straight arrows, and also rotate about the object being imaged. along the line between the source and detector when at position c in projection j. In practice one does not work with Ijc, but instead with Rjc given by K Rjc = -n(Ijc/Io) = t Pksjck (4.4) k=l where pk is the density in pixel k, Asjck is the length of line jc contained in pixel k, and K is the total number of pixels. The value of Asjck is nonzero only if line jc crosses pixel k and so the sum in equation (4.4) is equivalent to the integral in equation (4.3). Since the liquid and vapor are of the same material, the mass attenuation coefficient is a constant and is taken outside the integral. The average density, Pjc, along line Zjc can be computed from equation (4.4) as K Pjc = Rjc/{t Asjck} = Rjc/Ljc (4.5) k=l A backprojection image is made by forming a composite of the average line densities, Pjc' J C(j) bpk ={ I PjcASjck/ j=l c=l J C(j) 1 sjck j=l c=l (4.6) where bpk is the density of pixel k in the reconstruction, and Asjck is used as a weighting factor. An analogous relationship for computing a backpro- jection image in CST can also be developed. The quantity corresponding to Rjc in CST is the singly scattered photon flux values, Sc, given by K S= S PkAjkPSFjckAs k=l where p is the incoherent scattering coefficient, S is the uninteracted source flux, Pk is the density in pixel k, As is the length of one side of a (square) pixel, Ajk is an average quantity which takes into account all aspects of photon transport from the source to pixel k, and from pixel k to detector j, and PSFjck is the.fraction of the photons arriving at detector j from pixel k which have an energy in the range of channel c (Ec to Ec+l). Note that PSFjck is zero if it is not physically possible for a singly scattered photon from pixel k to arrive at detector j with an energy in the range Ec to Ec+1. So, the sum of equation (4.8) represents an integral over the isogonic segment c of detector or projection j and is equivalent to the sum of equation (4.5). Another point to (4.7) note about equation (4.7) is the that double subscript jc can be replaced by an equivalent single subscript g, K S = o PkAgkPSFgkAS (4.7*) k=l which is convenient in formulating the problem in matrix notation as in equation (4.1). The average density in an isogonic segment, Pjc, is computed as K jc = Sc /AsS I AjkPSFjck (4.8) k=l and similarly the backprojection image densities, bpk, for CST are computed as J C(j) I I PjcPSFjck bpk = j=l c=l .(4.9) J C(j) 1 l PSFjck j=l c=l A basic problem with equation (4.8) is that the density distribution must be known in order to calculate the Ajk factors. A reasonable solution is to use a homogeneous distribution with an average void fraction that is either assumed or gained from an independent measurement. The average void fraction of a CST backprojection image must then be checked for consistency with the initial void frac- tion. Any inconsistency can be resolved by some iterative refinement. Equations (4.5) and (4.6), in contrast, require no knowledge of the average void fraction and a one-step backprojection computation is possible in CT. Monte-Carlo Backprojection The motivation for using a statistical approach to form a backprojection image comes from three observations. First, the singly scattered photon flux, Sjc, is directly proportional to the average density in the isogonic segment jc. Second, any detected photons are due almost entirely to scattering in the liquid phase since the density ratio of liquid-to-vapor phase for water is usually 1000 to 1 in nuclear reactors. Third, the features of one detected spectrum can be cross correlated to the features of another spectrum detected at another location. These observations are summarized by the following example. If a single peak is observed in channel c of detector j and a similar peak is observed in channel c" of detector j', then the liquid phase is localized to the region where isogonic segment jc intersects segment j'c'. Monte-Carlo backprojection treats each detected spectrum as an independent one dimensional probability density function of where the liquid phase is located. A reconstruction of the liquid phase distribution is made by statistically cross correlating the detector spectra. The difference between Monte-Carlo backprojection and traditional backprojection is that traditional backprojec- tion ignores the cross correlation information. The initial step in Monte-Carlo backprojection is to form an independent probability density function (pdf) from the modified average isogonic segment densities of each detector PiC pdf (c) =. (4.10) C(j) SPjc c=l The average isogonic segment density is modified in the following manner: K P;c = S/ Asso X Ajk (4.8*) k=l in order to account for the non-uniform spatial width of the isogonic segments. Assume that a fine-mesh square pixel grid is employed and that a pixel is either liquid filled or vapor filled. Provided a fine enough mesh is used, most density distribu- tions can be reasonably represented by such a scheme. The number of pixels, N, which are liquid filled is easily computed from an estimate of the average void fraction, a, N = round((l ()-K) (4.11) where K is the total number of pixels. All the pixels are assumed vapor filled to start and the determination of which N pixels to designate as fluid filled proceeds as follows: 1. Given a set of J random numbers {r-}, an associated set of channel numbers {chj} is generated by inverting the pdf 's ch rj 1 pdfj(c); j=l, 2, ..., J (4.12) c=l There are J detectors and one channel from each detected spectrum is chosen by equation (4.12). 2. Associated with each channel chosen is a set of pixels, Set(ch.), which contains all the pixels within the corresponding isogonic segment. A new set is formed by taking the intersection of all J of these Set(chg)'s. The intersection set, I-Set, contains a t6tal of M pixels. Taking a new random number, rJ+1, the mth element m = integer(M-rj+1) + 1 (4.14) of I-Set is chosen as the pixel to be fluid filled. If M is zero, then no pixel is designated as being fluid filled. 3. Steps 1 and 2 are repeated sequentially until the necessary N pixels are turned on. Figure 4.2 compares the results of Monte-Carlo and traditional backprojection. The Monte-Carlo reconstruction is not as blurred and better represents the actual flow regime. Although the reconstructed image resembles the correct flow regime, a comparison of the calculated and detected singly scattered photon flux spectra shows that there is still room for improvement (Figure 4.3). (a) (b) (c) META ULFd .Detector Source * -.Detector Figure 4.2: Reconstruction of a 50% void fraction annular flow regime. a) Monte-Carlo backprojection; b) Simple backprojection; c) source and detector locations. Calculated ------- "Measured" (a) (b) Figure 4.3: Comparison of the calculated singly scattered photon flux for the Monte-Carlo backprojection reconstruction of Figure 4.2 to the "measured" spectra. a) Upper detector; b) Lower detector. The addition of a feedback mechanism to the pixel selection process is one possible improvement. Instead of randomly choosing a pixel from the intersection set, those pixels which improve the calculated spectra could be chosen preferentially. A more complete evaluation of the Monte- Carlo reconstruction method and possible variations is left for future work. Matrix Methods Generally matrix inversion methods are avoided in CT due to the large size (i.e. > 1000 x 1000) of the matrices involved. Even in CST where there are only a few projections and, with a coarser pixel mesh, the matrices can be on the order of 200 x 200. Inverting matrices of this size requires large amounts of storage, long computing times, and is prone to round off error problems. Whatever the difficulties in implementing a matrix method solution, the compact notation is useful. There are many ways of formulating the reconstruction problem as a series of matrix operations. The most basic approach is to treat the problem as an optimization problem. A search is made for the density distribution, p*, which minimizes, either directly or indirectly, the difference between the measured data, S1, and the calculated singly scattered photon flux, F, that is, t All data in this report, unless otherwise stated, are the result of Monte-Carlo photon transport calculations (i.e. numerical experiments) and not physical measurements. min IS1 FI or min 1S1 [T(p)]p| (4.15a) P P subject to the constraints Pvapor < Pk < Pfluid k=l,2,...,K. (4.15b) The notation I..I represents some norm functional. Any reconstruction algorithm in CST has a structure similar to that depicted in Figure 4.4. Note that the nonlinearity of the problem requires the recalculation of the transport matrix for each new density distribution, pm. Changes in the density distribution directly affect the attenuation of the source photons and scattered photons which is accounted for in the transport matrix. If the variation in pm is small so that [T] is approximately constant, then it may be possible to skip the recalculation step for a few iterations. In contrast, the corresponding matrix in CT is only a function of the system geometry and remains constant throughout the calculation. Another point at which the nonlinearity of CST needs to be addressed is in the calculation of Em+l. One method of calculating Tpm+l is to use the method of steepest descent pm+l = (Jacobian of P(pm))-1.(l p( m)) (4.16) The transport matrix can be treated as a constant or as a function of the density in the calculation of the Jacobian. Figure 4.4: Flow chart outlining the general structure of a reconstruction algorithm for Compton scatter tomography. The superscript m denotes the iteration index. Treating [T] as independent of p means a smaller step size, f, is necessary and that Tp does not point directly toward the path of steepest descent. Including the dependence of [T] on p results in a significant increase in the storage requirements and in computing time. Algebraic Reconstruction Technique ART is historically the first of the series expansion techniques applied to commercial CT devices [31]. There are two basic forms of ART algorithms: multiplicative and addative. Presented here is a variant of the additive algorithm. A single ART iterative step involves correcting the density vector by adding to it a correction vector 5m+l = pm + cm+l.wm+l (4.17) where cm+l is a constant correction factor and wm+l is a weighting vector. The correction factor is computed so that one of the calculated singly scattered photon flux values is consistent with the corresponding measured value, that is, S Fm+l (= Tk p ) (4.18) Subsequent iterations correct the density vector with respect to the other measured values, one at a time. The subscript g is incremented cyclically gm+l = (gm + 1) modula G + 1 with each iteration in order to cover all G data points. The correction factor can be derived by substituting equation (4.17) into (4.18) S1 VTm+Ipm cm+l 9 g k gm+l. (4.19) yTm+lwm+l gk The transport matrix elements are unknown until the new density vector is known which in turn requires the correction factor. One solution to this predicament is to use the transport matrix calculated during the previous iteration, that is, ql TM Mpi cm+l S gk P gm+l (4.19*) ITm m+l TgkWk This approach is valid whenever the transport matrix does not change much from one iteration to the next, as when close to the final solution. There are several possible choices of weighting vectors. The simplest choice is an identity vector, however, this results in all the pixel values being altered including pixels which do not affect the sum in equation (4.18). An identity vector modified such that m+l 1 if Tg 0m+ w, =0 gk f m+l k { 0 if Tgk 0 g = g (4.20a) sensibly corrects only those pixels contained in the isogonic segment corresponding to S1. Another reasonable weighting vector is w = Tk g = gm+1 (4.20b) which effectively weights a pixel according to the proba- bility that a scattered photon from that pixel is counted in S A third possible weighting vector is wm+1 = PSFgk; g=gm+l (4.20c) which is analogous to the fractional ray length weighting scheme of CT [32]. The weighting scheme chosen here and elsewhere in this report is that of equation (4.20b). Figure 4.5 graphically shows how an ART solution proceeds for a two data point, two pixel situation. The lines Lm and Lm correspond to the two simultaneous equations given by Lm S1 = Tml + Tm22 ; g = 1,2. (4.21) The lines change as the density vector is corrected since the transport matrix is a function of the current density distribution iterate. The lines tend to stabilize as the 1.0 L L= L L, 1 2 2 3 p2 04Li L- - L \ 0.8 L 3 22 L1 0.6 0 - P 1 0.4 - pp 0.2 0 0.2 0.4 0.6 0.8 1.0 P, Figure 4.5: Graphic demonstration of how an ART solution proceeds for a two data point, {Li}, two pixel, {pi}, solution. Superscripts identify the iteration number. final solution is approached, indicating that the transport matrix is only varying slightly. Simultaneous Iterative Reconstruction Technique SIRT is similar to ART except that the density vector is corrected using all the projection data, 91, simulta- neously. SIRT attempts to combine a complete cycle of ART iterations (i.e. G iterations) in calculating the correction vector. The correction vector elements are given by I (Tmk/a ){(S Fm)/( ~Tk p+ = g k (4.22) I (Tmk/og) g where a2 is the measurement variance. The correction vector g resembles a backprojection image of the error vector m = S1 F" (4.23) which can be seen by comparing equations (4.22) and (4.9). Reference 30 indicates that some method of damping may be necessary to avoid a solution which oscillates and fails to converge. Once the correction vector is obtained it is multiplied by a constant damping factor, f. The damping factor is calculated so that the variance weighted euclidean norm of the error vector is minimized min S1 [Tm](pm + fm+l1m+l) (4.24) f The solution to equation (4.24) is +1 (em/2 )(Tm Am+1 H {m+l g/g)(gkk)} f = _g k (4.25) I{(Tm m+l }2 { k TkAPk )/Gg} g k The proposed SIRT algorithm for CST is Pm+l = mid {Pvapor' Pk + fm+lApm+' Pfluid} (4.26) where the correction vector is given by equation (4.22) and the damping factor by equation (4.25). Iterative Least-Squares Technique The ILS technique computes a correction vector which minimizes the variance weighted euclidean norm of the error vector. This can be expressed as U = I{(Sl Fm TkAp+)/a }. (4.27) g K The correction vector which minimizes U is found by setting the derivative of U with respect to the correction vector equal to zero. This leads to the set of equations pm+l(TmkT /a) = iTm-em/aG; k'=l,...,K (4.28) k g g Solving the system of equations of (4.28) would violate the principal motivation of going to an iterative approach, which is introduced to avoid the computational overhead and pitfalls of matrix methods. The usual course of action is to assume that all the the correction vector elements are zero except the element corresponding to the current pixel k'. Equation (4.28) becomes Vrm em/a2 Apk+1 = g k k'=l,...,K (4.28*) (ITmk/og)2 g Equation (4.28*) overestimates the magnitude of the indivi- dual elements of the correction vector which necessitates the introduction of a damping factor as is done for SIRT. The form of the ILS damping factor, which is calculated according to equation (4.24) also, is 1{(em/a2)(Tm -pm+l m+l {(eg/ g)(TgkPk)} fm+l = -g k g. (4.29) Tm Apm+l)/7 }2 g k The proposed ILS algorithm for CST is p+1 = mid vaporo' m + fm+1lm+ Pfluid} (4.30) where the correction vector is given by equation (4.28*) and the damping factor by equation (4.29). Comparison of Reconstruction Techniques The ART, SIRT, and ILS reconstruction techniques are compared here in order to select one of them to use in a detailed evaluation of CST. Each of the techniques are used to reconstruct the two flow regimes shown in Figure 4.6. The flow regimes are analogous to the "hot spot" and "cold spot" phantoms used in CT studies. The criterion used to judge the best technique is the euclidean distance between the true density vector and the reconstructed density vector. The rate of convergence is considered as a secondary criterion. The projection data (spectra) used for each technique are generated by Monte-Carlo transport calculations. The singly scattered photon spectra calculated for the two detectors are shown in Figures 4.7 and 4.8. Table 4.1 lists the pertinent "experimental" conditions. The reconstructions of the hot-spot and cold-spot flow regimes are shown in Figures 4.9 and 4.10, respectively. All the reconstructions exhibit undesirable streaking artifacts. The artifacts are due to the lack of edge infor- mation from using only two views. The reconstructions of the hot-spot are better than those of the cold spot. Intui- tively this is expected by observing that the signal-to- noise ratio is higher for the hot-spot spectra. The ART results have more pronounced streaking artifacts than either the ILS or SIRT results. The poorer performance of ART is in keeping with similar CT phantom study tests done for Detector Source Collimator (dimensions in centimeters) K5 *- 17 of -- 17 Detector -.j x =6 15 .j Figure 4.6: Source and detector configuration used for the hot-spot and cold spot tests. Energy [keV] (a) 400 400 Energy [keV] (b) Figure 4.7: The measured singly scattered photon spectra for the hot-spot test. a) Spectrum from the upper detector; b) Spectrum from the lower detector. \ iI' 4J _P r i1 U) o 0 U U 0 I 0 0 0 400 0 400 Energy [keV] Energy [keV] (a) (b) Figure 4.8: The measured singly scattered photon spectra for the cold-spot test. The ordinate scale is the same as for Figure 4.7. a) Spectrum from the upper detector; b) Spectrum from the lower detector. TABLE 4.1 Data pertaining to the hot-spot and cold-spot experiment. Detectors: Number of detectors 2 Detector type planar Diameter 1 cm Active thickness 0.1 cm MCA channel width 1 keV Source: Isotope Cs-137 Physical shape cylinder Diameter 0.31 cm Height 0.63 cm Material Density: Hot-spot bubble density bulk density Cold-spot bubble density bulk density 1.155 g/cc 0.0 g/cc 0.0 g/cc 1.155 g/cc Y -x -.1 * 1. ILS t.2 SIRT ART Figure 4.9: C";: :'1 4 "- 7 Reconstructions of the hot-spot flow pattern by the ILS, SIRT, and ART algorithms. Y x ILS SIRT ART Figure 4.10: Reconstructions of the cold-spot flow pattern by the ILS, SIRT, and ART algorithms. medical applications. There is no discernible difference between the ILS or SIRT results. The rate of convergence and accuracy of the reconstruc- tion methods is shown in Figure 4.11. The discrepancy function is defined as the euclidean distance between the true flow regime and the reconstructed flow regime. The spectra difference function is defined as the euclidean distance between the measured spectra and the spectra calcu- lated using the reconstructed flow regime. Both the discre- pancy and spectra difference functions are normalized by their values for the initial density distribution, which is an all vapor condition. The convergence of all the recon- struction algorithms for the cold spot test is similar. The ILS and SIRT methods show better convergence properties than the ART method for the hot spot test as judged by the discrepancy function. The ILS and SIRT method are more desirable than the ART method due to their consistent performance. There is no significant difference in the performance of the ILS or SIRT methods based on examining the discrepancy or spectra difference functions. The choice of whether to use either the ILS method or the SIRT method appears to be one of personal preference. This chapter successfully demonstrates that CT recon- struction algorithms can be used with CST. As vital as the reconstruction algorithms are, they are not sufficient for implementing a CST system. The next chapter considers the effect of noise in the spectra on the reconstruction process Iteration * ILS o SIRT A ART 5 Iteration (a) Figure 4.11: Comparison of the accuracy and convergence of the ILS, SIRT, and ART algorithms, a) Hot-spot flow regime; b) Cold-spot flow regime. cu a) 0 0) 0- 0) 4- 4-1 **- Iteration OILS o SIRT A ART 5 10 Iteration Figure 4.11--continued. >1 u ) a4 0 Ul -H 0 |

Full Text |

TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT USING COMPTON SCATTERING OF GAMMA RAYS By DAVID E. BODETTE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1986 Copyright 1986 by David E. Bodette ACKNOWLEDGMENTS I thank Dr. Alan M. Jacobs, my adviser, for his patient guidance and careful attention to detail. This work was partially supported by a Liquid Metal Fast Breeder Reactor Technology Research Award from the Argonne National Laboratory under contract with the United States Department of Energy (DOE). The author was also supported by the DOE Energy Graduate Traineeship Program. iii TABLE OF CONTENTS Page ACKNOWLEDGMENTS iii LIST OF TABLES vi LIST OF FIGURES vii ABSTRACT xi CHAPTERS 1 INTRODUCTION 1 2 THE FORWARD PROBLEM 10 Collimated Point Source/Collimated Detector . . 11 Collimated Point Source/Uncollimated Detector 15 Uncollimated Point Source/Uncollimated Detector 17 Parallel Beam Source/Uncollimated Detector. . . 28 Discrete Forward Equations 29 3 METHOD OF SPECTRAL MOMENTS 3 5 Description of Method 35 Procedure and Results 42 Interpreting the Moments 44 A Detailed Example 46 Using the Moments 52 4 THE RECONSTRUCTION PROBLEM 54 General Solution Methods 56 A Closer Look 58 Backprojection 59 Monte-Carlo Backprojection 64 Matrix Methods 69 Algebraic Reconstruction Technique 72 Simultaneous Iterative Reconstruction Technique .... 76 Iterative Least-Squares Technique 77 Comparison of Reconstruction Techniques .... 79 4 COMPTON SCATTER TOMOGRAPHY 90 An Example 90 Accuracy 92 Resolution 106 Image Contrast 109 IV Page 6 SUMMARY AND CONCLUSIONS 118 Summary 119 A Practical Application 121 Conclusions and Future Research 126 APPENDICES I MONTE-CARLO TRANSPORT PROGRAM 130 II INSIGHTS AND FINER POINTS 133 Parallel Versus Diverging Source Beams 133 Multiple Scatter Versus Detector Position . . . 137 Source Placement 142 Convergence Criterion 143 III RECONSTRUCTION COMPUTER PROGRAMS 144 Running the Programs 144 Flow of Programs 145 IV MAXIMUM ENTROPY METHOD 146 V DATA 149 REFERENCES 191 BIBLIOGRAPHY 195 BIOGRAPHICAL SKETCH 196 v LIST OF TABLES Page 2.1 Multiply scattered photon spectra 22 3.1 Spectral moment results 45 4.1 Data for hot-spot and cold-spot tests 83 5.1 Test parameters and values 98 5.2 Summary of reconstruction results 105 5.3 Parameters affecting image resolution 107 vi LIST OF FIGURES Page 2.1 Collimated source/collimated detector configuration 12 2.2 Collimated source/uncollimated detector configuration 16 2.3 Uncollimated source/uncollimated detector configuration 19 2.4 Illustration of scattering angle blur 25 2.5 Illustration of an isogonic segment 27 2.6 Parallel beam source/uncollimated detector configuration 30 2.7 Overlaying of pipe cross section with an n x n pixel grid 30 2.8 Feathering of isogonic segments 34 3.1 Compton scatter densitometry configuration 36 3.2 Spectral moments of a single-scattered spectrum 39 3.3 Source-pipe-detector configuration 44 3.4 Illumination of stratified flow regime 44 3.5 Pipe illumination for various beam widths 48 3.6 Spectrum of upper detector versus beam width 49 3.7 Spectrum of lower detector versus beam width 50 4.1 Conventional CT scanner 60 4.2 Monte-Carlo backprojection 67 4.3 Spectra from Monte-Carlo backprojection 68 4.4 Flow chart of reconstruction algorithms 71 4.5 Example of an ART solution for two pixel image 75 4.6 Hot-spot and cold-spot test configuration 80 4.7 Hot-spot spectra 81 4.8 Cold-spot spectra 82 vii Page 4.9 Hot-spot reconstruction 84 4.10 Cold-spot reconstruction 85 4.11 Convergence of reconstruction algorithms 87 5.1 Calculated spectrum using incorrect source strength value 91 5.2 Reconstruction using incorrect source strength value 93 5.3 Experimental setup 95 5.4 Flow regime models 99 5.5 Singly scattered photon spectrum 100 5.6 Reconstruction with and without noise 101 5.7 Spectra with noise added 103 5.8 Detected spectrum with pipe wall effects 117 6.1 Spectrum from a 4 inch pipe 122 6.2 Detail of spectrum from the pipe 123 6.3 Effect of signal-to-noise ratio 125 6.4 Detector response versus bubble size 127 1.1 Example geometry of photon trajectory 132 11.1 Detector response for diverging source beam 135 11.2 Detector response for parallel source beam 136 11.3 Detector configuration for multiple scatter versus detector position 138 11.4 Spectra for multiple scatter versus detector position 139 V.l Spectra from annular flow 151 V.2 Spectra from annular flow . 152 V.3 Spectra from annular flow 153 V.4 Spectra from inverted annular flow 154 viii Page V.5 Spectra from inverted annular flow 155 V.6 Spectra from inverted annular flow 156 V.7 Spectra from bubbly flow 157 V.8 Spectra from bubbly flow 158 V.9 Spectra from bubbly flow 159 V.10 Spectra from stratified symmetric flow 160 V.ll Spectra from stratified symmetric flow 161 V.12 Spectra from stratified symmetric flow 162 V.13 Spectra from stratified asymmetric flow 163 V.14 Spectra from stratified asymmetric flow 164 V.15 Spectra from stratified asymmetric flow 165 V.16 Annular flow pattern 166 V.17 Reconstruction of annular flow 167 V.18 Reconstruction of annular flow 168 V.19 Reconstruction of annular flow 169 V.20 Inverted annular flow pattern 170 V.21 Reconstruction of inverted annular flow 171 V.22 Bubbly flow pattern 172 V.23 Reconstruction of bubbly flow 173 V.24 Reconstruction of bubbly flow 174 V.25 Reconstruction of bubbly flow 175 V.26 Stratified symmetric flow pattern 176 V.27 Reconstruction of stratified symmetric flow 177 V.28 Stratified asymmetric flow pattern 178 V.29 Reconstruction stratified asymmetric flow 179 V.30 Example of smooth convergence 180 IX Page V.31 Example of divergence 181 V.32 Example of divergence 182 V.33 Example of oscillating divergence 183 V.34 Example of non-convergence 184 V.35 Setup for laboratory experiments 185 V.36 Experimental spectrum of annular flow 186 V.37 Experimental and Monte-Carlo spectra 188 V.38 Experimental spectrum of stratified flow 189 V.39 Experimental spectrum of bubbly flow 190 x Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TOMOGRAPHIC TWO-PHASE FLOW MEASUREMENT USING COMPTON SCATTERING OF GAMMA RAYS By DAVID E. BODETTE December, 1986 Chairman: Dr. Alan M. Jacobs Major Department: Nuclear Engineering Sciences Using Compton scattered gamma-rays to measure local void fraction was first suggested by Kondic and Hahn in 1970. The Compton scattered gamma-ray densitometers they suggested employ a single, narrow beam source and either a well collimated detector or an uncollimated detector. The collimated detector configuration only gives the local void fraction measurement in the small volume where the source and detector collimators intersect. The uncollimated detector configuration is a more efficient design since the local void fraction along the source beam's path is measured in a single reading. A logical extension of the technique is to use wide beam illumination and uncollimated detectors to sample an even greater portion of the flow cross section. This report investigates and demonstrates several methods of inferring two-phase flow parameters using wide beam illumination coupled with two detectors placed xi symmetrically about the source and pipe. The spatial distribution of the fluid in a slice of the pipe is encoded with respect to energy in the singly scattered photon flux. Two basic techniques are detailed for decoding the spatial information: the method of spectral moments and the method of computed tomography (CT). Examination of the low-order moments of the singly scattered photon spectra from the two detectors provides sufficient information for flow regime identification. Flow asymmetries are revealed by comparison of the spectral moments of the two measured spectra. Further classification of the flow pattern is made on the basis of the first and second moments of the spectra. The real focus of this report is the adaptation and successful demonstration of CT techniques with Compton scattering. Three series expansion techniques are considered: the algebraic reconstruction technique (ART), the simultaneous iterative reconstruction technique (SIRT), and the iterative least squares (ILS) technique. Application of the modified ART, SIRT, and ILS algorithms to a hot-spot and a cold-spot reconstruction problem indicates that SIRT and ILS are more accurate than ART. A more extenÂ¬ sive testing of the ILS algorithm for a variety of model flow regimes demonstrates the potential of Compton scatter tomography as a quantitative measurement technique. xii CHAPTER 1 INTRODUCTION The case for developing improved non-intrusive measurÂ¬ ing technologies for sensing material/phase distribution in two-phase flow has been well argued by many authors. The problem is one which spans a broad number of fields incluÂ¬ ding manufacturing, engineering, security, and medicine. Each potential application imposes its own peculiar demands- which are often at odds with constraints in other applicaÂ¬ tions. This situation has resulted in many diverse techniques being pursued and often optimized for each specific problem. Examples of the technologies employed include ultrasonics, microwaves, lasers, optics, and irraÂ¬ diation with particles, X rays, or gamma rays. One specific problem encountered in the nuclear power industry and its associated research industry is to measure the void fraction and/or liquid distribution (flow regime) of two-phase flows inside steel walled pipes. Any intrusive instrument which disturbs the flow pattern is certainly not desirable, but an even more strict requirement often encounÂ¬ tered in these circumstances is that penetrations through the pipe wall are not permitted. Another requirement of any measuring instrument is that its output should be unambiÂ¬ guous. In the case of two-phase flow void fraction 1 2 measurement this means the entire flow cross section must be interrogated. There may also be physical limitations in terms of size and orientation, as well as access to the pipe. A perusal of the recent literature on measurements in two-phase flow [1, 2, 3, 4] reveals that the choice of techniques for making non-intrusive local void fraction and/or phase distribution measurements is limited to irraÂ¬ diation with particles, gamma rays, and laser light. When intrusion is not an issue, a variety of probe techniques are available for local property measurements. Reference [4] lists four basic probe methods which appear to be reliable enough for general use: electrical probes, optical probes, thermal anemometers, and micro-thermocouples. Laser scattering and interferometry is suggested as a very promising method when the flow boundary is not opaque. There are techniques for making average void fraction measurements which might be adaptable to the local void fraction problem. Ultrasonics is capable of measuring the average void fraction in very thick walled pipes, but only for low void fractions (i.e. < 0.20) [3], An array of ultrasonic sensors around a pipe might produce an image of the flow provided the voids are not too big or too many. Microwave resonant cavities are also capable of measuring average void fraction. Using microwaves to sense local void fraction would probably run into difficulties from the metal pipes [1]. Nuclear magnetic resonance (NMR) is currently 3 receiving a lot of attention in the medical field and certainly has potential for two-phase flow. The medical NMR machines involve a radio-frequency field and a large magnetic field which cause problems in nearby electronic equipment. The metal pipes seem to be another source of potential problems in implementing NMR for local . measurements; however, plastic pipes can be used to avoid such problems [5], Gamma-ray attenuation is the predominant method of choice for making average void fraction measurements and is often used as the calibration standard for other methods. Gamma-ray attenuation instrumentation and techniques are recognized as having reached the highest level of development. The reason for this dominance is that the method places a minimum of constraints on the user. The method is non-intrusive, fairly insensitive to flow pattern, and functions over a broad range of void fractions, flow regimes, and temperatures. There are many ways of utilizing gamma rays or neutrons to measure local void fraction and phase distribution. The basic methods of implementing tomographic two-phase flow measurement with gamma rays or neutrons can be grouped into the three following categories: 1) Measure the emission of gamma rays or neutrons from- activation products in the fluid, induced excitation, or injected radiotracers. 2) Measure the attenuation of narrow or wide source beams in traversing the pipe. 4 3) Measure the singly scattered gamma rays or neutrons from a primary source beam. Methods from category 1 could be roughly classified as emission computed tomography (ECT) techniques. An example of a category 1 method is a single photon ECT system which examines the ~6 MeV gamma rays from the activation-decay scheme pair 016(n,p)N16 n16 + 016 + Y + 0â€œ. Although recent advances in the analytical techniques of ECT and in instrumentation suggest this as a viable method, it appears that only cursory attention has been given to the technique [6]. An ECT system could also be employed using injected radiotracers; however, the quantity of radioÂ¬ nuclides required poses a serious drawback [7]. The most familiar tomographic instruments comprise category 2 and are generally referred to as computed axial tomography (CAT) scanners. The basic hardware setup of CAT scanners involves placing a source on one side of the pipe and a detector on the opposite side of the pipe. The whole source-detector assembly rotates in a controlled manner around the pipe while the attenuation of the source beam at different orientations is recorded. The rate at which CAT scanners can rotate about the pipe is a limiting factor. CAT scanners of the newest generation do not employ any 5 moving parts, but instead use a magnetically focused flying spot electron beam and a fixed circular detector array [8]. These fifth generation CAT scanners can acquire an image in approximately 1/10 second and are very promising where cine tomography is necessary. The capital investment in these new CAT scanners is significant and combined with the fact that they are not compact or portable means they are not a viable alternative in many applications. Category 3 techniques are closely related to ECT except that the emitted radiation does not come from radioactive material within the fluid, but from radiation scattered out of an external illuminating source beam. Gamma-ray and neutron sources are possible candidates, and methods based on both are actively being developed [9]. The earliest reported efforts with scattered radiation focused on gamma rays scattering and involved measuring the local density in a small volume [10]. A well collimated source illuminated an object and a well collimated detector measured the Compton scattered gamma rays coming from the small volume where the source and detector collimators' fields-of-view intersected. The output of the detector was proportional to the electron density in the small volume which in turn is proportional to the material density for a constant material composition. The Compton scattering method of local density measurement developed, at first, independently in the medical and industrial imaging fields. Lale published two articles, one in 1959 [11] and another in 6 1968 [12], which demonstrated how internal organs could be imaged by exploiting the Compton scattering of gamma rays. In 1970, Kondic and Hahn [13] published the first suggested engineering application of Compton scattering, that of making local density measurements in two-phase flow. They proposed both a collimated source/collimated detector (as Lale) and a collimated source/uncollimated detector arrangeÂ¬ ment. The uncollimated detector offered the improvement that local density values were given for points along the illuminated path of the source beam from a single measurement. Acting independently, Farmer and Collins [14, 15] proposed and demonstrated the collimated source/uncolliÂ¬ mated detector arrangement in determining anatomical cross sections. Their images were only qualitative since the calibration procedures used failed to completely remove the effects of multiple scattering or variations in attenuation. The analytic method for making quantitative measurements using this source-detector arrangement by deconvolving the substantial multiple scatter blur was developed by Anghaie [16, 17]. The use of Compton scattering to make local density measurements in two-phase flows was pursued earnestly in the late 70's and in to the 80's. Almost exclusively the focus was on the collimated source and either collimated or uncolÂ¬ limated detector arrangement. Tomographic images produced with either of these arrangements required scanning the source and/or detector to cover all the points in the plane of interest. Kondic [18] in 1978 first proposed using an uncollimated or wide beam source coupled with an uncolliÂ¬ mated detector to image an entire flow cross section. Such an arrangement made more efficient use of the source photons than the previous source-detector arrangements and required no moving parts. Both of these attributes made the arrangeÂ¬ ment more adapted to the high framing rates necessary for transient analysis. However, no critical study or viable implementation of the wide source beam/uncollimated detector arrangement has appeared in the literature. Imaging with scattered neutrons was developed in the late 7 0's and early 8 0's independently of and somewhat concurÂ¬ rently of the present work with scattered gamma rays. The use of neutrons to make local void fraction measurements was suggested by Banerjee and Lahey [19] and separately by other researchers. Neutron scattering tomography was developed by Hussein [20] as his Ph.D. thesis. He used a wide parallel beam of 14 MeV neutrons to illuminate a simulated two-phase flow and then measured the fluence and energy of the scattered neutrons using small organic scinti- lator detectors at different positions around the flow volume. The focus of the present dissertation is the developÂ¬ ment and demonstration of the algorithms for reconstructing the local density distribution of a two-phase flow using Compton scattered gamma rays. The specific experimental setup considered uses a wide-beam collimated gamma-ray 8 source to illuminate the flow and energy sensitive detectors to examine the emerging singly Compton scattered photon flux. To date, the reconstruction techniques proposed for this setup have been based on deterministic algebraic soluÂ¬ tions and have not addressed the problems of multiple scatÂ¬ tering, statistical error, image blur, and other errors. The contribution of this author is to apply the powerful reconstruction algorithms developed in the medical and optical imaging fields to this problem and to critically assay the effects ofâ€žmultiple scattering and statistical errors on the reconstruction. Tomographic imaging with Compton scattered gamma rays has two distinct features which distinguish it from previously encountered tomographic reconstruction problems. First, the imaging processes dealt with in the medical and optic fields have been linear or could be linearized by an r analytic technique. Imaging the phase distribution in a two-phase flow with Compton scattered gamma rays is a nonÂ¬ linear problem. Second, the combination of using a point source and the Compton scattering process introduces a unique curved geometry to the flow reconstruction process. The forward problem of predicting the emerging singly Compton scattered photon flux given the local density disÂ¬ tribution is discussed in Chapter 2. Chapter 3 details a method of inferring tomographic information about a two- phase flow by examining the spectral moments of the Compton scattered photons [21]. Chapter 4 develops several more sophisticated tomographic reconstruction techniques and compares their performance under ideal conditions using data generated by Monte-Carlo calculations. Chapter 5 evaluates Compton scatter tomography under more realistic conditions and discusses some of the considerations in designing a system. Chapter 6 summarizes the significant results of this study and offers directions for future work. CHAPTER 2 THE FORWARD PROBLEM Fundamental to the problem of reconstructing the phase or density distribution using Compton scattered gamma rays is being able to model the gamma ray interaction with the density field. The aim of this chapter is to describe the basic equations predicting the emerging scattered photon flux from a known density distribution illuminated by a point qamma-ray source. The basic concepts of using gamma- ray scattering to make local density measurements are introduced first with the collimated point source/collimated detector arrangement. Once the equations for analyzing this setup are developed, the analysis is extended to collimated point source/uncollimated detector arrangements by introducÂ¬ ing the equivalent approach of energy collimation. The analysis is extended further to the wide-beam collimated point source/uncollimÃ¡ted detector arrangement. Lastly, parallel source beam illumination is considered and the differences between it and the diverging wide-beam illumination are discussed. Throughout, a detailed analysis is only given to the singly scattered photon flux. The presence of the multiply scattered photon flux is included in all the equations and discussed in general terms. However, for the case of the wide-beam collimated source 10 11 /uncollimated detector configuration a specific method is given for removing the multiple scatter component. Collimated Point Source/Collimated Detector Figure 2.1 illustrates the arrangement of a collimated source and collimated detector for making local density measurements. Monoenergetic gamma rays leave the source and are collimated to a thin pencil beam which illuminates a chord of the pipe. A similar collimator on the detector is aligned so that its field-of-view intersects the source beam at the small volume denoted by V. The detector senses the source gamma rays scattered in the volume V. The detector response, D(E), is D(E) = eSÂ°puwAiA0Az + M(E) (2.1) where SÂ° is the source activity; p is the density of the fluid in the volume V; pw is the incoherent differential scattering cross section; Az is the path length of source photons in V; and AQ are the total attenuation functions, including interaction attenuation and the solid angles of view, along the incoming and outgoing paths respectively; e is the detector efficiency; and M(E) represents the contribution to the spectrum of multiply scattered photons. 12 Figure 2.1: A collimated source/collimated detector gamma- ray scattering configuration for measuring the density in the volume V. 13 The functional form of both attenuation functions is Ax = Cxâ€¢exp(-J pxp(z)dz) (2.1a) where Cx contains the geometric or field-of-view factors, Px is the mass attenuation coefficient, p(z) is the local density at z, and x is a dummy subscript representing i and o for the incoming and outgoing paths, respectively. The line integral in equation (2.1a) is carried out either along the incoming path, 1^, from the source to the volume V or along the outgoing path, lQt from V to the detector. The value of the integral is in units of optical path length. The predominant interaction of gamma rays with energies between 0.1 and 1.0 MeV in water is Compton scattering. The energy of the singly Compton scattered photons arriving at the detector is related to the polar scattering angle, Â©, by the kinematic relation Eq = EÂ±/{1 + (Ei/mec2)(l - cos 0)) (2.2) where EQ is the energy of the scattered photon, Ej_ is the energy of the source photon, and mec is the rest mass energy of an electron. The singly scattered photons, represented by the first expression on the right hand side (RHS) of equation (2.1), all arrive at the detector after scattering through 14 approximately the same angle and thus have approximately the same energy. Assuming the attenuation along the incoming and outgoing paths remains constant in time, the detector response to the singly scattered photons is proportional to the density in the small volume V. The multiply scattered photons arriving at the detector, unlike the singly scattered photons, have a broad range of energies since they can suffer any number of collisions (with concomitant changes in energy) prior to reaching the detector. The functional form of the multiply scattered photon component of the detector response is M(E) = ////$(r ',E ',Q ' )K(r ',E ',Q ' ;r^,E,Q)dr 'dE 'dfi 'dr^dft (2.1b) where $(..)is the angular photon flux and K (. -) is the transport kernel. The kernel can be viewed as the probabi- \ lity that a photon at r' going in direction ft' with an energy E' suffers a scattering interaction at r'; the scatÂ¬ tered photon goes in direction ft toward the detector with an energy E; the scattered photon reaches point r^ in the detector's active volume, interacts with the detector and is counted. The integral with respect to r' is carried out only over the detector collimator's field-of-view. The detector collimator serves not only to pick out the volume of interest, V, but to reduce the magnitude of the multiply scattered photon spectral component by restricting the 15 region of space where the photons can suffer their last collision prior to reaching the detector. Collimated Point Source/Uncollimated Detector The function of the detector collimator in picking out the small volume V can also be achieved by electronic colli- mation. An example is to use a high energy resolution detector coupled to a multi-channel analyzer (MCA) instead of the detector collimator. The point of scattering of a singly scattered photon is related to its energy by equaÂ¬ tion (2.2) and the polar scattering angle as determined from the measurement geometry. The singly scattered photons arriving at the detector now have a range of energies. The detector response becomes D(E) = SÂ°p(x) vi^ÃxjAQUJeÃE) + M(E) (2.3) where SÂ°, p, yw, Aj_, AQ/e , and M are as defined for equaÂ¬ tion (2.1) and x is the point of scattering. The x-axis is colinear with the source beam illumination path (Figure 2.2). Note that x and 9 can be expressed in terms of E using equation (2.2) as x = x^ - l - cos^e cos 9=1+ (mecz/E^)(1 - E^/E) (2.4) 16 Source Pipe Configuration of a collimated source/uncollimated detector gamma-ray scattering densitometer. Figure 2.2: 17 where and are the detector coordinates illustrated in Figure 2.2. The first term on the RHS of equation (2.3) is strictly a function of a single variable for a fixed geometry. The multiple scattering term in equation (2.3) has the same functional form as equation (2.1) except that the inteÂ¬ gration with respect to r' is over all space since there is no detector collimator. Therefore, the magnitude of the multiply scattered photon component of the spectrum is greater with an uncollimated detector than with a collimated detector. Offsetting this increased measurement structured noise is the fact that the singly scattered spectrum now contains the local density information for the entire source beam illumination path. The information which would have required many measurements with a collimated detector is available from a single measurement with an uncollimated detector. Uncollimated Point Source/Uncollimated Detector Just as eliminating the detector collimator increases the information contained in a single measurement, so does eliminating the narrow source collimator. An uncollimated or a wide angle-of-view collimated source can illuminate the entire flow cross section. The detector response then contains the local density information from every point in the flow cross section. The manner in which the local density information is encoded in the detector response, 18 however, complicates the data processing necessary to reconstruct the local density distribution. Figure 2.3 illustrates a wide beam source Compton scatÂ¬ tering densitometer setup. The curves drawn across the pipe represent lines of constant scattering angle, called isogonic lines by Kondic [18]. Primary source photons scattering at points along an isogonic line all scatter through the same angle in order to reach the detector. Using elementary geometry it is easy to show that the isogonic lines are all segments of circles which pass through both the source and the detector. The axis labeled Y in. Figure 2.3 is the loci of all the centers of these circles. The radius, r, and the location of the center- point, yc, of each isogonic circle as a function of the scattering angle, 9, are r = I SD sin Â© 2 yc = - -2SD cot Â© (2.5) where SD is the distance between the source and the detector. The detector response is now represented by an integral equation ME) D ( E ) = / D(x,y)dÂ£ + M(E) (2.6) 19 ELEVATION VIEW Figure 2.3: Configuration of an uncollimated source/ uncollimated detector gamma-ray scattering system for tomography. 20 where M(E) is the multiply scattered photon spectrum, D(x,y) is the contribution to the detected photon spectrum of singly Compton scattered photons scattered at point (x,y), and dÂ£ is the differential line element along the isogonic contour Â£(E). The first term on the RHS of equation (2.6) is essentially the integral of the first term on the RHS of equation (2.3). The second term, as before, represents the multiply scattered photon contribution to the spectrum and has the same functional form as equation (2.1b). The multiple scatter component can comprise approxiÂ¬ mately one-third of the total signal [22, 23] and must be removed. It may be possible to remove the multiple scatterÂ¬ ing component analytically since the transport kernel is a known function and the angular photon flux can be generated from an order-of-scattering calculation once a good estimate of the density distribution is obtained. However, this would be a cumbersome process and does not appear to be necessary as there appears to be an acceptable alternative approximation. A multiply scattered photon arriving at the detector with energy E could have come from any point within a large volume of the pipe. The structure of the fluid distribution imposed on the multiply scattered photon spectrum would tend to be averaged out because of the many possible paths that the photon could have taken prior to arriving at the detector with energy E. The result is that for wide beam source/uncollimated detector arrangements the multiply scattered photon spectrum should be smooth and a simple background subtraction procedure can be employed to remove it. The validity of this intuitive argument has been tested by a series of Monte-Carlo calculations (see Appendix I). The multiply scattered photon spectrums have been comÂ¬ puted for five flow regimes at both a high and low void fraction. The calculations assume that the flow is contained in a pipe whose thickness and composition are such that it can be neglected. The results, summarized pictorially in Table 2.1, show that to a first order approxÂ¬ imation equation (2.6) can be rewritten as D (E) = f D (x, y ) dÂ£ + B ( E) . 2- (E) where B(E) is some function of energy which is fitted to the smooth multiply scattered photon spectrum, and D(x,y) is the contribution to the detected photon spectrum of singly Compton scattered photons scattered at point (x,y). The function D(x,y) can be written more explicitly as D(x,y)- SÂ°e(E)yw(Ei,E0)p(x,y)Ai(x,y)AQ(x,y) (2.7) where e(E) is the detector efficiency, Mw is the incoherent differential scattering cross section for water, p(x,y) is the local fluid density at (x,y), 22 TABLE 2.1 Monte-Carlo multiple scatter calculation results. Flow Regime Void Fraction =0.7 Stratified â– D Bubbly â– D r Jy l Annular â– d Inverted Annular â– D Â¿i ,Â«,1-^..^^^^ 23 Table 2.1 continued. - - multiple scatter single scatter 24 A^ and AQ are as defined for equation (2.1), and SÂ° is the source strength. There are several implicit assumptions expressed in equations (2.6*) and (2.7). They are that a) the divergence and width of the source beam in the direction perpendicular to the imaging plane is small, b) the source is a point body, c) the detector is a point body, and d) the effects of the detector electronics and MCA are neglected. Careful construction of the source collimator can insure that a) is true. Including the effects of b), c), and d), poses no serious analytic hurdles as far as the forward problem of predicting the detector output given the density distribution. However, these details have serious implications on the resolution obtainable when trying to reconstruct the density distribution from the detector output. ,,Equation (2.4) provides a one-to-one correspondence between the energy of a detected singly scattered photon and the point of scattering. The relationship is true only if the source and detector are point bodies. If the so'^ce ana detector have finite extent, then there is some uncertainty associated with the scattering angle as depicted in Figure 2.4. This means that photons scattering at a point in the density field do not arrive at the detector with a single energy but with an energy over a small range, AE = (3E/3Â©)A0 25 Figure 2.4: Illustration of the scattering angle blur problem arising from the use of non-point source and detector. 26 This "blurring" effect is analogous to taking a picture with an out-of-focus camera. The amount of uncertainty in the scattering angle is not constant over the entire density field nor even along an isogonic line. Consequently, the blurring cannot be expressed as the convolution of some spread function with the spectrum and cannot oe removed by a straignr forward deconvolutional procedure. There is also blurring of the detected photon spectrum by the detector and electronics. Unlike the blurring discussed above, this blurring is generally energy invariant and can be removed by a deconvolutional procedure. Methods for removing the detector system blur and correcting for the energy dependent efficiency are well developed and it is assumed that the measured data are automatically corrected for these effects. Another aspect of the data collection process which affects image quality is the use of a multi-channel analyzer (MCA) to record the photon spectrum. The impact of collecÂ¬ ting the spectrum data in an MCA is that the spectrum is not available as a continuous function but as a discrete set of values. That is, D(E) becomes the set {Dc;c=l...I} where Dc = / D(E)dE. AEc It is now convenient to introduce the isogonic segment as shown in Figure 2.5. An isogonic segment is the region bounded by the isogonic curves which correspond to the Y Figure 2.5: Illustration of an isogonic segment. (Reproduced by permission. Copyright Â® ISA 1985 [21].) 28 energy limits of a single channel in the MCA. Any features of the spectrum smaller than the width of a single MCA channel are smoothed over and cannot be recovered. Generally this does not pose a problem since modern equipÂ¬ ment can collect data with very narrow channel widths. Parallel Beam Source/Uncollimated Detector The analysis so far has dealt exclusively with a small source and a diverging beam. It is relevant to extend the analysis to parallel source beams since work has been done using parallel beam illumination with neutron scattering [20]. Figure 2.6 illustrates one possible way of constructing a parallel beam source from an array of narrow beam sources. Alternatively, a parallel beam source can be seen as the result of moving a wide beam point source infinitely far away from the pipe. As the source is moved away, the source-to-detector distance goes to infinity and by equation (2.5) the radius of an isogonic curve also goes to infinity. The curvature (i.e. the reciprocal of the radius) goes to zero and the isogonic curves become lines or rays which originate at the detector. The exit paths of the scattered photons (or neutrons) are, therefore, colinear with the isogonic lines (see Appendix II). The same equations describing the detector response for the wide beam source/uncollimated detector case are valid provided the appropriate limits on the integrals are used. 29 Discrete Forward Equation This section deals with the practical aspects of expressing equation (2.7) in a form which enables computing the integral of equation (2.6*) using a digital computer. The details of how the computations are done are presented in Appendix III. The purpose hereis to show the relationÂ¬ ship between the discrete equations, which are used in reconstructing an image of a flow cross section, and the "exact" integral form. The measured data are available as a discrete set from which the singly scattered flux spectrum must be extracted. The extraction process of removing the multiply scattered flux component and correcting for the detector system efficiency can be expressed as sc = detector having suffered only a single Compton scattering collision. The notation denotes a vector with components S^. The object here is to be able to calculate given the density distribution in the illuminated cross section of the pipe. An expression for calculating is derived using equations (2.7) and (2.7*) = S Â°f dEj VnpA^dl . AEC ME) (2.9) 30 Beam Figure 2.6: A parallel beam illumination/uncollimated detector gamma-ray scattering system. Isogonic Segment Figure 2.7: An n x n pixel grid is overlayed the flow cross section in order to facilitate calculation or the singly scattered photon flux. Individual pixels are designated 2^Y a single number k ranging from 1 to K(=n ). 31 The first step in numerically evaluating equation (2.9) is to overlay the illuminated flow cross section with an n x n pixel grid as shown in Figure 2.7. Equation (2.9) is rewritten as a sum of integrals over the individual pixels si = SÂ° l / dEf v pkAÂ±A dl (2.10) kc aec Me) where pk(x,y) is equal to p(xfy) for points inside pixel k and is zero elsewhere. The notation kc indicates that the sum is to be carried out only over those pixels in the isogonic segment c. The integrals in equation (2.10) can be expressed equivalently by = SÂ°i Jdx/dyywpkAiA0PSFc(x,y) . (2.11) k The term PSFc(x,y) is the characteristic function of the isogonic segment c defined as PSFc (x'y) 0 if (x,y) (x,y) is inside isogonic segment c is outside isogonic segment c. Any real imaging system imposes a lower limit on the spatial resolution below which variations in the density can not be seen. If the size of the pixels are about the same as the spatial resolution limit, then it is reasonable (even necessary) to assume a constant density within each pixel. 32 The remaining integral in equation (2.11) is approximated by Sc = sÂ° I PkAckPSFck ' (2'12> k The factor AQk is an average quantity defined by Ack //dxdyuw(x,y)AÂ±(x,y)AQ(x,y) J/dxdy k (2.13) where the integrals are carried out only over the area of pixel k. The function PSFck is the fractional area of pixel k which is inside isogonic segment c or PSFck = //dxdy PSFc(x,y). (2.14) k An implicit assumption in equation (2.12) is that the attenuation factors do not vary significantly across a pixel. This assumption is valid except for pixels near an edge in the density distribution and simultaneously on the edge of an isogonic segment. The assumption also breaks down when the size of the pixels are on the order of a photon mean-free-path. In practice there are no sharp edges to the isogonic segments due to the scattering angle blur. The scattering angle blur causes the isogonic segments to "feather" off at the edges (and overlap) as shown in 33 Figure 2.8. The effect is that the assumption works in real applications where non-point source and detectors are used. The definitions for A^ and PSFC^ must be modified to account for the scattering angle blur. The new definition for Acj< specifies that the attenuation factors must not only be averaged over the area of the pixel but also over the range of possible energies with which the photons can reach the detector. The new definition of PSFC^ is the fraction of photons reaching the detector from pixel k with energies in the range of channel c multiplied by the pixel area. Equation (2.12) gives an expression for calculating the singly scattered photon spectrum for a detector looking at the scattered photon flux of a known density distribution illuminated by a wide beam mono-energetic gamma-ray source. Each channel of the singly scattered photon spectrum is seen to represent an integral measure of the density distribution within the corresponding isogonic segment. The goal of the / remainder of this report is to do the inverse operation of extracting information about an unknown density distribution given measurements of the singly scattered photon flux. Chapter 3 tackles this inverse problem by correlating the size and shape of spectra with the different possible flow regimes. Chapter 4 goes one step further in solving the inverse problem and looks at ways of inverting equation (2.12). 34 PSFch(x,y) y The scattering angle blur causes the isogonic segments to feather off and overlap. PSFch(x,y) is the fraction of those photons which scatter at point (x,y) toward the detector that have an energy in the spectrum channel ch. Figure 2.8: CHAPTER 3 METHOD OF SPECTRAL MOMENTS The basic premise of Compton scatter tomography is that the singly scattered photon spectrum contains information about the density distribution. If this is true then each flow pattern should impart some unique characteristics to the spectrum. Examination of the first few moments of the spectrum is a simple pattern recognition technique employed in this chapter to extract these unique characteristics. Description of the Method Consider a Compton scattering setup using a wide beam source and uncollimated detector as illustrated in Figure 3.1. Mono-energetic gamma rays leave the source and undergo Compton scattering interactions inside the pipe. Some of the scattered gammas are received by the detector and the energy spectrum of the detected gamma rays is accuÂ¬ mulated in an MCA. The total count in each MCA channel is proportional to the density of the material in the isogonic segment for the channel. Note that the source beam is stopped down slightly in order to reduce the amount of the pipe wall which is illuminated. The typical two-phase fluid problem consists of a fluid and its vapor. The vapor is, in general, considerably less 35 36 Figure 3.1: Compton scatter densitometry configuration. Scattered radiation from the illumination beam is received by the detector. (Reproduced by permission. Copyright Â® Instrument Society of America (ISA) 1985 [21].) 37 dense than the liquid. Therefore, an MCA channel represents an integral measure of the liquid in the isogonic segment. The singly scattered photon spectrum is a projection of the liquid distribution onto the energy axis. The energy axis can, by suitable transformation, be related to any spatial axis intersecting the isogonic lines. A convenient reference axis perpendicular to the isogonic lines is the Y-axis. The transformation relating the energy axis to the Y-axis is available from equation (2.5) Y(E) = r(0(E))+yc(6(E)) . (3.1) Note that each detector has a characteristic Y-axis. The fluid distribution in the pipe can be inferred from the one dimensional image (i.e. projection) by examining the moments of the spectrum with respect to the Y-axis. The zeroth moment is a measure of the amount of liquid phase present. The first moment is a measure of where the center- of-mass of the liquid phase is. The standard deviation is a measure of how the liquid phase region is distributed about the center-of-mass. The moments are not direct indicators of their respective quantities because of the attenuation factors. The zeroth moment is calculated in a straight forward manner as 38 c (3.2) where is a channel of the singly scattered photon spectrum [see equation (2.8)]. The first moment and standard deviation are calculated according to the formulae and c AY YdY c Yc+1 (3.3a) l si c AY, fY, Y - â– c+1 (3.3b) I si Figure 3.2 graphically depicts the zeroth moment, first moment, and standard deviation for annular and inverted annular flow at 50% void fraction. The two peaks in the spectrum for the annular flow correspond to the liquid in the illuminated portions of the annulus separated by the low density vapor core region. The peak heights are not the same because the source beam is attenuated before it reaches the far side of the annulus which results in fewer scatÂ¬ tering collisions taking place there. Geometric attenuation RELATIVE DETECTOR COUNT 39 Figure 3.2: Spectral moments of the singly scattered photon spectrum. a) Spectrum for an annular flow; b) Spectrum for an inverted annular flow. (Reproduced by permission. Copyright Â© ISA 1985 [21].) 40 1.00 o o cc o H O Ui H UJ O Ui > H < 0.75- 0.50 0.25 UJ CC < yO > - 0.340 reÃ. units r1> \ I' ,0* | \ aA I. I I, i \ i , V 10 15 Y (cm) 20 (b) Figure 3.2â€”continued 41 of the photon fluence (i.e. photons/area) for points further from the source also occurs. Similarly, these same attenuaÂ¬ tion terms are applicable with respect to the detector. Only a single peak for the liquid core region is seen in the inverted annular spectrum. The peak is not symmetric due to attenuation as the gamma rays pass through the liquid and to geometric attenuation. The zeroth moments are not identical even though the void fraction is the same. All of the liquid core region is illuminated in the inverted annular flow. Not all of the liquid annulus in the annular flow is illuminated. In addition, the different liquid distribution leads to differÂ¬ ences in attenuation, solid angle, and differential scattering cross section. Consequently, the zeroth moment of the inverted annular flow is larger than that of the annular flow. Both the annular and inverted annular flows are symmetrically distributed within the pipe and so, the first moment would be expected to lie at the center of the pipe. However, due to the attenuation effects the first moment is skewed closer to the source and detector. The standard deviation of the annular spectrum is larger than that of the inverted annular. This is expected since the annular spectrum consists of two widely separated peaks while the inverted annular spectrum is a more compact single peak. 42 Procedure and Results Monte Carlo transport calculations are carried out to test whether flow regime identification is possible by examining the first few moments of the singly scattered gamma ray spectrum. Several different flow regimes are analyzed; in each case the energy spectrum from single scattering interactions with liquid or vapor inside a pipe is calculated for two detectors. The detectors are assumed to be placed symmetrically with respect to the pipe and source as shown in Figure 3.3. The source-pipe-detector geometry is the same for all cases. The source and detector are treated as points for purposes of calculating the scattering angle. However, for purposes of calculating the solid angle view of each interaction, the detector has a radius of 0.5 cm. Detected scattering interactions within the pipe wall are not considered since the primary focus is to evaluate the gamma rays coming from the fluid. However, attenuation due to the pipe is taken into account. The flow regimes considered are annular, inverted annular, stratified, bubbly, and mist. The annular and inverted annular flows are modeled as a liquid/vapor annulus and a vapor/liquid central core concentric with the center of the pipe, respectively. The stratified flow is modeled as a pipe half-filled with liquid as illustrated in Figure 3.4. The bubbly and mist flows are simulated by choosing fifteen randomly sized (0.85 43 bubbles of vapor/liquid in a matrix of liquid/vapor, respecÂ¬ tively. New bubbles are periodically chosen during a calÂ¬ culation to simulate stochastic flow conditions. The void fraction is fixed at 50% for the annular, inverted annular, and stratified flow. Due to the random nature of choosing the bubbles, the bubbly flow has a void fraction of ~45% and the mist flow has a void fraction of ~58%. Four calculaÂ¬ tions are made for each flow regime where the width of the source beam in the plane of the pipe cross section varies from 4.5 to 28 degrees. The results of the moment calculations are tabulated in Table 3.1. Interpreting the Moments The moments are functions of the magnitude, shape, and width of the spectra. The spectra change as the source beam width varies since the amount of the flow cross section illuminated changes. Stating general trends of the spectra based on the changing amount of illuminated fluid is diffiÂ¬ cult due to the effects of attenuation. However, the following two observations can be made. First, the magniÂ¬ tude of the spectrum from a narrow source beam is smaller than that from a wider source beam. Assuming the source strength remains constant, the number of photons illuminaÂ¬ ting the pipe is less for the narrow source beam. The result is that the zeroth moment becomes smaller as the source beam narrows. Second, the width of the spectrum decreases as the source width is narrowed. As illustrated in Figure 3.5, the 44 Source Figure 3.3: One source, two symmetrically positioned detector configuration used for the Monte Carlo studies presented in this chapter. (Reproduced by permission. Copyright Â® ISA 1985 [21].) Figure Vapor Liquid 3.4: Source-detector configuration relative to the stratified flow regime case. (Reproduced by permission. Copyright Â® ISA 1985 [21].) TABLE 3.1 Zeroth moment, first moment, and standard deviation of the singly scattered photon spectra for the five flow regimes listed. Source Beam Width 2 8Â° 14Â° 8Â° 4.5Â° Flow Regime [relative] [cm] A 0 V A o >i V A Q M ANNULAR Detector 1 0.311 15 12 0.0935 15 15 0.0485 15 14 0.0267 15 14 Detector 2 0.313 15 12 0.0935 15 15 0.0484 15 14 0.0267 15 14 INV. ANNULAR Detector 1 0.337 16 3.7 0.194 16 3.7 0.108 16 3.5 0.0601 16 3.5 Detector 2 0.340 16 3.7 0.193 16 3.6 0.108 16 3.6 0.0600 16 3.5 STRATIFIED Detector 1 0.621 17 6.1 0.360 17 6.6 0.207 16 6.9 0.116 16 6.7 Detector 2 1.00 14 5.4 0.353 15 6.3 0.158 16 6.6 0.0789 16 6.4 BUBBLY Detector 1 0.434 15 8.4 0.204 16 7.6 0.108 16 8.1 0.0597 16 8.2 Detector 2 0.436 15 8.4 0.205 16 7.5 0.107 16 8.0 0.0602 16 8.2 MIST Detector 1 0.431 16 6.6 0.183 16 6.5 0.108 16 5.9 0.0614 16 5.3 Detector 2 0.423 16 6.9 0.184 16 6.7 0.112 16 5.7 0.0600 16 5.2 Note: The spectra are from a Cs-137 source illuminating Ã a 12 i cm diameter pipe with a 0 .75 cm thick steel wall . The beam divergence normal to the flow cross section is 5 degrees in all cases. The void fraction is nominally 50% for all cases. 46 region within the isogonic segments associated with the extreme channels of the spectrum is no longer illuminated as the beam narrows. The effect is more pronounced for the isogonic segments at low values of the Y-coordinate. If liquid is present in the affected regions, the result is that the spectrum narrows and the narrowing can be more pronounced at the low end of the Y-axis. Such an asymmetric change in the spectrum width means the first moment shifts to larger values of Y as the beam narrows. A narrower spectrum would also be expected to yield a smaller standard deviation. The magnitude of change in the first moment or standard deviation, or whether any change occurs at all, depends on any accompanying change in the shape of the spectrum. If there is no liquid in the affected regions then, of course, there is no change in the width of the spectrum. A Detailed Example '' The following is a discussion in physical terms of how the moments vary with changes in the source beam width. So that sufficient attention can be given to the changes in each moment, the discussion is confined to only the stratiÂ¬ fied flow regime. Enough insight should be gained to interpret the data in Table 3.1 for the remaining flow regimes. Figures 3.6 and 3.7 are the spectra for the upper and lower detector, respectively, viewing a stratified flow. These figures show the progressive changes in the magnitude, 47 shape, and width of the singly scattered photon spectrum as the source beam width narrows. The magnitudes of the spectra decrease as the beam narrows because there are fewer photons illuminating the flow. The magnitudes of the two figures should not be compared since the vertical scales are diffeÂ¬ rent. Instead the figures should be compared on the basis of the zeroth moments. The zeroth moment of the lower detector's spectrum is larger when the source beam is wide, while the zeroth moment of the upper detector's spectrum is larger when the source beam is narrow. The behavior of the zeroth moments can be explained by considering the effect of attenuation on the scattered gamma rays and the differences in solid angle each detector views of the scattering interactions. When the source beam is widest the average attenuation of the scattered gamma rays reaching each detector is about the same. The average position of scattering occurs in the middle of the liquid and the scattered gamma rays must pass through almost the same amount of liquid and steel to reach either detector. When the source beam is narrowed the average scattering interaction takes place toward the center of the pipe and the scattered gamma rays must pass through more liquid to reach the lower detector. Thus, attenuation of the scattered gamma rays explains why the zeroth moment of the upper detector is larger when the source beam is narrow. When the source beam is wide the zeroth moments should be the same since the average attenuation is the 48 Figure 3.5: Darkened regions are those isogonic segments which are eclipsed when the more restricted beam illumination is used. (Reproduced by permission. Copyright @ ISA 1985 [21].) RELATIVE DETECTOR COUNT 49 Figure 3.6: Detected photon energy spectra of once- scattered gamma rays for upper detector viewing a stratified flow as in Figure 3.4. (Reproduced by permission. Copyright Â® ISA 1985 [21].) RELATIVE DETECTOR COUNT 50 Figure 3.7: Detected photon energy spectra of once- scattered gamma rays for the lower detector viewing a stratified flow as in Figure 3.4. (Reproduced by permission. Copyright Â© ISA 1985 [21] . ) 51 same. The reason they are not is because the upper detector is further away and presents a smaller solid angle to each scattering interaction. The behavior of the first moments can be explained in the following manner. Observe that the spectrum for the widest source beam is peaked at about 13 cm for the lower detector and at about 17 cm for the upper detector. Due to the source-pipe-detector geometry the isogonic segments do not all intersect the liquid region. The projection of the liquid region onto the corresponding Y-axis for the two detectors is not symmetric. As the source beam narrows both spectra become more uniform since approximately the same amount of liquid is illuminated in each isogonic segment. The remaining slight peak is due to the non-uniform amount of attenuation and to solid angle considerations. The more uniform a spectrum the more the center-of-mass coincides with the center of the spectrum. The result is that the first moment shifts to higher Y values for the lower detector and shifts to lower Y values for the upper detector. Note also the preferential narrowing of the spectra which tends to enhance the changes in the first moments. The standard deviation for both detectors' spectra increases uniformly as the source beam narrows from 28 to 8 degrees and then decrease slightly as the beam narrows further to 4.5 degrees. When the source beam is widest both spectra are peaked near the center-of-mass which gives a 52 relatively small standard deviation. The standard deviation increases as the source beam narrows because the spectra become more uniform. The slight narrowing of the spectra leads to a smaller standard deviation. The two competing effects result in the observed behavior. Using the Moments The purpose of calculating the moments is to identify the flow regime automatically and without any a priori knowledge. Inter-comparing the moments of all the flow regimes at each particular source beam width suggests that using the widest source beam offers the greatest discriminaÂ¬ tion between flow regimes. A close agreement between the moments of the upper and lower detectors is seen in all the data for the annular, inverted annular, bubbly and mist flows where the fluid distribution is symmetric. The stratified flow is asymmetric, but only for the wider source beams (i.e. 14 degrees) is the asymmetry discernible from comparing the moments of the upper and lower detector. The annular flow is discernible from the other symmetric flows with the wide source illumination because of its large standard deviation. The inverted annular flow conversely is distinguishable because of its small standard deviation. Discriminating between the bubbly and mist flows is not possible on the basis of the first moment or standard deviation due to their similar nature. Bubbly flow is expected when the void fraction is low (i.e. < 50%) while mist is expected when the void fraction is high. The zeroth 53 moment is sensitive to the void fraction and could be used to determine which flow exists. The results presented so far show that flow regime identification can be made using the singly scattered photon spectra from two detectors sensing the emerging scattered photon flux from a flow cross section illuminated by a wide beam source. A tomographic image of the flow can be produced using a simple model of the flow cross section for each flow regime. The consistency between the actual flow regime and the flow regime as identified by the spectral moments can be checked using the solution to the forward problem [i.e. equation (2.12)]. The calculated singly scattered photon spectra for the two detectors based on the identified flow regime should match the measured spectra. If the calculated and measured spectra are not close enough, then the flow cross section model is adjusted accordingly. For example, if the spectral moments indicate that the flow is annular, then the size of the annulus is adjusted until the measured and calculated spectra match. The next chapter develops other methods of producing such tomographic images of the flow. These other methods are more systematic in adjusting the flow cross section image so that the measured and calculated singly scattered photon spectra match. CHAPTER 4 THE RECONSTRUCTION PROBLEM A general statement of the problem at hand is Given a set of singly scattered photon flux measureÂ¬ ments, {S1}, made using a wide-beam source and uncollimated detectors, reconstruct the liquid-vapor phase distribution {p} that produced {S^. Ideally the set of flux measurements is acquired with a single exposure. That is, the source and detector locations are fixed and there are no moving parts. Insight into this problem can be gained by considering the following formulation of the forward problem from equations (2.6*) and (2.7): Sl=[T]p (4.1) where [T] is the matrix which mathematically describes the scattering transport process, Â§â– *â– is the vector of singly scattered photon flux measurements, and p is the vector of pixel densities. The obvious solution of inverting [T], that is, p=[T] 1 S1 (4.2) 54 55 is not strictly possible since the elements of [T] are functions of p. The problem at hand is nonlinear. Until now the general approach has been to avoid dealing directly with the nonlinearity of [T], One often-proposed method (e.g. [24, 25]) is to use multiple source-detector locations and/or source energies in order to cancel out the nonlinear aspects of [T]. There are numerous practical limitations in implementing this approach for dynamic measurements. Another proposed method is to assume a density distribution based on prior experience and use this in computing [T]. Since [T] is a function of the density, a bias in p is likely for the high contrast reconstructions expected from two-phase flows. A characÂ¬ teristic of these methods is the idea of calculating the density distribution in an algebraic, once-through fashion. It is proposed here that the problem be viewed as a tomographic reconstruction problem where the measured data, s} represent projections of the phase distribution, p, onto non-physical energy axes. The projection data are incomplete since only a few detectors are employed. SimilaÂ¬ rly the elements of [T] are at best known only approximately since they depend on the unknown phase distribution p. When the available information is incomplete and/or imperfect, an algebraic solution may not exist or several solutions may exist. These are familiar situations in CT and the analytic techniques are already developed to handle them. The analoÂ¬ gous projection operator, [T], of current CT applications is 56 mathematically well-behaved (i.e. linear, single-valued, etc.). The CT reconstruction methods must, therefore, be modified to handle the nonlinearities of Compton scatter tomography (CST). Several CT reconstruction techniques are adapted to CST in this chapter. General Solution Methods The basic methods of CT for solving the reconstruction problem can be grouped into five categories: 1. Fourier transform: exponential Radon transform [26], filtered backprojection [27]. 2. Direct scanning [12, 14, 15]. 3. Matrix inversion: successive approximation [20], maximum entropy (see Appendix IV) [28], generalized inverse, and pseudo inverse [29]. 4. Series expansion [30]: algebraic reconstruction technique (ART), simultaneous iterative reconstrucÂ¬ tion technique (SIRT), and iterative least squares (ILS). 5. Statistical: Monte-Carlo backprojection. The transform techniques are recognized as being the more elegant, mathematical reconstruction methods. TransÂ¬ form methods in general require a large number of projecÂ¬ tions and do not show good noise immunity. The cost of detectors and the physical limitations in mounting and cooling the detectors imply a limit on the number of projecÂ¬ tions available. Practical limits on source size, limits on detector size, and desired framing rate are all going to 57 result in statistical noise in the measurements. These features of using Compton scatter to image dynamic two-phase systems make transform methods unsuitable as a reconÂ¬ struction technique. Direct scanning techniques inherently require a narrow collimated source beam and a moving source detector arrangeÂ¬ ment. A narrow collimated source is less photon efficient than a wide collimated source. This inefficiency coupled with the time necessary to translate the source and detecÂ¬ tors means that a direct scanning system is slower than a comparable wide-beam source, uncollimated detector system. Since one goal is to measure dynamic two-phase flow, the faster measurement system is the obvious choice. There is another reason for dismissing direct scanning techniques that involves accounting for the multiply scattered photon flux. The structure of the flow regime is more strongly imprinted on the multiply scattered photon flux when a narrow-beam source is used than when a wide-beam source is used. Any structure in the multiply scattered photon spectrum means that approximating the multiply scattered photon spectrum by a general function B(E) [see equation (2.6*)] is not possible. Removing the multiple scatter component in direct scanning systems, therefore, requires more computation and is more prone to error. Matrix inversion and series expansion reconstruction techniques are equivalent in that a solution is sought by iteratively improving an initial guess of p (and hence [T]) 58 until some integral transformation of the reconstructed phase distribution and the true phase distribution are judged close enough with respect to the measurement error. The difference between matrix methods and series expansion methods is that at each iterative step matrix methods update all the pixel values at once, whereas series expansion methods update the pixel values sequentially. Matrix methods usually exhibit better convergence rates than series expansion methods since they can take into account higher order effects in a more straightforward manner. Techniques from either of these categories appear to be well suited for use in a CST system. A Closer Look This section details how specific CT techniques are adapted for use in CST. The first technique to be considered is backprojection. Generally the initial step in a reconstruction algorithm is a backprojection operation which yields a very smooth representation of the true image. Much of the fundamental differences between traditional CT and CST are evident in the way the projection data are handled in order to perform a basic backprojection operaÂ¬ tion. The second technique considered is Monte-Carlo backprojection (an original concept by this author and Dr. E. E. Carroll at University of Florida). The result of Monte-Carlo backprojection is a higher contrast image than obtained with traditional backprojection. 59 The third part of this section is a general discussion of matrix inversion techniques. The purpose of this discussion is to introduce the basic iterative nature of any practical solution strategy imposed by the nonlinearity of CST. The fourth part of this section adapts three series expansion techniques (ART, SIRT, and ILS) for use in CST. Backprojection Consider a projection formed by the conventional CT scanner shown in Figure 4.1. The pencil beam source irraÂ¬ diates the object and the transmitted radiation is measured at points in the image plane. The source and detector are moved in the direction indicated to scan the entire object and form a single projection. Other projections are formed by rotating the source and detector around the object and taking another scan. The number of photons reaching the detector at each point in a projection is a function of the average optical path through the object. The functional relationship is Ijc = IQexp{-/ p(s)vit(s)ds} (4.3) Â£jc where I^c is the output of the detector at position c in projection j, IQ is the number of photons leaving the source, p(s) is the density at point s, Pt(s) is the total mass attenuation coefficient of the material at s, and indicates that the integral is to be carried out 60 Figure 4.1: Example of a conventional transmission CT scanner. The source and detector move parallel in the direction indicated by the straight arrows, and also rotate about the object being imaged. 61 along the line between the source and detector when at position c in projection j. In practice one does not work with I_:_, but instead with R â– J given by K Rjc = -ln(Ijc/I0) = yt l PkAsjck (4.4) k-1 where pk is the density in pixel k, As^ck is the length of line il^c contained in pixel k, and K is the total number of pixels. The value of ASjck is nonzero only if line &jc crosses pixel k and so the sum in equation (4.4) is equivalent to the integral in equation (4.3). Since the liquid and vapor are of the same material, the mass attenuation coefficient is a constant and is taken outside the integral. The average density, p^â€ž, along line SLj.â€ž can be computed from equation (4.4) as Jc = R DC K /{J-i-f- I Asjck} k=l ^jc^kjc (4.5) A backprojection image is made by forming a composite of the average line densities, p.;â€ž, J ^ 62 J C ( j ) bPk =i l l j=l c=l PjcAsjck }/ J I C( j ) l As_i j=l c=l jck (4.6) where bpk is the density of pixel k in the reconstruction, and ASjc]< is used as a weighting factor. An analogous relationship for computing a backpro- jection image in CST can also be developed. The quantity corresponding to R^_ in CST is the singly scattered photon flux values, S^, given by K sÂ°j MPkAjkPSFjckAs JC â€” -L (4.7) where u SÂ° Pk As Ajk PSF^ck is uieuractiun ui me pnuuona arrxvxng at detector j from pixel k which have an energy in the range of channel c (Ec to Note that PSFjck is zero if it is not physically possible for a singly scattered photon from pixel k to arrive at detector j with an energy in the range Ec to Ec+p So, the sum of equation (4.8) represents an integral over the isogonic segment c of detector or projection j and is equivalent to the sum of equation (4.5). Another point to is the incoherent scattering coefficient, is the uninteracted source flux, is the density in pixel k, is the length of one side of a (square) pixel, is an average quantity which takes into account all aspects of photon transport from the source to pixel k, and from pixel k to detector j, and +Â» U -P v" o /â– > 4Â» i y> y* -P 4Â» K tr\ 4Â» j*Â» v\ e~* 63 note about equation (4.7) is the that double subscript jc can be replaced by an equivalent single subscript g, sg = sÂ°J PkAgkPSFgkis (4'7 }C â€” J. which is convenient in formulating the problem in matrix notation as in equation (4.1). The average density in an isogonic segment, Pjc, is computed as pjc = K Dc /ASSÂ° 1 iAjkPSFjck (4.8) and similarly the backprojection image densities, bp^, for CST are computed as bPk = J C ( j ) I I P 1=1 c=l jcPSFjck J c { j ) l l PSF j=l c=l jck (4.9) A basic problem with equation (4.8) is that the density distribution must be known in order to calculate the Ajj. factors. A reasonable solution is to use a homogeneous distribution with an average void fraction that is either assumed or gained from an independent measurement. The average void fraction of a CST backprojection image must 64 then be checked for consistency with the initial void fracÂ¬ tion. Any inconsistency can be resolved by some iterative refinement. Equations (4.5) and (4.6), in contrast, require no knowledge of the average void fraction and a one-step backprojection computation is possible in CT. Monte-Carlo Backprojection The motivation for using a statistical approach to form a backprojection image comes from three observations. First, the singly scattered photon flux, Sjc, is directly proportional to the average density in the isogonic segment jc. Second, any detected photons are due almost entirely to scattering in the liquid phase since the density ratio of liquid-to-vapor phase for water is usually 1000 to 1 in nuclear reactors. Third, the features of one detected spectrum can be cross correlated to the features of another spectrum detected at another location. These observations are summarized by the following example. If a single peak is observed in channel c of detector j and a similar peak is observed in channel c' of detector j', then the liquid phase is localized to the region where isogonic segment jc intersects segment j'c'. Monte-Carlo backprojection treats each detected spectrum as an independent one dimensional probability density function of where the liquid phase is located. A reconstruction of the liquid phase distribution is made by statistically cross correlating the detector spectra. The difference between Monte-Carlo backprojection and 65 traditional backprojection is that traditional backprojec- tion ignores the cross correlation information. The initial step in Monte-Carlo backprojection is to form an independent probability density function (pdf) from the modified average isogonic segment densities of each detector pdf j (c) Dc C ( j ) I P c=l DC (4.10) The average isogonic segment density is modified in the following manner: , K SÂ±/ AsSÂ° l A4 J k=l jk (4.8*) in order to account for the non-uniform spatial width of the isogonic segments. Assume that a fine-mesh square pixel grid is employed and that a pixel is either liquid filled or vapor filled. Provided a fine enough mesh is used, most density distribuÂ¬ tions can be reasonably represented by such a scheme. The number of pixels, N, which are liquid filled is easily computed from an estimate of the average void fraction, a, N = round((l - a)*K) (4.11) 66 where K is the total number of pixels. All the pixels are assumed vapor filled to start and the determination of which N pixels to designate as fluid filled proceeds as follows: 1.Given a set of J random numbers {r } , an associated set of channel numbers {ch-} is generated by inverting the pdfj's chi r^ Ã¡ Â£J pdf.; (c) ; j = l, 2, J (4.12) J c=l J There are J detectors and one channel from each detected spectrum is chosen by equation (4.12). 2.Associated with each channel chosen is a set of pixels, Set(ch'), which contains all the pixels within the corresponding isogonic segment. A new set is formed by taking the intersection of all J of these Set(ch-:)'s. The intersection set, I-Set, contains a total of M pixels. Taking a new random number, rJ+^, the mth element m = integer(M*rJ+1) + 1 (4.14) of I-Set is chosen as the pixel to be fluid filled. If M is zero, then no pixel is designated as being fluid filled. 3.Steps 1 and 2 are repeated sequentially until the necessary N pixels are turned on. Figure 4.2 compares the results of Monte-Carlo and traditional backprojection. The Monte-Carlo reconstruction is not as blurred and better represents the actual flow regime. Although the reconstructed image resembles the correct flow regime, a comparison of the calculated and detectedt singly scattered photon flux spectra shows that there is still room for improvement (Figure 4.3). 67 (b) (c) Figure 4. ^.Detector --Detector 2: Reconstruction of a 50% void flow regime. a) Monte-Carlo b) Simple backprojection; c) locations. fraction annular backprojection; source and detector 68 Figure 4.3: Comparison of the calculated singly scattered photon flux for the Monte-Carlo backprojection reconstruction of Figure 4.2 to the "measured" spectra, a) Upper detector; b) Lower detector. 69 The addition of a feedback mechanism to the pixel selection process is one possible improvement. Instead of randomly choosing a pixel from the intersection set, those pixels which improve the calculated spectra could be chosen preferentially. A more complete evaluation of the Monte- Carlo reconstruction method and possible variations is left for future work. Matrix Methods Generally matrix inversion methods are avoided in CT due to the large size (i.e. > 1 000 x 1000 ) of the matrices involved. Even in CST where there are only a few projections and, with a coarser pixel mesh, the matrices can be on the order of 200 x 200. Inverting matrices of this size requires large amounts of storage, long computing times, and is prone to round off error problems. Whatever the difficulties in implementing a matrix method solution, the compact notation is useful. There are many ways of formulating the reconstruction problem as a series of matrix operations. The most basic approach is to treat the problem as an optimization problem. A search is made for the density distribution, p*, which minimizes, either directly or indirectly, the difference between the measured data, S^, and the calculated singly scattered photon flux, F, that is, All data in this report, unless otherwise stated, are the result of Monte-Carlo photon transport calculations (i.e. numerical experiments) and not physical measurements. t 70 min |3^ - F| P or min |- [T(p)]p| P (4.15a) subject to the constraints Pvapor < Pk < Pfluid k=l, 2 / â€¢ â€¢ â€¢ t K. (4.15b) The notation |..| represents some norm functional. Any reconstruction algorithm in CST has a structure similar to that depicted in Figure 4.4. Note that the nonlinearity of the problem requires the recalculation of the transport matrix for each new density distribution, pm. Changes in the density distribution directly affect the attenuation of the source photons and scattered photons which is accounted for in the transport matrix. If the variation in pm is small so that [T] is approximately constant, then it may be possible to skip the recalculation step for a few iterations. In contrast, the corresponding matrix in CT is only a function of the system geometry and remains constant throughout the calculation. Another point at which the nonlinearity of CST needs to be addressed is in the calculation of Spm+^. One method of calculating Aprn+'*- is to use the method of steepest descent 5pm+1 = (Jacobian of F(pm))~1â€¢(S1 - F(pm)) . (4.16) The transport matrix can be treated as a constant or as a function of the density in the calculation of the Jacobian. 71 Figure 4.4: Flow chart outlining the general structure of a reconstruction algorithm for Compton scatter tomography. The superscript m denotes the iteration index. 72 Treating [T] as independent of p means a smaller step size, f, is necessary and that Ap does not point directly toward the path of steepest descent. Including the dependence of [T] on p results in a significant increase in the storage requirements and in computing time. Algebraic Reconstruction Technique ART is historically the first of the series expansion techniques applied to commercial CT devices [31]. There are two basic forms of ART algorithms: multiplicative and addative. Presented here is a variant of the additive algorithm. A single ART iterative step involves correcting the density vector by adding to it a correction vector -m+! = -m cm+lâ€¢wm+1 (4.17) where cm+^ is a constant correction factor and wm+-*- is a weighting vector. The correction factor is computed so that one of the calculated singly scattered photon flux values is consistent with the corresponding measured value, that is, S 1 g T?nt+1 Fg I- yrrm+lrim+l 1~ Â¿Tgk Pk K (4.18) Subsequent iterations correct the density vector with respect to the other measured values, one at a time. The subscript g is incremented cyclically 73 gm+-*- = (gm + 1) modula G + 1 with each iteration in order to cover all G data points. The correction factor can be derived by substituting equation (4.17) into (4.18) It k m+lTTm+l gk wk The transport matrix elements are unknown until the new density vector is known which in turn requires the correction factor. One solution to this predicament is to use the transport matrix calculated during the previous iteration, that is, ,m+l _ el _ Yrpin _.m sg Â¿Tgk Pk I mm wm+l Tgkwk .m+1 (4.19*) This approach is valid whenever the transport matrix does not change much from one iteration to the next, as when close to the final solution. There are several possible choices of weighting vectors. The simplest choice is an identity vector, however, this results in all the pixel values being altered including pixels which do not affect the sum in equation (4.18). An identity vector modified such that 74 g = g m+1 (4.20a) sensibly corrects only those pixels contained in the isogonic segment corresponding to SÂ¿. Another reasonable y weighting vector is g = g m+1 (4.20b) which effectively weights a pixel according to the probaÂ¬ bility that a scattered photon from that pixel is counted in Sg. A third possible weighting vector is m+1 (4.20c) g=g which is analogous to the fractional ray length weighting scheme of CT [32]. The weighting scheme chosen here and elsewhere in this report is that of equation (4.20b). Figure 4.5 graphically shows how an ART solution proceeds for a two data point, two pixel situation. The lines L^1 and Lâ„¢ correspond to the two simultaneous equations given by (4.21) The lines change as the density vector is corrected since the transport matrix is a function of the current density distribution iterate. The lines tend to stabilize as the 1.0 p2 Figure 4. 0.2 U1 J 0.2 _J 0.4 0.6 0.8 1.0 Pi Graphic demonstration of how an ART solution proceeds for a two data point, {L^}, two pixel, {pi}, solution. Superscripts identify the iteration number. 76 final solution is approached, indicating that the transport matrix is only varying slightly. Simultaneous Iterative Reconstruction Technique SIRT is similar to ART except that the density vector is corrected using all the projection data, Sâ€™*-, simultaÂ¬ neously. SIRT attempts to combine a complete cycle of ART iterations (i.e. G iterations) in calculating the correction vector. The correction vector elements are given by m+l _ Pk i 9 tVal,((sg ' f5,/( i <â€™VÂ°Â§> (4.22) where at is the measurement variance. The correction vector y resembles a backprojection image of the error vector em = s1 - Fm (4.23) which can be seen by comparing equations (4.22) and (4.9). Reference 30 indicates that some method of damping may be necessary to avoid a solution which oscillates and fails to converge. Once the correction vector is obtained it is multiplied by a constant damping factor, f. The damping factor is calculated so that the variance weighted euclidean norm of the error vector is minimized min | S1 - [Tm](pm + fm+1Apm+1)|. (4.24) f 77 The solution to equation (4.24) is (4.25) g k The proposed SIRT algorithm for CST is where the correction vector is given by equation (4.22) and the damping factor by equation (4.25). Iterative Least-Squares Technique The ILS technique computes a correction vector which minimizes the variance weighted euclidean norm of the error vector. This can be expressed as (4.27) The correction vector which minimizes U is found by setting the derivative of U with respect to the correction vector equal to zero. This leads to the set of equations k'=l K (4.28) Solving the system of equations of (4.28) would violate the principal motivation of going to an iterative approach, 78 which is introduced to avoid the computational overhead and pitfalls of matrix methods. The usual course of action is to assume that all the the correction vector elements are zero except the element corresponding to the current pixel k\ Equation (4.28) becomes k '=1 (4.28*) K Equation (4.28*) overestimates the magnitude of the indiviÂ¬ dual elements of the correction vector which necessitates the introduction of a damping factor as is done for SIRT. The form of the ILS damping factor, which is calculated according to equation (4.24) also, is IÃ(Â£g/ag) 9 iKÃT^Ap^h/OgÃ2 The proposed ILS algorithm for CST is Pk+1 = raid ÃPvapor' Pk + fm+1Apg+1, Pfiuid^ (4'30) where the correction vector is given by equation (4.28*) and the damping factor by equation (4.29). 79 Comparison of Reconstruction Techniques The ART, SIRT, and ILS reconstruction techniques are compared here in order to select one of them to use in a detailed evaluation of CST. Each of the techniques are used to reconstruct the two flow regimes shown in Figure 4.6. The flow regimes are analogous to the "hot spot" and "cold spot" phantoms used in CT studies. The criterion used to judge the best techniquÃ© is the euclidean distance between the true density vector and the reconstructed density vector. The rate of convergence is considered as a secondary criterion. The projection data (spectra) used for each technique are generated by Monte-Carlo transport calculations. The singly scattered photon spectra calculated for the two detectors are shown in Figures 4.7 and 4.8. Table 4.1 lists the pertinent "experimental" conditions. The reconstructions of the hot-spot and cold-spot flow regimes are shown in Figures 4.9 and 4.10, respectively. All the reconstructions exhibit undesirable streaking artifacts. The artifacts are due to the lack of edge inforÂ¬ mation from using only two views. The reconstructions of the hot-spot are better than those of the cold spot. IntuiÂ¬ tively this is expected by observing that the signal-to- noise ratio is higher for the hot-spot spectra. The ART results have more pronounced streaking artifacts than either the ILS or SIRT results. The poorer performance of ART is in keeping with similar CT phantom study tests done for (dimensions in centimeters) x 00 o Figure 4.6: Source and detector configuration used for the hot-spot and cold spot tests. Energy [keV] (a) (b) oo The measured singly scattered photon spectra for the hot-spot test, a) Spectrum from the upper detector; b) Spectrum from the lower detector. Figure 4.7: (a) (b) Figure 4.8: The measured singly scattered photon spectra for the cold-spot test. The ordinate scale is the same as for Figure 4.7. a) Spectrum from the upper detector; b) Spectrum from the lower detector. 83 TABLE 4.1 Data pertaining to the hot-spot and cold-spot experiment. Detectors: Number of detectors Detector type Diameter Active thickness MCA channel width 2 planar 1 cm 0.1 cm 1 keV Source: Isotope Cs-137 Physical shape cylinder Diameter 0.31 cm Height 0.63 cm Material Density: Hot-spot bubble density 1.155 g/cc bulk density 0.0 g/cc Cold-spot bubble density 0.0 g/cc bulk density 1.155 g/cc SIRT Figure 4 V:Â«'f ART r'*'- â– s~. â– â– f .9: Reconstructions of the hot-spot flow pattern by the ILS, SIRT, and ART algorithms. 85 Figure 4.10: Reconstructions of the cold-spot flow pattern by the ILS, SIRT, and ART algorithms. 86 medical applications. There is no discernible difference between the ILS or SIRT results. The rate of convergence and accuracy of the reconstrucÂ¬ tion methods is shown in Figure 4.11. The discrepancy function is defined as the euclidean distance between the true flow regime and the reconstructed flow regime. The spectra difference function is defined as the euclidean distance between the measured spectra and the spectra calcuÂ¬ lated using the reconstructed flow regime. Both the discreÂ¬ pancy and spectra difference functions are normalized by their values for the initial density distribution, which is an all vapor condition. The convergence of all the reconÂ¬ struction algorithms for the cold spot test is similar. The ILS and SIRT methods show better convergence properties than the ART method for the hot spot test as judged by the discrepancy function. The ILS and SIRT method are more desirable than the ART method due to their consistent performance. There is no significant difference in the performance of the ILS or SIRT methods based on examining the discrepancy or spectra difference functions. The choice of whether to use either the ILS method or the SIRT method appears to be one of personal preference. This chapter successfully demonstrates that CT reconÂ¬ struction algorithms can be used with CST. As vital as the reconstruction algorithms are, they are not sufficient for implementing a CST system. The next chapter considers the effect of noise in the spectra on the reconstruction process Discrepancy o ILS â–¡ SIRT A ART Comparison of the accuracy and convergence of the ILS, SIRT, and ART algorithms. a) Hot-spot flow regime; b) Cold-spot flow regime. oo Figure 4.11: Figure 4.11â€”continued. 89 and looks at some of the other practical aspects of designing a complete CST system. The ILS algorithm is chosen as the reconstruction algorithm in the discussions of the next chapter and the remainder of the dissertation. CHAPTER 5 COMPTON SCATTER TOMOGRAPHY Three major issues of interest in designing a CST system or in evaluating the suitability of CST for a particular purpose are accuracy, image resolution, and image contrast. The nonlinearity of CST makes it difficult to quantify any of these for all values of void fraction, flow regime, position in space, etc. A complete analysis requires varying all pertinent parameters combinatorially. The following discussions do not attempt a complete quantiÂ¬ tative analysis, but seek to identify general trends and limitations of CST. An Example In order to illustrate how the non-linear aspects of CST affect the flow image reconstruction, consider the case of a small error in the source strength specification. The detector response is directly proportional to the source strength [see equation (2.6*)]. If the correct flow distriÂ¬ bution along with too low a value for the source strength is used in calculating the detector response [see equation (2.12)] then the calculated spectrum will be uniformly lower than the measured spectrum (Figure 5.1). The reconstruction algorithm is only able to compensate for the error in SÂ° by 90 Counts Correct Figure 5.1: Result of specifying too low a source strength to the reconstruction program. The calculated spectrum is based on the correct flow distribution, but with SÂ° 10% too low. 92 adjusting the fluid density in the individual pixels. Examining equations (2.1) and (2.6*) it is seen that the density enters the detector response equation as a linear factor and as an argument of the exponential attenuation factors. If the attenuation factors did not exist then the low value of SÂ° would result in a uniform increase in the density throughout the flow image thereby uniformly increasing the calculated spectrum to match the measured spectrum. However, the attenuation factors are present and as the density increases the attenuation increases (i.e. the magnitudes of the attenuation factors decrease) requiring more of an increase in the density than if the attenuation factors are not present. Since attenuation is not uniform over the flow cross section this leads to a non-uniform change in the density or, equivalently, a change in the density distribution. The constraint conditions placed on the minimum and maximum allowable density values [see equaÂ¬ tion (4.30)] can exacerbate the error in the density distriÂ¬ bution if some of the pixels are constrained from changing while others are not. Specifying a source strength value in the input to the reconstruction algorithm which is lower than actually used in making the measurements results in the reconstructed density values being too high and the density distribution being incorrect (Figure 5.2). Accuracy The topic of accuracy is confined here to examining the extent to which errors in the input data result in errors in 93 (a) (b) y x Figure 5.2 Reconstruction of an annular flow using an incorÂ¬ rect source strength value. a) SÂ° 10% too high; b) SÂ° correct; c) SÂ° 10% too low. 94 the reconstructed density distribution. Such an approach views error propagation as a concern in the operation of an imaging system rather than in its design. The design process is limited to minimizing the errors in the measured data due to the lack of analytic tools to handle nonlinearÂ¬ ities. In recognizing this limitation the approach taken here is to look at a specific system and its response to the expected types of error. The specific system chosen is designed to measure the two-phase distribution of water in a pipe. The system consists of two detectors and a source placed as shown in Figure 5.3. No actual experiments are done using this setup, instead the detector responses are generated by a 3-D Monte-Carlo photon transport computer code. The various components of the setup are modeled in the transport calcuÂ¬ lation in the following manner. The source is cesium-137 which has a single gamma ray emission energy of 662 KeV. The shape of the source is a cylinder 0.64 cm long and 0.32 cm in diameter. A perfectly absorbing shield collimates the source into a wedge beam which illuminates a thin cross section of the pipe. The two detectors are 1 cm in diameter and 0.1 cm thick. These high energy resolution detectors sense the scattered photon flux. There are a total of 132 channels each 1 keV wide in each spectrum and due to scattering angle blur it is estimated that this represents 126 independent channels. y Figure 5.3: The source-pipe-detector configuration used in evaluating the effects of flow regime, void fraction, and noise on CST. 96 The presence of the pipe wall is neglected in the Monte-Carlo calculations and in reconstructing the flow cross sections. In a real experiment the pipe would be a fixed scattering (and absorbing) material whose dimensions, location, and density are known and can be accounted for in the reconstruction process. Ignoring the pipe reduces the computer runtime and programming effort considerably. The most obvious source of error is the statistical noise in the spectra due to using radiation in the imaging process. Systematic errors may also be present in the spectra due to poor instrument calibration, dynamic bias [33], or spectra processing errors. The only other data required by the reconstruction algorithm are the physical dimensions of the system. Since there are no moving parts in the system, these dimensions can be measured quite accuÂ¬ rately and do not represent a source of significant error. The system as defined leaves only three operational parameters which can vary-- data acquisition time, flow regime, and void fraction. Flow regime and void fraction are the two quantities being measured while data acquisition time is the only parameter under control of the operator. The error in a spectrum is closely related to the data acquisition time. The shorter the data acquisition time the fewer the photons detected and the higher the statistical noise in the spectrum. The longer the data acquisition time the more the stochastic flow variations introduce dynamic bias and blur the spectra. The data acquisition time is 97 chosen to balance the statistical and systematic errors in order to minimize the resulting error in the reconstructed density distribution. In the following analysis the error in the reconÂ¬ structed density distribution due to errors in the spectra under different flow regime and void fraction conditions is considered. Table 5.1 lists the 3 parameters being varied and their range of values. The four flow regime models are illustrated in Figure 5.4 as they should appear when correctly reconstructed. These four regimes cover a broad range of two-phase flows and most regimes can be viewed as combinations or extreme cases of these basic types. Two particular cases of the stratified regime are considered; stratification symmetric and stratification asymmetric with respect to the two detectors. These two cases are chosen in order to judge the sensitivity of CST to asymmetries in the flow. Three separate void fraction cases (i.e. 0.3, 0.5, and 0.7) are looked at for each of the five flow regimes. Figure 5.5 is an example of the detected scattered photon spectra for the case of annular flow and void fraction of 0.5. Similar plots for all fifteen flow regime- void fraction cases are given in Appendix V. Figure 5.6a is the reconstructed image of the flow cross section using the spectra in Figure 5.5. A 15 x 15 pixel array is used and the number of dots is linearly proportional to the density. It is known a priori that 6 pixels in each corner are voided which leaves 201 unknown pixel densities to calculate. The 98 TABLE 5.1 The three parameters propagation analysis, being examined as part of the error Parameter Range of values Flow regime Annular Inverted annular Bubbly Stratified Void fraction 0.3 0.5 0.7 Spectrum noise No added noise 10% additive random noise and 25% magnitude decrease 30% additive random noise 99 Figure 5.4: The five flow regime cases at 50% void fraction as they should appear when correctly reconÂ¬ structed. a) Annular flow; b) Inverted annular flow; c) Bubbly flow; d) Symmetric stratified flow; d) Asymmetric stratified flow. 100 Energy [keV] Figure 5.5: Example of the singly scattered photon spectra used as input to the reconstruction program. This particular spectrum is from a 50% void fraction annular flow. 101 Figure 5.6: Reconstruction of the annular flow regime using the spectrum in Figure 5.5 with the three noise conditions. a) No added noise; b) 10% additive noise and 25% magnitude decrease; c) 30% addiÂ¬ tive noise and no magnitude decrease. 102 ratio of unknowns (pixels) to knowns (spectrum channels) is approximately 0.80. The ILS algorithm is employed in all the reconstrucÂ¬ tions in Figure 5.6 and this chapter. The reconstruction algorithm is an integral part of a CST imaging system and the sensitivity of CST to the experimental parameters is somewhat dependent on the algorithm chosen. In the most general system in which no assumptions are made about the flow regime the ILS algorithm should perform as well or better than any other algorithm. Three different noise or error cases are considered. The first case is the detected spectra with only the inherent noise from the Monte-Carlo calculations (less than 6%). The second case adds a normally distributed random noise component to each measured spectrum and introduces a systematic error by decreasing the magnitude by 25% (i.e. source magnitude error). The standard deviation of the noise component is 10% of the measured spectrum's magnitude. Similarly, the third noise case adds a 30% normal random noise component, but introduces no systematic error. The normally distributed random noise is generated by the method of adding five uniform random deviates, correcting the values to have zero mean, and adjusting the standard deviaÂ¬ tion to the appropriate value for the particular noise level. Figure 5.7 illustrates the results of introducing the second and third noise cases to the measured spectra for annular flow. Counts [relative] Figure 5 .7: Results of introducing the noise and systematic error to the spectrum in Figure 5.5. a) 10% noise and 25% magnitude decrease; b) 30% noise only. â€”1 400 103 104 Table 5.2 summarizes the results of reconstructing the flow cross section density distribution for the various combinations of flow regime, void fraction, and noise content. The reconstructions are rated according to two criteria: void fraction and image quality. The cross sectional average void fraction is computed using those pixels in the image which fall within the pipe boundary. The image quality is a subjective measure of how well the reconstructed image matches the true flow distribution. The error in the measured void fractions is consisÂ¬ tently larger for the low void fraction cases. The detector response is less sensitive to density variations at low void fractions than at high void fractions. Therefore, a small error in the detector response results in a larger error in the density at low void fractions than at high void fractions. Another observation from Table 5.2 is that CST is more sensitive to systematic error in the spectra than to random error. The robustness of the reconstruction process to random error is due to the fact that neighboring channels in a spectrum are not independent with respect to the indiviÂ¬ dual pixel densities. Due to the finite size of the source and detector even an abrupt spatial density change results in smooth spectral features. The best this reconstruction algorithm can do with a noisy spectrum is to find a density distribution that results in a spectrum which fits the average trends and not the individual channels. Systematic 105 TABLE 5.2 Results of reconstruction tests. Flow Reqime Actual Void Fraction Noise Case 1* Noise Case 2 Noise Case 3 Void Frac. Imaqe Void Frac. Imaqe Void Frac. Imaqe 0.30 0.22 9 0.68 0 0.26 9 Bubbly 0.50 0.45 9 0.75 0 0.47 9 0.70 0.71 0 0.85 0 0.72 0 0.30 0.24 9 0.70 0 0.28 9 Strat. 0.50 0.45 9 0.76 0 0.49 9 (sym.) 0.70 0.68 9 0.84 0 0.69 9 0.30 0.34 9 0.72 0 0.34 9 Strat. 0.50 0.48 9 0.77 0 0.52 9 (asym. 0.70 0.68 9 0.84 0 0.70 9 0.30 0.26 9 0.68 0 0.30 9 Annular 0.50 0.46 9 0.73 0 0.47 9 0.70 0.66 9 0.82 0 0.66 9 0.30 0.27 9 0.73 0 0.32 9 Invert 0.50 0.47 9 0.78 0 0.50 9 Annular 0.70 0.67 9 0.84 0 0.69 9 Key to table â€¢ Compares well with the actual distribution. 9 Matches the actual distribution except in a few places. 0 Flow is identifiable but densities are incorrect. 0 Flow is unidentifiable. * There is some inherent noise in the Monte-Carlo calculation, but this is less than 6% in all cases. 106 error, unlike random error, affects a spectrum in a more uniform or smoothly varying manner and does not "average out" in the reconstruction process. Since systematic error has a greater impact on the reconstructed density distribuÂ¬ tion it should be minimized. Resolution The only available information in CST about the density distribution comes from the scattered photon spectra. In order for an object to be resolved it must cause a statistiÂ¬ cally significant change in the spectra by its presence. This provides a convenient working definition of resolution and a useful perspective for optimizing system parameters with respect to resolution. The following discussion identifies those parameters which affect resolution and indicates in what manner a parameter can be altered to improve resolution. Table 5.3 lists the parameters which affect image resoÂ¬ lution and their ideal values. The ideal values are classified according to whether the value increases the spatial information content in the spectra or increases the number of photons detected thereby reducing the statistical noise. There exist practical considerations not listed in the table which force the parameter values in a real system to deviate from the ideal. A single detector is only able to give explicit spatial information in the direction normal to its characteristic isogonic lines. The energy resolution of the detector; the 107 TABLE 5.3 Some of the parameters which affect image resolution. Parameter Ideal value Detector size point (S) efficiency 100% (P) energy resolution <<0.001 (S) number of detector infinite (S)(P) detector-to- object distance (S)(P) Source size point (S) collimator perfect absorber (S) energy 300 keV to 4 MeV (P) source-to- object distance as close as possible (S)(P) Container (pipe) thickness as thin as possible (P) . density , < density of fluid (P) size (diameter) (P) S - parameter choice affects the spatial information. P - parameter choice affects the number of photons detected. 108 size of the detector and source; the source energy; and the relative locations of the pipe, source, and detector all combine to define an effective spatial isogonic line frequency (ESILF). The ESILF is the effective number of isogonic lines per centimeter measured in the direction normal to the isogonic lines. This number is not a constant, but varies over an image. Consider the case of a point detector and point source where the detector has a 1 KeV energy resolution. The resolution of the detector dicÂ¬ tates an isogonic line spacing no closer than 1 KeV. The 1 KeV energy spacing translates to an actual spatial spacing determined by the scattering geometry and source energy. If the point detector is replaced by a detector 1 cm in diameter then the upper bound on the ESILF is 1 per centimeter regardless of the detector resolution. This is due to the scattering angle blur. In general, the size of the source and detector determine the ESILF in the region close to the source and the detector energy resolution determines the ESILF in the regions further away. Ideally the smallest possible source extension and detector are used in order to enhance the ESILF. There are no drawbacks to making a source as small and compact as possible. The limiting factors are the specific activity of the source material and the source activity required to meet the statistical noise and framing rate criteria. A typical 10 Ci cesium-137 source in the shape of a ball can be made with a diameter of ~6 millimeters. 109 There are definite drawbacks to reducing the size of the detector, however. Making a detector small adversely affects its intrinsic efficiency. Also, the smaller a detector the fewer photons it counts due to field-of-view effects. The size of the detector and source requires balancing the competing needs of spatial resolution, time resolution, and the uncertainty in the density values. The choice of what type of detector to use is very restricted. Since all the spatial information of the image is encoded in the photon energies, detectors with the highest possible energy resolution are required. Liquid nitrogen cooled, semiconductor detectors are the only viable choice with resolutions on the order of 0.002 (i.e. the full-width-at-half-maximum per photon energy). Another detector feature to consider is the peak-to-Compton ratio. Counts in the Compton continuum are unusable due to the multiply scattered photons in the spectrum. In addition the tCompton continuum of the high energy photons may interfere with the low energy end of the singly scattered photon spectrum. A high peak-to-Compton ratio helps to avoid these problems. Image Contrast The discussion on resolution has concentrated on what might be called the spatial sampling rate. An object smaller than the spatial sampling rate may enter the image plane and still cause a statistically significant change in 110 the spectra. The object, however, is blurred in the reconÂ¬ structed image. It is also possible for an object to be larger than the spatial sampling rate and not cause a statistically significant change in the spectra. This can occur for two reasons. First, the object could have the same electron density as the rest of the surrounding material and be indistinguishable from it. Second, if the object is sufficiently shielded by the bulk material, then the detectors never see the scattered photons coming from the region of the object. The relationship between the size and density of an object which can be sensed is analogous to the relationship between sharpness and contrast in photography. The two factors affecting whether an object can be seen are its size and the difference between its electron density and that of the bulk material. Specifying the minimum size and density disturbance that can be sensed necessitates estimating what constitutes a statistically significant change in a spectrum. The two major sources of error which define statistical significance are measurement-data processing error and reconstruction error. The measurement- data processing error has three main components: counting statistics, background removal, and detector efficiency correction. Equation (2.8) gives an expression from which the error in these three components can be related to the error in the singly scattered photon spectrum. The relative Ill error, xcr in the single scatter counts for a particular channel, c, in a spectrum is given by J a - -Ã -t 'D 1 o + Si,2 a2 1/2 ] / (5.1) where aÂ£ is the standard deviation of the detector efficiency correction, aD is the standard deviation of the detected count in channel c, and aB is the standard deviation of the background count in channel c. The variance of the total detected counts, Dc, and the background counts, Bc, can be approximated by = Dc and aB = Bc . (5.2) The efficiency correction, e, does not represent a signifÂ¬ icant source of error. This is because the value of the efficiency correction is not as important as its functional form. Any error in the value of the efficiency correction can be compensated for by a calibration procedure, but only as long as the correction in one channel relative to the others is done correctly. The functional variation of detector efficiency with energy is well understood in the range of energies of interest in Compton scattering. The relative error in the efficiency correction is probably on 112 the order of 1% i.e. aÂ£/e =0.01. (5.3) The reconstruction error involves the error in the cross section library, the usual error in finite precision arithmetic, and the error associated in representing the two-phase flow with discrete pixels. Hubbell [34] estimates the error in the National Bureau of Standards published cross section library for the energy range of interest (200 keV to 4 MeV) to be on the order of 1%. The reconstruction algorithms developed here are not very sensitive to round off error, etc. and with the advent of microcomputers which use 64 bit real arithmetic means that this error can be neglected. The primary source of error in the reconstrucÂ¬ tion process is in the discrete representation of the flow. Experience has shown that most of this error is in calculaÂ¬ ting the attenuation factors. Using small pixels (~0.1 mean-free-path) and averaging the attenuation factors from several random points per pixel can reduce this error to the order of 5%. This error mainly affects those pixels on the edge of features in the flow. The overall result of the reconstruction error is that even when the correct density distribution (in discrete form) is used the difference between the calculated spectra and the measured spectra is about 5% aF/sl =0.05 . i. e. (5.4) 113 The change in the detected number of photons, ADC, which an object must cause in order to be seen is given approximately by ADc/ec > {a g + Qp}1/2 (5.5) 2 where ap is the variance of the reconstruction error. The value of ADC for a particular object can be estimated using a modified form of equation (2.1) ADC = SÂ°ApywAzAiA0ecPSFc (5.6) where Ap is the difference between the object density and the bulk material density and Az is the mean distance through the object as seen by the source photons. The attenuation function, A^, and the isogonic segment characÂ¬ teristic function, PSFC, are both functions of the object size, shape, and orientation. Even assuming a simple object shape, PSFC is not easily expressed as an analytic function of Az and, in general, PSFC must be calculated numerically. In practice, equation (5.6) is computed for several sizes of an object and the value of ADC which satisfies equation (5.5) is used to interpolate the results to find ApAz. The simplicity of equations (5.1) through (5.6) belie the problem of specifying a minimum size object and density variation that can be sensed. The attenuation factors in equation (5.4) introduce spatial, void fraction, and flow 114 regime dependence to ApAz. The detector efficiency adds an energy dependence and the other variables introduce yet other functional dependencies. Note that the choice of source energy is coupled to the problem of image contrast by the mass attenuation and scatÂ¬ tering coefficients. The mass interaction coefficients decrease as the source energy increases. At higher source energies the photons penetrate better and shielding by the bulk material is less of a problem. At lower source energies the photons interact better and there is better sensitivity to a given change in the electron density. There exists, for each particular imaging system, a source energy which optimizes the detector response by balancing the effects of attenuation and interaction. An optimum source energy can be found which maximizes the number of detected photons coming from a point in the flow. That is, a search is made for the source energy which satisfies 3D(x,y)/3E = 0 . (5.7) In the case of the experimental setup in Figure 5.3 the optimum source energy is close to that of cesium-137 for 0% void fraction and a point in the flow furthest from the source and detector. The optimum energy decreases with increasing void fraction and/or for points closer to the source and detector. 115 The choice of source energy must also take into account the practical considerations of finding an isotope with a mono-energetic gamma of the right energy and with a reasonable half-life. A dual energy source such as cobalt- 60 could be used with some care to separate the two singly scattered photon spectra. Another practical problem is source handling and shielding. High energy, high activity sources are impractical due to the amount of shielding material required. Image contrast is also affected by the presence of a dense scattering and attenuating object such as a steel pipe. The pipe wall introduces a threshold effect. If the pipe is thick enough and its electron density much greater than the bulk material being imaged, then the scattered spectrum of the pipe wall can dominate the detector response. The statistical variation in the pipe wall spectral component sets a threshold of statistical signifiÂ¬ cance which variations in the bulk material must overcome in order to be sensed. The counts in the singly scattered photon spectrum due to the bulk material, S^-^, must be greater than the expected Poisson statistical variation in the spectrum, i-e- Sbulk > aS (5*8) where ag is given by equation (5.1). The threshold effect is important at high void fractions where there are few 116 photons interacting in the bulk fluid compared to the pipe. Increasing the data acquisition time in order to increase the number of detected photons from the bulk fluid is one way of reducing the impact of the threshold. Figure 5.8 shows the spectrum from a 50% void fraction annular flow (same as Figure 5.5) with a 0.6 centimeter thick iron pipe taken into account. The peaks at the upper and lower ends of the single scatter spectrum are dominated by the pipe wall and it is more difficult to see changes in the bulk material density in the isogonic segments associated with these two peaks. The pipe also attenuates the singly scattered photons from the bulk fluid especially toward the low energy end of the spectrum. A CST instrument is a highly integrated system in which each design parameter is chosen in concert with the others and with respect to the characteristics of the object being imaged. The design process can be understood as optimizing three basic system specifications: accuracy, image resoluÂ¬ tion, and contrast resolution. Accuracy determines how close a reconstructed image matches the gross features of the object being imaged. Image resolution and contrast resolution determine the minimum size features of an object which can accurately be reconstructed in the image. 117 Energy [keV] Figure 5.8: Detected spectrum from a 50% void fraction annular flow inside a 0.6 cm thick walled steel pipe. CHAPTER 6 SUMMARY AND CONCLUSIONS This study investigates and demonstrates several methods of inferring the flow parameters of two-phase flows using wide beam illumination and few (i.e. 2) detectors. The first step is the recognition of the singly scattered photon spectra as projections of the phase distribution. The spatial distribution of the fluid-vapor phase in a slice of the pipe is encoded with respect to energy in the singly scattered photon flux. Two basic methods are detailed for decoding the spatial information: the method of spectral moments and the method of computed tomography. Previous authors have recognized that tomography of ' two-phase flow is possible using Compton scattering of gamma rays. The majority of their efforts employ elaborate colli- mation schemes to reduce the likelihood of multiply scattered photons reaching the detectors. This report demonstrates using 3-D Monte-Carlo photon transport calculaÂ¬ tions that the multiply scattered photon spectrum for a wide beam source and uncollimated detector setup lacks any structure which would interfere with extracting the singly scattered photon spectrum. A simple background subtraction procedure appears sufficient for removing the multiply scattered photon spectrum component. Not only is the 118 119 detector collimator unnecessary it is undesirable since it can reduce the sensing of singly scattered photons and thereby degrade system performance. Summary Examination of the low-order moments of the singly scattered photon spectra for a one source, two detector configuration provides sufficient information for flow regime identification. The two detectors are placed in symmetric positions with respect to the source and pipe. Flow asymmetries are revealed by differences in their spectral moments. Further identification of the actual symmetric or asymmetric flow pattern is made on the basis of the values of the first and second moments of the spectra. Quantitative measures of the flow may be possible using this method and a spectral calibration technique. The real'focus of this report is the adaptation and successful demonstration of computed tomography techniques with Compton scattering. Three series expansion techniques are chosen for adaptation because of their iterative natureâ€” ART, SIRT, and ILS. The adaptation process involves deriving the equations for each algorithm from the forward Compton scattering relations and altering the structure of the iterations to include a recalculation of the transport matrix. Recalculation of the transport matrix overcomes the nonlinearity associated with attenuation without having to include it explicitly in any derivatives. Application of the modified ART, SIRT, and ILS 120 algorithms to a hot-spot and a cold-spot reconstruction problem indicates that SIRT and ILS are more accurate than ART. A more extensive testing of the ILS algorithm for a variety of model flow regimes demonstrates the potential of CST as a quantitative measurement technique. Monte-Carlo backprojection is another reconstruction technique looked at in this report. This is a new reconÂ¬ struction method which operates by statistically cross correlating the spectra of a few detectors to produce a tomographic image. This technique appears well suited to situations where the spectra contain a lot of structure. One example would be the binary type flow patterns expected for high frame rate imaging of two-phase flows. Another application would be to identify and localize flow anomalies by looking at the difference between a time averaged spectrum and an instantaneous spectrum. The difference spectrum would show structure due to any flow anomalies. Based on the observations and results of Chapter 5 a CST system is best designed with specific measurement goals in mind. The design process involves striking a balance between image resolution, image accuracy, framing rate, and sensitivity. The balancing or optimization is such that there'is a tight coupling between a CST system and a particular application. Even with careful optimization it may not be possible to achieve set design goals due to practical limitations. 121 A Practical Application As a practical illustration and summary of the design concerns discussed in Chapter 5 consider the problem of imaging the phase distribution of a steam-water flow inside a 4 inch schedule 40 steel pipe. The source-pipe-detector configuration is assumed to be the same as shown in Figure 5.3. Figures 6.1 and 6.2 show the results of a Monte-Carlo calculation of the scattered photon spectrum for this setup when the void fraction is zero. The two quantities of interest to an experimenter are the source activity required to achieve a set framing rate and the minimum size bubble that can be sensed. Estimating the required source strength necessitates first calculating the minimum number of photons which must be counted in a single channel in order to overcome statistical noise in the measurements. The source strength is given by dividing the minimum number of photons for a single channel, Dc, by the number of photons detected in channel c per emitted source photon and multiplying by the framing rate. A relation for Dc can be obtained by substituting equations (5.1) and (5.2) into equation (5.8) rD(SNR)2(rD + rfi) Dc = 5 5â€” e(l - r2(ae/e)2) where rD is the ratio Dc/SÂ¿, rB is the ratio Bc/SÂ¿, rc is the ratio SÂ¿/SÂ¿, (6.1) Photons detected per source photon [photons/Ci-sec Figure 6.1: Results of a Monte-Carlo calculation of the scattered photon spectrum for a 4 inch schedule 40 pipe full of water. The calculation assumed a 100% efficient detector and a 1 keV MCA channel width setting. 122 123 Energy [keV] Figure 6.2: Detail of the enclosed region on Figure 6.1 showing the total number of photons detected, D , the multiply scattered photons, Bc, the singly c scattered photons, SÂ¿, and the singly scattered photons from the bulk fluid, SÂ¿. 124 SNR is the signal-to-noise ratio, and Dc, SÂ¿, Bc, and are defined pictorially in Figure 6.2. Equation (5.8) gives an inequality which must be satisfied in order to extract information about the flow from the photon spectra. Replacing the inequality by an equality is equivalent to specifying a working SNR of 1. Figure 6.3 illustrates some flow reconstructions obtained when the SNR is near 1. The ratios rD, rB, and rc for channel c estimated from Figure 6.2 are 2.79, 1.27, and 1.51, respectively. SubstiÂ¬ tuting these ratio values, equation (5.3), an assumed detector efficiency of 15%, and an assumed SNR of 2 into equation (6.1) indicates that 300 photons must be counted in channel c (i.e. Dc = 300) per measurement. The Monte-Carlo calculations indicate that 36 photons are detected in channel c per Curie of source activity per second for a 100% efficient detector. The required source activity for a 15% efficient detector is given by (300 photons/frame)(10 frames/sec) activity = = 556 Ci (0.15) (36 photons/Ci-sec) (6.2) The minimum size steam bubble that can be sensed with this CST system is found using equations (5.5) and (5.6). The RHS of equation (5.5) is evaluated using equations (5.1) and (5.4) along with the information from the source activity calculation which yields 125 Figure 6.3: Reconstruction of a 50% void fraction annular flow regime for three different SNR values of photon spectra. 126 ADC > 23 counts . (6.3) The same comment about the SNR with respect to equation (6.1) applies to this inequality and it is assumed that ADC needs to be twice as large as the RHS of equation (6.3) i.e. ADC = 46 or ADC/DC = 46/300 = 0.15 . (6.4) Figure 6.4 plots the results of a numerical evaluation of equation (5.6) for spherical bubbles located in the center of the 4 inch pipe. The peak in the curve arises because of the decreased attenuation of the scattered photons as they pass through the larger bubbles. The smaller bubbles do not significantly affect attenuation and simply result in less photons being scattered. A conservative estimate of the minimum size steam bubble is read from the graph as 3.5 cm. In summary a CST system to measure steam-water flow in a 4 inch schedule 40 pipe at 10 frames per second requires a 556 Curie cesium-137 source and is able to sense bubbles on the order of 3 cm in diameter. Conclusions and Future Research Looking specifically at the two-phase flow measurement problem there are several conclusions which can be drawn from the analysis of this report. It is possible to build a CST system to image two-phase flow inside steel pipes using a single source and two detectors. Such a CST system can be expected to function well in the presence of moderate ADC/D Figure 6.4: Results from the numerical evaluation of equation (5.6). 127 128 statistical noise, but may give erroneous results if systemÂ¬ atic errors exist in the data. The performance of the system is best at low void fractions and degrades with increasing void fraction. The pipe can place a practical upper limit on how high the void fraction can go and meaningful measurements still be made. If the pipe wall is too thick it may not be possible to design a functioning system at all. Some of these limitations can be over come by augmenting the CST design of this report with a few detectors to measure the attenuation of the transmitted source beam. The attenuation measurements complement the scatter measurements since they perform best at high void fractions. The potential advantages of a hybrid scatter- transmission system warrants further investigation and it is concluded that the next logical step would be the construcÂ¬ tion of a hybrid CST prototype system. A reasonable design goal is to image air/water of steam/water flows inside thin walled pipes with diameters in the range of 8 to 16 cm. Real-time data analysis and display could be achieved by integrating a super-micro computer into the system. Besides functioning as a useful laboratory instrument the prototype would provide a means for evaluating analytic techniques for boosting the effective system framing rate without requiring excessively large source strengths. Despite the evidence in this report that CST is a viable imaging technique, a prototype system is necessary to gain the acceptance of other researchers, particularly in the medical community. 129 The basic CST concepts discussed here require mono- energetic sources of gamma rays. This limits the choice of sources to only a few commercially available isotopes. Extension of the techniques to poly-energetic sources would allow closer optimization with respect to source energy and perhaps the use of X-ray machines. The higher photon outputs of X-ray machines appear necessary to achieve true stop action imaging. APPENDIX I MONTE-CARLO TRANSPORT PROGRAM The 3-D Monte-Carlo photon transport computer code is written in Microsoft FORTRAN for the IBM PC. This program is the one used to generate the "measured" data of this report. The program is written for a very specific geometry: a cylindrical pipe running parallel to the z-axis and illuminated by an external source located at the origin. The histories of a fixed number of photons are followed from birth through a fixed number of scattering interactions. It is assumed that only Compton scattering collisions and photoelectric absorption interactions take place. Three different source models are available: point, distributed (i.e. cylindrical), and pencil beam. All photons are born at the origin with the point source and pencil beam models. The distributed source model uniformly chooses the coordinates of a photon's birth from within a cylindrical source. The initial trajectory of all photons is chosen so that their paths cross the pipe. Any trajectories allowed by the source collimator which do not intersect the pipe are ignored. The point of intersection between a photon's initial trajectory and the pipe wall is used to start the photon tracking process. 130 131 The photon tracking process involves the repeated exeÂ¬ cution of three steps: 1) choose an interaction point inside the pipe along the current photon trajectory, 2) score the fraction of photons interacting at this point that would reach the detector, and 3) choose a new photon trajectory by sampling the Klein-Nishina differential scattering cross section. These three steps are executed once for each order of scattering. The importance of each photon is corrected in step 1 for the probability of the photon not escaping the pipe and for the probability that the photon is scattered instead of absorbed. The correction for the probability of not escaping is necessary since all photons are forced to have their next interaction inside the pipe. The scoring in step 2 takes into account the photon importance. It is often stated that Monte-Carlo does nothing more than trade the difficulties of solving a transport equation for an equally difficult geometry problem. Figure 1.1 describes a common geometry problem which must be solved by the program and should be helpful in understanding some of the logic and equations in the subroutines. An effort has been made to use variable names consistent with this figure. Two floppy disks containing this computer code are available by contacting the author through the Department of Nuclear Engineering Sciences, 202 Nuclear Sciences Center, Gainesville, Florida 32611. = R - (R sin9) P COS0 = V *v P o Figure I .1: Geometry and variable definitions used in the transport program. 132 APPENDIX II INSIGHTS AND FINER POINTS Parallel Versus Diverging Source Beams The reconstruction algorithms in Chapter 4 explicitly ignore the density terms in the arguments of the exponential attenuation factors. The expression for the number of singly scattered photons counted in a single MCA channel can be simplified to the following form Sc = pc*Constant*exp[ -p'*(uiLi + yQL0)] (II.1) where pc is the average density of the isogonic segment associated with channel c, p' is the average density in the flow cross section as seen by the incoming source and outgoing scattered photons, Pj_,U0 are the mass attenuation coefficients for the incoming source and outgoing scattered photons, respectively, and Lj_,LQ are the average distances traveled in the two- phase flow by the incoming and outgoing photons, respectively. In a fan beam source geometry pc and p' are not necessarily related since the paths of the incoming source photons and the outgoing scattered photons only traverse the isogonic segment c for relatively short distances (see Figure 2.5). The reconstruction algorithms take advantage of the 133 134 observation that when an individual pixel density is adjusted the associated attenuation factor for that pixel does not change. Obviously the attenuation factors for the other pixels do change if the incoming or scattered photons of those pixels cross the pixel whose density is changed. In the case of a parallel source beam the path of an outgoing scattered photon is colinear with the isogonic segment (see Figure 2.6). Equation (II.1),' rewritten for a parallel source beam, is Sc = Pc*Constant*exp[ -(p'y^L.^ + PCU0L0)] . (II.2) A reconstruction algorithm can more easily account for the linear-exponential dependence on pc in equation (II.2) than in equation (II.1). Figures II.1 and II.2 graphically depict equations (II.1) and (II.2) for an isogonic segment crossing the middle of the pipe in Figure 5.3. Equation (II.2) can exhibit a maximum with respect to pc if the pipe is large enough. The presence of a maximum can lead to an ambiguity in determining pc. The upshot of this is that the diverging source beam configuration leads to a more tightly coupled set of reconstruction equations than the parallel source beam. Figure II.1: Detector response for channel c versus the isogonic segment average density, p , and the average density seen by the source and scattered photons, p". Wide source beam illumination is used. 135 Figure II.2: Detector response for channel c versus the isogonic segment average density, p , and the average density seen by the source photons, p' Parallel sSurce beam illumination is used. 136 137 Multiple Scatter Versus Detector Position A convenient measure of the importance of the multiply scattered photons is the ratio of the singly scattered photons to multiply scattered photons averaged over the singly scattered photon spectrum. This ratio varies with detector position for a lot of reasons. The most important reason is that Compton scattering is non-isotropic. Figure II.3 shows the same experimental setup as Figure 5.3, but with three different detector positions. Figure II.4 shows the detected spectra for those positions. The single-to- multiple scatter ratio is highest for position 3 and lowest for position 2. Position 3 is favored for reducing the importance of multiple scattering, however, it may not be the best choice. Using position 3 would require a higher energy resolution detector in order to have the same ESILF as position 2 or 3, since the single scatter spectrum is narrower in position 3. Reference [35] discusses similar ways of optimizing the detector placement. Other factors to consider in placing the detector are detector efficiency versus angle of incidence, geometric attenuation, and spatial resolution. It is desirable to place the detector as close as possible to the object in order to increase the number of detected photons and the width (on the energy axis) of the singly scattered photon spectrum. The wider the spectrum the more spatial information is available for a given energy resolution (12,15) D (23,13) y A x, y) (4,9) Figure II.3: Three detector positions used to illustrate the variation of the multiple scatter component of the spectrum with detector position. 138 Counts [relative] Figure II.4: Spectra for the detectors in Figure II.3. 139 Figure II.4â€”continued. Detector 2 Single-to-multiple scatter ratio 7.67 Energy [keV] 600 140 Counts [relative] Figure II.4--continued. 141 142 detector assuming the detector and source sizes are negligible. The detector efficiency can be a strong function of a photon's angle of incidence to the detector. This is a problem when photons of the same energy (i.e. coming from the same isogonic segment) have widely ranging angles of incidence. The usual efficiency correction procedures assume all photons enter at the same angle and can not correct for such a variation. Moving the detector away from the object reduces the range of angle of incidence. Planar detectors are more sensitive to angle of incidence and less efficient than coaxial detectors. Therefore, it is recomÂ¬ mended that coaxial detectors be used in CST. Source Placement There are some practical matters to consider in the placement of the source. Photon economy is directly related to the source-to-object distance. The closer the source is to the object the larger the ratio of photons illuminating the object to the photons emitted. The trade-off in moving the source closer to the object is some loss in flexibility in the placement of the detectors. The detectors must be placed (or shielded) such that they do not have a line-of- sight view of the source. Another point to consider is backscatter from any surÂ¬ rounding structures. The wide source beam associated with small source-to-object distances is more likely to cause 143 problems by illuminating nearby objects. Scatter from these nearby objects could interfere with the measurements. Convergence Criterion The choice of the criterion used to end the reconstrucÂ¬ tion iterations can be tricky. Reference [36] looks at a similar situation where there are only a few projections available and considers two types of criteria: a) criteria based on properties of the reconstruction, and b) criteria based on the residual errors in the projections. Their conclusion based on tests is that type b) methods are better when there are only a few views. The type b) convergence criteria for CST corresponds to examining the difference between the measured spectra and the calculated spectra for the reconstruction. Appendix V contains' plots of t'he spectrum difference function and the discrepancy function (see the discussion of Figure 4.11) versus iteration number. These plots show that there is a very weak correlation between the rate of change in the spectrum difference function and how close the reconstrucÂ¬ tion is to the actual flow pattern as measured by the disÂ¬ crepancy function. This author prefers to use a fixed number of iterations instead of a convergence criterion. APPENDIX III RECONSTRUCTION COMPUTER PROGRAMS The reconstruction algorithms derived in Chapter 4 are implemented here in Turbo Pascal for an IBM PC. There is one main program for each of the four algorithms: Monte-Carlo Backprojection (MCB), ART, SIRT, and ILS. There are also four separate procedures which are the actual implementation of these algorithms. All remaining procedures used are common to the four main programs and are included at compile time. Running the Programs Prior to running the programs for a particular source- pipe-detector configuration a Pixel Spread Function (PSF) file must be generated for that particular setup. The PSF file contains the values of the PSF matrix introduced in Chapter 4. The program Make_PSF is used to generate this matrix using a Monte-Carlo technique and to store it in the proper format. The programs are designed to be run interactively. The first thing the programs do is read in the input data file and the PSF file. Reading in the PSF file takes approxiÂ¬ mately 3 minutes. The programs then run through 10 144 145 iterations of their respective algorithms (~30 minutes). The programs then prompt for the name of an output file. The input data file contains the information about the system geometry, descriptions of the source and detector, and the measured spectrum for each detector. The data fields are free format, however, data types must be adhered to. The file SAMPLE.INP contains enough descriptive comments that it can be used as a model for other problems. Three floppy disks are available containing these programs plus sample data by contacting the author through the Department of Nuclear Engineering Sciences, 202 Nuclear Sciences Center, Gainesville, Florida 32611. Flow of Programs The ART, SIRT, and ILS programs all have the same basic internal structure. The first statement reserves space on the heap for the PSF function. The next block of statements initialize all the global variables, loop counters, etc. Following the initialization code block is a repeat-until block which executes the particular reconstruction algorithm a fixed number of iterations. The final statements save the reconstructed density distribution and other results in an output file. APPENDIX IV MAXIMUM ENTROPY METHOD The maximum entropy method begins by redefining the fluid density distribution as a probability density distribution. The density distribution, p, usually sought is replaced by the probability function, E, where it k = i and (IV.l) fk i 0. The relationship between f and p is fk = Pk/(P*K) (IV.2) where p is the average density and K is the total number of pixels. The average density may not be known exactly and changing the normalization of f to an inequality allows using the fluid density in place of the average density. The entropy of a probability distribution is H(f) = -Ifkln(fkK) . (IV.3) The derivation given here is based on the work of Minerbo [28]. 146 147 The maximization is carried out by introducing a Lagrange multiplier XCT for each constraint condition (i.e. the y measured data) and forming the Lagrangian 'F(f) = H (f ) - Â£xg* (Sg - s| ( f ) ) (IV.4) Taking the functional derivative with respect to f while assuming that the transport matrix [equation (4.1)] is independent of f and setting it equal to zero yields H-l - ln(fkK) + (l/PfiuidK)Â£TgkV = Â° â€¢ (IV*5) k g Now it is argued that in order to have K linearly indepenÂ¬ dent density values requires that the term in brackets be zero for each value of k. Solving for fk yields fk = exp{ 1 - UgTgk}/K . (IV. 6) g The Xg's are found by substituting the density values given by equations (IV.6) and (IV.2) into the constraint conditions [i.e. equation (4.1)]. This gives G equations of the form Sg ' = Pfluid* Â£Tgk*exP{1 â€œ IXgTgk} k g where the - values are the measured data. y (IV.7) 148 It is interesting to note here the relationship between the maximum entropy solution and the least squares solution [see equation (4.28)]. When the number of pixels, K, is less than the number of measured data points, G, the lambdas are undefined. When the number of pixels is the same as the number of data points the maximum entropy solution is idenÂ¬ tical to the least squares solution with the weights set equal to one. When the number of pixels is larger than the number of data points the maximum entropy solution exists. The inforÂ¬ mation needed to define the extra pixels is provided by the entropy condition. The problem is that a nonlinear set of equations must be solved for the lambdas. The experience of this author with using the maximum entropy method under these conditions is that the solution fails to converge. Even for K only slightly greater than G leads to an ill- conditioned set of equations. The problem seems to be that the f^ values are too sensitive to small changes in the lambdas coupled with the fact that the transport matrix elements are functions of the density. An even more rigorous derivation of the maximum entropy solution can be made by explicitly including the dependence of the transport matrix on the density distribuÂ¬ tion. The equation analogous to equation (IV.5) arrived at by this route leads to a nonlinear relation between the f's and lambdas which just complicates the solution even more. APPENDIX V DATA This appendix is divided into two parts. The first part contains a representative selection of the data used in compiling Table 5.2. Figures V.l through V.15 are some of the spectra generated by the Monte-Carlo transport code. Figures V.16 through V.29 are plots of the reconstructed flow distributions for various cases of flow regime, void fraction, and noise content. Figures V.30 through V.34 are plots of the spectrum difference and discrepancy functions showing the different types of convergence behavior. The second part contains plots of the spectra from actual laboratory experiments. This data is not usable because of a lack of available sources to calibrate the detector and electronics. Figure V.35 shows a sketch of the experimental setup used. Figure V.36 illustrates the steps of correcting the spectrum for the detector efficiency and removing the multiply scattered photon spectrum. The multiple-scatter "background" is removed by fitting a 2nd order equation to the two small spectrum segments on either side of the single scatter peak. Also, shown in Figure V.37 is the spectrum generated by the Monte-Carlo photon transport code for the same annular setup. The slight misalignment of the two spectra is 149 150 probably due to poor energy calibration in the MCA and some difficulty in determining the detector efficiency. A three point energy calibration would be better than the two point procedure actually used. The detector efficiency was determined using sources of unknown pedigree. Although energy calibration of an MCA and detector efficiency deterÂ¬ mination are routine procedures in most counting laboratoÂ¬ ries this figure illustrates the need for care in performing these tasks. Counts [relative] <1) > â€¢H -p nJ râ€”I Q) U W -P a d o u Figure V.l: Detected spectra from a 30% void fraction annular flow regime. a) Spectrum of upper detector; b) Spectrum of lower detector. 151 Counts [relative] Figure V.2: Detected spectra from a 50% void fraction annular flow regime. a) Spectrum of upper detector; b) Spectrum of lower detector. 152 Counts [relative] Figure V.3: Detected spectra from 70% void fraction annular flow regime. a) Spectrum from the upper detector; b) Spectrum from the lower detector. 153 Counts [relative] Figure V.4: Detected spectra from a 30% void fraction inverted annular flow regime, a) Spectrum from the upper detector; b) Spectrum from the lower detector. 154 Counts [relative] 0.1 a) > â€¢H -P m iâ€”i 0 w 4J Â£ Â£ O O 200 Energy [keV] (b) "1 400 Figure V.5: Detected spectra from a 50% void fraction inverted annular flow regime, a) Spectrum of upper detector; b) Spectrum of lower detector. 155 Counts [relative] 1 400 Figure V.6: Detected spectra from a 70% void fraction inverted annular flow regime, a) Spectrum of upper detector; b) Spectrum of lower detector. 156 Counts [relative] 0.1 (U > â€¢H -P rtf râ€”4 o W -P C 3 O U (b) â€œ1 400 Figure V. 7 : Detected spectra from a 30% void fraction bubbly flow regime. a) Spectrum of upper detector; b) Spectrum of lower detector. 157 Figure V.8: Detected spectra from a 50% void fraction bubbly flow regime- a) Spectrum of upper detector; b) Spectrum of lower detector. 400 158 Counts [relative] 0 Figure .9: Detected spectra from a 70% void fraction bubbly flow regime. a) Spectrum of upper detector; b) Spectrum of lower detector. 159 Counts [relative] 0.1 CD > +j Ctf râ€”I CD M W 4J G G O U 0 (b) Figure V.10: Detected spectra from a 30% void fraction stratified symmetric flow regime, a) Spectrum from upper detector; b) Spectrum from lower detector. 400 160 Counts [relative] -p (0 iâ€”i 0) u m -p c 3 o o Figure V.11: Detected spectra from a 50% void fraction stratified symmetric flow regime, a) Spectrum from upper detector; b) Spectrum from lower detector. 161 Counts [relative] Figure V.12: Detected spectra from a 70% void fraction stratified symmetric flow regime, a) Spectrum from upper detector; b) Spectrum from lower detector. 162 Counts [relative] Figure V.13: Detected spectra from a 30% void a) Spectrum from upper detector; fraction stratified asymmetric flow regime, b) Spectrum from lower detector. 163 Counts [relative] Figure V.14: Detected spectra from a 50% void fraction stratified asymmetric flow regime, a) Spectrum of upper detector; b) Spectrum of lower detector. 164 Counts [relative] Figure V.15: Detected spectra from a 70% void fraction stratified asymmetric flow regime, a) Spectrum of upper detector; b) Spectrum of lower detector. 165 166 Annular flow regine as it should appear when correctly reconstructed, a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. Figure V.16: 167 Figure V.17: Reconstruction of annular flow for noise case 1. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 168 Figure V.18: Reconstruction of annular flow for noise case 2. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 169 Figure V.19: Reconstruction of annular flow for noise case 3. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 170 Inverted annular flow regime as it appear when correctly reconstructed void fraction; b) 50% void fraction void fraction. should . a) 30% ; c ) 7 0 % Figure V.20: 171 (c) Figure V.21: Reconstruction of the inverted annular flow for noise case 1. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 172 Figure V.22: Bubbly flow regime as it should appear when correctly reconstructed, a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 173 Figure V.23: Reconstruction of the bubbly flow for noise case 1. a) 30% void fraction; b) 50% void fraction; c) 70 % void fraction. 174 (c) .BBT"' Jig* ty'.:;'. â€¢?AÃf. V-1 UM.W â€¢;H{ M?S_. . .... . Â¡?T;r, Ã/Ã¼tVfVÃ ;â€¢: \y-C...Â«5 ,:..\/ â€¢., .ÃXÃ* t,â€™. â€¢'â– ?'â€¢ vÃ â€¢. vÂ«Ã â€˜ Vy-CS-A *S T.W>\>. Figure V.24: Reconstruction of the bubbly flow for noise case 2. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 175 Figure V.25: Reconstruction of the bubbly flow for noise case 3. a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. 176 Stratified symmetric flow regime as it should appear when correctly reconstructed, a) 30% void fraction; b) 50% void fraction; c) 70% void fraction. Figure V.26: 177 (b) (c) Figure V.27: Reconstruction of the stratified symmetri for noise case 1. a) 30% void fraction; void fraction; c) 70% void fraction. c flow b) 50% 178 Stratified asymmetric flow regime as i appear when correctly reconstructed, void fraction; b) 50% void fraction; c void fraction. t should a) 30% ) 70% Figure V.28: 179 (b) (c) Figure V.29: Reconstruction of for noise case 1. void fraction; c) the stratified asymmetric flow a) 30% void fraction; b) 50% 70% void fraction. Figure V.30: Example of smooth convergence. 180 1 cu o tÃ cu tÃ O 4-4 4-4 â€¢r-Ã Q e d tÃ -p o (U Ch cn 0 Iteration â€œI 10 Figure V.31: Example of overiteration with resulting slight divergence. 181 Discrepancy Figure V.32: Example of overiteration with resulting divergence. 182 Figure V.33: Example of overiteration with oscilating divergence. 183 Figure V.34: Example of failure to converge. 184 Coax Detector Figure V.35: Experimental source-pipe-detector configuration. 185 6- 5- (a) Figure V.36: Experimental Spectrum from the annular plastic phantom. a) raw data; b) detector efficiency corrected data with dashed line indicating fit to multiple-scatter background. 186 Figure V.36â€”continued. 1 i 1 1 "I 400 600 Energy [KeV] (b) 187 Figure V.37: Comparison of Monte-Carlo generated spectrum to experimental spectrum from the annular flow phantom. 188 7 Figure V.38: Experimental spectrum from a stratified flow phantom. 189 Counts 7 Energy [KeV] Figure V.39: Experimental spectrum from a bubbly flow phantom. 190 REFERENCES 1. Hewitt, G.F., Measurement of Two-Phase Flow Parameters, Academic Press, London (1978). 2. Heiselmann, H.W., Davis, L.L., "Measurement Techniques for Transient Two Phase Coolant Flow During Blowdown," Proceedings of Two-Phase Flow and Heat Transfer Symposium-Workshop, Fort Lauderdale, Florida, (October 18-20, 1976), (distributed by Pergamon Press). 3. Delhaye, J.M., Cognet, G., "Measuring Techniques in Gas- Liquid Two-Phase Flows," Proceedings of International Union of Theoretical and Applied Mechanics Symposium, Nancy, France, July 5-8, 1983, Springer-Verlag, Berlin (1984). 4. Bergles, A.E., Collier, J.G., Delhaye, J.M., Hewitt, G.F., Mayinger, F., Two-Phase Flow and Heat Transfer in the Power and Process Industries, Hemisphere Publishing Company, Washington, D.C. (1981). 5. Bryant, D.T., Payne, J.A., Firmin, D.N., Longmore, D.B., "Measurement of Flow with NMR Imaging Using a Gradient Pulse and Phase Difference Technique," Journal of Computer Assisted Tomography, Vol. 8, pg. 588 (1984). 6. Jacobs, A. M., "Nonintrusive transaxial tomography technique for velocity profile measurement," Proceedings of the Society of Photo-Optical Instrumentation Engineers, Vol. 182 Imaging Applications for Automated Industrial Inspection and Assembly, pp. 165-170 (1979). 7. Jacobs, A. M., "First Report on the Evaluation of Gamma- Ray Scattering Techniques for the Tomographic Measurement of Water Two-Phase Density Distribution," EG&G Consultant Agreement No. K-8118-TKP-70-79, EG&G Idaho, Idaho Falls, Idaho, July (1979). 8. Peschmann, K.R., Napel, S., Couch, J.L., Rand, R.E., Alei, R., Ackelsberg, S.M., Gould, R, Boyd, D.P., "High- Speed Computed Tomography: Sytems and Performance," Applied Optics, Vol. 24, No. 23, pp. 4052-60 (1985). 9. Kondic, N. N., Jacobs, A. M., Anghaie, S., "NonÂ¬ intrusive Determination of Mass Density Distribution by Gamma-Ray Scattering: Historical Survey and 191 192 Perspectives," American Society of Mechanical Engineers- Fluids Engineering Conference, Cavitation and Multiphase Flow Forum, Albuquerque, New Mexico, June (1985). 10. Odeblad, E., Norhagen, A., "Electron Density in a Localized Volume by Compton Scattering," Acta RadiolÃ³gica, Vol. 45, p. 161 (1959). 11. Lale, P. G., "The Examination of Internal Tissues Using Gamma-Ray Scatter with a Possible Extension to Megavolt Radiography," Physics in Medicine and Biology, Vol. 4, p. 159 (1959). 12. Lale, P. G., "Examination of Internal Tissues by High- Energy Scattered X-Radiation," Radiology, Vol. 90, p.510 (1968) . 13. Kondic, N. N., Hahn, O. J., "Theory and Application of the Parallel and Diverging Radiation Beam Method in Two- Phase Systems," Paris, France, 4th International Heat Transfer Conference, Vol. 7, MT-1.4 (1970). 14. Farmer, F. T., Collins, M. P., "A New Approach to the Determination of Anatomical Cross-sections of the Body by Compton Scattering of Gamma-Rays," Physics in Medicine and Biology, Vol. 16, No. 4, p. 557 (1971). 15. Farmer, F. T., Collins, M. P., "A Further Appraisal of the Compton Scattering Method for Determining Anatomical Cross-sections of the Body," Physics in Medicine and Biology, Vol. 19, No. 6, p. 808 (1974). 16. Anghaie, S., Jacobs, A. M., Kondic, N. N., "Application of Gamma Scattering Compton Profile for Two-Phase Systems Densitometry," invited paper, American Nuclear Society Winter Annual Meeting, San Francisco, California (November, 1981). 17. Anghaie, S., "Accuracy Improvement for Gamma-Ray Techniques in Two-Phase Measurements," Ph.D. Thesis, Pennsylvania State University, Nuclear Engineering (1982). 18. Kondic, N., "Density Field Determination by an External Stationary Radiation Source Using a Kernal Techniques," Symposium on Polyphase Flow Measurement, ASME Winter Annual Meeting, San Francisco (December, 1973). 19. Banerjee, S., Lahey, R. T., Jr., "Advances in Two-Phase Flow Instrumentation," Advances in Nuclear Science and Technology, Vol. 13, Plenum Publishing Co., New York, New York, pp. 227-414 (1981). 193 20. Hussein, E., "A Fast Neutron Scattering Technique for Measurement of Local Void Fraction Distribution in Two- Phase Flow," Ph.D. Thesis, McMaster University (1983). 21. Bodette, D.E., Jacobs, A.M., "Tomographic Two-Phase Flow Distribution Measurement Using Gamma-Ray Scattering," Proceedings of the Symposium on New Technology in Nuclear Power Plant Instrumentation and Control, Washington, D.C., November 28-30, 1984 (distributed by Instrument Society of America). 22. Battista, J., SantÃ³n, L., Bronskill, M., "Compton- Scatter Imaging of Transverse Sections: Corrections for Multiple-Scatter and Attenuations," Physics in Medicine and Biology, Vol. 22, p. 229 (1977). 23. Battista, J., Bronskill, M., "Compton-Scatter Tissue Densitometry: Calculation of Single and Multiple-Scatter Photon Fluences," Physics in Medicine and Biology, Vol. 23, p. 1 (1978). 24. Kondic, N. N., Lassahn, G., "Nonintrusive Density Distribution Measurement in Dynamic High Temperature Systems," 24th International Instrumentation Symposium, Albuquerque, New Mexico (May, 1978). 25. Kondic, N.N., Jacobs, A.M., Ebert, D., "Three- Dimensional Density Field Distribution by External Stationary Detectors and Gamma Sources Using Selective Scattering," invited paper, 2nd International Topical Meeting on Nuclear Reactor Thermal-Hydraulics, Santa Barbara, California (January 11-14, 1981). 26. Tretiak, O. J., Metz, C., "The Exponential Radon Transform," SIAM Journal of Applied Mathematics, Vol. 39, No. 2, pp. 341-354 (1980). 27. Gullberg, G. T., Budinger, T. F., "The Use of Filtering Methods to Compensate for Constant Attenuation in Single-Photon Emission Computed Tomography," IEEE Transactions on Biomedical Engineering, Vol. BME-28, No. 2, pp. 142-157 (1981) . 28. Minerbo, G., "MENT: A Maximum Entropy Algorithm for Reconstructing a Source from Projection Data," Computer Graphics and Image Processing, Vol. 10, pg. 48 (1979). 29. Llacer, Jorge, "Tomographic Image Reconstruction by Eigenvector Decomposition: Its Limitations and Areas of Applicability," IEEE Transactions on Medical Imaging, Vol. MI-1, No. 1, pg. 34 (July, 1982). 194 30. Herman, G. T., Image Reconstruction from Projections The Fundamentals of Computerized Tomography, Academic Press, New York (1980). 31. Gilbert, P., "Iterative Methods for the Three- Dimensional Reconstruction of an Object from ProjecÂ¬ tions," Journal of Theoretical Biology, Vol. 36, pp. 105-17 (1972). 32. Joseph, P.M., "An Improved Algorithm for Reprojecting Rays Through Pixel Images," IEEE Transactions on Medical Imaging, Vol. MI-1, No. 3, pg. 192 (November, 1982). 33. Harms, A.A., Forrest, V.F., "Dynamic Effects in Radiation Diagnosis of Fluctuating Voids," Nuclear Science and Engineering, Vol. 46, pp. 408-413 (1971). 34. Hubbell, J.H., "Photon Mass Attenuation and Energy- absorption Coefficients from 1 keV to 20 MeV," International Journal of Applied Radiation and Isotopes, Vol. 33, pp. 1269-1290 (1982). 35. Anghaie, S., Finlon, K.D., Jacobs, A.M., "Optimizing and Compton Profile Approach to Two-Phase Distribution Measurement," Third Multiphase Flow and Heat Transfer Symposium Workshop, Miami, Florida (April 18-20, 1983). 36. Minerbo, G.N., Sanderson, J.G., "Reconstruction of a Source from a Few (2 or 3) Projections," LASL informal report LA-6747-MS (March, 1977). BIBLIOGRAPHY Abramowitz, M., Stegun, I. A., Handbook of Mathematical Functions, ninth edition, Dover Publications, Inc., New York (1972). Barth, K., Deimling, M., Fritschy, P., Mueller, E., Reinhardt, E.R., "Flow Measurement by Magnetic Resonance Imaging," Application of Optical Instrumentation in Medicine XIII - Medical Image Production, Processing, and Display, ed. by R.H. Schneider and S.J. Dwyer, Published by Society of Photo-Optical Instrumentation Engineers, Vol. 535, pp. 261-267 (1985). Gullberg, G. T., Budinger, T. F., "Three-Dimensional Reconstruction in Nuclear Medicine Emission Imaging," IEEE Transactions on Nuclear Science, Vol. NS-21, p. 2 (1974). Heinzerling, J., Schlindwein, M., "Non-Linear Techniques in Multi-Spectral X-Ray Imaging," IEEE Transactions on Nuclear Science, Vol. NS-27, No. 2, pg. 961 (April, 1980). Holt, R.S., "Compton Imaging," Endeavour, New Series, Vol. 9, No. 2, pg. 97 (1985). Hussein, E., Banerjee, S., Meneley, D.A., "A New Fast Neutron Scattering Technique for Local Void Fraction Measurement in Two-Phase Flow," Thermal Hydraulics of Nuclear Reactors, Vol. 2 (1983). Hussein, E.M.A., "Numerical Method for Radiation Scatter Imaging," IEEE/Seventh Annual Conference of the Engineering in Medicine and Biology Society (1985). King, J.D., "Coal Flowmeter Development," (CONF-840694) From AR and TD and Suface Coal Gasification InstrumentaÂ¬ tion Projects Contractors' Meeting, Morgantown, WV, (June 13, 1984). Rousseau, J.C., Czerny, J., Riegel, B., "Void Fraction Measurement by Neutron Absorption or Scattering Methods," OELD/NEA Specialists Meeting on Transient Two- Phase Flow, Toronto, Canada (August, 1976). 195 BIOGRAPHICAL SKETCH David E. Bodette was born on September 14, 1958 on Midway Island, Hawaii. In June 1976, he graduated from Eau Gallie High School in Melbourne, Florida. In December 1980, he received a Bachelor of Science degree in nuclear engineering from the University of Florida. He was awarded an Institute of Nuclear Power Operations Fellowship in September 1982 and in December 1982 received a Master of Science degree in nuclear engineering from the University of Florida. In 1983, David Bodette was awarded a Department of Energy Graduate Traineeship. As part of this traineeship he spent three months at EG&G Idaho conducting experiments with two-phase flow densitometry. He expects to receive a Doctor of Philosophy from the University of Florida in December 1986. 196 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Alan M. Jacobs Professor of Nuclear Engineering Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of-Etoctorof Philosopi Edward E. Carroll, Professor of Nuclear Sciences engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. â€¢/ ' //â€¢â€¢ r/ â– -v - /-â– â– ' Edward T. Dugan Associate Professor of Nuclear Engineering Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Robert J. Professor nrahan Chemistry I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scppÃ©^and quality, as a dissertation for the degree of Doctor/Ã³f Philosophy, Vernon P. Roan, Jr. Professor of Mechanical Engineering This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1986 /ibhiÃ¼jf Ct' Dean, College of Engineering Dean, Graduate School Internet Distribution Consent Agreement In reference to the following dissertation: AUTHOR: Bodette, David TITLE: Tomographic two-phase flow measurement using Compton scattering of gamma rays / (record number: 940973) PUBLICATION DATE: 1986 y)A\/rd l< , as copyright holder for the aforementioned dissertation, hereby I, _ grant specific and limited archive and distribution rights to the Board of Trustees of the University of Florida and its agents. I authorize the University of Florida to digitize and distribute the dissertation described above for nonprofit, educational purposes via the Internet or successive technologies. This is a non-exclusive grant of permissions for specific off-line and on-line uses for an indefinite term. Off-line uses shall be limited to those specifically allowed by "Fair Use" as prescribed by the terms of United States copyright legislation (cf, Title 17, U.S. Code) as well as to the maintenance and preservation of a digital archive copy. Digitization allows the University of Florida to generate image- and text-based versions as appropriate and to provide and enhance access using search software. This grant of permissions prohibits use of the digitized versions for commercialese or profit. Signature of Copyright Holder Â¿I l^oDzTr^ Printed or Typed Name of Copyright Holder/Licensee Personal information blurred / f ZoQÂ£> Date of Signature Please print, sign and return to: Cathleen Martyniak UF Dissertation Project Preservation Department University of Florida Libraries P.O. Box 117007 Gainesville, FL 32611-7007 7/1/2008 |