International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP) (12th : 2015)

Use of kriging to surrogate finite element models of bonded double cantilever beams Sessa, Salvatore; Valoroso, Nunziante 2015-07

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata


53032-Paper_439_Sessa.pdf [ 296.15kB ]
JSON: 53032-1.0076210.json
JSON-LD: 53032-1.0076210-ld.json
RDF/XML (Pretty): 53032-1.0076210-rdf.xml
RDF/JSON: 53032-1.0076210-rdf.json
Turtle: 53032-1.0076210-turtle.txt
N-Triples: 53032-1.0076210-rdf-ntriples.txt
Original Record: 53032-1.0076210-source.json
Full Text

Full Text

12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015Use of Kriging to Surrogate Finite Element Models of BondedDouble Cantilever BeamsSalvatore SessaResearch Fellow, Department of Structures for Engineering and Architecture, Universityof Naples Federico II, Naples, ItalyNunziante ValorosoAssociate Professor, Department of Engineering, University of Naples Parthenope,Naples, ItalyABSTRACT: An algorithm based on kriging statistical interpolation for computing the surrogate re-sponse of a Finite Element model is presented. The interpolation model is calibrated via computationof Finite Element responses at a set of random occurrences of a material parameter. Such random gen-eration concentrates at locations where the numerical model requires a higher amount of data to get thedesired accuracy. As a model problem a standard fracture propagation test is analyzed. The proposedprocedure is shown to be robust and accurate since responses obtained via a direct computation and useof the surrogate model turn out to be undistinguishable.1. INTRODUCTIONWe discuss a kriging-based surrogate model ap-plied to the inverse identification of mode-I fractureparameters on a bonded Double Cantilever Beam(DCB) specimen.Finite Element(FE) based inverse identificationof fracture parameters from a DCB test has beenrecently assessed by Valoroso et al. (2013); the ba-sic idea developed therein amounts to simulatingthe test via a numerical model depending on a setof (unknown) material parameters whose value iscomputed via minimization of a least-squares resid-ual between the computed response and the experi-mental one. Objective of the present work is to de-velop a kriging-based surrogate model able to com-pute the response of a finite element model for theDCB depending on varying material parameters.The theoretical basis of kriging interpolation, orGaussian process regression, has been set up byMatheron (1973) based on the research of Krige(1951). The theory has acquired increasing pop-ularity due to many successful applications devel-oped mainly in spatial analysis and geostatisticalestimation, see e.g. Cressie (1993) for a completeoverview on the subject. The basic idea of themethod consists in predicting the value of a func-tion at a given point as a weighted average of ob-served data, whose weights are defined by means ofa stochastic model related to the observations cross-covariance. The main appeal of kriging interpola-tion consists in its capability to compute fast esti-mations of the function values at unknown place-ments regardless of the complexity of the observeddata providing at the same time the estimation of aconfidence interval.The surrogate model discussed in this paper usesthe responses of a set of FE simulations as the ob-served data set of the kriging interpolation. Specif-ically, a set of random points is generated in a suit-able admissibility range of the material parametersparameters to be identified and, at each point, aFE analysis is performed in order to get the cor-responding observed response. The kriging modelbased on these observations allows evaluating es-timations of the surrogate response and a corre-sponding confidence interval. A peculiar featureof the presented algorithm consists in the gener-112th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015ation procedure of the random occurrences whichis governed by the confidence interval amplitude,whereby new observations will cluster at those lo-cations where the surrogate model is less precise.It is worth emphasizing that the phase of datageneration remains computationally costly sinceeach point of the data set requires a full directFE analysis. Nevertheless, in case of experimentaltest campaigns, the proposed procedure provides isexpected to provide improved performances com-pared to the procedure presented in Valoroso et al.(2013) since specimens and protocols are usuallyfixed. Hence, the surrogate model should be gen-erated only once and its redundant employmentamounts to a set of linear computations, wherebyit is faster than direct computations.The outline of the paper is as follows. The testprotocol and the relevant FE analysis is summa-rized in section 2. Section 3 describes the krig-ing interpolation procedure. The random gener-ation algorithm and the calibration of the krigingmodel are presented in section 4; herein practi-cal applications and numerical results are presentedthat show the capabilities of the procedure. Specif-ically, the model is employed for the computationof the response surrogate of cohesive bonding in-terface models presented in Fedele et al. (2012) andValoroso et al. (2013) belonging to a Mode–I frac-ture test campaign. Finally, section 5 summarizesadvantages and drawbacks of the proposed algo-rithm and future research directions.2. MODE-I FRACTURE TEST AND FI-NITE ELEMENT ANALYSISThe symmetric DCB test depicted in Figure 1 isthe standard test employed for obtaining the mode-I fracture energy GIc of bonded assemblies. Thespecimen analyzed by Valoroso et al. (2013) con-sists of two adherends made of aluminum alloyAl 2024-T351 with length l = 200mm, width b =20mm and thickness h = 8mm; each adherend isconnected to stainless steel load block.Tests were carried out on an electromechanicalmaterial testing system shown in figure Figure 2following the guidelines of standard ISO 25217(2009); load and displacements data were recordedduring the test.Figure 1: DCB specimen geometry.Figure 2: DCB delamination test configuration.Figure 3: DCB specimen. Geometry and FE mesh usedfor identification.The FE mesh used in computations is shown inFigure 3; it consists of 3580 4-node Enhanced As-sumed Strain elements for the bulk material and125 2-node interface elements. The left end of themodel is free while the right part presents two sim-212th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015ple supports at the centroids of the load blocks withan increasing vertical displacement prescribed tothe upper block, see e.g. Valoroso et al. (2013) fora full detail.0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6020406080100120140160180δ (mm) P (N)PresentISOExperimentFigure 4: DCB specimen. Experimental VS computedload-deflection curves.Measured and computed load-deflection curvesare given in Figure 4. Here on the x-axis is re-ported the displacement δ [mm] of the upper con-straint and on the y-axis the corresponding loadP [N]. The experimental load–deflection is plottedin blue; the green and the red curves are two nu-merical responses computed using different valuesof GIc, which is the target parameter of the presentidentification exercise.3. KRIGING INTERPOLATIONKriging interpolation performs a weighted averageof the known values of a function in the neighbor-hood of the point of interest; it is a form of Bayesianinference of observed, or measured, data since av-erage weights are obtained from a Gaussian regres-sion of data. The function of interest is modeledas a Gaussian process which is defined by a priornormal distribution; then, a set of values associatedwith a spatial location is observed and the covari-ance of these evidences is computed. This defines,for each of the observations, a Gaussian likelihoodfunction which is combined with the prior one inorder to get an updated, posterior Gaussian pro-cess. The recondite action of the weighted average,concealed in the definition of the weights, consistsin computing the expected value of the posteriorGaussian process at the point of interest.Let r = [r1 . . .rn]T a vector containing the n ob-served values of the function of interest R(x), x =[x1 . . .xn]T the vector containing their spatial loca-tion and xˆ the point of interest in the domain of thefunction R(x). The weighted average consists in thefollowing:R(xˆ)≃ rˆ = wT r (1)where w is the weight vector and rˆ the estimatedvalue of R.3.1. Computation of the average weightsAverage weights are obtained via Bayesian updat-ing of the prior Gaussian process by inference ofthe cross-correlation of the observed data. From anoperational point of view, it is more convenient toemploy a function derived by the cross-correlationnamed the variogram, which is defined as the vari-ance of the difference between process values attwo locations (xi and x j) across realizations of astochastic process or a random field:γi j = γ(∣∣xi − x j∣∣)= Var[R(x j)−R(xi)] (2)the computation of the variogram at the n observedpoints defines the n+1×n+1 matrix Sd as:Sd =γ11 · · · γ1n 1............γn1 · · · γnn 11 · · · 1 0(3)while the variogram at the point of interest xˆ definesthe n+1 vector s0 as:s0 =[γ (|xˆ− x1|) · · · γ (|xˆ− xn|) 1]T (4)Average weights are obtained by the combinationof Sd and s0:[wλ]= S−1d s0 (5)where w is the weight vector to be used in Eq. (1)and λ is an error parameter. Finally, the estimationof R(xˆ) and its confidence interval amplitude τ are:R(xˆ)≃ rˆ = wT rτ =√wT s0+λn(6)312th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015where n is the number of observations and r is thevector containing the observed values.The computation of Eq. (3) and Eq. (4) requiresa proper definition of the variogram function. Aslong as the theoretical variogram can be computed,the described procedure is quite straightforward;nonetheless, it is frequent that the observed func-tion is not available in mathematical form. Insuch circumstances approximated or empirical var-iograms can be used, see e.g. Cressie (1993). Inthe present work the Gaussian variogram has beenused; it is defined as:γ(xi,x j)=0; xi = x jc0 +(cs − c0)(1− e−3|xi−x j|a2); xi 6= x j(7)where c0 is the nugget parameter, i.e. the jumpof the variogram function at xi − x j = 0; cs is thesill parameter, i.e. the limit of the variogram for∣∣xi − x j∣∣→ ∞ and a is the range parameter, i.e. thevalue of∣∣xi − x j∣∣ where the difference between thevariogram and its sill becomes negligible. Here-after, the nugget parameter will be set c0 = 0 sincethe FEM responses are not affected by any localrandomness.4. MODEL CALIBRATION ALGORITHMThe kriging interpolation procedure summarizedin the previous section requires a set of observeddata in order to completely define its mathemati-cal model. While in traditional applications exper-imental evidences are used, in this work the ob-served data are defined as the responses of finite el-ement analysis. This section presents an algorithmwhich aims to calibrate a kriging model for the esti-mation of the surrogate response of a finite elementDCB test simulation.Let then RFE (x) be a response of interest com-puted via a finite element analysis. Such a responsedepends on a generic parameter x with xl and xuas lower and upper bound respectively. Please notethat the procedure described in this section is gen-eralized and it can be used to describe any of theresponses which are usually computed by finite el-ement analysis (such as displacements, stresses, in-ternal forces etc). Analogously, the parameter x canbe any of the input variables of the FEM model. Itis worth to emphasize that, for simplicity purposes,the presented formulation considers a single, scalarparameter only; though, an extension to multiple–parameters cases can be easily implemented but itgoes beyond the aims of this work.In order to define the set of observed data, the al-gorithm randomly generates occurrences of the pa-rameter x following a Cumulative Probability Func-tion (CDF). At the stage of the initialization of theprocedure, the uniform CDF F0X (x) = x− xl/xu − xlis considered so that a set of m occurrences ofthe parameter is generated and, for each one ofthem, the finite element analysis computes the cor-responding observed value RFE :rFEi = RFE (xi) ; i = 1 · · ·m (8)Such a set of observations rFEi and its location inthe parameter space xi are capable to define a sur-rogate response by kriging interpolation as long asa variogram is defined. For this purpose, the Gaus-sian variogram of Eq. (7) is considered; however,its parameters c0, cs and a need to be convenientlycalibrated taking into account the observations rFEi .This can be easily done by a minimization symplexalgorithm applied to the optimization problem:[cˆ0 cˆs aˆ]== arg min[c0,cs,a∣∣∣∣∣m∑i=1(rFEi −w(c0,cs,a)rFEj 6=i)2](9)where rFEi is the ith observation and w(c0,cs,a)rFEj 6=iis the kriging estimate of rFEi provided by the re-maining observations. In this sense, the minimizingargument is the mean square of the residual com-puted at each one of the observation points.Once that the variogram is calibrated, a first–tentative surrogate model is easily defined by com-puting S0d , s00 and w0 by Eq. (3), Eq. (4) and Eq. (5)respectively, while the surrogate response at thepoint xˆ comes from Eq. (6).The accuracy of the surrogate model depends onthe number of observations m and on their distribu-tion in the parameter domain: a larger value of mwould provide better results, however, the random412th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015generation above does not consider any informa-tion about the features of the actual response. It isconvenient to perform the generation of the randomlocations xi by the following recursive procedure:after the first m generations, Eq. (6) provides, forany point of the domain of x, the surrogate responseand its confidence interval τ (x) which is obviouslylarger where the surrogate response is less accurate.The function τ (x) can be used as the probabilitydensity function of the random generation, if conve-niently normalized. Specifically, an updated CDFcan be defined as:F1X (x) =∫ xxlτ0 (ξ ) dξ∫ xuxlτ0 (ξ ) dξ (10)so that the random generation will provide out-comes with higher probability where the surrogatemodel lacks of precision. Note that the denomina-tor of Eq. (10) is a measure of the accuracy of thewhole surrogate model since it is its cumulative er-ror. This allows one to iterate such procedure withsubsequent generations of m observation points aslong as the overall error is greater than a tolerancevalue.The iterative procedure for generating the sur-rogate model can be summarized in the followingsteps:1. generation m new values of xi with CDF Fk−1X ;2. computation of rFEi by Finite Element analy-sis;3. calibration of the variogram γk and identifica-tion of the parameters cˆk0 cˆks aˆk through theoptimization problem:[cˆk0 cˆks aˆk]= arg min [c0,cs,a |∣∣∣∣∣(k+1)m∑i=1(rFEi −w(c0,cs,a)rFEj 6=i)2](11)4. computation of Skd , sk0 and wk by Eq. (3),Eq. (4) and Eq. (5);5. evaluation of the overall error εk =∫ xuxlτ0 (ξ ) dξ and convergence check;6. updating of the CDF:FkX (x) =∫ xxlτk (ξ ) dξ∫ xuxlτk (ξ ) dξ (12)while convergence is not reached, then the proce-dure is iterated updating k = k+1 and returning tostep 1.4.1. Numerical ApplicationIn order to provide a clearer idea about the proposedalgorithm, a numerical application is presented inthis section. It concerns the finite element simu-lation of a DCB test whose specifications are pro-vided in Section 2, with maximum displacementof the upper load block δ = 2.4mm. The criticalfracture energy release rate GIc has been used asvarying parameter x while the total energy release(i.e. the integral of the constraint reaction amongthe whole test) is the response of interest.Considering the experimental results presentedin Valoroso et al. (2013), the upper and lowerbounds of GIc have been set to 75 and 750 J/m2.Finally, the number of generated random observa-tions at each iteration of the algorithm has been setat m= 8. This is in order to take advantage of paral-lel computing: the algorithm has been implementedso that several FEM analysis run separately on dif-ferent cores and a commitment of eight cores is theoptimal choice for the employed hardware.0 100 200 300 400 500 600 700 8001015202530354045GIc [J/m2]RFigure 5: Surrogate response (solid line) and observa-tions (diamonds) at the 1st iteration.Figure 5 shows the surrogate response at the firstiterations with the parameter GIc reported on thex–axis and the surrogate response R(GIc) on the y–axis. The blue diamonds represent the responsescomputed by the FEM algorithm at each one of the512th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 20150 100 200 300 400 500 600 700 80000. [J/m2]τFigure 6: Normalized amplitude of the confidence in-terval after the 1st iteration (solid curve) and randomvalues of the parameter generated at the 2nd iteration(blue dashed lines).first eight randomly generated values of the param-eter. The surrogate response has been subsequentlycomputed by the kriging procedure described aboveand it is plotted as the red solid line. The corre-sponging amplitude of the confidence interval τ isplotted as a black line in Figure 6 where the bluevertical lines correspond to the eight values of theparameter generated at iteration 2. It is worth toemphasize how the new random observations clus-ter where the error results higher. After five iter-ations (40 FEM analysis) the algorithm provides asufficiently accurate surrogate model shown in Fig-ure 7.0 10h 20h 30h 40h 50h 60h 70h 80h10152025303540455055GIcRFigure 7: Surrogate response (solid line) and observa-tions (diamonds) at the 5th iteration.Once that the observation set corresponding tothe total energy release is defined, the calibrationof a kriging model can be performed for any otheroutput of the FEM analysis.The peculiar surrogate model has been also com-puted for the load–displacement response of theconsidered DCB specimen; specifically, each time–step of the DCB test simulation can be treated asa single, scalar response so that a kriging surrogatemodel can be defined. The whole set of surrogatemodels corresponding to each one of the responses–in–time leads to the definition of a complete load–displacement curve surrogate.Figure 8 shows a comparison between the out-put of the FEM analysis (black curves) and the cor-responding surrogate response (colored diamonds)for three values of the parameter GIc. A good fit ofthe surrogate model can easily be appreciated sincethe corresponding plots are not distinguishable; themaximum error of the surrogate model, for this spe-cific case, results about 0.1%.0 0.5 1 1.5 2 2.5050100150200250300δ [mm]P [N]  GIc=0.1 KrigGIc=0.1 FemGIc=0.23 KrigGIc=0.23 FemGIc=0.75 KrigGIc=0.75 FemFigure 8: Comparison between surrogate responses(colored curves) and FEM analysis (black curves).5. DISCUSSION AND CONCLUSIONSA surrogate model of finite element simulationsbased on a kriging interpolation procedure is pre-sented in this work. kriging interpolates observeddata of a function of interest in order to predict the612th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12Vancouver, Canada, July 12-15, 2015function’s value at an assigned point in the param-eter domain. The presented procedure uses as ob-servations the outcomes of finite element analysisperformed at randomly generated values of the pa-rameter. The peculiar capability of kriging inter-polation to provide a confidence interval leads toan efficient random generation of the observationsso that the generated points cluster where the surro-gate model lacks of precision. The numerical appli-cation shown in Section 4.1 shows the consistencyof the presented algorithm: the maximum error ofthe surrogate model results about 0.1% of the cor-responding finite element response.The presented algorithm is not strictly limited toapplications similar to the one shown in Section 4.1.Since the procedure processes numerical responses,it is possible to perform it for any kind of numer-ical analysis or structural models. Even the lim-itation about the parameter domain can be easilyovertaken: from a mathematical point of view, aconvenient definition of the CDFs could generateunbounded random values of the parameter. Up-per and lower bounds introduced in Section 4 areaimed at obtaining an easy implementation of thealgorithm.An interesting consideration concerns the com-putational times: the calibration of the krigingmodel required about 40 finite element analysis,thus, the pre–processing stage turns out quite com-putationally demanding. Although, once that thekriging model is calibrated, the presented numeri-cal application requires about 120 sec. for a com-plete finite element analysis while the surrogatemodel runs in about 2 seconds. This enormousadvantage makes such method appealing wheneverseveral analysis are required, as in case of inverseanalysis in experimental test campaigns. In general,it results convenient as long as redundant analysisare required.The algorithm high speed is due to the fact that,once the kriging model is calibrated, the surro-gate response prediction consists in linear opera-tions with matrices and vectors. This aspect con-cerns a straightforward improvement of the pro-cedure which has been omitted for lack of space.Specifically, the algorithm can be extended with aDirect Differentiation routine in order to get the sur-rogate response derivatives. This is a further ap-pealing benefit since several application cases re-quire finite differences for the computation of theirderivatives as in the case of the parameter identifi-cation presented in Valoroso et al. (2013) and con-cerning the same structural model employed in thiswork. The availability of a surrogate model provid-ing the response derivative would permit the em-ployment of gradient algorithm which are faster andmore reliable than gradient–free algorithms. Anexhaustive focus on direct differentiation would bebeyond the purposes of this work but are one of thefuture research direction.Further improvements concern the case of multi-ple parameters analysis. The procedure can be eas-ily extended in order to be performed in a multi–dimensional parameter space; nevertheless, in thiscase the peculiar relationship between the parame-ters strongly influences the definition of the auto-correlation functions, thus, specific investigationsare required. This further extension will also be atopic of the upcoming research.6. REFERENCESCressie, N. (1993). Statistics for spatial data, revisededition. John Wiley & Sons, Inc.Fedele, R., Sessa, S., and Valoroso, N. (2012). “Im-age correlation–based identification of fracture pa-rameters for structural adhesives.” TECHNISCHEMECHANIK, 32(2–5), 195–204.ISO 25217 (2009). “ISO 25217. Adhesives – Deter-mination of the mode I fracture energy of structuraladhesive joints using double cantilever beam and ta-pered double cantilever beam specimens.” Report no.,Geneva, Switzerland.Krige, d. g. (1951). A Statistical Approach to Some MineValuation and Allied Problems on the Witwatersrand.Matheron, G. (1973). “The intrinsic random functionsand their applications.” Advances in Applied Proba-bility, 5(3), pp. 439–468.Valoroso, N., Sessa, S., Lepore, M., and Cricrì, G.(2013). “Identification of mode-i cohesive parametersfor bonded interfaces based on dcb test.” EngineeringFracture Mechanics, 104, 56–59.7


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items