Reliability-based design optimization of shells with uncertain geometry using adaptive Kriging metamodels
RReliability-based design optimization of shells with uncertaingeometry using adaptive Kriging metamodels
V. Dubourg , J.-M. Bourinet , and B. Sudret Phimeca Engineering, Cournon d’Auvergne, France Universit´e Clermont Auvergne, Sigma Clermont, Institut Pascal, Clermont Ferrand, France CNRS, UMR 6602, Institut Pascal, Aubi`ere, France Institute of Structural Engineering, ETH Z¨urich, Z¨urich, Switzerland
Abstract
Optimal design under uncertainty has gained much attention in the past ten years due tothe ever increasing need for manufacturers to build robust systems at the lowest cost.Reliability-based design optimization (RBDO) allows the analyst to minimize some cost function whileensuring some minimal performances cast as admissible failure probabilities for a set of perfor-mance functions. In order to address real-world engineering problems in which the performanceis assessed through computational models (e.g. finite element models in structural mechanics)metamodeling techniques have been developed in the past decade. This paper introduces adap-tive Kriging surrogate models to solve the RBDO problem. The latter is cast in an augmentedspace that “sums up” the range of the design space and the aleatory uncertainty in the designparameters and the environmental conditions. The surrogate model is used (i) for evaluatingrobust estimates of the failure probabilities (and for enhancing the computational experimentaldesign by adaptive sampling) in order to achieve the requested accuracy and (ii) for applyinga gradient-based optimization algorithm to get optimal values of the design parameters. Theapproach is applied to the optimal design of ring-stiffened cylindrical shells used in submarineengineering under uncertain geometric imperfections. For this application the performance ofthe structure is related to buckling which is addressed here by means of a finite element solutionbased on the asymptotic numerical method.
Keywords : reliability-based design optimization – Kriging surrogates – shell buckling –geometric imperfections – asymptotic numerical method a r X i v : . [ s t a t . M E ] A p r Introduction
Shell structures occupy a predominant part of our landscape (see e.g.
Ramm and Wall, 2004, fora review of their applications). They owe this predominance to their curvature which allow themto withstand large transverse loading by a membrane-dominated stress state. As a result, theycan be used for building large-span shelters such as roofs, fuselages or boat and submarine hullswithout requiring too many intermediate supports such as stiffening beams or rims. Nonetheless,as many optimized and therefore slender structures, the strength of thin shells also exhibits asignificant sensitivity with respect to geometrical, material and other environmental conditionswhich are typically unknown to some extent.Early work on the elastic stability of slender structures (such as beams, plates or shells) isoften attributed to Euler in 1744, although most of the theoretical concepts in practice today forshells are due to Lorentz, Timoshenko and Southwell in early nineties. In parallel to theoreticaladvances, experimental studies revealed embarrassing discrepancies between the predicted buck-ling loads and those obtained from real tests. Koiter (1945) was certainly the first researcherto point out that these discrepancies are mostly explained by the imperfect geometry, boundaryconditions and material properties of the experimental specimens. This premise is now fullyacknowledged by the whole community of engineers and scientists in structural mechanics in thelight of other studies by Arbocz and Babcock (1969); Singer et al. (1971); Singer and Abramovich(1995) amongst others. The reader may refer to (Baˇzant, 2000) for a review of works in the fieldof stability of structures with an emphasis on anelastic structures and to the recent paper fromElishakoff (2012) for a detailed history of works on elastic stability of shells.A key aspect of these imperfections though is that they are extremely varying in termsof shape and amplitude. Hence, for the sake of structural safety, designers have to accountfor extreme and fortunately unlikely imperfections. A common practice is to assume a givenshape in the calculations corresponding to the worst case structural strength and then resort toadvanced numerical schemes in order to justify the design. However, this approach, referred toas the worst case approach in the sequel, introduces an unknown degree of conservatism whichmay not suit the safety requirements imposed by stakeholders.As early suggested by Bolotin (1962) in his pioneering works, it is argued that a better so-lution may be obtained by means of statistical methods and that the design of imperfect shellsnecessarily falls under a probabilistic formulation. Several imperfection surveys were later car-ried out in order to assess the statistical properties of imperfections present in both small andlarge-scale shells. These statistics such as those gathered in the imperfection data bank (Arbocz,1982) were then introduced into stochastic buckling analysis by researchers. Elishakoff (1979)was the first researcher to use random initial imperfections of compressed cylindrical shells ina Monte Carlo analysis. The initial imperfections were expanded in double Fourier series and he Fourier coefficients were considered as Gaussian random variables. Such a representation ofthe random imperfections is also used in further studies carried out by the same authors, seee.g. (Elishakoff and Arbocz, 1982; Elishakoff et al., 1987; Arbocz and Hol, 1995; Arbocz andHilburger, 2005). This double Fourier representation was later investigated by other researcherswith an effort to reduce the number of random Fourier coefficients in order to alleviate the cost ofthe stochastic analysis (Noirfalise, 2009; Kriegesmann et al., 2010, 2011). In the works of Schenkand Schu¨eller (Schenk and Schu¨eller, 2003, 2007), the geometric imperfection was modeled as atwo-dimensional univariate non-homogeneous Gaussian random field by means of a KarhunenLo`eve expansion. This representation is also used in references (Craig and Roux, 2008; Noir-falise, 2009; Dubourg et al., 2009). As an alternative technique, the spectral representationmethod was applied in several works for modeling two-dimensional univariate random fields,see e.g. the work by Stefanou and Papadrakakis (2004) which models homogenous Gaussianrandom fields representing spatially-varying material properties (Young’s modulus and Poisson’sratio) and thickness imperfections. The spectral representation method is also used in con-junction with an autoregressive moving average technique with evolutionary power spectra in(Papadopoulos and Papadrakakis, 2005) for modeling non-homogenous Gaussian random fieldsrepresenting spatially-varying Young’s modulus in addition to geometric and thickness imperfec-tions. A similar approach based on non-Gaussian translation fields is adopted by Papadopoulosand Papadrakakis (2005), which puts the emphasis on the influence of the Gaussianity/non-Gaussianity assumption on the results. A recourse to evolutionary power spectra estimated bymeans of a moving window averaging technique is also found in (Broggi and Schu¨eller, 2011) fornon-homogenous random fields representing geometric and thickness imperfections. Most of theworks reported in the literature focused on metallic shells with geometric imperfections possiblycombined with spatially-varying material properties and thickness imperfections. A few researchstudies were also conducted on anisotropic composite shells (Chryssanthopoulos and Poggi, 1995;Kriegesmann et al., 2010; Broggi and Schu¨eller, 2011; Kriegesmann et al., 2011) characterizedby larger imperfections due to their complex manufacturing processes. For a more realistictreatment of imperfections, some other sources of random imperfections were additionally incor-porated in probabilistic buckling studies such as those arising from a non-uniform distributionof the axial loading (Papadopoulos and Iglesis, 2007) or those coming from the application ofuncertain boundary conditions (Schenk and Schu¨eller, 2007). It is of importance to mention thatmost of the models used for describing the geometric imperfections in the above cited referenceswere identified from the experimental data available in the imperfection data bank. For thicknessimperfections and spatially-varying material properties, the parameters of the random fields areassumed to have specified values and they are sometimes varied in a parametric analysis.In early probabilistic studies, buckling loads were computed by means of Koiter’s theory Elishakoff, 1979; Elishakoff and Arbocz, 1982). More refined numerical solutions based on amultimode analysis were later used by Elishakoff et al. (1987) and Arbocz and Hilburger (2005).A recourse to a nonlinear finite element (FE) model is often advocated for an enhanced accuracyon the limit loads. FE-based probabilistic approaches were carried out with STAGS in (Arboczand Hol, 1995; Schenk and Schu¨eller, 2003, 2007), ABAQUS in (Kriegesmann et al., 2010; Broggiand Schu¨eller, 2011; Kriegesmann et al., 2011) and LS-DYNA in (Craig and Stander, 2007;Craig and Roux, 2008) (note that this last code is specifically used in the context of dynamicbuckling). In other studies, some other researchers implemented specific shell elements for thepurpose of their probabilistic buckling analysis. The TRIC triangular shell element as describedin (Argyris et al., 2002) is used in all the works carried out at the National Technical Universityof Athens (Stefanou and Papadrakakis, 2004; Papadopoulos and Papadrakakis, 2005; Lagarosand Papadopoulos, 2006; Papadopoulos and Iglesis, 2007; Papadopoulos et al., 2009). In theirworks, Noirfalise (2009) and Dubourg (2011) made a recourse to the B¨uchter and Ramm 8-node shell element (B¨uchter et al., 1994) for their reliability analysis and reliability-based designoptimization. Most of the probabilistic buckling analyzes found in the literature assume linearelasticity of metallic shells, which is perfectly appropriate for the studied thin shells taken fromthe imperfection data bank such as the so-called A-shells. A nonlinear behavior of the shellmaterial is addressed in (Papadopoulos and Papadrakakis, 2005; Lagaros and Papadopoulos,2006; Dubourg et al., 2009). The Young’s modulus is considered as a random field in all thesereferences. The yield strength is taken as an additional random field independent of the Young’smodulus one by Dubourg et al. in (Dubourg et al., 2009).In several studies, the probabilistic analysis consisted in constructing limit load histogramswith comparison to those obtained experimentally and analyzing the second-order statistics ofthese loads. Note that the numbers of samples used in the FE-based Monte Carlo simulationsof the reported references are most often in the order of a few hundreds. Another directionfollowed by researchers has consisted in performing a structural reliability analysis of shells byimposing that their limit loads should be greater than a prescribed service load. The earliestoccurrences of such studies were based on the first-order second-moment method (Elishakoffet al., 1987; Arbocz and Hol, 1995; Arbocz and Hilburger, 2005). The first-order reliabilitymethod (FORM) was later used in many works based on analytical models not listed here forthe sake of brevity or on FE models, see e.g. (Bourinet et al., 2000; Dubourg et al., 2009). For thespecific case of imperfections modeled by random fields, a recourse to subset simulation method(Au and Beck, 2001) is considered as the most suitable solution as investigated by Noirfalise(2009) and Dubourg et al. (2009). For the purpose of improved designs, the optimization of shellsubject to buckling based on FE models has also been of interest in recent years. In (Lagarosand Papadopoulos, 2006), the weight of shells with random geometric imperfections and space- ariant Young’s modulus and thickness is minimized in the framework of reliability-based designoptimization. This work is based on a (5 + 5)-evolution strategy optimization algorithm and theprobabilistic constraint is assessed by means of a crude MC with 1 ,
000 samples. In (Craig andRoux, 2008), shells with stochastic imperfections are optimized in a dynamic buckling context.Two optimization studies are performed: the first one aims at minimizing the weight of the shellwith constraints on the average peak force and average internal energy, the second one aims atincreasing the robustness w.r.t. the variations of the normal peak force with the same constraints.The strategy used by Craig and Roux is to construct a quadratic polynomial response surfacewith a 96-sample MC simulation at each point of the design of experiments.In this article, it is proposed to address the optimization under uncertainty of the weightof ring-stiffened cylindrical shells representative of those used in submarine pressure hulls. Thegeometric imperfection of a single bay is considered as random in the shape of a few selectedcritical buckling modes as previously addressed in (Dubourg et al., 2008). The uncertainty inmaterial properties and thicknesses is accounted for by means of a random variable approach.This uncertainty is therefore not modeled as space-variant random fields as reported in some ofthe previously cited works. The parameters of the random models are chosen in accordance torequirements imposed by standards or general recommendations, if available, or arbitrarily setup to prescribed values. The limit loads are accurately assessed by FE with a non-incrementalnon-iterative method known as the asymptotic numerical method (ANM). A nonlinear elasticbehavior of the material is assumed in the FE analysis with random material properties (Young’smodulus, yield strength and ultimate strength). An approach based on semi-numerical modelsis also proposed as an alternative for computing the limit loads. The weight optimization ofthe single bay design is carried out under the constraint that the failure probability of theimperfect shell w.r.t. buckling is kept lower than a prescribed target failure probability. Theproposed approach is therefore a reliability-based design optimization (RBDO) as known in theliterature, see e.g. (Tsompanakis et al., 2008). Since the limit loads are predicted by means ofa costly-to-evaluate FE analysis, the following constraints are taken into account: • the proposed approach shall be applicable within a few hundred runs of the FE modelonly: this implies resorting to metamodeling techniques. A recourse to an adaptive Krigingmethod developed in a so-called augmented space is made (Dubourg et al., 2011) . • the probability of failure shall be evaluated precisely at each step of the optimizationalgorithm and it is likely to be small (e.g. less than 10 − for the forecast applications).Thus subset simulation (Au and Beck, 2001) is used.This paper is organized as follows. Section 2 first introduces the fundamental concepts ofnonlinear stability analysis in structural mechanics. It also defines the assumptions that areaccounted for in the subsequent FE-based probabilistic buckling analyses. FE-based reliability nalyses are known to be too-computationally-demanding for being applicable to industrial con-cerns. In this respect, Section 3 exposes the key concepts of the Kriging surrogate-based strategy.In Section 4, RBDO is applied to a single bay of a submarine pressure hull with random geomet-ric imperfections. This probabilistic design philosophy is opposed to the state-of-the-art worstcase approach. Buckling is a structural instability phenomenon triggered by some excessive load that needs tobe identified. This load will be referred to as the critical buckling (or collapse) load in the sequel.In practice it is determined by applying a so-called load proportionality factor (LPF) λ which isinitialized to zero and then incrementally increased until collapse is observed. In continuum mechanics, the equilibrium state of a conservative mechanical system is character-ized by a zero elementary variation of its total potential energy denoted by E t . This fundamentalprinciple leads to the establishment of the following so-called variational formulation of equilib-rium states: δE t = E t , u ( λ, u ) δ u = 0 , (1)where u denotes any admissible displacement of the structure and E t , u is the first-order func-tional derivative of the total potential energy that depends on the LPF λ . The infinite set ofvalues of λ and u satisfying the latter equation is known as the equilibrium path of the struc-ture. This path is usually constructed incrementally from a known initial state λ (0) , u (0) , e.g.the reference state of the unloaded structure for which λ = 0.For stable structures, the only state of interest corresponds to a unit value of the LPF λ . Unstable (resp. stable) equilibrium states are characterized by a negative (resp. strictlypositive) second-order functional derivative of the total potential energy E t , u u , meaning thatthey correspond to local maxima (resp. minima) of this energy.There exists basically two kinds of instabilities, both potentially leading to buckling and/orpremature plastic collapse: bifurcation points and limit points, see Figure 1 for illustration.Regarding bifurcation points, the structure may lose its stability along the equilibrium path,resulting in sudden and large displacements which often lead to collapse. Regarding limit points,this occurs when the structure is no longer able to withstand loads due to nonlinear geometricaland/or material effects. For many structures including shells, these two kinds of points generallyinteract in a joint manner, one triggering the other and conversely. egular pointBifurcationLimit load point u Bifurcated path Imperfection w/increasing amplitude
Figure 1: Equilibrium paths and stability
Practical detection of these instabilities is a non-trivial task and it involves the resolution of aperturbed equilibrium problem along with the resolution of Eqn. (1). The present study focuseson imperfect structures which are commonly assumed to fail at their limit load. Indeed, Koiter(1945) showed that the presence of initial imperfections in the structure smooths the equilibriumpaths. As a result, imperfect structures regularly feature smooth limit load points rather thansharp bifurcation points which are only observed on perfect structures (see Figure 1). Therefore,in the sequel, the detection of singular points along the equilibrium path will be restricted tolimit load points characterized by horizontal tangents of the equilibrium path. In a static analysis, the total potential energy of a structure of volume V is given by: E t ( λ, u ) = (cid:90) V W int ( ε ) d v − λ W ext ( u ) , (2)where W int is the strain energy density in the structure, W ext is the work of external forces andd v is the infinitesimal volume. ε stands for the Green-Lagrange strain tensor which is definedas: ε = ε ( u ) = 12 (cid:0) ∇ u + ∇ T u (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) ε l ( u ) + 12 ∇ u ∇ T u (cid:124) (cid:123)(cid:122) (cid:125) ε nl ( u , u ) , (3)where ε l ( u ) (resp. ε nl ( u , u )) denotes the linear (resp. symmetric quadratic bilinear) term of ε and ∇ the gradient operator. Assuming linear elasticity, the strain energy density W int reduces The initial study was carried out on cylinders under a compressive axial load with modal imperfections, but thisgeneralizes to other imperfect structures. o the following quadratic form: W int ( ε ) = 12 ε : D : ε , (4)where D is the elasticity tensor of the material and the symbol : denotes the double contractionof tensors.Introducing the second Piola-Kirchhoff stress tensor S = D : ε , the variational formulationof the equilibrium equation rewrites as the following set of equations: δE t = (cid:90) V S : δ ε d v − λ W ext ( δ u ) = 0 S = D : ε , (5)where δ ε = ε l ( δ u ) + ε nl ( u , δ u ). The nonlinear problem in Eqn. (5) is usually solved by means of so-called incremental iterativemethods such as the Newton-Raphson algorithm. The present work is based on an originalalternative known as the asymptotic numerical method (ANM) proposed by Damil and Potier-Ferry (1990) and Cochelin (1994).
It is first proposed to rewrite the nonlinear problem in Eqn. (5) into the following convenientquadratic form: R ( Y , λ ) = L ( Y ) + Q ( Y , Y ) − λ F = 0 , (6)where R is a vector of residuals, L is a linear operator, Q is a bilinear quadratic operator, F isa known vector and Y T = ( u T , S T ) groups the unknowns of the problem.A key idea of the ANM then consists in expanding the unknowns Y and λ over a uniquepath parameter denoted by a in the form of the following polynomial series expansions: Y ( a ) = Y + a Y + a Y + . . . + a N Y N λ ( a ) = λ + a λ + a λ + . . . + a N λ N , (7)where ( Y , λ ) describes the initial state of the system, supposedly known. In this study, thepolynomial expansions are truncated after N = 30 terms.Introducing these expansions into Eqn. (6) and grouping the terms with the same power of a then yield the following succession of linear systems for orders p = 1 , . . . , N : L t ( Y ) = λ F L t ( Y ) = λ F − Q ( Y , Y )... L t ( Y N ) = λ N F − N − (cid:80) p =1 Q ( Y p , Y N − p ) , (8) here L t ( • ) = L ( • ) + 2 Q ( Y , • ) is the tangent operator, which is the same at all orders.At this stage, the problem involves one more unknown than the number of available equations,namely the parameter a . Similarly to a classical incremental iterative method, the ANM uses apseudo arc-length technique by setting: a = ( Y − Y ) T Y + ( λ − λ ) λ , (9)which completes the system of equations in Eqn. (8).Hence, it can be seen that the initial nonlinear problem in Eqn. (6) has been genuinelytransformed into a set of N linear systems by rejecting all nonlinearities to the right-hand side ofEqn. (8). In addition, the N linear systems composing Eqn. (8) feature a single linear operator L t which is the same at all orders. When switching to the discrete form of the problem (by meansof a classical FE displacement formulation), the resolution of the N linear systems requires onlyone decomposition of the tangent stiffness matrix K t which is the discrete counterpart of L t .This latter remark makes the ANM very efficient as the tangent stiffness matrix K t is large inpractice.Eventually, the ANM provides a continuous representation of the equilibrium path for anyvalue of λ thanks to the series expansion in a . This is an interesting property with respect tothe incremental iterative methods that need to solve the problem for each value of λ . Due to the use of finite expansions in Eqn. (7), the solution becomes invalid for large values of a . Thus, it is proposed to truncate the solution below a maximum value of a denoted by a max .This maximum value is based on a study of the norm of the residual R ( a ) = R ( Y ( a ) , λ ( a )).Cochelin (1994) proved that it is reasonable to approximate this quantity by the norm of thefirst omitted term in the expansion, so that: (cid:107) R ( a ) (cid:107) ≈ (cid:13)(cid:13) a N +1 R N +1 (cid:13)(cid:13) . (10)Based on this approximation, Cochelin then came up with the following expression for a max : a max = (cid:18) (cid:15) (cid:107) F (cid:107)(cid:107) R N +1 (cid:107) (cid:19) N +1 , (11)where (cid:15) is the maximum tolerance on the norm of the residual. This tolerance is usually setequal to a small value (here 10 − ) thanks to the normalization of the residual with respect tothe right-hand side (cid:107) F (cid:107) of Eqn. (6).The description of the whole equilibrium path is therefore made piecewise by repeating theprocedure incrementally, i.e. by resetting the initial state of the system ( Y , λ ) to ( Y ( a max ) , λ ( a max )).It is worth pointing out that the ANM remains more computationally efficient than its incremen-tal iterative counterparts because it is incremental only. Indeed, incremental iterative methods eed to iterate within the increments in order to remain on the equilibrium path, which involvesexpensive decompositions of the tangent stiffness matrix during the iterative process. The determination of the limit-load carrying capacity exploits the parametric approximationof the load proportionality factor. Indeed, limit load points are characterized by an horizontaltangent on the equilibrium path thus meaning that the derivative of the load proportionalityfactor with respect to a equals zero at the critical limit load. Hence, the limit LPF is defined as λ limit = λ ( a limit ) where: a limit = min (cid:26) a ∈ [0; a max ] : d λ d a = 0 (cid:27) . (12)Thanks to the chosen polynomial series expansion for the LPF, determining the limit load simplyconsists in finding the roots of a polynomial of order ( N −
1) and retaining the lowest positiveroot that is less than a max , provided it exists. In the present work, the ANM is applied to geometric nonlinearities extended to large rotationsbased on the work of Zahrouni et al. (1999). It can be shown that the corresponding equilibriumequations under such an assumption conveniently fit in the quadratic formulation of the ANMgiven in Eqn. (6). Two additional sources of nonlinearity are explicitly accounted for within theANM.The first source of nonlinearity is due to follower forces resulting from the hydrostatic pressureexerted on the shell. Accounting for the specific effects of follower forces in a buckling analysiscould be of utmost importance for some structural components such as the single-bay of thesubmarine pressure hull with an overall geometrical imperfection studied in Section 4. From acomputational viewpoint, this additional assumption introduces a dependence of the virtual workof external forces on the LPF (see e.g. Noirfalise, 2009, pp. 81–86). This results in additionalterms appearing in the right-hand side forces of Eqn. (8), on the one hand, and a nonsymmetrictangent operator L t (nonsymmetric tangent stiffness matrix K t in the matrix formulation), onthe other hand.The second source of nonlinearity accounted for in this work is due to the assumption of anonlinear behavior of the constitutive material of the shell. A nonlinear elastic Ramberg-Osgoodconstitutive law characterized by the following stress-strain relationship is considered within theANM (Zahrouni et al., 1998): E ε = (1 + ν ) S d − (1 − ν ) P I + 32 α (cid:18) S eq σ y (cid:19) n − S d , (13) here E is the Young’s modulus, ν is the Poisson’s ratio, σ y is the yield strength, P = − S : I = − tr ( S ) is the hydrostatic pressure, S d = S + P I is the deviatoric part of the stress tensor S , S eq = (cid:113) S d : S d is the von Mises equivalent stress, α and n are the two Ramberg-Osgoodparameters.Note that plasticity is not taken into account in the present analysis. Even though, it isargued that the structure under concern here does not present any significant local unloadinguntil the collapse load of interest is reached. In such a case, nonlinear elasticity represents afairly accurate model.As additional details, it is important to mention that large rotations and shear strains in thethin-walled shell structure of interest are accounted for by means of the shell FE formulation.The present approach resorts to a three-dimensional seven-parameter shell formulation proposedby B¨uchter et al. (1994). This formulation based on the enhanced assumed strain (EAS) conceptdisables the usual locking problems featured by shell elements. However, it introduces anotherset of nonlinear compatibility equations (see e.g. Baguet, 2001, pp. 43–48). Reliability-based design optimization (RBDO) aims at obtaining an optimal design which guar-antees a chosen reliability level w.r.t. various performance criteria. More specifically a costfunction c is to be minimized by selecting an optimal set of design parameters denoted by d ∗ while fulfilling probabilistic constraints: d ∗ = arg min d ∈D d c ( d ) : b i ( d ) ≤ , i = 1 , . . . , n c P [ g l ( X ( d )) ≤ ≤ p f l , l = 1 , . . . , n p (14)In this equation D d is the design space, { b i ( d ) ≤ , i = 1 , . . . , n c } are deterministic feasibilityconstraints on the parameters (also called soft constraints such as bounds) and { P [ g l ( X ( d )) ≤ ≤ p f l , l = 1 , . . . , n p } are the reliability constraints meaning that the retained design d ∗ shouldlead to failure probabilities smaller than the requirements { p f l , l = 1 , . . . , n p } . The notation X ( d ) is used to emphasize that the input random vector modeling the uncertainty (of prescribedPDF f X ( x | d )) contains two types of variables:- variables whose mean values are the design parameters d : their variability correspond tounavoidable scattering in the manufacturing of the system;- environmental parameters (e.g. loads) and possibly material properties that are uncertainin nature. s argued in the introduction, the complex structural behaviors that are to be addressed bythe proposed approach cannot usually be dealt with by approximate reliability methods (such asFORM/SORM) due to the complex shapes of the limit-state surfaces { g l ( X ( d )) = 0 , l = 1 , . . . , n p } .In contrast, simulation methods such as Monte Carlo or more advanced approaches such as sub-set simulation (Au and Beck, 2001) are unaffordable when using FE models. This naturallyleads to introducing metamodels for both computing the failure probabilities and optimizing thesystem. In the present paper, Kriging surrogates (Santner et al., 2003) are used in an augmentedspace which enables solving the two issues in one single step. Let us consider a single performance function g ( x ) , x ∈ D x ⊆ R n . Kriging is a statisticallearning technique that comes from geostatistics and that is now used in computer experiments(Sacks et al., 1989). A Kriging metamodel is an analytical function ˜ g that is inexpensive toevaluate (w.r.t. a performance function that involves running a large FE model) and that maybe built from an experimental design denoted by X = { x i ∈ D x , i = 1 , . . . , m } and the associateperformances gathered in y = { g ( x i ) , i = 1 , . . . , m } .Kriging assumes that the function g of interest is a realization of a Gaussian process denotedby Y ( x ) , x ∈ D x which is defined as follows: Y ( x ) = f ( x ) T a + Z ( x ) (15)In this equation f ( x ) T a is the mean of the process, which is represented by a set of basisfunctions { f i , i = 1 , . . . , P } (e.g. polynomial functions) and Z ( x ) , x ∈ D x is a stationaryzero-mean Gaussian process with variance σ Y and autocorrelation function: C Y Y ( x , x (cid:48) ) = σ Y R ( x − x (cid:48) , θ ) , ( x , x (cid:48) ) ∈ D X × D X (16)In the above equation θ gathers all the parameters defining C Y Y . In practice, square exponentialmodels are generally postulated (Lophaven et al., 2002): R ( x − x (cid:48) , θ ) = exp (cid:34) n (cid:88) k =1 − (cid:18) x k − x (cid:48) k θ k (cid:19) (cid:35) (17)The Kriging estimator at a test point x ∈ D x is obtained by the best linear unbiased estimator(BLUE) which is the linear combination of the observations (i.e. the points in X ) that provides aminimal variance. By construction the Kriging estimator is a Gaussian random variable (cid:98) Y ( x ) ∼N (cid:0) µ (cid:98) Y ( x ) , σ (cid:98) Y ( x ) (cid:1) whose mean value is the metamodel of interest:˜ g ( x ) ≡ µ (cid:98) Y ( x ) = f ( x ) T ˆ a + r ( x ) T R − ( y − F ˆ a ) (18) n this equation the following notations r , R et F are used: r i ( x ) = R ( x − x i , θ ) , i = 1 , . . . , m (19) R ij = R ( x i − x j , θ ) , i, j = 1 , . . . , m (20) F ij = f j ( x i ) , i = 1 , . . . , m, j = 1 , . . . , P (21)The Kriging variance σ (cid:98) Y ( x ) whose expression can be found in Santner et al. (2003) reflectsthe precision of the estimator. It is a measure of the epistemic uncertainty that comes from thelimited information on g gathered in the observations {X , y } . A large value of σ (cid:98) Y ( x ) means thatthe metamodel cannot be trusted at point x whereas a small value guarantees a good accuracyof the metamodel at this point.Finally the parameters σ Y , a , θ are estimated by means of the maximum likelihood principleusing the observations in y (see e.g. Marrel et al., 2008; Dubourg, 2011, for more details) or bycross-validation techniques (see e.g. Bachoc, 2013). The RBDO problem may be solved through an iterative optimization algorithm in which thefailure probability associated to the current values of the design parameters d k shall be computedfor each iteration k . In a naive approach a new metamodel could be built at each iterationconsistently with the domain of variation of the current probabilistic model X ( d k ). In order toavoid such a computational burden Dubourg et al. (2011) propose to build a single global Krigingsurrogate in an augmented space that “mixes” both the (supposedly bounded) design space D d (of volume V D d ) and the randomness in X ( d ). This corresponds to defining an augmenteddistribution h ( v ): h ( v ) = (cid:90) D d f X ( x | d ) π ( d ) d d (22)where π ( d ) denotes a prior probability distribution on the design space, which is usually auniform distribution when no specific design is preferred a priori. In this case one gets: h ( v ) ∝ (cid:90) D d f X ( x | d ) d d (23)The space on which the global Kriging surrogate is built is constructed as a confidence domain D V that is sufficiently large to allow an accurate reliability estimation whatever the currentprobabilistic model X ( d k ). In the case when the variables associated with the design param-eters are independent, bounds may be computed for each single component X ( j ) ( d ( j ) k ) and theconfidence domain D V is the hyperrectangle obtained by tensor products of these ranges (seeDubourg et al., 2011; Dubourg, 2011, for more details). ( x ) = µ c Y ( x ) = Figure 2: Example of a true limit-state surface, its Kriging surrogate and the bounds of the marginof uncertainty.
As shown in the previous paragraph the Kriging surrogate shall be precise enough in order tocompute the failure probability for various values of the current design parameters d k . TheKriging variance provides a measure of “how good the surrogate is” at each point. However in areliability analysis, the quantity of interest is the limit-state surface i.e. the contour { g ( x ) = 0 } .Hence from the metamodel ˜ g = µ (cid:98) Y and the Kriging variance Dubourg et al. (2011) define amargin of uncertainty in the vicinity of this contour: M = (cid:8) x : − k σ (cid:98) Y ( x ) ≤ µ (cid:98) Y ( x ) ≤ + k σ (cid:98) Y ( x ) (cid:9) (24)where k is a “number of standard deviations” of the epistemic uncertainty. Heuristically it isa domain in which there is a high probability (e.g. 95% for k = 1 .
96) that the true limit-statesurface lies in, up to the validity of Kriging assumptions (see Figure 2).The use of this margin (which is a subset of D x ) is twofold:- Its boundaries define two domains which yield two estimates of the probability of failurefor the current domain. The closeness of these two values denoted by p − f and p + f will be anindicator of the precision of the reliability estimates and in most cases the true value of p f lies between these bounds although there is no formal proof (see also Sch¨obi et al. (2016)for more details).- A probabilistic classification function is defined by: π ( x , t ) = P (cid:104) ˆ Y ( x ) ≤ t (cid:105) = Φ (cid:18) t − µ (cid:98) Y ( x ) σ (cid:98) Y ( x ) (cid:19) (25) ote that in this equation, P [ • ] denotes the Gaussian probability measure associated withthe epistemic uncertainty of Kriging and not the aleatoric uncertainty in X (the latterbeing denoted by P [ • ]).The probability that a point x belongs to the margin of uncertainty M is readily obtained by: C ( x ) = P (cid:104) (cid:98) Y ( x ) ∈ M (cid:105) = π ( x , k σ (cid:98) Y ( x )) − π ( x , − k σ (cid:98) Y ( x )) (26)This quantity C ( x ) is interpreted as a sampling PDF whose maximal values correspond to beingin the margin of uncertainty, i.e. close to the zero-value of the limit-state function and inregions where there is a lack of accuracy in the surrogate. By using this PDF (defined up toits normalizing constant) and a Markov chain Monte Carlo sampling technique such as slicesampling (Neal, 2003), one can draw a large set of points that lie in the margin M , i.e. that arepotential candidates for enriching the experimental design X . A technique of K -means clusteringhelps reducing this number of points to a set of limited size (say K enrich = 5 to 10) which areadded to the experimental design. A new Kriging metamodel is then built from this updatedexperimental design.In the present approach, the Kriging surrogate is refined by repeating the above procedureuntil the log-ratio of the probabilities of failure obtained from the bounds of M is smaller thanan admissible value, i.e. log( p + f /p − f ) ≤ ε p f . The proposed RBDO algorithm relies upon the gradient-based Polak-He algorithm (see Polak,1997; Dubourg, 2011). Two types of convergence criteria shall be checked for:- the Kriging surrogate shall be accurate enough to evaluate the failure probability in eachiteration, as previously described;- the optimal design of the structure should be converged.As a summary of the above sections, the proposed RBDO algorithms reads as follows:1.
Initialization : Define the ranges of the design parameters and an initial value d . Definethe conditional probabilistic model X ( d ). Compute the confidence domain D V . Define aninitial experimental design X on this domain and build the initial Kriging surrogate.2. Optimization iteration k • Estimate the current failure probability p f ( d k ) and associated bounds obtained fromthe Kriging margin of uncertainty. Possibly add new points to X and refine the Krigingsurrogate until the error on p f ( d k ) is considered acceptable. Improve the current design by one step of the gradient-based optimization algorithm(i.e. d k +1 = PolakHe( d k )). This requires the evaluation of the gradient of p f w.r.t. d . • Check the convergence of the optimization.3.
Check the reliability constraints associated with the final design. A Kriging-based impor-tance sampling scheme proposed by Dubourg et al. (2013) is used for this purpose. Thislatter technique can also be used for the analytical computation of sensitivity measures thatare used in the gradient-based optimization algorithm, see (Dubourg and Sudret, 2014).
Submarine pressure hulls are mainly composed of structural components such as ring-stiffenedcylinders, cones, elliptical or spherical ends, internal diaphragms, bulkheads and deep frames. Ata diving depth I , these structures are subjected to an external hydrostatic pressure p = ρ water g I where ρ water is the sea water density (set here equal to 1 ,
000 kg/m ) and g ≈
10 m/s isthe gravitational constant. Such a loading induces a compression stress state that is mostlymembrane dominated. Buckling therefore constitutes a critical failure mode for submarines.The design practice is usually based on specific standards and design codes such as the BritishStandard 5500 (BS5500) or the more recent Eurocode 3, possibly along with finite-element-basedsimulations. It often makes use of long-term-experience-based safety factors at various designstages, which eventually implies an unknown degree of conservatism. Hence, structural reliabilitymethods reveal a promising tool for investigating the safety margins attached to the currentsubmarine design practices (see e.g. Faulkner and Das, 1990; Pegg, 1995; Groen and Kaminski,1996; Bourinet et al., 2000).Another major challenge for the designer consists in finding an optimal ratio between theweight of the resistant structure and the buoyancy of the submersible. The latter point fallsunder the reliability-based design optimization (RBDO) formulation. The work presented in thesequel is based on preliminary studies published by Dubourg et al. (2008, 2011). The present study does not consider the submarine pressure hull as a whole. It focuses insteadon a single bay reference structure consisting of a shell cylinder with a single inner T-section ringstiffener and whose length is equal to the stiffener spacing. The dimensions of this elementarystructure are shown in Figure 3. In the following, the outer cylinder is referred to as the shell z p R² p L s = 600 mm e w = 10 mm h w = 156 mm e = 24 mm R = 2,488 mm u r u w f = 120 mm e f = 24 mm Figure 3: Single bay reference structure and initial design. (A) Overall mode (B) Interframe mode (C) Frame tripping mode
Figure 4: Schematic representation of the most critical buckling patterns of a ring-stiffened shellcylinder. plating, and the web (resp. the flange) designates the vertical (resp. horizontal) part of theT-section ring stiffener. This simplified model with well-chosen boundary conditions is supposedrepresentative of the behavior of a central bay of a sufficiently long pressure hull compartment(infinite-length in the present study).The linear elastic stability analysis of this ring-stiffened shell exhibits some typical bucklingpatterns. The three most critical kinds of buckling patterns are known as overall buckling, inter-frame buckling and frame tripping and they are basically illustrated in Figure 4. Actual struc-tures exhibit some unavoidable shape imperfections due to the manufacturing process (mostlycold-bending- and welding-based) and heavy loads connected to the hull (e.g. the nuclear reac-tor). These initial imperfections may trigger buckling or premature plastic collapse at pressurefar below those corresponding to elastic buckling, even if these imperfections are of moderateamplitude due to the stringent tolerances used in fabrication.Predicting the collapse pressure for any given imperfect geometry is not straightforward hough because the structure may feature a considerable degree of interaction between the afore-mentioned buckling modes. For solving the buckling problem at hand, the designer may resortto closed-form solutions or other semi-numerical methods available in the codes of practice (e.g.the BS5500). Another alternative that is investigated here consists in using an appropriate FEmodel. Note that the present analysis is restricted to the effects of overall and interframe shape im-perfections. The collapse due to frame tripping is avoided here by imposing some conservativerules taken from BS5500 regarding the proportions of the stiffener web and flange during theoptimization. The overall (resp. interframe) radial geometric imperfection is given by: ζ n ( z, θ ) = A n cos ( n θ ) , (27) ζ m ( z, θ ) = A m cos (cid:18) π zL s (cid:19) cos ( m θ ) , (28)where n (resp. m ) is the number of circumferential waves that typically ranges from 2 to 6(resp. 10 to 20), A n (resp. A m ) denotes the amplitude of the the overall (resp. interframe)radial imperfection, and 0 ≤ θ < π , 0 ≤ z ≤ L s . In the present application, only two modesare considered: n = 2 and m = 14. These two modes correspond to the most critical bucklingpatterns of the initial design. A finer study would consist in considering a larger spectrum ofimperfections depending on the current design at each iteration of the optimization process. It is proposed here to compute the collapse pressure by means of the asymptotic numericalmethod, accounting for material and geometric nonlinearities. The steel that constitutes thepressure hull is assumed to follow a nonlinear elastic Ramberg-Osgood constitutive law as de-scribed in Section 2.3.4. Follower forces are taken into account for the hydrostatic pressure field p so that it is always exerted normally with respect to the deformed structure.Rigid body modes are eliminated in three nodes as illustrated in Figure 5(A):- in A, the three translations are set equal to zero,- in B, the translation along the z -axis is set equal to zero,- in C, the translations along the y - and z -axes are set equal to zero.The orthoradial rotations of the two circular ends of the cylinder are set equal to zero in orderto fulfill the assumption of repeated adjacent bays. As an additional hypothesis, these two endsare supposed to remain plane and normal to the z -axis during the whole loading process, i.e. thenodes of each end cross-section undergo a constant but unknown overall axial displacement. yxy A B C
DD D-D (A) Nodes for the elimination of rigid bodymodes (B) Perfect mesh
Figure 5: Finite element modeling of the ring-stiffened shell cylinder. (A) Overall imperfection (B) Interframe imperfection (C) Both imperfections
Figure 6: Ring-stiffened shell cylinder with amplified imperfections.
In addition to the hydrostatic pressure exerted on the outer cylinder, an axial membranecompressive stress of amplitude p π R is applied as indicated in Figure 3. This additional loadis due to the hydrostatic pressure exerted on both ends of the pressure hull.The structure is meshed with 1,540 B¨uchter and Ramm 8-node shell elements featuring about40,000 degrees of freedom: 70 ×
10 elements for the outer cylinder, 70 × × In the sequel the designs obtained with the ANM-based FE model are compared with thoseobtained from approximate semi-numerical solutions available in the submarine pressure hulldesign codes of practice (see Dubourg et al., 2008, for a review). These approximations suppose n general a geometrical imperfection of a given modal shape and, as a consequence, they arenot able to account for the possible interactions between buckling modes in case of multimodal(therefore more general) imperfections. The model for predicting the overall plastic collapsepressure p n pl is based here on the Bryant formula embedded in the BS5500. The one used for theinterframe plastic collapse pressure p m pl resorts to an interpolated table of numerical solutionsderived by the Krylov Shipbuilding Research Institute (KSRI). These two models assume anoverall (resp. interframe) modal imperfection of amplitude A n (resp. A m ).The final semi-numerical model yielding the plastic collapse pressure of an infinite lengthring-stiffened cylinder with both an overall and an interframe imperfections is approximated asfollows: p critical ( A n , A m ) = min ( p n pl ( A n ) , p m pl ( A m )) . (29) In this section two design philosophies are opposed. The first one resorts to the so-called worstcase approach that consists in designing for an extreme configuration specified by experts. Theother one uses a more comprehensive probabilistic model and eventually falls under the RBDOformulation.
First, the objective of the design optimization is to find the set of parameters defining the ge-ometry of the structure d = ( e, h w , e w , w f , e f ) T that minimizes the ratio between the structuralweight and the weight of the displaced water. The latter ratio reads as follows: c ( d ) = ρ steel V steel ( d ) ρ water π ( R + e/ L s , (30)where V steel is the volume of steel composing the ring-stiffened bay and ρ steel = 7 ,
650 kg/m isthe density of steel.The admissible design space is bounded with the following constraints: (i) Since the semi-numerical model lacks consideration of the frame tripping collapse mode, itis proposed to resort to the following conservative safety criteria prescribed in the BS5500: h w ≤ . (cid:115) Eσ y e w , (31) w f ≤ (cid:115) Eσ y e f . (32)These two constraints actually bound the slenderness ratios of the stiffener components. (ii) The stiffener flange should not be too large with respect to the interframe distance:445 mm ≤ L s − w f . (33) iii) The design space is bounded with the following reasonable values: p Rσ y ≤ e ≤
50 mm , (34) w f ≤ h w ≤ w f , (35)5 mm ≤ e w ≤
25 mm , (36)70 mm ≤ w f ≤
150 mm , (37)15 mm ≤ e f ≤
50 mm . (38)The first lower constraint on the hull thickness e means that the circumferential stress inthe equivalent non-stiffened cylinder should not exceed the yield strength.At last, the predictive models for the collapse pressure (namely the FE model and the semi-numerical solutions) are used for guaranteeing that collapse does not occur at some prescribedaccidental diving depth I acc . Therefore, it leads to the establishment of the following last con-straint: I acc ρ water g ≤ p critical ( d ) . (39)It is assumed that the present submarine is designed for an accidental diving depth I acc of 250 m. The worst case approach basically consists in setting all the demand (resp. capacity) variables totheir highest (resp. lowest) possible values and to find the optimal design for this worst scenario.In the present context of shell design, this resorts to (i) prescribed maximum imperfectionamplitudes and (ii) a destruction diving depth I des that is significantly larger than the accidentaldiving depth I acc .Here, the maximum overall imperfection amplitude is taken from the BS5500 recommenda-tions and is set equal to A = 5 R/ , A
14 max = L s / Arguing that the previous worst case approach introduces an unknown degree of conservatism,it is proposed to resort to a more comprehensive probabilistic model for describing the possibleconfigurations of the hull. This probabilistic model is specified in Table 1.Since no data is available, the probabilistic model for the material properties is built from therecommendations available in the JCSS probabilistic modeling code (Vrouwenvelder, 1997). Thiscode also prescribes a linear correlation between the yield strength σ y and the ultimate stress σ u in the form of a Pearson correlation coefficient ρ = 0 .
75, which is taken into account in thepresent analysis. The right-skewed probabilistic model for the amplitudes of the imperfections
Variable Distribution Mean C.o.V. E (MPa) Lognormal 200 ,
000 0 . σ y (MPa) Lognormal 390 0 . σ u (MPa) Lognormal 570 0 . e (mm) Lognormal µ e . h w (mm) Lognormal µ h w . e w (mm) Lognormal µ e w . w f (mm) Lognormal µ w f . e f (mm) Lognormal µ e f . A (mm) Lognormal
13 5 R , . A (mm) Lognormal L s . was built with an empirical coefficient of variation of 50% and the mean is such that the previousworst imperfections A and A
14 max matches the 99.5%-quantile of the present probabilisticmodel. This thus leads approximately to set the mean value equal to one third of the latterworst imperfection amplitudes as indicated in Table 1.Given this probabilistic model, the original deterministic design optimization problem istransformed into a reliability-based design problem where safety is measured by means of thefollowing failure probability: p f ( d ) = P [ p critical ( d , X ) ≤ I acc ρ water g ] , (40)where X is the random vector that collects all the random variables of the probabilistic model.The optimization is performed w.r.t. the means of the random design variables e , h w , e w , w f and e f . The single probabilistic constraint ( n p = 1) reads as follows: p f ( d ) ≤ Φ( − β ) , (41)where β = 6 in the present application (i.e. p f 0 ≤ − ). The deterministic design optimization problem underlying the worst case approach is solved hereby means of the Polak-He gradient-based optimizer. It uses the two proposed mechanical modelsfor the buckling strength of the structure, namely the semi-numerical (SN) and the ANM-basedfinite element (FE) models. he reliability-based design optimization problem underlying the probabilistic approach issolved with the metamodel-based RBDO strategy presented in Section 3. Again, two designsare computed with either of the mechanical models.Once the four optimal designs are found, a reliability analysis is performed in order to computethe safety level of the optimally designed structures at both the accidental and the destructiondiving depth using the probabilistic model of Table 1. Since the FE model is expensive toevaluate, this resorts to the metamodel-based importance sampling technique (Dubourg et al.,2013) with a 5% target coefficient of variation on the failure probability. For the less expensivesemi-numerical model, it is proposed to resort to direct subset simulation in order to computethe whole CDF of the critical pressure which yields a relationship between the failure probabilityand the diving depth in a single run for each design. The results are given in Table 2 and the corresponding designs are illustrated in Figure 7. First,it should be noticed that the FE-based design is always more cost-optimal than its SN-basedcounterpart. Actually, this confirms the initial intuition as the semi-numerical solutions involvea set of built-in safety factors that eventually lead to an important (although unknown) degreeof conservatism. In the worst case approach, the relative gain in using a FE model w.r.t. theSN-cost is only 2%, whereas it reaches 17% in the RBDO approach.It should also be noticed that the SN-based design always features a more slender stiffenerweb than the FE-based designs. This is because the SN-solution lacks an explicit considerationof the frame tripping buckling mode. This lack is such that in the deterministic worst caseapproach, the BS5500 safety constraint regarding this mode and defined in Eqn. (31) is activeat the optimal design. Indeed, in this case, the stiffener web is clearly too slender as illustratedin Figure 7(A).As expected, the worst case approach offers a significant degree of safety at the accidentaldiving depth and it even remains a little margin at the destruction diving depth although thefailure probability is much greater there ( p f ≈ − ). The probabilistic approach enables anexplicit control of the safety level at the accidental diving depth. Due to the important targetedlevel of safety ( p f < Φ( − ≈ − ), the reliability-based optimal designs are of course lessoptimal than their worst case counterparts.The relationship between the diving depth and the failure probability is illustrated in Figure 8.The subset sampling technique applied with the semi-numerical model enables a reconstruction ofthe full CDF. The metamodel-based importance sampling applied on the expensive-to-evaluatefinite-element model only yields the failure probability estimates at the two diving depths ofinterest. It can be seen from Figure 8(B) that the failure probability matches the maximum Worst case approach RBDO ( β = 6 )Method FE-based SN-based FE-based SN-based e (mm) 21.99 26.56 28.65 35.85 h w (mm) 186.01 a e w (mm) 19.47 a w f (mm) 119.57 101.22 130.62 146.18 e f (mm) 23.97 24.53 29.68 32.77Cost 0.1960 0.2004 0.2356 0.2847 β ( I acc ) 4.99 3.81 6.06 6.11 β ( I des ) 1.40 2.00 4.42 4.99 a The frame tripping constraint is active.
FE-basedSN-based (A) Worst case approach
FE-basedSN-based (B) RBDO ( β = 6) Figure 7: Comparison of the optimal designs for the imperfect infinite-length ring-stiffened shellcylinder. 24
50 200 250 300 350 400 450 500 550 600 650 700
Diving depth (m) -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 F a il u r e p r o b a b ili t y SN-basedFE-based (A) Worst case approach
150 200 250 300 350 400 450 500 550 600 650 700
Diving depth (m) -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 F a il u r e p r o b a b ili t y SN-basedFE-based (B) RBDO ( β = 6) Figure 8: Relation between the diving depth and the failure probability for the imperfect infinite-length ring-stiffened shell cylinder. tolerance set here equal to p f = Φ( − < − .Convergence of the metamodel-based RBDO strategy is obtained within 850 calls to thebuckling strength models. Note that it is of utmost importance for the FE-based applicationdue to the important numerical effort required by a single FE analysis (about 10 minutes ofCPU time on our computer). The paper recalls the principles of an adaptive Kriging method for efficiently approximating thelimit-state surfaces that appear in RBDO problems. The Kriging variance, which is a nativeestimation of the surrogate precision, enables (i) an adaptive enrichment of the experimentaldesign and (ii) the computation of approximate confidence bounds on the probabilities of fail-ure. The definition of the augmented space that sums the input uncertainties and the range ofvariation of the design parameters enables the construction and the adaptive enrichment of aunique surrogate for the whole RBDO procedure.RBDO is applied to the design of an imperfect submarine pressure hull prone to buckling. Thesafety margin associated with the current worst case design methodology has been quantifiedin the form of a failure probability. It reveals that this common practice yields a significantlevel of safety although it is not truly mastered. Second, in order to address this latter remarkit is proposed to explicitly account for the uncertainties in the optimization problem. Thiseventually falls under the so-called RBDO formulation which is commonly identified to be toocomputationally demanding for application to real industrial problems. In this context, the roposed metamodel-based RBDO strategy truly reveals interesting to come up with a solutionwithin less than a thousand runs of the FE model. Note that in the case of larger target failureprobabilities, the RBDO problem can be transformed so as to evaluate quantiles of the limit-state function at each step of the optimization, as recently shown in (Moustapha et al., 2016) inthe context of car design against crashworthiness.Finally, the method could be improved so as to address geometrical uncertainties that aremodeled by random fields. This type of problems shows a larger stochastic dimension andrequires algorithms more advanced that the one presented in this paper to fit the surrogates. Acknowledgements
The first author was funded by a CIFRE grant subsidized by the French national agency forresearch and technology (ANRT) through the company Phimeca Engineering.
References
Arbocz, J. (1982).
Buckling of Shells: Proceedings of a State-of-the-Art Colloquium, Universit¨atStuttgart, Germany, May 6–7, 1982 , Chapter The imperfection data bank, a mean to obtainrealistic buckling loads, pp. 535–567. Springer Berlin Heidelberg.Arbocz, J. and C. D. Babcock (1969). The effect of general imperfections on the buckling ofcylindrical shells.
Journal of Applied Mechanics 36 (1), 28–38.Arbocz, J. and M. W. Hilburger (2005). Toward a probabilistic preliminary design criterion forbuckling critical composite shells.
AIAA Journal 43 (8), 1823–1827.Arbocz, J. and J. M. A. M. Hol (1995). Collapse of axially compressed cylindrical shells withrandom imperfections.
Thin-Walled Structures 23 (1–4), 131–158.Argyris, J., M. Papadrakakis, and G. Stefanou (2002). Stochastic finite element analysis of shells.
Computer Methods in Applied Mechanics and Engineering 191 (4142), 4781–4804.Au, S.-K. and J. L. Beck (2001). Estimation of small failure probabilities in high dimensions bysubset simulation.
Probabilistic Engineering Mechanics 16 (4), 263–277.Bachoc, F. (2013). Cross validation and maximum likelihood estimations of hyper-parameters ofGaussian processes with model misspecification.
Computational Statistics & Data Analysis 66 ,55–69.Baguet, S. (2001).
Stabilit´e des structures minces et sensibilit´e aux imperfections par la M´ethodeAsymptotique Num´erique . Ph. D. thesis, Universit´e d’Aix-Marseille II, Marseille, France. aˇzant, Z. P. (2000). Stability of elastic, anelastic, and disintegrating structures: a conspectusof main results. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift f¨urAngewandte Mathematik und Mechanik 80 (11-12), 709–732.Bolotin, V. V. (1962).
Statistical method in the nonlinear theory of elastic shell . NASA TT F-85.Bourinet, J.-M., N. Gayton, M. Lemaire, and A. Combescure (2000). Reliability analysis ofstability of shells based on combined finite element and response surface methods. In M. Pa-padrakakis, A. Samartin, and E. O˜nate (Eds.),
Proc. 4 th International Colloquium on Com-putation of Shell & Spatial Structures (IASS-IACM 2000), Chania, Crete, Greece, June 5–7,2000 .Broggi, M. and G. I. Schu¨eller (2011). Efficient modeling of imperfections for buckling analysisof composite cylindrical shells.
Engineering Structures 33 (5), 1796–1806.BS5500 (1997). Unfired fusion-welded pressure vessels. British Standard Institutions.B¨uchter, N., E. Ramm, and D. Roehl (1994). Three-dimensional extension of non-linear shell for-mulation based on the enhanced assumed strain concept.
International Journal for NumericalMethods in Engineering 37 (15), 2551–2568.Chryssanthopoulos, M. K. and C. Poggi (1995). Stochastic imperfection modelling in shellbuckling studies.
Thin-Walled Structures 23 (1–4), 179–200.Cochelin, B. (1994). A path-following technique via an asymptotic-numerical method.
Computers& Structures 53 (5), 1181–1192.Craig, K. J. and W. J. Roux (2008). On the investigation of shell buckling due to randomgeometrical imperfections implemented using karhunenlove expansions.
International Journalfor Numerical Methods in Engineering 73 (12), 1715–1726.Craig, K. J. and N. Stander (2007). Optimization of shell buckling incorporating karhunen-lo`eve-based geometrical imperfections.
Structural and Multidisciplinary Optimization 37 (2),185–194.Damil, N. and M. Potier-Ferry (1990). A new method to compute perturbed bifurcations: Ap-plication to the buckling of imperfect elastic structures.
International Journal of EngineeringScience 28 (9), 943–957.Dubourg, V. (2011).
Adaptive surrogate models for reliability analysis and reliability-based designoptimization . Ph. D. thesis, Universit´e Blaise Pascal, Clermont Ferrand, France. ubourg, V., J.-M. Bourinet, B. Sudret, and M. Cazuguel (2011). Reliability-based design op-timization of an imperfect submarine pressure hull. In M. H. Faber, J. K¨ohler, and K. Nishi-jima (Eds.), Proc. 11 th International Conference on Applications of Statistics and Probabilityin Civil Engineering (ICASP 11), Z¨urich, Switzerland, August 1–4, 2011 , pp. 703–711. CRCPress.Dubourg, V., F. Deheeger, and B. Sudret (2013). Metamodel-based importance sampling forstructural reliability analysis.
Probabilistic Engineering Mechanics 33 (0), 47–57.Dubourg, V., C. Noirfalise, and J.-M. Bourinet (2008). Reliability-based design optimization:an application to the buckling of imperfect shells. In P. K. Das (Ed.),
Proc. 4 th InternationalASRANet Colloquium, Athens, Greece, June 25–27, 2008 .Dubourg, V., C. Noirfalise, J.-M. Bourinet, and M. Fogli (2009). FE-based reliability analysis ofthe buckling of shells with random shape, material and thickness imperfections. In H. Furuta,D. M. Frangopol, and M. Shinozuka (Eds.),
Proc. 10 th International Conference on StructuralSafety and Reliability (ICOSSAR 2009), Osaka, Japan, September 13–17, 2009 . CRC Press.Dubourg, V. and B. Sudret (2014). Meta-model-based importance sampling for reliability sen-sitivity analysis.
Structural Safety 49 , 27–36. Special Issue In Honor of Prof. Wilson H.Tang.Dubourg, V., B. Sudret, and J.-M. Bourinet (2011). Reliability-based design optimization usingkriging surrogates and subset simulation.
Structural and Multidisciplinary Optimization 44 (5),673–690.Elishakoff, I. (1979). Simulation of space-random fields for solution of stochastic boundary valueproblems.
The Journal of the Acoustical Society of America 65 (2), 399–403.Elishakoff, I. (2012). Probabilistic resolution of the twentieth century conundrum in elasticstability.
Thin-Walled Structures 59 , 35–57.Elishakoff, I. and J. Arbocz (1982). Reliability of axially compressed cylindrical shells withrandom axisymmetric imperfections.
International Journal of Solids and Structures 18 (7),563–585.Elishakoff, I., S. van Manen, P. G. Vermeulen, and J. Arbocz (1987). First-order second-momentanalysis of the buckling of shells with random imperfections.
AIAA Journal 25 (8), 1113–1117.Faulkner, D. and P. Das (1990). Application of reliability theory to structural design and as-sessment of submarines and other externally pressurized cylindrical structures. In
Proc. 4thInternational Symposium on Integrity of Offshore Structures, Glasgow, UK , pp. 199–230. roen, H. and M. Kaminski (1996). Optimisation of pressure vessels under reliability constraints.In Proc. 15th International Conference on Ocean, Offshore & Arctic Engineering (OMAE’96) ,pp. 177–185.Koiter, W. (1945).
On the stability of elastic equilibrium (in Dutch) . Ph. D. thesis, DelftUniversity.Kriegesmann, B., R. Rolfes, C. H¨uhne, and A. Kling (2011). Fast probabilistic design procedurefor axially compressed composite cylinders.
Composite Structures 93 (12), 3140–3149.Kriegesmann, B., R. Rolfes, C. H¨uhne, J. Tessmer, and J. Arbocz (2010). Probabilistic designof axially compressed composite cylinders with geometric and loading imperfections.
Interna-tional Journal of Structural Stability and Dynamics 10 (4), 623–644.KSRI (1998). Etablissement de la m´ethode de calcul de la stabilit´e des coques. Technical report,Krylov Shipbuilding Research Institute.Lagaros, N. D. and V. Papadopoulos (2006). Optimum design of shell structures with randomgeometric, material and thickness imperfections.
International Journal of Solids and Struc-tures 43 (22–23), 6948–6964.Lophaven, S., H. Nielsen, and J. Søndergaard (2002). DACE, A Matlab Kriging Toolbox. Tech-nical report.Marrel, A., B. Iooss, F. Van Dorpe, and E. Volkova (2008). An efficient methodology for mod-eling complex computer codes with Gaussian processes.
Computational Statistics & DataAnalysis 52 , 4731–4744.Moustapha, M., B. Sudret, J.-M. Bourinet, and B. Guillaume (2016). Quantile-based opti-mization under uncertainties using adaptive kriging surrogate models.
Structural and Multi-disciplinary Optimization . Special issue on Physical, Model, and Statistical Uncertainty inStructural and Multidisciplinary Optimization (Accepted for publication).Neal, R. (2003). Slice sampling.
Annals of Statistics 31 , 705–767.Noirfalise, C. (2009).
Analyse fiabiliste de la stabilit´e des coques minces avec imperfections par lam´ethode asymptotique num´erique . Ph. D. thesis, Universit´e Blaise Pascal, Clermont Ferrand,France.Papadopoulos, V. and P. Iglesis (2007). The effect of non-uniformity of axial loading on thebuckling behaviour of shells with random imperfections.
International Journal of Solids andStructures 44 (18–19), 6299–6317. apadopoulos, V. and M. Papadrakakis (2005). The effect of material and thickness variabilityon the buckling load of shells with random initial imperfections. Computer Methods in AppliedMechanics and Engineering 194 (1216), 1405 – 1426. Special Issue on Computational Methodsin Stochastic Mechanics and Reliability Analysis.Papadopoulos, V., G. Stefanou, and M. Papadrakakis (2009). Buckling analysis of imperfectshells with stochastic non-Gaussian material and thickness properties.
International Journalof Solids and Structures 46 (14–15), 2800–2808.Pegg, N. G. (1995). The application of structural reliability methods to submarine pressurevessels. Technical memorandum DREA-TM-95-223, Defence Research Establishment Atlantic,Dartmouth NS, Canada.Polak, E. (1997).
Optimization algorithms and consistent approximations . Springer.Ramm, E. and W. A. Wall (2004). Shell structures - a sensitive interrelation between physicsand numerics.
International Journal for Numerical Methods in Engineering 60 (1), 381–427.Sacks, J., W. Welch, T. Mitchell, and H. Wynn (1989). Design and analysis of computer exper-iments.
Statistical Science 4 , 409–435.Santner, T., B. Williams, and W. Notz (2003).
The Design and Analysis of Computer Experi-ments . Springer, New York.Schenk, C. A. and G. I. Schu¨eller (2003). Buckling analysis of cylindrical shells with randomgeometric imperfections.
International Journal of Non-Linear Mechanics 38 (7), 1119–1132.Schenk, C. A. and G. I. Schu¨eller (2007). Buckling analysis of cylindrical shells with cutoutsincluding random boundary and geometric imperfections.
Computer Methods in Applied Me-chanics and Engineering 196 (35–36), 3424–3434.Sch¨obi, R., B. Sudret, and S. Marelli (2016). Rare event estimation using polynomial-chaoskriging.
ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: CivilEngineering . D4016002.Singer, J. and H. Abramovich (1995). The development of shell imperfection measurementtechniques.
Thin-Walled Structures 23 (1–4), 379–398.Singer, J., J. Arbocz, and C. D. Babcock (1971). Buckling of imperfect stiffened cylindricalshells under axial compression.
AIAA Journal 9 (1), 68–75. tefanou, G. and M. Papadrakakis (2004). Stochastic finite element analysis of shells withcombined random material and geometric properties. Computer Methods in Applied Mechanicsand Engineering 193 (1–2), 139–160.Tsompanakis, Y., N. Lagaros, and M. Papadrakis (Eds.) (2008).
Structural design optimizationconsidering uncertainties . Taylor & Francis.Vrouwenvelder, T. (1997). The JCSS probabilistic model code.
Structural Safety 19
Computer Methods in Applied Mechanics and Engineer-ing 175 (1–2), 71–85.Zahrouni, H., M. Potier-Ferry, H. Elasmar, and N. Damil (1998). Asymptotic numerical methodfor nonlinear constitutive laws.
Revue Europ´eenne des ´El´ements Finis 7 (7), 841–869.(7), 841–869.