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.
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, Figure 1. 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].
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, Figure 1. 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 Section 2, 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 , where refers to factors and 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 .
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 .
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 , where 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 .
Central Composite Design (CCD) is a combination of a two-level FFD and a SED. Therefore, the sample size is . 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 .(ii)Central Composite Faced (CCF) and Inscribed (CCI): if .(iii)Central Composite Scaled (CCS): if is arbitrarily selected to be either or .
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 -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 . Sample size is , 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 , in which the subscripts and 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, , 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, , is divided into equal intervals . 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, , is an unknown function of the input parameters, , then the response surface is an approximation of this function:where is the error in the estimated response.
The outcome of a DOE with experiments or simulations can be collected as couples in which any EDP is associated with a point in the design space. The response surface is said to be interpolating if for each sample point holds, or approximating if [16]. There are different techniques in order to approximate a response surface, e.g., least squares method, optimal RSM, shepard and -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). The model function is defined as , where is the vector of unknown coefficients to be found and is the vector of input parameters. The data set consists of pairs, , where is the input parameters of the th simulation, whose EDP is . In LSM, the coefficients, , , are estimated by minimizing the function (i.e., the sum of squared residuals at the points in the data set):where the residuals are the difference between the actual responses and the predicted ones at the locations and can be written as . The minimum of can be easily found by setting the gradient equal to zero:
Least 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 simulations/experiments and parameters, the metamodel function takes the following form:or the form of matrix notation,where
Subsequently, (2) can be written as
Deriving (7), equating to zero, and solving in yieldand the EDP of the estimated (fitted) metamodel is
Nonlinear problems should be solved iteratively. In this method, first the initial values for the coefficients, , are chosen. Then, it is updated iteratively (it is known as Gauss-Newton algorithm):where refers to the shift vector. This vector can be updated using an iterative model by approximation to a first-order Taylor series expansion about where is a Jacobian matrix of with respect to . This equation can be written in the form of matrix notation
Subsequently, derivative of with respect to takes the following form:
Solving in yields
Both linear and nonlinear problems can be used in the context of the complete or incomplete polynomials. The number of required coefficients, , for a th-degree polynomial, with variables, can be calculated asand, subsequently, the general expression for the th-degree polynomial can be written as
Majority of the real-world civil engineering problems can be estimated using one of the following expressions for the linear and quadratic forms:
The 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
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.
Figure 2 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.
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.
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.
Table 1 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 .
Having four RVs, many of the DOE techniques in Section 2 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), (ii), (iii), (iv),
In 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, , instead of the actual values, . In the case of two-level designs, they are referred to “”, and in the case of three-level designs, they are “”, “0” and “+1”. Any desired value in the range of (min, max) of the actual values can be converted to the coded one:where the superscripts and 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), Figure 5, where the samples are at the center and middle of the edges. Such a plot is not possible for in the Cartesian coordinate system; however, the simulations can be presented in the form of Table 2. In this table, to correspond to , , , and , respectively. The lower bound, mean, and the upper bound for each RV can be found in Table 1.
6. Results
The coupled system is excited using the ground motion shown in Figure 6(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.
(a) Applied acceleration
(b) Upstream/downstream crest displacement
(a) First crack
(b) Second crack
(c) Third crack
(d) Crush
Next, the BBED analyses are performed as illustrated in Table 3. Note that 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, (and ) has no effect on the static displacement.(iv)Per this table, increasing , , and reduces the displacement. However, sensitivity of is higher than the others, and has the minimum impact.
On 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:
Note 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 that . 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) varies from 0.8 to 24.4%, varies from 0 to 9.5%, varies from 0 to 7.0%, and 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.
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: where the hat sign represents the estimated EDP.
Five EDPs are used in (21), i.e., two displacements, and three cracking-based DIs. The predicted values are shown in Figure 12. To make the comparison easy, the absolute logarithmic values are shown. Large negative values for present the lower contribution of that parameter in overall response surface metamodel. For example, based on Figure 12(a), contribution of (i.e., ) is negligible in static analyses; however, it plays an important role in dynamic part (Figure 12(b)).
(a) Static displacement
(b) Dynamic displacement
(c) First crack
(d) Second crack
(e) Third crack
In addition, the probability of exceedance of response parameters is shown in the last column of Figure 12. 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, , is 0.96, 0.92, 0.91, 0.90, and 0.88 for , , , , and , respectively. Smaller 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, , for the five EDPs is 0.24, 3.89, 0.65, 0.43, and 40.
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.
Data Availability
All the raw data used to support the findings of this study (resulting from numerical simulations) are included within the article. However, the finite element models are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The first author would like to express his sincere appreciation to his former advisor (and the current PI), Professor Victor E. Saouma, at the University of Colorado Boulder for his enthusiastic guidance and advice throughout this research.