# Response Surface Method for Material Uncertainty Quantification of Infrastructures

Authors: Mohammad Amin Hariri-Ardebili; S. Mahdi Seyed-Kolbadi; Mohammad Noori
Journal: Shock and Vibration
 (2018)

---

## Abstract

Recently, probabilistic simulations became an inseparable part of risk analysis. Managers and stakeholders prefer to make their decision knowing the existing uncertainties in the system. Nonlinear dynamic analysis and design of infrastructures are affected by two main uncertainty sources, i.e., epistemic and aleatory. In the present paper, the epistemic uncertainty is addressed in the context of material randomness. An old ultra-high arch dam is selected as a vehicle for numerical analyses. Four material properties are selected as random variables in the coupled dam-reservoir-foundation system, i.e., concrete elasticity, mass density, compressive (and tensile) strength, and the rock modulus of elasticity. The efficient Box-Behnken experimental design is adopted to minimize the required simulations. A response surface metamodel is developed for the system based on different outputs, i.e., displacement and damage index. The polynomial-based response surface model is subsequently validated with a large number of simulations based on Latin Hypercube sampling. Results confirm the high accuracy of proposed technique in material uncertainty quantification.


---

## Body

## 1. Introduction

With the recent advances in computational tools, the probabilistic numerical simulations became an important aspect in risk analysis and risk management. Nowadays, decision-making is based on the uncertainties in the system and not a deterministic simulation. Nonlinear dynamic response of concrete dams can be affected by two main uncertainty sources: epistemic and aleatory [1].The aleatoric uncertainty stems from intrinsic randomness of a phenomenon. Uncertainty in the seismic hazard (e.g., intensity, time, and return period) is the most dominant one. On the other hand, the epistemic uncertainty is due to lack of knowledge. Uncertainty in material characteristics (e.g., modulus of elasticity and strength) is the main source in this category. Both these uncertainties can be incorporated in the numerical simulations, which results in the uncertainty propagation though the model, Figure1. Ground motion record-to-record variability is usually presented as a fragility function. A comprehensive state-of-the-art review on the fragility analysis of concrete dam can be found in [2]. Material uncertainties were also studied in few cases [3–5].Figure 1
General framework for uncertainty quantification.Having all the tools for the probabilistic simulations, still it is computationally expensive to perform a large set of nonlinear dynamic analyses on 3D model of dam-foundation-reservoir coupled system. Thus, it is important to look for methods which efficiently reduce the sample size. In the case of fragility curves, such a method is proposed by Azarbakht and Dolšek [7] which reduces the total number of required ground motions for an incremental dynamic analysis (IDA). In the case of material and modeling uncertainty, response surface method (RSM) is found to be a good solution [8]. The computational cost in RSM is significantly reduced compared to the crude Monte Carlo Simulation (MCS). Moreover, the explicit nature of response surface makes is useful for future assessments.The primary goal of this paper is to develop a response surface to calculate the probabilistic nonlinear seismic response of the arch dams, Figure1. Although the RSM is a well-established method, it has been applied to limited number of geostructures such as rock-fills [9], concrete-faced rock-fill dams [10], slope stability reliability analysis [11–13], and concrete gravity dams [14, 15].The main objective of this paper is to focus on the epistemic random variables (RVs). First, the concept of experimental design is studied in Section2, followed by response surface method, Section 3. Fundamentals of smeared crack in the mass concrete are reviewed in Section 4. Next, the case study and its finite element model are discussed in Section 5. Finally, the results of metamodeling are provided in Section 6 and validated.


## 2. Design of Experiment

A Design of Experiment (DOE) refers to a statistical procedure that systematically defines the efficient number of sampling data points to optimize the computed responses [16]. A DOE includes a “Factor” and “Level”. Factor is a parameter over which the designer or analyzer has direct control on an experiment. Moreover, level is the number of different values a factor can be assigned based on its discretization. In this section some DOE techniques are briefly presented.Randomized Complete Block Design (RCBD) is a procedure based on blocking and is used when the analyzer wants to focus on one or more particular factors. The number of required experiments for a RCBD is NDOE=∏i=1kLi, where k refers to factors and L represents levels.Latin Square Experimental Design (LSED) is similar to RCBD while it requires less samples. It performs single experiment in each block and requires some conditions that are needed to implement this technique. The sample size is NDOE=L2.Full Factorial Design (FFD) is a widely used technique in developing the metamodels and is usually constructed on a two-level or three-level basis. The samples are given by all the possible combination of the factors. Therefore, the sample size is NDOE=Lk.Fractional Factorial Design (FrFD) is a subset of FFD which reduces the total number of experiments/simulations. In this method, the sample size is usually one-half, one-third, or one-quarter of FFD. It is important to select the samples in a way that they be balanced. The sample size is NDOE=Lk-p, where p is related to the fraction of design.Star Experimental Design (SED) contains two axial points on the axis of each factor with a physical distance δ from the center. Also, it includes the center point with all the factors in their mean value. The sample size is NDOE=2k+1.Central Composite Design (CCD) is a combination of a two-level FFD and a SED. Therefore, the sample size is NDOE=k2+2k+1. This method is capable of estimating the curvature of the design space. Depending on location of star points (δ), CCD can take different forms:(i)
Central Composite Circumscribed (CCC): ifδ=1.(ii)
Central Composite Faced (CCF) and Inscribed (CCI): ifδ=1/k.(iii)
Central Composite Scaled (CCS): ifδ is arbitrarily selected to be either >1 or <1/k.Box-Behnken Experimental Design (BBED) is an incomplete three-level FFD. It is supposed to reduce the sample size as the number of parameters grows; however, it is still sufficient to estimate the coefficients of a 2nd-degree least squares polynomial. Specific tables are available to construct the BBED models [17].Plackett-Burman Experimental Design (PBED) includes a set of economical designs with the run number a multiple of four and they are appropriate for screening purposes. It is constrained with the condition that NDOE=k+1. Sample size is NDOE=k+4-MOD(k,4), in which MOD is modulo operation.Taguchi Experimental Design (TED) is based on distinction between the controllable and noise factors to reduce the sensitivity of the problem to the variations in uncontrollable factors. Sample size is NDOE=Lki+ko, in which the subscripts i and o refer to the inner and outer factors, respectively.Random Experimental Design (RED) relies on different techniques for filling uniformly the design space. RED is not based on the concept of levels and does not require discretization. Sample size, NDOE, is selected independently.Quasi Random Sequences (QRS) are generated from a completely deterministic, low-discrepancy process and possess no inherent statistical properties. Discrepancy is a metric for the degree of nonuniformity of numbers in a sequence [18]. Two main sequences are Halton and Sobol.Latin Hypercube Sampling (LHS) reduces the variance in the crude Monte Carlo Simulation (MCS) [19]. In LHS, first the given range, [0,1], is divided into N equal intervals 1/N. Then, a point is randomly selected from each interval.


## 3. Response Surface Method

Response surface method (RSM) has been used for a variety of problems in structural reliability and optimization [8, 20–27]. The fundamental idea is to use the results of a DOE to create an approximation of the response parameters. This approximation, which is called the response surface or “metamodel” (model of the model), can be constructed for different engineering demand parameters (EDPs).Since the response surface provides an analytical or explicit function, the further operations (e.g., reliability or optimization) on the system will be very fast and do not require extra experiments or simulations. Assuming that the response variable,y, is an unknown function of the input parameters, x, then the response surface y^ is an approximation of this function:(1)y=fx=f^x+ϵx⟹y^=f^xwhere ϵx is the error in the estimated response.The outcome of a DOE withNDOE experiments or simulations can be collected as NDOExi,yi couples in which any EDP yi is associated with a point xi in the design space. The response surface is said to beinterpolating if for each sample point yi=f^xi holds, orapproximating if ϵxi≠0 [16]. There are different techniques in order to approximate a response surface, e.g., least squares method, optimal RSM, shepard and K-Nearest, Kriging, Gaussian processes, radial basis functions, and artificial neural networks.In this paper, the least squares method (LSM) is used to construct the response surface. This method is originally developed by Gauss [28]. It is based on adjusting the coefficients in the response surface metamodel so that it best fits an observed data set (from experiments or simulations). Themodel function is defined as f^x,β, where β=β1,…,βmT is the vector of m unknown coefficients to be found and x=x1,…,xkT is the vector of k input parameters. The data set consists of xi,yi pairs, i=1,…,N, where xi is the input parameters of the ith simulation, whose EDP is yi. In LSM, the coefficients, βj, j=1,…,m, are estimated by minimizing the S function (i.e., the sum of squared residuals at the points in the data set):(2)S=∑i=1NDOEϵi2where the residuals are the difference between the actual responses and the predicted ones at the locations xi and can be written as ϵi=yi-f^xi,β,i=1,…,N. The minimum of S can be easily found by setting the gradient equal to zero:(3)∂S∂βj=2∑i=1NDOEϵi∂ϵi∂βj=-2∑i=1NDOEyi-f^xi,β∂f^xi,β∂βj=0,j=1,…,mLeast squares problems are divided into two groups: linear and nonlinear. The following provides insight into each solution technique.Linear problems have a closed-form (analytical) solution; however, they are not quit accurate and they only provide the general trends of the EDPs over the design space. For a problem with NDOE simulations/experiments and k parameters, the metamodel function takes the following form:(4)f^x,β=β0+β1x1+⋯+βkxkf^xi,β=β0+∑j=1kxi,jβjor the form of matrix notation,(5)y=Xβ+ϵ,where(6)y=y1y2⋮yN,X=1x1,1⋯x1,k1x2,1⋯x2,k⋮⋮⋱⋮1xN,1⋯xN,k,β=β1β2⋮βN,ϵ=ϵ1ϵ2⋮ϵNSubsequently, (2) can be written as(7)S=ϵTϵ=y-XβTy-Xβ=yTy-2βTXTy+βXTXβDeriving (7), equating to zero, and solving in β yield(8)∂S∂β=-2XTy+2XTXβ=0⟹β=XTX-1XTyand the EDP of the estimated (fitted) metamodel is(9)y^=XβNonlinear problems should be solved iteratively. In this method, first the initial values for the coefficients, β(1), are chosen. Then, it is updated iteratively (it is known as Gauss-Newton algorithm):(10)βk+1=βk+Δβkwhere Δβ(k) refers to the shift vector. This vector can be updated using an iterative model by approximation to a first-order Taylor series expansion about β(k)(11)f^xi,βk+1=f^xi,βk+∑j=1m∂f^xi,βk∂βjβjk+1-βjk︸Ji,jkΔβjkwhere J is a N×m Jacobian matrix of f^ with respect to β. This equation can be written in the form of matrix notation(12)y=y^k+ϵk=y^k+1+ϵk+1=y^k+JkΔβk+ϵk+1Subsequently, derivative ofS with respect to β takes the following form:(13)∂S∂β=-2JkTϵk+2JkTJkΔβkSolving inΔβ(k) yields (14)Δβk=JkTJk-1JkTϵkBoth linear and nonlinear problems can be used in the context of the complete or incomplete polynomials. The number of required coefficients,m, for a nth-degree polynomial, with k variables, can be calculated as(15)m=k+nn=k+n!n!k!and, subsequently, the general expression for the nth-degree polynomial can be written as(16)f^nthx,β=f^n-1th+nth  order  termsMajority of the real-world civil engineering problems can be estimated using one of the following expressions for the linear and quadratic forms:(17)f^1stx,β=β0+∑i=1kβixif^2ndx,β=f^1st+∑i=1k∑j=1iβi,jxixjThe goodness-of-fit (GOF) in an approximation can be estimated by regression parameters which varies between 0.0 and 1.0. The higher the GOF is, the better the model is expected to be. A widely used GOF is(18)R2=1-∑i=1Nyi-y^i2∑i=1Nyi-y¯2;y¯=1N∑i=1Nyi


## 4. Smeared Crack Model

In concrete dam engineering, the nonlinearity is originated mainly from two sources: continuum crack model and discrete crack model. The latter one is usually used when the location and direction of a potential crack or joint are already known. This model has been successfully used for modeling the contraction joints in arch dams [6]. The former one itself can be divided into two major groups, i.e., the damage mechanics approach [29–31] and the smeared crack approach [32–34].The nonlinear response of the case study dam in this paper is modeled by smeared crack approach. Thus, this section briefly reviewed the main formulation of this technique. Smeared cracks are convenient when the crack orientations are not known beforehand. It does not require a remeshing or new degrees of freedom. A coaxial rotating model is used to simulate the concrete crack under the dynamic loading. In this model, the precracked constitutive relation is replaced by the cracked one where the reference axis is aligned with the fracture direction. The main assumptions are as follows: (1) the concrete is initially linear isotropic until it reaches the ultimate strength, (2) the concrete modulus of elasticity is taken as the average instead of the linear actual one, and (3) during the softening phase, an anisotropic modulus matrix is considered for material.Figure2 shows the smeared crack-based precracked and cracked constitutive relationships. Cracking occurs when the principal stresses, in any direction, lie outside the failure surface. In this model, cracking is permitted at each Gaussian point and in three orthogonal directions.Figure 2
Constitutive model for the concrete based on coaxial rotating smeared crack model; adopted from [6].Next, a failure criterion (e.g., yielding, load carrying capacity, and initiation of cracking) is required to quantify the ultimate strength surface. In the present paper, Willam and Warnke [35] failure criterion is used to identify the initiation and propagation of the cracks in mass concrete. Since this formulation is not the main scope of this paper, only a brief formulation is illustrated in Figure 3.Figure 3
Formulation of failure criterion for concrete; adopted from [6].


## 5. Numerical Example

An ultra-high arch dam is selected as a vehicle for the numerical examples [6]. The finite element model of the dam, the foundation, and the reservoir are prepared in ANSYS [36] software, Figure 4. The model consists of 792 and 3,770 solid brick elements for dam and foundation, and 3,660 fluid brick elements for the reservoir domain. All the nodes on the far-end boundary of the foundation are restricted in three transitional directions. Nonlinearity stems only from the concrete cracking and the contraction joints are modeled in this study. This is to keep the number of overall variables in a reasonable range to perform the response surface metamodeling.Figure 4
Finite element model of dam, reservoir, and foundation system.Table1 lists the most important material properties that were used to develop the finite element model. Four material properties are assumed to be random variables (RVs). Three of them belong to concrete, i.e., elasticity, mass density, and compressive strength, and the remaining one belongs to the foundation, i.e., modulus of elasticity. These parameters are among those most sensitive RVs identified for concrete dams in [37]. Furthermore, it is noteworthy to say that a full correlation is assumed among the concrete tensile and compressive strength through a deterministic relationship ft′=0.1fc′.Table 1
Material properties for the finite element model.
Quantity Symbol Unit DOE UQ Mean (0) Variation Lower (-1) Upper (+1) Model COV Truncation Concrete modulus of elasticity E c MPa 40,000 15% 34,000 46,000 Normal 0.08 [ 34 , 000   46 , 000 ] Concrete Poisson’s ratio ν c - 0.2 - - - - - - Concrete mass density ρ c kg/m3 2,400 10% 2,160 2,640 Normal 0.06 [ 2 , 160   2 , 600 ] Concrete compressive strength f c ′ MPa 35.0 20% 28.0 42.0 Normal 0.10 [ 25.0   40.0 ] Concrete tensile strength f f ′ MPa 3.5 20% 2.8 4.2 Normal 0.10 [ 2.5   4.0 ]  Foundation modulus of elasticity E f MPa 20,000 25% 15,000 25,000 Normal 0.10 [ 15 , 000   25 , 000 ] Foundation Poisson’s ratio ν f - 0.25 - - - - - - Foundation mass density ρ f kg/m3 2,500 - - - - - -  Water mass density ρ w kg/m3 1,000 - - - - - - Wave velocity in water C w m/s 1,440 - - - - - - Wave reflection coefficient at the bottom α w - 0.80 - - - - - - Sediment/silt mass density ρ s kg/m3 1,360 - - - - - - DOE: Design of Experiment; UQ: Uncertainty Quantification; COV: Coefficient of Variation.Having four RVs, many of the DOE techniques in Section2 can be adopted. For example, two-level and three-level FFD require 16 and 81 simulations, respectively. A SED needs 9 and a CCD requires 25 simulations. However, in the present study, the BBED technique is used. The main reason can be attributed to the fact that it is an incomplete three-level FFD. Thus, it can capture the potential curvature of the response surface, while it requires less simulations compared to FFD. Box and Behnken [17] determine the required number of simulations as(i)
k=3, NDOE=12+1=13(ii)
k=4, NDOE=24+1=25(iii)
k=5, NDOE=40+1=41(iv)
k=6, NDOE=48+1=49In two-level designs (e.g., FFD and SED), lower and upper bounds are assumed for experimental design. However, in three-level designs, another level is assumed in between which is referred to mean. It is convenient to use the “coded” (standardized) factors,xicoded, instead of the actual values, xiactual. In the case of two-level designs, they are referred to “±1”, and in the case of three-level designs, they are “-1”, “0” and “+1”. Any desired value in the range of (min, max) of the actual values can be converted to the coded one:(19)xicoded=xiactual-xiactualU+xiactualL/2xiactualU-xiactualL/2where the superscripts U and L present the upper and lower bounds of the actual values, respectively.It is possible to graphically show the BBED with three factors (as a 3D cube), Figure5, where the samples are at the center and middle of the edges. Such a plot is not possible for k=4 in the Cartesian coordinate system; however, the simulations can be presented in the form of Table 2. In this table, x1 to x4 correspond to Ec, ρc, fc′, and Ef, respectively. The lower bound, mean, and the upper bound for each RV can be found in Table 1.Table 2
Tabulated simulations in a Box-Behnken experimental design withk=4.
Simulations Coded values x 1 x 2 x 3 x 4 1-4 ± 1 ± 1 0 0 5-8 ± 1 0 ± 1 0 9-12 ± 1 0 0 ± 1 13-16 0 ± 1 ± 1 0 17-20 0 ± 1 0 ± 1 21-24 0 0 ± 1 ± 1 25 0 0 0 0Figure 5
Graphical representation of Box-Behnken experimental design fork=3.


## 6. Results

The coupled system is excited using the ground motion shown in Figure6(a). This ground motion is scaled to maximum design level (MDL) for the dam site Hariri-Ardebili and Kianoush [6]. Displacement response of the pilot dam subjected to this signal is also shown in Figure 6(b). Pilot dam is referred to the one with all the material properties at their mean value, Table 1. The pilot model is also called Sim-25 in Table 2. Note that there is an initial displacement towards downstream due to hydrostatic pressure. Figure 7 illustrates the cracking and crushing of pilot dams. As discussed already, cracking is allowed in three orthogonal directions and they are presented as first, second, and third cracks. Obviously, the damage distribution is higher in Figure 7(a) than Figures 7(b) and 7(c). Furthermore, there is no concrete crushing for this pilot model, Figure 7(d). Note that color variation presents the degree of element damage.Figure 6
Structural analysis of pilot model.
(a)
Applied acceleration (b)
Upstream/downstream crest displacementFigure 7
Damage distribution in upstream and downstream faces of pilot model.
(a)
First crack (b)
Second crack (c)
Third crack (d)
CrushNext, the BBED analyses are performed as illustrated in Table3. Note that ft′ always is assumed to be 10% of the compressive strength. The displacement response is considered in two levels: (1) preseismic and (2) maximum seismic. The major observations based on static displacement are as follows:(i)
In general, the ratio of maximum dynamic displacement to the static one is 5-7 for all the 25 simulations.(ii)
There is no cracking in the dam under the static loads.(iii)
Based on the results,fc′ (and ft′) has no effect on the static displacement.(iv)
Per this table, increasingEc, ρc, and Ef reduces the displacement. However, sensitivity of Ec is higher than the others, and ρc has the minimum impact.Table 3
Summary of the BBED simulations.
Sim E c ρ c f c ′ E f Δ m a x s t a [mm] Δ m a x d y n [mm] DIvolcr1[%] DIvolcr2[%] DIvolcr3[%] DIvolcrush[%] 1 46 2640 35 20 19.7 117.6 17.6 6.3 3.2 0.1 1 46 2160 35 20 19.8 95.1 11.5 2.1 0.7 0.0 3 34 2640 35 20 25.4 168.7 14.2 5.8 3.7 0.3 4 34 2160 35 20 25.4 137.6 11.6 3.0 1.2 0.0 5 46 2400 42 20 19.8 99.5 8.4 1.0 0.3 0.0 6 46 2400 28 20 19.8 99.5 8.4 1.0 0.3 0.0 7 34 2400 42 20 25.4 141.6 8.9 1.3 0.3 0.0 8 34 2400 28 20 25.4 177.0 16.8 7.8 5.6 1.0 9 46 2400 35 25 19.0 98.2 13.5 3.3 1.3 0.0 10 46 2400 35 15 21.1 113.7 19.8 4.3 1.7 0.1 11 34 2400 35 25 24.7 146.8 12.3 4.4 1.9 0.0 12 34 2400 35 15 26.7 159.5 15.0 5.2 2.1 0.1 13 40 2640 42 20 22.1 130.2 0.8 0.0 0.0 0.0 14 40 2640 28 20 22.1 149.0 22.4 9.5 6.5 1.9 15 40 2160 42 20 22.2 108.8 5.5 0.3 0.1 0.0 16 40 2160 28 20 22.2 151.0 18.7 9.2 7.0 2.4 17 40 2640 35 25 21.4 133.0 14.4 5.8 2.7 0.0 18 40 2640 35 15 23.3 151.0 19.1 6.3 3.2 0.2 19 40 2160 35 25 21.4 109.1 10.2 1.9 0.8 0.0 20 40 2160 35 15 23.5 122.4 14.1 3.1 1.0 0.1 21 40 2400 42 25 21.4 114.8 7.8 1.3 0.4 0.0 22 40 2400 42 15 23.4 127.0 10.9 1.4 0.2 0.0 23 40 2400 28 25 21.4 145.0 17.4 8.4 5.9 1.3 24 40 2400 28 15 23.4 168.0 24.4 8.8 6.1 1.7 25 40 2400 35 20 22.2 126.3 14.7 4.5 1.9 0.0On the other hand, 5 quantities are reported for the seismic simulations, i.e., maximum displacement, cracking in three directions, and crushing. Except the displacement, the others are presented in terms of the volumetric damage index, DI:(20)DIvol=Damaged  volumeTotal  volumeNote that all the solid elements in the present formulation have eight Gaussian points in which the damage may occur in these points. Subsequently, status of each element can be presented with a number as 0, 1/8,  ..., 7/8, 1, in which zero means no damage in the element and one shows the full damage. It is worth mentioning thatDIvolcr3⊂DIvolcr2⊂DIvolcr1. Furthermore, Figures 8, 9, 10, and 11 show the damage distribution results from first and third cracks on upstream (US) and downstream (DS) faces of dam, respectively. The major observations are the following:(i)
Dynamic displacement varies from 95 to about 177 mm in different models.(ii)
DIvolcr1 varies from 0.8 to 24.4%, DIvolcr2 varies from 0 to 9.5%, DIvolcr3 varies from 0 to 7.0%, and DIvolcrush is limited to 2.4%.(iii)
As seen, the main cracking is distributed in central upper parts of the dam in vicinity of crest. In some cases, there are limited cracking at the lower parts, near the dam-foundation interface.(iv)
Cracking in the third direction is more concentrated on the downstream face, which shows its vulnerability.Figure 8
Distribution of first crack on upstream face; top left: Sim-1; first row: 1 to 6; bottom right: Sim-24.Figure 9
Distribution of first crack on downstream face; top left: Sim-1; first row: 1 to 6; bottom right: Sim-24.Figure 10
Distribution of third crack on upstream face; top left: Sim-1; first row: 1 to 6; bottom right: Sim-24.Figure 11
Distribution of third crack on downstream face; top left: Sim-1; first row: 1 to 6; bottom right: Sim-24.In order to quantify the sensitivity of all material trampers, and also derive a response surface metamodel, a second-order polynomial including the interaction terms is used:(21)EDP^=β1+β2Ec+β3ρc+β4fc′+β5Ef+β6Ec2+β7ρc2+β8fc′2+β9Ef2+β10Ecρc+β11Ecfc′+β12EcEf+β13ρcfc′+β14ρcEf+β15fc′Efwhere the hat sign represents the estimated EDP.Five EDPs are used in (21), i.e., two displacements, and three cracking-based DIs. The predicted βi values are shown in Figure 12. To make the comparison easy, the absolute logarithmic values are shown. Large negative values for log⁡βi present the lower contribution of that parameter in overall response surface metamodel. For example, based on Figure 12(a), contribution of fc′ (i.e., β4) is negligible in static analyses; however, it plays an important role in dynamic part (Figure 12(b)).Figure 12
Validation of the predicted response surface using LHS method.
(a)
Static displacement (b)
Dynamic displacement (c)
First crack (d)
Second crack (e)
Third crackIn addition, the probability of exceedance of response parameters is shown in the last column of Figure12. In each plot, the empirical cumulative distribution function (CDF) of the FEM and BBED is plotted, as well as the fitted curves using the log-normal (LN) distributional model. LN is the best model to fit the data points among others (e.g., normal, beta). In all cases, these are a good match between the empirical and fitted curves. For all five responses, the BBED model fluctuates around the FEM, and two fitted curves cross at the probability of 0.50-0.65.The estimated response surfaces should be validated using an independent data set. For this purpose, 250 samples for each of four RVs are computed using the LHS technique [38]. No correlation is assumed among the RVs. Mean, COV, and the upper and lower bounds are already provided in Table 1. The resulting samples are shown in Figure 13. Subsequently, 250 extra nonlinear transient analyses are performed and the EDPs are computed. Real versus estimated values are also shown in Figure 12. As seen the individual data points are close to equity line. The estimated coefficient of determination, R2, is 0.96, 0.92, 0.91, 0.90, and 0.88 for Δmaxsta, Δmaxdyn, DIvolcr1[%], DIvolcr2[%], and DIvolcr3[%], respectively. Smaller R2 for the third crack can be attributed to this fact that there is a smaller variation among the percent cracking of the models. Moreover, the root mean square error, RMSE=1/N∑i=1N‍EDPi-EDP^i2, for the five EDPs is 0.24, 3.89, 0.65, 0.43, and 40.Figure 13
LHS-based sampling for validation of metamodels.


## 7. Summary

This paper presented the results of probabilistic analysis of an arch dam with uncertainty in the material properties. A detailed finite element model of the dam-foundation-reservoir coupled system is developed. The fluid-structure interaction is modeled by Eulerian approach. The nonlinearity in the system is originated from the smeared crack modeling in the mass concrete. The coupled system is excited using an earthquake record in the maximum design level of the dam site.Over ten different designs of experiment techniques are critically reviewed and the Box-Behnken experimental design (BBED) is chosen to be used in this paper. Advantage of BBED over the other techniques is that it is an incomplete three-level full factorial design and thus it can capture the nonlinearity of the design space, while it requires small number of simulations.For an engineering problem with four random variables, BBED requires only 25 simulations. The outputs of these analyses are reported in terms of the preseismic and seismic displacements, as well as the volumetric damage index in three orthogonal cracking directions. Then, a response surface metamodel is developed for each of the engineering demand parameters. The metamodels are provided in the explicit form based on the quadratic polynomial including the interaction terms.Next, it is important to validate the applicability of the developed analytical models through another set of independent simulations. Since nonlinear dynamic analysis of 3D arch dams with foundation and reservoir is computationally expensive, Latin Hypercube Sampling is used with 250 simulations. No correlation is assumed among the generated random numbers.Results show a good consistency between the finite element based outputs and those estimated by response surface metamodel. Coefficient of determination is about 0.9 in all cases. Furthermore, a discussion is provided about the cracking pattern of arch dams under the material uncertainty.Results of this study are important since they provide a framework for uncertainty quantification of the material properties in the structural and infrastructural system subjected to varying environmental conditions. There are a large number of dams, bridges, and nuclear power plants which are aging and they need to be reevaluated regularly based on the updated (or degraded) material properties. Once such metamodels are developed, they can be used without any extra cost to evaluate the reliability of the system under the new condition. This saves a considerable amount of money and time and facilitates the risk-based decision-making.



---

*Source: 1784203-2018-07-05.xml*