The High Efficiency fuel Channel (HEC) consists of a pressure tube, a ceramic insulator, a liner tube, and fuel bundles. Figure 5 shows a 3-D view of HEC. The outer surface of the pressure tube is exposed to a moderator. The moderator could be a liquid moderator such as heavy-water or a solid moderator. The purpose of using an insulator is to reduce the operating temperature of the pressure tube and heat losses from the coolant to the moderator. Low operating temperatures of the pressure tube would allow for the use of available materials such as Zr-2.5%Nb, which has low absorption cross-sections for thermal neutrons (Chow and Khartabil, 2008).

The proposed material for the ceramic insulator is Yttria Stabilized Zirconia (YSZ) (Chow and Khartabil, 2008). YSZ has a low neutron absorption cross-section, low thermal — conductivity and high corrosion resistance in exposure to water at supercritical conditions (Chow and Khartabil, 2008). These properties make YSZ a good candidate as an insulator. The liner, which is a perforated tube and made of stainless steel, intends to protect the ceramic insulator from being damaged during operation or possible refuelling due to stresses introduced by fuel bundles and from erosion by the coolant flow.

In this simulation the input parameter value (H) is attributed to arbitrary Keff. So if the favorite Kgff is: 1 then the value of H will be defined: 100 and it is also for Set point (reference Keff). So Set Point: 50 means the reference K, ff is: 0.50 and in this situation this K, ff is enumerated as the arbitrary and favorite Keff that stability of reactor in this situation is based on it. This arbitrary Keff with output of Zero-Order Hold block is compared for performance of simulation by SIMULINK software. Also for example velocity: 3 means the velocity of control rod in this simulation is: 3 units per second (for example: 3mm/ s). In this state the velocity of control rod is increased comparing to the last stage which was: 1 unit per second. The velocity of control rod belongs to Speed block which is an input to Tsp Sum block in the block diagram.

By changing the values of the cited parameters (which were: Set Point, the velocity of control rod (v), H and the delay time), the different states of the graphs can be shown according to Figs.4, 5, 6 and 7:

The Figs.4, 5, 6 and 7 show if the velocity of control rod for upward and downward movements is increased then the Keff will be pendulous surrounds of defined Set Point (reference Keff) and also the oscillation amplitude will be more than lower velocity situations. For low velocities the oscillation amplitude is slight and acceptable. Another effective factor is inherent delay (such as derived delay of control rod mechanism) and its inordinate increasing can cause unstable states.

3. Conclusion

By this simulation the best response and operation which a reactor can have from stability aspect, according to its control rod velocity is derived.

According to Figs.4-8 the best status in the Fig.4 is observed in which there is no vibration in response. In this simulation the stability of reactor depends on either velocity or delay time values directly because delay time plays a key role. Therefore in this simulation the admissible ranges of velocity and delay time which can be caused to stable the reactor are respectively: low velocity of control rod around 1-3 units per second and short delay time (10ms). However in this case reach critical state (Keff=1) for nuclear reactor will be taken more than modes Figs.5-8. As in all the figures is observed, so can deduce the velocity of control rod plays the more important role than delay time of mechanism in stability of nuclear reactor. Whereas for minor and major changes of reactivity and shut down of reactor in emergency states, there are some kinds of control rods at nuclear reactors cores such as regulating rods, safety rods and shim rods, therefore this simulation can be applied for each control rod in either LWR nuclear reactor or research reactor cores which have vertical control rods. Also this simulation can be applied for each batch control rods which

act in the same way as the cluster at fuel assembly at core of nuclear power reactors though

the moving speed of regulating rod is much less than the safety rod.

4. Acknowledgment

This book chapter is related to a research project entitled: "The Simulation of a Model by

SIMULINK of MATLAB for Determining the Best Ranges for Velocity and Delay Time of

Control Rod Movement in LWR Reactors" that by financial supporting the Islamic Azad

University-South Tehran Branch has been carried out.

Kinetic Monte Carlo models and simulations have been employed to study cascade ageing and defect accumulation at low doses using the input from atomistic simulations of cascades and defect migration properties (Barashev, Bacon et al. 1999; Gao, Bacon et al. 1999; Barashev, Bacon et al. 2000; Heinisch, Singh et al. 2000; Heinisch and Singh 2003; Domain, Becquart et al. 2004; Becquart, Domain et al. 2005; Caturla, Soneda et al. 2006). These have mostly focused on intrinsic defect (interstitial and vacancy) migration, clustering and annihilation under irradiation conditions. In this example, we focus on extrinsic gas atoms such as helium or hydrogen migration in irradiated materials.

For the Fe-He system, modeling of helium clusters has been performed by Morishita, Wirth and co-workers (Morishita, Sugano et al. 2003a; Morishita, Sugano et al. 2003b) and Bringa, Wirth et al. (2003)using semi-empirical potentials. In one paper, Morishita, Sugano et al. (2003a) performed molecular dynamics (MD) calculations to evaluate the thermal stability of helium-vacancy clusters in Fe In another paper, Morishita, Sugano et al. (2003b)have looked at dissolution of helium-vacancy clusters as a function of temperature increase using the empirical potentials for the Fe-He system. Wirth and Bringa (Wirth and Bringa 2004) have simulated the motion of one single 2He-3Vac cluster at 1000 K using the same potential system. First principles calculations of helium atoms in interstitial and substitutional sites has been performed by Fu and Willaime 2005 and can be used to provide an input parameter set for kinetic Monte Carlo calculations.

The kinetic Monte Carlo model(Deo, Srivilliputhur et al. 2006; Deo, Okuniewski et al. 2007; Deo, Srinivasan et al. 2007) consists of helium interstitials on the octahedral sublattice and vacancies on the bcc iron lattice. The migration of the free (not clustered) helium is shown in figure 5(a) and that of the vacancies and self interstitial atoms is shown in Figure 5(b). The
lowest-energy migration path of the SIA corresponds to a nearest-neighbor translation- rotation jump in the <111> direction. The rates of migration of the point defect entities are calculated as

i Л

rmigration V migrationexp

where the superscript i refers to the helium, self interstitial atoms and the vacancy point defect entities. The rate of migration of the point defect entity is rmigration, the attempt frequency is v’migration, the migration barrier is Emigration, while kB and T are the Boltzmann constant and the temperature respectively. Two point defect entities are considered to be in a cluster when the distance between them is less than a0, which is the lattice constant of bcc iron. Interstitial atoms are not allowed to dissociate from clusters. Dissociation of the helium and the vacancy from the cluster is described in Figures 2(c) and 2(d) respectively. The rate
of dissociation of a point defect entity (i = helium or vacancy) from a cluster into the bulk lattice is considered to be thermally activated and is calculated as:

( E Л

ri = v exp dissociation (4)

fdissociation udissociationг І і г-,-, I ^ ‘

L kBT J

where Гdissociation is the rate or dissociation, V’dissocation is the attempt frequencies, Edissociation is the energy of dissociation. The dissociation energy Edissociation of a point defect from a cluster

is taken to be the sum of the energy to bind a point defect entity to the cluster and E‘migration.

Small bubbles migrate with a Arrhenius migration rate parameterized by Table 1. Larger bubbles migrate by surface diffusion at the bubble-matrix interface (Cottrell 2003; Evans 2004),

D, = Ds(f^), (5)

where DB and Ds are the bubble and surface diffusivities respectively, Q is the atomic volume and r the radius of the bubble.

Morishita et al. (Morishita et al. 2003) and Fu et al (Fu and Willaime 2005) have calculated the migration energies of helium and vacancies as well as the binding energies of some helium — vacancy clusters. We employ these migration barriers to calculate the rate of migration of the point defect entities. Parameters used to calculate these quantities are described in Table 1. These parameters are used to calculate the rates migration and dissociation events (Equations 1-3) in the system and build the event catalog for the kMC simulation.

The event catalog is generated by calculating the rates of migration or dissociation of the point defect entities using Equations (3-6). The kMC event catalog consists of the migration, clustering and dissociation of the point defect entities, helium, self interstitial atoms and vacancies. The transition probability of each event is proportional to the rate of event occurrence, calculated by the Equations (3-6). We follow the well established kMC simulation algorithm (Bortz, Kalos et al. 1975; Fichthorn and Weinberg 1991) which is a stochastic, atomic — scale method to simulate the time-evolution of defects and nano/ microstructural evolution that focuses on individual defects and not on atomic vibrations.

Reaction pathways available in the system are tabulated in table 1. Rates of migration of events are calculated at each kMC step. Parameters are obtained both from literature using first principles calculations (Fu, Willaime et al. 2004; Fu and Willaime 2005) and also from molecular statics calculations using semi-empirical potentials as employed by Morishita et al (Morishita et al. 2003) and Wirth and Bringa (Wirth and Bringa 2004) (using the Ackland Finnis-Sinclair potential, the Wilson-Johnson potential and the Ziegler-Biersack-Littmark — Beck potential for describing the interactions of Fe-Fe, Fe-He and He-He, respectively). Helium atoms are introduced at random positions at the beginning of the simulation.

Self interstitial atoms (SIA) are produced in the simulation along with vacancies as Frenkel pairs as cascade debris. Self interstitial atoms are mobile and cluster. Self interstitial clusters up to size 5 migrate one dimensionally. Self interstitials of higher size are stationary. Mono — and clustered vacancies are mobile. Vacancies can dissociate from a vacancy cluster as well as a helium-vacancy cluster. Helium atoms migrate on the octahedral sublattice as well as part of helium-vacancy bubbles. A substitutional helium is considered as a 1-1 He1V1 and is mobile. If a bubble has a helium-vacancy ratio greater than 5, it emits a self interstitial atom. Small bubbles migrate according to Equation 3 and Table 1 while large bubbles migrate according to Eq. 3. Self interstitial atoms and vacancies annihilate when they meet either as point defects or in a cluster. The boundary of the simulation cell acts as a sink for the point defect entities.

At each kMC step, the system is monitored to identify a clustering event. When two point defect entities (helium-helium, vacancy-vacancy, helium-vacancy) are in a cluster the simulation creates a mapping between the entities and the cluster such that for each cluster there are at least two entities associated with the cluster. The event catalog is updated with the new rates of event occurrence and the transition probabilities for the next kMC event are calculated using Equations (4-6) using the parameters from Table 1.

Simulations were performed for a damage energy of 100keV. A production rate of randomly distributed cascades is assumed such that the damage is introduced at a rate of 10-6 dpa/ s. The simulation cell size is 400aox400a0x400a0/ where a0 is the lattice parameter of iron. The number of Frenkel pairs introduced in the simulation cell are calculated by the Norgett- Robinson-Torrens relationship,

displacements per cascade = (6)

where Ed is the energy of incident particles while Ed is the threshold energy (40eV for iron (1994)). Initial concentration of helium is a parameter in the simulation and is varied from 1 to 25 appm/ dpa. Damage measured in dpa is another parameter in the system and determines the length of the simulation. The incident energy is also a simulation control parameter that determines the number of defects introduced in a single cascade. The simulation is performed until sufficient damage is accumulated in the simulation cell.

Simulations are performed for values of the He concentrations varying from 1 to 25 appm/ dpa (appm = atomic parts per million). Initial concentrations of both point defect entities are calculated from the NRT formula (Equation 4) using an incident energy of 100keV. Overall damage is introduced at the rate of 10-6 dpa/ s and the simulation is carried out until the required amount of damage has accumulated in the system.

Figure 3 shows the concentration of bubbles as a function of He concentration (appm/ dpa) after damage equal to 0.1 and 1 dpa is introduced. The simulation temperature is 0.3Tm. The bubble density increases with increasing He/dpa ratio. The bubble density can be described by a power law expression,

CB = K (cHe ) , CB = K (cHe) (7)

where is the bubble density, is the helium concentration expressed as appm He/dpa and K and m are constants determined by the kinetic Monte Carlo simulations. We find that the exponent m is approximately 0.5. In Figure 6, the solid lines are fit to the data assuming a square root dependence of the bubble density on helium concentration. Thus the bubble density increases as the square root of the He/dpa ratio. While experimental evidence of this variation is difficult to find, such dependence has been suggested by rate theory calculations (SINGH and TRINKAUS 1992) for the case of cold helium implantation annealing. The square root dependence is also found for equilibrium bubbles containing an ideal gas(Marksworth 1973).

The bubble density increases with damage (expressed as dpa) as damage is accumulated till 1 dpa. At higher dpa ratios more vacancies are produced that may serve as nucleation sites for embryonic bubbles by trapping helium. Thus the bubble density scales directly with increasing damage. That bubble density increases with accumulating dpa is observed in experiments (Sencer, Bond et al 2001; Sencer, Garner et al. 2002; Zinkle, Hashimoto et al. 2003). These experiments suggest much smaller values of bubble density than those calculated in the present simulations. Embryonic bubbles are submicroscopic and difficult to estimate from experimental observations; while the kMC simulations have a large fraction of nanometer size bubbles; making comparison with experiment difficult.

This example demonstrates the use of the kinetic Monte Carlo method to simulate defect diffusion in irradiated materials. The method depends on the ability to calculate rates of migration events of individual defects. The advantage of the method is that all atoms do not necessarily need to be simulated and large time scales are accessible to the simulation provided the underlying physics remains invariant.

In the 99Mo production system with the solution irradiation method, a flowing target solution with a high concentration is in contact with the structural material of the capsule and the pipes, and then it is important to investigate compatibility between the flowing target solution and the structural material. In the previous tests (Inaba et al., 2009), the circulating solution test was carried out under y-ray irradiation. However, the SUS304 specimen used in the test was only immersed in the bottom of an irradiation container with a volume of 2,000 cm3, and the specimen had no influence of the circulating solution flow, and then the compatibility was not cleared. Therefore, the compatibility test between the flowing target solution and the structural material was carried out, and the corrosivity of the flowing target solution for the structural material as well as the chemical stability of the solution was investigated.

An aqueous K2MoO4 solution, which was the first candidate of the irradiation target, was used in the test. The purity of K2MoO4 used in the test was over 98%.

The temperature profile of the coolant along the heated length of the fuel channel can be calculated based on the heat balance. Equation (16) was used to calculate the temperature profile of the coolant. The NIST REPFROP software Version 8.0 was used to determine the thermophysical properties at a bulk-fluid temperature corresponding to each one-millimeter interval.

hi+1 = hj + ЕЛ*.. Дх (16)

In Eq. (16), qx is the axial heat flux value, which is variable along the heated length of the fuel channel if a non-uniform Axial Heat Flux Profile (AHFP) is used. In the present chapter, four AHFPs have been applied in order to calculate the fuel centerline temperature in fuel channels at the maximum channel thermal power. These AHFPs are cosine, upstream — skewed cosine, downstream-skewed cosine, and uniform. The aforementioned AHFPs were calculated based on power profiles listed in Leung (2008) while the downstream-skewed AHFP was determined as the mirror image of the upstream-skewed AHFP. A local heat flux can be calculated by multiplying the average heat flux by the corresponding power ratio from Fig. 13.

It should be noted that there are many power profiles in a reactor core. In other words, the axial heat flux profile in each fuel channel differs from those of the other fuel channels. This variation in power profiles is due to the radial and axial power distribution, fuel burn-up, presence of reactivity control mechanisms, and refuelling scheme. Thus, a detailed design requires the maximum thermal power in the core, which can be determined based on neutronic analysis of the core which is beyond the scope of this chapter. However, the four examined AHFPs envelope a wide range of power profiles.

Detectors for the routine monitoring of reactor power in a PWR are located outside the reactor pressure vessel and are characterized by the following typical environmental conditions: neutron flux up to 1011 n cm-2 s-1, gamma irradiation rates up to 106 R h-1, and temperatures of approximately 100 °C. Out-of-core sensors are the usual basis of reactor control and safety channels in a PWR. In choosing specific detector types, consideration must be given to the expected neutron signal level compared with noise sources, the speed of response of the detector, and the ability to discriminate against gamma-induced signals.

Each of these criteria assumes different importance over various ranges of reactor power, and as a result multiple detector systems are usually provided, each designed to cover a specific subset of the power range (Knoll, 2000). Figure 3 illustrates a typical scheme for a PWR in which three sets of sensors with overlapping operating ranges are used to cover the entire power range of the reactor.

The lowest range, usually called the source start-up range, is encountered first when bringing up reactor power from shut-down conditions. This range is characterized by conditions in which the gamma flux from the fission product inventory in the core may be large compared with the small neutron flux at these low power levels. Under these conditions, good discrimination against gamma rays is at a premium. Also, the expected neutron interaction rates will be relatively low in this range. Pulse mode operation of either fission chambers or BF3 proportional counters is therefore possible, and the required gamma-ray discrimination can be accomplished by accepting only the much larger amplitude neutron pulses. As the power level is increased, an intermediate range is encountered in which pulse mode operation is no longer possible because of the excessive neutron interaction rate. In this region the gamma-ray-induced events are still significant compared with the neutron flux, and therefore simple current mode operation is not suitable. The MSV mode of operation can reduce the importance of the gamma-ray signal in this range, but a more common method used in PWRs is to employ direct gamma-ray compensation using a compensated ionization chamber (CIC). A third range of operation corresponds to the region near the full operating power of the reactor. The neutron flux here is usually so large that gamma-ray-induced currents in ion chambers are no longer significant, and simple uncompensated ion chambers are commonly used as the principal neutron sensor. Because these instruments are often part of the reactor safety system, there is a premium on simplicity that also favors uncompensated ion chamber construction.

The diffusion equation is one of the most widely used reactor physical models. It describes the neutron balance in a volume V, the neutron energy may be continuous or discretized (multi group model). The multi group version is:

where the processes leading to energy change are collected in Tkk’:

where subscripts k, k’ label the energy groups, Vk is the speed of neutrons in energy group

k, Yk(r) is the space dependent neutron flux in group k, and keff = 1. In general, the cross-sections Dk, Etk, Ek’^k, Efk’ are the space dependent diffusion constant, the total cross-section, the scattering cross-section, and the fission cross-section. Xk is called the fission spectrum. Equation (4.1) is a set of partial differencial equations, to which the initial condition Yk(r,0),r Є V and a suitable boundary condition, e. g. Yk(r, t),r Є 9V are given for every energy group k and every time t. The boundary conditions used in diffusion problems are of the type

(Vn)Yk(r) + bk(r)Yk(r) = hk(r) k = 1,…, G. (4.3)

for r Є dV. Here bk(r) depends on the boundary condition and may contain material properties, for example albedo.

The diffusion equation is a relationship between the cross-sections in V and the neutron flux Yk(r, t). The equation is linear in Yk(r, t). The main variants of equation (4.1) that are of interest in reactor physics are: l.

where the eigenvalue keff introduced as a parameter in Tkk’ thus allowing for a non-trivial solution Yk(r). That usage is typical in core design calculations.

2. Time dependent solution allowing time dependence in some cross-sections. A typical application is transient analysis.

3. Equation (4.1) is homogeneous but it is possible to add an external source and to seek the response of V to the source.

The structure of the diffusion equation is simple. Mathematical operations, like summation

and differentiation, and multiplication by material parameters (cross-sections) are applied

to the neutron flux. In such equations the symmetries are mostly determined by the space dependence of the material properties. In the next subsection we investigate the possible symmetries of equation (4.4) and the exploitation of those symmetries.

When the solutions (r), k = 1,…, G are known, not only the reaction rates, and net-

and partial currents can be determined, but also matrices can be created to transform these quantities into each other. From diffusion theory it is known that the solution is determined by specifying the entering current along the boundary dV. Thus the boundary flux is also determined. But the given boundary flux also determines the solution everywhere in V. The solution is given formally by a Green’s function as follows:

t G

Yk(r) = E Gk0,k(r0 ^ r)fk0 (r0)dr0. (4.5)

JdVk0=1

Here Gk0,k (r0 ^ r) is the Green’s function, it gives the neutron flux created at point r in energy group k by one neutron entering V at r0 in energy group k0; and fk0 (r0) is the given flux in energy group k0 at boundary point r0. Similarly the net current is obtained as

r G

Jnk(r) = — DkV E Gk0,k(r0 ^ r)fk0 (r0)dr0 (4.6)

JdVk0 = 1

where the V operator acts on variable r.

After applying the oscillation acceleration, the void fraction distribution fluctuated with the same period as the oscillation acceleration. Figure 12 shows the isosurfaces of the void fraction at t = 0.8 s. In the whole area where boiling occurs, the void fraction in the center of the subchannel was relatively low, and the void fraction concentrated in the positive directions along the X and Y axes was high.

Fig. 12. Isosurface distribution of void fraction

red circles in Fig. 13(a). At t = 0.85 s, a high void fraction moved in a direction opposite to the oscillation acceleration as shown in Fig. 13(b). At t = 0.9 s, a high void fraction could be seen near the fuel-rod surface in the regions marked by red circles in Fig. 13(c). This indicates that the magnitude of void fraction fluctuation at Z = 2.3 m is particularly large near the fuel-rod surface. This tendency of void fraction fluctuation is the same between t = 0.9 s and 1.0 s, when the oscillation acceleration acted in a direction opposite to that of the black arrow.

Figure 14 shows the time variation in the void fraction at Z = 3.4 m in the downstream region of Fig. 12. The black arrow shows the direction in which the oscillation acceleration acts. At t = 0.78 s, the vapor phase moved in a direction opposite to the black arrow and concentrated in the regions marked by red circles, shown in Fig. 14(a). At t = 0.8 s, the void fraction in the region marked by the red circles in Fig. 14(b) increased. In addition, a high

void fraction could also be seen away from the fuel rod surface, as shown by the blue circles in Fig. 14(b). The high void fraction in the regions marked by the blue circles in Fig. 14(b) split, and the high void fraction in the blue circles in Fig. 14(c) was formed at t = 0.82 s. The high void fraction regions represented by red and blue circles in Fig. 14(c) moved in a direction opposite to the black arrow as shown Fig. 14(d). While the void fraction regions indicated by the red and blue circles in Fig. 14(d) decreased as shown in Fig. 14(e), high void fraction was concentrated in the regions marked by red circles; high void fraction could also be seen in the regions away from the fuel rod surface, such as the regions indicated by the blue circles at t = 0.9 s, as shown in Fig. 14(f). Near the fuel rod surface, void fluctuation with a different period to that of the oscillation acceleration was seen while the magnitude of the void fraction was relatively small.

Figure 15 shows the time variation in vapor velocity at Z = 3.4 m. The black arrow shows the direction in which the oscillation acceleration acts at each time. At t = 0.78 s, the vapor velocity acted in the direction indicated by red arrows in Fig. 15(a); this direction is opposite, but not parallel to, the black arrow. Between t = 0.8 s and t = 0.82 s, in spite of the changing direction of the oscillation acceleration, the vapor velocity decreased but still acted in the direction of the red arrows, shown in Fig. 15(b) and Fig. 15(c), because of the effect of inertia. At t = 0.8 s, the high void fraction indicated by blue circles in Fig. 14(b) was moved by the vapor velocity. This caused the high void fraction shown in Fig. 15(b) to split, and the high void fraction represented by blue circles in Fig. 14(c) and Fig. 14(d) was formed.

Figure 17 shows the time variation in the Eotvos number at Z = 3.4 m and also shows a range of Eotvos number from 4 to 10 for which the effect of bubble deformation upon the lift force is dependent upon Eotvos number, as shown in Eq. (7). The black arrow shows the direction in which the oscillation acceleration acts. The red and blue circles in Fig. 17 correspond to regions where the magnitude of the lift force was large; the lift force acted in a direction facing away from the fuel rod surface, as shown in Fig. 16. In these regions, the effect of bubble deformation on the lift force was dominant because the Eotvos number exhibited high values. Near the fuel rod surface, the Eotvos numbers less than 4 and greater than 10 were mixed, indicating that the magnitude and direction of the lift force were not uniform near the fuel rod surface.

Figure 18 shows the variation in bubble diameter with time at Z = 3.4 m. The black arrow shows the direction in which the oscillation acceleration acts. Bubble diameters greater than 7 mm are distributed in the region where the Eotvos number is greater than 10, as shown in Fig. 17. The bubble diameter distribution shown in Fig. 18 is strongly inhomogeneous and physically invalid because large bubble diameters are mainly observed in small regions in the subchannel, while small bubble diameters of less than 3 mm are observed in the center of the subchannel. This strongly inhomogeneous bubble diameter distribution resulted in locally high Eotvos numbers and fluctuation in the direction of the lift force vectors.

The region where large bubble diameters are seen corresponds to the region of high void fraction, as shown in Fig. 14. According to Eq. (8), the bubble diameter is significantly
dependent on the void fraction, and a local high void fraction results in a local large bubble diameter. Thus, a strongly inhomogeneous bubble diameter distribution results from void fraction fluctuation.

It is necessary to adequately evaluate the influence of the void fraction upon bubble diameter in order to avoid a strongly inhomogeneous bubble diameter distribution under oscillation conditions.

According to our results, void fraction fluctuation in the downstream region is significantly dependent on the lift force caused by a strongly inhomogeneous bubble diameter distribution.

2. Conclusion

A new external force term, which can simulate the oscillation acceleration, was added to the momentum conservation equations in order to apply the three-dimensional two-fluid model analysis code ACE-3D under earthquake conditions.

A boiling two-phase flow excited by applying vertical and horizontal oscillation acceleration was simulated in order to confirm that the simulation can be performed under oscillation conditions. It was confirmed that the void fraction fluctuation with the same period as that of the oscillation acceleration could be calculated in the case of both horizontal and vertical oscillation acceleration.

The influence of the oscillation period of the oscillation acceleration on the boiling two — phase flow behavior in a fuel assembly was investigated in order to evaluate the highest frequency necessary for the improved method to be consistent with the time-series data of oscillation acceleration and the shortest period of oscillation acceleration for which the boiling two-phase flow shows quasi-steady time variation. It was confirmed that a boiling two-phase flow analysis consistent with the time-series data of oscillation acceleration and with a time interval greater than 0.01 s, can be performed. It was also shown that an effective analysis can be performed by extracting an earthquake motion of about 1 s at any time during the earthquake.

The three-dimensional behavior of boiling two-phase flow in a fuel assembly under oscillation conditions was evaluated using a simulated fuel assembly excited by oscillation acceleration. On the basis of this evaluation, it was confirmed that void fraction fluctuation
in the downstream region is significantly dependent on the lift force caused by a strongly inhomogeneous bubble diameter distribution and that it is necessary to adequately evaluate the influence of void fraction on bubble diameter in order to avoid strongly inhomogeneous bubble diameter distribution under oscillation conditions.

3. Acknowledgment

The present study includes the result of "Research of simulation technology for estimation of quake-proof strength of nuclear power plant" conducted by the University of Tokyo as Core Research for Evolutional Science and Technology (CREST). This research was conducted using a supercomputer of the Japan Atomic Energy Agency.

Some developed codes and simulators for improving the power monitoring will be reviewed in this section. For example, MCNP (monte-carlo n-particle transport code) is developed for neutron detector design, or modeling a fission chamber to optimize its performance

7.1 Computational tools to conduct experimental optimization

Research reactors need a handy computational tool to predict spatial flux changes and following power distribution due to experimental requirements. Therefore it is important to get accurate and precise information ahead of any modifications. To meet this demand, flux measurements were conducted in case that a typical flux trap inside the core to be allocated. In TRR, one of standard fuel boxes, in position D6 in core configuration of the year 1999, was taken out of the core and a water trap was formed in its place. With the aid of miniature neutron detector (MND) using standard procedure, thermal neutron flux is measured inside the water trap. To calculate the flux and power theoretically, two different computational approaches such as diffusion and Monte Carlo methods were chosen. Combination of cell calculation transport code, WIMS-D5, and three-dimensional core calculation diffusion code, such as CITATION, were used to calculate neutron flux inside the whole core either in two or five energy groups. However, MCNP-4, as a Monte Carlo code, was used to calculate neutron flux again inside the whole core as well as inside the trap (Khalafi and Gharib, 1999). Figure 14 shows axial thermal flux distribution along the D6 position by measurement and computation.

It is obvious from the figure; the both calculation codes are satisfactory and a good agreement exists between detector measurements and code computations. However, diffusion method is a rational choice especially for survey calculation where the Monte Carlo approach is more time demanding. For some consideration, in order to measure spectrum, a fixed point on the midplane along D6 axis was chosen. A variety of foils of different material was selected as measuring windows to determine differential fluxes at specified energy bins. Metal foils such as Ti, Se, Mg, Ni, Al, Co, Au, In, and Fe were selected as energy windows. These foils are sensitive to a part of neutron energy spectrum starting form high energies and ending to thermal energies. Induced activity of each foil is measured based on gamma spectroscopy using high purity Germanium (HPGe) detector. By providing raw counts to SAND-II computer code, neutron energy spectrum was calculated. The measured and calculated spectrum using neutron detector, MCNP and WIMS codes is shown in Figure 15.

Spectrum calculations were also checked against measurements. Monte Carlo shows a better prediction while WIMS provides a fair result. It is notable that combination of WIMS/CITATION would be sufficient for neutron flux calculations while Monte Carlo technique should be reserved for the final stages of simulation. A good choice of computational tools would save time a lot in this respect and one is encouraged to perform a comprehensive simulation ahead of design and construction of irradiation facility.

A comparison between the three 99Mo production methods (the fission method, the solid irradiation method and the solution irradiation method) is shown in Table 1, assuming the 99Mo production in JMTR, which is a tank-type reactor.

1.1 Fission method ((n, f) method)

In the conventional fission method ((n, f) method), high-enriched uranium targets are irradiated with neutrons in a testing reactor, and 99Mo is produced by the 235U (n, f) 99Mo reaction. Most of the world supply of 99Mo is produced by the (n, f) method since 99Mo with a high-level specific activity of 370 TBq/ g-Mo is obtained. However, the method has problems about the nuclear nonproliferation and the generation of a significant amount of radioactive waste including Fission Products (FPs) and Pu. Caused by the radioactive waste, the separation process of 99Mo is too complex, and 99Mo production with the (n, f) method needs expensive facilities and extreme care to avoid contamination with FPs. The 99Mo production cost by this method achieves 57 US$/37 GBq (Boyd, 1997), and it is too expensive.