Development of a new method to determine the axial void velocity profile in BWRs from measurements of the in-core neutron noise

Determination of the local void fraction in BWRs from in-core neutron noise measurements requires the knowledge of the axial velocity of the void. The purpose of this paper is to revisit the problem of determining the axial void velocity profile from the transit times of the void between axially placed detectors, determined from in-core neutron noise measurements. In order to determine a realistic velocity profile which shows an inflection point and hence has to be at least a third order polynomial, one needs four transit times and hence five in-core detectors at various axial elevations, whereas the standard instrumentation usually consists only of four in-core detectors. Attempts to determine a fourth transit time by adding a TIP detector to the existing four LPRMs and cross-correlate it with any of the LPRMs have been unsuccessful so far. In this paper we thus propose another approach, where the TIP detector is only used for the determination of the axial position of the onset of boiling. By this approach it is sufficient to use only three transit times. Moreover, with another parametrisation of the velocity profile, it is possible to reconstruct the velocity profile even without knowing the onset point of boiling, in which case the TIP is not needed, although at the expense of a less flexible modelling of the velocity profile. In the paper the principles are presented, and the strategy is demonstrated by concrete examples, with a comparison of the performance of the two different ways of modelling the velocity profile. The method is tested also on velocity profiles supplied by system codes, as well as on transit times from neutron noise measurements.


Introduction
Ever since early work in the mid-70's on the in-core neutron noise in BWRs revealed that direct information on the local two-phase flow fluctuations can be obtained through the local component of the neutron noise (Wach and Kosály, 1974;Behringer et al., 1979), it was thought that such measurements could also be used to determine the (axially) local void fraction in the core. Such attempts were also made quite early by putting forward suggestions on how to extract the local void fraction from in-core neutron noise measurements (Kosály et al., 1975;Kosály, 1980). However, as described also recently in Hursin et al. (2020), the suggested methods were either incomplete and required either calibration, or several auxiliary conditions, whose fulfilment was unclear and rather uncertain, or both. Hence, to date no routine method exists for extracting the local void fraction in BWRs from measurements.
calculations which do not require detailed knowledge of the two-phase flow structure.
As known from previous work, by using the signals of two detectors in the same string, one can determine the transit time of the twophase flow between the two detectors. The transit time can be obtained either from the maximum of the cross-correlation function (CCF) of the two signals, or from the slope of the phase of the cross-power spectral density (CPSD), as a function of the frequency (Pázsit and Demazière, 2010;Hursin et al., 2020). One further possibility is to use the frequencies where the phase is equal to zero or ± (Pázsit and Glöckler, 1994;Pázsit and Demazière, 2010;Yamamoto and Sakamoto, 2021), or the impulse response function (Czibók et al., 2003;Hursin et al., 2020). It is generally assumed that since the neutron noise is induced by the void fluctuations, the transit time is related to the steam velocity (this point will be returned to below). Hence, in principle, there is a possibility to determine the void velocity at the detector positions from the transit times.
However, this in itself is not sufficient to determine the void fraction. To make progress, in a series of works we investigated how the knowledge of the void velocity at the detector positions could be utilised to determine the void fraction (Pázsit et al., 2011;Dykin and Pázsit, 2013;Dykin et al., 2014). One possibility was to use a relationship between the void fraction and the void velocity through mass conservation and a known slip ratio. This method has the disadvantage that is not based on purely neutronic measurements, rather it assumes knowledge of the flow properties. The other possibility was based on the fact that the neutron noise is induced by the passage of the fluctuating two-phase flow structure through the so-called detector field of view, determined by the range of the local component of the neutron noise (Wach and Kosály, 1974;Dykin and Pázsit, 2013). This range is determined by an exponent ( ), which is a spatial decay constant (with dimensions of inverse length), describing the spatial decay of the local component. It is a function of the void fraction (hence also of the axial elevation ), but independent from all other thermal hydraulic parameters such as the flow regime or the slip ratio. For any given void fraction, it can be calculated by reactor physics methods. Due to the existence of the local component, the auto power spectral density (APSD) of the neutron detectors will then have a break frequency at An illustration of the break frequency between 10 and 15 Hz is shown in Fig. 1 Since the break frequency can be obtained in a straightforward way from the detector signal, then, if also the void velocity can be extracted from the neutron noise measurements, the range ( ) of the local component can also be obtained and, from the calculated correlation between the void fraction and , the void fraction can also be extracted. Although this method is not purely based on measurements, nevertheless similarly to the neutron energy spectrum method, it has the advantage that it does not require any knowledge on the thermal hydraulic conditions in the core.
A pilot study on the feasibility of both methods (mass conservation and break frequency methods) was investigated in simplified models. A bubbly flow was generated through a Monte-Carlo simulation, with the help of which the performance of the two methods could be investigated in model cases (Pázsit et al., 2011;Dykin and Pázsit, 2013). Regarding the first method, a slip ratio equal to unity was assumed. This is not valid in realistic cases, it was just used for test purposes. Regarding the break frequency method, the dependence of the range of the local component on the void fraction was calculated in a simple analytical model. Similar numerical studies were performed also by other groups with the purpose of investigating the possibilities of determining the transit time from in-core neutron noise measurements, as well as to calculate the local component of the neutron noise in realistic cases (Yamamoto, Fig. 1. APSDs of simulated in-core detector signals induced by a bubbly flow at various elevations (Pázsit et al., 2011). 2014; Yamamoto and Sakamoto, 2016;Yamamoto, 2018;Yamamoto and Sakamoto, 2021). What regards the flow simulation, even in these works a bubbly flow was simulated. However, the local component of the neutron noise, or rather the full neutronic transfer function, was calculated with a frequency dependent Monte-Carlo method with complex weights. Therefore, the range of the local component is fully realistic in these calculations, and can be used in practical applications, such as the evaluation of pilot experiments in research reactors (Yamamoto and Sakamoto, 2021).
A note on the usage of the word ''void velocity'' is in order here. In the neutron noise community it is tacitly assumed that the transit time deduced from the neutron noise measurements corresponds to the transit time of the void (steam). However, the neutron noise is induced by the temporal fluctuations of the reactor material around its mean value. In a binary or dichotomic medium (fluid-void), it is represented by the fluctuations of the minority component. At low void fraction, such as a sparse bubbly flow, the neutron noise is indeed generated by the fluctuations represented by the void, and hence the transit time obtained by the noise measurement corresponds to the transit time of the void. At high void fractions, the fluid becomes the minority component, hence the neutron noise is generated by the water droplets/mist, and/or the propagating surface waves of the water in an annular flow regime. Therefore, in the upper part of the core, i.e. between the uppermost two detectors, it is more correct to talk about the transit time of the perturbation.
It is a great advantage of the break frequency method that the break frequency depends on the transit time of the perturbation through the detector field of view. Hence, it is completely independent of whether the fluctuations of the void or the fluid generate the detected neutron noise. This means that the range of the local component will be determined correctly, and hence also the break frequency method will supply the void fraction correctly, since this latter is extracted from the dependence of on the void fraction. Because of this fact, and because full realistic calculations of the detector field of view are underway (Yamamoto, 2018), the break frequency method appears to be more suitable and effective for determining the void fraction in operating BWRs.
The above also means that it would be more correct to refer to ''perturbation velocity'' rather than ''void velocity'' in this paper. However, this would be cumbersome and even confusing, and for practical reasons, we will use the terminology ''void velocity'' or ''steam velocity'' throughout, on the understanding that in the upper part of the core, it actually means the velocity of the perturbation, which may differ from that of the void (steam). Whichever method is used for recovering the void fraction, the void/perturbation velocity is needed at the detector positions. Elaboration of a method of how the void velocity can be extracted from the detector signals is the sole objective of this paper. Namely, in incore noise measurements only the transit times of the void between two axially displaced neutron detectors can be obtained. The transit times are integrals of the inverse of the velocity, which is not constant between the detectors. The relationship between the void velocity at the detector positions, and the transit time between the detector pairs, is hence rather involved.
Determination of the void velocity at the detector positions is therefore only possible if a functional form is assumed for the velocity profile, which depends only on a few parameters, which then can be determined from the available transit time data. Since the axial dependence of the velocity has an inflection point, it has to be described by a non-linear function. The simplest such function, which was also suggested by Pázsit et al. (2011) and Dykin and Pázsit (2013), and which is the only one tested so far, is a third order polynomial.
However, a third order polynomial has four parameters. To determine these, one would need four independent transit times, hence access to five detectors. The standard instrumentation of BWRs comprises only 4 detectors axially at one radial core position, thus one has access only to three transit times between the three detector pairs.
To solve this problem, it was suggested that one could use, in addition to the four standard LPRMs (Local Power Range Monitors), an additional TIP detector (Transverse In-core Probe), by placing the TIP at an axial position either between the four LPRMs, or outside these, i.e. in a position different from those of the LPRM positions, and determine the transit time between the TIP and the nearest LPRM. This approach was tried in measurements, performed in the Swedish Ringhals-1 BWR (Dykin et al., 2014). Unfortunately, as is also described in Dykin et al. (2014), the attempt was unsuccessful. Due to the fact that, for obvious reasons, the data acquisition for the LPRMs and the TIP detectors belong to completely separate measurement chains, the data acquisition is made separately, which made synchronising the two data acquisitions with a sufficient accuracy impossible. Hence the transit time between the TIP and any of the LPRMs was not reliable. The conclusion in Dykin et al. (2014) was that the application of the TIP detector for acquiring a fourth transit time is not feasible.
Therefore, here we suggest a different strategy. First, we realise that there is no need for a fourth transit time to determine four parameters of the velocity profile, either a polynomial or some other form, if both the axial point of the onset of the boiling, as well as the steam velocity in this point are known. The onset point of the boiling can be determined with a TIP detector alone, from the amplitude of its root mean square noise (RMS) or its APDS, or, alternatively, from the coherence between the TIP and the lowermost LPRM, if these are determined as a function of the axial position of the TIP. Neither of these requires a perfect, or any, synchronisation between the two data acquisition system. At the onset of the boiling the void velocity can be assumed to be equal to the inlet coolant velocity, which is known. Thus, knowledge of these two quantities reduces the number of unknowns of the axial velocity profile to be determined from four to three.
Second, there exist non-linear functions with an inflection point, which represent an even higher order non-linearity than a third order polynomial, but which nevertheless can be parametrised with only three parameters instead of four. Examples are certain trigonometric or sigmoid functions. For simplicity these profile types will be referred to as 'trigonometric''. In this case not even the onset point of the boiling needs to be known; determination of the void profile is then possible based on solely of the three measured transit times with the standard instrumentation, without the need for using a TIP detector at all.
In the following, the principles, as well as the applicability of both types of velocity profile forms (trigonometric and polynomial) will be investigated in conceptual studies. Various types of velocity profiles, both trigonometric and polynomial, will be assumed as the ''true'' profiles as a starting point. From the true profiles, the three transit times between the four detectors can be calculated, and then the inversion procedure applied and its accuracy investigated. Since we do not have access to measurement data with four transit times or three transit times plus the knowledge of the axial onset point of the boiling, we need to make some assumptions when investigating the applicability of the polynomial form. The sensitivity of the results on the accuracy of these assumptions will be illustrated in model calculations. We will also investigate the flexibility of the two forms of the velocity profiles to reconstruct various types of axial velocity dependences, and the sensitivity of the reconstruction on the correct assumption on the form of the profile (i.e. starting with a trigonometric as ''true'' and performing the reconstruction with the polynomial form, and vice versa), as well as taking a true polynomial profile with a given onset point of the boiling and making the reconstruction with a different onset point as an incorrect guess. A test will be also made with transit times corresponding to void velocity calculations with a thermal hydraulic system code. Finally, an attempt will be made to reconstruct the (unknown) velocities at the detector positions from a measurement at Ringhals-1, both with the trigonometric and the polynomial velocity forms. The focus of the investigation is to see which method can reconstruct the known transit times better, and which inversion method is more robust and convergent.

Characteristics of the velocity profile
As is general in reactor noise diagnostics problems, when only a limited number of measurements is available, obtained from detectors in a few specified spatial positions, this is not sufficient to reconstruct the full spatial dependence of the noise source. Inevitably, one needs to make an assumption on the space dependence of the noise source in an analytical form, which contains only a limited number of free parameters. These can then be determined from the limited number of measurements (Pázsit and Demazière, 2010).
This strategy is easy to follow for localised perturbations, such as a local channel instability or the vibrations of a control rod, since the perturbation can be simplified to a spatial Dirac-delta function, either with a variable strength, or with a variable position. All these cases can be described by a few parameters, whose physical meaning is obvious, and the guess on the analytical form is rather straightforward.
What regards the reconstruction of the velocity profile, the case is more complicated. Here a whole profile (the axial dependence of the velocity) needs to be reconstructed, and it is not obvious how to parametrise it. One main difficulty is that, for obvious reasons, no directly measured velocity profiles are available, which would give a definite hint on a functional form with only a few parameters to be determined. Only qualitative information is known either from calculations with system codes, from common sense considerations, or from simulations.
An inventory of the available knowledge yields the following. What regards results from calculations with system codes, there are some data available from calculations with the system codes TRACE (USNRC, 2008a,b,c) and RAMONA (Wulff, 1984;RAMONA, 2001). Fig. 2 shows a few profiles from calculations with TRACE, where account was taken for the fact that the boiling does not start at the inlet, rather at a higher elevation (Hursin et al., 2017). In Fig. 3, calculations with RAMONA of the steam velocity in Ringhals 1 in a few selected channels are shown. RAMONA5 is capable of treating a non-homogeneous two-phase flow, and the vapour velocity is related to the liquid velocity through the expression is the slip factor, and 0 is the bubble rise velocity (the vapour velocity relative to stagnant liquid) using the notations in the I. Pázsit et al.  RAMONA5 user manual (RAMONA, 2001). The slip parameter is calculated using the option for the Bankoff-Malnes correlation. From the RAMONA result files, the nodal vapour velocity (also referred to as steam velocity) of each channel can be extracted.
In Fig. 3, the discontinuity at around 2.5 m is due to the fact that the fuel assemblies, in which the calculations were made, contain partial length fuel rods, with a different length in one of the channels. At the elevation of the end of the partial length, there is an abrupt change in the void/fuel ratio, hence the sudden change in the void velocity. The effect of such a discontinuity on the proposed method of the velocity reconstruction will be assessed in Section 6.
Another possibility is to use results from simulations of a bubbly flow in a heated channel, which were performed by an in-house Monte Carlo code. This code was developed earlier and was used in previous work (Pázsit et al., 2011;Dykin and Pázsit, 2013). Some profiles, resulting from these simulations, are shown in Fig. 4.
What these figures tell us is that the velocity increases monotonically in the channel from the inlet, first in a quadratic manner, then the increase slows down, either leading to an inflection point, or to a linear increase towards the core exit. One has though to keep in mind that these are all calculated/simulated values, and a direct measurement of the velocity profiles inside a BWR core will never be available.Moreover, as was mentioned in the Introduction, in the upper part of the core, the void velocity can differ from the velocity of the perturbation (which is just as impossible to measure directly). Hence, the simulated/calculated velocity profiles are mostly used to get a hint on the possible axial velocity profiles to be investigated. On the other hand, regarding the void velocity profiles, based on the ever improving reliability and accuracy of the system codes, one may rely on the type of the velocity profiles which these codes predict. In this respect one can draw the conclusion that all the profiles shown above can be satisfactorily approximated by a 3rd order polynomial. Trying to fit data from real measurements to such profiles will reveal which out of the above two types, close to linear increase in the upper part of the core or a strong inflection point is more likely. A fit to the profiles obtained from RAMONA calculations will be made in Section 6. A fit to real measurements will be made in Section 7.

Possible analytical forms
In our previous works (Pázsit et al., 2011;Dykin and Pázsit, 2013;Dykin et al., 2014) a third order polynomial was assumed for the axial dependence of the void velocity: This form has found to have some disadvantageous properties: partly that the integral of −1 ( ) w.r.t. does not exist in an analytical form, and partly that it assumes that the boiling starts at the inlet, i.e. at = 0, which is not true in practical cases. The first of these disadvantages does not represent a significant difficulty, since the unknown parameters -can also be determined by numerical unfolding methods, as it will be shown below. Even for the trigonometrical profile, where the same integral exists in analytical form, the numerical unfolding method is more effective than root finding of a highly not-transcendental analytical function in several variables.
The second property poses somewhat larger problems. Accounting for the fact that the onset of the boiling is at = ℎ where ℎ is an unknown, would increase the number of parameters to be determined to 5. However, as suggested in this work, if the onset point = ℎ of the boiling is known from measurements, then the third order polynomial form of (2) can be written in the form where ( ) is the unit step function, and 0 is the (known) inlet coolant velocity. This form contains only three unknowns, which can be determined from the three measured transit times. This procedure is suggested for future use, such that the onset point of boiling is determined by measurements with movable TIP detectors. Since such measurements are not available at this point, the test of the polynomial form will be performed by making a qualified guess on the onset point. The uncertainty of the unfolding procedure with a polynomial profile can be assessed with respect to the error in the estimation of the position of the onset point of boiling, which will be performed in Section 5.2.2.
In addition we propose also to investigate another path. The essence is the recognition that there exist non-linear functions other than a third-order polynomial which have an inflection point, and which contain only three free adjustable parameters. These include trigonometric functions, such as ⋅ atan ( ( − )), where , and are constants, or the so-called ''sigmoid'' function, used in the training of artificial neural networks (ANNs). In the continuation we will refer to such profiles as ''trigonometric''. For such profiles the onset point of the boiling does not need to be known. In the next section such a model is proposed, and a procedure for its use for the unfolding of the velocity profile is suggested.
Of course, the price one has to pay for the convenience of only needing to fit three parameters instead of four is that the structure of the profile is more ''rigid'' than that of the more general polynomial form, hence its flexibility of modelling and reconstructing a wide range of velocity profiles is reduced as compared to the polynomial fitting. If the onset point of boiling was known, then clearly the polynomial I. Pázsit et al. profile would be recommended. It the onset point is not known, it is not clear whether the use of a trigonometric form, or that of the polynomial form used with a guess for the axial point of the onset of the boiling yields better results.

Construction of a simple non-polynomial velocity profile
In order to obtain a velocity profile with an inflection point, which can be described by only three parameters, we shall assume a very simple phenomenological model based on simple considerations. The model does not have any deep physical meaning, or justification. One of its advantages, besides its simplicity, is that since it is based on a physical model, whatever coarse it is, it makes it simpler to estimate the possible range of the model parameters (which is useful in the inversion process), and in particular it is more straightforward to find initial guesses of the parameters included to the numerical inversion procedure than for the polynomial model. Although, the comparative investigations made later on in this chapter will show that this latter advantage is not significant in the sense that the polynomial model is much less sensitive to the correct choice of the starting guess of the sought parameters than the non-polynomial model.Because of its simplicity, such a model of assuming the void fraction being proportional to the integral of the heat generation rate was used also in other works (Yamamoto and Sakamoto, 2016).
Assume that the core boundaries lie between = 0 and = in the axial direction with a static flux ( ). Assuming that the boiling starts at the axial elevation = ℎ, and that there is a simple monotonic relationship between void fraction and void velocity, and that the latter at point is proportional to the accumulated heat production between the boiling onset and the actual position, gives the form where ℎ and are unknown constants. 1 Assume now, for simplicity, a simple cosine flux shape as In reality, the axial flux shape in a BWR deviates quite appreciably from a cosine-shaped profile, and moreover that profile is known from in-core fuel management calculations. Hence, the assumption of the simple cosine flux profile could be replaced with a more realistic one, although presumably at the price that the simplicity of the model, and hence its advantages, would be lost. In Eq. (5) it is not assumed that = ∕ , rather is kept as an independent (unknown) parameter. By allowing < ∕ , the effect of the reflector, represented as an extrapolation length as an independent parameter, can be accounted for.
With this choice, after integration, the velocity profile is obtained in the simple form with and 1 = A qualitative illustration of a typical velocity profile provided by this model, referred to as the trigonometric profile, is given below. To this order, geometrical as well as inlet and outlet velocity data are taken from the Ringhals-1 plant. The geometrical arrangement is depicted on Fig. 5 In the next section, the unfolding procedure (the algorithm for the reconstruction of the velocity profile from the transit times) will be briefly described. The velocity reconstruction method will then be first tested on various trigonometric profiles supplied by the above model, together with those of the polynomial profile (referred to as synthetic velocity profiles). Thereafter the reconstruction of the velocity profiles will be tested on the data given by calculations with RAMONA, shown in Fig. 3. Finally, an attempt will be made to reconstruct the velocity profile from a Ringhals measurement.

The unfolding procedure
First we tried to use the velocity profile given in Eq. (6), since it depends only on three parameters, hence three transit times, derived from four LPRM signals, should be sufficient for reconstructing the velocity profile. Eq. (6) has the further property that its inverse is analytically integrable, thereby giving a possibility to express the transit time 1,2 of the void between the detector positions 1 and 2 , with 1 < 2 , as analytical functions of the unknown parameters 1 , 1 and . For practical reasons we will number the detector positions such that 1 corresponds to the lowermost detector, LPRM 4, and the transit times between the detector pairs will be indexed by the position of the lower detector, i.e. 1,2 ≡ 1 etc. With these notations, one has

I. Pázsit et al.
Our expectation was that in possession of the analytical expressions for ( 1 , 1 , ), = 1, 2, 3 in the above form, and having access to given values of the three measured transit times , = 1, 2, 3, the unknown parameters 1 , 1 , can be determined as the roots of the non-linear equation system This strategy was tested by choosing detector positions, core size, as well as inlet coolant velocity and the same value for 1 which were used in calculating the profile in the right hand side of Fig. 6. Having the analytical form of ( ), the concrete transit times , = 1, 2, 3 can be numerically evaluated and used in (9), with the given in the analytical form (8). For the numerical solution of this non-linear equation system, the numerical root finding routine NSolve of Mathematica 12.1.1.0 was used (Wolfram Research). However, the root finding did not converge, even if quite accurate starting values were specified. It appears that the NSolve routine is primarily designed for treating polynomial equations, rather than transcendental ones. Therefore, another path was followed to unfold the parameters of the void profile from the transit times. Instead of using Nsolve, a kind of fitting procedure was selected by searching for the minimum of the penalty function as functions of 1 , 1 and . First the FindMinimum routine of Mathematica, was used. This procedure worked well and was able to reproduce the input parameters of the velocity profile. Initially the analytical form (8) was used for the ( 1 , 1 , ). However, it turned out that defining these latter as numerical integrals with free parameters 1 , 1 , worked much faster and with better precision, showing also that for the unfolding, it is not necessary that the transit times are given in an analytical form. Consequently, the modified polynomial form of ( ) in Eq. (3) can also be used, despite that −1 ( ) is not integrable analytically.
The unfolding procedure was tested using both the trigonometric velocity profile given in (6), as well as with the polynomial profile of Eq. (3). Tests were made with various values of the parameters, also with combinations that yielded velocity profiles similar to those in Fig. 2. These extended numerical tests were made by Matlab. The minimum of the penalty function (10) was found by own MATLAB scripts. In addition, unlike for the case with Mathematica, for the trigonometrical profile, using the parameter values obtained from the minimisation process as starting values to the routine fsolve helped to successfully solve also the nonlinear system of Eqs. (8), to get the velocity profile.

Trigonometric profile
Two tests will be shown for illustration, with two different profiles. We used a more curled and flatter profile, respective, with the following  The true (=starting) and reconstructed profiles for these two cases are shown in Fig. 7. The solid red line represents the true profile, and the broken blue the reconstructed one. It is seen that the inversion algorithm was able to reconstruct the original profiles in both cases quite well.
Tests made on a large variety of different profiles revealed that finding the minimum of the penalty function, Eq. (10), with the Matlab routine fsolve, the procedure in some cases did not converge to the true parameters. In some cases the minimum searching ended up by providing complex numbers for the searched parameters, even if quite accurate starting values and searching domains were specified. This lack of convergence is a reason for concern, since in a real application one does not know the searched parameters and hence cannot specify good starting values.
However, the fact of sometimes obtaining complex values of 1 , 1 and gave the idea of taking only the real part of the search function, such that the minimisation was performed on the modified penalty function With this, the convergence problems experienced previously ceased, and in all cases the minimisation procedure found the correct parameters for the reconstruction of the initial profile. An illustration of the performance of the method with a large selection of different profiles is shown in 8.

Polynomial profile
For the tests with the polynomial profile, the form (3) was used. This form has four fitting parameters. It was tested in two different ways. First, we assumed that the correct axial position ℎ of the onset of boiling is known (e.g. from a tip measurement). In that case, there are only the three parameters , and to be fitted. Second, we assumed that the correct value of ℎ is not known, rather it was guessed incorrectly, with a certain error. The interesting question was then to see how large an error this incorrect estimate causes in the reconstruction process.

Reconstruction with a known boiling onset point ℎ
In this case the unfolding worked always correctly and promptly, without the need of taking the real value of the penalty function. Finding correct initial values of the parameters for the minimisation process was easy, by taking a qualified guess of the void velocity at the outlet, the velocity gradient at the axial position of the inflection point of the profile and the void velocity at the position of the second detector. It seemed that handling the polynomial profile was more efficient than that of the trigonometric profile. One case of a successful reconstruction is shown in Fig. 9, where the following data were used: H = 3.68 m; h = 0.2 m; 0 = 2 m/s; ( = ) = 12.1 m/s. Fig. 9 that the reconstruction is completely successful.

Reconstruction with an unknown boiling onset point ℎ
In this case we assumed a value for the onset point in the reconstruction procedure which was different from the true one. In a practical case, when no information on the boiling onset point is available, it is a reasonable choice to assume the position at the boiling onset halfway between the core inlet and the position of the lowermost detector, because this minimises the error of the guess. Since the lowermost detector position is at 0.66 m, we selected ℎ = 0.33 m. Two reconstructions were made, one by taking ℎ = 0.45 m for the true onset point, and another by taking ℎ = 0.15 m for the true onset point.
The results of the reconstruction are seen in Fig. 10. It is seen that, as expected, the reconstruction will not be perfect, especially in the lower section of the core. However, as it is also seen in the figure, the only difference between the true and the reconstructed profiles is at the lowermost part of the core, and the incorrect reconstruction affects I. Pázsit et al.

Significance of choosing the right type of profile
One might also be interested to know the significance of choosing the right type of profile. In other words, to check the performance of the reconstruction procedure when the true profile is trigonometric, and the reconstruction is attempted by using a polynomial form, and vice versa.
The results of such a test are shown in Fig. 11. In the left hand side figure the true profile is trigonometric, whereas the reconstruction is made by the assumption of a polynomial form. In the right hand side figure the opposite case is shown, i.e. when the true profile is polynomial, whereas the reconstruction is made by the assumption of a trigonometric form.
It is seen that the use of the polynomial profile is more flexible than that of the trigonometric profile. It can very well reconstruct a true trigonometric profile throughout the whole axial range. It has though to be added, that here it was assumed that the onset point of the boiling was known. The figure also shows that when the true profile is polynomial, the trigonometric form has a slight error in the reconstruction both at low and at high elevations.
The overall conclusion of these model tests is that use of the polynomial profile is preferred to be used in the reconstruction over the simpler trigonometric profile.

Test with velocity profiles obtained from RAMONA
Another test of the unfolding procedure can be made by using the velocity profiles generated by the RAMONA calculations, shown in Fig. 3. This is an interesting exercise, even if, as mentioned earlier, the true void velocity does not agree with the velocity of the perturbation in the higher upper part of the core, because it uses non-analytical (non-synthetic) profiles, rather numerically calculated ones. It also represents a challenge, due to the discontinuity in the velocity profiles, which arises because of the presence of partial length rods. In such a case one can count on that the velocity of the perturbation will also be affected the same way, i.e. it will experience a discontinuity at the top of the partial-length rod. Hence this exercise will give information on the possibilities or the reconstruction of the velocity profile for such cases.
For the test, first the transit times between the detector positions had to be determined. Since the RAMONA calculations give the velocity in a number of discrete points (26 positions), for the accurate determination of the transit time, first a piece-wise continuous function was fitted to the calculated profiles. From the core inlet up to the lower end of the discontinuity, as well as from the top end of the discontinuity to the core outlet, a polynomial fit was made. The discontinuity, which occurs between two adjacent RAMONA points, was represented by a linear fit. The result of this fitting is shown in Fig. 12.
Thereafter, the transit times were calculated by an integration of the inverse velocity from the fitted curves, and used in the unfolding procedure. Due to its larger flexibility, a polynomial fit was used. The onset point of the boiling, and the steam velocity at this point was taken from the RAMONA data. The results of the reconstructed profiles are shown in Fig. 13, and the reconstructed velocities at the detector positions are listed in Table 1. Fig. 13 shows that, for trivial reasons, the reconstructed profiles cannot display any discontinuity. However, they reconstruct the RAMONA velocity profiles quite accurately up to the discontinuity, after which there is a significant deviation between the true and the reconstructed values. The reconstructed velocity in this section, i.e. in the uppermost part of the core, overestimates the true velocity. Accordingly, the steam velocity values are reproduced quite accurately in the lower three detectors, whereas there is an error between 5%-10% in the uppermost detector (Table 1). Regarding this latter detector, one has to add that it is quite close to the discontinuity, which means an abrupt change in    the velocities of all phases (fluid and steam). In such a case the concept of ''local velocity'' and ''local void fraction'' becomes problematic, so reproducing the local void fraction in that position is not a prime priority. In view of the above, it is quite encouraging that despite the discontinuous character of the velocity profile, the true velocity values were correctly reproduced at 3 of the 4 detectors, and with an overestimation of the true velocity by only 5%-10% in the uppermost detector. It can be added that, as is seen in Eq. (1), an overestimation of the velocity leads to an underestimation of the detector field of view . Due to the inverse relationship between the field of view and the void fraction, this also means an overestimation of the void content. This way, one can claim that in cores containing partial length fuel with characteristic length up to the uppermost detector, the reconstruction procedure yields a correct value for three lower detectors, and supplies an upper limit on the void fraction at the position of the uppermost detector.

Test with Ringhals-1 data
It might be interesting to test the procedure with pure measurement data, where the true values of the flow profile parameters are not known. This has the disadvantage, that in such a case the validity of the reconstructed velocity profile cannot be verified, but it is a test of whether the unfolding procedure works when one cannot give a qualified guess of the starting values for the search of the minimum. To this end we took real measurement data from Ringhals-1 (Dykin et al., 2014). In this particular measurement campaign, four identical measurements were taken, while a TIP detector was placed at the four LPRM positions, respectively. Since the position of the TIP does not influence the thermal hydraulic conditions, the four transit times, obtained from the fitting of the phase curves, can be used for a rough estimate of the uncertainty of the transit time estimation. The three transit times for the four measurements, together with the mean values and the relative standard deviations are given below in Table 2. It is Table 2 Transit times from Ringhals-1 (from Dykin et al. (2014)). All times are in seconds. seen that the uncertainty of the transit time estimation is about 1%. For the velocity reconstruction, the mean value of the transit times was used.

Measurement number
Since in this case neither the true character of the profile, nor the values of the corresponding parameters are known, the only assurance of the successful reconstruction is that the reconstructed values at least reproduce the transit times properly. One could expect that the task is underdetermined, i.e. that several void velocity profiles can be constructed which all reproduce the proper transit times, but are otherwise different, and supply therefore different values for the velocities at the detector positions.
Both the trigonometric and the polynomial forms were used in the attempt of reconstructing the velocity profile. For the initialisation, the parameters 0 = 2 m/s and ℎ = 0.3 m were used. It was seen that the profiles, either trigonometric or polynomial, which were able to reconstruct the measured transit times, resembled much more to the TRACE simulations in Fig. 2 without a marked inflection point, rather than to the more ''curved'' profiles in Figs. 3-4. The reconstructed profile, which yielded the best agreement with the measured transit times, is shown in Fig. 14. There is no calculation of the void velocity by either TRACE or RAMONA available for this measurement, and moreover it is not practical either, in view of the difference between the void velocity and the velocity of the perturbations as mentioned before. Hence a comparison between the profile reconstructed from the measurement with curve fitting, to the simulated profile from a system code, is not practical.
When comparing the performance of the two methods, i.e. the trigonometric vs the polynomial profile, it was once again found that assumption of the polynomial profile in the reconstruction performed better. This is because it has more free parameters that can be fitted, and hence this model is more flexible. Finding the proper parameter values which must be fixed for the search for the minimum of the penalty function took more trial and error, but also it made possible to find a better fit in the end.
Thus it turned out that the original concern that the case is underdetermined and one may obtain multiple solutions, was not valid for this case. This is not a proof that this should be the case in all other measurements, but at least it is reassuring. Significantly more cases need to be investigated to get a confirmation of the validity of the procedure, and validation against calculated/simulated values is desired. Unfortunately, there is no possibility to validate the method against explicit measurements of void velocity profiles. However, there is a database of noise measurements available, made in Swedish BWRs by GSE Power Systems, which at least yield a wide database of transit time data between four detectors, on which the method can be further tested (Bergdahl, 2002).

Conclusions
The results obtained by both simulations as well as to data from a system code and from a single application to a real case are promising, but further work is required in several areas. There is a thorough need for verification of the method, which in turn requires access to realistic void velocity profiles. Due to the lack of direct void velocity measurements, the closest possibility for validation is to make measurements with four LPRMs plus one movable TIP detector to obtain three transit times and the axial onset position of the boiling, and at the same time generate high-fidelity realistic void velocity profiles by system codes. These could be obtained either from dedicated measurements in critical assemblies or research reactors, or, more likely, from instrumented fuel assemblies at operating BWRs, such as all three Forsmark reactors, or Oskarshamn 3. Such validation work is planned in the continuation, as well as using the validated model for the next step, i.e. to determine the local void fraction from in-core noise measurements.