- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

Science and Technology of Nuclear Installations

Volume 2011 (2011), Article ID 693245, 12 pages

http://dx.doi.org/10.1155/2011/693245

## Experimental Validation of RELAP5 and TRACE5 for Licensing Studies of the Boron Injection System of Atucha II

^{1}Nuclear Regulatory Authority of Argentina, Buenos Aires, C1429BNP, Argentina^{2}School of Nuclear Engineering, Purdue University, West Lafayette, IN 47907, USA

Received 4 August 2010; Accepted 16 November 2010

Academic Editor: Alejandro Clausse

Copyright © 2011 Alejandro I. Lazarte et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This paper presents an experimental validation of RELAP5 and TRACE5 for licensing studies of the Atucha II-PHWR nuclear power plant. A scaled experimental facility, representing the boron injection system of Atucha II, was built. The system has a fundamental importance for loss of coolant accidents (LOCA) and anticipated transients without scram (ATWS). The experiment consists of the discharge of a tank that represents the boron tank filled with air or a mixture of air-water onto a discharge tank that represents the moderator tank. Both tanks are connected by a pipe which includes a valve and an orifice plate to model the pressure losses due to the fittings in the real system. The pressure and water level measured in the tanks are compared with the RELAP5 and TRACE5 predictions. The codes predict the pressure in the tanks accurately. However, both codes overpredict the heat transfer in the boron tank air-water interface which produces a greater expansion of the air which leads to a small discrepancy in the boron tank level prediction.

#### 1. Introduction

The objective of this paper is to benchmark and compare the accuracy of RELAP5 MOD 3.3 and TRAC/RELAP (TRACE5) Patch 1 to a high-pressure injection experiment. RELAP5 and TRACE5 are the two computer codes based on the one-dimensional two-fluid model approved by the USNRC to perform nuclear reactor licensing calculations of LOCA events and other transients. During the last three decades, RELAP5 has been compared, assessed, and validated against several separate effects and integral test facilities [1] as well as TRACE5 [2] and its predecessor, TRAC. Because of this, the codes are suitable, among other things, to be used for modeling this safety system with good accuracy.

The 1/2 scale experiment was performed at Purdue University to simulate several aspects of the boron injection (JDJ) system of the Atucha II NPP. The prototypic JDJ system consists of a high-pressure air tank connected to two boron tanks at atmospheric pressure by a fast opening valve [3, 4]. The boron tanks inject boron into the reactor moderator tank through a long pipe (approximately 30 m). A rupture disk in the pipe is designed to burst when the forward pressure difference exceeds 31.5 bar (3.15 MPa). The components of the JDJ system are shown in Figure 1.

The experiment was designed to validate the performance of the JDJ system for the license of Atucha II [5]. The licensing calculations were performed in two steps. The 1/2 scale experiment was used to validate the computational models (i.e., RELAP5 or TRACE5). These results are presented in this paper. The second step was the scaleup of the validated computational models to prototypic dimensions and conditions. The second step differs from traditional scaling, where the experimental results are used directly, because the computational model allows the correction of the scaling distortions and the simplifications of the experiment.

This paper is divided into four sections. The first is the description of the scaling of the experiment of the JDJ system. To begin the design of the experiment it was necessary to identify the relevant phenomena including any effect that may change its performance. Then, scaling laws based on the identified phenomena were used to design a 1/2 scale experiment with a pressure difference of 100 bar (10 MPa) between the tanks. Some components which were not simulated are the air tank and the fast opening depressurization valve because the exact valve was not available. This valve needs to be tested separately under dynamic conditions. The piston in the boron tank, separating air and water, was not simulated either. This is followed by a description of the experiment and the experimental results. The experiment was conducted with pressure and level measurements to determine the boron injection speed. The data was collected with a PC-based data acquisition system. Next the numerical models of the experiment are described for both RELAP5 and TRACE5. The model results were compared with the data and to each other to determine whether RELAP5 and/or TRACE5 are satisfactory for this application. The last part is a discussion of the results followed by the conclusions.

#### 2. Scaling

The experiment scaling is based on the approach used at Purdue University PUMA facility for the safety systems of the SBWR (i.e., Simplified Boiling Water Reactor) [6]. However, the task is made much simpler because of the physical boundary in the JDJ system separating the phases. Furthermore, the experimental procedure of the simplified system is such that only single-phase relevant phenomena need to be considered. The data used for the scaling was obtained from the JDJ system report [3].

For scaling, since the time of the boron injection delay is very short, the flow is assumed to be adiabatic. Furthermore, the water flow is assumed to be single-phase incompressible. These two assumptions make the conservation of mass and conservation of energy scaling trivial. Then the model is derived from a simple integral momentum balance of the water in the boron tank and pipe that is derived in Appendix A: where is the liquid density. The first term on the RHS is the pressure difference between the tanks, and the last term is the frictional and other irreversible losses in the pipe, entrance contraction, and exit expansion.

An integral model is described in Appendix B. The model consists of three first-order differential equations: conservation of air mass in the boron tank, conservation of energy in the boron tank, and conservation of momentum in the pipe joining the boron and discharge tanks. The pressure in the discharge tank is assumed to be constant. The *Mathematica* [7] model written to solve the equations is shown in Appendix B. The pressure difference between the two tanks remains approximately constant during the transient corresponding to the boron injection time, and the *Mathematica* model shows that this pressure difference is approximately
Then,
This equation may be nondimensionalized as follows,
Finally defining , one obtains
If and are the same, the experiment is similar to the prototype. The definition gives flexibility to the design of the experiment. In fact it is possible to scale down the size and pressure of the experiment significantly. *However, in order to consider unforeseen effects, it was decided to scale down the facility only by a linear geometric factor of 1/2 and to maintain the pressure difference between the tanks equal to the JDJ system value (i.e., 100 bar). *Since is the same, this implies that the liquid velocities in the experiment and the JDJ system are approximately the same, and since is 1/2, similarity of means that .

The scaling based on an integral equation does not imply that every component will be simulated exactly. Detailed phenomena like turbulence will be slightly different because the Reynolds number of the experiment is 1/2 of the prototypical value. Also many fittings like elbows and tees are different. Finally, the rupture disk assembly is different. An orifice plate is used to compensate these differences and may be adjusted to satisfy the integral scaling criterion.

#### 3. Experiment

The JDJ system schematic is shown in Figure 1. The air tank is pressurized to about 200 bar (20 MPa), while the boron tank is at 1 bar (0.1 MPa). The volumes are 1 m^{3} for the air tank and 0.2 m^{3} for each boron tank. The latter is divided into two parts: one containing borated water (0.162 m^{3}) and the other containing air (0.038 m^{3}) separated by a piston (stainless steel). Once the signal is given for opening the quick opening air valve, the air tank discharges onto the boron tanks, increasing the pressure. When the forward pressure difference between the latter and the moderator tank is greater than 31.5 bar (3.15 MPa), the bursting disc ruptures and borated water injects into the moderator.

The experimental schematic is shown in Figure 2. The tanks used for the experiment are shown in Figure 3. The experiment is simplified to one boron tank, where the volume of the tank is 1/8 (i.e., a linear scale of 1/2) of the volume of the two boron tanks of the JDJ system. The second major simplification is that the air tank subsystem of the JDJ system is not simulated in the experiment. Instead the air volume at the top of the boron tank is slowly pressurized with an air supply system. The quasistatic pressurization eliminates mixing in the boron tank, and therefore no piston is necessary. These simplifications were performed because the exact fast opening valve was not available and to reduce the cost of the experiment.

The 103.5 bar (MPa) is 53.5 L (0.0535 m^{3}). Boron tank is instrumented with a 140 bar (14 MPa) pressure and a 0.14 bar (14,000 Pa) differential pressure transducers. The flow rate of water between the boron tank and discharge tank is obtained from the differential pressure sensor that measures the water level in the boron tank. Additionally, a 276 bar (27.6 MPa) analog pressure gauge is used for instantaneous pressure readings. A 103.5 bar (10.35 MPa) ASME pressure relief valve has been installed to the boron tank to meet safety regulations. The outlet below the boron tank contains a tee for three purposes: to fill the boron tank with water, to drain the experiment of water, and to direct the water flow towards the discharge tank during the experimental runs. The air supply system consists of a 150 bar (15 MPa) standard air supply cylinder containing 8.8 m^{3} of air at NTP. The air flow rate to the boron tank is regulated by the 103.5 bar (10.35 MPa) pressure regulator.

The diameter of the pipe between the boron tank and the discharge tank (i.e., moderator tank) is 1–1/2 (i.e., 38.6 mm) which is approximately 1/2 of the prototypic value (i.e., 68.9 mm). The length of the pipe to the tip of the lance is 20.7 m which is roughly 1/2 of the length of the prototype (i.e., 36 m). All tanks and piping are made of steel. The maximum pressure difference between the boron tank and discharge tank is 102 bar (10.2 MPa), which matches the prototypic value. A rupture disk is located in the pipe to simulate the original rupture disk design (i.e., size). However, the disk was not used for the cases presented in this paper. Instead, the ball valve upstream of the rupture disk was quickly opened to start the transient. The irreversible loss in the ball valve is much smaller than the loss in the rupture disk. The results presented are characterization tests and do not account for the loss of the disk. These tests are not similar to the prototype; however, other cases involving the disk were. For these characterization tests, we have accounted for this difference in the computational models. The 13.5 bar (MPa) 907 L (0.907 m^{3}) discharge tank is instrumented with a 20 bar (2 MPa) pressure transducer. Additionally, a 27.5 bar (2.75 MPa) analog pressure gauge has been added for visual pressure readings. A 13.5 bar (1.35 MPa) ASME pressure relief valve has been installed to meet safety regulations. The discharge tank has two view ports to perform flow visualization. A video camera was placed on a viewport of the discharge tank to check the incoming water flow for cavitation. The National Instruments USB-6008 data acquisition system, connected to a PC, was used to collect the data. Data was collected from the three pressure transducers at 1 kHz.

Shakedown tests were initially performed to test the experimental equipment, procedure, and instrumentation performance. One significant phenomenon that was monitored during shakedown was cavitation in the pipe. Since no cavitation occurs in the prototype, it is necessary to avoid cavitation in the experiment. The discharge tank was designed to suppress cavitation by pressurizing it with air. All tests were performed at an initial discharge tank pressure of 8 bars (0.8 MPa) or higher to prevent cavitation.

Another important phenomenon concerning the performance of the experiment is the presence of the rarefaction (depressurization) wave at the beginning of the experiment which is captured by the DP transducer. The wave appears as a dip shown in Figure 7 for the valve tests. The rarefaction wave is only seen in the DP measurement as this was the only sensor with taps on *both* sides of the semireflective interface.

The data tests presented in this paper summarized in Table 1 and are of two types: (i)one air-only tests to check the air models of RELAP5 and TRACE5. These tests were performed with an orifice plate to specify the exact point at which the flow chokes in the calculations;(ii)one air-water tests without a rupture disk to test the complete models and to characterize the pipe. The idea was to make tests without the uncertainty of the rupture disk area. The tests were initiated by opening a manual ball valve in the pipe. Because of safety concerns, these tests were performed below 25 bar.

#### 4. RELAP5 Model

The RELAP5 nodalization is shown in Figure 4. The dimensions of the RELAP5 model were kept as close to the actual experiment as possible with distance and elevation measurements of each section of pipe. The aspect ratio of the computational nodes was chosen to be order 1 as recommended in the manual [1]. All pipe nodes have an inside diameter of 38.1 mm which is taken from manufacturer’s data. The nodes are grouped into “pipes” as shown in Figure 4. For example, “pipe” 107 consists of 90 nodes, 3.81 cm each. The node lengths have a range between 3.81 and 7.37 cm (i.e., 1 to 2 pipe diameters). By comparison with experimental results, a wall roughness of 50 *μ*m was selected for all nodes. Junctions representing 90° and 45° elbows have a Reynolds-independent loss factor of 0.32 and 0.16, respectively. The overall was measured. The roughness is not the same as the prototypic value, so this is a scaling distortion. The Ransom-Trapp critical flow model was specified, although no two-phase choking should occur before the boron tank empties. Choking was deactivated at all junctions except the boron tank outlet, the orifice plate, and the discharge tank inlet. At these junctions the abrupt area change model was turned on. The mixture level tacking model was specified in the boron tank and the first vertical pipe below the boron tank. This model helped in the level calculations when the level dropped out of the boron tank. Dry air was specified in the boron tank and saturated (humid) air was specified in the discharge tank. The simulation begins at 0.5 seconds when a trip activates the valve.

The 14 bar (1.4 MPa) air test comparisons with RELAP5 are shown in Figures 5 and 6. This test was performed with an orifice plate, where the flow chokes during most of the transient. Two models were tested, without and with heat structures to represent the walls of the tanks and the pipe. Due to the short characteristic time of the transient, negligible heat transfer is assumed to occur with the environment. Therefore, an insulated, or adiabatic, boundary condition was used on the outer wall of the pipe and tanks. It is seen that including the heat structures has a significant effect on the pressures and that the results become more accurate. The RELAP5 predictions match the data closely indicating that the dry air assumption is satisfactory.

The 9 bar (0.9 MPa) air-water test comparisons are shown in Figures 11 to 13. The test was performed with either orifice or rupture disk. The figures show results for three different pipe roughness. Good agreement is obtained with a roughness of 50 *μ*m. There is some discrepancy in the level, specifically the asymptotic value, and also a slight difference in the pressure of the boron tank.

The reason why the RELAP5 pipe model cannot predict the level and the boron tank pressure more accurately is that the interfacial heat transfer coefficient between the air and the liquid surface in the boron tank is too high. In fact it is so high that the RELAP5 pipe model practically predicts thermodynamic equilibrium in the control volumes, where the interface is located. This overprediction of the energy transfer between the water at room temperature and the cold air causes the air to expand more than it should, and therefore the asymptotic water level in the boron tank predicted by the model is lower than the data. This effect is small for the low-pressure difference cases presented here. However, for higher pressures, such as the prototypic case, this effect becomes significant. Another consequence of this discrepancy is that the RELAP5 overpredicts the water flow rate in the pipe. However, this error only has a small effect in the initial part of the transient corresponding to the boron injection time.

The velocity in the orifice is show in Figure 9, where the chocked flow during the first instant of the transient can be seen. Also, in the same figure, the velocity for the TRACE5 model is illustrated, obtained using the TPR exported file and the TRACE5 nodalization. Clearly, the RELAP5 and TRACE5 results are similar, showing in the TRACE5 the choked flow during roughly 1.5 s after valve opening.

#### 5. TRACE5 Model

The TRACE5 model was built using the RELAP5 nodalization by means of the SNAP tool. For the sake of simplicity, a nodalization is not shown since it is similar as it was illustrated in Figure 4 without junctions between components. The dimensions of the model are exactly equal to the RELAP5 model as well as the aspect ratio in each control volume, considering that junctions at close end pipes have a null area. The initial conditions were set equal to RELAP5 model. The optional multipliers for choked flow (ICFLOW = 2) were activated using the default values in the same junctions as it was done in RELAP5. At these contraction/expansion junctions, the abrupt area change was also selected. At all other junctions, the friction option FRIC = 1 was selected, except at closed end of the pipe, where FRIC = 0 was activated. All the irreversible -factor pressure losses were input using forward and reverse KFAC as it is recommended in user’s guideline for new nodalizations [2]. Dry and saturated (humid) air was also specified as it was done in RELAP5. In order to implement the same boundary condition, the pressures for fluid and noncondensable were extracted from the exported TPR file during conversion.

The material selected for heat structure modeling was carbon steel. The equation of state used was the IAPWS for all the cases since the curve fitting table required a lower maximum time step to avoid code failures without having a beneficial effect. The numerical scheme used in TRACE is the default option (SETS).

The 14 bar (1.4 MPa) air test comparisons with TRACE5 are shown in Figures 7 to 9. A description of this test had been already given in the previous subsection. In this case, the model with heat structure representing the tank wall is more accurate and even has better agreement with data than the RELAP5 model. However, a slightly negative slope is seen at the end of the transient. Also a comparison with results obtained with TRACE5.0 is shown. The main discrepancy is that the model underestimates the asymptotic values for the pressure. Comparing the pressure in moderator tanks for the RELAP5 and TRACE5 models, it could be seen that RELAP5 overpredicts while TRACE underpredicts the asymptotic pressures.

The 9 bar (0.9 MPa) air-water test comparisons are shown in Figures 13 to 16. Figures 13 to 15 show results for three different pipe wall roughnesses. Similar to the RELAP5 model, best agreement is obtained with a roughness of 50 *μ*m. However, the RELAP5 model prediction is more accurate than the TRACE5 model for the pressures, but the level prediction is slightly better. In these cases, the heat structures in all pipes were considered although the results are not shown because the heat structures model for the air-water experiments has a small effect compared with air-only test, where expansion is more abrupt and faster.

Figure 16 shows the fluid discharge velocity into the moderator tank. This is one of the more relevant parameters of this safety system. RELAP5 and TRACE5 are very similar; the first peak value is almost equal, and a small time shift may be seen.

#### 6. Discussion

Data was collected to test the RELAP5 and TRACE5 models for the air expansion in the boron tank and the quick injection of boron into the moderator tank. The comparisons of the RELAP5 and TRACE5 models with the data show that both codes provide a satisfactory prediction for pressure and level.

The air-only test indicates that the effect of the heat structures (i.e., tank and pipe walls) is considerable and they need to be included to model the air expansion more accurately.

The air-water models show good general agreement with the data. However, when the experiment was run at much higher pressures, not shown in this paper, the level prediction diverged more significantly from the data. The reason for this is believed to be an overprediction of the heat transfer at the boron tank air-water interface. When the interface is in a cell, the interface-to-gas, , and direct, , heat transfer coefficients are of the order 10^{9} and 10^{6} W/m^{3}K, respectively, in RELAP5. This results in gas temperatures being very near to the liquid temperature, when the level is in the cell, as shown in Figure 18. Figure 19 shows that TRACE5 predicts thermodynamic equilibrium in the cell, where the interface is located. Although the cell that contains the interface is in near thermodynamic equilibrium (i.e., isothermal condition), once the interface leaves the cell, the air temperature is allowed to decrease during the expansion. The expansion of the air then becomes a polytropic process controlled by the heat transfer between the air and the wall of the tank which the codes predict correctly. This problem may be reduced, but not resolved, by refining the mesh.

In order to study the air gas thermodynamic behavior in the boron tank discharge, a PV diagram is shown in Figure 17. The experimental curve was built with the measured pressure in the upper part of the boron tank and with air volume as function of the level from the DP transducer. Also, the figure shows an isothermal curve and an adiabatic curve (with ). Both curves were built with TRACE5 assuming an isothermal or adiabatic quasistatic evolution using the level (i.e., volume) as the independent variable and calculating the pressure. The first point in the diagram corresponds to the instant previous to the valve opening. The curve from the experimental data shows the depressurization wave during the first 300 ms as measured by the DP detector that can also be seen in Figure 12. After the pressure wave, the experimental results and the RELAP5 and TRACE5 show a very similar behavior. The results show that a polyprotic evolution with a ratio of specific heats () fits reasonably well.

#### 7. Conclusions

The codes predict the pressure in the tanks accurately. However, both codes overpredict the heat transfer in the air-water interface in the boron tank which produces a greater expansion of the air which leads to a small discrepancy in the boron tank level prediction. Furthermore, the good agreement between the integral model and the data (see Appendix B) shows consistency at all levels (i.e., measurements, integral model, and differential models) and provides even more validation. Therefore, both codes may be used for licensing studies and to test any minor modification of the boron injection system of Atucha II.

#### Appendices

#### A. Scaling Equation

The integral momentum equation for a component of uniform cross-sectional area is The LHS is the fluid inertia, the first term on the RHS is the convection, the second term is the pressure difference across the component multiplied by the cross-sectional area, and the last term is the irreversible losses in the pipe. The irreversible losses are comprised of the friction factor and the form losses. The effect of gravity is neglected. The elevation change in the experiment is about 1 m, and the prototype is slightly lower than 8.5 m. Two control volumes of uniform cross-section are considered: the pipe and the boron tank (Figure 20).

For the boron tank the control volume upper side moves with the interface. Then applying Leibnitz integral rule to the LHS (A.1) becomes In the first term in the RHS since the top moves with the water surface. Furthermore, the last term in the RHS may be neglected. Then (A.2) may be written as For incompressible flow so the first term in the RHS is zero. Dividing by , or The contraction may be treated as a singularity: For incompressible flow in the pipe the so the convection term in (A.1) is zero and the LHS is simplified since is constant. Then dividing by the pipe cross-sectional area, or The value in (A.8) stands for the accessories in the system, basically, the valve. The expansion between the pipe and the discharge tank may be treated as a singularity: Finally adding (A.5), (A.6), (A.8), and (A.9) yields where . This equation may be simplified considering that and inserting and in ,

#### B. Integral Model

The model consists of three first-order differential equations:(1)conservation of air mass in boron tank,(2)conservation of energy in boron tanks, and(3)conservation of momentum in the pipe.

The air quantities in the boron tank are labeled 1, and the liquid properties are labeled .

The conservation energy equation in the boron tank is written in terms of the pressure because the flow is calculated in terms of the pressure difference. The relationship between the time derivatives of the pressure and the density is obtained from the polytropic relationship: where is the constant ratio of the specific heats. Then Performing the derivative by parts yields which is simplified as follows: The conservation equation for the air in tank 1 is where is the volume of the air. Then combining (B.4) and (B.5) results in This equation considers the change in the air volume of the boron tanks as the water level decreases. The equation for the air volume in the boron tanks is simply Combining (B.6) and (B.7) yields Finally, combining (B.4) and (B.8) one obtains Equations (B.7) and (B.9) represent the desired conservation of mass and energy equations.

The conservation of momentum equation is derived from an integral momentum balance (see Appendix A) applied to the control volume corresponding to the pipe from the boron tank to the moderator tank: The LHS is the fluid inertia, the first term on the RHS is the convection, and the second term is the pressure difference between the tanks multiplied by the area of the pipe, the last term is the irreversible losses in the pipe. The irreversible losses are comprised of the friction factor and the form losses. For an incompressible flow the above equation may be simplified as Then, This then is the equation for the boron injection velocity, and the mass flow in the piping is .

A *Mathematica* model has been developed to solve these equations. Figures 21 and 22 show the pressure and velocity during the initial part of the transient corresponding to the boron injection time, which is the interval of most interest for licensing. However, the data are inaccurate during this time phase and are not shown. In addition, no mass flow or velocity measurement were done. The change in the pressure during this time is relatively small so the scaling assumption given in (2) is justified.

The long-term results for the pressure and the level in the boron tank are shown in Figures 23 and 24. Very good agreement with the experimental data is demonstrated in the figures. These results were calculated with an air polytropic coefficient of 1.25 that was experimentally determined (see Figure 17). This demonstrates that the physical mechanisms identified and modeled in this appendix are sufficient.

#### Acknowledgments

This work was funded by the Nuclear Regulatory Authority of Argentina. The authors would like to thank Deoras Prabhudharwadkar, Kevin Mueller, Matthew Kennedy, and Victor Kovach for their work on the experimental test section and Dr. M. Ishii and Dr. V. Ramson’s expert opinions in the experiment and the RELAP5 model, respectively.

#### References

- RELAP5/MOD3.3 Code Manual volume I: Code Structure, System Models, and Solution Methods, Nuclear Safety Analysis Division, Information Systems Laboratories Inc., December 2001.
- Trace V5.0 Theory Manual and User’s Manuals. Division of System Analysis Office of Nuclear Regulatory Research U. S. Nuclear Regulatory Commission Washington, DC, USA 20555-0001.
- G. Salom, “CNA II: Modelo Hidraulico Del Sistema De Inyeccion De Boro De Emergencia (JDJ),” Informe NA-SA LN-18-2007, 2007.
- Blueprint of the system JDJ10. KWU R34/92/1154. S322-11-93462.
- “CSAU-BEPU statistical methods for sensitivity and uncertainty analysis of atucha ii loss of coolant accidents using relap5,” Purdue Report OZ622A, Contract, February 2009.
- M. Ishii, S. T. Revankar, T. Leonardi et al., “The three-level scaling approach with application to the Purdue University Multi-Dimensional Integral Test Assembly (PUMA),”
*Nuclear Engineering and Design*, vol. 186, no. 1-2, pp. 177–211, 1998. View at Scopus - S Wolfram,
*The Mathematica Book*, Wolfram Media/Cambridge University Press, New York, NY, USA, 4th edition, 1999.