- Email: [email protected]
A novel X-ray computed tomography method for fast measurement of multiphase flow
Chemical Engineering Science 62 (2007) 4325 – 4335 www.elsevier.com/locate/ces
A novel X-ray computed tomography method for fast measurement of multiphase flow Changning Wu, Yi Cheng ∗ , Yulong Ding, Fei Wei, Yong Jin Department of Chemical Engineering, Beijing Key Laboratory of Green Chemical Reaction Engineering and Technology, Tsinghua University, Beijing 100084, PR China Received 12 May 2006; received in revised form 18 January 2007; accepted 25 April 2007 Available online 6 May 2007
Abstract A new approach based on an improved genetic algorithm (GA) was proposed to implement the image reconstruction when using X-ray computed tomography (XCT) for the application of fast measurement of multiphase flow dynamics. Instead of directly using a traditional XCT, we pursued to develop a different discrete tomography (DT) method, aiming to achieve a high resolution in time during the measurements with only limited projection data. The proposed method assumed that the interested multiphase flow can be simplified as having distinct dense and dilute phases so that the local phase concentration can be binary-coded, e.g., 0 or 1 in a gas bubbling system. The mathematical problem under these circumstances is strongly ill-posed, and thus tackled with an optimization approach, i.e., a GA incorporated with the underlying physics as some constraints. The numerical simulations mimicking the physical measurements demonstrated the feasibility of the new approach, namely GA-XCT, especially with high robustness to the noise. Experiments were performed to simulate a transient measurement on the gas bubbles in water, with a portable X-ray tube and a 2D plane detector as the hardware and a static object rotating in between. The results further provided the validation of the GA-XCT being superior to the conventional algorithm, e.g., filtered back-projection (FBP) technique, in dealing with the tomography of multiphase system with binary local density field. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Multiphase flow measurement; X-ray computed tomography; Discrete tomography; Genetic algorithm; Limited projection data
1. Introduction Multiphase flow is prevalent in petroleum refining, synthesis gas conversion to fuels and chemicals, bulk commodity chemicals, fine chemicals, manufacture of polymers, biotechnology and pollution control (Dudukovic, 2002). However, the complexity of the underlying dynamics due to the turbulence and nonlinear contact between phases hinders the theoretical understanding on the multiphase flows. Subsequently, the design, control and scaling-up of multiphase flow equipment are still
Abbreviations: ART, algebraic reconstruction technique; DT, discrete tomography; ECT, electrical capacitance tomography; FBP, filtered backprojection; GA, genetic algorithm; -CT, gamma computed tomography; IGP, integer programming problem; LDT, limited data tomography; PT, process tomography; XCT, X-ray computed tomography ∗ Corresponding author. Fax: +86 10 62772051. E-mail address: [email protected] (Y. Cheng). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.04.026
strongly empirical or semi-empirical. Measurement techniques for multiphase flows, especially the non-invasive techniques, always play an important role in promoting the development of multiphase flow theories and engineering applications (Chaouki et al., 1997). With the rapid development of fundamental studies, more stringent demands are brought out for measurement techniques, such as non-invasive measurement with 2D/3D filed information, high resolution in space and/or in time, etc. Under these circumstances, process tomography (PT) finds their unique role in studying the multiphase flow phenomena encountered in chemical engineering fields, involving a number of tomographic techniques, such as X-ray, -ray, optical, ultrasonic, electrical, and nuclear magnetic resonance imaging (Beck and Williams, 1996; Du et al., 2004; Dyakowski et al., 1997; George et al., 2000; Kumar et al., 1997; Rowe, 1971; Seville et al., 1986; Warsito and Fan, 2003a,b,c; Warsito et al., 1999). Among them, X-ray computed tomography (XCT) and electrical capacitance
4326
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
tomography (ECT) are regarded as the two types of the most attractive techniques (Warsito and Fan, 2003c; Werther, 1999). Compared with the ECT technique, XCT (and/or -CT) provides higher spatial resolution up to 1% of a column diameter and even more details at the scale of 200.400 m. Especially, XCT applied in the medical systems has been a mature technique both in theory and in hardware design, also acknowledged as a useful and important technique in the study of multiphase flows. The initial application of the basic principle of CT to the multiphase flow was done by the researchers at Shell Oil Co. Inc. (Bartholomew and Casagrande, 1957). In that work, -ray was used to measure the suspension density in an industrial riser. Different from the current CT method, an empirical density function was firstly assumed, and then the parameters in the function were correlated using a lot of data measured. Similar empirical method to implement a tomography measurement was also used by Weinstein’s group (Weinstein et al., 1984). Medical XCTs started to be applied to study multiphase flows in 1990s (Banholzer et al., 1987; Kantzas, 1994; Kantzas and Kalogeraki, 1996; Kantzas et al., 1997; Holoboff et al., 1995; Grassler and Wirth, 1999, 2000). However, time-averaged measurement results from XCTs were mostly reported in the literature as above. This common feature of XCTs as the low temporal resolution is not applicable in studying the time evolving phenomena of multiphase flows. More detailed researches about applying medical/commercial CT scanners to multiphase flow systems can be found in the review articles (Chaouki et al., 1997; Werther, 1999). Essentially, limited temporal resolution of XCTs comes from a basic request to the image reconstruction for a medical CT measurement, that is, a complete data set for the image reconstruction is inevitably demanded considering the unknown density distribution in human bodies and requested spatial resolution. This consequently causes the rotating action of either X-ray source(s) or the detectors or both of them around the object, which increases the data acquisition time with decreased transient scanning capability for one cross-section imaging. Developed to the fifth generation CT, the hardware has no mechanically moving parts at all, but the electron beams scan around the annular target. In this way, the CT scanner can obtain the data per slice very fast, for instance, in 10 ms (Tian et al., 2003). Regarding the few applications of fast XCT in multiphase flow measurement, the goal can be achieved using more fixed X-ray sources and detectors. Harvel et al. (1996, 1999) built a high speed XCT system to determine the flow regime and the 2D void fraction distribution for a vertical annulus two-phase flow. To increase the time resolution, 18 X-ray sources and 122 CdWO4 detectors were arranged around the test channel. The cross-section sampling time for each slice was 4.0 ms (i.e., 250 frames/s) with a spatial resolution of 2.0 mm. Kai et al. (2000) used eight electron guns and 122 detectors in the X-ray scanner to improve the time resolution. The maximum scan rate of this system could reach 250 slices per second, too. The measurement area was limited to 50 mm in diameter for the requirement of spatial resolution.
For a better understanding of the dynamics in multiphase flows, there is obviously a current need/challenge in improving the response time of the measurement while keeping the high resolution for instantaneous observation. While, compared with a medical XCT, the XCT applied in process engineering is of the limitation in the number of sensors arranged, and consequently in the relative lower spatial resolution due to the limited integral measurement data obtained. The need for high-speed imaging must prompt the number of the sensors to be limited as few as possible, but with sufficient accuracy for process quantification from the limited experimental integral data without sacrificing the time resolution for reliable real time measurements (Warsito and Fan, 2003c). Improving the resolution in time but maintaining the high resolution in space together with not much cost is therefore very challenging to further develop the CT technique for industrial applications. Different from the original idea of an XCT applied in medical purposes, it would be realized that there is often no need to study the gradient change of the phase holdup or the density in many multiphase systems. Distinguishing the dense and dilute phases would be enough. In this sense, some instrumentation/computing power could be saved to increase the process speed, but the resolution would be still high. Such an approach is expected to have significant impact on the multiphase study, even serving to the relevant industry. From this view point, we aim to develop a special fast XCT technique in this work, particularly applied in certain multiphase flow systems, where the local phase distribution shows 0 or 1 feature at this stage. The gas bubbling phenomena in bubble columns and slurry reactors would be the appropriate candidates for the potential applications. Owing to the simplified model system, limited projection data would be sufficient in meeting the request of high spatial resolution and temporal resolution as well. That is to say, limited data tomography (LDT) is concerned in this work, which refers to the situation where ideal tomography is not possible due to insufficient or inappropriately placed parallel-beam/fanbeam projections (for 2D inverse problem). Since the conventional CT algorithms are not applicable when the projection data are very limited, developing a new methodology to efficiently reconstruct the cross-sectional images of bubbling phenomena becomes the major task in this work. 2. Mathematical method: GA-XCT for LDT 2.1. Problem re-definition The principle of the radiation attenuation is the first foundation of XCT. Assuming that the X-rays are monoenergetic from a parallel beam, the residual intensity of X-rays is attenuated by the materials that the X-rays pass through: I = I0 e −
L m ds
,
(1)
where m is the mass attenuation coefficient (m2 /kg) mainly due to the photoelectric, Compton and pair production effects; is the density of the object (kg/m3 ); ds represents infinitesimal
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4327
thin layer of the object (m). The corresponding expression of the projection data is I0 bL = ln (2) = m ds. I L For air–water binary system, the attenuation of X-rays in air can be neglected due to the small density of air compared with water. The accumulated transmission thickness of water by X-ray along line L (Ll ) can be calculated by I0 bL = ln (3) = m ds = m,l l Ll . I L The image reconstruction techniques use the measured projection data as the input to calculate the density distribution of the desired cross-section of the investigated sample as the output. Accordingly, the 2D image of the desired cross-section can be obtained. The mostly applied reconstruction techniques are divided into two categories, the back-projection reconstruction and the iterative reconstruction. As an example of the back-projection reconstruction technique, filtered back-projection (FBP) technique is used in most commercial medical scanners and has proved to be extremely accurate and amenable to fast implementation. This technique can give a rather straightforward intuitive rationale because each projection represents a nearly independent measurement of the object. FBP is applied as the conventional method to be compared with the new approach proposed in this work in the context. Algebraic reconstruction technique (ART) is an entirely different approach for tomographic imaging assuming that the cross-section consists of an array of unknowns and the setting up algebraic equations for the unknowns in terms of the measured projection data (Kak and Slaney, 1988). ART is capable for the situations such as the limited data of projections and the uneven angle-distributed projections. All ART methods are iterative procedures. As shown in Fig. 1, a 2D object is placed between the X-ray source and the detector. The investigated area is discretized by rectangular girds or triangular ones. The mathematical expression of the tomographic problem is given by the equations: b = Ax,
i.e., bk =
M×N
akj xj , k = 1, 2, . . . , K,
(4)
j =1
where K is the total number of the projections, M × N is the total number of pixels across the area, akj means the weighting factor for the contribution of pixel j to the projection element bk . The goal is to determine the unknown image (x) when the experimental projections (b) are available. The weighting factor, akj , is calculated based on the geometrical considerations as the intersection length/area of the kth projection ray with the image pixel j. They depend only on the geometry of the problem and therefore need to be calculated only once. Radon (1917) demonstrated that the mathematical problem in tomography is an inverse problem and that a unique solution is guaranteed provided that data from an infinite number of
Fig. 1. Illustration of the principle of XCT.
views be available. However, for LDT, the problem is ill-posed; we re-defined it to an optimization problem, to find the best solution for the following function: min {Ax − b22 }.
xj ∈R
(5)
For the specific application in this work, xj ∈ {0, 1}, where “0” represents air and “1” represents water. Though the model system is simplified greatly, the resulted nonlinear integer programming problem (IGP) offers a bigger computational complexity, in general with a greater computational cost, which even makes the solution unfeasible in many cases. Caballero et al. (2002) claimed that when solving a complex IGP problem, genetic algorithm (GA) can be an efficient option from the point of view of computational cost as well as the quality of the solution. 2.2. Improved GA GA is a global search and optimization algorithm based on the idea of evolution theory that individuals having a high value of quality will survive to the next generation with a greater probability. In general, GA uses guided random search method, instead of calculating the derivatives or gradients. Therefore, GA shows very robust to the different problems. In a GA approach, an “individual” represents a variable for our problem, i.e., x in formula (5), with a corresponding encoded genome for each individual. Genomes are vectors of genes and the components of a genome are closely related to the parameters of the given problem. A set of genomes is called a population. Genes are sets of real numbers or bits. A GA works with the encoded genome instead of the variables themselves to minimize the objective function, which is often called fitness function, dependent on the GA. The fitness function generates an output of a real value from each input individual, decoded from its corresponding genome. The following operations are defined for a population: the selection
4328
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
operation selects the best genomes from the modified population for the next generation; the crossover operation exchanges a number of genes of two genomes of a population; the mutation operation changes two randomly selected genes of a genome. During each iteration, genomes are evaluated and some are selected as parents based on their fitness value. Pairs of parents are recombined by means of genetic operators and the current population is improved in the sense that some members of the population are replaced by individuals with better fitness values obtained as offspring. The process continues until some stopping condition is satisfied. The algorithm proposed in this work is problem-oriented. The fitness function for GA approach is stated as formula (5), which belongs to the problem of an integer goal programming due to the interested binary value system. A binary coding is thus utilized to represent the parameters of the fitness function. That is to say, each gene, xj , can only be 0 or 1. And, an individual or a genome is just a set of solution to the problem re-defined. A common selection approach assigns a probability of selection, Pj , to each individual j based on its fitness value. Ranking selection operator, only requiring the fitness function to map the solutions to a partially ordered set, is chosen in this paper. Ranking methods assign Pj based on the rank of solution j when all solutions are sorted. Normalized geometric ranking defines Pj for each individual by (Joines and Houck, 1994): P [Selecting the j th individual] = q (1 − q)r−1
(6)
where q is the probability of selecting the best individual, r is the rank of the individual (where 1 is the best), and q = q/(1 − (1 − q)NP ), NP is the population size. The crossover operator takes two individuals and produces two new individuals while the mutation operator alters one individual to produce a single new solution. The application of these two basic types of operators and their derivatives depends on the genome representation used. The one-point crossover works by selecting a common crossover point in each pair of randomly selected two parents and then swapping the corresponding subtrees like standard crossover, with the generation of two genomes of the offspring. If there are two parents S1 = (S11 , S12 , . . . , S1N ) and S2 = (S21 , S22 , . . . , S2N ), their children will be (S11 , . . . , S1m , S2,m+1 , . . . , S2N ) and (S21 , . . . , S2m , S1,m+1 , . . . , S1N ). That is, a head and a tail of an offspring genome will be taken from different parents. For the mentioned problem-oriented nature of the proposed methodology, some prior knowledge has been incorporated into the GA operators. Firstly, a node-based mutation operator is proposed to consider the dependence of the gene mutation rate on the status of its neighbors. The detailed implementation is shown in Fig. 2a. The numbers of “0” and “1” elements are counted around a randomly chosen node. If all the elements around the node are 0s or 1s, the mutation rate will be set to zero, meaning no change for the next iteration/generation; in other cases, the mutation rate of each element, p(0) or p(1) will be assigned according to the proportion of 0s and 1s, which
is given as N (1) ; [N (0) + N (1)]N (0) N (0) , p(1) = [N (0) + N (1)]N (1) p(0) =
(7)
where N (0) and N (1) represent the number of elements 0s and 1s around the selected node, respectively. This treatment imposes the spatial correlations among the status of the neighbored grids, which can be regarded as the continuity assumption of the gas inside each bubble. Compared with the other iterative techniques in the image reconstruction application, the above method incorporates some intelligence into the algorithm. In this way, the GA is no longer a random search method, but performs a directed search to recognize the “bubbles” or maintain the bubble existence, which in turn speeds up the iteration procedure. Secondly, a swapping-based mutation operator is used to increase the speed of convergence. As illustrated in Fig. 2b, to implement a status change from A to C, the iteration procedure of a GA generally steps from A to B then to C. Correspondingly, the fitness function first increases and then decreases, which costs additional computation time. The proposed swappingbased mutation operator allows for the right evaluation on this situation and makes the iteration stepping from A to C directly by simply swapping the corresponding elements. This becomes an auxiliary criterion to help the algorithm being guided. 2.3. Implementation of the GA-XCT method The GA approach to the XCT measurement for multiphase flows is simply named as GA-XCT in this work, which includes four modules in the programming: (1) finite element gridding of the cross-section (i.e., triangular grids); (2) computation of the weighting factor matrix, A, according to the geometry; (3) computation/measurement result of the projections, b; (4) image reconstruction using the improved GA approach until convergence. The GA-XCT algorithm is coded in C + +, in which a GALib developed by MIT researchers (http://lancet.mit.edu/ga) is employed. GA-Lib provides a framework and the basic functions for general applications of GA. The new modifications described above are superimposed on the GA-Lib successfully. Besides, multi-population GA approach is applied to speed up the calculation and convergence of the IGP problem instead of a common GA method. To start the GA-based image reconstruction, an initial population must be provided for GA, which is usually generated in a random manner for the entire population. However, since the algorithm can iteratively improve existing solutions, the beginning population can be seeded with potentially good solutions, with the remainder of the population being randomly generated.
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4329
Fig. 2. Illustration of GA operators and zooming operation: (a) the node-based mutation operator; (b) the swapping-based mutation operator; and (c) the grid refinement.
2.4. Method to improve the spatial resolution Based on the triangle meshing method, the GA approach can be conveniently applied to higher resolution image restoration by refining the grids in the way shown in Fig. 2c. It is simply named as “zooming operation” in this work for easy expression. The solutions under the coarser scale grids are taken as a better initial guess to the problem solutions under the finer scale grids. By inheriting the solutions from the coarser meshes, it can be found in following context that the image reconstruction results under finer resolutions in space can be achieved stably with good accuracy. 3. Validation by numerical simulations 3.1. The test phantom and simulation details To evaluate the performance of the proposed GA-XCT method, numerical simulations are carried out using various simulated distribution functions (i.e., phantoms) as the references. Some (near-) round objects with different sizes are arranged randomly in a bigger round area, representing the situation of the cross-section in a gas–liquid bubbling column. At this stage, we start from few bubbles cases in the measurement plane, such as four bubbles for the simpler case and ten bubbles for the more complicated case. The arrangement of the X-ray sources and the detectors in the numerical experiments is illustrated in Fig. 3. Simultaneous projections from fan beam X-rays through the model system is assumed from limited view angles, e.g., 3-angle projections shown in Fig. 3. Since the model system is known, the projection data can be directly calculated based on the geometry and the hardware properties. If more projection angles are needed, the X-ray sources and corresponding detectors will be arranged around the object evenly. The fast XCT is assumed to be achievable if the transient cross-sectional image could be reconstructed with the simultaneous projection data (or projection data in a very fast sequence) from the limited views. In this way, the sampling rate of the CT (in frame per second) is dependent on the re-start frequency of either the X-ray source or the data collection speed.
Fig. 3. Arrangement of X-ray sources and detectors in the simulation.
Provided that the column is filled with water and air bubbles, the projection data are achieved by calculating the accumulated penetration area of water by each X-ray. The width of each X-ray is considered as 1 mm due to the size of the X-ray source. Corresponding to each fan beam X-ray source, there are 196 pixels in the linear array of each detector. The number of pixels in each detector is an assumption based on the current X-ray source and the detector system in our laboratory for the following validation using experimental data. It should be noted that the up-to-date detectors have much better resolution than this. Hence, better performance can be expected with improved hardware in use. The model system is divided into 128 triangular elements as the start. For better resolution in space, each element is further refined to four elements as shown in Fig. 4. In the following text, 512, 2048 and even 8192 triangular elements are employed in the image reconstruction.
4330
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
Fig. 4. Reconstruction results of bubble objects with ideal/non-ideal shapes in the simulations under different level of noise.
Table 1 Parameters used in the GA method for image reconstruction
3.2. Image reconstruction for bubble objects with ideal shapes
Parameters
Value
Number of population Population size Crossover rate Maximum number of generations Number of grids Number of nodes
20 30 0.8 10,000 128 (mesh unrefined) 77 (mesh unrefined)
The bubble objects with ideal shapes mean that the objects can be perfectly described by the grids, e.g., the triangular ones in this work. In this case, the error in the image reconstruction should approach to zero ideally. As shown in Fig. 4a, four such bubbles are placed in a round area randomly. Using three angles of fan beam X-ray projection data, the image can be reconstructed exactly the same as the original one even if there is a certain amount of noise (up to 40%) superimposed on the projection data. This also demonstrates the robustness of the GA method, not sensitive to the noise.
In each numerical simulation, the initial population for GA iterations is generated in a random manner. The parameters adopted in the GA method are listed in Table 1. These parameters are determined on the basis of numerous numerical experiments in the simulation of various phantoms and reconstruction geometry for faster convergence and better quality image reconstruction as well as the consideration of memory storage requirement.
3.3. Image reconstruction for bubble objects with non-ideal shapes In general, a real bubble object has smooth boundary around instead of the ideally described ones by the triangular grids. Accordingly, the discretization error may cause the defect in the
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4331
Fig. 5. Comparison of reconstruction results under different number of view angles using GA-XCT/FBP method in the simulations.
image reconstruction, especially for the ill-posed mathematical problem in this work. Similar to the arrangement of bubbles in Fig. 4a and b plots the original image with four round bubbles and the corresponding results from image reconstruction based on the three-angle projection data without noise. When the total elements are 128, the bubbles are captured almost accurately both in location and in size, but the image reconstruction cannot go deeper to smooth the bubbles. Further refining the grids to 512 and to 2048, the bubbles are finely characterized with better and better resolution in space. This procedure can be conducted endlessly as long as the computation time is not the limit. This means that for such kind of cases, the resolution in space can reach very high, though the problem in mathematics is still ill-posed. Fig. 4c shows the reconstructed results from projections with random noise of 20% level. As can be seen in Fig. 4c, the image reconstructed at 128 grids is rather coarse compared with the original one. However, using this as the initial guess, the refining operation on the grids results in much better images reconstructed at finer resolution, similar to the above discussion. This further demonstrates the robustness of the proposed GAXCT method to the noise. It is clearly shown from Figs. 4b and c that the attractive feature of “zooming operation” allows the improvement of the overall performance of image reconstruction, especially of the spatial resolution. The effect of noise on the reconstructed results is not apparent even under the noise level of 20%. Generally speaking, the solution to the algebraic equations depends on the good matching between the number of equations and the unknown variables. In the present problem, the number of the unknowns is more than the number of the equations, that is, the so-called ill-posed problem. An optimization problem is thus defined here instead of a deterministic problem in the conventional CT. When applying the “zooming operation”, the number of equations is not changed for the fixed measurement hardware, but the number of unknowns seems to increase a lot. While in fact, when a better initial guess to the image reconstruction is achieved at a coarser scale, each refining operation
Fig. 6. Reconstruction results under three view angles at 8192 grids using GA-XCT in the simulations.
only generates the uncertainty to some extent, where most of the “solutions” in pixels are already the good solutions to the problem. In other words, a good initial guess may probably reduce the illness of the underlying problem, guiding the optimization to the right places, i.e., around the bubble objects to be reconstructed in more detail. 3.4. Image reconstruction under different number of view angles It is very important to determine the minimum number of the view angles for reconstructing the specific investigated objects with the expected resolution. Smaller minimum number
4332
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
Fig. 7. Reconstruction results for ten bubble objects at different grids using GA-XCT in the simulations.
of view angles means greater reduction of the requirement of projections, lower cost of CT hardware, and faster scanning speed for one slice reconstruction. A series of tests are carried out using 2048 grids for different number of view angles, 3, 4, 6, 8, 12 and 24, respectively. Four bubbles with different sizes are randomly arranged, which is also used as the example for the followed validation using real experiments. Fig. 5 shows the reconstructed results by GA-XCT method, with the comparison to the ones by the FBP method that is well acknowledged in the commercial CT algorithm. In the results by present method, the positions of bubbles are well confirmed, and the shapes are well described even under the situation of just three view angles. Based on the three-angle reconstruction result, much better result can be obtained when refining the grids to 8192 (shown in Fig. 6). With more view angles, the results become better and better. As the comparison, the results by FBP method show different performances. Even in the situation of eight view angles, the bubbles cannot be identified. Given more angles data, the bubbles can be recognized in the images with the help of human mind. In mathematics, the bubbles reconstructed still cannot be defined even for the best case using 24-angle projections. This comparison clearly shows that the proposed GA-XCT method is superior to the conventional CT algorithm when dealing with the specific applications. In other words, the new method has potential to reconstruct the phase distributions in the binary systems, and serves as a good solution for the LDT problem. 3.5. Image reconstruction for more complex system Fig. 7 illustrates the image reconstruction of ten bubble objects on a plane at different grids, showing a more complex system than the above ones. For such cases, finer resolution in space is needed to be able to “recognize” the bubbles. At these simulations, at least four view angles are needed to well reveal the positions and sizes of the ten bubbles in the reconstructed images. On one hand, this further demonstrates the strong capability of current GA-XCT method in handling more complex systems. On the other hand, there is still much work to do to
Fig. 8. Reconstruction results for bubble objects with different shapes.
fully understand the know-how about the new method to the X-ray CT applications, especially in the mathematical foundation and the physics. 3.6. Image reconstruction for bubble objects with different shapes All the above validations use round bubbles. While in reality, gas bubbles in liquids may show different shapes. Generally, the bubbles moving up freely maintain three types of shapes: sphericity, spheroidicity and spherical-cap, with observed shapes in a 2D plane (i.e., a cross-section of a column) as circle, ellipse and ring. To show the dependence of the current GA-XCT method on the object shapes, Fig. 8 further discusses the image reconstruction of the bubble objects with different shapes. The results do not show any dependence of the newly proposed method on the bubble shapes. 3.7. Computation time for the image reconstruction Fig. 9 shows the computation time for 1000 generations iteration varying with the number of grids and the number of
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
populations, using a PC with a CPU of Pentium 2.0 GHz. It can be found clearly that the computation time increases linearly with the complexity in the number of grids and in the
4333
population size. The linearity of the present algorithm is favorable to reconstruct the images under finer meshes with not much increase of the computation time. In particular, this algorithm can be easily paralleled to speed up the calculation. 4. Validation by experiments In practical measurements, the projection data would be different from the numerical ones mentioned above due to the existence of the systematic and random errors together. Following the above validation in simulations, this part is to testify the GA-XCT method using the measured projection data in real experiments. 4.1. Experimental setup
Fig. 9. Variation of computation time with the number of grids and parallel populations.
The measurement system consists of a Model 150-I Mobile X-ray Unit operated at the potential of 150 kV, a MiniX-1 X-ray Flat Plate Detector, and a data acquisition system. The experimental model is built by a perspex vessel of 90 mm i.d. and four test-tubes inserted inside. The test-tubes are arranged in the vessel as the bubbles and the positions are the same as drawn in Fig. 5. The vessel is filled with water, while keeping the test-tubes empty. The static object is rotated every 15◦ so that the total number of the view angles is 24 around the object. At each angle, the projection data are sampled more than 20 times to obtain better averaged results. The background correction is considered carefully in the experimental design to eliminate the effect of the walls of the perspex vessel and the test-tubes. This is the essential treatment on the original data from measurements. Thus the details are not discussed here. 4.2. Image reconstruction based on the measured data
Fig. 10. Comparison of measured projections and expected ones at three view angles of projection data.
Limited by the hardware in our measurement system, there are some errors between the measured projections and the expected ones. Fig. 10 shows an example of the comparison results for three view angles’ projection data. It can be seen that
Fig. 11. Comparison of reconstruction results under different number of view angles using GA-XCT/FBP method in the experiments.
4334
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
the noise is not really a random one. Some positions give larger errors, which may somewhat influence the image reconstruction. Fig. 11 shows the reconstructed results by GA-XCT method based on the measured projection data, with the comparison to the ones by the FBP method. Clearly, the GA-XCT generates much better results than the FBP method, and is more robust to the measurement noise.
5. Conclusions A novel measurement technique named GA-XCT was proposed to implement fast measurements of hydrodynamics in multiphase flows using simultaneous, limited angle projection data from X-rays. The model system in the multiphase flows was firstly simplified to have binary-value feature at the local density field exampled as a gas bubbling system. LDT was employed in order to meet the request of high temporal resolution in the X-ray measurement and accordingly to reduce the cost of X-ray CT facilities. This ill-posed problem for the image reconstruction in mathematics was re-defined to an optimization problem, which is a discrete tomography (DT) method, different from the medical CT algorithms. Strengthened relationship among the neighbor elements in the discretized space was considered in the improved GA method to guide the iterations in the image reconstruction. Large amount of numerical simulations and real experiments demonstrated the feasibility and capability of the GA-XCT method in the application of “bubbles” reconstruction using just a few angles of projection data. The proposed method also showed much better performance than the conventional CT reconstruction algorithm (FBP) when dealing with the LDT. Meanwhile, the GA based image reconstruction was robust to the noise, and not dependent on the bubble shapes. In spite of the successful demonstration of GA-XCT in some case studies, it must be reminded that the method is a newly developed idea, still having large room to be further developed both in the theoretical part and in the practical application. Further work is to develop more appropriate criterion(s) to guide the iterations and then judge the final convergence. In this sense, a multi-criteria optimization methodology would be beneficial as successfully applied in ECT’s work (Warsito and Fan, 2001). Fuzzy logic incorporated GA method may be helpful to enhance the capability of this approach to the application of a gas-solids bubbling system.
Acknowledgments The authors would like to thank the Chinese Ministry of Education for the financial support to Dr. Yi Cheng, i.e., A Foundation for the Author of National Excellent Doctoral Dissertation of PR China (No. 200245). Fruitful discussions with Dr. L.-S. Fan at Ohio State University and Dr. P. Cai at Dow Chemical are also greatly acknowledged.
References Banholzer, W.F., Spiro, C.L., Kosky, P.G., Maylotte, D.H., 1987. Direct imaging of time-averaged flow patterns in a fluidized reactor using X-ray computed tomography. Industrial & Engineering Chemistry Research 26, 763. Bartholomew, R.N., Casagrande, R.M., 1957. Measuring solids concentration in fluidized systems by Gamma-ray absorption. Industrial & Engineering Chemistry 49, 428. Beck, M.S., Williams, R.A., 1996. Process tomography: a European innovation and its application. Measurement Science & Technology 7, 215. Caballero, R., Molina, J., Gómez, T., Luque, M., Torrico, A., 2002. A genetic algorithm to solve an integer goal programming model for the higher education. Workshop on Multiple Objective Metaheuristics, Paris. Chaouki, J., Larachi, F., Dudukovic, M.P., 1997. Noninvasive tomographic and velocimetric monitoring of multiphase flows. Industrial & Engineering Chemistry Research 36, 4476. Du, B., Warsito, W., Fan, L.-S., 2004. ECT studies of the choking phenomenon in a gas–solid circulating fluidized bed. A.I.Ch.E. Journal 50 (7), 1386. Dudukovic, M.P., 2002. Opaque multiphase flows: experiments and modeling. Experimental Thermal and Fluid Science 26, 747. Dyakowski, T., Edwards, R.B., Xie, C.G., Williams, R.A., 1997. Application of capacitance tomography to gas solid flows. Chemical Engineering Science 52, 2099. George, D.L., Torczynski, J.R., Shollenberger, K.A., O’Hern, T.J., Ceccio, S.L., 2000. Quantitative tomographic measurements of opaque multiphase flows. Sandia Report SAND2000-0441, Sandia National Laboratories, Livermore, CA. Grassler, T., Wirth, K.E., 1999. X-ray computed tomography in mechanical engineering—a non-intrusive technique to characterize vertical multiphase flows. Proceedings of the Computerized Tomography for Industrial Applications and Image Processing in Radiology, BB 67-CD, Berlin, p. 203. Grassler, T., Wirth, K.E., 2000. X-ray computer tomography—potential and limitation for the measurement of local solids distribution in circulating fluidized beds. Chemical Engineering Journal 77, 65. Harvel, G.D., Hori, K., Kawanishi, K., Chang, J.S., 1996. Real-time crosssectional averaged void fraction measurements in vertical annulus gasliquid two-phase flow by neutron radiography and X-ray tomography techniques. Nuclear Instruments and Methods in Physics Research Section A 371 (3), 544. Harvel, G.D., Hori, K., Kawanishi, K., Chang, J.S., 1999. Cross-sectional void fraction distribution measurements in a vertical annulus two-phase flow by high speed X-ray computed tomography and real-time neutron radiography techniques. Flow Measurement and Instrumentation 10 (4), 259. Holoboff, J., Kantzas, A., Kalogerakis, N., 1995. Utilization of computerassisted tomography in the determination of the spatial variability of voidage in a gas–solid fluidized bed. In: Fluidization, vol. VIII. Preprints Tours, France, Engineering Foundation, New York, p. 295. Joines, J., Houck, C., 1994. On the use of non-stationary penalty functions to solve constrained optimization problems with genetic algorithms. In: IEEE International Symposium on Evolutionary Computation, Orlando, FL, p. 579. Kai, T., Misawa, M., Takahashi, T., Tiseanu, I., Ichikawa, N., Takada, N., 2000. Application of fast X-ray CT scanner to visualization of bubbles in fluidized bed. Journal of Chemical Engineering of Japan 33 (6), 906. Kak, A.C., Slaney, M., 1988. Principles of Computerized Tomographic Imaging. New York, IEEE Press, p. 275. Kantzas, A., 1994. Computation of holdups in fluidized and trickle beds by computer-assisted tomography. A.I.Ch.E. Journal 40 (7), 1254. Kantzas, A., Kalogeraki, N., 1996. Monitoring the fluidization characteristics of polyolefin resins using X-ray computer assisted tomography scanning. Chemical Engineering Science 51 (10), 1979. Kantzas, A., Wright, I., Kalogerakis, N., 1997. Quantification of channeling in polyethylene resin fluid beds using X-ray computer assisted tomography (CAT). Chemical Engineering Science 52 (13), 2023.
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335 Kumar, S.B., Mosleman, D., Dudukovic, M.P., 1997. Gas-holdup measurements in bubble columns using computed tomography. A.I.Ch.E. Journal 43, 1414. Radon, J., 1917. Ueber die bestimmung von funktionen durch ihre integralwerte laengs gewisser manningfaltigkeiten. Berichte ueber die Verhandlungen der Saechsischen Akademie der Wissenschaften zu Leipzig, Mathematisch–Physische Klasse 69, 262. Rowe, P.N., 1971. Experimental properties of bubbles. In: Davidson, J.F., Harrison, D., (Eds.), Fluidization. Academic Press, New York, p. 121. Seville, J.P.K., Morgan, J.E.P., Clift, R., 1986. Tomographic determination of the voidage structure of gas fluidized beds in the jet region. In: ]stergaard, V.K., SZrensen, A., (Eds.), Fluidization. Engineering Foundation, New York, p. 87. Tian, J., Bao, S.L., Zhou, M.Q., 2003. Medical Image Processing and Analysis: Principles and Methods. Beijing, Publishing House of Electronics Industry. Warsito, W., Fan, L.-S., 2001. Neural network based on multi-criteria optimization image reconstruction technique for imaging two- and threephase flow systems using electrical capacitance tomography. Measurement Science & Technology 12, 2198.
4335
Warsito, W., Fan, L.-S., 2003a. 3D-ECT velocimetry for flow structure quantification of gas–liquid–solid fluidized beds. Canadian Journal of Chemical Engineering 81 (3–4), 875. Warsito, W., Fan, L.-S., 2003b. ECT imaging of three-phase fluidized bed based on three-phase capacitance model. Chemical Engineering Science 58 (3–6), 823. Warsito, W., Fan, L.-S., 2003c. Neural network multi-criteria optimization image reconstruction technique (NN-MOIRT) for linear and non-linear process tomography. Chemical Engineering and Processing 42 (8–9), 663. Warsito, W., Ohkawa, M., Kawata, N., Uchida, S., 1999. Cross-sectional distributions of gas and solid holdups in slurry bubble column investigated by ultrasonic computed tomography. Chemical Engineering Science 54, 4711. Weinstein, H., Shao, M., Wasserzug, L., 1984. Radial solid density variation in a fast fluidized bed. A.I.Ch.E. Symposium Series 80, 117. Werther, J., 1999. Measurement techniques in fluidized beds. Powder Technology 102, 15.
A novel X-ray computed tomography method for fast measurement of multiphase flow Changning Wu, Yi Cheng ∗ , Yulong Ding, Fei Wei, Yong Jin Department of Chemical Engineering, Beijing Key Laboratory of Green Chemical Reaction Engineering and Technology, Tsinghua University, Beijing 100084, PR China Received 12 May 2006; received in revised form 18 January 2007; accepted 25 April 2007 Available online 6 May 2007
Abstract A new approach based on an improved genetic algorithm (GA) was proposed to implement the image reconstruction when using X-ray computed tomography (XCT) for the application of fast measurement of multiphase flow dynamics. Instead of directly using a traditional XCT, we pursued to develop a different discrete tomography (DT) method, aiming to achieve a high resolution in time during the measurements with only limited projection data. The proposed method assumed that the interested multiphase flow can be simplified as having distinct dense and dilute phases so that the local phase concentration can be binary-coded, e.g., 0 or 1 in a gas bubbling system. The mathematical problem under these circumstances is strongly ill-posed, and thus tackled with an optimization approach, i.e., a GA incorporated with the underlying physics as some constraints. The numerical simulations mimicking the physical measurements demonstrated the feasibility of the new approach, namely GA-XCT, especially with high robustness to the noise. Experiments were performed to simulate a transient measurement on the gas bubbles in water, with a portable X-ray tube and a 2D plane detector as the hardware and a static object rotating in between. The results further provided the validation of the GA-XCT being superior to the conventional algorithm, e.g., filtered back-projection (FBP) technique, in dealing with the tomography of multiphase system with binary local density field. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Multiphase flow measurement; X-ray computed tomography; Discrete tomography; Genetic algorithm; Limited projection data
1. Introduction Multiphase flow is prevalent in petroleum refining, synthesis gas conversion to fuels and chemicals, bulk commodity chemicals, fine chemicals, manufacture of polymers, biotechnology and pollution control (Dudukovic, 2002). However, the complexity of the underlying dynamics due to the turbulence and nonlinear contact between phases hinders the theoretical understanding on the multiphase flows. Subsequently, the design, control and scaling-up of multiphase flow equipment are still
Abbreviations: ART, algebraic reconstruction technique; DT, discrete tomography; ECT, electrical capacitance tomography; FBP, filtered backprojection; GA, genetic algorithm; -CT, gamma computed tomography; IGP, integer programming problem; LDT, limited data tomography; PT, process tomography; XCT, X-ray computed tomography ∗ Corresponding author. Fax: +86 10 62772051. E-mail address: [email protected] (Y. Cheng). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.04.026
strongly empirical or semi-empirical. Measurement techniques for multiphase flows, especially the non-invasive techniques, always play an important role in promoting the development of multiphase flow theories and engineering applications (Chaouki et al., 1997). With the rapid development of fundamental studies, more stringent demands are brought out for measurement techniques, such as non-invasive measurement with 2D/3D filed information, high resolution in space and/or in time, etc. Under these circumstances, process tomography (PT) finds their unique role in studying the multiphase flow phenomena encountered in chemical engineering fields, involving a number of tomographic techniques, such as X-ray, -ray, optical, ultrasonic, electrical, and nuclear magnetic resonance imaging (Beck and Williams, 1996; Du et al., 2004; Dyakowski et al., 1997; George et al., 2000; Kumar et al., 1997; Rowe, 1971; Seville et al., 1986; Warsito and Fan, 2003a,b,c; Warsito et al., 1999). Among them, X-ray computed tomography (XCT) and electrical capacitance
4326
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
tomography (ECT) are regarded as the two types of the most attractive techniques (Warsito and Fan, 2003c; Werther, 1999). Compared with the ECT technique, XCT (and/or -CT) provides higher spatial resolution up to 1% of a column diameter and even more details at the scale of 200.400 m. Especially, XCT applied in the medical systems has been a mature technique both in theory and in hardware design, also acknowledged as a useful and important technique in the study of multiphase flows. The initial application of the basic principle of CT to the multiphase flow was done by the researchers at Shell Oil Co. Inc. (Bartholomew and Casagrande, 1957). In that work, -ray was used to measure the suspension density in an industrial riser. Different from the current CT method, an empirical density function was firstly assumed, and then the parameters in the function were correlated using a lot of data measured. Similar empirical method to implement a tomography measurement was also used by Weinstein’s group (Weinstein et al., 1984). Medical XCTs started to be applied to study multiphase flows in 1990s (Banholzer et al., 1987; Kantzas, 1994; Kantzas and Kalogeraki, 1996; Kantzas et al., 1997; Holoboff et al., 1995; Grassler and Wirth, 1999, 2000). However, time-averaged measurement results from XCTs were mostly reported in the literature as above. This common feature of XCTs as the low temporal resolution is not applicable in studying the time evolving phenomena of multiphase flows. More detailed researches about applying medical/commercial CT scanners to multiphase flow systems can be found in the review articles (Chaouki et al., 1997; Werther, 1999). Essentially, limited temporal resolution of XCTs comes from a basic request to the image reconstruction for a medical CT measurement, that is, a complete data set for the image reconstruction is inevitably demanded considering the unknown density distribution in human bodies and requested spatial resolution. This consequently causes the rotating action of either X-ray source(s) or the detectors or both of them around the object, which increases the data acquisition time with decreased transient scanning capability for one cross-section imaging. Developed to the fifth generation CT, the hardware has no mechanically moving parts at all, but the electron beams scan around the annular target. In this way, the CT scanner can obtain the data per slice very fast, for instance, in 10 ms (Tian et al., 2003). Regarding the few applications of fast XCT in multiphase flow measurement, the goal can be achieved using more fixed X-ray sources and detectors. Harvel et al. (1996, 1999) built a high speed XCT system to determine the flow regime and the 2D void fraction distribution for a vertical annulus two-phase flow. To increase the time resolution, 18 X-ray sources and 122 CdWO4 detectors were arranged around the test channel. The cross-section sampling time for each slice was 4.0 ms (i.e., 250 frames/s) with a spatial resolution of 2.0 mm. Kai et al. (2000) used eight electron guns and 122 detectors in the X-ray scanner to improve the time resolution. The maximum scan rate of this system could reach 250 slices per second, too. The measurement area was limited to 50 mm in diameter for the requirement of spatial resolution.
For a better understanding of the dynamics in multiphase flows, there is obviously a current need/challenge in improving the response time of the measurement while keeping the high resolution for instantaneous observation. While, compared with a medical XCT, the XCT applied in process engineering is of the limitation in the number of sensors arranged, and consequently in the relative lower spatial resolution due to the limited integral measurement data obtained. The need for high-speed imaging must prompt the number of the sensors to be limited as few as possible, but with sufficient accuracy for process quantification from the limited experimental integral data without sacrificing the time resolution for reliable real time measurements (Warsito and Fan, 2003c). Improving the resolution in time but maintaining the high resolution in space together with not much cost is therefore very challenging to further develop the CT technique for industrial applications. Different from the original idea of an XCT applied in medical purposes, it would be realized that there is often no need to study the gradient change of the phase holdup or the density in many multiphase systems. Distinguishing the dense and dilute phases would be enough. In this sense, some instrumentation/computing power could be saved to increase the process speed, but the resolution would be still high. Such an approach is expected to have significant impact on the multiphase study, even serving to the relevant industry. From this view point, we aim to develop a special fast XCT technique in this work, particularly applied in certain multiphase flow systems, where the local phase distribution shows 0 or 1 feature at this stage. The gas bubbling phenomena in bubble columns and slurry reactors would be the appropriate candidates for the potential applications. Owing to the simplified model system, limited projection data would be sufficient in meeting the request of high spatial resolution and temporal resolution as well. That is to say, limited data tomography (LDT) is concerned in this work, which refers to the situation where ideal tomography is not possible due to insufficient or inappropriately placed parallel-beam/fanbeam projections (for 2D inverse problem). Since the conventional CT algorithms are not applicable when the projection data are very limited, developing a new methodology to efficiently reconstruct the cross-sectional images of bubbling phenomena becomes the major task in this work. 2. Mathematical method: GA-XCT for LDT 2.1. Problem re-definition The principle of the radiation attenuation is the first foundation of XCT. Assuming that the X-rays are monoenergetic from a parallel beam, the residual intensity of X-rays is attenuated by the materials that the X-rays pass through: I = I0 e −
L m ds
,
(1)
where m is the mass attenuation coefficient (m2 /kg) mainly due to the photoelectric, Compton and pair production effects; is the density of the object (kg/m3 ); ds represents infinitesimal
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4327
thin layer of the object (m). The corresponding expression of the projection data is I0 bL = ln (2) = m ds. I L For air–water binary system, the attenuation of X-rays in air can be neglected due to the small density of air compared with water. The accumulated transmission thickness of water by X-ray along line L (Ll ) can be calculated by I0 bL = ln (3) = m ds = m,l l Ll . I L The image reconstruction techniques use the measured projection data as the input to calculate the density distribution of the desired cross-section of the investigated sample as the output. Accordingly, the 2D image of the desired cross-section can be obtained. The mostly applied reconstruction techniques are divided into two categories, the back-projection reconstruction and the iterative reconstruction. As an example of the back-projection reconstruction technique, filtered back-projection (FBP) technique is used in most commercial medical scanners and has proved to be extremely accurate and amenable to fast implementation. This technique can give a rather straightforward intuitive rationale because each projection represents a nearly independent measurement of the object. FBP is applied as the conventional method to be compared with the new approach proposed in this work in the context. Algebraic reconstruction technique (ART) is an entirely different approach for tomographic imaging assuming that the cross-section consists of an array of unknowns and the setting up algebraic equations for the unknowns in terms of the measured projection data (Kak and Slaney, 1988). ART is capable for the situations such as the limited data of projections and the uneven angle-distributed projections. All ART methods are iterative procedures. As shown in Fig. 1, a 2D object is placed between the X-ray source and the detector. The investigated area is discretized by rectangular girds or triangular ones. The mathematical expression of the tomographic problem is given by the equations: b = Ax,
i.e., bk =
M×N
akj xj , k = 1, 2, . . . , K,
(4)
j =1
where K is the total number of the projections, M × N is the total number of pixels across the area, akj means the weighting factor for the contribution of pixel j to the projection element bk . The goal is to determine the unknown image (x) when the experimental projections (b) are available. The weighting factor, akj , is calculated based on the geometrical considerations as the intersection length/area of the kth projection ray with the image pixel j. They depend only on the geometry of the problem and therefore need to be calculated only once. Radon (1917) demonstrated that the mathematical problem in tomography is an inverse problem and that a unique solution is guaranteed provided that data from an infinite number of
Fig. 1. Illustration of the principle of XCT.
views be available. However, for LDT, the problem is ill-posed; we re-defined it to an optimization problem, to find the best solution for the following function: min {Ax − b22 }.
xj ∈R
(5)
For the specific application in this work, xj ∈ {0, 1}, where “0” represents air and “1” represents water. Though the model system is simplified greatly, the resulted nonlinear integer programming problem (IGP) offers a bigger computational complexity, in general with a greater computational cost, which even makes the solution unfeasible in many cases. Caballero et al. (2002) claimed that when solving a complex IGP problem, genetic algorithm (GA) can be an efficient option from the point of view of computational cost as well as the quality of the solution. 2.2. Improved GA GA is a global search and optimization algorithm based on the idea of evolution theory that individuals having a high value of quality will survive to the next generation with a greater probability. In general, GA uses guided random search method, instead of calculating the derivatives or gradients. Therefore, GA shows very robust to the different problems. In a GA approach, an “individual” represents a variable for our problem, i.e., x in formula (5), with a corresponding encoded genome for each individual. Genomes are vectors of genes and the components of a genome are closely related to the parameters of the given problem. A set of genomes is called a population. Genes are sets of real numbers or bits. A GA works with the encoded genome instead of the variables themselves to minimize the objective function, which is often called fitness function, dependent on the GA. The fitness function generates an output of a real value from each input individual, decoded from its corresponding genome. The following operations are defined for a population: the selection
4328
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
operation selects the best genomes from the modified population for the next generation; the crossover operation exchanges a number of genes of two genomes of a population; the mutation operation changes two randomly selected genes of a genome. During each iteration, genomes are evaluated and some are selected as parents based on their fitness value. Pairs of parents are recombined by means of genetic operators and the current population is improved in the sense that some members of the population are replaced by individuals with better fitness values obtained as offspring. The process continues until some stopping condition is satisfied. The algorithm proposed in this work is problem-oriented. The fitness function for GA approach is stated as formula (5), which belongs to the problem of an integer goal programming due to the interested binary value system. A binary coding is thus utilized to represent the parameters of the fitness function. That is to say, each gene, xj , can only be 0 or 1. And, an individual or a genome is just a set of solution to the problem re-defined. A common selection approach assigns a probability of selection, Pj , to each individual j based on its fitness value. Ranking selection operator, only requiring the fitness function to map the solutions to a partially ordered set, is chosen in this paper. Ranking methods assign Pj based on the rank of solution j when all solutions are sorted. Normalized geometric ranking defines Pj for each individual by (Joines and Houck, 1994): P [Selecting the j th individual] = q (1 − q)r−1
(6)
where q is the probability of selecting the best individual, r is the rank of the individual (where 1 is the best), and q = q/(1 − (1 − q)NP ), NP is the population size. The crossover operator takes two individuals and produces two new individuals while the mutation operator alters one individual to produce a single new solution. The application of these two basic types of operators and their derivatives depends on the genome representation used. The one-point crossover works by selecting a common crossover point in each pair of randomly selected two parents and then swapping the corresponding subtrees like standard crossover, with the generation of two genomes of the offspring. If there are two parents S1 = (S11 , S12 , . . . , S1N ) and S2 = (S21 , S22 , . . . , S2N ), their children will be (S11 , . . . , S1m , S2,m+1 , . . . , S2N ) and (S21 , . . . , S2m , S1,m+1 , . . . , S1N ). That is, a head and a tail of an offspring genome will be taken from different parents. For the mentioned problem-oriented nature of the proposed methodology, some prior knowledge has been incorporated into the GA operators. Firstly, a node-based mutation operator is proposed to consider the dependence of the gene mutation rate on the status of its neighbors. The detailed implementation is shown in Fig. 2a. The numbers of “0” and “1” elements are counted around a randomly chosen node. If all the elements around the node are 0s or 1s, the mutation rate will be set to zero, meaning no change for the next iteration/generation; in other cases, the mutation rate of each element, p(0) or p(1) will be assigned according to the proportion of 0s and 1s, which
is given as N (1) ; [N (0) + N (1)]N (0) N (0) , p(1) = [N (0) + N (1)]N (1) p(0) =
(7)
where N (0) and N (1) represent the number of elements 0s and 1s around the selected node, respectively. This treatment imposes the spatial correlations among the status of the neighbored grids, which can be regarded as the continuity assumption of the gas inside each bubble. Compared with the other iterative techniques in the image reconstruction application, the above method incorporates some intelligence into the algorithm. In this way, the GA is no longer a random search method, but performs a directed search to recognize the “bubbles” or maintain the bubble existence, which in turn speeds up the iteration procedure. Secondly, a swapping-based mutation operator is used to increase the speed of convergence. As illustrated in Fig. 2b, to implement a status change from A to C, the iteration procedure of a GA generally steps from A to B then to C. Correspondingly, the fitness function first increases and then decreases, which costs additional computation time. The proposed swappingbased mutation operator allows for the right evaluation on this situation and makes the iteration stepping from A to C directly by simply swapping the corresponding elements. This becomes an auxiliary criterion to help the algorithm being guided. 2.3. Implementation of the GA-XCT method The GA approach to the XCT measurement for multiphase flows is simply named as GA-XCT in this work, which includes four modules in the programming: (1) finite element gridding of the cross-section (i.e., triangular grids); (2) computation of the weighting factor matrix, A, according to the geometry; (3) computation/measurement result of the projections, b; (4) image reconstruction using the improved GA approach until convergence. The GA-XCT algorithm is coded in C + +, in which a GALib developed by MIT researchers (http://lancet.mit.edu/ga) is employed. GA-Lib provides a framework and the basic functions for general applications of GA. The new modifications described above are superimposed on the GA-Lib successfully. Besides, multi-population GA approach is applied to speed up the calculation and convergence of the IGP problem instead of a common GA method. To start the GA-based image reconstruction, an initial population must be provided for GA, which is usually generated in a random manner for the entire population. However, since the algorithm can iteratively improve existing solutions, the beginning population can be seeded with potentially good solutions, with the remainder of the population being randomly generated.
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4329
Fig. 2. Illustration of GA operators and zooming operation: (a) the node-based mutation operator; (b) the swapping-based mutation operator; and (c) the grid refinement.
2.4. Method to improve the spatial resolution Based on the triangle meshing method, the GA approach can be conveniently applied to higher resolution image restoration by refining the grids in the way shown in Fig. 2c. It is simply named as “zooming operation” in this work for easy expression. The solutions under the coarser scale grids are taken as a better initial guess to the problem solutions under the finer scale grids. By inheriting the solutions from the coarser meshes, it can be found in following context that the image reconstruction results under finer resolutions in space can be achieved stably with good accuracy. 3. Validation by numerical simulations 3.1. The test phantom and simulation details To evaluate the performance of the proposed GA-XCT method, numerical simulations are carried out using various simulated distribution functions (i.e., phantoms) as the references. Some (near-) round objects with different sizes are arranged randomly in a bigger round area, representing the situation of the cross-section in a gas–liquid bubbling column. At this stage, we start from few bubbles cases in the measurement plane, such as four bubbles for the simpler case and ten bubbles for the more complicated case. The arrangement of the X-ray sources and the detectors in the numerical experiments is illustrated in Fig. 3. Simultaneous projections from fan beam X-rays through the model system is assumed from limited view angles, e.g., 3-angle projections shown in Fig. 3. Since the model system is known, the projection data can be directly calculated based on the geometry and the hardware properties. If more projection angles are needed, the X-ray sources and corresponding detectors will be arranged around the object evenly. The fast XCT is assumed to be achievable if the transient cross-sectional image could be reconstructed with the simultaneous projection data (or projection data in a very fast sequence) from the limited views. In this way, the sampling rate of the CT (in frame per second) is dependent on the re-start frequency of either the X-ray source or the data collection speed.
Fig. 3. Arrangement of X-ray sources and detectors in the simulation.
Provided that the column is filled with water and air bubbles, the projection data are achieved by calculating the accumulated penetration area of water by each X-ray. The width of each X-ray is considered as 1 mm due to the size of the X-ray source. Corresponding to each fan beam X-ray source, there are 196 pixels in the linear array of each detector. The number of pixels in each detector is an assumption based on the current X-ray source and the detector system in our laboratory for the following validation using experimental data. It should be noted that the up-to-date detectors have much better resolution than this. Hence, better performance can be expected with improved hardware in use. The model system is divided into 128 triangular elements as the start. For better resolution in space, each element is further refined to four elements as shown in Fig. 4. In the following text, 512, 2048 and even 8192 triangular elements are employed in the image reconstruction.
4330
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
Fig. 4. Reconstruction results of bubble objects with ideal/non-ideal shapes in the simulations under different level of noise.
Table 1 Parameters used in the GA method for image reconstruction
3.2. Image reconstruction for bubble objects with ideal shapes
Parameters
Value
Number of population Population size Crossover rate Maximum number of generations Number of grids Number of nodes
20 30 0.8 10,000 128 (mesh unrefined) 77 (mesh unrefined)
The bubble objects with ideal shapes mean that the objects can be perfectly described by the grids, e.g., the triangular ones in this work. In this case, the error in the image reconstruction should approach to zero ideally. As shown in Fig. 4a, four such bubbles are placed in a round area randomly. Using three angles of fan beam X-ray projection data, the image can be reconstructed exactly the same as the original one even if there is a certain amount of noise (up to 40%) superimposed on the projection data. This also demonstrates the robustness of the GA method, not sensitive to the noise.
In each numerical simulation, the initial population for GA iterations is generated in a random manner. The parameters adopted in the GA method are listed in Table 1. These parameters are determined on the basis of numerous numerical experiments in the simulation of various phantoms and reconstruction geometry for faster convergence and better quality image reconstruction as well as the consideration of memory storage requirement.
3.3. Image reconstruction for bubble objects with non-ideal shapes In general, a real bubble object has smooth boundary around instead of the ideally described ones by the triangular grids. Accordingly, the discretization error may cause the defect in the
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
4331
Fig. 5. Comparison of reconstruction results under different number of view angles using GA-XCT/FBP method in the simulations.
image reconstruction, especially for the ill-posed mathematical problem in this work. Similar to the arrangement of bubbles in Fig. 4a and b plots the original image with four round bubbles and the corresponding results from image reconstruction based on the three-angle projection data without noise. When the total elements are 128, the bubbles are captured almost accurately both in location and in size, but the image reconstruction cannot go deeper to smooth the bubbles. Further refining the grids to 512 and to 2048, the bubbles are finely characterized with better and better resolution in space. This procedure can be conducted endlessly as long as the computation time is not the limit. This means that for such kind of cases, the resolution in space can reach very high, though the problem in mathematics is still ill-posed. Fig. 4c shows the reconstructed results from projections with random noise of 20% level. As can be seen in Fig. 4c, the image reconstructed at 128 grids is rather coarse compared with the original one. However, using this as the initial guess, the refining operation on the grids results in much better images reconstructed at finer resolution, similar to the above discussion. This further demonstrates the robustness of the proposed GAXCT method to the noise. It is clearly shown from Figs. 4b and c that the attractive feature of “zooming operation” allows the improvement of the overall performance of image reconstruction, especially of the spatial resolution. The effect of noise on the reconstructed results is not apparent even under the noise level of 20%. Generally speaking, the solution to the algebraic equations depends on the good matching between the number of equations and the unknown variables. In the present problem, the number of the unknowns is more than the number of the equations, that is, the so-called ill-posed problem. An optimization problem is thus defined here instead of a deterministic problem in the conventional CT. When applying the “zooming operation”, the number of equations is not changed for the fixed measurement hardware, but the number of unknowns seems to increase a lot. While in fact, when a better initial guess to the image reconstruction is achieved at a coarser scale, each refining operation
Fig. 6. Reconstruction results under three view angles at 8192 grids using GA-XCT in the simulations.
only generates the uncertainty to some extent, where most of the “solutions” in pixels are already the good solutions to the problem. In other words, a good initial guess may probably reduce the illness of the underlying problem, guiding the optimization to the right places, i.e., around the bubble objects to be reconstructed in more detail. 3.4. Image reconstruction under different number of view angles It is very important to determine the minimum number of the view angles for reconstructing the specific investigated objects with the expected resolution. Smaller minimum number
4332
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
Fig. 7. Reconstruction results for ten bubble objects at different grids using GA-XCT in the simulations.
of view angles means greater reduction of the requirement of projections, lower cost of CT hardware, and faster scanning speed for one slice reconstruction. A series of tests are carried out using 2048 grids for different number of view angles, 3, 4, 6, 8, 12 and 24, respectively. Four bubbles with different sizes are randomly arranged, which is also used as the example for the followed validation using real experiments. Fig. 5 shows the reconstructed results by GA-XCT method, with the comparison to the ones by the FBP method that is well acknowledged in the commercial CT algorithm. In the results by present method, the positions of bubbles are well confirmed, and the shapes are well described even under the situation of just three view angles. Based on the three-angle reconstruction result, much better result can be obtained when refining the grids to 8192 (shown in Fig. 6). With more view angles, the results become better and better. As the comparison, the results by FBP method show different performances. Even in the situation of eight view angles, the bubbles cannot be identified. Given more angles data, the bubbles can be recognized in the images with the help of human mind. In mathematics, the bubbles reconstructed still cannot be defined even for the best case using 24-angle projections. This comparison clearly shows that the proposed GA-XCT method is superior to the conventional CT algorithm when dealing with the specific applications. In other words, the new method has potential to reconstruct the phase distributions in the binary systems, and serves as a good solution for the LDT problem. 3.5. Image reconstruction for more complex system Fig. 7 illustrates the image reconstruction of ten bubble objects on a plane at different grids, showing a more complex system than the above ones. For such cases, finer resolution in space is needed to be able to “recognize” the bubbles. At these simulations, at least four view angles are needed to well reveal the positions and sizes of the ten bubbles in the reconstructed images. On one hand, this further demonstrates the strong capability of current GA-XCT method in handling more complex systems. On the other hand, there is still much work to do to
Fig. 8. Reconstruction results for bubble objects with different shapes.
fully understand the know-how about the new method to the X-ray CT applications, especially in the mathematical foundation and the physics. 3.6. Image reconstruction for bubble objects with different shapes All the above validations use round bubbles. While in reality, gas bubbles in liquids may show different shapes. Generally, the bubbles moving up freely maintain three types of shapes: sphericity, spheroidicity and spherical-cap, with observed shapes in a 2D plane (i.e., a cross-section of a column) as circle, ellipse and ring. To show the dependence of the current GA-XCT method on the object shapes, Fig. 8 further discusses the image reconstruction of the bubble objects with different shapes. The results do not show any dependence of the newly proposed method on the bubble shapes. 3.7. Computation time for the image reconstruction Fig. 9 shows the computation time for 1000 generations iteration varying with the number of grids and the number of
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
populations, using a PC with a CPU of Pentium 2.0 GHz. It can be found clearly that the computation time increases linearly with the complexity in the number of grids and in the
4333
population size. The linearity of the present algorithm is favorable to reconstruct the images under finer meshes with not much increase of the computation time. In particular, this algorithm can be easily paralleled to speed up the calculation. 4. Validation by experiments In practical measurements, the projection data would be different from the numerical ones mentioned above due to the existence of the systematic and random errors together. Following the above validation in simulations, this part is to testify the GA-XCT method using the measured projection data in real experiments. 4.1. Experimental setup
Fig. 9. Variation of computation time with the number of grids and parallel populations.
The measurement system consists of a Model 150-I Mobile X-ray Unit operated at the potential of 150 kV, a MiniX-1 X-ray Flat Plate Detector, and a data acquisition system. The experimental model is built by a perspex vessel of 90 mm i.d. and four test-tubes inserted inside. The test-tubes are arranged in the vessel as the bubbles and the positions are the same as drawn in Fig. 5. The vessel is filled with water, while keeping the test-tubes empty. The static object is rotated every 15◦ so that the total number of the view angles is 24 around the object. At each angle, the projection data are sampled more than 20 times to obtain better averaged results. The background correction is considered carefully in the experimental design to eliminate the effect of the walls of the perspex vessel and the test-tubes. This is the essential treatment on the original data from measurements. Thus the details are not discussed here. 4.2. Image reconstruction based on the measured data
Fig. 10. Comparison of measured projections and expected ones at three view angles of projection data.
Limited by the hardware in our measurement system, there are some errors between the measured projections and the expected ones. Fig. 10 shows an example of the comparison results for three view angles’ projection data. It can be seen that
Fig. 11. Comparison of reconstruction results under different number of view angles using GA-XCT/FBP method in the experiments.
4334
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335
the noise is not really a random one. Some positions give larger errors, which may somewhat influence the image reconstruction. Fig. 11 shows the reconstructed results by GA-XCT method based on the measured projection data, with the comparison to the ones by the FBP method. Clearly, the GA-XCT generates much better results than the FBP method, and is more robust to the measurement noise.
5. Conclusions A novel measurement technique named GA-XCT was proposed to implement fast measurements of hydrodynamics in multiphase flows using simultaneous, limited angle projection data from X-rays. The model system in the multiphase flows was firstly simplified to have binary-value feature at the local density field exampled as a gas bubbling system. LDT was employed in order to meet the request of high temporal resolution in the X-ray measurement and accordingly to reduce the cost of X-ray CT facilities. This ill-posed problem for the image reconstruction in mathematics was re-defined to an optimization problem, which is a discrete tomography (DT) method, different from the medical CT algorithms. Strengthened relationship among the neighbor elements in the discretized space was considered in the improved GA method to guide the iterations in the image reconstruction. Large amount of numerical simulations and real experiments demonstrated the feasibility and capability of the GA-XCT method in the application of “bubbles” reconstruction using just a few angles of projection data. The proposed method also showed much better performance than the conventional CT reconstruction algorithm (FBP) when dealing with the LDT. Meanwhile, the GA based image reconstruction was robust to the noise, and not dependent on the bubble shapes. In spite of the successful demonstration of GA-XCT in some case studies, it must be reminded that the method is a newly developed idea, still having large room to be further developed both in the theoretical part and in the practical application. Further work is to develop more appropriate criterion(s) to guide the iterations and then judge the final convergence. In this sense, a multi-criteria optimization methodology would be beneficial as successfully applied in ECT’s work (Warsito and Fan, 2001). Fuzzy logic incorporated GA method may be helpful to enhance the capability of this approach to the application of a gas-solids bubbling system.
Acknowledgments The authors would like to thank the Chinese Ministry of Education for the financial support to Dr. Yi Cheng, i.e., A Foundation for the Author of National Excellent Doctoral Dissertation of PR China (No. 200245). Fruitful discussions with Dr. L.-S. Fan at Ohio State University and Dr. P. Cai at Dow Chemical are also greatly acknowledged.
References Banholzer, W.F., Spiro, C.L., Kosky, P.G., Maylotte, D.H., 1987. Direct imaging of time-averaged flow patterns in a fluidized reactor using X-ray computed tomography. Industrial & Engineering Chemistry Research 26, 763. Bartholomew, R.N., Casagrande, R.M., 1957. Measuring solids concentration in fluidized systems by Gamma-ray absorption. Industrial & Engineering Chemistry 49, 428. Beck, M.S., Williams, R.A., 1996. Process tomography: a European innovation and its application. Measurement Science & Technology 7, 215. Caballero, R., Molina, J., Gómez, T., Luque, M., Torrico, A., 2002. A genetic algorithm to solve an integer goal programming model for the higher education. Workshop on Multiple Objective Metaheuristics, Paris. Chaouki, J., Larachi, F., Dudukovic, M.P., 1997. Noninvasive tomographic and velocimetric monitoring of multiphase flows. Industrial & Engineering Chemistry Research 36, 4476. Du, B., Warsito, W., Fan, L.-S., 2004. ECT studies of the choking phenomenon in a gas–solid circulating fluidized bed. A.I.Ch.E. Journal 50 (7), 1386. Dudukovic, M.P., 2002. Opaque multiphase flows: experiments and modeling. Experimental Thermal and Fluid Science 26, 747. Dyakowski, T., Edwards, R.B., Xie, C.G., Williams, R.A., 1997. Application of capacitance tomography to gas solid flows. Chemical Engineering Science 52, 2099. George, D.L., Torczynski, J.R., Shollenberger, K.A., O’Hern, T.J., Ceccio, S.L., 2000. Quantitative tomographic measurements of opaque multiphase flows. Sandia Report SAND2000-0441, Sandia National Laboratories, Livermore, CA. Grassler, T., Wirth, K.E., 1999. X-ray computed tomography in mechanical engineering—a non-intrusive technique to characterize vertical multiphase flows. Proceedings of the Computerized Tomography for Industrial Applications and Image Processing in Radiology, BB 67-CD, Berlin, p. 203. Grassler, T., Wirth, K.E., 2000. X-ray computer tomography—potential and limitation for the measurement of local solids distribution in circulating fluidized beds. Chemical Engineering Journal 77, 65. Harvel, G.D., Hori, K., Kawanishi, K., Chang, J.S., 1996. Real-time crosssectional averaged void fraction measurements in vertical annulus gasliquid two-phase flow by neutron radiography and X-ray tomography techniques. Nuclear Instruments and Methods in Physics Research Section A 371 (3), 544. Harvel, G.D., Hori, K., Kawanishi, K., Chang, J.S., 1999. Cross-sectional void fraction distribution measurements in a vertical annulus two-phase flow by high speed X-ray computed tomography and real-time neutron radiography techniques. Flow Measurement and Instrumentation 10 (4), 259. Holoboff, J., Kantzas, A., Kalogerakis, N., 1995. Utilization of computerassisted tomography in the determination of the spatial variability of voidage in a gas–solid fluidized bed. In: Fluidization, vol. VIII. Preprints Tours, France, Engineering Foundation, New York, p. 295. Joines, J., Houck, C., 1994. On the use of non-stationary penalty functions to solve constrained optimization problems with genetic algorithms. In: IEEE International Symposium on Evolutionary Computation, Orlando, FL, p. 579. Kai, T., Misawa, M., Takahashi, T., Tiseanu, I., Ichikawa, N., Takada, N., 2000. Application of fast X-ray CT scanner to visualization of bubbles in fluidized bed. Journal of Chemical Engineering of Japan 33 (6), 906. Kak, A.C., Slaney, M., 1988. Principles of Computerized Tomographic Imaging. New York, IEEE Press, p. 275. Kantzas, A., 1994. Computation of holdups in fluidized and trickle beds by computer-assisted tomography. A.I.Ch.E. Journal 40 (7), 1254. Kantzas, A., Kalogeraki, N., 1996. Monitoring the fluidization characteristics of polyolefin resins using X-ray computer assisted tomography scanning. Chemical Engineering Science 51 (10), 1979. Kantzas, A., Wright, I., Kalogerakis, N., 1997. Quantification of channeling in polyethylene resin fluid beds using X-ray computer assisted tomography (CAT). Chemical Engineering Science 52 (13), 2023.
C. Wu et al. / Chemical Engineering Science 62 (2007) 4325 – 4335 Kumar, S.B., Mosleman, D., Dudukovic, M.P., 1997. Gas-holdup measurements in bubble columns using computed tomography. A.I.Ch.E. Journal 43, 1414. Radon, J., 1917. Ueber die bestimmung von funktionen durch ihre integralwerte laengs gewisser manningfaltigkeiten. Berichte ueber die Verhandlungen der Saechsischen Akademie der Wissenschaften zu Leipzig, Mathematisch–Physische Klasse 69, 262. Rowe, P.N., 1971. Experimental properties of bubbles. In: Davidson, J.F., Harrison, D., (Eds.), Fluidization. Academic Press, New York, p. 121. Seville, J.P.K., Morgan, J.E.P., Clift, R., 1986. Tomographic determination of the voidage structure of gas fluidized beds in the jet region. In: ]stergaard, V.K., SZrensen, A., (Eds.), Fluidization. Engineering Foundation, New York, p. 87. Tian, J., Bao, S.L., Zhou, M.Q., 2003. Medical Image Processing and Analysis: Principles and Methods. Beijing, Publishing House of Electronics Industry. Warsito, W., Fan, L.-S., 2001. Neural network based on multi-criteria optimization image reconstruction technique for imaging two- and threephase flow systems using electrical capacitance tomography. Measurement Science & Technology 12, 2198.
4335
Warsito, W., Fan, L.-S., 2003a. 3D-ECT velocimetry for flow structure quantification of gas–liquid–solid fluidized beds. Canadian Journal of Chemical Engineering 81 (3–4), 875. Warsito, W., Fan, L.-S., 2003b. ECT imaging of three-phase fluidized bed based on three-phase capacitance model. Chemical Engineering Science 58 (3–6), 823. Warsito, W., Fan, L.-S., 2003c. Neural network multi-criteria optimization image reconstruction technique (NN-MOIRT) for linear and non-linear process tomography. Chemical Engineering and Processing 42 (8–9), 663. Warsito, W., Ohkawa, M., Kawata, N., Uchida, S., 1999. Cross-sectional distributions of gas and solid holdups in slurry bubble column investigated by ultrasonic computed tomography. Chemical Engineering Science 54, 4711. Weinstein, H., Shao, M., Wasserzug, L., 1984. Radial solid density variation in a fast fluidized bed. A.I.Ch.E. Symposium Series 80, 117. Werther, J., 1999. Measurement techniques in fluidized beds. Powder Technology 102, 15.