Extensions of the fast Fourier transform based method for quantitative assessment of code calculations
Andrej Prošek, Matjaž Leskovar
Jožef Stefan Institute, Jamova cesta 39, SI-1000 Ljubljana, Slovenia E-pošta: andrej.prosek@ijs.si, matjaz.leskovar@ijs.si
Abstract. In the last years, the number of quantitative comparisons between the experimental data and the calculated data in the area of nuclear technology is increasing. The most widely used method for the accuracy assessment of thermal-hydraulic code calculations was the fast Fourier transform based method (FFTBM). Recent applications of the original FFTBM showed the need for further improvements. The purpose of the present study was therefore to extend the FFTBM. The most important improvement was the introduction of an index for the indication of the time shift between the compared signals. Also was proposed to properly prepare the time domain signals before the fast Fourier transform is performed. Namely, little can be done in this respect in the frequency domain. For oscillations due to measuring reasons (noise), the moving average was proposed, for phenomena of the logarithmic nature the calculation of the logarithm of the time signals was proposed, etc. To validate the extended FFTBM, the results of the Organisation for Economic Co-operation and Development (OECD) and the International Atomic Energy Agency (IAEA) international studies were used. The results show that the extended FFTBM correctly indicates the time shift and that it addresses the accuracy in a consistent way.
Keywords: FFTBM, time-shift indication, accuracy, quantitative assessment, thermal-hydraulics
Razširitev metode na podlagi hitre Fourierjeve transformacije za kvantitativno ocenitev računalniških izračunov
Povzetek. V zadnjih letih na področju jedrske tehnologije čedalje bolj narašča število kvantitativnih primerjav med eksperimentalnimi podatki in rezultati izračunov. Najbolj uporabljena je bila metoda na podlagi hitre Fourierjeve transformacije (FFTBM), ki je namenjena za določanje natančnosti termohidravličnih izračunov. Dosedanje uporabe FFTBM so pokazale potrebo po nadaljnjih izboljšavah. Zato je bil namen študije razširiti metodo FFTBM. Najpomembnejša izboljšava je uvedba kazalca za časovni premik med primerjanima signaloma. Hkrati smo predlagali, da se kvantitativna ocenitev izboljša z ustrezno pripravo časovnih signalov, preden uporabimo hitro Fourierjevo transformacijo, ker se v frekvenčnem prostoru ne da veliko narediti. Za izločitev oscilacij zaradi meritev (šum) smo predlagali drseče povprečje, za pojave, katerih narava je logaritemska, smo predlagali logaritmiranje časovnih signalov ipd. Za potrjevanje FFTBM so bili uporabljeni rezultati mednarodnih študij v okviru Organizacije za gospodarsko sodelovanje in razvoj (OECD) in Mednarodne agencije za atomsko energijo (IAEA). Rezultati kažejo, da FFTBM pravilno zazna časovni premik in da konsistentno izračuna natančnost.
Ključne besede: FFTBM, indikacija časovnega premika, natančnost, kvantitativna ocenitev, termohidravlika
1 Introduction
The fast Fourier transform based method (FFTBM) was proposed in 1990 [1]. Since then it has been used for
accuracy quantification of code calculations [2], [3] and it became a standard tool for nodalization qualification
[4].	FFTBM shows the measurement-prediction discrepancies - accuracy quantification - in the frequency domain. It is incorporated in a modified way also in the automated code assessment program (ACAP)
[5],	developed to provide quantitative comparisons between nuclear reactor systems code results and experimental measurements. ACAP is intended for single variable comparison only, while with FFTBM the total code accuracy can also be assessed. Significant improvements of FFTBM have been done in the last years. FFTBM with new accuracy measures is described in detail in [6] and was applied to the main coolant pump trip transient in the Mochovce nuclear power plant [7]. Later, FFTBM with the capability to calculate time dependent code accuracy was developed and was used for quantifying code calculations of the large-break loss-of-coolant accident in the RD-14M facility [8]. Recently, it was observed that FFTBM favours certain trends when an edge (difference) is present in the signal between the first and the last data point of the investigated time signal. Namely, the discrete Fourier transform used for the code accuracy calculation views the time domain signal as an infinite periodic signal. The problem of the edge effect was resolved in a unique way by signal mirroring [9].
Recent applications of FFTBM have shown the need for a few further extension of the method. For example, in the original FFTBM, the time shift is determined in the qualitative analysis. To eliminate subjectivity, a quantitative tool for time shift indication was developed. The time shift of the results is very important in safety analyses, as the transient progression mostly depends on the primary pressure. Many safety systems actuate on set pressure setpoints. The discrepancies in the pressure can therefore cause time shifting of other variables, e.g. the calculated fuel cladding temperature, while the calculated trend may be good. Other improvements towards making a more objective judgement based on the results of FFTBM are also presented. The first is elimination of the noise (if present) from the signal by the moving average, and the second is consideration of the logarithmic nature of the phenomenon. Once the time shift is indicated, shifting must be made in the time domain. Similarly, selecting the units (e.g. K or °C) may affect the results because of the zero-frequency component (the mean value of the time-domain signal).
In Section 2, the FFTBM method is described first. In Section 3, the results of the extended FFTBM validation are presented for selected examples. Finally, in Section 4, some conclusions are drawn.
2 Extended FFTBM method description
The methodology of the code-accuracy assessment consists of three steps: a) selection of the test case (experimental or plant measured data to compare), b) qualitative analysis, and c) quantitative analysis. The qualitative analysis is a prerequisite to perform the quantitative analysis. The qualitative analysis, including visual observation of plots, is done by evaluating and ranking the discrepancies between the measured and calculated variable trends. The quantitative analysis is meaningless unless all the important phenomena are predicted.
For brevity reasons, only the basic quantitative accuracy measures with no derivations will be described together with proposed extensions. For more details, the reader can refer to refs. [3], [6], [9].
2.1 Basic accuracy measures
The accuracy measures are defined in the frequency domain and the discrete Fourier transform of time record data is computed with the fast Fourier transform (FFT). To calculate discrepancies, the experimental signal (Fexp(t)) and error function are needed. The error function in the time domain AF(t) is defined:
AF (t) = Fa() - FeJ),	(1)
where Fcal(t) is the calculated signal. The code accuracy quantification for an individual calculated variable is
based on amplitudes of the discrete experimental
l^xpft^j and error signal |2i F(fn) obtained by FFT at
frequencies (n, where n=0,1.....2m (N=2m) and m is the
exponent (m=(8, 9, 10, 11)). These spectra of amplitudes are used for the calculation of the average amplitude (AA) that characterizes the code accuracy:
2m
X| ~F(fn\
AA = -n=°-.	(2)
2m
X Fexp((j
n=0
AA is the relative magnitude of the discrepancy coming from the comparison between the calculation and the corresponding experimental variable time history. When the calculated and the experimental data are equal, then the error function is zero (average amplitude is also equal to zero), characterizing perfect agreement.
The overall picture of the accuracy for a given code calculation is obtained by defining the total weighted average amplitude (total accuracy):
Nvar
AAot = X (AA\.(wf )i,	(3)
i=1
with
Nvar
X (Wf ) = 1,	(4)
i=1
where Nvar is the number of the variables analysed, and (AA)i and (wf)i are the average amplitude and the weighting factors for the i-th analyzed variable, respectively. Detailed information on weighting factors can be found in ref. [6]. The acceptability factor for the calculation is 0.4, while for pressure, the criterion for AA is to be below 0.1.
2.2 Index for time shift indication
It should be noted that the AA figure of merit (Eq. (2)) was not obtained by comparing the experimental and calculated magnitude spectra, but by calculating the magnitude spectrum of the difference signal. Nevertheless, due to the Fourier transform properties the magnitude spectrum of the difference signal can also be obtained by adding the experimental and calculated signal magnitude spectra (actually subtraction); they must be converted into a rectangular notation, added, and then reconverted back to a polar form. When the spectra are in a polar form, they cannot be added by simply adding the magnitudes and phases. The error
function amplitude spectrum DF(fn) can be expressed as:
\^F(fn\ = \Fef - FJfn\ =
= ^(Re( Ff - FJfJ)2 + (lm( FeJfn) - FcJfJ))2 = = yjM2 + M22 - 2M1M2(cos jcos j2 + sin jsin j2),
(5)
where	Fep = M1 cos j + iM1 sin j	and
Fml = M2 cos j2 + iM2 sin j2 (rectangular form).
This example shows that to calculate the difference magnitude spectrum we need both the magnitude and the phase of the experimental and calculated spectra. Information about the shape of the time domain signal is contained in the magnitude and in the phase. In other words, comparing the shapes of the time domain signals is done through calculating the difference signal magnitude spectrum. At the time of the development of the original FFTBM [1] it was mentioned that a possible improvement of the method could involve "the development of the procedure taking into account the information represented by the phase spectrum of the Fast Fourier Transform in the evaluation of accuracy". As we can see from Eq. (5), the difference signal magnitude inherently includes the magnitude and the phase of the experimental and calculated signal. The finding that both the magnitude and the phase of the experimental and calculated frequency spectra are contained in AA by making the Fourier transform of the difference signals is very important as this gives the possibility to compare the shapes of the signals. We agree with the authors of [10] that it is difficult to imagine which information is contained in the phase spectrum of the difference signal, since the experimental and calculated phase cannot be simply added. Therefore, to the authors opinion the difference signal phase information is not applicable for the comparison of two signals. In the following, AA will be referred to
as AAMj, since it contains magnitude M and phase j information.
In the FFTBM package there is an option of time shifting of the data trends to analyze separately the effects of delayed or anticipated code predictions concerning some particular phenomena or systems interventions. It is a Fourier transform property that a shift in the time domain corresponds to a change in the phase. This property was used to identify the signals which differ in the time shift. Namely, the magnitudes of such signals are the same and only their phases are different. Therefore, the following expression, not taking into account the phase, is proposed:
M _ n=0
AAM =
where:
2
E |F(-\~Jfn\
(6)

Fexp(fn\ - \FcJfn\ (Re(~pO)2 +(lm(FpO)2 (Re( KO) + (lm(Fca/a):
(7)
= |Mj - M2| = 7(Mj - M2 )2.
When j =j2, the Eq. (5) is equal to Eq. (7). This
means that expression AAM is a measure containing information from magnitudes M only. It is known that when two signals are only time shifted, the magnitude spectra are the same and the value of AAM is consequently zero. It is very unlikely that a calculated signal, which is not shifted, would have a shape giving the same magnitudes as the experimental signal, as the predictions are required to be qualitatively correct. Therefore, AAM can be used to establish the value by which AAMj ° AA is increased due to the time shift contribution. In AAMj, the information from both, the shape of the time domain signal and the time shift, is provided, while in AAM , only the time invariant information of the time domain signal is provided, what can be regarded to a certain degree as the shape of the time domain signal. Therefore, the difference AAj = AAMj - AAm gives the information about the time shift contribution. This difference is further normalized to:
I =
AAMj - AAm AAj
AAM
AAM
(8)
where the indicator I tells how the compared time signals are shifted, and is therefore called the time shift indicator. The larger the value of the time shift indicator I, the larger is the contribution of the time shift to AAMj of the difference signal. A large value of I (I > 1) indicates that the compared signals are maybe shifted in time. When I > 2, we can be quite confident into time shift.
2.3 Calculation of AA
M
As for the AAM calculation, the sum is made from the difference of the magnitudes, there was a need to
m
increase the number of samples in the frequency domain (to make the frequency spectrum more continuous) to correctly calculate the magnitude differences. This can be done by padding the time signal with zeros before taking the discrete Fourier transform [10]. In the FFTBM program [9], the zeros can be padded to make the signal length 16-times the original time signal length. It was verified that when the ratio of the samples between the padded and the non-padded signal is 4 or higher, the amplitude spectrum is properly described.
In order not to distort the amplitude spectrum, the last point in the time domain should be zero. In the opposite, the spectrum of the step function would be added to the spectrum of the original signal. In the FFTBM program this was solved so that the last value of the signal was subtracted from the original signal. In terms of FFT this only changes the mean value of the signal while the frequency spectrum remains the same. After the transformation the mean value was corrected with the subtracted value. AA remains practically the same, using the original number of points or the increased number of points.
2.4 Other extensions
To be able to make a reliable accuracy assessment, it is important that the compared signals are presented appropriately already in the time domain, since it is very difficult to make any adjustments in the frequency domain. For example, if a time shift is needed it should be done in the time domain. In the case of numerous oscillations, AA is very high as normally the oscillations of the experimental and the calculated signals are not in phase. If the oscillations are due to measuring reasons (noise), they may be filtered. In the original FFTBM, there is a parameter, called cut frequency, which filters the spurious contribution of higher frequency components from the accuracy measure in the frequency domain. However, the selection of the cut frequency may be subjective. Our goal is to make the quantitative assessment independent of the cut frequency, i.e. to consider the whole spectrum. Therefore, it was proposed to use the moving average of the time domain signal, which is equivalent to low pass filtering. The advantage is that the user can see the result of smoothing the signal, i.e. eliminating noise from the signal. When we are interested in the accuracy of temperatures, it makes sense always to treat all the signals by the same unit, as using different temperature units (K, °C) produces different results. When the nature of a phenomenon is logarithmic (logarithmic plots), it makes sense to calculate the logarithm of the time signals before making the difference signal. In the next section it is shown on some examples from the literature how these extensions in general lead to a more objective judgement with FFTBM.
3 Results
In validation of the extended FFTBM, results of the international studies under Organisation for Economic Co-operation and Development (OECD) or International Atomic Energy Agency (IAEA) were used.
The first two examples are taken from the LOFT L2-5 test calculations, performed in the frame of OECD BEMUSE (Best-Estimate Methods Uncertainty and Sensitivity Evaluation) programme. The first example is the rod surface temperature, the most important safety parameter. In Figure 1, the UPI calculation, for which a possible time shift (/=1.43) was detected among 14 calculations, is shown. For all plots of all calculations the reader can refer to Figure 24 of the BEMUSE Phase II final report [11]. This time shift can be easily verified by visual observation. It is around 1 s. When shifting the signal, the accuracy measure AA changes from 0.324 to 0.214, indicating a very good agreement. The duration of a heatup was correctly predicted, while the rod started to heatup a little bit too early. From the safety point, the maximum temperature and duration of heatup are more important than the exact heatup starting time.
1200
400
20
40
60 Time (s)
80
100
120
Figure 1. Rod surface temperature for the LOFT L2-5-test.
The second example is break flow. Figure 2 shows the break flow rate in the cold leg for the first 20 s. It can be seen that the measured flow is delayed compared to calculations.
800
600
et400
• 200
	-EXP GID GRS KAERI -•- KFKI	-
		
		
		
10 Time (s)
15
Figure 2. Break flow rate for the LOFT L2-5-test.
0
0
0
5
As can be seen from Table 1, for the above shown four calculations (GID, GRS, KAERI and KFKI) a possible time shift was indicated. Further investigation showed that it was claimed that in the experiment the break opened in 0 s (Table 9 of ref. [11]), while from the measured data it can be concluded that the break opened at 0.3 s. It seems that the participants were not aware of this fact, but this time shift has a significant influence on the quantitative results. The indicator for time shift indications may help to discover such inconsistencies. For the time shift indication all frequencies were used while in the case of the BEMUSE study only half of the frequencies were considered. This is equivalent to some sort of moving average, leading to lower AA values (last column in Table 1).
Table 1. Calculated AA for time window (0 - 100 s)
ID	AA	I	AA*
CEA	0.407	0.5	0.302
GID	0.664	3.8	0.457
GRS	1.160	2.0	1.065
IRSN	0.573	0.6	0.498
JNES	0.770	0.4	0.716
KAERI	0.661	2.9	0.444
KFKI	0.627	2.1	0.443
KINS	0.353	0.4	0.238
NRI-K	0.351	0.3	0.236
NRI-M	0.446	0.3	0.347
PSI	0.487	0.4	0.394
TAEK	0.420	0.2	0.316
UPC	0.748	0.5	0.648
UPI	0.478	0.3	0.384
g 0.02
0
0
5000
10000 15000 Time (s)
Figure 3. H2 mass flow rate for the Phebus FPT1 test.
20000
Without considering the time shift, the calculated AA was very poor (AA=1.095). When considering the time shift of 281 s, the calculated AA was 0.442, indicating a good calculation.
The fourth example is taken from the IAEA study on heavy water reactors [8]. In the frame of quantitative code assessment, there was a discussion on primary pressure. Figure 4 shows a comparison between the experiment and the best calculated pressure. The calculation seems good; however the quantitative analysis did not confirm this. Namely, the pressure criterion 0.1 for AA was not satisfied (AA=0.117).
0
100
300
400
cut frequency used as in BEMUSE study (see Table 14 of [11])
The third example is taken from the Phebus FPT1 test, selected for the OECD international standard problem no. 46 (ISP-46). Figure 3 shows the experimental and calculated hydrogen mass flow rate which was judged as one of the best calculations in the frame of the ISP-46 experiment study [12].
0.12 0.1
3
is 0.08
J 0.06
s
!S 0.04
200 Time (s)
Figure 4. Header 7 pressure for the RD-14M test.
Figure 4 also shows that the experimental pressure signal has several small amplitude oscillations (noise). In the past, such problems were not encountered as the original FFTBM [1] was limited to 1000 data points and reducing the data has a similar effect as smoothing the data. When moving average was used for the experimental signal (curve EXP-ma), the pressure criterion 0.1 for AA was satisfied (AA=0.079). The advantage of using the moving average is that all frequencies are used and the user does not need to evaluate the higher frequency component contribution, what can be more demanding and time consuming than using the moving average. Also, the user can plot the smoothed signal and compare it to the original signal in order to check that the trend of the smoothed signal is correct.
The last example shows the Iodine (I) release at the end of the OECD ISP-46, Phebus FPT1 test. In Figure 5, the I deposition is shown as a function of the core level. The measured spikes indicate elevations of grids, which prevented the measurements. When the plots are made in a logarithmic scale, it seems reasonable to compare the logarithms of the signals. In the opposite case the accuracy would be judged very poor. Besides, the moving average was used to eliminate the contribution of spikes and the so calculated AA was 0.40, indicating a good calculation. Namely, for parameters other than pressure, the criterion for acceptable prediction is AA below 0.7. Without
considering the moving average for the experimental signal the obtained AA was 0.80.
-2.5	-1.5	-0.5	0.5
Level (m)
Figure 5. Deposition of I for the Phebus FPT1 test.
4	Conclusions
In the paper, extensions of FFTBM, intended for quantitative assessment of thermal-hydraulic code calculations, are presented. A new measure for indication of the time shift between the experimental and the calculated signal is proposed. It is also suggested to make all operations in the time domain, as it is very difficult to make adjustments in the frequency domain (e.g. logarithmic scale, moving average). For the validation of the extended FFTBM, results from OECD and IAEA international studies have been used. They show that considering the extensions, the quantitative assessment is more objective and helps the analyst to get an insight into the results when several calculations are performed by different participants. Finally, the quantitative results are compared with the evolutionary conclusions of the expert panel evaluating the calculated results of the considered test and the quality of the codes used. They agree well and with the extended FFTBM a few disagreements, resulting from improper use of FFTBM, are resolved.
Acknowledgements
The authors acknowledge the support of the Ministry of higher education, science and technology of the republic of Slovenia within the P2-0026 program.
5	References
[1]	W. Ambrosini, R. Bovalini, F. D'Auria, Evaluation, of accuracy of thermal-hydraulic code calculations, Energia Nucleare, 1990, Vol. 7, pp. 5-16.
[2]	B. Mavko, A. Prosek, F. D'Auria, Determination of code accuracy in predicting small-break LOCA experiment, Nuclear Technology, 1997, Vol. 120, pp. 1-19.
[3]	A. Prosek, F. D'Auria, B. Mavko, Review of quantitative accuracy assessments with fast Fourier transform based method (FFTBM), Nuclear Engineering and Design, 2002, Vol. 217, 1-2, pp. 179-206.
[4]	F. D'Auria, W. Giannotti, Development of a code with the capability of internal assessment of uncertainty, Nuclear Technology, 2000, Vol. 131, pp. 159-196.
[5]	R.F. Kunz, G.F. Kasamala, J.H. Mahaffy, C.J. Murray, On the automated assessment of nuclear reactor systems code accuracy, Nuclear Engineering and Design, Vol. 211, PP. 245-272, 2002.
[6]	A. Prošek, B. Mavko, B., A tool for quantitative assessment of code calculations with an improved fast Fourier transform based method, Electrotechnical Review, 2003, Vol. 70(5), pp. 291-296. Available: http://ev.fe.uni-lj.si/5-2003/prosek.pdf.
[7]	A. Prošek, B. Kvizda, B. Mavko, T. Kliment, Quantitative assessment of MPC trip transient in a VVER, Nuclear Engineering and Design, 2004, Vol. 227, pp. 85-96.
[8]	A. Prošek, F. D'Auria, D. J. Richards, B. Mavko, Quantitative assessment of thermal-hydraulic codes used for heavy water reactor calculations, Nuclear Engineering and Design, 2006, Vol. 236, pp. 295-308.
[9]	A. Prošek, M. Leskovar, B. Mavko, Quantitative assessment with improved fast Fourier transform based method by signal mirroring, Nuclear Engineering and Design, 2008, Vol. 238(10), pp. 2668-2677.
[10]	Smith S. W., The Scientist and Engineer's Guide to Digital Signal Processing, Second Edition, 1999, California Technical Publishing San Diego, California.
[11]	BEMUSE PHASE II REPORT: Re-Analysis of the ISP-13 Exercise, Post Test Analysis of the LOFT L2-5 Test Calculation, NEA/CSNI/R(2006)2. Available at http://www.nea.fr/html/nsd/docs/2006/csni-r2006-2.pdf.
[12]	M. Leskovar, B. Mavko, Simulation of the Phebus FPT1 severe accident experiment with the MELCOR computer code, Stroj. vestn., 2006, Vol. 52(3), pp. 142-160.
Andrej Prošek received his diploma degree from the Faculty of Electrical Engineering, University of Ljubljana in 1987. The same year he joined the Jožef Stefan Institute (JSI), Reactor Engineering Division and began his postgraduate studies. He received his M.Sc. and Ph.D. degrees in nuclear engineering from the Faculty of Mathematics and Physics, University of Ljubljana, in 1992 and 1999, respectively. Currently he is senior research associate at JSI. His research interests include nuclear safety, thermal-hydraulic (TH) safety analyses, uncertainty evaluation and accuracy quantification of TH codes.
Matjaž Leskovar received his diploma degree from the Faculty of Mathematics and Physics, University of Ljubljana, in 1993. The same year he joined the Jožef Stefan Institute (JSI), Reactor Engineering Division and began his postgraduate studies. He received his M.Sc. and Ph.D. degrees in nuclear engineering from the Faculty of Mathematics and Physics, University of Ljubljana, in 1996 and 2000, respectively. Currently he is senior research associate at JSI. His research interests include nuclear safety and severe accidents, focusing on the modelling of fuel-coolant interactions.