Gauss-Kronrod Integration For Rectangular Closed Conduits: P. Gruber, T. Staubli, F. Fahrni
Gauss-Kronrod Integration For Rectangular Closed Conduits: P. Gruber, T. Staubli, F. Fahrni
Gauss-Kronrod Integration For Rectangular Closed Conduits: P. Gruber, T. Staubli, F. Fahrni
Abstract
There exist cases of closed conduits applications, where the performance of the chosen configuration of the ATT flow
measurement is unsatisfactory. One method of improving the situation is to add more paths to the existing installations In
case of circular closed conduits and Gauss-Jacobi integration, there exists the nice feature, that by interlacing new paths in
between the old paths, the positions of the old paths do not change without losing in accuracy of the discharge calculation.
That means that polynomial deviations from the assumed weighting function W=(1-x2)1/2, of order 2(2N*+1)-1 =4N*+1
(equal the degree of exactness of the quadrature formula) can be integrated without error. N* corresponds to the number
of original paths. In case of an original 4-path in one plane configuration extended to a 9-path in one plane configuration,
N*=4, N=9, degree of exactness=17. This has an enormous advantage, because the already implemented path positions do
not have to be changed. In case of rectangular closed conduits, this is unfortunately no more the case. For the optimal
interlacing arrangement, all positions change. For rectangular conduits, in theory the positioning is less of a problem due
to the fact, that all the transducers have the same pill angle independent on the transducer positions. Nevertheless, from the
installation effort point of view, it is often a request to not move the old positions. Kronrod [1] derived optimal formulas
for the interlacing positions and all weights under the constraints of the old N* positions. Due to the constraints, the degree
of the quadrature formula is 3N*+1 (N* even, e.g. N*=4, degree=13) or 3N*+2 (N* odd). The paper gives the
corresponding positions and weights for the interesting cases N*=2, 3 and 4, and compares the results with the known
choice of positions and weights for assumed and simulated flow profiles. As the flow profiles are in reality not rational
functions, the loss in polynomial accuracy is not always clear.
1. Introduction
The performance of an ATT flow measurement installation is in applications, where the hydraulic situation is difficult or
not well known, hard to evaluate. In recent years, improved CFD flow simulations for hydraulic difficult situations allow
a reasonable estimate of the integration error of the measurement. The question of measurement location, orientation and
path configuration (number of layers, crossed, not crossed) are important parameters that can be varied for finding an
optimal solution by CFD simulations. Figure 1 shows two examples of crossed 8-path configurations (2 vertical planes and
4 layers each)
After CFD simulation studies or from experience of experts, the measurement is finally parametrized. A common example
is the 8-path (2 times 4 paths in 2 planes) configurations with Gauss-Jacobi positions and weightings for circular closed
conduits. After the installation, commissioning and putting into operation of such a measurement, the individual path
velocities respectively axial and transverse velocities, are recorded and monitored. If unusual profiles are observed and
installation errors can be excluded, one way to improve the accuracy of the flow measurement is increasing the number of
paths. Marushchenko, Gruber [2] showed that with an increase from 4 paths to 5 paths (centre path) in one plane, substantial
improvements can be achieved. The future revised version of IEC60041 standard [3] will open up the number of paths up
Table 1: general table for position & weights for variable number of paths in one plane
The table is based on the original Gaussian quadrature method with an additional weighting function W(z). In this
application of the method, the area flow function (AFF) F(z) corresponds to the weighting function W(z).
F(z) describes the distribution of the partial flow rates for each height z and is expressed by the product between the
averaged axial velocity and the width b(z) of the conduit at height z (Staubli & al [5])
At the height of the path position 𝑧𝑖 , the width 𝑏(𝑧𝑖 ) is equal to the projected path length Li,proj :
Integration of F(z) over the whole cross section yields the discharge Q
𝐷
𝑄 = ∫ 2𝐷 𝐹(𝑧)𝑑𝑧 (3)
−
2
The different integration weights of the Gauss-Jacobi and the OWICS method for circular section and Gauss-Legendre
and OWIRS methods have their origin in different assumption on the shape of the theoretical reference area flow func-
tion respectively velocity distribution. They can be distinguished by a single parameter . The assumed AFF is therefore
given by
z2 𝑚2
𝐹𝑟𝑒𝑓 (𝑧) = 𝐶 · (1 − 𝐷 2
) [ ] (4)
( ) 𝑠
2
𝑚2
𝐹𝑟𝑒𝑓 (𝑧) = 𝐶 · (1 − 𝜁 2 ) [ ] (5)
𝑠
With this assumed AFF Fref and the proper choice of , the positions and weights wi of Table 1 can be determined by the
Gaussian quadrature.
𝐷
𝑄= · ∑𝑁
𝑖=1 𝑤𝑖 · 𝑣̅𝑎𝑥 (𝑧𝑖 ) · 𝑏(𝑧𝑖 ) (6)
2
Gauss-Jacobi positions & Gauss-Jacobi positions & Gauss-Legendre positions & Gauss-Legendre positions &
Gauss-Jacobi weights OWICS weights Gauss-Legendre weights OWIRS weights
Table 2: augmented table for positions and weights proposed for IEC 60041 revision,
The original Gaussian N quadrature rule for the approximation of the integral of a function f with a weight function W is
given for an interval [-1, 1] by the following formula:
1
∫−1 𝑊(𝜁)𝑓 (𝜁)𝑑𝜁 = ∑𝑁 𝐺
𝑖=1 𝑤𝑖 𝑓(𝜁𝑖 ) + 𝑅𝑁 (𝑓) (7)
𝜁𝑖 and wi are the positions and the positive weights which can be determined such that the degree of exactness is [9]
𝑑𝑁𝐺 = 2𝑁 − 1 (8)
That means that a polynomial function f of order 2N-1 can be integrated by the weighted sum of equation (7) without
error:
Polynomial function of higher order or transcendental functions can only integrated with a remaining error term. With
increasing the number N, the approximation of equation (7) gets better for transcendental functions.
The position and weights wi and 𝜁𝑖 , i=1,…N are the ones listed in Table 1 for four different kind of weight function W,
where only the nonnegative positions are listed due to the fact that the positions are symmetrical to the 𝜁-axis.
The Gauss-Kronrod quadrature formula, extending equation (7), has the form [9]
1 ∗ ∗ 𝐾
∫−1 𝑊(𝜁)𝑓 (𝜁)𝑑𝜁 = ∑𝑁 ∗ 𝑁 +1 ∗
𝑖=1 𝑤𝑖 𝑓(𝜁𝑖 ) + ∑𝑙=1 𝑤𝑙 𝑓(𝜁𝑙 ) + 𝑅𝑁∗ (𝑓) (10)
Equation (10) means that the approximation of the integration is as follows: At the N* old positions 𝜁𝑖 , i=1,…N*, the N*
weights 𝑤𝑖∗ are newly computed. The old function values 𝑓(𝜁𝑖 ) are then added up weighted by 𝑤𝑖∗ . The second sum
considers the N*+1 new interlacing positions 𝜁𝑙 , l=1,…,N*+1. Kronrod derived the new positions 𝜁𝑙 , l=1,…,N*+1 and
the new weights 𝑤𝑖∗ , i=1,…,N* and 𝑤𝑙∗ , l=1,…,N*+1, such that despite the constraints of keeping the old positions fixed
the degree of exactness is
𝑑𝑁𝐾 = 3𝑁 ∗ + 1 for N* even
In case of an odd N* one degree can be gained because the centre position 0 is not changed. If we compare this exactness
with the one of the pure Gaussian quadrature formula where all N=2N*+1 positions and the corresponding weights are
determined, a reduction of exactness of
is inevitable.
For the interesting cases in for ATT applications, the following table results
which corresponds to the AFF for of equation (5). It turns out that for this case the Gauss-Kronrod formula is the
the N=2N*+1 Gauss-Jacobi formula for this weight. Therefore, in case of circular closed conduits and Gauss-Jacobi
integration, there exists the nice feature, that by interlacing N*+1 new paths in between the old N* paths to total of
N=2N*+1 paths (e.g. N*=4, 2N*+1=9), the positions of the old paths do not change without losing in accuracy [10]. That
means, that polynomial deviations from the assumed weighting function 𝑊 (0.5,0.5) (𝜁) = (1 − 𝜁 2 )0.5 of order 2(2N*+1)-1
=4N*+1 (equal the degree of the quadrature formula, e.g. 17) can still be integrated without error. This has the additional
benefit, that the already implemented path positions do not have to be changed. This can for instance be seen in Table 1 if
one goes up in the case of Gauss-Jacobi from N=2 to N=5, the position 0.5 remains the same.
The positions in bold figures indicate the Gauss-Legendre positions for the old path numbers N*=2, 3 and 4
N=5
Gauss-Legendre positions & weights OWIRS positions & weights Gauss- Legendre positions &
OWIRS weights
di/(D/2) wi di/(D/2) wi di/(D/2) wi
0 0.568889 0 0.554092 0 0.562705
+/-0.538469 0.478629 +/-0.525989 0.470657 +/-0.538469 0.485402
+/-0.906180 0.236927 +/-0.893646 0.245772 +/-0.906180 0.228094
Fig. 3: path positions and weights for G-L (Gauss-Legendre), G-K (Gauss-Kronrod), OWIRS and Legendre positions
and OWIRS weights
4. Analytical examples
Table 6: Discharge measurement of the four cases for the five polynoms
- Gauss-Legendre: best performance. Up to 9th order polynomials, the error is zero. For higher order a slow increase
of the error can be observed.
- Gauss-Kronrod: second best performance. Up to 7 th order polynomials, the error is zero. For higher order, the
error rises up to 2-3 times larger values as for the Gauss-Legendre case.
- G-L positions & OWIRS weights: Even for small order polynomials, the error is as high as for the G-L case for
high order polynomials. With increasing order of the polynomials, the error is comparable to the Gauss-Kronrod
case.
- OWIRS: The worst performance because the method is optimized for the OWIRS profile. For higher order
polynomials, the error is ~constant (0.15%).
The four parameters are chosen, such that a realistic flow profile can be obtained, which is similar to the Gersten Herwig
profile [11,12]. The maximal value and the value for the representative wall distance were defined accordingly. The inner
exponential decay rate T2, ymax and were kept constant, while the wall exponential decay rate was varied between
1.4
1.2
0.8
0.6
0.4
0.2
0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
Again, the test is carried out with the approximation order of N=5.
Table 8: Integral (discharge) of the four cases for five exponential profiles
Table 9: measurement error of the four cases for five exponential profiles, green: best, red: worst
- Contrary to the polynomial profiles, this choice of exponential profiles never has Gauss-Legendre as winner. In
two cases, Gauss-Legendre even shows the worst performance.
- Surprisingly enough, the Gauss-Kronrod approximation is in three out of 5 cases the best. The reason for these
results is not yet clear. However, it seems to be that the interlacing positions in between the old Gauss-Legendre
positions can better cope with the exponential behaviour of the profile.
- OWIRS behaves badly in 3 out of 5 cases. This is not surprising, as the OWIRS positions and weights are not
optimized for an exponential profile.
5. CFD simulations
CFD simulations were carried out for two types of cross sections and three different flow velocities. The two cross
sections are:
- Square shape with side length a=b=1.4m
- Rectangular shape with 1.5m width a=1.5m and height b=4.5m.
The velocities v are chosen as 0.2m/s, 2m/s and 10m/s. This results for the Reynolds number
2𝑎𝑏
𝑅𝑒 = 𝑣 106 (18)
𝑎+𝑏
The simulated square cross section is shown in Figure 6. For the rectangular cross section the length of the conduit is also
6m.
Fig. 6: simulated square cross section, 1.4m x 1.4m x 6m, 2E10P (two planes with five paths in each plane)
Figure 7 shows contour plot of the axial velocity (x-axis) for the square cross section. It is clearly visible that at the centre
of the square the maximal velocity is attained. Figure 8 visualizes the secondary flow components. For each simulated
Fig. 7: axial velocity contour for square and 2m/s Fig. 8: secondary flow for square and 2m/s
The simulations were performed with ANSYS CFX 18. The mesh of the conduit is a manually generated structured
hexahedral mesh. The inlet boundary comprises about 22500 hexahedral elements. Minimum and maximum angles were
90° and the maximum volume ratio were 1.2. The dimensionless wall distance (y+) is on average about 0.01 (for 0.2 m/s),
4 (for 2m/s) and 20 (for 10 m/s). The simulations were performed in steady state and with translational periodic boundary
conditions. The solutions are based on the SST (shear stress transport) turbulence model. The calculations are solved with
the high resolution advection scheme and automatic timescale. The RMS residuals were lower than 10-12, the maximum
residuals were at 1.1·10-11, and the imbalance at 4.8·10-12. The simulations required about 300 iterations to reach this
convergence. The wall roughness was chosen as the one of sand/grain and set equal to 0.1mm For steel the corresponding
value would be ten times smaller.
In Figure 9 the 3-dimensional plots of the normalized velocity profiles for the three velocities are shown. From these it
seems that the shape of the profiles do not change much, especially for larger velocities. This can also be seen in the three
area flow functions of Figure 10.
Fig. 9: normalized velocity profiles: left 0.2m/s, middle 2m/s, right 10m/s
Fig. 10: Area flow functions (AFF): left 0.2m/s, middle 2m/s, right 10m/s
Figure 11 summarizes all the obtained relative error results. The following remarks can be added:
c) d)
e) f)
- With increasing number of paths all error curves decrease strongly and slow down for higher number of paths.
- The Legendre positions & OWIRS case and the pure OWIRS case behave very similar and are clearly superior
to the other two cases. Only in the case a): velocity of 0.2m/s and square cross section, the error is with ~0.1% for
all path configurations constant and relatively high. For low velocities and square geometry, the velocity profile
seems less to be of the OWIRS type as for higher velocities. For this situation it could be interesting to optimize
the parameter as was suggested by Gruber & al [13].
- It is interesting to note that the Gauss-Kronrod case always performs better than the Gauss-Legendre case. To be
specific, for N= 5 paths the improvement is ~0.1%,<~0.01% for N=7 and finally ~0,02% for N=9. The reason for
the obtained type of decrease of improvement is not clear at the moment.
6. Conclusions
In situations, where the accuracy of the ATT multipath discharge measurement in rectangular conduits seems to be not
reliable enough due to unknown or not removable adverse effects, the method of increasing the path number is a viable
option. The strategy should be as follows:
- If all path positons can easily been shifted then an incremental increase from for instance 4 to 5 path with OWIRS
positions and weights is for realistic cases the best choice.
- If all path positons can easily been shifted and in cases of a purely constant weighting function and multiplicative
polynomial deviations the Gauss-Legendre is clearly the best method.
References
[1] A. S. Kronrod: Nodes and weights of quadrature formulas. Sixteen-place tables, New York 1965: Consultants
Bureau (Authorized translation from the Russian)
[2] S. Marushchenko, P. Gruber: “Comparative Study of 4(8)-path and 5(10)-path configurations for ATT flow
measurements in circular conduits”, IGHEM 2014, Itajuba, Brasil
[3] IEC60041 norm: Field acceptance tests to determine the hydraulic performance of hydraulic turbines, storage pumps,
storage pumps and pump-turbines, 1991
[4] ASME PTC-18 norm: Hydraulic turbines & pumps, 2016
[5] T. Staubli, P. Gruber, F. Fahrni: “Diagnosis of Acoustic Transit Time Data based on the Area Flow Function”,
IGHEM 2018, Beijing
[6] T. Tresch, B. Lüscher, T. Staubli, P. Gruber: “Presentation of optimized integration methods and weighting correction
for the acoustic discharge measurement”, IGHEM 2008, Milano
[7] B. Lüscher, T. Staubli, T. Tresch, P. Gruber: “Optimizing the acoustic discharge measurement for rectangular conduits”,
IGHEM 2008, Milano
[8] https://en.wikipedia.org/wiki/Gauss–Kronrod_quadrature_formula, cited June 2018
[9] S. E. Notaris: “Gauss-Kronrod Quadrature Formulae – Survey of Fifty Years of Research, Electronic Transactions on
Numerical Analysis, vol. 45, pp. 371-404, 2016
[10] I. P. Mysovskih: “A special case of quadrature formulae containing preassigned nodes”, (Russian), Vesci Acad. Navuk
BSSR Ser. Fiz.-Techn. Navuk, 4 (1964), pp.125-127
[11] Arbeitsgruppe für laseroptische Strömungsdiagnostik, 2009
[12] Gersten, K., Herwig, G.: Strömungsmechanik. Grundlagen der Impuls-, Wärme- und Stoffübertragung aus asymp-
totischer Sicht. 1. Aufl., Vieweg-Verlag, Braunschweig Wiesbaden 1992.
[13] P. Gruber, T. Staubli, T. Tresch, F. Wermelinger: Optimization of the ADM by Adaptive Weighting for the
Gaussian Quadrature Integration, IGHEM 2010, Roorkee, India