Full text
Abstract— This article presents a stochastic computational model for the analysis of the reliability of a drawn steel bar. The whole distribution of the limit state function is studied using global sensitivity analysis based on Cramér-von Mises distance. The algorithm for estimating the sensitivity indices is based on one loop of the Latin Hypercube Sampling method in combination with numerical integration. The algorithm is effective due to the approximation of resistance using a threeparameter lognormal distribution. Goodness-of-fit tests and other comparative studies demonstrate the significant accuracy and suitability of the three-parameter lognormal distribution, which provides better results and faster response than sampling-based methods. Global sensitivity analysis is evaluated for two load cases with proven dominant effect of the long-term variation load action, which is introduced using Gumbel probability density function. The Cramér-von Mises indices are discussed in the context of other types of probability-oriented sensitivity indices whose performance has been studied earlier. Keywords— Global sensitivity analysis, steel, probability, failure, reliability, random sampling. I. INTRODUCTION eliability is described as the ability of a system or component to function under stated conditions for a specified period of time [1, 2]. Structural reliability is assessed using methods of probabilistic analysis [3]. Structural sensitivity analysis is a suitable complement to probabilistic reliability analysis [4]. The basic measure of reliability is the probability that failure of a load-bearing structure does not occur [5]. The most serious failure is loss of load-carrying capacity of a component or member within a structure or of the structure itself, see for e.g. [6]. Failure of a structure occurs when the material in the structure is stressed to its strength limit [7]. Although stresses cannot be directly investigated experimentally, limit states can be investigated by measuring permanent deformations or observing fractures [8]. Load-bearing structures must be designed with such geometric and material characteristics that the probability of failure is very low. For reliable design, probability of 7.2E-5 or less is commonly required [5]. A specific type of failure of slender steel structures is loss of stability, which causes the structure to collapse before reaching its material strength limit [9]. In stochastic systems, stability often means insensitivity or low sensitivity of their output characteristics to the shapes of some input distributions [10]. In operations research, sensitivity analysis is developed as a method of critical assessment of decisional variables, which have interrelations among them, and is capable of identifying those sensitive variables that have an influence on the final desired result [11-12]. One of the principal problems in system reliability studies is the problem of their sensitivity to the shape of input distributions [13-15]. The desired result may be the validation of the attributes of static [16] or dynamic [17] models that require simultaneous optimization, but also other performance characteristics associated with the structure, see for e.g. [18-19]. Assessment of decisional variables can be based on MCDM [20] and computer experiments [21]. The basic calibration tools of reliability according to standards for structural design are the partial safety factors of load and resistance [22-25]. The verification of reliability of standard design is commonly performed using probabilistic methods [26, 27], but only deterministic types of sensitivity studies, see for e.g. [28, 29]. The task of sensitivity analysis in Limit states of structures and global sensitivity analysis based on Cramér-von Mises distance Z. Kala, Department of Structural Mechanics, Faculty of Civil Engineering, Brno University of Technology Veveří Str. 95, Brno, ZIP 602 00 Czech Republic Received: May 5, 2020. Revised: June 27, 2020. Accepted: July 3, 2020. Published: July 7, 2020. R INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 107
stochastic models oriented to structural reliability is to quantify the effects of input variables on the output characteristics, the choice of which is specific to the given type of structure and computational model, see e.g. [30-35]. According to conventional definition, reliability sensitivity analysis focuses on measuring the effects of uncertainty or distribution parameters of input variables on failure probability [36]. Sensitivity analysis subordinated to contrasts [37] appears to be an appropriate reliability-oriented method, which can be used to analyse both the probability of failure [38-39] and design quantiles [40-42]. Another method discussed is global sensitivity analysis based on Cramér-von Mises distance [43], which is focused on the distribution function of the response. Other types of sensitivity analysis based on measuring the absolute difference between the unconditional and conditional failure probability were published in [44-45]. Borgonovo moment independent importance measure [46] was not developed directly for the analysis of failure probability but can be mentioned as one of the few methods that permits the analysis of statistically dependent input random variables. The analysis and comparison of different techniques of sensitivity analysis and a discussion on their strengths and weaknesses has been published in [47]. If an analytical solution is not available, the failure probability can be estimated using methods such as Monte Carlo, Importance Sampling [48], Adaptive Sampling [49], Response Surface [50], artificial neural networks [51] or polynomial chaos expansions [52]. Unfortunately, a plan tailored to estimate a particular type of sensitivity index, or set thereof, may be suboptimal for other indices. In this paper, the effects of input random variables on the distribution function of the limit state are studied using Cramér-von Mises indices. One of the objectives of this paper is to develop a stochastic computational model of a load-bearing element that will facilitate the efficient evaluation of global sensitivity analysis and which will be used in further research. II. STOCHASTIC MODEL AND RELIABILITY ANALYSIS The aim of the study is to analyse the theoretical reliability of a steel structure, which is at the design stage. This is primarily ensured by design reliability conditions of standard EN 1990 [5]. A. Load Action Let a drawn bar be subjected to a combination of permanent load action G and long-term leading single variable load action Q. Permanent load action G can be considered with mean value µ G, which is equal to the characteristic value Gk [53]. The random deviations of the permanent load action from the mean value should be symmetrical, which corresponds well with the Gauss probability density function, see Fig. 1. The standard deviation σ G is selected under the assumption that approximately 50% of random realizations of G are found in the interval ± σ G. Fig. 1 Density probability functions of permanent load actions The leading variable load action is introduced with Gumbel probability density function [54], see Fig. 2, which is the characteristic function for extreme value distributions [55, 56]. The cumulative distribution function of Gumbel distribution is ( ) − −−= β α x xF expexp , (1) for −∞ < x < ∞, −∞ < α < ∞ and β > 0. Parameters α and β can be expressed from βγαµ += Q , (2) 6 πβ σ = Q , (3) where γ=0.5772156649 is the Euler–Mascheroni constant. The mean value µ Q and standard deviation σ Q are introduced as functions of the characteristic value Qk using the equations µ Q=0.6 Qk, σ Q=0.21 Qk, which were used in [53]. Fig. 2 Density probability functions of variable load actions It can be noted that Gauss and Gumbel probability density functions are distributions recommended in standard EN1990 [5] and can be used in the studies presented here, see Table I. INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 108
TABLE I INPUT RANDOM VARIABLES OF LOAD ACTION Characteristic Pdf Mean value St. deviation Permanent load G Variable load action Q Gauss Gumbel-max G k 0.6 ⋅ Qk 0.1⋅G k 0.21⋅Qk However, other types of functions may also be used in justified cases or specific load cases. A detailed description of these and other less frequently used theoretical models of random variables can be found in International Standards ISO 13822, ISO 12491, ISO 3534-1 and specialized literature [57]. B. Resistance The drawn bar is made from a steel sheet of length L, width b and thickness t. Width and thickness are the dimensions of the rectangular cross-section. The resistance R is expressed as the product of the yield strength fy and cross-sectional area A of the rectangular cross-section. AfR y ⋅= , (4) where btA ⋅= . (5) The yield strength, which is generally a random variable, is an important material characteristic of steel. In the past, numerous experimental tests on the material and geometric characteristics of structural steel have been performed [58], which provide important background information for the analysis of resistance. The input random variables are listed in Table II. TABLE II INPUT RANDOM VARIABLES OF RESISTANCE PARAMETERS Characteristic Pdf Mean value St. deviation Yield strength f y Specimen thickness t Specimen width b Gauss Gauss Gauss 395.68 MPa 16 mm 80 mm 25.13 MPa 0.734 mm 0.8 mm Input variable A has an asymmetric chi square probability density function. Resistance R is the product of the yield strength and the cross-sectional area given by Eq. (5). If we multiply three random variables fy⋅t⋅b with Gauss probability density functions, the computed resistance R does not have a Gauss probability density function (has an asymmetric distribution) and it is necessary to seek the most suitable distribution for the approximation of the distribution of R. It can be noted that the probability density distribution of R can be analysed by the numerical (but not analytical) integration of a triple integral or by approximating the Gauss probability density function from Table II using such functions that can be integrated analytically. This approach would have to be used if none of the known types of probability density functions proved to be suitable for approximating the probability density function of R. A suitable approximation of the probability density function of random variable R can be selected from numerous types of probability density functions among which certain types can be fitted more closely to the observed frequency of the resistance than others, depending on the characteristics of the phenomenon and of the distribution. One of the possibilities of approximating R in the case study analysed herein is by using the three-parameter lognormal distribution. The software Statrel [59] refers to the three-parameter lognormal distribution, which will be used below, as the shifted lognormal probability density function. ( ) ( ) ( )( ) ⋅ −− − ⋅⋅− = 2 2 0 0 2 ln exp 2 1 u u u S mxx Sxx xf π , (6) for x0 < x < ∞, where xx tx σµ ⋅+= 00 , (7) += 2 0 1 1ln t S u , (8) 2 ln 2 0 u xu S tm −⋅= σ , (9) where x x xxxx a a azaz t1 2 22 3 32 32 0 − ⋅ ++−−− = , (10) where ax is standardized skewness and 4 2+= xxx aaz . (11) C. Probability Analysis of Resistance The selection of the suitable probability density function of resistance was verified numerically by comparing the results of two probabilistic studies. It was tested whether the shifted lognormal probability density function gives a close fit and leads to good predictions. The first probability IPfR is calculated by numerical integration using the shifted lognormal probability density function. The second probability IIPfR is calculated by numerical integration of the function R=fy ⋅ t ⋅ b in threedimensional space. The value IIPfR can be described as a very good approximation of the “real” probability that could theoretically be achieved if the integration step approaches zero. ( ) ifR rRPP <= , (12) The suitability of the shifted lognormal probability density function is verified by comparing the values of IPfR and IIPfR. The inaccuracies of the estimates of both IPfR and IIPfR are INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 109
negligible because numerical integration was performed with 7000000 (IPfR) and 70003 steps, see the variable max in the algorithm below. Throughout the article, numerical integration is performed on the interval mean value plus minus ten standard deviations of the integrated variable. The introduction of wider intervals has a negligible effect on increasing the accuracy of the results of the numerical integrations. Three nested loops are used in the algorithm of the numerical integration for the computation of IIPfR. The basic form of the numerical integration algorithm can clearly be written in the programming language Pascal as: function NormR2(x,m,S:extended):extended; begin NormR2:=exp(-sqr(x-m)/(2*sqr(S)))/(sqrt(2*pi)*S); end; {Three nested loops numerical integration} Function Pf(R:Extended):Extended; var i1,i2,i3,max :integer; a1,a2,a3,b1,b2,b3 :extended; PfR,Sm10,R04,R01,R02,R03 :extended; m1,m2,m3,S1,S2,S3 :extended; begin max:=7000; // steps of numerical integration m1:=395.68;m2:=16;m3:=80; S1:=25.13;S2:=0.734;S3:=0.8; R:=200E3; // R is value of resistance Sm10:=10; R01:=m1-Sm10*S1; R02:=m2-Sm10*S2; R03:=m3-Sm10*S3; b1:=R01; PfR:=0; // PfR is probability of R for i1:=1 to max do begin a1:=b1;b1:=R01+(m1+Sm10*S1-R01)*i1/max; b2:=R02; for i2:=1 to max do begin a2:=b2;b2:=R02+(m2+Sm10*S2-R02)*i2/max; R04:=(b1-a1)*(b2-a2); b3:=R03; for i3:=1 to max do begin a3:=b3;b3:=R03+(m3+Sm10*S3-R03)*i3/max; if (a1+b1)*(a2+b2)*(a3+b3)/8<R then PfR:=PfR+R04*(b3-a3)*NormR2((a1+b1)/2,m1,S1)* NormR2((a2+b2)/2,m2,S2)*NormR2((a3+b3)/2,m3,S3); end; end; end; Pf:=PfR; // Pf is probability of R end; The computational complexity of the algorithm is 1012 steps of mainly arithmetic operations. The algorithm runs for 70 minutes on a single core CPU Intel Core i7-3740QM processor. It can be noted that a similar algorithm is applied in the DOProC method [60], which, however, has input random variables defined by histograms. The assumption of using three nested loops to estimate IIPfR is demanding on CPU time even though it is a relatively simple algorithm. The estimation of IIPfR is more numerically demanding the higher the value of IIPfR. Therefore, for IIPfR>0.5 it is possible to speed up the algorithm by modifying the condition into the form “(a1+b1)*(a2+b2)*(a3+b3)/8>=R“ and “Pf:=1-PfR;“, which is a key optimization. Additional optimizations were also performed. The computation of IPfR using the shifted lognormal probability density function is not so numerically demanding and was therefore performed using seven million steps. The computational complexity of the calculation of IPfR is approximately one hundred thousand times lower than the calculation of IIPfR. Table III contains a comparison of IIPfR values with IPfR values. From a numerical point of view, it is apparent that the estimate of IPfR is sufficiently accurate because the differences between the estimates of IPfR and IIPfR are minimal, with the exception of very small probabilities. Conformity in exponents with the power of E-19 is fully sufficient in terms of engineering probability of failures. Possible limitations pertain to probabilities with exponents lower than E-19, which, however, is irrelevant in engineering applications. TABLE III COMPARISON OF PROBABILITY COMPUTATIONS Resistance ri IPfR Shifted lognormal pdf IIP fR Numerical integration using three nested loops 0 100 kN 200 kN 300 kN 400 kN 500 kN 600 kN 700 kN 2.6215E-65 2.7419E-37 1.3556E-19 5.8430E-09 2.7029E-03 4.4285E-01 9.8810E-01 9.99996E-01 0 3.1258E-47 1.1362E-19 5.2958E-09 2.6880E-03 4.4288E-01 9.8812E-01 9.99996E-01 The resistance must be positive. This is one of the arguments that sometimes calls into question the suitability of approximating resistance using the shifted lognormal probability density function in the event that it is also defined for negative values. The shifted lognormal probability density function used in this paper is defined on the interval (-579.377 kN; ∞). It is evident from Table III that the probability that the resistance has a negative value is IPfR=2.6215E-65. This is a negligible value. If we introduce a trimmed shifted lognormal probability density function, the value of the obtained probability is practically the same. This implies that trimming is not necessary. The shifted lognormal probability density function can be used as it is. Table IV contains a comparison of IIPfR with probability IIIPfR obtained by approximation using Gauss probability density function, which is not considered appropriate. It is apparent that the estimate IPfR (Shifted lognormal probability density function) is significantly more accurate in comparison INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 110
with IIIPfR, see Table IV. Estimates of IIIPfR obtained by approximation using Gauss probability density function have discrepancies from the target IIPfR for small probabilities. TABLE IV COMPARISON OF PROBABILITY COMPUTATIONS Resistance ri IIIPfR Gauss pdf IIP fR Numerical integration using three nested loops 0 100 kN 200 kN 300 kN 400 kN 500 kN 600 kN 700 kN 5.4609E-37 1.5929E-24 9.6066E-15 1.2498E-07 3.9107E-03 4.3580E-01 9.9027E-01 9.999999E-01 0 3.1258E-47 1.1362E-19 5.2958E-09 2.6880E-03 4.4288E-01 9.8812E-01 9.99996E-01 The comparative studies confirmed that the shifted lognormal probability density function is a suitable approximation of resistance computed from Eq. (4). The following chapter validates this conclusion using goodness-offit tests. D. Sampling Based and Statistical Analysis of Resistance The goodness-of-fit test can be used to test if sample data fit a distribution from a random observation. The random observations were obtained using the Latin Hypercube Sampling Method (LHS) [61, 62]. One million steps of the LHS method were used to simulate the random realizations of R. The use of the LHS method for such a large number of samples has proven appropriate since the results obtained using the LHS method are always better than the results obtained using the simple Monte Carlo method, including the large numbers of runs, which is useful in global sensitivity analyses of the limit state function distribution [38]. Fig. 3 Approximation of the histogram of resistance The results of the statistical analysis of resistance R obtained using the LHS method are shown in Table V. Three variants of the approximation of the histogram of resistance computed using the input random variables listed in Table II are shown in Fig. 3. Arithmetic mean µ R, standard deviation σ R and standard skewness aR are important results of the statistical analysis, which are also applicable as parameters of the tested shifted lognormal probability density function. The values of µ R, σ R and aR can be computed analytically from the function Y 321 XXXY ⋅⋅= , (13) of three random variables X1, X2, X3 with Gauss probability density functions. Each random variable Xi is characterized by its arithmetic mean µ i and standard deviation σ i. Arithmetic mean µ Y, standard deviation σ Y and standard skewness aY of random variable Y are then expressed by the following equations: 321 µµµµ ⋅⋅= Y , (14) ( ) γσσγσµγσσµµσ ⋅⋅+⋅⋅+⋅+⋅⋅= 2 2 2 1 2 1 2 2 2 2 2 3 2 2 2 1Y , (15) where 2 3 2 3 σµγ += , (16) ( ) .4 6 2 3 2 2 2 1 2 3 2 2 2 1 2 3 2 2 2 1 2 3 2 2 2 1 3 σσσ µσσσµσσσµ σ µ ⋅⋅⋅+ ⋅⋅+⋅⋅+⋅⋅⋅⋅= Y Y Y a , (17) The statistical characteristics of resistance R Eq. (4) computed in this study using Eq. (14), (15), (17) and Table II are µ R=506.47 kN, σ R=40.031 kN and aR= 0.11065. The comparison of µ R, σ R and aR with the results in Table V shows perfect agreement between the analytical solution and the statistical analysis. TABLE V STATISTICAL CHARACTERISTICS OF HISTOGRAM IN FIG. 3 Characteristic Value Valid observations Minimum Maximum Range Median Arithmetic mean Geometric mean Mean square Variance Stand. deviation Coef. of variation Third moment Stand. skewness Fourth moment Stand. kurtosis Variance of mean Var. of variance Var. of 3. moment Var. of 4. moment 1000000 323.62 710.63 387.01 505.73 506.47 504.89 1602.4 1602.4 40.030 0.79037E-01 7041.1 0.10977 0.77691E+07 3.0257 0.16023E-02 0.60313E+08 26005. 0.68249E+09 INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 111
E. Goodness-of-fit Tests of Resistance Goodness-of-fit tests are statistical tools that can reliably reject approximations that are not perfectly adapted to data obtained using the LHS method. Measures of goodness-of-fit typically summarize the discrepancy between observed values and the values expected under the model in question. Goodness-of-fit tests may reject unsuitable types of approximations with a selected probability, but they cannot confirm the appropriate types of probability density functions. Thus, the aim of the study can only be to reject inappropriate types of probability density functions for the approximation of the resistance. The tested models in question are shifted lognormal, Gauss, Lognormal, Hermite, Gamma and Beta probability density functions. Chi-square, KolmogorovSmirnov and Anderson-Darling goodness-of-fit tests are used, see Table VI. Testing for the shifted lognormal probability density function is depicted in Fig. 4, Fig. 5 and Fig. 6. The tests were performed using the program Statrel 3.10. Fig. 4 Chi-square goodness-of-fit test [59] of shifted Lognormal pdf Fig. 5 Kolmogorov-Smirnov goodness-of-fit test [59] of shifted Lognormal pdf Based on the results of three goodness-of-fit tests, it can be concluded that the shifted lognormal probability density function cannot be rejected for the approximation of the probability density function of R. Fig. 6 Anderson-Darling goodness-of-fit test [59] of shifted Lognormal pdf It is evident from the results in Table VI that the conclusions of Chi-square and Kolmogorov-Smirnov tests are the same and are consistent with the results of the studies in the previous chapters. The Anderson-Darling test did not reject the hypothesis in any case, even though, for example, in Table IV it was clearly shown that the choice of Gauss probability density function is not suitable for the approximation of the probability density function of R. In addition to the shifted lognormal probability density function, four-parameter Hermite probability density function could also be considered. TABLE VI GOODNESS-OF-FIT TEST OF HISTOGRAM IN FIG. 3 Probability density function Chi-square test KolmogorovSmirnov test AndersonDarling test Shifted Lognormal Gauss Lognormal Hermite (Statrel) Gamma Beta Not rejected Rejected Rejected Not rejected Rejected Rejected Not rejected Rejected Rejected Not rejected Rejected Rejected Not rejected Not rejected Not rejected Not rejected Not rejected Not rejected III. PROBABILITY ASSESSMENT OF RELIABILITY Reliability can be understood as the ability of a structure or a structural member to fulfil stated requirements for which it has been designed. Reliability is generally expressed in terms of probability. The ultimate limit states are defined as states that are associated with collapse or with similar forms of structural failure. Loads higher than resistance can be understood as random failure, which occurs in dependence on random material and geometric characteristics of structural members and load actions. The probability of failure can be written as: ( ) RQGPP f >+= , (18) where G, Q and R are random variables defined in the previous chapter. The statistical characteristics of random variables G and Q are defined in Table I and are functions of INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 112
parameter δ as mentioned below, see also [53]. The statistical characteristics of the resistance function R=fy·t·b are defined in Table I. Practically Eq. (13) can be evaluated in the form ( ) 0<−= QDPP f , (19) where distance D=R-G is considered as a random variable with shifted lognormal probability density function Eq. (6). Arithmetic mean µ D, standard deviation σ D and standard skewness aD are expressed by the following equations: GRD µµµ −= , (20) 22 GRD σσσ += , (21) 3 3 D R RD aa σ σ ⋅= , (22) The use of the shifted lognormal probability density function for DRG is justified similarly as for R in chapter II. It can be shown that if the shifted lognormal probability density function Eq. (6) is suitable for R, then it is also suitable for D, because a lower value of skewness (if σ D> σ R then aD<aR) leads to a better approximation according to Eq. (6). The probability of failure Eq. (19) can be computed numerically from the integral: ( ) ( ) ∫ ∞ ∞− Φ= yyyP QDf d ϕ , (23) where ϕ Q(y) is Gumbel probability density function, ΦD(y) is cumulative distribution function of the shifted lognormal distribution and y denotes a general point of the observed force variable, through which both variables D and Q are expressed. Integration in Eq. (23) is performed numerically by Simpson’s rule using more than ten thousand integration steps. Numerical integration Eq. (23) allows relatively accurate computation of very small values of Pf in situations where computation using Monte Carlo numerical simulation methods is extremely numerically demanding or impossible. The methodology for the verification of standard design reliability is described in detail, for example, in [53]. Mean values µ G, µ Q and standard deviations σ G, σ Q of random variables G and Q are functions of characteristic values Gk and Qk according to Eq. (24) and Eq. (25), where δ is the parameter of the study. kk k QG Q + = δ , (24) dkk RQG =⋅+⋅ 5.135.1 , (25) where Rd=355E6·16E-3·80E-3=454.4E3N is the design resistance according to EUROCODE 3. The values 1.35 and 1.5 are partial safety factors of permanent and long-term variable load action [53]. Table VII contains the results of probabilistic analysis computed in relation to parameter δ , which has a step of 0.1. In Table VII minimum Pf occurs for δ =0.2 and maximum for δ =1.0. It is evident that the proportions between permanent and long-term variable load action strongly influence Pf. One way to optimize Pf for different load combinations is the calibration of partial safety factors. However, this would require a more complex standard assessment, which is debatable. TABLE VII RESULTS OF FAILURE PROBABILITY ANALYSIS δ µ G σ G µ Q σ Q P f 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 336.6 299.6 263.4 228.0 193.4 159.4 126.2 93.7 61.8 30.6 0.0 33.7 30.0 26.3 22.8 19.3 15.9 12.6 9.4 6.2 3.1 0.0 0.0 20.0 39.5 58.6 77.3 95.7 113.6 131.2 148.4 165.2 181.8 0.0 7.0 13.8 20.5 27.1 33.5 39.8 45.9 51.9 57.8 63.6 0.000439 0.000068 0.000022 0.000029 0.000065 0.000135 0.000248 0.000406 0.000604 0.000839 0.001105 IV. GLOBAL SENSITIVITY ANALYSIS The aim of global sensitivity analysis is to quantify the effects of input random variables on the limit state function taking into account the whole distribution of random variables. One approach to estimating these effects is through the application of Cramér-von Mises indices [44]. The adaptation of Cramér-von Mises indices to the analysis of whole distribution of random variables can be performed using Eq. (19) or Eq. (23). The estimation of the first-order indices is based on the computation of probability Φ DQ( τ ) and conditional probability ( ) τ DQ Φ when input Xi is fixed [44]. Input random variables G, Q, fy, t2, b are assumed to be statistically independent. ( ) ( ) ττ ≤−=Φ QDP DQ , (26) ( ) ( ) i DQ XQDP ττ ≤−=Φ , (27) ( ) ( ) ( ) [ ] ( ) ( ) ( ) ( ) τ ττ ττ DQ RDQDQ DQ DQ i E SΦ Φ−Φ Φ−Φ =∫d 1 2 . (28) where the Cramér-von Mises distance is ( ) ( ) ( ) [ ] ( ) τττ DQ R DQ DQCVM ED ΦΦ−Φ= ∫ d 2 . (29) Estimates of Eq. (26) and Eq. (27) are computed by numerical integration. The mean value E[⋅] is estimated using INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 113
500 steps of the LHS method. Higher-order sensitivity indices Sij, Sijk, etc., are computed similarly, however, by fixing two, three, etc., input variables [43]. If two of the three random variables of R Eq. (4) are fixed, then the approximation using the shifted lognormal probability density function is not needed and the Gauss probability density function is used. Sensitivity indices Si, Sij, Sijk, etc., are based on Hoeffding decomposition and thus the sum of all indices must be equal to one [43]. It can be noted that the numerical complexity of the estimates of Cramér-von Mises indices can be reduced by using polynomial chaos expansions [52] as presented in [63]. In the study presented here, five input random variables yield five first-order sensitivity indices, and this leads to a total of 25-1 sensitivity indices; therefore, 31 sensitivity indices are estimated. Computer parallel processing can be employed to evaluate sensitivity indices using the LHS method. The computation of a sensitivity index is a random realization of this index when the LHS method is used. The use of new (pseudo-) random numbers (new random realizations of input random variables) gives a new random realization of the sensitivity index. The theoretically accurate value of the sensitivity index can be obtained (theoretically) either (i) using an infinite number of steps of the LHS method or (ii) by evaluating an infinite number of random realizations of the sensitivity index (using a finite number of steps of the LHS method to compute one realization of the sensitivity index). In approach (ii), the precise value of the sensitivity index can be obtained as the arithmetic mean (theoretical) of an infinite number of random realizations of the sensitivity index. Although we do not have the possibility to perform an infinite number of random realizations of sensitivity indices, a finite number can be realized, where each index can be evaluated on one CPU core. The advantage of approach (ii) is that the random realizations of sensitivity indices can be computed parallelly on the computer. In practice, the estimate of each index is computed using eight random realizations, which are computed on eight CPU cores for eight different sets of (pseudo-) random numbers (LHS runs). The resulting estimate of sensitivity indices is computed as the arithmetic mean of these eight estimates, resulting in higher accuracy. The sum of all numerically estimated indices is 1.000664, which is one of the indicators of sufficient accuracy of the results. Cramér-von Mises sensitivity indices are computed for δ =0.5 and δ =0.7, see Fig. 7 and Fig. 8. The first-order indices account for more than 50% of the influence on the output distribution in both cases. The results show the dominant effect of long-term variable load action Q and yield strength fy. The long-term variable load action Q has a dominant effect even in the case of δ =0.5, Gk =Qk = 159.4 kN, see Fig. 7. The cause of the dominance is the variation coefficient of Q, which has a value of 0.35, while G has a variation coefficient of 0.1. Another cause of the dominance of Q is Gumbel probability density function, which is used in probability theory and statistics to model the distribution of the maximum number of samples of different distributions [55, 56]. In the case of δ =0.7, the dominance of Q is more pronounced in both the first-order effects and in interactions with other variables, see Fig. 8. The second most influential variable is the yield strength, see Fig. 7 and Fig. 8 Conversely, the specimen width b has practically no significant influence. Introducing b as a deterministic variable anywhere in the domain would have a much smaller effect on limit state distribution than if we did the same with any other input random variable. Fig. 7 Cramér-von Mises indices for δ =0.5 Sobol’s decomposition is a proven and widely used method of global sensitivity analysis with many useful applications, see e.g. [64-67]. The sensitivity analysis results have been used as the basis for decision making on the modification of technological procedures for manufacturing beams [68]. Although these are methods of different types, the comparison of quantile-oriented sensitivity analysis and Sobol' sensitivity analysis showed that there are no large differences between the order of importance of the input variables, but certain differences in the size of the indices exist [42]. It can be discussed whether Cramér-von Mises indices are suitable fsor studying structural reliability. Estimates of Cramér-von Mises indices are based on integration over dΦ DQ( τ ) using Eq. (26), (27), i.e., for all τ belonging to R. As a result of integration over all τ , change in the mean value of D or Q is not reflected in the values of the indices. The Cramér-von Mises sensitivity analysis takes into account the whole distribution of random variables. Thus, Cramér-von Mises indices do not examine the effects on one Pf value, but examine the effects on all Pf values from all τ . INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 114
Fig. 8 Cramér-von Mises indices for δ =0.7 Contrast-based sensitivity indices [37], which examine the effects of input variables on one target Pf and are directly related to the assessed reliability, do not have this shortcoming [38-39]. Contrast-based sensitivity indices belong to the category of so-called reliability-oriented sensitivity analysis [69]. It can be noted that contrast-based sensitivity indices [37-39] can be written in a very transparent form Eq. (30), but only if the values of Pf, PfXi etc. are very small. It can be noted that Eq. (30) has no integration over τ , respectively, t is still zero. ( )( ) ( ) ( )( ) ( ) ( )( ) ,smallveryforonly ,1,1 1 11 0 0 2 0 0 0 f XZ Z Z XZ Z f iff i PEcorr E EEE P XPEP C i i < < < < <= − = − ≈ . (30) where Z=D-Q, if Z<0 then failure occurs, otherwise there is no failure, and corr is Pearson correlation coefficient. Numerical experience with the LHS method shows that the convergence of Cramér-von Mises indices [43] is much better than the convergence of probability-oriented contrast sensitivity indices [37]. The sequence in Eq. (29) is strongly convergent as Monte Carlo runs go to infinity and is asymptotically Gaussian. In contrast, the convergence of Ci Eq. (30) and all other higher order contrast-based sensitivity indices is very slow and requires an extreme number of sampling runs or direct integration calculations. It can be noted that algorithms for effective estimations of contrast-based sensitivity indices [37] using polynomial chaos expansions [52] or neural networks [70] would be desirable so that more complex problems could be solved. V. CONCLUSION The methodology of stochastic reliability analysis, which consists of the probabilistic analysis of failure and global sensitivity analysis, is described in this article. Global sensitivity analysis is based on the computation of the Cramérvon Mises indices whilst taking into account the whole distribution of random variables. In the presented case study, resistance is a function of the product of three random variables with Gauss probability density functions (yield strength, specimen thickness, specimen width), load action is a function of the sum of two random variables (permanent load action, variable action). Arithmetic mean, standard deviation, and standard skewness of resistance are calculated analytically using the mean values and standard deviations of yield strength, specimen thickness and specimen width. It is shown that the use of the shifted lognormal probability density function is suitable for the approximation of the relative frequency of random resistance. This was confirmed by comparing the results of statistical and probabilistic analyses with the results of numerical integration and goodness-of-fit tests. Eight deterministic values of the loadcarrying capacity are selected and the probability that the random resistance is less than the selected values is computed eight times. The probability obtained using the shifted lognormal probability density function is in perfect agreement with the result of numerical integration, which from an engineering point of view can be considered as a practically accurate (but also very numerically demanding) solution. Gauss probability density function did not show such agreement. Goodness-of-fit tests have shown that the use of the shifted lognormal probability density function for the approximation of the probability density of resistance is better than using the Gauss Lognormal, Gamma or Beta probability density function. Alternatively, the Hermite probability density function can also be considered, which could be explored in the future. These conclusions were made for the case study presented here. Computation of the probability of failure is based on the approximation of resistance using a three-parameter lognormal distribution, which is referred to as the shifted lognormal probability density function. The presented solution is highly effective for the case study presented in this paper. In the presented case study, the advantages of the presented method are the very fast and accurate estimation of the probability of failure. The stochastic model based on the shifted lognormal probability density function and Gumbel probability density function gives very good accuracy of the estimate of the probability of failure and at the same time has relatively low demands on CPU time of the computer. Both of these properties are desirable in reliability-oriented sensitivity INTERNATIONAL JOURNAL OF MECHANICS DOI: 10.46300/9104.2020.14.14 Volume 14, 2020 ISSN: 1998-4448 115