Sensitivity analysis using contribution to sample variance plot: Application to a water hammer model

Reliability Engineering and System Safety 99 (2012) 62–73

Contents lists available at SciVerse ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Sensitivity analysis using contribution to sample variance plot: Application to a water hammer model S. Tarantola a, V. Kopustinskas b,n, R. Bolado-Lavin b, A. Kaliatka c, E. Uˇspuras c, M. Vaiˇsnoras c a

European Commission, Institute for Protection and Security of Citizen, Ispra, Italy European Commission, Institute for Energy and Transport, Petten, Netherlands c Lithuanian Energy Institute, Kaunas, Lithuania b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 29 March 2011 Received in revised form 13 October 2011 Accepted 16 October 2011 Available online 25 October 2011

This paper presents ‘‘contribution to sample variance plot’’, a natural extension of the ‘‘contribution to the sample mean plot’’, which is a graphical tool for global sensitivity analysis originally proposed by Sinclair. These graphical tools have a great potential to display graphically sensitivity information given a generic input sample and its related model realizations. The contribution to the sample variance can be obtained at no extra computational cost, i.e. from the same points used for deriving the contribution to the sample mean and/or scatter-plots. The proposed approach effectively instructs the analyst on how to achieve a targeted reduction of the variance, by operating on the extremes of the input parameters’ ranges. The approach is tested against a known benchmark for sensitivity studies, the Ishigami test function, and a numerical model simulating the behaviour of a water hammer effect in a piping system. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Sensitivity analysis Contribution to sample variance plot Variance reduction Water hammer RELAP model Ishigami function

1. Introduction Let us consider a simulation model represented by an input/ output function Y ¼G(X), where Y is a scalar model output and X ¼(X1,X2,y,Xn) defines the generic vector of n input parameters. The values of the model parameters are not perfectly known; in other words, they are affected by some uncertainty. Therefore, although it is not the case, each input can be considered as a random variable characterized statistically by its probability density function pi(Xi). The model output Y can also be thought as a random variable; the estimation of its pdf is the objective of uncertainty analysis and is carried out by evaluating the function G(X) on a sample of points generated from the pdf pi(Xi). The analysis described here is in the context of deterministic models. Probabilities have only been introduced to represent imprecision about model parameter values. Graphical sensitivity tools provide valuable information on the relationship between uncertain model inputs and model outputs. In 1989, Sacks et al. proposed the use of scatter plots, i.e. projections on a two-dimensional plane of the hyper-surface describing the input/output mapping [1]. A number of other graphical techniques are discussed in [2–5]. In 1993, Sinclair [6] introduced the contribution to the sample mean plot (CSM) which was further developed by Bolado-Lavin et al. [7].

n

Corresponding author. Tel.: þ39 0332 786257. E-mail address: [email protected] (V. Kopustinskas).

0951-8320/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2011.10.007

The idea behind CSM is to use a given random sample of the input parameters—that is generally used for uncertainty analysis, to draw conclusions about the sensitivity of the model output. This paper presents an extension of CSM plots, called ‘‘contribution to sample variance plots’’ (CSV). When investigating input–output relationship, CSV contains a considerable higher amount of information than that provided by standard global sensitivity indices [8]. Once the most important input has been detected, global sensitivity indices do not inform the analyst about how to act operatively in order to reduce the range of uncertainty of the important input for a given target reduction of the output variance. Contrarily, CSV does give us the amount of the variance reduction that would be achieved for any arbitrarily chosen restriction of the input uncertainty range. The CSV measure is defined and its properties are described in Section 2. The interpretation of CSV is given in terms of change of variance that can be achieved by trimming the range of the input parameters in Section 3. The CSV plot can use the same sample points utilized for the CSM, for the scatter-plots, and for uncertainty analysis in general. Therefore, new information can be obtained without additional model evaluations. As the CSM is a powerful tool to identify local regions of the input space that contribute substantially to the mean value of the model output, the CSV is powerful to localize areas where the contribution to the variance of the model output is considerable. The theoretical results are shown on the Ishigami test function in Section 4. The CSV technique is also applied to a numerical model simulating the behaviour of a water hammer in a pipeline.

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

Table 1 Notation. n N Xi Y V(  ) E(  ) p(  ) F(  ) G(  ) CSMXi(  ) CSVXi(  ) q

Number of parameters Sample size Model parameter i Model output value Variance Mean value Probability density function Cumulative distribution function Model function Contribution to sample mean for parameter Xi Contribution to sample variance for parameter Xi quantile

The model is described in Section 5 and the results of CSV and the comparison with the results from CSM are shown in Section 6. Table 1 presents notation used throughout the article.

63

computational cost, i.e. using the simulations already carried out for the scatter-plot (and the CSM). Given a set of N sample points (xj1,xj2,y,xjn), j¼1,2,y,N and the corresponding model output values yj, the CSV for input Xi can be computed using the following procedure: 1. Compute the output sample mean Ym and transform each observation (y1,y2,y,yN) by subtracting the mean value Ym: ytj ¼yj  Ym, j ¼1,2,y,N. The transformed output Yt has zero mean value. (2) (N) 2. Sort in ascending order the sample of input Xi (x(1) i ,xi , y, xi ) and the corresponding set of the transformed outputs ytj, obtaining a series (yt(i,1), yt(i,2), y, yt(i,N)). 3. The CSV at q quantile for the input parameter Xi is computed as PbqNc j¼1

yt2ði,jÞ

j¼1

yt2ði,jÞ

CSV X i ðqÞ ¼ PN

,

q A ½0,1

ð3Þ

where qN is the largest integer not greater than qN. The CSV for input Xi is obtained by plotting the CSVXi(q) against the cumulative distribution function of Xi. The plot consists of pairs (Fi(xi), CSVXi(Fi(xi)) for each sample point xi of Xi.

2. The contribution to the sample variance plot Let us recall the definition of CSM for a given input Xi [7]: 1 EðYÞ

CSMX i ðqÞ ¼

EðYÞ ¼

Z

Z

1

1

Z

F 1 i ðqÞ

... 1

1

1

n Y

pi ðxi ÞGðx1 ,x2 ,. . .,xn Þ

dxi dx1 . . .dxi1 dxi þ 1 . . .dxn Z 1 Y n ... pi ðxi ÞGðx1 ,. . .,xn Þ dx1 . . .dxn

Z

3. Interpretation of the CSV plot

i¼1

1

1

ð1Þ

1 i ¼ 1

where qA[0,1], E(Y) is the mean value of the model output, Fi 1(q) is the inverse cumulative distribution of Xi at quantile q. The multiple integral in (1) is computed on the range [  N,N] for all input parameters except for Xi, for which the range is [  N,Fi 1(q)]. CSMXi(q) is plotted on the [0,1]2 space: q is a point on x-axis representing a fraction of distribution range Xi, and CSMXi (q) is a fraction of the output mean corresponding to the values of Xi smaller or equal than its q-quantile. By definition, CSMXi (0) ¼0 and CSMXi (1)¼1. The CSV for the variable Xi (CSVXi) is defined similarly to CSM: CSV X i ðqÞ ¼

VðYÞ ¼

Z

1 VðYÞ

Z

Z

1

1

Z

F 1 i ðqÞ

... 1

1

1

n Y

pi ðxi ÞðGðx1 ,x2 ,. . .,xn Þ

i¼1

EðYÞÞ2 dxi dx1 . . .dxi1 dxi þ 1 . . .dxn Z 1 Y n 1 ... pi ðxi ÞðGðx1 ,. . .,xn ÞEðYÞÞ2 dx1 . . .dxn

1

ð2Þ

Both CSM and CSV are tools to estimate the contribution to the sample mean or sample variance of a particular range of input parameter values. If CSM or CSV are close to diagonal, it indicates that the contribution to the mean or to the variance is equal throughout the range of the input parameter. The CSV value at point q¼0.1 provides an estimate of the model output variance due to the 10% of the smallest values of the input parameter. Where the CSV curve is steep, the contribution to the sample variance is large (variance is larger than on average). Where CSV curve is flat, contribution to the sample variance is small (variance is smaller than on average). Similar considerations are valid also for the CSM plot when referring to the contributions to the sample mean. The CSV plot is a useful tool to analyse the effect of reduced range of an input parameter to the variance of the model output. Further, we will develop a relationship between input parameter range reduction and change of the model output variance. Let us consider a standard setting with input parameter Xi having probability density function pi(xi) over [ N,N]. In the new setting, the range of Xi is reduced to [u,z],  Nou ozoN, having probability density function pni (xi). It can be shown that

1 i ¼ 1

where qA[0,1], V(Y)—variance of the model output, Fi 1(q) is the inverse cumulative distribution of Xi at quantile q. The multiple integral in (2) is computed on the range [ N,N] for all input parameters except for Xi, for which the range is [ N,Fi 1(q)]. CSVXi is also plotted on the [0,1]2 space, where CSVXi(q) is a fraction of the output variance corresponding to the values of Xi smaller or equal than its q quantile. Also for the CSV holds CSVXi(0) ¼0 and CSVXi(1)¼ 1. It is important to note that the CSV is defined in (2) as a contribution to variance with respect to constant mean E(Y) over the full range of all parameters. The CSM allows the analyst to identify the effect of each model input on the average of the model output for given percentiles within its uncertainty range; the CSV does the same for the variance. Graphical tools provide useful insights into the I/O relationship using reasonably low number of model simulations. CSV plots can be placed side by side to scatter-plots and to CSM plots at no extra

pni ðxi Þ ¼ R z u

pi ðxi Þ pi ðvÞ dv

ð4Þ

The variance V(Yn[u,z]) in the new setting (range of the input parameter Xi is [u,z]) is computed as Z 1 Z 1Z z ... p1 ðx1 Þp2 ðx2 Þ. . .pni ðxi Þ. . .pn ðxn ÞðGðx1 ,. . .,xn Þ VðY n½u,z Þ ¼ 1

1

u

EðYÞÞ2 dxi dx1 . . .dxi1 dxi þ 1 . . .dxn Z 1 Z 1Z z 1 ¼ Rz ... p1 ðx1 Þp2 ðx2 Þ. . .pi ðxi Þ. . .pn ðxn ÞðGðx1 ,. . .,xn Þ 1 u u pi ðvÞ dv 1 EðYÞÞ2 dxi dx1 . . .dxi1 dxi þ 1 . . .dxn

ð5Þ

The variance (5) defines variance of the model output with the reduced range of the parameter Xi, but with respect to constant mean E(Y) over the full range of all parameters. Further

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

combining (5) with (2) and noting that Fi 1(q)¼z and u ¼  N: Z 1 Z 1 Z F 1 n i ðqÞ Y 1 ... pi ðxi ÞðGðx1 ,x2 ,. . .,xn Þ CSV X i ðqÞ ¼ VðYÞ 1 1 1 i¼1

GðX 1 ,X 2 ,X 3 Þ ¼ sin X 1 þ 7sin2 X 2 þ 0:1X 43 sin X 1

ð6Þ

where V(Yn[  N,z]) is the variance of the model output over the reduced range [ N,z] of parameter Xi. Eq. (6) provides direct interpretation of the CSV plot: CSV value for the input parameter Xi at each quantile q, qA[0,1], is proportional to the ratio between the variance over the reduced range [ N,Fi 1(q)] of Xi and the total variance V(Y). The proportion Rz constant 1 pi ðvÞ dv is derived from (4) and reflects changes in the distribution function due to the reduction of the range. Another important result can be obtained by analysing the difference between two points in the CSV plot: CSVXi(q2)– CSVXi(q1). In this case Eq. (6) transforms as follows by noting Fi 1(q1)¼u and Fi 1(q2) ¼z: CSV X i ðq2ÞCSV X i ðq1Þ Z 1 Z 1 Z F 1 n i ðq2Þ Y 1 ... pi ðxi ÞðGðx1 ,x2 ,. . .,xn Þ ¼ 1 VðYÞ 1 1 F i ðq1Þ i ¼ 1 EðYÞÞ2 dxi dx1 . . .dxi1 dxi þ 1 . . .dxn Z z VðY n½u,z Þ ¼ pi ðvÞ dv VðYÞ u

ð7Þ

Eq. (7) transformed into the following form: CSV X i ðq2ÞCSV X i ðq1Þ VðY n½u,z Þ ¼ Rz VðYÞ u pi ðvÞ dv

ð8Þ

provides important practical information: by choosing q1 and q2, reduction of the variance due to reduced range of the input parameter Xi from [0,1] to [q1,q2], 0 oq1oq2 o1 can be estimated. As q1 and q2 can be any from the interval [0,1] satisfying the condition q1oq2, a 3D plot can be created (shown in the next section) indicating the region where the variance reduction can be most effectively achieved. The classical sensitivity analysis does not tell the analyst how to reduce the range of an important input parameter, but CSV does. Eq. (8) holds under the condition that the variance in (5) is defined as a variance over the reduced range of the parameter Xi, but with respect to constant mean E(Y) over the full range of all parameters. This is important to note because if in the reduced range of Xi the mean value changes significantly, the variance as defined in (5) increases as well. Following similar reasoning as for derivation of Eqs. (6)–(8), the following relationship is also valid for the CSM plot: CSMX i ðq2ÞCSM X i ðq1Þ Z 1 Z 1 Z F 1 i ðq2Þ 1 ... Gðx1 ,x2 ,. . .,xn Þ ¼ 1 EðYÞ 1 1 F i ðq1Þ

Z ¼

u

20

pi ðvÞ dv

EðY n½u,z Þ EðYÞ

10 0 -10 -4

-3

-2

-1

0

1

2

3

4

1

2

3

4

1

2

3

4

Parameter X1 20 10 0 -10 -4

-3

-2

-1

0 Parameter X2

pi ðxi Þ dxi dx1 . . .dxi1 dxi þ 1 . . .dxn

i¼1 z

ð10Þ

and input parameters follow uniform distribution U(  p, p). The scatter plot for (10) is shown in Fig. 1 by using sample size of 1000 points. The analytical values of the sensitivity indices can be computed for (10). The first-order sensitivity indices Si ¼V(E(Y9Xi))/ V(Y) are as follows: S1 ¼0.3138; S2 ¼0.4424; S3 ¼0. This indicates that X2 is the most important parameter, while X3 is not important to the mean of the model output in terms of the first order effect. The higher order indices and total indices can be computed analytically and are presented in [10]. As shown in [7] similar results can be obtained by using CSM plot technique with simple random sample. Fig. 2 shows the CSM plot for Ishigami function (10) by using analytic computations. By definition the CSM indicates importance of the input parameter to the mean value of the model output, so it can be directly compared to the first-order indices Si. Note that the CSM curve for X3 is exactly diagonal. The more the CSM curve deviates from the diagonal, the more important is the parameter. The CSM curve for X1 has a maximum distance from the diagonal. In general, the CSM curve can cross the diagonal several times and in this case [7] suggests using the sum of maximum distances from the diagonal (in absolute values). According to the criterion, the most important factor is found to be X2, as it has the largest sum of the maximum distances of the CSM curve to the diagonal (Table 2). This result is consistent with the ranking provided by the first order indices. The CSV plot for Ishigami function (10) is shown by using analytic computations (Fig. 3) and random sample of size 1000 (Fig. 4). The plot is different from the CSM plot. The CSV indicates

20

ð9Þ

The interpretation of the result in (9) is that by reducing the range for the input parameter Xi from [0,1] to [q1,q2], the reduction in the model output mean can be estimated. This result might be important if the model output represents a critical value for the integrity of the system (pressure, temperature, etc.), and its reduction has important practical implications.

G(X1,X2,X3)

n Y



The Ishigami function [9] is a frequently used test function in sensitivity analysis, having three input parameters:

G(X1,X2,X3)

EðYÞÞ2 dxi dx1 . . .dxi1 dxi þ 1 . . .dxn Z z VðY n½1,z Þ ¼ pi ðvÞ dv VðYÞ 1 CSV X i ðqÞ n½1,z VðY VðYÞ Þ ¼ Rz 1 pi ðvÞ dv

4. CSM/CSV results on Ishigami test function

G(X1,X2,X3)

64

10 0 -10 -4

-3

-2

-1

0 Parameter X3

Fig. 1. Scatter plot of Ishigami function (10), N ¼1000.

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

Ishigami function

1

0.8

0.7

0.7

0.6

0.6

0.5 0.4

0.5 0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

X1 X2 X3 Diagonal

0.9

CSV(q)

CSM(q)

0.8

Ishigami function

1

X1 X2 X3 Diagonal

0.9

65

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Cumulative fraction of input parameter range, q

0

1

Fig. 2. The CSM plot for Ishigami function (10) computed analytically.

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Cumulative fraction of input parameter range, q

1

Fig. 4. The CSV plot for Ishigami function (10), random sample N ¼ 1000.

Table 2 The maximum distance of the CSM curve from the diagonal and their sum for Ishigami function.

Sum

CSMX1(q1)

q2

CSMX2(q2)

q3

CSMX3(q3)

0.5

0.07

0.125 0.375 0.625 0.875

0.02

0

0

0.07

0.08

Ratio of variances: Fx1(q1,q2)

Values

q1

0

Ishigami function

1

2 1.5 1 0.5

0.8

0.6

0.2

X1 X2 X3 Diagonal

0.9

1 0.8

0 0

0.6 q1 v aria ble

0.2

0.8 1

0.7

ble

0.4

0.4

ria

a 2v

q

0

Fig. 5. The variance ratio function Fx1(q1,q2) for the input parameter X1.

CSV(q)

0.6 In this test case when all parameters are distributed uniformly  U(  p, p), Eq. (4) is reduced to

0.5 0.4

pi ðxi Þ ¼

0.3

ð11Þ

and correspondingly Eq. (7) transforms into

0.2

CSV X i ðq2ÞCSV X i ðq1Þ ¼ ðq2q1Þ

0.1 0

zu n p ðx Þ 2p i i

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Cumulative fraction of input parameter range, q

1

Fig. 3. The CSV plot for Ishigami function (10) computed analytically.

that input parameter X3 deviates from the diagonal most severely. This result can be explained as the output variance is the most unevenly distributed within the range of X3. By looking at the scatter plot (Fig. 1) it is evident that output variance is large at both ends of X3 range and small in the middle of the range, while the mean is constant throughout the range. This information can be only qualitatively assessed from the scatter plot; however, the CSV plot allows to quantify the effect.

VðY n½pð2q11Þ, pð2q21Þ Þ VðYÞ

ð12Þ

where Fi 1(q)¼ p(2q  1). Let us define the variance ratio function Fxi(q1,q2) as follows: CSV X i ðq2ÞCSV X i ðq1Þ VðY n½pð2q11Þ, pð2q21Þ Þ , ¼ ðq2q1Þ VðYÞ 0 rq1r q2 r1:

Fxi ðq1,q2Þ ¼

ð13Þ

The function Fxi(q1,q2) represents the ratio of two variances—variance of the model output with reduced range [q1,q2] with respect to the mean over the full range for input parameter Xi as defined in (5) and the total variance of the model calculated over the full range [0,1]. The variance ratio function for the input parameters X1–X3 of (10) are shown in Figs. 5–7. The computations were performed analytically by using Mathematicas. As evident from those figures, the most effective variance

66

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

1.5 1 0.5 0 0

1 0.8

0.2 q1

0.6

0.4 va

ria

0.4

0.6 ble

0.2

0.8

le

ab

ari

v q2

Ratio of variances (reduced/full)

Ratio of variances: Fx2(q1,q2)

1.6

0

1

X1 X2 X3

1.4

1.2

1

0.8

0.6

0.4

0

Fig. 6. The variance ratio function Fx2(q1,q2) for the input parameter X2.

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 q

Fig. 8. The effect on variance by symmetric reduction of the range [q,1  q] for all input parameters.

X1 X2 X3

4.5 4 3 2 1 1

0 0

0.8 0.2

0.6 0.4 q1 va

riab

0.4

0.6 le

0.2

0.8 1

q2

ble

ria

va

0

Ratio of variances (fixed/full)

Ratio of variances: Fx3(q1,q2)

5

4 3.5 3 2.5 2 1.5 1 0.5

Fig. 7. The variance ratio function Fx3(q1,q2) for the input parameter X3.

0

reduction can be achieved by reducing the range of X3. The 3D variance ratio plot enables to derive a number of interesting observations: – The effect on variance change can be estimated for any range of q1 and q2 satisfying condition q1rq2. – The counter-diagonal line (from (0,1) to (0,5,0,5)) represents symmetric reduction of the range from both sides of the range. This can be plotted on a separate two-dimensional plot as well and illustrated in Fig. 8. – The diagonal line (from (0,0) to (1,1)) represents the situation q1¼q2 ¼q or in derivative of CSV: Fxi(q,q)¼ (dCSVXi(q)/dq) and means the variance change by fixing parameter Xi at quantile q. This can be plotted on a separate two-dimensional plot as shown in Fig. 9. Similarly, the ratio of means can be investigated from Eq. (9). This will show the effect of the reduction in the range of input parameters on the model output mean. Fig. 8 illustrates the power of CSV plot to quantify the reduction of variance by reducing symmetrically ranges of input parameters. This plot indicates that by reducing the range of X3 from [  p, p] to [  0.8p, 0.8p] which corresponds to 10% reduction from both sides of cumulative function, the variance is reduced by 40% (variance ratio of 0.6). The plot also indicates that by further reducing the range of X3 the variance can be

0

0.1

0.2

0.3

0.4

0.5 q

0.6

0.7

0.8

0.9

1

Fig. 9. The effect on the variance by fixing parameters at q for all input parameters.

reduced up to maximum 52%. The plot also indicates that by fixing X1 at q ¼0.5 even lower value of variance can be achieved. Fig. 8 indicates that the variance increases when the range of X2 is symmetrically reduced towards the zero point (q¼0.5). This is due to the fact that the variance in (5) is defined with respect to the mean over the full range. In this case, as clear from the scatter plot (Fig. 1) the mean around zero decreases compared to the global mean over the full range of X2 and this causes the variance as defined in (5) to exceed the total variance. Fig. 9 shows the effect on the variance by fixing each parameter at different values of q, qA[0,1]. The mathematical interpretation of it is derivative of the CSV. In case of parameter X3, fixing it at the point –p (q¼ 0) or p (q¼1), the variance increases by nearly five times, but by fixing it at the point 0 (q¼ 0.5) the variance is reduced more than twice (52%). Similar interpretation for the other parameters can be derived from Fig. 9.

5. UMSICHT test facility experimental case For the present case study, water hammer experimental test performed at Fraunhofer Institute for Environmental, Safety and

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

closing valve in a steady-state liquid flow. A centrifugal pump produces steady-state flow into the circuit from the pressurized vessel (B2) into the test pipe section back to the vessel. When at t¼ 0 s the valve closes rapidly while the pump is still running, pressure waves are induced in the whole pipe system and measured by fast pressure transducers (P01–P23). During the first phase of the transient, a rarefaction wave is travelling inside the pipe towards the downstream reservoir. As a consequence, cavitation occurs downstream the valve and a vapour bubble is formed. The generated pressure wave oscillates in the pipeline between the vessel and the vapour bubble until the condensation, inducing a cavitational hammer. The objective of the Pilot plant pipework test is to investigate pressure surges in fast closing processes and to provide experimental data for the validation of calculation software for water supply and chemical plants [12].

Energy Technology (UMSICHT) in Oberhausen was studied [11]. The test facility was modelled by using RELAP5/Mod3.3 thermal hydraulic code and the modelling results were compared with the experimental data. The UMSICHT test facility enables various operations and transport of compressible and incompressible liquid due to a modular construction system. Using modern high speed measurement (frequency 1–10 kHz) the local phase distribution, the system pressure, the fluid velocity as well as the effective force on the pipe restraints can be measured and calculated. The configuration of the Fraunhofer UMSICHT test loop Pilot Plant Pipework (PPP) is shown in Fig. 10. The piping length measurements start from the closure valve which initiates water hammer transient (located between P02 and P03). The piping lengths downstream the closure valve are marked positively. The piping lengths upstream the closure valve are marked negatively. The circulation loop is closed at the pressurized vessel (point B2). The total test pipe length (from FP1 to FP3) is 137 m in order to simulate typical power plant piping. Within the experimental test, the fast acting valve was fitted in the horizontal part of the 110 mm inner diameter piping system. The test was carried out with water at room temperature (20 1C) and initial fluid velocity v0 ¼4 m/s. Fig. 10 shows the UMSICHT test loop length marks. The experiments were conducted using the dynamic behaviour of

5.1. UMSICHT facility modelling with RELAP5 code The UMSICHT test facility model was developed by RELAP5/ Mod3.3 code and its nodalization scheme is shown in Fig. 11. Two time-dependent volumes (components ‘‘500’’ and ‘‘650’’) with the specified constant pressures and temperatures to obtain steady

P09

44.4 m

-8.7 m

P23

P01

P03

P06

WM 0.2 m

-0.2 m

P15

67.0 m

60.8 m

34.5 m

81.6 m

145.5 m 139.4 m

84.6 m

FP 2

P18 88.7 m

B2

50.9 m

0m

P02 FP 1

146.8 m

P12

bridge

closure valve

-14.5 m

67

77.5 m

142.9 m

67.9 m

FP 3 -18.2 m

149.4 m

137.0 m

90.7 m

75.5 m

Fig. 10. Configuration of the Fraunhofer UMSICHT test loop PPP with pressure measurement points: B2—the pressurized vessel, FP1–3—forces on pipe supports measurement points, P01–P23—pressure transducers, WM—wire mesh sensor.

44.4 m

Flow direction VALVE 754

PIPE 753

-18.2 m

p1

SINGLE JUNCTION 752

PIPE 755

WM

0.0 m

650 TMDPVOL 2

67.0 m

0.2 m

500 TMDPVOL 1 p2

50.9 m

BRIDGE

P03

84.6 m

SINGLE JUNCTION 759

139.4 m

137.0 m

PIPE 758

NEW PIPE BRIDGE 2 81.6 m

90.7 m

88.7 m

67.9 m

SINGLE JUNCTION 757

148.0 m 75.5 m

Fig. 11. Pilot plant pipework (UMSICHT) RELAP5/Mod3.3 code model nodalization scheme.

68

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

state liquid velocity were simulated in the model. This approach was used to avoid modelling of the pump, which operates in the actual facility. The facility piping from the tank upstream the shut-off valve (component ‘‘754’’) was simulated using the pipe component ‘‘753’’. The segment of actual facility piping with the bridge downstream the valve was modelled using the pipe component ‘‘755’’. The last segment of the piping with the pipe bridge 2 was modelled by employing the pipe component ‘‘758’’. Adapting the developed model for water hammer analysis, high attention must be given to selection of the cell (control volume) size, valve model, time step of calculation and other parameters. Adoption process of UMSICHT test facility model is presented in a number of articles [14,15]. RELAP5 modelling of water hammer experiment in UMSICHT test facility has shown that the first (and the highest) calculated

5

Relap5 calculated UMSICHT measured

Pressure, MPa

4 3

6. Selection and quantification of the most important input parameters The main task of the RELAP5 calculations was validation of the software code, therefore evaluation of sensitivity and uncertainties is very important. The sensitivity analysis started with the selection of parameters influencing the calculation (pressure behaviour downstream the closed valve) results. There are two types of such parameters: (1) parameters describing initial condition of the system, and (2) model parameters. The following two parameters specify the initial condition of the system: X1 X2

2

water pressure in the pump header (pos. P01 in Fig. 10); water temperature in the system.

The following model parameters may contribute to the model output uncertainties:

1 0 0

2

4

6

8

10

Time, s

X3 X4 X5–X7

valve closing rate; pipe wall roughness; flow energy loss coefficients in different piping segments.

0.4

Table 3 provides the list of the selected input parameters and properties of their distributions for the sensitivity analysis. The uncertainty range of the input parameters was estimated from the UMSICHT measurement data, technical data sheets and handbooks. The random sample for sensitivity analysis was generated by using truncated normal distribution, accounting for the range of double standard deviation to both directions from the mean, corresponding to 2.5% truncation from both sides.

0.2

7. Results of CSM and CSV plots on the water hammer model

1

Relap5 calculated epswms (UMSICHT measured)

0.8 Void fraction, -

pressure peak matches very well the measured value of pressure. Comparison of the basic case calculation with experimental data is presented in Fig. 12. The prediction of the first cavitation hammer is the most important, because the attained value during the first pressure peak is the most dangerous in comparison with the successive pressure peaks and can lead to damages of the equipment (valves, pumps, pipe bends) in the piping. Therefore the further analysis was carried out only for investigation of the first pressure increase peak.

0.6

0 0

2

4

8

6

10

Time, s Fig. 12. Pressure history (a) and simulated average void fraction (b) downstream the closed valve (P03).

The sample size of 500 was generated and RELAP model was computed on the sampled input parameters. As an example of the model results, pressure peak time plots of 20 random RELAP model runs are shown in Fig. 13. The sensitivity analysis by using CSM and CSV plots was performed by using MATLAB script. Two types of model outputs were considered: maximum value of

Table 3 Parameters selected for the sensitivity analysis and their distribution parameters. #

Parameter

Range of values

Mean value

Standard deviation

Probability distribution type

Min

Max

Initial conditions X1 Pressure at the pump header (Pa) X2 Water temperature (K)

3.871  105 290.7

4.029  105 296.6

3.95  105 293.65

3.95  103 1.475

Truncated normal Truncated normal

Assumptions X3 Valve closing rate (s  1) X4 Wall roughness (m) X5 Junction loss coefficient in pipe component ‘‘753’’ X6 Junction loss coefficient in pipe component ‘‘755’’ X7 Junction loss coefficient in pipe component ‘‘758’’

17.52 2.0  10  5 0.04476 0.01459 6.95  10  3

30.48 3.0  10  5 0.08313 0.02709 0.0129

24 (according to UMSICHT data) 25.0  10  6 (according to UMSICHT data) 0.06395 0.02084 9.92  10  3

3.24 2.5  10  6 9.59  10  3 3.13  10  3 1.49  10  3

Truncated Truncated Truncated Truncated Truncated

normal normal normal normal normal

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

pressure attained the time range of the modelling (1) and pressure at each time moment of the modelling within the first 4 s after closure of the valve (2). For the first type of model output only one model output value was computed for each model run. For the second type of model output in total 89 time moments (0– 4 s) were considered each representing a separate model output.

4.5 4

Pressure, MPa

3.5 3 2.5 2 1.5 1

69

7.1. Maximum pressure value as a model output In this case only one model output value was computed for each model run. The CSM and CSV plots are shown in Figs. 14 and 15, respectively. The CSM plot indicates that all parameters equally contribute to the mean of the model output, i.e. all might be fixed to point values without any significant effect on the mean of the model output. The CSV plot shows that variance of model output due to X6 (junction loss coefficient in pipe component ‘‘755’’) is not uniform and it can be effectively reduced by reducing the range of X6 values. The same effect, but less significant can be achieved with parameter X4 (wall roughness). All the other parameters are invariant and point estimates can be used with no effect on the model output mean and variance values. The CSM/CSV analysis was performed also for the time moment of the maximum pressure value as a model output, giving very similar results as for the pressure value discussed above. Another study [13] reported the first order indices and total indices computed by extended FAST (EFAST) method and are presented in Table 4. An interesting finding is that the CSV detects importance of X6 and X4 as well as EFAST, although CSM does not indicate any particular effect for this parameter.

0.5 0

7.2. Time dependent pressure value as a model output

3

3.1

3.2

3.3 3.4 Time, sec

3.5

3.6

3.7

1

1

0.8

0.8

0.8

0.6 0.4 0.2

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X1, q

0

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X2, q

0

1

1

0.8

0.8

0.8

0.6 0.4 0.2 0

CSM(q)

1

CSM(q)

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X4, q

0

1

0

0.2

0.4 0.6 0.8 Parameter X3, q

1

0

0.2

0.4 0.6 0.8 Parameter X6, q

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X5, q

1

0

0.2

0.4 0.6 0.8 Parameter X7, q

1

0

1 0.8 CSM(q)

CSM(q)

0

CSM(q)

1

CSM(q)

CSM(q)

Fig. 13. Results of 20 random model runs in the time range of 3–3.7 s.

Analysis of this type of model output is more complicated because of large number time points (89 in total) each being treated as an independent model output. In addition, the CSM and

0.6 0.4 0.2 0

Fig. 14. CSM plot for the maximum pressure value as model output. Definition of the parameters X1–X7 is presented in Table 3.

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

1

1

0.8

0.8

0.8

0.6 0.4 0.2

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X1, q

0

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X2, q

0

1

1

0.8

0.8

0.8

0.6 0.4 0.2 0

CSV(q)

1

CSV(q)

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X4, q

0

1

0

0.2

0.4 0.6 0.8 Parameter X3, q

1

0

0.2

0.4 0.6 0.8 Parameter X6, q

1

0.6 0.4 0.2

0

0.2

0.4 0.6 0.8 Parameter X5, q

1

0

0.2

0.4 0.6 0.8 Parameter X7, q

1

0

1 0.8 CSV(q)

CSV(q)

0

CSV(q)

1

CSV(q)

CSV(q)

70

0.6 0.4 0.2 0

Fig. 15. CSV plot for the maximum pressure value as model output. Definition of the parameters X1–X7 is presented in Table 3.

Table 4 The first order and total indices computed by EFAST as reported in [13]. Parameter

First order index Si

Total index STi

Interactions STi–Si

X1 X2 X3 X4 X5 X6 X7

0.13 0.01 0.01 0.18 0.02 0.54 0.02

0.25 0.11 0.10 0.28 0.11 0.64 0.08

0.12 0.10 0.10 0.10 0.09 0.11 0.06

Table 5 Significant parameters throughout the time range of 0–4 s of simulation. Definition of the parameters X1–X7 is presented in Table 3. Time range (s)

Most significant parameters by CSM

Most significant parameters by CSV

Qualitative comparison

0–2 3.0–3.17 3.18–3.20 3.21–3.24 3.25–3.28 3.29–3.30 3.31–3.34 3.35–3.36 3.37–3.40 3.41–3.43 3.44 3.45–3.50 3.6–4.0

None None Low importance of X1, X4, X6 All X1, X4, X6 X1,X4,X6 None X1,X4,X6 X1,X4,X6 X1, X2, X4, X6, X7 X1,X4,X6 Low importance of X1,X4,X6 None

X2 X1, X4, X6 All All X1, X4, X6 None None None X1,X4,X6 X1, X2, X4, X6, X7 All X1, X4, X6, X7 X2

X2 is important in the CSV X1, X4, X6 are important in the CSV CSV detects all to be important The same for both The same for both None important in the CSV The same for both None is important in the CSV The same for both The same for both CSV detects all to be important X7 is important in the CSV X2 is important in the CSV

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

CSV plots change significantly throughout the time range of the simulation (0–4 s). Analysis of each time point is summarized in Table 5. T = 0.1501sec

1

The results presented in Table 5 show the complexity of the model and very heterogeneous effects throughout the simulation time interval of 0–4 s after the closure of the valve 754 (Fig. 11).

T = 3.11sec

1 0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

T = 3.16sec

1

0

0.2

0.4

0.6

0.8

1

T = 3.185sec

1

0

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0.2

0.4

0.6

0.8

1

T = 3.2351sec

1

0

0

0.2

0.4

0.6

0.8

1

T = 3.2601sec

1

0

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

T = 3.3101sec

1

0

0.2

0.4

0.6

0.8

1

T = 3.3351sec

1

0

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0.2

0.4

0.6

0.8

1

T = 3.3851sec

1

0

0

0.2

0.4

0.6

0.8

1

T = 3.41sec

1

0

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

T = 3.4601sec

1

0

0

0.2

0.4

0.6

0.8

1

T = 3.485sec

1

0

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.2

0.8

1

0

0.6

0.8

1

0.4

0.6

0.8

1

0.8

1

0.8

1

0.8

1

0.8

1

T = 3.2851sec

0

0.2

0.4

0.6

T = 3.36sec

0

0.2

0.4

0.6

T = 3.4351sec

0

0.2

0.4

0.6

T = 3.7sec

1

0.8

0

0

1

0.8

0.4

T = 3.2101sec

1

0.8

0

0.2

1

0.8

0

0

1

0.8

0

T = 3.1351sec

1

0.8

0

71

0

0.2

0.4

0.6

Fig. 16. CSV plot for the parameter X6 at different time moments. The value q is plotted on x-axis and value of CSV(q) is plotted on y-axis.

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

The results indicate that during steady state conditions just before the water hammer caused pressure peak (the first 2 s of the simulation) and after it (later than 3.6 s) CSM does not detect any parameter to be important, but CSV detects X2 (water temperature) to be important in reducing the model output (pressure) variance. The simulation period of 3.0–3.6 s is full of heterogeneous effects as in this period the water hammer caused pressure peak is achieved in all simulations. In general view, parameters X1, X4 and X6 are most often appear among the most important. On the other hand, the rest of the parameters X2, X3, X5 and especially X7 are most often important in the CSV plot. Although the findings of the CSM and CSV are the same in certain periods of the simulation, the overall conclusion is that none of the parameters can be considered as unimportant. Dynamics of the CSV plot for the parameter X6 at different time moments is shown in Fig. 16. Not all 89 time moments are represented in the plot, but only each 5th time moment. The behaviour of the CSV changes from non-significant to significant and vice versa, the most important time periods being around 3.2 and 3.4 s. A number of effects are shown in Figs. 17–22 that present a snapshot of the time-dependent evolution. Fig. 17 presents a

0.12

0.1 Model output: pressure, MPa

72

0.08

0.06

0.04

0.02

0 3.86

3.88

3.9

3.92

3.94

3.96

3.98

4

4.02

4.04 x105

Parameter X1, Pa Fig. 19. Scatter plot of X1 at time point T¼ 3.20 s.

1 0.9

0.016

0.8 0.7 CSM(q), CSV(q)

Model output: pressure, MPa

0.014

0.012

0.01

0.6 0.5 0.4 0.3

0.008 0.2 CSM CSV diagonal

0.1

0.006

0

0.004 0.014

0.016

0.018

0.02

0.022

0.024

0.026

0

0.1

0.2

0.3

0.028

0.4

0.5

0.6

0.7

0.8

0.9

1

2.9

3

Parameter X1, q

Parameter X6

Fig. 20. CSM and CSV plot of X1 at time point T¼ 3.20 s.

Fig. 17. Scatter plot of X6 at time point T¼ 3.15 s.

4.5

1 4

CSM CSV diagonal

0.9

3.5 Model output: pressure, MPa

0.8

CSM(q), CSV(q)

0.7 0.6 0.5 0.4 0.3

2.5 2 1.5 1

0.2

0.5

0.1 0

3

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Parameter X6, q Fig. 18. CSM and CSV plot of X6 at time point T¼ 3.15 s.

0.9

1

0

2

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

Parameter X4, m Fig. 21. Scatter plot of X4 at time point T¼3.35 s.

x 10-5

S. Tarantola et al. / Reliability Engineering and System Safety 99 (2012) 62–73

The theoretical study of Ishigami function showed some new findings which are not visible by applying mean based methods like CSM or first order FAST indices, e.g. importance of the parameter X3 is neglected as having zero effect on the mean of Y, but it has a significant effect on the variance of Y, which can be quantified by using the CSV. The CSM and CSV techniques were applied to water hammer model case study. The CSV plots indicate a number of new findings compared to the CSM plots and provide the most effective ways to reduce variance of the model output. The results were also compared to a sensitivity study by using EFAST method. In a number of cases the results of this complex model study indicate that the CSM and CSV tools complement each other and both should be studied in a single framework.

1 CSM CSV diagonal

0.9 0.8 0.7 CSM(q), CSV(q)

73

0.6 0.5 0.4 0.3 0.2 0.1 0

Acknowledgements 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Parameter X4, q Fig. 22. CSM and CSV plot of X4 at time point T¼ 3.35 s.

scatter plot of the parameter X6 at time point of 3.15 s. The scatter plot indicates that variance is not equally distributed along the range of X6. This result is shown quantitatively by the CSV plot (Fig. 18), indicating that the highest 10% values of the X6 range increases the variance by 30%. This also shows the way how the variance can be reduced and how much (see Eq. (8). Note that the CSM plot does not indicate any particular effect for the parameter X6. Fig. 19 presents a scatter plot of the parameter X1 at time point of 3.20 s. The CSM plot (Fig. 20) does not show any particular importance of X1, but the CSV plot indicates that the variance is the largest for the lowest 10% values of the X1 range. The CSV plot more clearly indicates different behaviour within different intervals of the X1 range. Fig. 21 presents a scatter plot of the parameter X4 at time point of 3.35 s. In this case the CSV plot (Fig. 22) does not indicate any effect of X4, as the variance is almost equally distributed along the X4 range. However, the CSM plot shows that the model output mean is higher in the range of low X4 values and is lower in the range of higher X4 values. This shows the importance of X4 to the model output that cannot be observed by using only the CSV plot.

8. Conclusions and discussion The CSV technique was developed as an extension of the CSM plot. The definition of the CSV proposed in this article is able to provide information how much the variance as defined in (5) with respect to the full range mean can be reduced and what is the most effective way to do it.

This research was funded by Research Council of Lithuania (Grant no. ATE-10/2010).1 References [1] Sacks J, Welch WJ, Mitchell TJ, Wynn HP. Design and analysis of computer experiments. Statistical Science 1989;4(4):409–23. [2] Wegman EJ. Hyperdimensional data analysis using parallel coordinates. Journal of American Statistical Association 1990;85(411):664–75. [3] Kurowicka D, Cooke R. Uncertainty analysis with high dimensional dependence modelling. John Wiley & Sons; 2006. [4] Theus M, Urbanek S. Interactive graphics for data analysis. Principles and examples. CRC Press; 2009. [5] Saltelli A, Tarantola S, Camplongo F, Ratto M. Sensitivity analysis in practice. A guide to assessing scientific models. John Wiley & Sons; 2004. [6] Sinclair J. Response to the PSACOIN Level S exercise. PSACOIN Level S intercomparison. Nuclear Energy Agency, OECD; 1993. [7] Bolado-Lavin R, Castaings W, Tarantola S. Contribution to the sample mean plot for graphical and numerical sensitivity analysis. Reliability Engineering and System Safety 2009;94:1041–9. [8] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global sensitivity analysis. The primer. John Wiley & Sons; 2008. [9] Ishigami T, Homma T. An importance quantification technique in uncertainty analysis for computer models. In: Proceedings of the ISUMA’90: first international symposium on uncertainty modelling and analysis, p. 398–403. [10] Homma T, Saltelli A. Importance measures in global sensitivity analysis of nonlinear models. Reliability Engineering and System Safety 1996;52:1–17. [11] Dudlik A. Pilot plant pipework test facility and water hammer benchmark experiments at Fraunhofer UMSICHT, Oberhausen hydraulic systems. In: Proceedings of the NURETH11 conference; 2005. ¨ [12] Dudlik A, Schonfeld SBH, Hagemann O, Fahlenkamp H. Water hammer and cavitation hammer in process plant pipe systems. Kerntechnik 2003;68:91–6. [13] Kaliatka A, Kopustinskas V, Vaiˇsnoras M. Water hammer model sensitivity study by the FAST method. Energetika 2009;55(1):13–9. [14] Kaliatka A, Uˇspuras E, Vaiˇsnoras M. Uncertainty and sensitivity analysis of parameters affecting water hammer pressure wave behaviour. Kerntechnik 2006;5–6:270–8. [15] Kaliatka A, Uˇspuras E, Vaiˇsnoras M. Benchmarking analysis of water hammer effects using RELAP5 code and development of RBMK-1500 reactor main circulation circuit model. Annals of Nuclear Energy 2007;34(1–2):1–12.

1

Regards only authors from Lithuanian Energy Institute.