A response surface method for stochastic dynamic ... - Semantic Scholar

Reliability Engineering and System Safety 126 (2014) 44–53

Contents lists available at ScienceDirect

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

A response surface method for stochastic dynamic analysis Umberto Alibrandi n National University of Singapore, Civil and Environmental Engineering, 1 Engineering Drive 2, Singapore 117576, Singapore

art ic l e i nf o

a b s t r a c t

Article history: Received 4 June 2012 Received in revised form 20 September 2013 Accepted 5 January 2014 Available online 17 January 2014

A Response Surface (RS) strategy is presented for the evaluation of the response statistics of dynamic systems subjected to stochastic excitation. The proposed approach adopts a strategy based on the High Dimensional Model Representation (HDMR), which gives a Gaussian Model (GM) of the response. The GM requires only a reduced number of analyses which can be adopted for all the degrees of freedom of a MDOF dynamic system and it can be successfully adopted for weakly nonlinear dynamic systems. For more strongly nonlinear systems a Non-Gaussian approximation may be necessary for the highest response thresholds. In this paper this issue is accomplished through the FORM solution, and the design point is obtained by using a response surface method recently proposed by the author and Der Kiureghian to this aim. The latter response surface is based on a variant of the Model Correction Factor Method (MCFM), which is here applied by using as a starting model the GM itself. In many applications of engineering interest, both the input and the response processes are stationary, so that the stochastic excitation through the Fourier series can be modeled in terms of the underlying Power Spectral Density (PSD). In these cases, it is seen that the dynamic computations required by the proposed approach can decrease significantly. The application to SDOF and MDOF hysteretic systems shows the effectiveness of the presented method. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Stochastic dynamics FORM Response surface Linearization Bouc–Wen

1. Introduction The aim of the stochastic dynamic analysis is the evaluation of the response of dynamic systems subjected to stochastic input. If the stochastic excitation follows a Gaussian distribution and the dynamic system is linear, then the response follows a Gaussian distribution, too, and it is relatively easy to be determined. Conversely, the response of a nonlinear dynamic system follows a Non-Gaussian distribution, and its determination is a very complicated task to be accomplished. In the last decades a lot of research has been devoted to this topic, however most approaches are not easily applicable to general non-linear MDOF systems, and so they are difficult to apply in practice. These drawbacks are not shared by the Equivalent Linearization Method (ELM) [1], which has gained wide popularity because of its versatility and applicability, see among others, [2,3]. The basic idea is to replace the original non-linear system by an equivalent linear one, whose determination is performed by minimizing the difference between the two systems in some statistical sense. The ELM exhibits different forms based on the adopted probability density function for the evaluation of the coefficients that appear in the

n

Tel.: þ 65 6516 6498; mob.: þ 65 8670 5478. E-mail addresses: [email protected], [email protected]

0951-8320/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ress.2014.01.003

linearized system. Usually the Gaussian distribution is adopted, which allows to approximate the second order moments of the response. Against relatively little numerical efforts, unfortunately it gives accurate results for weakly non-linear systems only [1,4]. Moreover, the ELM generally cannot approximate adequately the distribution probability of the response, especially in the tail region. Therefore some response statistics as crossing rates and first passage probability will be inaccurate at high thresholds for systems strongly nonlinear. The most robust procedure is given by the Monte Carlo Simulation (MCS), which is however strongly demanding in its crude form. For this reason, recently some smart simulation techniques have been proposed, among the other we recall the subset simulation [5], line sampling [6], asymptotic sampling [7]. Promising results are given from the application to the nonlinear stochastic dynamic analysis of the analytical methods of structural reliability, particularly the first-order reliability method (FORM) [8,9]. At first, the stochastic input is discretized into a large number of standard normal random variables. The tail probability is defined as the probability that the response of the dynamic system is greater than a chosen threshold x at fixed time instant t. In this way, for given x and t the dynamic problem may be solved by using FORM. The knowledge of the design point is of great importance, because: (i) it corresponds to the most likely realization of the stochastic input that gives rise to the tail exceedance event, and therefore it defines a critical excitation for the system,

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

(ii) it gives the FORM solution, which has been shown to give very good approximations of the tail probability in many cases of practical interest [10–12], (iii) it allows the application of the recently developed tail equivalent linearization method (TELM) [13–17]. However, in high-dimensional spaces, as the one here analyzed, the evaluation of the design point is a challenging task. Indeed, the design point is obtained as a solution of a constrained optimization problem [18,19]. For some material models, the problem is not numerically very smooth, and the finite element code fails to produce a result due to lack of numerical convergence, so that some suitable remedying strategies have to be adopted [20]. Moreover, for the solution of the optimization problem gradient-based procedures are usually adopted, which require repeated evaluation of the response gradient. Typically, the finite elements codes do not provide these gradients, which can be approximated by the finite differences method (FDM). By using FDM, each gradient at each iteration of the iterative procedure requires n þ1 nonlinear dynamic computations, n being the number of random variables, and consequently the cost of computation may be excessive. Moreover the selection of the perturbation parameter in FDM is questionable, and in particular in high-dimensional spaces the accuracy of the response gradient may be lost. Consistency, accuracy and efficiency of the gradients can be achieved by using the direct differentiation method (DDM) [21–23], which involves analytical differentiation of the discretized response equations. However, the DDM requires alterations at the finite element (FE) code level, and it is necessary to use a DDM-enabled software, such as OpenSees [24]. To reduce the computational cost and using an existing FE code as a “black-box”, an alternative strategy is given from the Response Surface Methodology (RSM), which builds a surrogate model of the target Limit State Function (LSF), defined in a simple and explicit mathematical form [25–28]. Once the Response Surface (RS) is built, it is possible to substitute the RS for the target LSF, and then it is no longer necessary to run demanding finite element analyses. In this paper we used a novel response surface strategy based on the High Dimensional Model Representation (HDMR) [29]. The HDMR is a set of analysis tools for capturing the high-dimensional relationships between sets of input and output model variables. The existing HDMR-based response surfaces [30] require for the first-order representation of the limit state function on average 5–7 points for each random variable, so for the problem at hand a huge computational effort arises, since we have a different limit state function for each response threshold. In this paper we present a suitable chosen response surface requiring only 2 points for each random variable, reducing significantly the dynamic computations. The proposed application of the HDMR gives a Gaussian Model (GM) for the system response and it gives quite good approximations for weakly nonlinear dynamic systems. For systems with strong nonlinearities it is likely that especially for the highest thresholds the FORM solution outperforms the GM. In any case, the design point of the GM is likely to be close to the design point of the original problem, so it can represent a good starting point for FORM. It is expected that in general the algorithm converges with only a few iterations. In this paper, to reduce further the computational cost, we adopted the “improved Model Correction Factor Method” (iMCFM), developed by the author with A. Der Kiureghian in recent papers [11,12]. The paper is organized as follows: in Sections 2 and 3 we present the nonlinear stochastic dynamic problem by using the proposed response surface strategy giving rise to the GM, in Section 4 its improvement through the MCFM is discussed, and

45

finally in Section 5 the method is applied to hysteretic systems, showing its accuracy and efficiency.

2. Discretization of the stochastic input The proposed approach requires the preliminary discretization of the stochastic input into a set of standard normal random variables. Several formulations for this purpose are available, see Der Kiureghian [9] for a brief review. For a zero-mean, Gaussian excitation process, all representations have the form nðτÞ

Fðt; uÞ ¼ ∑ sði τÞ ðtÞui ¼ sðτÞ ðtÞ  u

ð1Þ

i¼1

where n ¼ nðτÞ is a measure of the resolution of the discretization, u ¼ fu1 ; u2 ; …; un gT is an n-vector of standard normal random ðτÞ variables, sðτÞ ðtÞ ¼ f s1 ðtÞ

sð2τÞ ðtÞ

…

sðnτÞ ðtÞ g

T

is an n-vector of

deterministic shape functions dependent on the covariance function of the process. In (1) the superscript “τ” refers to the discretization of the input in the time domain [8–14,16,17]. In the special case of sði τÞ ðtÞ ¼ s  Δt  δðt  t i Þ, where δðtÞ is the Dirac delta function and t i ¼ iΔt, i ¼ 1; 2; …; n, are a set of equally spaced time points, then, Fðt; uÞ represents a discretized band-limited white noise excitation of intensity SW , of variance s2F ¼ ð2π SW Þ=Δt and with cut-off frequency ωc ¼ π =Δt. If the excitation is a stationary process defined by the power spectral density (PSD) function SF ðωÞ, through the discrete Fourier series the stochastic input can be expressed as [15–17] nðωÞ

nðωÞ

k¼1

k¼1

Fðt; uÞ ¼ ∑ F k ðt; uÞ ¼ ∑

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2SF ðωk ÞΔω½ cos ðωk tÞu0k þ sin ðωk tÞu″k 

nðωÞ

ωÞ ðtÞu″k ¼ sðωÞ ðtÞ  u ¼ ∑ s0ðk ωÞ ðtÞu0k þ s″ð k

ð2Þ

k¼1

where n ¼ 2nðωÞ is a measure of the resolution of the discretization, u ¼ fu01 ; u″1 ; u02 ; u″2 ; …; u0n=2 ; u″n=2 gT is an n-vector of standard normal random variables, sðωÞ ðtÞ ¼ f s1 ðωÞ 0

s″ð1ωÞ

…

s0n=2 ðωÞ

ωÞ s″ðn=2 g

T

is an n-vector of deterministic shape functions dependent on the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PSD of the process, being s0k ðωÞ ¼ 2SF ðωk ÞΔω cos ðωk tÞ and ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p s″ðkωÞ ¼ 2SF ðωk ÞΔω sin ðωk tÞ; here the superscript “ω” refers to the discretization of the input in the frequency domain. To obtain a good approximation of Fðt; uÞ the frequency step is Δω r ð2π Þ=t, being t the time instant when the input process reaches stationarity; if t is equal to the last time instant of (1) and choosing Δω ¼ ð2π Þ=t so that nðωÞ ¼ nðτÞ =2, then the (2) corresponds to (1), in the sense that the PSD SF ðωÞ is sampled at the frequencies 0; ω1 ; ω2 ; …; ωn=2 that we would obtain by applying the discrete Fourier transform (DFT) to each sample of (1). Under these conditions the variance of a band-limited white noise excitation of intensity SW by using (2) is s2F ¼ ð2π SW Þ=Δt.

3. An alternative Gaussian approximation of the system response Consider the response of a dynamical system to the excitation in (1)or (2). Owing to the random variables u, the response is stochastic and we denote it as Xðt; uÞ. For a specified threshold x and time t, we define the tail probability as Prob½Xðt; uÞ Z x. To apply the tools of the structural reliability theory, we define the limit state function (LSF) gðx; t; uÞ ¼ x  Xðt; uÞ

ð3Þ

46

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

so that the failure probability with respect to the limit state P f ¼ Prob½gðx; t; uÞ r0 is equal to the tail probability. To simplify the computational effort required by the evaluation of the tail probability, a response surface approach is here adopted, based on a simple application of the High Dimensional Model Representation (HDMR) to the stochastic dynamic analysis.

retaining only the terms up to the first-order, and choosing as a reference point the origin of the standard normal space n

n

i¼1

i¼1

Xðt; uÞ ¼ X 0 ðtÞ þ ∑ X i ðt; ui Þ þ ⋯ ffi ∑ X i ðt; ui Þ

ð7Þ

where taking account of (5) and (6) X 0 ðtÞ ¼ Xðt; 0; 0; …; 0Þ ¼ 0;

3.1. The High‐Dimensional Model Representation (HDMR)

X i ðt; ui Þ ¼ Xðt; 0; …; 0; ui ; 0; …; 0Þ

The HDMR is a general set of quantitative tools for capturing the relationships between the output function f ðxÞ and the input T variables x ¼ f x1 x2 … xn g . The HDMR expansion reads as

where

It is noted that the first-order components X i ðt; ui Þ give the response of the system at the time instant t when only ui a 0, while uk ¼ 0, with k a i. The X i ðt; ui Þ have different physical meaning if we consider the excitation in terms of random pulses (1) or harmonic components (2). In the first case, X i ðt; ui Þ is the response of the system at the time instant t when only the random pulse at the time instant t i acts. It is noted that on a large time scale and for systems containing damping, the most important terms of the series expansion (7) correspond to the random pulses acting in time instants closest to t. This implies that it is possible consider a reduced number n^ o n of components X i ðt; ui Þ without affecting significantly the accuracy of the solution. Consider now the stationary case, where we represent the stochastic excitation in terms of Gaussian harmonic components. In this case X i ðt; ui Þ is the response of the system at the time instant t pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi when only the harmonic component 2S ð ωi ÞΔω cos ðωi tÞ or F pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2SF ðωi ÞΔω sin ðωi tÞ acts. The most important terms of the series expansion (7) correspond to the lowest frequencies and closest to the ðωÞ first modal frequencies of the system. The choice of n^ ¼ 2n^ o2nðωÞ components implies that we are neglecting the content of the highest frequencies. It is clear that using the representation (2) we can easily detect the range of the frequencies of interest. Note moreover that for the symmetry of the problem, X i ðt; ui Þ is an odd function of ui and X i ðt; ui Þ ¼ X i ðt;  ui Þ so that it is enough to consider only points ui 4 0. At first, consider the linear case. The first-order components X i ðt; ui Þ are linear in ui , and they are exactly interpolated by the knowledge of only one point

ðxi ; xi0 Þ ¼ ðx0;1 ; :::; x0;i  1 ; xi ; x0;i þ 1 ; :::; x0;n Þ

X^ i ðt; ui Þ ¼ ci ðtÞui

n

f ðxÞ ¼ f 0 þ ∑ f i ðxi Þ þ i¼1

∑

1riojrn

f ij ðxi ; xj Þ þ ⋯ þ f 12…n ðx1 ; x2 ; …xn Þ

ð4Þ

where f 0 is a constant term representing the zeroth-order component function or the mean response of f ðxÞ, f i ðxi Þ is a first-order term expressing the effect of the variable xi acting alone, although nonlinearly, upon the output f ðxÞ, f ij ðxi ; xj Þ is a second-order term which describes the cooperative effects of the variables xi and xj upon the output f ðxÞ, and so on. The expansion (4) contains a finite number of terms and it is always exact. Usually the HDMR expansion converges very rapidly, and the expansions up to second-order are often enough to describe the outputs of many realistic systems. Of interest is the cut-HDMR expansion which is an exact representation of f ðxÞ, in the hyperplane passing through a reference point in the variable space. With the cut-HDMR, first a reference point x0 is selected. The optimal component functions of cut-HDMR in (4) possess the following structure [29] f 0 ¼ f ðx0 Þ f i ðxi Þ ¼ f ðxi ; xi0 Þ  f 0 f ij ðxi ; xj Þ ¼ f ðxi ; xj ; xij0 Þ  f i ðxi Þ  f j ðxj Þ  f 0 :::

ðxi ; xj ; xij0 Þ ¼ ðx0;1 ; :::; x0;i  1 ; xi ; x0;i þ 1 ; :::; x0;j  1 ; xj ; x0;j þ 1 ; :::; x0;n Þ

ð5Þ

ð6Þ

The cut-HDMR component functions are defined along some cut lines, planes, subvolumes, etc. across the reference point T x0 ¼ f x01 x02 ::: x0n g . This is the origin of the name. The choice of the reference point x0 does not affect the convergence of the HDMR expansion, however, in practice, it is chosen in the neighborhood of interest. In a typical response surface based on HDMR (or quite similar formulations), the series expansion (4) is applied directly to the limit state function (LSF) gðuÞ where the reference point is chosen as the origin of the normal standard space or the design point. The first-order components are obtained through interpolation or regression, and a good estimate of the failure probability and of the reliability index may be obtained. However, this approach is hardly applicable to nonlinear stochastic dynamic analysis, since we need 5–7 points for axes, which would imply only for the definition of the first-order cut components at least 5n analyses for each threshold, being n the number of basic random variables. The task can be accomplished using a very simple trick that will be described in detail in the following subsection. 3.2. The proposed Gaussian Model (GM) The computational effort required by the HDMR in nonlinear stochastic dynamic analysis is greatly simplified if instead of applying (4) to the target LSF (3), it is applied directly to Xðt; uÞ,

ð8Þ

ð9Þ

where we chose as an interpolation point u^ i ¼ 1. Consider now a nonlinear system. A good interpolation of X i ðt; ui Þ is given by a polynomial model of degree d Z 2, so that we would need an increased number of points for each random variable. However, numerical experimentation has shown that an adequate approximation of the response may be obtained using the linear model (9) and keeping as an interpolation point u^ i ¼ 1. Therefore, taking account of (7) and (9) the proposed approach gives a linear response surface n

X G ðt; uÞ ¼ ∑ ci ðtÞui ¼ cðtÞ  u

ð10Þ

i¼1

corresponding to a Gaussian approximation of the response, of variance s2G ðtÞ ¼ jjcðtÞjj2 . In Fig. 1 we represented the curves exact and fitted of X 330 ðt; u330 Þ and X 475 ðt; u475 Þ relative to a hysteretic oscillator, with the input discretized in terms of random pulses, and correspondent to variables far (Fig. 1a) and near (Fig. 1b) the considered time instant t. In the figures we represented the target first-order components X i ðt; ui Þ (continuous line), the interpolated linear functions X^ i ðt; ui Þ (dashed line), and the corresponding first-order components X lin i ðt; ui Þ of the linear system (dotted line). It is possible to see that the greater contribution to the response is given by the components X i ðt; ui Þ with greater i, as expected. Therefore, a reduced number of analyses can be adopted; these analyses can be used for all the degrees of freedom of a MDOF system, so that there is no computational increase for

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

the HDMR-based series expansion, see (5) and (6). In this paper the improvement of the solution given by GM is achieved by using the FORM solution, as described in the next section.

-3

3

x 10

Exact Fitted Linear

2

47

4. Form solution of stochastic dynamic problems X(t,u330)

1

Let us consider again the target LSF (3), Xðt; uÞ being the response of the dynamical system. If we have a linear dynamic system, using the Duhamel integral, the response is given as a linear combination of standard normal random variables Z t n Xðt; uÞ ¼ hðt–tÞFðt; uÞdt ¼ ∑ ai ðtÞui ¼ aðtÞ  u ð12Þ

0

-1

-2

-3 -5

0

-4

-3

-2

-1

0

1

2

3

4

5

u330 0.02 Exact Fitted Linear

0.015 0.01

X(t,u475)

0.005 0

un ðx; tÞ ¼ arg min fjjujj : gðx; t; uÞ ¼ 0g

-0.005 -0.01 -0.015 -0.02 -5

i¼1

where hðtÞ denotes the Impulse‐Response Function (IRF) of the system for the particular input–output pair and the deterministic shape functions ai ðtÞ, are given as convolution integrals of the shape function (sði τÞ ðtÞ or sði ωÞ ðtÞ) of the input. In this simple case the limit state is an hyperplane, the reliability index and the corresponding design point are given as βðx; tÞ ¼ x=jjaðtÞjj and un ðx; tÞ ¼ x aðtÞ=jjaðtÞjj2 . It is here also noted that in this case the shape functions ai ðtÞ of (12) are coinciding with the ci ðtÞ of (10), and therefore P f ;FORM ¼ Φð  βÞ ¼ P f ;GM , coinciding with the exact solution. Conversely, if we have a nonlinear dynamic system, then the limit state is not linear, and the design point is given as a solution of a constrained optimization problem

-4

-3

-2

-1

0

1

2

3

4

5

u475

Fig. 1. First-order components of the HDMR of the system response of a Bouc– Wen SDOF.

the evaluation of the response of a MDOF system with respect to a simple oscillator, except the greater cost of each deterministic analysis. 3.3. Evaluation of the tail probability of the Gaussian Model After the building of the Gaussian model (10) it is possible to determine the surrogate models of the LSF, one for each threshold g G ðx; t; uÞ ¼ x–X G ðt; uÞ ¼ x–cðtÞ  u n

ð11Þ

The design point uG of the surrogate LSF does not require any dynamic computation. In this simple case the surrogate limit state is an hyperplane, then its reliability index and the corresponding design point are known in closed form, namely β G ðx; tÞ ¼ x=jjcðtÞjj and unG ðx; tÞ ¼ x cðtÞ=jjcðtÞjj2 . Note that the tail probability of the Gaussian Model P f ;GM ¼ Prob½X G ðt; uÞ Z x ¼ Φð  βG Þ gives the exact solution for dynamic linear systems, and an approximation for non-linear systems. Numerical experimentation has shown that the GM gives good approximations of the response for weakly non-linear systems or for the lowest thresholds of systems with greater non-linearity. Unfortunately, the accuracy of the response is not guaranteed for the highest thresholds of strongly non-linear systems. This result could be achieved considering the second-order components of

ð13Þ

being the design point the point belonging to the limit state closest to the origin of the standard normal space. The reliability index is given as β ¼ jjun jj and the FORM solution gives a firstorder approximation of the tail probability. It is possible to see that there exists a correspondence one-toone between the design point and the FORM hyperplane [13] aðx; tÞ ¼ x

un ðx; tÞ jjun ðx; tÞjj2

ð14Þ

From (12) and (14) it is seen that the response of a linear system follows a Gaussian distribution with variance equal to jjaðtÞjj2 , while the response of a nonlinear system may be approximated for each response threshold x by a Gaussian distribution with variance equal to jjaðx; tÞjj2 . Consequently, FORM gives a NonGaussian approximation of the system response, giving rise to an improved solution with respect to any Gaussian linearization technique, especially at the tails of the distribution. The main challenge of FORM is however the evaluation of the design point of (3) for each threshold, indeed the optimization problem (13), especially in very high dimensional spaces requires very high computational effort. The gradient-free optimization algorithms, like the Cross–Entropy algorithm [31–33], or genetic algorithms [34] are indeed too demanding. The most efficient procedures to accomplish this task are gradient-based [18,19], which at each step of the iterative procedure require the gradient evaluation of the response. Typically, the finite elements codes do not provide these gradients, which can be approximated by the finite differences method (FDM). By using FDM, each gradient at each iteration of the iterative procedure requires n þ1 nonlinear dynamic computations, n being the number of random variables, and consequently the cost of computation may be excessive. Moreover the value of the perturbation parameter in FDM is a crucial choice, and often numerical noise may be caused by wide tolerance limits in the nonlinear Finite Element (FE) analysis, or by the low number of significant digits in the in the numerical value

48

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

Fig. 2. Description of MCFM for the evaluation of the design point.

of the response [35]. Under these circumstances, the FDM gives inaccurate values of the gradients, which lead to imprecise search directions, and may prevent the algorithm to converge to the design point. It is however worth underlining that increasing the number of significant digits or choosing narrow tolerance criteria in the FE code may be detrimental for the convergence, too; then the analyst has to choose carefully both the perturbation parameter and the tolerances of the FE code. This shortcoming can be solved through the Direct Differentiation Method [21–23] which adopts analytical differentiation of the discretized response equations. In DDM we are guaranteed about consistency, accuracy and efficiency of the gradients. The sensitivities have to be evaluated with an approximation consistent to response itself, and this is obtained with DDM by differentiating the finite element response equations discretized in time and space. The accuracy necessary to the convergence toward the design point is obtained in DDM since the exact derivatives of the response are computed. The efficiency is achieved because the sensitivities are obtained from a linear equation upon convergence of the response of the FE analysis. The main drawback of the DDM is the need to modify the finite element code in order to incorporate for each kind of finite element the exact derivatives. It is therefore necessary to adopt a software, such as OpenSees [24], where the DDM is implemented. However, inside a reliability framework, it may be desirable to adopt the FE code as a “blackbox” which gives the response of the system for each realization of the random variables; this may become necessary when the analyst has to solve problems which require specialized FE codes or elements, where the DDM has not been developed yet. Moreover usually the commercial software packages do not allow to access and modify the code to the finite element level, so that it is necessary to apply the FDM, with the drawbacks above described. It is also noted that to achieve the convergence to the design point the LSF has to be smooth, which means that its derivative with respect to all the random variables is continuous in the space of the random variables. It has been noted that when the material model has a sudden change of slope, such as the uniaxial bilinear stress–strain model, then it is likely that the discontinuities in the gradient will arise, too; in such case the convergence can be obtained smoothing the material model [20]. Unfortunately, using finite element models the detection of discontinuities of the gradient is not an easy task, and so the requisite of smoothness case may be easily violated. Another difficulty for obtaining the FORM solution is that for some realizations of the random variables the FE code fails during the evaluation. Indeed, during the iterative procedure, the

optimization algorithm can require the response of a system in the non-convergence zone of the FE analysis. In this paper the convergence to the design point or to its close approximation is obtained with low computational cost using a suitable response surface strategy, presented in [11,12] and based on the Model Correction Factor Method (MCFM) [36–38]. Instead of selecting a more or less complicated mathematical form, the RS is chosen by idealizing (or over-idealizing) the system under study. If the simplified model is able to capture the main features of the complicated model, and it should do so because both models represent the same physical model, it is expected that the design point usð1Þ of the idealized (or simple) model is not too far from the target design point un . Then, the simple model is moved through an iterative procedure such that the design points of the two models are as close as possible. A natural choice for the simple model is clearly represented by a linear SDOF, to which corresponds as a mathematical model a hyperplane of orientation rðtÞ in the space of the normal standard random variables. Under this hypothesis, the iterative algorithm of the MCFM simply moves the design point uð1Þ s of the simple model along the direction e^ ðtÞ ¼ rðtÞ=jrðtÞj onto the LSF gðx; t; uÞ ¼ 0 of the realistic model, obtaining so the point uns which is at the same time the design point of the “shifted simple model” and a point belonging to the target LSF, implying gðx; t; uns Þ ¼ 0. The procedure is illustrated in Fig. 2. This approach solves the above cited drawbacks, because (i) it doesn0 t require the evaluation of any gradient, (ii) the simple model is simply shifted onto the target LSF, which can therefore be non-smooth, (iii) if the simple model is well chosen, it is likely that the algorithm require to the FEM solver to run analyses only in the convergence zone. 4.1. Determination of the simple model A major issue in an MCFM-based approach is clearly the choice of the simple model. It is noted that the design point is the point belonging to the limit state closest to the origin of the standard normal space, therefore the reliability index of the simple model βs ¼ jjuns jj is an upper bound to the target reliability index β ¼ jjun jj, that is

β r βs

ð15Þ

Taking account of (15), in [11,12] it has been adopted as a simple model a linear dynamic system, whose parameters ϑ ¼ fϑ1 ; ϑ2 ; …; ϑm gT are chosen in order to minimize βs ðϑÞ. It is however seen that the method requires a relatively high computational cost if a good choice of the simple model is well represented by an MDOF linear system [12]; moreover, the computations have to be repeated for each quantity of interest. Conversely, in this paper we adopt as a simple model a nonparametric one, born from the Gaussian model previously described. If the behavior of the nonlinear system is not too different from linear, it is likely that unG is not too far from un , being unG the design point of the GM. It is also noted that it is enough to build only one GM for all the degrees of freedom of the nonlinear MDOF system, reducing further the computational cost. The simple model is written as X s ðt; uÞ ¼ rðtÞ  u

ð16Þ

with rðtÞ ¼ cðtÞ þ wðtÞ

ð17Þ

where cðtÞ is the gradient vector of the Gaussian model (10), respectively, while wðtÞ is a correction term. To obtain the hyperplane ^ ðkÞ , k¼1,2,…,m, mrn, X s ðt; uÞ we select a set of sampling points u located in each direction at a small distance from unG . At each of these

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

49

points we match the real and the idealized limit state surfaces, which ^ ðkÞ Þ ¼ X s ðt; u ^ ðkÞ Þ. Note that the number m of support implies Xðt; u points can be less than the number n of random variables. Assuming unG is not too far from the target design point un , the most important m random variables evaluated at unG are approximately the same as the most important random variables at un . Thus, after evaluating the importance measures of the random variables according to the Gaussian model, and subsequent sorting, the support points are chosen along the m most important directions. 4.2. The MCFM based on FORM After evaluating the simple model X s ðt; uÞ it is possible to build the corresponding LSF g s ðx; t; uÞ ¼ x–X s ðt; uÞ ¼ x–rðtÞ  u ð1Þ s

ð18Þ uð1Þ s

¼ x r=jjrjj2 . with reliability index β ¼ x=jjrjj and design point In general, uð1Þ does not belong to the realistic limit state s gðx; t; uÞ ¼ 0, of course. Define now the “exactly corrected simple model” g~ s ðx; t; uÞ ¼ x  νðuÞX s ðt; uÞ

ð19Þ

where νðuÞ is known as “effectivity factor” [36–38]. It is always possible to choose νðuÞ such that the LSF (19) is equal to (3), but is also as difficult. However, since we are interested in approximating only a narrow region around the design point, it is reasonable to substitute the effectivity factor with a constant νðuÞ ffi νn , giving rise to the “corrected simple model” g 0s ðx; t; uÞ ¼ x–νn X s ðt; uÞ ¼ x–νn ½rðtÞ  u

ð20Þ

being νn the “model correction factor”. The latter is determined through an iterative procedure such that the design point uns of the corrected simple model g 0s ðx; t; uÞ belongs to the limit state of the realistic model gðx; t; uÞ. It is noted that since X s ðt; uÞ is a linear ð1Þ n n model, then βs ¼ β s =νn and uns ¼ uð1Þ s =ν , where ν is determined n such that g s ðx; t; us Þ ¼ 0. If the simple model X s ðt; uÞ is not too far from the target one Xðt; uÞ, then it is expected that uns is quite close to un , see Fig. 2. Note that in view of (15) the MCFM gives an upper bound of the reliability index and correspondingly we obtain a lower bound (i.e. unsafe) of the FORM approximation of the tail probability. The quality of the approximation can be further improved building a new simple model X s ðt; uÞ centered in uns . Numerical experimentation has shown that usually a few iterations suffice.

Fig. 3. Sample loop hysteresis (a) and tail probability (b) of a Bouc–Wen SDOF (α ¼ 0:05, r ¼ 1, δ ¼ γ ¼ 0:5) with the excitation discretized in terms of random pulses.

corresponding cut-off frequency of this stochastic excitation is

ωc ¼ π =Δt ¼ 157:08 rad/s.

5. Numerical investigation 5.1. Bouc–Wen SDOF system subjected to white noise—Excitation expressed in terms of random pulses Consider a Bouc–Wen [39] hysteretic oscillator defined by the differential equations mX€ þ cX_ þk½αX þ ð1  αÞZ ¼ FðtÞ r1 _ Z_ ¼  δjXjjZj Z  γ jZjr X_ þAX_

ð21Þ

where m, c and k are the mass, damping and stiffness coefficients, α controls the degree of hysteresis, and r, A, γ and δ are parameters defining the shape of the hysteresis loop. The excitation is defined as FðtÞ ¼  mag ðtÞ, where ag ðtÞ denotes the base acceleration modeled as a white-noise process of intensity SW ¼ 0:0117 m2/s3. Let s0 ¼ m2 π SW =ðckÞ the mean-square response of the linear system, obtained setting α ¼ 1. The excitation has been modeled using representation (1), the time step Δt ¼ 0:02 s is used and the response at time t ¼ 10 s is considered; consequently, for the case under exam, we have n ¼ nðτÞ ¼ 500 random variables. The

The result is calculated for 6 normalized threshold values x=s0 ranging from 0.5 to 3 with increments of 0.5. It is seen that the most important n^ ¼ 200 first-order components of the HDMR series expansion are enough, so that for all the considered cases the Gaussian Model required only 200 dynamic computations. At first we chose

α ¼ 0:05;

r ¼ 1;

A ¼ 1;

δ ¼ γ ¼ 0:5

ð22Þ

In Fig. 3(a) the hysteresis loop of the oscillator for a sample of the stochastic input is shown. In Fig. 3(b) we represented in semilogarithmic scale the tail probability of the linear system, obtained by setting α ¼ 1 (dash-dot line with diamond markers), together with the GM (dashed line). It is seen that the latter agrees very well with the MCS with 100,000 samples (circle markers). This happens because the hysteresis loop holds approximately a constant slope both for low and high thresholds. Then we considered a case where FORM solution is needed. We chose

α ¼ 0:5;

r ¼ 3;

A ¼ 1;

δ ¼ γ ¼ 1=ð2sr0 Þ

ð23Þ

50

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

In Fig. 4(a) we represented the hysteresis loop of the oscillator for a sample of the stochastic input, while in Fig. 4(b) it is shown in semilogarithmic scale the tail probability of: (i) the linear system (dash-dot line with diamond markers), (ii) the GM (dashed line), (iii) the FORM solution (continuous line), (iv) the MCFM with m ¼ 150 analyses per threshold (dotted line), (v) the MCS with 100,000 samples (circle markers). The MCFM required a total amount of 200 þ 150  6 ¼ 1100 analyses. 5.2. Bouc–Wen SDOF system subjected to white noise—Excitation expressed in terms of harmonic components If one is interested in stationary responses of the system, the discretization (2) can be used, which requires a less computational cost. Consider now the Bouc–Wen SDOF with parameters (23) By choosing t ¼ 10 s, we have Δω r ð2π Þ=t ¼ 0:628 rad/s. If we select Δω ¼ 0:628 rad/s, correspondingly we have nðωÞ ¼ 500=2 ¼ 250 harmonic components, and the cut-off frequency ωc ¼ nðωÞ Δω ¼ 157:08 rad/s. Under these conditions, the SDOF Bouc–Wen with the discretization (2) requires n ¼ 2nðωÞ ¼ 500 random variables, likewise the discretization (1). However, in this example the stationarity is reached after t ¼ 5=6 s, so it is possible to choose a wider frequency step, namely Δω ¼ 1 rad/s. Moreover, we are not required to consider all the range of frequencies, since the fundamental frequency of the linear SDOF system is ω0 ¼ 8:36 rad/s. So it is possible to choose ωmax ¼ 16 rad/s, to which correspond only n^ ¼ 2n^ ðωÞ ¼ 32 random variables. The results are shown in Fig. 5, here we mention that the GM required 32 dynamic computations, while the MCFM 32 analyses per threshold, for a total amount of 32 þð6  32Þ ¼ 224 analyses. 5.3. Bouc–Wen MDOF system subjected to filtered white noise

Fig. 4. Sample loop hysteresis (a) and tail probability (b) of a Bouc–Wen SDOF (α ¼ 0:5, r ¼ 3, δ ¼ γ ¼ 0:5=sr0 ) with the excitation discretized in terms of random pulses.

The method can be easily applied to MDOF systems, indeed the n^ analyses required for the evaluation of the GM can be used for all the dofs of the system. As an application to a multi-degree-of-freedom (MDOF) system, we considered a 6-story shear-building model, see Fig. 6, subjected to a stochastic base acceleration, and presented in [12]. The hysteretic behavior in each story is idealized by a Bouc–Wen model so that the equations of motion for the ith story are described as follows: ! i N € € mi ∑ Q þ U g þ ∑ C ij Q_ þr i  r i þ 1 ¼ 0; i ¼ 1; 2; …; N ð24Þ j¼1

Fig. 5. Tail probability of a Bouc–Wen SDOF (α ¼ 0:5, r ¼ 3, δ ¼ γ ¼ 0:5=sr0 ) with the excitation discretized in terms of harmonic components.

j

j¼1

j

r i ðQ i ; Q_ i ; Z i Þ ¼ αi ki Q i þ ð1  αi Þki Z i

ð25Þ

Z_ i ¼  δi jQ_ i jjZ i jni  1 Z i  γ i jZ i jni Q_ i þ Ai Q_ i

ð26Þ

In the above, N is the number of stories, mi are floor masses, C ij are viscous damping coefficients, r i are story resisting forces, Q i are interstory displacements, Z i are hysteretic components of the interstory displacements, and U€ g ðtÞ denotes the base acceleration, which is modeled as a filtered white-noise process with filter parameters ωf ¼ 5π rad/s, ζ f ¼ 0:6 and the intensity of the underlying white noise as SW ¼ 0:00332 m2/s3. At first the discretization in terms of random pulses is adopted. The time step Δt ¼ 0:02 s is used and the response at time t ¼ 10 s is considered. Thus, n ¼ 501 standard normal random variables define the stochastic input. For the GM we used the most important n^ ¼ 200 random variables; the results shown in Fig. 7 indicate that the GM with only 200 dynamic computations gives very good approximations for both the first and the sixth interstory drift. The FORM solution required only further 50 analyses per

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

51

Fig. 6. MDOF Bouc–Wen system.

each threshold, while the target solution, represented from the MCS, needed 100,000 samples. Likewise for the SDOF system, we can obtain a saving of computational cost, by expressing the stochastic input in terms of harmonic components. Let the frequency step be Δω ¼ ð2π Þ=t ¼ 0:628 rad/s; the first two frequencies of the elastic system are ω1 ¼ 10:84 rad/s and ω2 ¼ 26:54 rad/s. So we built three different discretizations for the evaluation of the GM, the first one ð2Þ X ð1Þ G ðt; uÞ with ω1; max ¼ 15 rad/s, the second one X G ðt; uÞ with ð3Þ ω2; max ¼ 40 rad/s, and the last one X G ðt; uÞ with ω3; max ¼ 100 rad/ s. The Gaussian Models required a total amount of n^ 1 ¼ 46, n^ 2 ¼ 126 and n^ 3 ¼ 318 dynamic computations, respectively. The second and the third Gaussian Model give the same tail probability, which differs from the first one. The results in Fig. 8 show that the only the model X ð1Þ G ðt; uÞ does not give a safe approximation; moreover it has been noted that the target FORM solution can be obtained by MCFM only with a number of support points m Z n^ 2 . This means that the system is governed by first two modes, which confirms previous results obtained for this numerical application [12].

6. Concluding remarks In this paper we have adopted a novel response surface strategy for nonlinear stochastic dynamic analysis. At first, a

Fig. 7. Tail probability of the first (a) and the sixth (b) interstory drift of a Bouc– Wen MDOF system with the excitation discretized in terms of random pulses.

response surface based on the High Dimensional Model Representation is adopted to build a good Gaussian approximation of the response, called Gaussian Model (GM). Numerical experimentation has shown that it works quite well for weakly nonlinear systems, and for the lowest thresholds of strongly nonlinear systems. The GM requires a low computational effort, as shown from the application to hysteretic systems, and summarized in Table 1. The developed examples have also shown that if one is interested to the stationary response of the system, the discretization in the frequency domain allows to reduce significantly the computational effort, especially if we have a first-mode dominated structure. For the highest threshold of structures with higher degree of nonlinearity the FORM solution may be necessary. This happened for the Bouc–Wen example of parameters(23) A good approximation of the design point is obtained by using a variant of the Model Correction Factor Method (MCFM), and choosing as a starting model the GM. The joint use of GM with MCFM can be then considered also as a response surface method for the evaluation of the design point excitation, as alternative to the procedure presented in [12]. However, the computational effort required by the proposed approach is clearly less, see Table 1.

52

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

not possible to know in advance the level of error of the achieved approximation. Research to overcome these limitations is in progress [40]. Another issue which deserves further study is represented by the application of the GM to real-world engineering MDOF systems, with different degree and kind of nonlinearities.

Acknowledgements The author wishes to thank Prof. Armen Der Kiureghian of the University of California at Berkeley for his suggestion to explore the discretization of the input in terms of harmonic components. References

Fig. 8. Tail probability of the first (a) and the sixth (b) interstory drift of a Bouc– Wen MDOF system with the excitation discretized in terms of harmonic components.

Table 1 Results of the SDOF and MDOF systems. System

Bouc–Wen Eq. (23) Bouc–Wen Eq. (24) MDOF

Number of analyses GM (τ)

GM (ω)

200 200 200

32 32 126

GM (τ) þ MCFM

GM (ω) þ MCFM

MCFM (Ref. 12)

1100 700

224 626

– 1890 4647

In this paper the focus was devoted to the point-in-time distribution of the response, expressed in terms of its Cumulative Distribution Function (CDF). A topic of interest is the time-variant reliability analysis, however the adoption of linear response surfaces allows the straightforward application of the recently developed Tail Equivalent Linearization Method; further research will be devoted to this topic. The reliability approach here presented is based on FORM, and it shares with FORM the main shortcomings: (i) in some strongly nonlinear problems the accuracy of FORM is not excellent, (ii) it is

[1] Roberts JB, Spanos PD. Random vibration and statistical linearization. Chichester: Wiley; 1991. [2] Casciati F, Faravelli L. Hysteretic 3-dimensional frames under stochastic excitation. Res Mech 1989;26(3):193–213. [3] Pradlwarter HJ, Schuëller GI, Schenk CA. A computational procedure to estimate the stochastic dynamic response of large non-linear FE-models. Comput Meth Appl Mech Eng 2003;192(7-8):777–801. [4] Lutes LD, Sarkani S. Random vibrations: analysis of structural and mechanical systems. Burlington (MA): Elsevier; 2004. [5] Au SK, Beck JL. Estimation of small failure probabilities in high dimensions by subset simulation. Probab Eng Mech 2001;16:263–77. [6] Pradlwarter HJ, Schueller GI, Kotsourelakis PS, Charmpis DC. Application of line sampling simulation method to reliability benchmark problems, 2007, 29: 208-221. [7] Bucher C. Asymptotic sampling for high-dimensional reliability analysis. Probab Eng Mech 2009;24:504–10. [8] Der Kiureghian A, Li C-C Nonlinear random vibration analysis through optimization. In: Proceedings of the seventh IFIP WG 7.5 conference on optimization of structural systems, Boulder, Colorado; 1996, 197–206. [9] Der Kiureghian A. The geometry of random vibration and solutions by FORM and SORM. Probab Eng Mech 2000;15:81–90. [10] Koo H, Der Kiureghian A, Fujimura K. Design-point excitation for non-linear random vibration. Probab Eng Mech 2005;20:136–47. [11] Alibrandi U, Der Kiureghian A. The model correction factor method in NonLinear Stochastic Dynamic Analysis. In: Sixth computational stochastic mechanics conference, June 13–16; 2010. Rhodos, Greece. [12] Alibrandi U, Der Kiureghian A. A gradient-free method for determining the design point in nonlinear stochastic dynamic analysis. Probab Eng Mech 2012;28:2–10. [13] Fujimura K, Der Kiureghian A. Tail-equivalent linearization method for nonlinear random vibration. Probab Eng Mech 2007;22:63–76. [14] Der Kiureghian A, Fujimura K. Nonlinear stochastic dynamic analysis for performance-based earthquake engineering. Earthquake Eng Struct Dyn 2009;38:719–38. [15] Garrè L, Der Kiureghian A. Tail equivalent linearization method in frequency domain and application to marine structures. Mar Struct 2010;23:322–38. [16] Broccardo M, Der Kiureghian A. Multi-component nonlinear stochastic dynamic analysis using tail-equivalent linearization method. In: Proceedings of the 15th world conference on earthquake engineering, Lisbon, Portugal; 2012. [17] Broccardo M, Der Kiureghian A. Non stationary stochastic dynamic analysis by tail-equivalent linearization. In: 11th International conference on structural safety and reliability, New York, June 16–20; 2013. [18] Liu P-L, Der Kiureghian A. Optimization algorithms for structural reliability. Struct Saf 1991;9(3):161–77. [19] Zhang Y, Der Kiureghian A. Two improved algorithms for reliability analysis. In: Reliability and optimization of structural systems, proceedings of the sixth IFIP WG 7.5 working conference on reliability and optimization of structural systems, Assisi, Italy, Chapman & Hall; 1995. [20] Haukaas T, Der Kiureghian A. Strategies for finding the design point in nonlinear finite element reliability analysis. Probab Eng Mech 2006;21:133–47. [21] Liu P-L, Der Kiureghian A. Finite element reliability of geometrically nonlinear uncertain structures. J Eng Mech 1991;117(8):1806–25. [22] Zhang Y, Der Kiureghian A. Dynamic response sensitivity of inelastic structures. Comput Meth Appl Mech Eng 1993;108:23–36. [23] Haukaas T, Der Kiureghian A. Finite element reliability and sensitivity methods for performance-based earthquake engineering, peer report 2003/ 14, Pacific Earthquake Engineering Research Center, University of California, Berkeley, CA; 2004. [24] Mckenna F, Fenves GL, Scott MH. Open system for earthquake engineering simulation. Pacific Earthquake Engineering Research Center, University of California, Berkeley, CA, 〈http://opensees.berkeley.edu〉; 2003. [25] Breitung K, Faravelli L. Log-likelihood maximization and response surface in reliability assessment. Nonlinear Dyn 1994;5(3):273–85.

U. Alibrandi / Reliability Engineering and System Safety 126 (2014) 44–53

[26] Faravelli L. Response-surface approach for reliability analysis. J Eng Mech 1989;115:2763–81. [27] Bucher C, Burgound U. A fast and efficient response surface approach for structural reliability problems. Struct Saf 1990;7:57–66. [28] Rajashekhar MR, Ellingwood BR. A new look at the response surface approach for reliability analysis. Struct Saf 1993;12:205–20. [29] Li G, Rosenthal C, Rabitz H. High dimensional model representation. J Physical Chem A 2001;105:7765–77. [30] Chowdury R, Rao BN, Prasad AM. High-dimensional model representation for structural reliability analysis. Commun Numer Methods Eng 2009;25:301–37. [31] de Boer PT, Kroese DP, Mannor S, Rubinstein RY. A tutorial on the crossentropy method. Ann Oper Res 2005;134(1):19–67. [32] Rubinstein RY. The cross-entropy method for combinatorial and continuous optimization. Methodol Comput Appl Probab 1999;1:127–90. [33] Valdebenito MA, Pradlwarter HJ, Schueller GI. The role of design point for calculating failure probabilities in view of dimensionality and structural nonlinearities. Struct Saf 2010;32:101–11.

53

[34] Koh CG, Hong B, Liaw CY. Substructural and progressive structural identification methods. Eng Struct 2003;25:1551–63. [35] Koduru SD, Haukaas T. Feasibility of FORM in finite element reliability analysis. Struct Saf 2010;32:145–53. [36] Ditlevsen O, Arnbjerg-Nielsen T. Model correction factor method in structural reliability analysis. J Eng Mech 1994;120(1):1–10. [37] Ditlevsen O, Johannesen JM. System reliability analysis by model correction factor method. In: ICASP8: eighth international conference on applications of statistics and probability, Sydney, Australia; 1999. [38] Ditlevsen O, Madsen HO. Structural reliability methods. Wiley; 1999. [39] Wen JK. Approximate method for non-linear random vibration. J Eng Div, ASCE, 101; 1975; 389–401. [40] Alibrandi U, An alternative linear response surface improving FORM solution for stochastic dynamic analysis. In: Reliability and optimization of structural systems, IFIP WG 7.5 working conference; 24–27 June 2012. Yerevan, Armenia.