University of Southern California Dissertations and Theses

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

(USC Thesis Other)

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Efficient Inverse Analysis with Dynamic and Stochastic Reductions for Large-Scale Models of Multi-Component Systems by Xiaoshu Zeng A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL ENGINEERING) May 2022 Copyright 2022 Xiaoshu Zeng Acknowledgements First and foremost, I would like to thank my advisor Prof. Roger Ghanem for his continuous guidance, support, and enlightenment on my thesis and my Ph.D. research in general, and his patience, kindness, motivation, and immense knowledge. Prof. Ghanem’s help is invaluable in both my research and career. Besides my advisor, I would like to thank the rest of my thesis committee: Prof. Bora Gencturk, Prof. Erik A Johnson, Prof. Sami F Masri, and Prof. Assad Oberai, for their encouragement, insightful comments, and helpful questions. I would also like to thank my colleagues with whom I have been closely working: Prof. Olivier Ezvan, Prof. Bora Gencturk, and Dr. Gianluca Geraci for their generous help, informative discussions, and technical support. My sincere thanks also go to Dr. Bert Debusschere, Dr. Michael S. Eldred, Dr. Gianluca Geraci, Dr. John Davis Jakeman, and Dr. Habib Najm for offering me two internship opportunities in their groups and leading me to work on diverse, exciting projects. I thank my fellow labmates and friends Ruda Zhang, Philippe Hawi, Zhiheng Wang, Panagiotis Tsilifis, Ziad Ghauch, Loujaine Mehrez, Xiaohui Tu, Kelli McCoy, Meghna Kiran, Zheming Gou, and Zhengtao Yao for the inspiring discussions and all the fun we have had in the past several years. Special thanks to Xiaoying Pan, who has given me a lot of support and encouragement. We also had numerous discussions on my research topic, and she always inspired me to continue the good work. Last but not least, I would like to thank my family: my parents Hongliang Zeng and Shenghua Zeng, for their unreserved love and support; my sisters, Xiaowei Zeng, Xiaoling Zeng, and my brother Xiaoyang Zeng, for their care and encouragement. Without them, I could not achieve anything. ii Table of Contents Acknowledgements ii List Of Tables v List Of Figures vi Abstract xiv Chapter 1: Introduction 1 1.1 Research context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Objectives, methodologies, and contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4.2 Methodologies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Chapter 2: Automatic mode selection and global reduced-order model 20 2.1 Modal analysis in structural dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Automatic selection of the dominant vibration modes . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Construction of global reduced-order basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.1 Reduced kinematics bases on scale separation . . . . . . . . . . . . . . . . . . . . . . . 29 2.3.2 Reduced kinematics without scale separation . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.3 Combined kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.4.1 Simplified nuclear fuel assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4.2 Boiling water reactor fuel assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Chapter 3: Efficient modal analysis of the fully-loaded spent nuclear fuel canister 44 3.1 Fundamental tools for structural dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.1.1 Shift-Invert Lanczos approach for eigenvalue analysis . . . . . . . . . . . . . . . . . . 44 3.1.2 Domain decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.1.3 Craig-Bampton sub-structuring technique . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Craig-Bampton implementation suitable for high modal density . . . . . . . . . . . . . . . . . 48 3.2.1 Craig-Bampton implementation for Shift-Invert Lanczos . . . . . . . . . . . . . . . . . 48 3.2.2 Craig-Bampton technique adapted to the fully-loaded spent nuclear fuel canister . . . 50 3.2.3 Error quantification of the Craig-Bampton model of the fully-loaded spent nuclear fuel canister . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3 Further reduced-order model by using dominant modes of sub-structures . . . . . . . . . . . . 54 3.4 Application on the fully-loaded spent nuclear fuel canister . . . . . . . . . . . . . . . . . . . . 58 3.4.1 Reference Craig-Bampton model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.2 Two-level nested Craig-Bampton model . . . . . . . . . . . . . . . . . . . . . . . . . . 61 iii 3.4.3 Reduced-order model based on dominant fuel assembly modes . . . . . . . . . . . . . . 64 Chapter 4: Accelerated basis adaptation method in the framework of polynomial chaos expansion 67 4.1 Basis adaptation in polynomial chaos expansion . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.1.1 Polynomial chaos expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.1.2 Basis adaptation and dimension reduction . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.2 Acceleration of the convergence of the basis adaptation method . . . . . . . . . . . . . . . . . 74 4.2.1 Correction of the mean and Gaussian coefficients . . . . . . . . . . . . . . . . . . . . . 75 4.2.2 Updating the rotation matrix by higher dimensional adaptation . . . . . . . . . . . . . 78 4.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.1 Borehole function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.3.2 Space structure model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Chapter 5: Stochastic reduction for the high-dimensional quantity of interest problems with parametric and nonparametric uncertainties 94 5.1 Karhunen-Loève expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.2 Integration of Karhunen-Loève expansion and basis adaptation . . . . . . . . . . . . . . . . . 96 5.3 Incorporation of nonparametric stochastic analysis . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3.1 Fundamentals of the nonparametric stochastic analysis . . . . . . . . . . . . . . . . . . 100 5.3.2 Incorporation of parametric and nonparametric stochastic analysis . . . . . . . . . . . 104 5.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.4.1 The fuel assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.4.2 The fully-loaded spent nuclear fuel canister model . . . . . . . . . . . . . . . . . . . . 113 Chapter 6: Efficient Bayesian inference for inverse analysis 127 6.1 Bayesian method for inverse problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.2 Efficient formulation in Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.3 Markov chain Monte Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.3.1 Classical Metropolis-Hasting algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.3.2 Block update Metropolis-Hasting algorithm . . . . . . . . . . . . . . . . . . . . . . . . 132 6.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.1 The fuel assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.4.2 The fully-loaded spent nuclear fuel canister model . . . . . . . . . . . . . . . . . . . . 141 Chapter 7: Multi-fidelity uncertainty quantification for models with dissimilar parameter- ization 154 7.1 Multi-fidelity sampling uncertainty quantification . . . . . . . . . . . . . . . . . . . . . . . . . 154 7.2 Enhancing multi-fidelity UQ via basis adaptation . . . . . . . . . . . . . . . . . . . . . . . . 156 7.2.1 Important directions as shared common space among models . . . . . . . . . . . . . . 157 7.2.2 Multi-fidelity adapted basis estimator: generalities . . . . . . . . . . . . . . . . . . . . 158 7.2.3 An illustrative example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 7.2.4 Using the adapted PCE to build bias and correlation metrics . . . . . . . . . . . . . . 162 7.2.5 Strategy 1: Construction of an unbiased MFAB estimator . . . . . . . . . . . . . . . . 169 7.2.6 Strategy 2: Construction of a biased MFAB estimator . . . . . . . . . . . . . . . . . . 171 7.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 7.3.1 Analytical test problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 7.3.2 An acoustic application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 Chapter 8: Conclusions and future prospects 193 8.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 8.2 Future prospects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 Bibliography 198 iv List Of Tables 4.1 Uncertain parameters and distributions for the borehole function. . . . . . . . . . . . . . . . 81 7.1 Uncertain parameters and distribution for the Acoustic problem. Speakers are ordered coun- terclockwise with the first speaker located on the right vertical edge of the octagon. . . . . . . 185 v List Of Figures 1.1 GE 14 boiling water reactor fuel assembly [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Basic layout [2] and cross-section [3] of the HI-STAR 100 rail transportation spent fuel cask. . 3 1.3 Possible failure modes under cask drop in horizontal orientation, SAND90-2406 [4]. . . . . . . 3 1.4 Control assembly models (PWR) within a three-dimensional slice of the global structural model (with MPC-24), EPRI 1011817 [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.5 End snapshot view at 20 ms of Control Assembly at the bottom – no displacement magnifi- cation, EPRI 1011817 [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Multiscale mode of a heterogeneous plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2 Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.29) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.3 Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.33) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.34) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Projected multiscale mode in Figure 2.1 by the reduced kinematics of exact in the stiff part and flexible part in static equilibrium (with minimum residual kinetic energy) . . . . . . . . 32 2.6 Orthogonal projection of the multiscale mode in Figure 2.1 given by the piece-wise constant approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.7 Orthogonal projection of the multiscale mode in Figure 2.1 given by the piece-wise rigid-body approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.8 Orthogonal projection of the multiscale mode in Figure 2.1 given by a degree-5 polynomial approximation over the whole structure domain . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.9 Simplifiedfuelassembly: undeformedconfiguration(left), aglobalmode(center), alocalmode as seen from two points of view (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.10 Control FRF 1 calculated with the global ROM of dimension n g = 76=5441 obtained by piece-wise constant kinematics (constant by slice) and for which = 5:6 dB . . . . . . . . . . 38 vi 2.11 Control FRF 2 calculated with the global ROM of dimension n g = 76=5441 obtained by piece-wise constant kinematics (constant by slice) and for which = 51:0 dB . . . . . . . . . . 38 2.12 Control FRF 1 calculated with the global ROM of dimension n g = 113=5441 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 4:1 dB . . . . . . . . 39 2.13 Control FRF 2 calculated with the global ROM of dimension n g = 113=5441 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 3:7 dB . . . . . . . . 39 2.14 Control FRF 1 calculated with the global ROM of dimension n g = 113=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 4:2 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.15 Control FRF 2 calculated with the global ROM of dimension n g = 113=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 3:7 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.16 Control FRF 1 calculated with the ROM of the n d = 113=5441 dominant elastic modes obtained by the automatic mode selection procedure and for which = 1:4 dB . . . . . . . . . 40 2.17 Control FRF 2 calculated with the ROM of the n d = 113=5441 dominant elastic modes obtained by the automatic mode selection procedure and for which = 1:3 dB . . . . . . . . . 40 2.18 Fuel assembly details: view of the top (top), and views of the midspan (middle) and of the bottom (bottom) without the channel thin faces . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.19 Fuel assembly: undeformed configuration with and without the channel thin faces (left) and, without the channel thin faces, a global mode (center) and a local mode as seen from two points of view (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.20 Control FRF 1 calculated with the global ROM of dimension n g = 151=6605 obtained by piece-wise constant kinematics (constant by slice) and for which = 6:1 dB . . . . . . . . . . 42 2.21 Control FRF 2 calculated with the global ROM of dimension n g = 151=6605 obtained by piece-wise constant kinematics (constant by slice) and for which = 11:6 dB . . . . . . . . . . 42 2.22 Control FRF 1 calculated with the global ROM of dimension n g = 217=6605 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 4:1 dB . . . . . . . . 42 2.23 Control FRF 2 calculated with the global ROM of dimension n g = 217=6605 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 8:2 dB . . . . . . . . 42 2.24 Control FRF 1 calculated with the global ROM of dimension n g = 681=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 4:2 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.25 Control FRF 2 calculated with the global ROM of dimension n g = 681=6605 obtained by combined kinematics (rigid body by fuel rod slice and exact elsewhere) and for which = 3:7 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.26 Control FRF 1 calculated with the ROM of the n d = 681=6605 dominant elastic modes obtained by the automatic mode selection procedure and for which = 2:2 dB . . . . . . . . . 43 vii 2.27 Control FRF 2 calculated with the ROM of the n d = 681=6605 dominant elastic modes obtained by the automatic mode selection procedure and for which = 2:6 dB . . . . . . . . . 43 3.1 FE model of frame sub-structure (canister + basket), represented without the canister thick top lid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Location of the nodes associated with N I = 1581 DOFs of interest on the exterior of the canister bottom plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3 Top of a longitudinal mid-section of the canister loaded with one fuel assembly. The structure is fixed at the top of the canister lid through 2 attachments. . . . . . . . . . . . . . . . . . . . 59 3.4 Bottom of a longitudinal mid-section of the canister loaded with one fuel assembly. The external loads (red) and dynamic response monitoring (blue) are restricted to the exterior of the canister bottom plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5 Location of the nodes associated withN I = 63 DOFs of interest on the exterior of the canister bottom plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.6 Histogram of the error ij for the 20 kHz Craig-Bampton model with respect to the exact modes, both restricted to the n d = 200 most dominant modes. Red line: global error = + 4 = 0:46 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.7 FRFcomparisonbetweenthe20kHzCraig-Bamptonmodel(red)andtheexactmodel(black), both restricted to the n d = 200 most dominant modes. Error: ij = 0:46 dB. . . . . . . . . . 61 3.8 Convergence of the Craig-Bampton model error = +4 with respect to the cutoff frequency. 62 3.9 Histogram of the 4 kHz Craig-Bampton model error ij . Global error = + 4 = 0:38 dB (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.10 FRF comparison between the 4 kHz Craig-Bampton model (red) and the 20 kHz reference Craig-Bampton (black). Error: ij = 0:38 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.11 Fuel assembly domain decomposition: all the sub-structures are represented as disconnected from the others. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.12 Histogram of the error ij of the two-level nested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner truncation. Global error = + 4 = 0:41 dB (red line). . . . . 64 3.13 FRFcomparisonbetweenthe20kHzreferenceCraig-Bampton(black)andthetwo-levelnested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner truncation (red). Error: ij = 0:41 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.14 Importance value (dB) of all the fuel assembly modes (sorted). . . . . . . . . . . . . . . . . . 64 3.15 For each fuel assembly, importance value (dB) of its modes (sorted). . . . . . . . . . . . . . . 64 3.16 Convergence of the global error = + 4 with respect to the percentage of dominant fuel assembly modes preserved. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.17 Histogram of the error ij for the proposed reduced-order model with 75% removal of the fuel assembly modes below 4 kHz frequency. Global error = + 4 = 0:45 dB. . . . . . . . . . 65 viii 3.18 FRF comparison between the 20 kHz Craig-Bampton reference (black) and the proposed reduced-order model with 75% removal of the fuel assembly modes below 4 kHz frequency (red). Error: ij = 0:45 dB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.19 Histogram of the error ij for the two-level nested Craig-Bampton model removal of the one most dominant fuel assembly mode. Global error = + 4 = 7:43 dB. . . . . . . . . . . . 66 3.20 FRFcomparisonbetweenthe20kHzCraig-Bamptonreference(black)andthetwo-levelnested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner truncation from which the one most dominant fuel assembly mode is removed (red). Error: ij = 7:43 dB. . . . . . . 66 4.1 Workflow of the sequentially optimized adaptation method. . . . . . . . . . . . . . . . . . . . 81 4.2 (a) PDFs of classical basis adaptation of dimensions from 1 to 5; and (b) PCE coefficients comparison of 4d and 5d adaptations in the space spanned byf A5 ;2J 5;p g . . . . . . . . 82 4.3 (a) PDFs of basis adaptation of dimensions up to 4 with updated coefficients (UCs); and (b) PCE coefficients comparison of 3d and 4d adaptations with UCs in the space spanned by f A4 ;2J 4;p g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.4 (a) PDFs of 1d and 2d adaptations by B (1) ; and (b) PCE coefficients comparison of 1d and 2d adaptations byB (1) in the space spanned byf A2 ;2J 2;p g . . . . . . . . . . . . . . . . 84 4.5 Comparison of absolute values of the first rows of rotationsB (1) andB (2) . . . . . . . . . . . 85 4.6 (a) PDFs of 2d and 3d adaptations by B (1) ; and (b) PCE coefficients comparison of 2d and 3d adaptations byB (2) in the space spanned byf A3 ;2J 3;p g . . . . . . . . . . . . . . . . 85 4.7 (a) PDFs obtained from Monte Carlo method with different samples and (b) a zoomed in view in the right tail region; the full dimensional PCE of order 3 is served as reference model. . . . 86 4.8 (a)Finiteelementmodelofthespacestructure; (b)Finiteelementmodelofthespacestructure without the outer shell and inner absorption block. . . . . . . . . . . . . . . . . . . . . . . . 87 4.9 (a) Time history of the impulse force and (b) time history of acceleration in X direction of the critical point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.10 PDFs of the maximum X acceleration of the critical point obtained by full dimensional PCE and classical basis adaptations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.11 (a)PDFcomparisonoffirstorderexpansionsbetweenfulldimensionalPCEand9dadaptation; (b) and Gaussian coefficient comparison between full dimension PCE and 9d adaptation in the full dimensional space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.12 (a) PDFs of the maximum X acceleration of the critical point obtained by basis adaptation of dimensions up to 3 with updated coefficients (UCs); and (b) PCE coefficients comparison of 2d and 3d adaptations with UCs in the space spanned byf A3 ;2J 3;p g. . . . . . . . . . . . 91 4.13 (a) PDFs of the maximum X acceleration of the critical point obtained by 1d and 2d adapta- tions byB (1) ; and (b) PCE coefficients comparison of 1d and 2d adaptations byB (1) in the space spanned byf A2 ;2J 2;p g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.14 Comparison of absolute values of the first rows of rotationsB (1) andB (2) for the space model. 92 ix 4.15 (a) PDFs of the maximum X acceleration of the critical point obtained by 2d and 3d adapta- tions byB (2) ; and (b) PCE coefficients comparison of 2d and 3d adaptations byB (2) in the space spanned byf A3 ;2J 3;p g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.1 Probability density function of theBe(5:0; 1:5) distribution with scale factor 1.85 kN/mm. . . 108 5.2 Mean FRF of the FA obtained from the level 1 quadrature. . . . . . . . . . . . . . . . . . . . 108 5.3 Eigenvalues of the covariance matrix of the FRF of the FA approximated by level 1 quadrature.109 5.4 The first 8 eigenvectors of the covariance matrix of the FRF of the FA approximated by level 1 quadrature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.5 Probability density functions of KL parameters obtained from PCEs of order up to 1 (in blue), up to 2 (in orange), and up to 3 (in green) for the FA. . . . . . . . . . . . . . . . . . . . . . . 110 5.6 Plot of 500 MC samples of the first 8 KL parameters for the FA with X-axis the values obtained from the physical model evaluations and Y-axis the values obtained from the PCE models; the blue lines in the figures are lines with unit slop. . . . . . . . . . . . . . . . . . . . 111 5.7 Confidence intervals of the probabilistic model of the FRF of the FA obtained from PCE via KLE; the green solid line presents the mean FRF. . . . . . . . . . . . . . . . . . . . . . . . . 112 5.8 Comparison of the FRF of FA at the frequencies around peaks; the blue lines are obtained from PCE of the FRF and the orange lines are obtained from PCE of the FRF via KLE. . . 113 5.9 The95%confidenceinterval(pinkshades)oftheprobabilisticmodeloftheFRFobtainedfrom nonparametric analysis; the reference deterministic FRF (in black) and the FRF obtained by the most dominant 14; 500 modes (in blue) are also presented. . . . . . . . . . . . . . . . . . . 115 5.10 The 95% confidence interval (pink shades) of probabilistic model of the FRF obtained from nonparametric analysis with the most 327 system modes randomized; the reference determin- istic FRF (in black) and the FRF obtained by the most dominant 14; 500 modes (in blue) are also presented. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.11 The 95% confidence interval (pink shades) of probabilistic model of the FRF obtained from nonparametricanalysiswiththemost327systemmodesrandomizedanddispersionparameter C = 0:005. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.12 Mean FRF obtained from the level 1 quadrature. . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.13 Eigenvalues of the covariance matrix of the FRF approximated by level 1 quadrature. . . . . 118 5.14 Thefirst20eigenvectorsofthecovariancematrixoftheFRFapproximatedbylevel1quadrature.118 5.15 Probability density functions of KL parameters obtained from classical 2d (in orange) and 3d (in green) adaptations, and 204-dimensional order 1 PCEs (in blue). . . . . . . . . . . . . . . 120 5.16 Probability density functions of KL parameters obtained from accelerated 2d (in orange) and 3d (in green) adaptations, and 204-dimensional order 1 PCEs (in blue). . . . . . . . . . . . . 121 5.17 PCE coefficients of the accelerated 3d adaptation models of the KL parameters. . . . . . . . . 122 x 5.18 Confidence intervals of the probabilistic model of the FRF obtained from PCE via KLE; the green solid line represents the mean FRF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.19 PDFs of the FRFs at selected frequencies obtained from classical 2d (in blue) and 3d (in orange) adaptations, and from 3d adaptation via KLE (in dashed green). . . . . . . . . . . . 124 5.20 PDFs of the FRFs at selected frequencies obtained from accelerated 2d (in blue) and 3d (in orange) adaptations, and from 3d adaptation via KLE (in dashed green). . . . . . . . . . . . 125 6.1 Comparison of the FRF obtained from the damage scenario I to the mean and confidence intervals of the prior distribution of the FRF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.2 MCMC chains of 1 and 2 for damage scenario I. . . . . . . . . . . . . . . . . . . . . . . . . 136 6.3 MCMC chains of 3 and 4 for damage scenario I. . . . . . . . . . . . . . . . . . . . . . . . . 136 6.4 Posterior distributions (in blue color) of the input parameters of the damage scenario I against their corresponding prior distributions (in red color). . . . . . . . . . . . . . . . . . . . . . . . 137 6.5 Comparison of the mean and confidence intervals of the posterior distribution of the FRF to the FRF obtained from damage scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 6.6 Comparison of the FRF obtained from the damage scenario II to the mean and confidence intervals of the prior distribution of the FRF. . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 6.7 MCMC chains of 1 and 2 for damage scenario II. . . . . . . . . . . . . . . . . . . . . . . . . 139 6.8 MCMC chains of 3 and 4 for damage scenario II. . . . . . . . . . . . . . . . . . . . . . . . . 139 6.9 Posteriordistributions(inbluecolor)oftheinputparametersofthedamagescenarioIIagainst their corresponding prior distributions (in red color). . . . . . . . . . . . . . . . . . . . . . . . 140 6.10 Comparison of the mean and confidence intervals of the posterior distribution of the FRF to the FRF obtained from damage scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.11 Canister bottom plate with locations of the fuel assemblies. . . . . . . . . . . . . . . . . . . . 142 6.12 Comparison of the FRF generated by the damage scenario I to the probability distribution of the FRF generated by the forward probabilistic model. . . . . . . . . . . . . . . . . . . . . . . 143 6.13 MCMC chain of spring constant of type I connections. . . . . . . . . . . . . . . . . . . . . . . 144 6.14 MCMC chain of spring constant of type II connections. . . . . . . . . . . . . . . . . . . . . . 144 6.15 MCMC chain of spring constant of type III connections. . . . . . . . . . . . . . . . . . . . . . 144 6.16 Posterior distributions (in blue color) of the type I connection of the damage scenario against their corresponding prior distributions (in red color) for the FLSNFC model. . . . . . . . . . 149 6.17 Posterior distributions (in blue color) of the type II connection of the damage scenario against their corresponding prior distributions (in red color) for the FLSNFC model; the posterior mean of the presenting connections has decreased from the prior mean by more than 15%. . . 150 xi 6.18 Posterior distributions (in blue color) of the type III connection of the damage scenario against their corresponding prior distributions (in red color) for the FLSNFC model; the posterior mean of the presenting connections has decreased from the prior mean by more than 15%. . . 151 6.19 Posterior distribution of the type I connection associated with damaged FAs obtained from posterior samples from iterations between 7; 500 and 12; 500 (in blue color), and between 12; 500 and 17; 500 (in green color); the associated prior distribution is also presented (in red color). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 6.20 Posterior distribution of the type II connection associated with damaged FAs obtained from posterior samples from iterations between 7; 500 and 12; 500 (in blue color), and between 12; 500 and 17; 500 (in green color); the associated prior distribution is also presented (in red color). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 6.21 Posterior distribution of the type III connection associated with damaged FAs obtained from posterior samples from iterations between 7; 500 and 12; 500 (in blue color), and between 12; 500 and 17; 500 (in green color); the associated prior distribution is also presented (in red color). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.22 Comparison of the confidence intervals of the reconstructed FRF and the FRF obtained from the damaged model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 7.1 Scatter plot of f and g in original coordinates ( 2 = 0:08), where the green line indicates the perfect linear correlation between the models. . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 7.2 (a) Scatter plot of f and g and (b) the variation of the functions f and g along the adapted dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 7.3 Normalized histogram for 1000 realizations of several estimators for the expected value . . . . 161 7.4 Histogram of different estimators; HF and LF are 1d and 1d adapted in MFAB . . . . . . . . 162 7.5 Histogram of different estimators; HF and LF are 2d and 1d adapted in MFAB . . . . . . . . 162 7.6 Scatter plot for f and g in the original coordinates. . . . . . . . . . . . . . . . . . . . . . . . . 174 7.7 (a) Scatter plot in 1d and 1d adapted spaces, and (b) responses with respect to the adapted variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 7.8 MSE measure for the exponential problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 7.9 Bias measure based on first order information for the exponential problem. . . . . . . . . . . 176 7.10 Cumulative to total different ratio of the MSE measure of the exponential problem. . . . . . . 177 7.11 Cumulative to total different ratio of the first order information measure of the exponential problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 7.12 Squared correlation between adapted HF and adapted LF computed by pilot samples. . . . . 177 7.13 Squared correlation between projected HF and adapted LF for the unbiased strategy. . . . . . 178 7.14 Scatter plot of f and g with adapted dimensions being (a) 1 and 1, and (b) 4 and 4; the unbiased sampling strategy is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 xii 7.15 Probability density functions for 500 realizations of several expected value estimations where the unbiased strategy is used in the MFAB estimator. Dimensions of HF and LF models in MFAB are 4 and 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 7.16 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB and MFAS are 1 and 1. . . . . . . . . . . . . . . . 182 7.17 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 3 and 3. . . . . . . . . . . . . . . . . . . . . . . 182 7.18 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 5 and 5. . . . . . . . . . . . . . . . . . . . . . . 183 7.19 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 6 and 5. . . . . . . . . . . . . . . . . . . . . . . 183 7.20 (a) Example pressure plot of the speaker cabinet. The QoI is the average pressure inside the red rectangle. (b) Scatter plot of HF and LF in the original coordinates of the acoustic application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 7.21 (a) Scatter plot in 1d and 1d adapted spaces, and (b) responses with respect to the adapted variable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 7.22 MSE measure for the acoustic application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 7.23 Bias measure based on first order information for the acoustic application. . . . . . . . . . . . 187 7.24 Cumulative to total difference ratio of the MSE measure for the acoustic application. . . . . . 187 7.25 Cumulative to total difference ratio of the first order information measure for the acoustic application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 7.26 Squared correlation between HF and LF of the acoustic application computed by the pilot samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 7.27 Squared correlation between HF and LF of the acoustic application computed by the unbiased strategy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 7.28 Probability density functions for 500 realizations of several expected value estimations where the unbiased strategy is used in the MFAB estimator. Dimensions of HF and LF models in MFAB are 2 and 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 7.29 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 1 and 1. . . . . . . . . . . . . . . . . . . . . . . 191 7.30 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 2 and 2. . . . . . . . . . . . . . . . . . . . . . . 191 7.31 Probability density functions for 500 realizations of several expected value estimations. Di- mensions of HF and LF models in MFAB are 3 and 2. . . . . . . . . . . . . . . . . . . . . . . 192 xiii Abstract This work has two main goals: the primary goal of dealing with inverse analysis for large-scale models of multi-component systems and the second goal of dealing with the multi-fidelity uncertainty quantification (UQ) for models with dissimilar parameterization. The primary goal involves efficient structural dynamic analysis, probabilistic modeling, and inverse analysis for the following reasons: 1. The context is the integrity assessment of the internals of a complex multi-scale structure, a fully- loaded spent nuclear fuel canister (FLSNFC), with the assessment based on dynamic signals on the exterior surface of the FLSNFC before and after transportation. 2. Since the observed data is inevitably subjected to errors, the inverse analysis usually requires an accurate probabilistic model to capture the uncertainty propagated to the quantities of interest (QoI). 3. An efficient forward model, a dynamic model, is essential to construct a probabilistic model for complex systems. A first attempt to build an efficient dynamic model is through constructing a global reduced-order basis (ROB). Usually, the complex multi-scale structures are characterized by numerous local vibration modes (or elastic modes) and the usual long-wavelength global vibration modes. Accordingly, a methodology that does not require the computation of the numerous local modes builds a global reduced-order model (ROM) by constructing a global ROB. In this method, the kinematics of the structure is modified to filter the local vibrations. Moreover, the reduced kinematics is combined with the idea of static condensation to achieve higher accuracy. Due to the high accuracy requirement, a second attempt for dynamic modeling is carried out. For the FLSNFC, a honeycomb basket is placed inside the cylindrical canister, and a fuel assembly (FA) that holds xiv nearly 100 fuel rods is inserted in each of the 68 basket cells. A multi-level nested Craig-Bampton (CB) sub-structuring method with shift-invert Lanczos (SIL) eigenvalue solver and filtering of the local vibration modes of the substructures is proposed. This method is adapted to the multi-scale nature and localized connections between the substructures. The CB sub-structuring technique takes advantage of the limited degrees of freedom (DOF) of internal boundary and is applied to modal analysis for two structural levels, the system and the FA levels. As a result, the integrated method achieves a computational gain of four orders of magnitude for the FLSNFC at the expense of negligible errors. For probabilistic modeling, polynomial chaos expansion (PCE) is an efficient method, but it suffers from the curse of dimensionality. However, a basis adaptation method proposed by Tipireddy and Ghanem (2014) [6] can reduce the dimension of the problem by rotating the input Gaussian random variables such that the quantity of interest (QoI) in the new space has concentrated representation. In this study, we proposed two novel approaches that can accelerate the convergence of the basis adaptation method. In the first approach, the mean and Gaussian coefficients in the adapted space are corrected by information obtained from a pilot PCE. The second approach updates the rotation matrix by taking advantage of the probabilistic information embedded in the higher dimensional adaptation gleaned from an initial adaptation. As a result, both approaches achieve accelerated convergence of the basis adaptation method with negligible additional costs. The basis adaptation is adequate to reduce the dimension for the scalar QoI problems. To deal with UQ problems with high dimensional QoI and high dimensional parameter space, the integration of Karhunen- Loève expansion (KLE) and basis adaptation is proposed. The KLE first approximates the QoI to reduce the dimension of the QoI to a limited number of KL terms. Then, for each KL term, an adapted PCE is built with the accelerated basis adaptation method to reduce the dimension of the input variables. The PCE models of the KL terms then can be substituted back to the KLE of the QoI to obtain a surrogate probabilistic model of the QoI. With the accuracy of the surrogate model verified, this parametric approach is combined with the nonparametric approach to account for uncertainty induced by modeling errors. The Bayesian method will be used for the inverse analysis of parameter inference given observed data of a possibly damaged model. The process usually involves evaluating the forward model numerous times, xv which is intractable for the complex system considered in the present study, even if the dynamic ROM is used. Thanks to the accurate surrogate probabilistic model of the QoI, the physical model required to be evaluated in Bayesian analysis can be replaced by the surrogate model to achieve several orders of computational efficiency. Nevertheless, generating posterior samples of high dimensional parameters can still be challenging. Thus, a block-update Markov Chain Monte Carlo (MCMC) method is applied to address this issue. By appropriately designing the inverse problem, the location and damage types of the internals can be identified based on dynamic signals on the exterior surface. Thesecondgoaloftheworkistoproposeanovelmulti-fidelityUQmethodfordissimilarparameterization models. In multi-fidelity UQ, credible prediction and analyses of high-fidelity (HF) models are obtained by leveraging evaluations of a large number of efficient low-fidelity (LF) models. The efficacy of the technique relies heavily on the correlation of the HF and LF models. We propose using the basis adaption method in the multi-fidelity technique to independently identify the important directions (or adapted variables) for each model, andtheimportantdirectionsassembleacommonlower-dimensionalspace. Sinceimportantdirections have concentrated information of the QoI and are aligned for different models, the samples generated on the common lower-dimensional manifold have enhanced correlations. Thus, the proposed method can increase the performance of the multi-fidelity technique, especially for models with dissimilar parameterization. The two main goals in the thesis are relevant since both involve stochastic reduction for complex multi- component models. xvi Chapter 1 Introduction 1.1 Research context This work has two main goals, with the primary goal dealing with parameter identification that belongs to the scope of inverse analysis for large-scale models of multi-component systems. To achieve the primary goal, dynamic and stochastic reductions of the models are required. The intention is to propose a probabilistically- informed methodology that involves using the non-destructive evaluation (NDE) techniques to determine the extent of potential damage or degradation of the internals of a fully-loaded spent nuclear fuel canister (FLSNFC) before and after the normal condition of transportation and hypothetically accident condition. The assessment is based on dynamic signals on the exterior surface of the FLSNFC, which could be high- dimensional. Thus, several analyses have to be carried out: dynamic analysis adapted to the large-scale model of multi-component systems, an efficient probabilistic analysis model for high-dimensional input and output problems, and inverse analysis for expensive forward models with high-dimensional parameter space. These analyses deal with models of high fidelity (HF). However, when low fidelity (LF) models are available and the quantity of interest is some statistics of the HF model, the multi-fidelity uncertainty quantification (MFUQ) techniques can be applied. Therefore, the second goal of the work is to deal with the MFUQ problem for models with dissimilar parameterization. The dissimilar parameterization issue is common yet challenging in the MFUQ field. The spent fuel [7] is the nuclear fuel that has undergone fission in the reactor and is no longer serviceable in the present form due to the depletion of fissile material, poison buildup, or radiation damage. The spent 1 fuel removed from the reactor is boiling and radioactive, and hence, storage is crucial for public safety. There are two main methods for the storage, wet and dry, which are distinguished by the medium for the cooling process. The wet storage usually happens at the reactor sites. The dry storage is often stored off-site, in which the spent fuel container is shipped away from the reactor and stored in, for example, the drywell/tunnel for permanent disposal. Currently, dry storage methods are more popular, especially metal cask storage and concrete cask storage methods [8, 7]. Due to a potential high-level waste repository being built, the number of shipments of the spent fuel containers by road and rail is expected to increase [9]. The spent fuel cask will receive mechanical vibration during transport. As the spent fuel’s structural integrity depends on the vibrations and can be damaged, especially in accident transport conditions, the dynamic analysis of the spent fuel cask is required. The spent fuel rods are bundled by two types of fuel assemblies, the boiling water reactor (BWR) and the pressurized water reactor (PWR). The BWR, more specifically the GE 14 assembly (see Figure 1.1), is used for investigation in the current study. The HI-STAR 100 cask model is used to place the BWR fuel assemblies in this study, in which 68 fuel assemblies are placed into a honeycomb basket structure inside a canister and is a computational benchmark model for BWR [10]. The HI-STAR 100 is an all-steel cask made for rail transport of fuel in an inner welded multi-purpose canister (MPC), the layout and cross-section of which can be seen in Figure 1.2. Specifically, the MPC-68 is used to place 68 BWR fuel assemblies in this study, and the thick layers of materials outside the MPC that are built to shield the radiations are not modeled. Figure 1.1: GE 14 boiling water reactor fuel assembly [1]. The best example of the analytical evaluation of the effects of hypothetical accident conditions on the structural integrity of the fuel assemblies and fuel rods is the well-known study by Sandia National Labora- tories [4, 5]. Their study performs detailed structural analyses and failure evaluations to quantify fuel rods and assemblies’ failure frequency and potential reconfiguration. A two-step process is proposed, namely, (i) 2 Figure 1.2: Basic layout [2] and cross-section [3] of the HI-STAR 100 rail transportation spent fuel cask. the global structural analysis of cask drop event (9-meter drop of a cask to model the accident condition) of the most dangerous orientation (side-drop) to determine the dynamic forces acting on the fuel rods and (ii) the application of global dynamic forces on the detailed fuel rods for failure analysis [5, 11]. Three possible failure modes of the fuel rods, a pinhole failure, a breakage mode failure, and an axial split failure, are depicted in Figure 1.3. SAND90-2406 demonstrated that the side drop resulted in the highest failure Figure 1.3: Possible failure modes under cask drop in horizontal orientation, SAND90-2406 [4]. probabilities for all failure modes, with 2E-4 for the pinhole failure, 5E-5 for the breakage failure, and 2E-5 for the axial split failure [4]. A global model that consists of a three-dimensional axial slice of the cask body, the honeycomb basket structure, and fuel assemblies is used to determine the dynamic forces that act on the fuel rods, see in Figure 1.4. In this example, 24 PWR fuel assemblies are placed in the MPC, different from 3 Figure 1.4: Control assembly models (PWR) within a three-dimensional slice of the global structural model (with MPC-24), EPRI 1011817 [5]. the configuration in this study; however, the methodology works for both cases. The axial slice includes the spacer grid plane with two half spans. One of the fuel assemblies is modeled as “Control Assembly” in detail, while the others are modeled as super beams with equivalent mass and stiffness. The fuel rods in the Control Assembly are modeled as beams. The Control Assembly is placed in different positions, top, bottom, at an angle to side drop orientation, to capture the variations of dynamic forces. After a 9-meter side-drop simulation, a snapshot end view of the deformations of the cask, basket, and Control Assembly (bottom) is presented in Figure 1.5. This figure shows the fuel rods’ rod-to-rod compaction, possibly caused by the Figure 1.5: End snapshot view at 20 ms of Control Assembly at the bottom – no displacement magnification, EPRI 1011817 [5]. connection degradation between the rods and spacer grids or the collapse of the spacer grids. In addition, the 4 connections between the fuel assemblies and the basket are loosed. Also, there is deformation in the basket plates, but it is important to note that they did not collapse. The detailed fuel rod deformation patterns and local fuel rod response due to the global dynamic forces can be seen in [5, 11]. In this analysis, the model is only an axial slice of the whole structure, and the assemblies are only detailed for the Control Assembly; in addition, the overall mesh size is not refined. Thus, several engaging scenarios are not investigated; for example, it is hard to know the geometry of the damages, the ensemble dynamical behavior of the fuel rods is not presented, and the global bending of the cask is also hard to predict. Therefore, a detailed model for the whole structure that can assess the structural integrity of the fuel assemblies and fuel rods is not available. In addition, the damages of the internal components introduce significant uncertainties of their behaviors. Various sources of uncertainty also lead to an inverse problem of high dimensions, which makes obtaining posterior distributions of the parameters intended to infer more challenging. The challenges of the inverse problem also arise from the need for repeated evaluations of the dynamic model. The second goal of the work is to obtain statistics with acceptable accuracy for the UQ problem of large- scale problems. MFUQ techniques are one of the tools to obtain credible statistics of the QoI from HF and computationally demanding models and have received significant attention in the UQ field in recent years. The main idea is to overcome the growth of computational power needed for the UQ of HF computations by leveraging numerous evaluations obtained through approximated solvers, the LF models. The main challenges of MFUQ include discovering and exploiting relationships among the HF and LF models and the optimal computational resource allocation among the models. To address these challenges is the second goal of the thesis. 1.2 Challenges Thedynamicanalysis,probabilisticmodeling,andinverseanalysisofthelarge-scalemodelofmulti-component systems involved in the primary goal of integrity assessment of the internals of the FLSNFC pose significant challenges, which some of them have already mentioned in Section 1.1. In addition, the second goal of 5 proposing credible statistics for HF models also exhibits challenges to improving computational efficiency. In this section, these challenges will be described further in more detail. The FLSNFC is sealed to block radiations. Consequently, the structural health of the internals is eval- uated based on the signatures detected on the exterior surface of the structure. The frequency response function (FRF) fully characterizes the linear dynamics of the structure through a frequency domain analysis in a broad frequency band. The complex multi-scale structures usually have numerous substructures (con- sisting of many small-scale, flexible sub-components) attached to the main body or frame structure through the internal boundaries. Examples include cars, nuclear fuel containers (like the FLSNFC), naval ships, and aerospace structures. In the FLSNFC, a honeycomb basket is placed inside the thick cylindrical canister, and a fuel assembly (FA) that holds 92 fuel rods is inserted in each of the 68 basket cells. The fuel rods (about 1-centimeter diameter) in each FA are guided and bundled by eight spacer grids along with height (about 4 meters). A clear separation of scales can be observed in the FLSNFC, in which the “lower scale” contains numerous fuel rods, the “middle scale” includes the 68 fuel assemblies, and the “upper scale” (also referred to as the structural frame) is constituted of the honeycomb basket and the thick canister. About the structural integrity of the internals, we are particularly interested in the damaged fuel rods, the loss of connection between fuel rods and spacer grids, and the loss of connections between fuel assemblies and the structural frame. The dynamics of the upper scale assess the lower and middle scales’ structural integrity. Therefore, a high-fidelity finite element (FE) [12] computational model that describes all three structural scales is built, through which it can discover the signatures transmitted from fuel rods and fuel assemblies to the exterior surface of the structural frame. Due to the multi-scale nature and the fine mesh of the model, there are approximately 140 million degrees of freedom (DOF) in the FE model. Moreover, in addition to the long-wavelength global vibration modes (or elastic modes), numerous local vibration modes associated with thousands of the relatively flexible fuel rods are presented in the broad frequency band. It is challenging for structures with high dimensionality and modal density to perform an accurate structural dynamics analysis. The possible unknown damages of the internals before and after transportation introduce uncertainties to the input space, and these uncertainties propagate to the FRF through the physical model. The damages in the internals of the structure could change the dynamic behaviors completely. However, it is usually hard to 6 know what kind of damage scenarios are occurred; for instance, in the FLSNFC, any of the fuel rods can be damaged, any of the connections between the fuel rods and spacer grids can be loosed, and the severity of the damages also varies due to different transportation conditions. The UQ analysis for risk assessment is thus performed. Due to a large number of fuel rods, enormous connections between different structural scales, and the unpredictability of various damage scenarios, the dimension of the forward UQ problem is high (on a scale of hundreds). Though the Monte Carlo (MC) [13] method for stochastic analysis is independent of the dimension of the problem, it has slow convergence. The well-known polynomial chaos expansion (PCE) [14] method is an efficient stochastic analysis method, but it suffers from the curse of dimensionality. The number of required samples to construct accurate PCE grows geometrically with dimension. Nowadays, the complexity of large-scale models people are dealing with adds more limits on the number of affordable samples for probabilistic analysis. Thus, efficiently carrying out stochastic analysis for complex systems with high dimensional random inputs remains a significant challenge at the forefront of UQ. Besides the high-dimensional input space, the FRF in a wide frequency range as the quantity of interest (QoI) is also of high dimension. Thus, the direct application of some dimension reduction technique that is adapted to scalar QoI, for example, the so-called basis adaptation method [6], is restrained. For the inverse analysis, we focus on Bayesian approaches [15, 16, 17]. The output of Bayesian inference is usually a probability distribution that contains all the available probabilistic information of the interesting parameters. However, generating samples from the posterior distribution usually requires numerous eval- uations of the forward model (the dynamic model in our case). Moreover, the computational cost might be unaffordable for a large-scale model and a high-dimensional input space inherited from the probabilistic model. The Markov Chain Monte Carlo (MCMC) algorithms [18, 19, 20] is one of the most popular tools to generate samples from multi-dimensional posterior distributions. However, the sampling can be significantly slow as the dimension of the parameters to be identified is significantly high, and the convergence rate of MCMC algorithms typically decreases with the increase of parameter dimension. Thus, proposing efficient sampling strategies is still an active research domain in inverse analysis. Though many advancements of the MFUQ technique can be found in the literature, some challenges re- main in the field that decreases the applicability [21]. The thesis focuses on the issue exhibited in dissimilar 7 or heterogeneous parameterization characteristics among the HF and LF models. In practical applications, the HF model typically has refined spatial/temporal resolutions with detailed physical models. In order to increase the efficiency of MF estimators, it is common in practice to resort to models with much-simplified physics to enable significant cost separation. However, altering the physics could result in dissimilar pa- rameterization among the models; for example, each model might have a set of parameters that are not shared with other models, or the variation of the QoI concerning the shared parameters is different among the models. The presence of dissimilar parameterization thus is common and is an obstacle for almost all the MFUQ techniques as it decreases the correlation among the HF and LF models, and as a result, a large number of HF samples are required for credible statistics. Therefore, the efficiency of the MF estimator is a function of correlation among the HF and LF models and the cost ratios of the LF models to the HF model. 1.3 Literature review A literature review of the research related to the challenges we are trying to address provides insight into where the current research stands and will be given in this section. A critical feature of the dynamics of the complex multi-scale structure is a mixture of global modes and numerous local modes. The global modes are associated with vibrations of the structural frame or stiff part of the structure, while the local modes are associated with vibrations of the flexible sub-components of the substructures with the structural frame remain still. Though with a high-fidelity FE model, the dynamics of the small-scale substructures can be well represented, the interest is to characterize the global dynamics of the structural frame. Therefore, to approximate the global dynamics, it is unnecessary to compute all the modes, especially the local modes that have insignificant contributions to the global dynamics. In addition, the prediction of the local dynamics is different and usually more complicated than the prediction of global dynamics. One of the most popular methods to tackle this difficulty is to select the dominant modes by the modal participation factor (or mass participation factor) [22, 23]. This approach has been used in many applications to remove the insignificant modes from the structural dynamics analysis, for example, the application on rail vehicles [24], bridges [25], and PWR fuel assemblies [8]. The modal participation factor is dependent on mass normalization and independent of the external forces. In addition, it is often 8 the case that the cumulative sum of the mass participation factors of the low-frequency modes practically reaches the total sum of all the modes; that is, this method is more suitable for cases when high-frequency dynamics is not crucial. Recently, a general methodology that addresses this problem has been proposed by constructing a global reduced order model (ROM) based on a global reduced-order basis (ROB) [26]. The ROB’s main idea is to alter the kinetic energy (represented by kinematics) to filter the local modes while the global dynamics can still be well approximated. The methodology has been used in cars [27] and PWR fuel assemblies [28], and in addition, it has been extended to multi-level reduction [29, 30]. The global ROM is highly dependent on the reduced kinematics, and hence, it has to be sophisticated enough to filter the insignificant local modes on the one hand and well represent the global dynamics that might have coupled contribution from the local vibrations on the other hand. Modal analysis based on the Rayleigh-Ritz method [23, 22] is a widely-used model reduction method that permits an approximation of the dynamics by a reduced set of elastic modes. In which the reduced set considers the elastic modes that have frequencies less than some pre-specified value. The ROM obtained can predict the low frequency (LF) dynamics covered by the modal analysis. The modes in the LF range are usually of long-wavelength type and well separated, and thus, the modal analysis can be compact and accurate. The vibration modes are the eigenvectors of a generalized eigenvalue problem (GEP) associated with the free vibration of the structure. In structural dynamics, some of the most popular eigensolvers are shift-invert Lanczos (SIL) algorithm [31, 32], the subspace iteration method (SIM) [33, 34], Craig-Bampton (CB)sub-structuringtechnique[35,23], andautomatedmultilevelsub-structuring(AMLS)[36,37]. Allthese methods are based on Rayleigh-Ritz analysis with the solution of the GEP obtained by solving the projected system to a reduced subspace. The SIL and SIM are iterative methods that converge to the true solution, while the CB and AMLS give approximated eigensolutions. By taking advantage of the eigenvalue shifts that enable one to solve the GEP by many slices (spectrum slicing), the SIL reduces the number of eigenvectors required to maintain the orthogonality. A proper selection of starting vectors can improve the efficiency of SIM. The CB method and AMLS are based on domain decomposition. The structure is partitioned into many substructures, and the approximation is in the space spanned by a subset of component vibration modes of the substructures and the static constrained modes associated with the internal boundaries of the 9 substructures. The accuracy of these two methods depends on the modal truncation for the substructures. The AMLS is an extended version of the CB method, with automatic recursive decomposition and interface reduction. The sub-structuring techniques are well suited for the problem with high modal density since the expensive computation of maintaining orthogonality between numerous modes is avoided. The UQ analysis is crucial for problems where significant uncertainties are presented in the input space. Although formally, the convergence of standard MC simulation is independent of the number of random parameters, it is a function of the unknown variance and distribution of the solution being sought, which is likely to grow as the number of these parameters increases. The probability density evolution method [38, 39, 40] decouples physical systems by combining the principle of preservation of probability and random event descriptionandajudiciouschoiceofrepresentativepointsinparameterspace[41]. EfficientUQmethodshave also been introduced that leverage the topology of dynamical systems with localized uncertainties, yielding orders of magnitude speed-ups over standard MC sampling [42, 43, 44]. For more general structures, an appealing idea to attenuate the computational burden is to find a surrogate model that is easy to solve to replacethephysicalmodel. PCE[14]hasbeenappliedinmanydomainsofscienceandengineering[45,46,47] toconstructsurrogatemodelsbycharacterizingquantitiesofinterest(QoI)bytheircoordinatesrelativetothe spacespannedbyabasisofmultivariateHermitepolynomials. ExtensionsofthestandardPCEhaveincluded hybrid bases [48], multi-resolution bases using wavelets [49], and generalized PCE [50] using the whole range of orthogonal polynomials [51, 52]. In the PCE context, the coefficients in the polynomial representation are called the PCE coefficients and are obtained by either intrusive [53, 54, 52] or non-intrusive methods [14, 55, 56]. For the former, minimizing the error in the governing equation yields a system of coupled equations for the PCE coefficients. These generally require the development of newly adapted solvers. In the non-intrusive methods, the coefficients are evaluated as generalized Fourier coefficients in the form of integrals that require function evaluations of the original problem. The numerical quadrature is usually used to compute the PCE coefficients; however, in the case of full tensor product implementation, the number of quadrature points increases geometrically with the dimension [56]. The sparse quadrature construction [57,58], suchasSmolyakalgorithm[59], canbeusedtotemperthecurseofdimensionalitydrastically. Further attempts, including the least angle regression [60, 61], compressive sensing [62, 63], and other techniques 10 [64], have been made to represent the QoI by only a selection of significant polynomials. While these methods achieved measurable success, the challenge remains as the stochastic dimension increases or the space-time discretization is refined. A key driver for the challenge described above is that the quantities being characterized (i.e., the solution of a differential equation) are very complex objects, namely stochastic processes, whose probability measures, and therefore relevant computational characteristics, are functions of the random parameters. In many applications of interest, the quantity of interest (QoI) is a much simpler object, namely a random variable characterizing failure or cost. The idea of basis adaptation in homogeneous chaos space was proposed [6] recently to leverage the simplicity of the QoI and has successfully applied to practical engineering problems [65, 66]. Accordingly, the Gaussian input variables are rotated by an isometry to obtain a new orthogonal basis such that the PCE of the QoI based on the new basis has concentrated representation in the first several dimensions. The dimension of the UQ is thus reduced, and a small number of physical model evaluations are required. This method is extended from scalar QoI to random fields [67] and combined with compressive sensing to obtain efficient schemes [68]. Furthermore, a partial least squares (PLS) based adaptation is proposed to obtain the basis transformation by taking advantage of the covariance of the input parameters [69]. Though the basis adaptation can converge to an accurate representation of the QoI in the first several dimensions for many applications, it can have slow convergence for some particular problems. Hence, extensive dimensional adaptation is required to achieve high accuracy. The active subspace (AS) [70, 71] is another method to detect directions of large variability by the eigenvalue decomposition of the covariance matrix of the gradient, and the most important directions are used to construct a response surface of low dimensions. Although the AS method can provide appreciable dimension reduction, the calculation of the gradient of the model can be challenging, and it may introduce bias at the same time [72]. Both the basis adaptation and AS methods are adapted to scalar QoI, developing dimension reduction techniques for multi-dimensional QoI remains an open research domain. For the uncertainty induced from modeling errors, a nonparametric approach was introduced in [73] for linear structural dynamics. The method first constructs a ROM from the dynamic model. Then, a stochastic ROM (SROM) is constructed by randomizing the linear ROM’s underlying matrices. The distributions of the random matrices are obtained by the Maximum Entropy principle [74, 75, 76], where each of the matrices 11 has a dispersion hyper-parameter to control the level of uncertainty. The method has been applied to various applications [77, 78] and has extended to nonlinear models [79]. The inverse analysis involves an indirect identification of the parameters given some observed data. Dynamic system identification methods can be found in [80, 81, 82, 83]. In practical applications, the observation data is inevitably subjected to noise and uncertainty. Also, the computational model might have limitations to precisely predict the response given the parameters. Thus, quantifying the uncertainties in the input parameters is essential for the identification and decision-making process. Bayesian approaches [15, 16, 17] provide a foundation for inverse analysis where the observation data is polluted. The output of Bayesianinferencemethodsisusuallyposteriordistributionsoftheparametersthatarenotinanalyticalform and are propositional to the product of the likelihood of the data and the prior probability of the parameters. The posterior distribution contains complete probabilistic information of the parameters given the observed data. However, since no analytical forms of the posterior are given, the distributions are often approximated by samples generated by numerical approaches, for example, the famous MCMC algorithms [18, 19, 20]. The sampling methods require various evaluations of the physical model at different input parameters, which makestheBayesianmethodsprohibitedwhenthecomputationalmodelisextensivelyexpensive. Attemptsto attenuate the computational burden have been made in the literature, for example, in [84, 85, 86, 87, 88, 89], which are based on reductions or surrogates modeling of the computational model. Besides the difficulty embedded in the numerous physical model evaluations, another major challenge for MCMC methods is the high-dimensionality of the parameter space that arises from various sources of uncertainty or discretization of the infinite-dimensional parameter fields. The Markov chain must explore the high-dimensional parameter space to find regions that have a high posterior probability. The high-dimensionality also impacts the model reduction or surrogate modeling of the forward model. To increase the acceptance probability in MCMC and thusacceleratethesamplingforhigh-dimensionalparametersspaces, theadaptiveMCMCmethodswereused in [90, 91, 92], the sequential MCMC was used in [93], the likelihood-informed MCMC was applied in [94], and an accelerated MCMC by AS was proposed in [95]. The more recent research [96, 97] uses Generative Adversarial Network (GAN) to map from latent low-dimensional variables to the high-dimensional input variables and thus address the high-dimensionality issue. However, the research that can reduce the cost of 12 the expensive forward model and reduce the dimension of the parameter space simultaneously or reduce the forward model with high-dimensional parameter space is limited. MFUQ techniques are designed to obtain credible statistics of the QoI from HF and computationally demanding models and have recently received significant attention in the UQ field. The main idea is to overcome the growth of computational power needed for the UQ of HF computations by leveraging numerous evaluations of low-cost LF models. The combination of information from various levels of approximation, for example, different levels of spatial/temporal discretization and/or altered physical models, can be based on either sampling or surrogate approaches [98, 99, 100, 101, 102, 103]. The thesis will focus on sampling- based methods since they are independent of the dimension and are more robust. The prototype algorithm in the family of MFUQ sampling strategies is the so-called Multilevel Monte Carlo (MLMC) which was first proposed by Heinrich [104] and, successively, further developed and generalized by Giles in the seminal contribution [105]. A detailed presentation of MLMC with an exhaustive review is presented in [106]. Several extensions have been proposed, for instance, [107]. The idea of an ordered sequence of models has been leveraged and extended in the so-called multi-fidelity MC (MFMC), in which a control variate is recursively applied to the sequence of models. Combinations of control variate and MLMC led to other approaches like [108, 109, 110]; however, all these strategies are based on hierarchical relationships among models, i.e., an ordered sequence of models is constructed according to their correlation and cost properties. Gorodetsky et al. [111] demonstrated that MLMC and MFMC could be interpreted as a particular instance of general construction, called Approximate Control Variate (ACV), which is virtually capable of exploiting better non-hierarchical relationships among models in order to provide a larger variance reduction, especially in the case of limited resources available for the highest fidelity models. Geraci et al. provided an example of its application in the context of UQ for computer networks [112]. Though many advancements of the MFUQ technique can be found in the literature, some challenges remain in the field that decreases the applicability [21]. The thesis focuses on the issue exhibited in dissimilar or heterogeneous parameterization characteristics among the HF and LF models. As described in Section 1.2, to increase the efficiency of the MFUQ technique, the LF models usually need to alter the physics of the HF model, which in turn reduced the correlation among the models and, consequently, a need for a more significant number of HF runs. 13 By resorting to a dimension reduction step, [113] proposed to use AS to enhance the correlation among the models with dissimilar parameterization. However, same as in the UQ analysis, the AS required the computation of the model’s gradient, which requires many evaluations of the physical model to obtain an accurate estimation, especially in the case with high-dimensional parameter space. Also, if the HF model is not accurately captured by the active dimensions discovered by the method, a bias in the statistics will be introduced. 1.4 Objectives, methodologies, and contributions The standpoint of the research related to the challenges we are trying to address can glance from the literature review, based on which the objectives, technical methodologies, and main contributions of the thesis are presented in the current section. 1.4.1 Objectives As described before, the structures we are interested in for this study are large-scale models of multi- component systems with intertwined local elastic modes and global elastic modes in a broad frequency band. The particular application in the thesis is the integrity assessment of the structural internals of the FLSNFC, analyzed by a FE model of high DOFs. The first objective of this study is to provide an efficient vibration analysis method with high accuracy for structures with numerous local modes. The dynamic model serves as a forward physical model. Then, the second objective is to provide an efficient forward UQ analysis adapted to the high-dimensional QoI, the FRF on the exterior surface of the model, and the high-dimensional random input space. The third objective is to perform an efficient inverse analysis of the FLSNFC with high- dimensional input space to identify the probability distributions of the parameters given observed data from some damage scenarios. A successful accomplishment of these three objectives will enable us to infer the internal damages by dynamic signals from the exterior surface of the canister. The last objective of the thesis is to propose robust strategies to construct efficient MFUQ estimators that can deal with models with dissimilar parameterization. 14 1.4.2 Methodologies Using a global ROM to represent the global dynamics [26] while avoiding the calculation of local modes is first improved and extended in this study, and the intention is to develop a representative small dimension subspace to describe the global dynamics of the complex multi-scale structures. An automatic procedure for selecting the most significant elastic modes is introduced to build an efficient global ROM that is not affected by the local vibration modes. An alternative methodology that does not require the computation of the numerous elastic modes is introduced, in which a global ROB for the construction of global ROM is proposed. In order to take into account the influence of the local vibrations on the global vibrations, several reduced kinematics are introduced. The use of simple reduced kinematics enables one to have a decent approximation of the dominant global subspace of the complex fuel assembly. Details of the methodologies can be found in [114]. However, we will show later that the global ROM method cannot obtain a dynamic solution that is accurate enough for the inverse problem (since the inverse analysis is usually sensitive to errors). Therefore, a second attempt on dynamic reductions will be carried out. In the FE model of the FLSNFC, the connections between the fuel assembly and the structural frame are localized. The fuel assembly is connected to the bottom plate of the canister through a nose-piece structure with its extremity node, and it is connected to the honeycomb basket nearly the top by several springs. Therefore, the CB method is well-suited for the FLSNFC, in first that the structure is of high modal density, second, the connections between substructures are localized such that the internal boundaries DOFs are small, and lastly, the fuel assemblies are identical and the cost to calculate the substructure modes is reduced. However, the high modal density leads to an extensive CB model, and the direct eigensolver for the CB GEP is prohibitive. Therefore, the SIL and spectrum slicing are used to overcome the difficulty of solving the CB GEP, enabling parallel and efficient solving. In addition to applying the CB technique at the system level, the modes of the fuel assembly can also be calculated with the CB method. The benefits of this lower level CB are simultaneously reducing the computational cost and enabling an efficient re-analysis with localized modifications of the fuel assembly (this is especially useful for the UQ analysis). By doing so, a two-level nested CB method is proposed for the FLSNFC. The most dominant modes of the substructures are used in the CB model to decrease the cost further. In which an efficient method to determine the 15 importance of each of the sub-structural modes is introduced. The accuracy of the CB model depends on the sub-structural truncation and is fully controllable. Details of the methodologies for this part can be found in [115, 116]. For the UQ analysis adapted to scalar QoI, we propose two novel methods to accelerate the convergence of the original basis adaptation approach, thus expanding its applicability while also providing insight into its performance. Ingeneral, toconstructtherotationmatrixforbasisadaptationmethods, apilot(i.e., aninitial guess) first-order PCE in the original (non-rotated) space is first built. Our first method corrects the mean and Gaussian coefficients of the adapted PCE by using the available information embedded in the pilot PCE, increasing the convergence rate by eliminating the need to enforce first and higher-order accuracy separately. The basis adaptation intends to find a converged adapted PCE by sequentially increasing the adapted dimension. The second acceleration method proposed in this work takes advantage of the probabilistic information encoded in higher dimensional adaptations to correct the rotation matrix such that additional probabilistic concentration is achieved in the first few dimensions of the new rotation. The updating process is framed as an optimization problem. By continuously updating the rotation matrix using the highest dimensional adapted PCE available, the wealthiest probabilistic information on hand will be concentrated in the first few dimensions of the updated rotation, and we are more likely to find a low-dimensional subspace to capture the probabilistic behavior of the QoI accurately. The second method can be combined with the first method to achieve even better performance. Details of these two methods can be found in [117]. For stochastic reduction of multi-dimensional or even high-dimensional QoI, we propose a method that integrates the Karhunen-Loève expansion (KLE) [14] and the adapted PCE. In the first step, the multi- dimensional QoI is expanded by the KLE and further approximated by a truncated sequence of weighted KL basis. The KL basis is constructed by eigenvectors of the covariance matrix of the QoI, and the associated eigenvalue characterizes the contribution of each basis to the covariance. Since the eigenvalue commonly decays fast in most applications, a limited number of KL terms is sufficient to represent the covariance presented in the QoI. With the KL basis including a small number of terms built, in the second step, the adapted PCE model of the KL coefficients associated with these KL bases can be constructed via the accelerated basis adaptation method. The KL coefficient associated with each KL basis is a scalar quantity 16 that enables the application of the basis adaptation method for dimension reduction, and the construction of adapted PCE is only required for a handful of KL coefficients. Lastly, the probabilistic model of the KL coefficients can be substituted back to the KLE of the QoI to form a probabilistic surrogate model of the multi-dimensional QoI. This parametric approach can be combined with the nonparametric approach to account for uncertainties induced by modeling errors (or background noises). The integrity assessment of the internals is carried out through Bayesian inference, where the posterior distributions of the parameters are obtained with given data from a damaged model. The damage scenario underlying the data is intended to be identified from the posterior distributions of the parameters. To decrease the computational burden from the expensive forward physical model, the surrogate model of the QoI obtained from the integration of KLE and adapted PCE is used in the likelihood computation. Furthermore, to address the high-dimensionality of the parameters in posterior sampling, the block-update MCMC increases the acceptance probability. We propose incorporating the basis adaptation method into the MF techniques to deal with MFUQ problems with dissimilar parameterization among the models. The basis adaptation is used to explore important directions for HF and LF models independently, and these directions form a shared reduced space where the QoI has concentrated information on this space. Since the directions are aligned for different models, the correlation among the models when sampled on the shared reduced space is enhanced. Based on this fundamental idea, two strategies are proposed to construct efficient MF estimators, one is unbiased, and the other has balanced bias and variance. 1.4.3 Contributions The contributions of the thesis are listed as follows. 1. The method of constructing a global ROM by a global ROB [26] that alters the kinetic energy (repre- sented by kinematics) to represent global modes while filtering local modes is improved and extended; in particular, the scale separation of the model is exploited, and several reduced kinematics are intro- duced. 17 2. A multi-level nested CB method that combines domain decomposition, SIL eigenvalue solver, dominant mode selection, and CB sub-structuring technique is introduced for efficient dynamic analysis and is adapted to large-scale models with multiple structural levels. 3. Two novel methods are proposed to accelerate the dimension reduction technique in UQ, the basis adaptation method, in the framework of PCE. The methods can reduce the high-dimensional input space of UQ problems to only a few. 4. A method that integrates the KLE and the basis adaptation of PCE is proposed to address the di- mension reduction for high-dimensional UQ problems with also high-dimensional output space. The temporal/spatial dependent process is an example of the high-dimensional output. 5. A unified process is proposed to account for parameter uncertainties by parametric stochastic analysis, such as PCE, and uncertainties induced by modeling errors by nonparametric stochastic analysis. 6. An efficient Bayesian inference method is proposed to replace the physical model with the surrogate model obtained from the integration of KLE and the basis adaptation of PCE in generating posterior samples of the parameters. 7. The applicability and efficiency of the MFUQ techniques for problems with dissimilar parameterization areimprovedbyincorporatingthebasisadaptationmethod. Theimprovementisachievedbyleveraging the basis adaptation to enhance the correlation among models with different fidelity. 1.5 Outline The thesis is organized as follows. In Chapter 2, the modal analysis in structural dynamics is recalled, the automaticmodeselectionproceduretoselectthemostdominantelasticmodesisintroduced, theconstruction of the ROB with several reduced kinematics is presented, and the applications of these methods on the simplified and detailed fuel assemblies are described in the last part of this chapter. Chapter 3 described the two-level nested CB method adapted to the FLSNFC. In this chapter, the SIL, domain decomposition, and CB sub-structuring techniques are first recalled, followed by the derivations of how the CB method is adapted 18 to the FLSNFC. Then the further reduction by using the dominant modes of the substructures in the CB model is proposed. The last part of this chapter applies the adapted CB method to the FLSNFC. In Chapter 4, the PCE and basis adaptation method is introduced, followed by introducing two novel approaches to accelerate the basis adaptation. Then, the accelerated methods are applied to a borehole function and a space structure model. Chapter 5 introduces KLE, the integration of KLE and the adapted PCE, and the nonparametric stochastic analysis. Then, the methods are applied to the detailed FA and the FLSNFC to build surrogate models of the FRFs. The efficient inverse analysis is presented in Chapter 6, in which the Bayesian methods are recalled, the efficient algorithm is described, and the block-update MCMC is introduced for posterior sampling. Finally, the methods are applied to the detailed FA and the FLSNFC for damage identifications. Chapter 7 presents the incorporation of basis adaptation into MFUQ techniques to deal with MFUQ problems with dissimilar parameterization. First, the MFUQ techniques are recalled, the idea of the incorporation is introduced, and two strategies are proposed to develop an efficient unbiased estimator and an estimator with bias and variance balanced. Then, these strategies are applied to an acoustic problem. Finally, the conclusion and future perspectives are given in Chapter 8. 19 Chapter 2 Automatic mode selection and global reduced-order model In this chapter, the modal analysis of the linear structural dynamics is recalled in section 2.1. Then we present the automatic mode selection procedure to find the most dominant elastic modes in section 2.2. This is followed by section 2.3 that introduces the methodology of constructing a global ROM that can filter the local vibrations associated with the flexible part of the structure while capturing the contributions of these local vibrations to the global dynamics at the same time. Finally, section 2.4 presents the applications of these two methods to the simplified and detailed fuel assembly. 2.1 Modal analysis in structural dynamics This work focuses onthedynamic characterization in the frequency domain. For all! in a specifiedfrequency bandB = 2]0;f u ] with f u the upper bound in Hz, the unknown displacement field U(!) is given by the following matrix equation ! 2 [M] +i! [D] + [K] U(!) =F(!); (2.1) where [M], [D], and [K] are the FE mass, damping, and stiffness matrices, which, without loss of generality, are considered to be symmetric and positive definite; F(!) is the discretized external force at frequency !; i denotes the imaginary unit. The dimension of the nodal DOF is denoted as N, which can be very large for complex structures. However, the FE matrices are also very sparse. Solving Eq (2.1) directly for many frequency points and external forces can be computationally prohibitive. 20 A traditional way to reduce the computational cost for solving Eq (2.1) is to construct a ROM using vibration eigenmodes. The eigenmodes are obtained by solving the following GEP [K]' =[M]'; (2.2) where ' is an eigenmode or eigenvector and is the associated eigenvalue. The eigenvalue in Hz can be obtained by p =2. There are N eigenpairs, but not all are necessary to capture the accurate response within the frequency bandB. In practice, the modes that have eigenfrequencies below some f c > f u are considered. The cutoff frequency f c is usually determined by convergence analysis such that the significant spread contributions of the modes that have eigenfrequencies larger than f u are considered. The number of modes in the frequency bandB c = 2]0;f c ] is denoted as n e with n e N. Eq (2.2) can be written in matrix form as [K][] = [M][][]; (2.3) where [] = [' 1 ;:::;' ne ] is the ROB, which contains the modes that have eigenfrequencies in the frequency band B c = 2]0;f c ]; and [] is a diagonal matrix with the diagonal entries the first n e eigenvalues. The eigenmodes are orthogonal with each other and can be normalized arbitrarily. In this work, mass normalization is used and satisfies [] T [M][] = [I ne ]; [] T [K][] = []; (2.4) where [I ne ] is the identity matrix of dimensionn e . Though the projection of the damping matrix, [] T [D][], is not diagonal in general, the modal damping model is used in this paper such that [] T [D][] = 2[][ ], where [] = diag( 1 ;:::; ne ) is the diagonal matrix of the modal damping ratios and [ ] = p []. The classical ROM is obtained by the following approximation [U](!) ne X =1 q (!)' = []q(!); (2.5) 21 where q(!) = (q 1 (!);:::;q ne (!)) is the generalized coordinates and is obtained by solving the following equation ! 2 [I ne ] + 2i![][ ] + [] q(!) =F(!); (2.6) whereF(!) = [] T F(!). Eq (2.6) can be solved very efficiently for numerous frequencies and external forces. As a result, the computational cost of the classical ROM mainly depends on the cost to build the ROB of Eq (2.2), and the dimension of the ROB, n e . Several structural scales are presented for interested complex structures in this study, and numerous local vibration modes are mixed with the usual global vibration modes. These local vibration modes can significantly increase the modal density in the frequency band of analysis but have an insignificant influence on the global displacement of the structure. Hence, the dimension of ROB n e will be large, and the cost to build the ROB is high. This chapter addresses the problems of (i) finding a small dimensional ROM with acceptable accuracy and (ii) avoiding the calculation of numerous local vibration modes. 2.2 Automatic selection of the dominant vibration modes A way to rank the importance of the modes is presented. The importance is regarding the modal contribution to the global dynamic response of the structure so that one can remove the insignificant modes and still get a good representation of the global behavior by superposition of the remaining modes. Denoting the dynamic stiffness matrix as [A(!)] = ! 2 [M] +i![D] + [K] , then the generalized coordinatesq(!) in Eq (2.6) can be written as q(!) = [] T [A(!)][] 1 [] T F(!): (2.7) Define the matrix [H(!)] := [] T [A(!)][] 1 , which has components h (!) for = 1;:::;n e that are given by h (!) = ! 2 + 2i! ! +! 2 1 : (2.8) Let U ij (!) denotes the response at DOF i results from the unit force at DOF j, where i; j = 1;:::;N, then the linear system Eq (2.1) can be fully described by the FRF matrix [U(!)] whose components are 22 fU ij (!); i;j = 1;:::;Ng. The displacement induced by a unit load is called compliance. Eq (2.5) and (2.7) implies the approximation [U(!)] [][H(!)][] T , and the approximation is exact when n e is sufficiently converged in the frequency bandB. Hence, U ij (!) = ne X =1 U ij (!); (2.9) where U ij (!) is the contribution of the mode to the compliance and is obtained by U ij (!) = i h (!) j ; (2.10) in which k = [] k for k = 1;:::;N and = 1;:::;n e . The FRFs are usually in decibel scale which are given by u ij (!) = 20 log 10 jU ij (!)j: (2.11) Let ~ U ij (!) be an approximation of U ij (!), then the FRF error measure associated with these two FRFs is defined as ij =(u ij ; ~ u ij ) = s 1 B Z B (u ij (!) ~ u ij (!)) 2 d!: (2.12) The error measure is independent of the magnitude of the response; that is, it is normalized. Also, the error measure is the same for displacement, velocity, and acceleration responses. The importance of one mode to the FRFU ij (!) can be quantified by the FRF error of Eq (2.12) induced by its sole removal. Note that the interest is to quantify the importance of the mode to the global dynamics of the whole structure rather than one specific FRF. Thus, many combinations of observationi’s and excitation j’s have to be considered; we refer to all the i’s and j’s considered the “excitation and observation pool”. Since the purpose is to find a ROM that can accurately represent the global dynamics, the excitation and observation pool contains DOFs that are important for global responses but not sensitive to local dynamics of the flexible sub-components. For example, the representative locations on the stiff skeleton of the structure (or structural frame). Assuming that there are N O interested observations and N E interested excitations, the total output and input combinations are N O N E . Fori = 1;:::;N O andj = 1;:::;N E , the responses in 23 decibel scale,u ij (!) and ~ u ij (!) =u ij (!)u ij (!) (for = 1;:::;n e ), are first computed, then the associated (u ij ; ~ u ij ) is computed by Eq (2.12). For each mode, there are N O N E error measures(u ij ; ~ u ij ) obtained. There are two proposed measures to quantify the importance of one mode regarding the structure’s global dynamics. The first one is to use the maximum FRF error in the excitation and observation pool, = max( ij ; i = 1;:::;N O ; j = 1;:::;N E ): (2.13) Let = 1= (N O N E ) P N O i=1 P N E j=1 (u ij ; ~ u ij ) and 2 = 1= (N O N E 1) P N O i=1 P N E j=1 ((u ij ; ~ u ij ) ) 2 be the sample mean and variance, respectively, the second measure is defined as = +k ; (2.14) where k is a positive constant that controls the acceptable percentage of FRFs for which the FRF error is greater than . Note that in the criterion of Eq (2.13), 0% of the FRFs have an FRF error greater than . Their importance regarding the global dynamics sorts the modes. The ROB contains n d modes with the most significant importance , where n d is the minimum dimension such that u ij ;u (n d ) ij < tol ; 8i = 1;:::;N O ; and8j = 1;:::;N E ; (2.15) whereu (n d ) ij is the FRF in decibel scale of output i and inputj computed by the ROM that contains the n d most dominant modes. Note that although Eq (2.15) use the criterion of Eq (2.13) that requires the FRF error of all the cases in the excitation and observation pool to be less than a tolerance, one can also choose the criterion of Eq (2.14) such that a large percentage, for instance, 95%, of the combinations have FRF errors satisfying the tolerance requirement. This automatic mode selection procedure addresses the problem of finding a small dimension ROM. However, to obtain the importance of the mode, the n e elastic modes need to be computed, which is of high cost. Note that the more generic mode importance (or participation) indicator, the modal participation factor [23, 22], uses the effective mass of each mode to rank the importance. The modes are assumed to competently 24 represent the underlying dynamics when the cumulative mass ratio reaches one. However, this metric is dependent on mass normalization and independent of the external forces. In addition, it is often the case that the cumulative sum of the mass ratios of the low-frequency modes reaches the total sum of all the modes, which is not suitable for our problem where the high-frequency modes also have significant contributions. For the same reason, using a simple truncation at low-frequency is also not applicable to our problem. A similar (yet different) indicator proposed in [118] uses the trace of matrix integral of the transfer function (which is the so-calledH 2 norm) to rank the modal contribution to the FRF. However, it considers all the output FRFs, and the size of the transfer function in our case will be too large. Nevertheless, the indicators presented in Eqs (2.13) and (2.14) are adapted directly to the interested FRFs. 2.3 Construction of global reduced-order basis This section intends to construct a small dimension ROB without calculating the insignificant local vibration modes. The goal is to find a set of basis vectors whose shapes filter the local deformations but have a good representation of the global deformations. Such basis vectors are referred to as “global ROB”. To achieve this, the eigenvalue problem Eq (2.2) is altered such that the numerous local modes will not be present. Before introducing global ROB, the static [23] and Guyan condensations [23, 119] are recalled. Static and Guyan condensations. Suppose the DOFs of the structure are split into two sets, the master and slave DOFs. Applying the partition to Eq (2.2), we can get " K mm K ms K sm K ss # ' m ' s ! = " M mm M ms M sm M ss # ' m ' s ! (2.16) where the subscript m and s denote the master and slave DOFs. The second row of the equation yields ' s = ([K ss ][M ss ]) 1 ([K sm ][M sm ])' m = [t sm ()]' m ; (2.17) 25 where [t sm ()] := ([K ss ][M ss ]) 1 ([K sm ][M sm ]). Then,' can be written as' = [T ()]' m where the transformation matrix is given by [T ()] = " I Nm t sm () # ; (2.18) where N m is the number of master DOFs. The eigenvalue problem Eq (2.16) can be reduced to the space spanned by the columns of [T ()] by Galerkin projection, which yields [K R ()]' =[M R ()]'; (2.19) where [K R ()] = [T ()] T [K][T ()], [M R ()] = [T ()] T [M][T ()], andthesubscriptRdenotesthatthematrix is reduced. Different choices of [T ()] yield different condensations. The Guyan condensation assume that the dynamic effects on the slaves DOFs are zero, which corresponds to [t sm ()] = [t sm (0)], then, the transformation matrix becomes [T G ] = " I Nm t sm (0) # = " I Nm [K ss ] 1 [K sm ] # : (2.20) The reduced stiffness and mass matrix are [K G ] = [K R (0)] = [K mm ] [K ms ][K ss ] 1 [K sm ] and [M G ] = [M R (0)] = [T G ] T [M][T G ] = M mm M ms K 1 ss K sm K ms K 1 ss M sm +K ms K 1 ss M ss K 1 ss K sm . The eigenvalue problem in the reduced space is [K G ]' m =[M G ]' m : (2.21) The eigensolution of Eq (2.16) is approximated by' = [T G ]' m . Note that' s =[K ss ] 1 [K sm ]' m , that is the maser DOFs and slave DOFs are in static equilibrium. The static condensation assumes that the mass associated with the slave DOFs is zero; that is, [M ss ], [M ms ], and [M sm ] are all replaced by zeros in Eq (2.16). Since the mass elimination altered the mass matrix in Eq (2.16), the associated eigenmodes are denoted as instead of'. Plug the altered mass matrix into 26 Eq (2.17), and we obtain the same transformation matrix as Guyan condensation. The first row of Eq (2.16) yields the reduced eigenvalue problem for static condensation, [K G ]' m =[M mm ]' m ; (2.22) different from Eq (2.21). Let [ S ] denote the (Nn S ) real matrix whose columns are the first n S < m approximated eigenmodes obtained by static condensation of Eq (2.16). Projecting Eq (2.2) to the space spanned by the columns of [ S ] yields the following eigenvalue problem [ S ] T [K][ S ] & = [ S ] T [M][ S ] &: (2.23) Then the eigenmodes ' approximated by static condensation can be obtained by ' = [ S ]&. Note that Eq (2.23) is different from Eq (2.22) and is referred to as a “Rayleigh-Ritz” step. The condensation can reduce the memory and computation cost for solving Eq (2.2) when m N. Nevertheless, the goal is to filter the insignificant local modes, where, in such a case, the master DOFs are chosen such that the local modes will not be present. This can be achieved by separating the DOFs into stiff parts associated with the structural frame and flexible parts. When the stiff part has a large portion of the total DOFs, mN is invalid. General methodology to construct the global modes. Let [B] be a ROB of dimension (Nr), and S r be the subspace ofR N spanned by the columns of [B]. Assume that by performing Galerkin projection of Eq (2.2) to [B], the reduced representation has exact elastic energy and reduced kinetic energy. The general framework for finding the global modes for any arbitrary [B] is presented below. Suppose that, in this chapter, the inner product of any two vectors u; v2R N is given by < u;v> W = v T [W]u; where [W] is (NN) a real weight matrix, and is symmetric positive-definite. For any u2R N , defining u r :=Pu, which is the orthogonal projection of u on the subspaceS r , and P = [B] [B] T [W][B] 1 [B] T [W]: (2.24) 27 is the (NN) orthogonal projection matrix [120]. Let _ u be the time-dependent velocity associated with the displacement field u, and then the kinetic energy is E k [ _ u] = 1 2 _ u T [M] _ u. The kinetic energy associated with the orthogonal projection _ u r is E r k [ _ u] = E k [ _ u r ] = 1 2 ( _ u r ) T [M] _ u r = 1 2 _ u T [M r ] _ u, in which the (NN) low-rank mass matrix [M r ] is given by [M r ] = [P] T [M][P]: (2.25) By using the exact stiffness matrix [K] and projected mass matrix [M r ], the following GEP is introduced [K] = [M r ] ; (2.26) which is different from Eq (2.2). Since the rank of [M r ] is r, there are r eigenpairs in this GEP. Similar to static condensation, an additional Rayleigh-Ritz step is required to obtain the global modes. Let [()] = [ 1 ;:::; ] with <r, the following eigenvalue problem is solved [()] T [K][()] r = g [()] T [M][()] r ; (2.27) for = 1;:::;. The global mode and associated eigenvalue are ' = [()]r and g , respectively. Let [ g ] = [' 1 ;:::;' ng ] be the global ROB with n g < and [ g ] = diag( g 1 ;:::; g ng ) be the diagonal matrix of the first n g global eigenvalues, the mass normalization can be verified, that is [ g ] T [M][ g ] = [I ng ] and [ g ] T [K][ g ] = [ g ]. Note that in Eq (2.26) is different from g in Eq (2.26), and the difference depends on how well the reduced kinematics can represent the global dynamics near the frequency f g . In cases where the reduced kinematics are satisfactory, we can choose such that > 2f c . The dimension of the global ROM can be chosen in the same manner, such as g ng > 2f u . In the following section, several reduced kinematics based on the scale separation are proposed. Each of the kinematics is defined by the ROB [B] and weight matrix [W]. Note that in Eq (2.26), [M r ] is multiplied by a vector on the right, and hence the assembly of [M r ] is not necessary. The same logic also applies to the orthogonal matrix [P]. 28 2.3.1 Reduced kinematics bases on scale separation Assume a clear separation between the stiff and flexible parts responsible for global and local dynamics, respectively. The sub-domains of stiff part and flexible part are denoted as 1 and 2 such that = 1 [ 2 and 1 \ 2 =;. To filter the vibrations inside 2 , the following ROB [B] is used [B] = " I N1 0 # ; (2.28) where N 1 is the number of DOFs in 1 ; [0] is of dimension (N 2 N 1 ), where N 2 is the number of DOFs in 2 . Assume that the matrices of the structure are partitioned accordingly. In this case, the dimension of the reduced subspaceS r is r =N 1 . Figure 2.1 depicts a multiscale vibration mode of a heterogeneous plate with two scales: a stiff skeleton and 12 flexible panels to present the filtering effects of different reduced kinematics. In which there is global deformation at the stiff part but also local deformations at the flexible panels. Figure 2.1: Multiscale mode of a heteroge- neous plate Figure 2.2: Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.29) Exact in stiff part and zero in the flexible part. Consider using weight matrix [W] = [M], and in this case, the kinetic energy difference (E k ( _ u)E r k ( _ u)) after the reduction is minimized [114] for any _ u2R N . Plugging Eq (2.28) into Eq (2.24), the orthogonal projection matrix is obtained as [P] = " I N1 M 1 11 M 12 0 0 # ; (2.29) 29 from which the projected mass matrixM r can be obtained as [M r ] = " M 11 M 12 M 21 M 21 M 1 11 M 12 # : (2.30) As mentioned before, it is not necessary to build and store [M r ] explicitly but rather compute the product of [M r ] with a vector on the right. For a given _ u2R N , the productp = [M r ] _ u can be given by p = M 11 _ u 1 +M 12 _ u 2 M 21 _ u 1 +M 21 M 1 11 M 12 _ u 2 ! ; (2.31) in which one needs to solve a linear system of dimension N 1 . Figure 2.2 presents the projection of the multiscale mode as seen in Figure 2.1 by orthogonal projection matrix [P] defined in Eq (2.29), from which we can see that the flexible part has zero deformation while the stiff part is kept exact everywhere except at the interface with the flexible part. Mass elimination or static condensation. In this case, [B] is also given by Eq (2.28), for the weight matrix, instead of using the mass matrix of the structure, the following is utilized [W] = " M 11 0 0 0 # : (2.32) It can be derived that [P] = " I N1 0 0 0 # ; [M r ] = " M 11 0 0 0 # ; p = [M r ] _ u = M 11 _ u 1 0 ! : (2.33) Which is the same as in static condensation, where the master DOFs are the DOFs belonging to 1 . Figure 2.3 presents the projection of the multiscale mode as seen in Figure 2.1 by orthogonal projection matrix [P] defined in Eq (2.33), from which we can see that the flexible part has zero deformation while the stiff part is kept exact. The difference between Figures 2.3 and 2.2 is due to the different choices of weight matrices. Note that for both cases in Figures 2.3 and 2.2, though there is a discontinuity between stiff and 30 flexible parts in the projected space, the eigensolutions of Eq (2.26) satisfy the static equilibrium at the interface, and thus, the interface is continuous. Figure 2.3: Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.33) Figure 2.4: Orthogonal projection of the multiscale mode in Figure 2.1 by projection matrix [P] defined in Eq (2.34) Exact in the stiff part and flexible part in static equilibrium, or Guyan condensation. Consider using the Guyan transformation matrix [T G ] (Eq. 2.20) as the ROB [B], and the weight matrix defined in Eq (2.32), and the following can be derived [P] = " I N1 0 K 1 22 K 21 0 # ; [M r ] = " M G 0 0 0 # ; p = [M r ] _ u = M G _ u 1 0 ! ; (2.34) where the subvectorp 1 can be written in detail by substituting the expression of [M G ] as p 1 = M 11 M 12 K 1 22 K 21 K 12 K 1 22 M 21 +K 12 K 1 22 M 22 K 1 22 K 21 _ u: Which involves solving two linear systems of dimension N 2 . Solving the GEP in Eq (2.26) by substituting [M r ] defined in Eq (2.34) is equivalent to solving GEP in Eq (2.21) and obtaining the displacement associated with DOFs in 2 by static equilibrium. Thus, this approximation is equivalent to Guyan condensation. Note that in this case, the Rayleigh-Ritz step in Eq (2.27) is not necessary since the resulting is already in the subspaceS r . Figure 2.4 presents the projection of the multiscale mode in Figure 2.1 by orthogonal projection matrix [P] defined in Eq (2.34), from which we can see that the stiff part is kept exact. In contrast, the flexible part deforms accordingly such that the static equilibrium between the stiff part and the flexible part is satisfied. 31 Exact in the stiff part and flexible part in static equilibrium, with minimum residual kinetic energy. This approximation is the same as in the previous case except that the weight matrix uses [W] = [M]. Please refer to [114] for detailed formulas. Figure 2.5 presents the projected multiscale mode in Figure 2.1, from which we can see the local deformation is reproduced, and the deformation of the stiff part is changed to satisfy the static equilibrium. This implies that the choice of [W] = [M] may not be suitable to filter the local deformations. The difference between Figures 2.5 and 2.4 is due to the different choices of weight matrices. Figure 2.5: Projected multiscale mode in Figure 2.1 by the reduced kinematics of exact in the stiff part and flexible part in static equilibrium (with minimum residual kinetic energy) The static condensation and Guyan condensation described in this section are applied to the hetero- geneous plate. We changed the mass of the flexible panels for several cases and found out that: (i) the static condensation and Guyan condensation give similar global ROM; (ii) the static condensation can well approximate the global dynamics of the plate with a small dimensional global ROM for the case when the mass of the stiff part is dominant; (iii) with the increase of mass of the flexible panels, both the complexity of the dynamics and the contribution of flexible modes increase, the dimension of the global ROM increases too; (iv) in the case where the mass of the flexible panels is dominant, it is hard to find a small dimensional ROM by static condensation that can represent the global dynamics well. 32 2.3.2 Reduced kinematics without scale separation This section considers structures with numerous local vibrations but no substructure supporting these local behaviors, for example, very complex structures like cars. It is a continuation of the work done by Soize and Batou [26, 28, 121]. The structure is first partitioned into N S non-overlapping subdomains i such that = Ns [ i=1 i ; 8i6=j i \ j =;: (2.35) The second step is to introduce the orthogonal projection operator [26],h r , by fh r (u)g (x) = Ns X i=1 1 i (x) 1 m i Z i (x 0 )u(x 0 )dx 0 ; (2.36) where u(x) is the displacement field; (x) is the mass density at the pointx 0 ; m i = R i (x)dx is the mass of the subdomain i , and 1 i is the indicator function. The discretization [H r ] ofh r is a special case of [P]. To obtain [H r ] by Eq (2.24), letting the weight matrix [W] = [M], and [B] be the following [B] = 2 6 6 6 6 6 4 S 1 0 0 0 S 2 . . . . . . . . . . . . . . . 0 0 0 S Ns 3 7 7 7 7 7 5 with [S i ] = 2 6 6 6 6 4 T i (X i;1 ) T i (X i;2 ) . . . T i (X i;N i ) 3 7 7 7 7 5 ; (2.37) where [T i ()] (6r i ) is the nodal transformation matrix, with r i the dimension of kinematics in subdomain i ;X i;j is the coordinate of node j in subdomain i ;N i is the number of nodes in subdomain i . If the displacement is constant in each subdomain, [T i ()] is independent of the space and the subdomain. We call this case piece-wise constant kinematics, in which the nodal transformation matrix has the following, [T i (X i;j )] = [T ] for all i = 1;:::;N S and j = 1;:::;N i , where [T ] = 2 6 6 6 6 6 6 6 6 4 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 3 7 7 7 7 7 7 7 7 5 : (2.38) 33 Note that it is assumed that the DOFs within every node are sorted byTX; TY; TZ; RX; RY; RZ, which are translation and rotation DOFs along X; Y; Z. It can be deduced that the dimension ofS r for piece- wise constant kinematics is r = 3N S . The local vibrations within each subdomain are not represented; the domain decomposition determines these filtering effects. In some cases, it is useful to use enriched kinematics within the subdomains. For the piece-wise rigid-body kinematics where each subdomain can have rigid-body displacement, the nodal transformation matrix is given by [T i (X )] = 2 6 6 6 6 6 6 6 6 4 1 0 0 0 Z Y 0 1 0 Z 0 X 0 0 1 Y X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 7 7 7 7 7 7 7 7 5 ; (2.39) where (X ;Y ;Z ) is the coordinate or node X in subdomain i . The dimension of S r for piece-wise constant kinematics is r = 6N S . The nodal transformation matrix for piece-wise linear kinematics is [T i (X )] = h 1[T ] X [T ] Y [T ] Z [T ] i : (2.40) Similarly, a piece-wise polynomial approximation has the following nodal transformation matrix, [T i (X )] = h p 1 (X )[T ] p 2 (X )[T ] ::: p ni (X )[T ] i ; (2.41) where n i is the number of polynomials considered for the reduced kinematics of subdomain i and p k with k = 1;:::;n i is a polynomial function of space. Figure 2.6 presents the shape of the orthogonal projection of the multiscale mode of Figure 2.1 by the piece-wise constant approximation, in which the local displacements are filtered and that the global deformation shape is reproduced very roughly. Figure 2.7 presents the shape of the orthogonal projection of the multiscale mode of Figure 2.1 given by the piece-wise rigid- body approximation, in which the local displacements are filtered and that the global deformation shape is approximately reproduced. Figure 2.8 presents the shape of the orthogonal projection of the multiscale 34 mode of Figure 2.1 given by a degree-5 polynomial approximation over the whole structure domain. The local displacements are filtered and the global deformation shape seems to be reproduced fairly accurately. Figure 2.6: Orthogonal projection of the multiscale mode in Figure 2.1 given by the piece-wise constant approximation Figure 2.7: Orthogonal projection of the multiscale mode in Figure 2.1 given by the piece-wise rigid-body approximation Figure 2.8: Orthogonal projection of the multiscale mode in Figure 2.1 given by a degree-5 polynomial approximation over the whole structure domain 2.3.3 Combined kinematics Thestructuresinterestedinthisstudyhaveclearscaleseparationsandthestrategiesinsection2.3.1and2.3.2 canbothbeusefulinthiscase. Morespecifically, thestrategyinsection2.3.1cankeepanexactrepresentation of global deformation while zeroing out the local displacements. However, it may eliminate the flexible parts’ inertia, which can significantly affect the global dynamics. For example, the static condensation in Eq (2.33) eliminates the inertia of the flexible parts, and the coupled effects between stiff parts and flexible parts are not considered in the resulting global modes. In addition, the eliminated mass can cause the shift of global eigenfrequencies, in which the shifting is significant when the flexible part has a large portion of the total mass. The strategy in section 2.3.2 takes the inertia of the flexible part into account, but the approximation 35 of the global deformations is dependent on the local fluctuations. The averaged displacement within each subdomain can be dominated by the local vibrations when their amplitude is large and cannot cancel out. The strategies, particularly the static condensation in section 2.3.1, and the reduced kinematics in section 2.3.2 are combined, in which the reduced kinematics are only applied to the flexible substructures. With proper choice of reduced kinematics, the combined strategy can (i) take into account the global dynamics due to the coupling of the flexible substructures and the stiff substructures, (ii) avoid representing the numerous local modes due to the local vibrations of the flexible substructures, and (iii) correct the shifting of the global eigenfrequencies. In this case, the reduced basis [B] is constructed as [B] = " I N1 0 0 b # ; (2.42) in which [b] is an (N 2 r 2 ) reduced basis for subdomain 2 (the flexible substructures). The representation of subdomain 1 (the stiff skeleton) is kept exact. The strategies presented in section 2.3.2 are used to construct [b] such that the global dynamics of the flexible substructures (which emerge by coupling with the stiff skeleton) is represented. By using the weight matrix [W] = [M] and further simplifying it by letting [W 12 ] = [0] and [W 21 ] = [0] (the error of the simplification is small since the coupling terms of the mass matrix are insignificant), the projection matrix [P] and reduced mass matrix [M r ] can be obtained. Again, in practice, [P] and [M r ] are not assembled, and the productp = [M r ] _ u is instead computed in preparation for solving the GEP Eq (2.26). The lumped mass approximation is also used, introducing negligible errors when the FE model uses a fine mesh. The GEP Eq (2.26) is solved very efficiently using the multiple shift-inverted Lanczos (SIL) algorithm [122, 123], which will be detailed in later chapters. 2.4 Applications This section presents the application of the methods described in sections 2.3.1, 2.3.2, and 2.3.3 to two models, namely, the simplified and detailed models of a fuel assembly. For both applications, the automatic mode selection is used and compared with other ROMs. 36 2.4.1 Simplified nuclear fuel assembly The simplified nuclear fuel assembly is a slender structure composed of a 10 10 regular array of fuel rods and ten uniformly placed tie plates that hold the rods together. The tie plates are modeled by shell elements and the fuel rods by beam elements. The FE model has N = 368; 892 DOFs. The total mass of the fuel rods is 262.0 kg, whereas the total mass of the structure is 276.4 kg. There are 5422 elastic in the interested frequency bandB = 2]0; 1000] rad/s, and 5441 elastic modes below 1200 Hz. More specifically, there is a limited number of global modes and numerous local modes due to the bending of the flexible fuel rods. Figure 2.9 presents the fuel assembly, a global mode, and a local mode. The stiff part consists of ten tie Figure 2.9: Simplified fuel assembly: undeformed configuration (left), a global mode (center), a local mode as seen from two points of view (right) plates for this structure, while the fuel rods support the flexible part. However, one should note that the 100 fuel rods constitute the skeleton and this beam-like structure still exhibits global modes. The reference ROM uses the firstn e = 5441 elastic modes as its basis. In addition, four additional ROMs are obtained. We are interested in the global dynamics of this structure. For comparison purposes, two control FRFs are considered; for both cases, the external force is applied at the edge of one intermediate tie plate inx direction (z is the longitudinal direction), for all! inB. The force is applied off-center of the edge such that the torsion modes will be present. Control FRF 1 is the acceleration in x direction of one node in the opposite direction of another intermediate tie plate, also off-center. Control FRF 2 is the acceleration in the y direction of the same observation location. The automatic mode selection procedure is applied to find dominant elastic modes. A subset of 240 DOFs based on 80 randomly selected nodes from the ten tie plates are served as the excitation and observation pool 37 to determine the importance of the elastic modes. Three global ROMs are constructed based on a domain decomposition that the fuel assembly is cut into 100 transverse slices uniformly. The first is constructed by using the piece-wise constant kinematics for all the subdomains. The second one is constructed by using piece-wise rigid-body kinematics. The third global ROM is based on the combined kinematics, in which the domain decomposition is only applied to the fuel rods while kinematics of the tie plates is kept exact. Finally, Piece-wise rigid-body kinematics is applied to the fuel rods. Note that, since the fuel rods represent about 94.8% of the total structure mass, the global ROM by static condensation does not give satisfactory results, which is not shown for brevity. The global ROM by piece-wise constant kinematics has n g = 79 global modes inB. The prediction of control FRF 1 is shown in Figure 2.10, in which the overall prediction is satisfactory though some peaks are missed in the low- and middle-frequency range, and a noticeable discrepancy is presented in the high- frequency domain. The prediction of control FRF 2 is presented in Figure 2.11, which shows that this global ROM fails to predict control FRF 2. This can be explained that the piece-wise constant kinematics is unable to represent torsion modes. 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact piecewise constant Figure 2.10: Control FRF 1 calculated with the global ROM of dimensionn g = 76=5441 obtained by piece-wise constant kinematics (constant by slice) and for which = 5:6 dB 0 200 400 600 800 1000 Frequency (Hz) -200 -150 -100 -50 0 Accelerance (dB) exact piecewise constant Figure 2.11: Control FRF 2 calculated with the global ROM of dimensionn g = 76=5441 obtained by piece-wise constant kinematics (constant by slice) and for which = 51:0 dB The global ROM by piece-wise rigid-body kinematic has n g = 113 global modes inB. The predictions of control FRF 1 and FRF 2 are shown in Figures 2.12 and 2.13. We can see that they can accurately predict the LF responses, and the middle frequency (MF) responses are quite satisfactory for both cases, though the high frequency (HF) predictions are not as good. The predictions of the global ROM by combined kinematics give similar results as that of using piece-wise rigid-body kinematics, as can be seen in Figures 2.14 and 2.15, and note that the number of global modes 38 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact piecewise rigid body Figure 2.12: Control FRF 1 calculated with the global ROM of dimension n g = 113=5441 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 4:1 dB 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact piecewise rigid body Figure 2.13: Control FRF 2 calculated with the global ROM of dimension n g = 113=5441 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 3:7 dB inB is also n g = 113. Therefore, the remaining discrepancy (particularly in the frequency band ]800, 1000] Hz) is due to the reduced kinematics of the fuel rods. 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact combined kinematics Figure 2.14: Control FRF 1 calculated with the global ROM of dimension n g = 113=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact else- where) and for which = 4:2 dB 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact combined kinematics Figure 2.15: Control FRF 2 calculated with the global ROM of dimension n g = 113=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact else- where) and for which = 3:7 dB The predictions given by the ROM of the 113 auto-selected dominant elastic modes are shown in Figures 2.16 and 2.17, which are nearly perfect everywhere except for some low amplitudes. Regarding the numerical cost, thecomputationoftheglobalROMtookabout3minutes, whereasthecomputationoftheelasticmodes took about 1 hour and 40 minutes. From this application, we can see that even though the mass of the flexible part is dominant, we can still represent the global dynamics excellently by a small dimensional ROM and reduce the computational cost a lot. 39 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact dominant modes Figure 2.16: Control FRF 1 calculated with the ROM of the n d = 113=5441 dominant elastic modes obtained by the automatic mode selection procedure and for which = 1:4 dB 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact dominant modes Figure 2.17: Control FRF 2 calculated with the ROM of the n d = 113=5441 dominant elastic modes obtained by the automatic mode selection procedure and for which = 1:3 dB 2.4.2 Boiling water reactor fuel assembly The vibration analysis of the detailed fuel assembly that is used in boiling water nuclear reactors (BWRs) is considered in this application. The fuel assembly follows the GE14 design. The fuel assembly is a slender structure composed of a 1010 regular array of fuel rods held together by eight uniformly distributed spacer grids and the upper and lower tie plates. Eight fuel rods are removed to place two hollow water rods near the center. There is also a channel that surrounds the fuel rods and spacer grids from outside. The connections between the fuel rods and spacer grid are modeled by springs in the transverse plane. The channel sits on the lower tie plate and is attached to the upper tie plate by a post, and there are springs between the channel and the spacer grids. The total mass of the fuel assembly is 279.2 kg, in which the fuel rods account for 82.5%. In the FE model, the fuel rods are modeled by beam elements, same as in section 2.4.1, the spacer grids, water rods, and channel are modeled by shell element, and solid element model the lower and upper tie plates. The total number of DOFs is N = 1; 912; 506. The extremity of the bottom nose piece is fixed, and on each side of the channel, two nodes near the top are constrained to only rotation displacement. There are 5946 elastic modes in the interested frequency bandB = 2]0; 1000] rad/s, and 6605 elastic modes below 1200 Hz. The reference ROM uses the first n e = 6605 elastic modes as its basis. Figure 2.18 shows three local views of the fuel assembly model, while Figure 2.19 presents the fuel assembly, a global mode, and a local mode. We are interested in the global dynamics of this structure. For comparison purposes, two control FRFs are considered; for both cases, the external force is applied at one node of the upper part of the lower 40 Figure 2.18: Fuel assembly de- tails: view of the top (top), and views of the midspan (mid- dle) and of the bottom (bottom) without the channel thin faces Figure 2.19: Fuel assembly: undeformed configuration with and without the channel thin faces (left) and, without the channel thin faces, a global mode (center) and a local mode as seen from two points of view (right) tie plate in the x direction. The force is not applied at the center. Control FRF 1 is the acceleration in the x direction of one node in one adjacent face. Control FRF 2 is the acceleration in the y direction of the same observation node. The study in this section is part of the work on damage identification of the canister internals, especially the integrity of the fuel rods, by a non-intrusive approach. Hence, only the global behavior of the fuel assembly is concerned. The automatic mode selection procedure is applied to find dominant elastic modes. A subset of 240 DOFs based on 80 randomly selected nodes from the upper and lower tie plates are served as the excitation and observation pool to determine the importance of the elastic modes. Three global ROMs are constructed: the global ROMs, using piece-wise constant, piece-wise rigid-body, and combined kinematics. These global ROMs are based on a domain decomposition that the fuel assembly is uniformly cut into 100 transverse slices. For the global ROM based on combined kinematics, the flexible part is the fuel rods (similarly to section 2.4.1), and the kinematics of the water rods, channel, spacer grids, and stiff lower and upper tie plate structures is kept exact. Again, the reduced kinematics for each of the subdomains of the fuel rods is rigid-body. The global ROM by piece-wise constant kinematics has n g = 151 global modes inB. The prediction of control FRF 1 is shown in Figure 2.20, in which we can see that the dynamics of the response is more 41 significant than that of in section 2.4.1. The prediction accuracy is reasonable though it misses many resonance peaks in the low and middle frequency, and the amplitude in the high-frequency range is not well represented. For control FRF 2, see Figure 2.21, the error is even more significant, especially in the low-frequency range. 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact piecewise constant Figure 2.20: Control FRF 1 calculated with the global ROM of dimension n g = 151=6605 obtained by piece-wise constant kinematics(constantbyslice)andforwhich = 6:1 dB 0 200 400 600 800 1000 Frequency (Hz) -100 -80 -60 -40 -20 Accelerance (dB) exact piecewise constant Figure 2.21: Control FRF 2 calculated with the global ROM of dimension n g = 151=6605 obtained by piece-wise constant kinematics(constantbyslice)andforwhich = 11:6 dB The global ROM by piece-wise rigid-body kinematic has n g = 217 global modes inB. The predictions of control FRF 1 and FRF 2 are shown in Figures 2.22 and 2.23. We can see that the accuracy is improved for both cases, especially in the high frequency of the control FRF 1 and the low frequency of the control FRF 2. The prediction of FRF 2 is again less satisfactory than that of FRF 1. 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact piecewise rigid body Figure 2.22: Control FRF 1 calculated with the global ROM of dimension n g = 217=6605 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 4:1 dB 0 200 400 600 800 1000 Frequency (Hz) -100 -80 -60 -40 -20 Accelerance (dB) exact piecewise rigid body Figure 2.23: Control FRF 2 calculated with the global ROM of dimension n g = 217=6605 obtained by piece-wise rigid-body kinematics (rigid body by slice) and for which = 8:2 dB The global ROM by combined kinematic has n g = 681 global modes inB. The predictions of control FRF 1 and FRF 2 are shown in Figures 2.24 and 2.25, from which we can see that the accuracy for both cases is significantly improved. The prediction for FRF 1 is nearly perfect. The prediction error of control FRF 2 is still more significant; this is explainable since it has more complexities than FRF 1. 42 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact combined kinematics Figure 2.24: Control FRF 1 calculated with the global ROM of dimension n g = 681=5441 obtained by combined kinematics (rigid body by fuel rod slice and exact else- where) and for which = 4:2 dB 0 200 400 600 800 1000 Frequency (Hz) -100 -80 -60 -40 -20 Accelerance (dB) exact combined kinematics Figure 2.25: Control FRF 2 calculated with the global ROM of dimension n g = 681=6605 obtained by combined kinematics (rigid body by fuel rod slice and exact else- where) and for which = 3:7 dB The predictions given by the ROM of the 681 dominant elastic modes are shown in Figures 2.26 and 2.27, which are satisfactory everywhere except for an amplitude shift in FRF 1, which is due to the static contribution of the discarded modes. Regarding the numerical cost, the computation of the global ROM took about 3 hours and 20 minutes, whereas the computation of the elastic modes took about 11 hours. For both applications in sections 2.4.1 and 2.4.2, a quad-core CPU on a 64-GB machine is used. 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 Accelerance (dB) exact dominant modes Figure 2.26: Control FRF 1 calculated with the ROM of the n d = 681=6605 dominant elastic modes obtained by the automatic mode selection procedure and for which = 2:2 dB 0 200 400 600 800 1000 Frequency (Hz) -100 -80 -60 -40 -20 Accelerance (dB) exact dominant modes Figure 2.27: Control FRF 2 calculated with the ROM of the n d = 681=6605 dominant elastic modes obtained by the automatic mode selection procedure and for which = 2:6 dB Weseethateventhoughthemassoftheflexiblepartisdominant, wecanstillobtainagoodrepresentation of the global dynamics of the detailed fuel assembly by a small dimensional ROM. By introducing simple reduced kinematics, the global dynamics is decently approximated, and the irrelevant local vibrations are filtered. However, there is still a discrepancy due to the reduced kinematics of the flexible part, the fuel rods. A better global ROM would require better kinematics to represent the fuel rods. It is challenging for such a complex structure to develop such kinematics to represent the local vibrations coupled with the global resonances and not to represent the local vibrations independent of the global dynamics. 43 Chapter 3 Efficient modal analysis of the fully-loaded spent nuclear fuel canister This section proposes an efficient vibration analysis adapted to accurately describe the global dynamics of the sealed FLSNFC with high dimension and model density. The chapter is organized as follows. First, the fundamental tools for structural dynamics are described, including the SIL, CB, and domain decomposition. Then the efficient implementation of the CB technique adapted to the high modal density exhibited by the FLSNFC is introduced, followed by the methodology to efficiently determine the dominant sub-structural modes. Finally, the numerical results for the FLSNFC are presented. 3.1 Fundamental tools for structural dynamics 3.1.1 Shift-Invert Lanczos approach for eigenvalue analysis In section 2.1, the modal analysis of structural dynamics is introduced. In order to characterize the dynamics of the structure in the interested frequency bandB, the GEP Eq (2.2) is solved in a larger frequency band B c . For the problems when there are a large number of modes in the interested frequency band, an efficient strategy is to divide the frequency band into slices and use the Lanczos solver for the slices. A spectral shift is introduced to the GEP [K]' =[M]', which leads to a new GEP ([K][M])' = ()[M]'. Rewriting the GEP as a standard eigenvalue problem (SEP) ([K][M]) 1 [M]' = 1 ', then the largest 44 eigenvalues of this SEP correspond to the eigenvalues that are closest to the shift in the GEP. Introducing the symmetric indefinite matrix [H()] = [K][M], the GEP Eq (2.2) is replaced by the SEP [H()] 1 [M] ' = 1 ': (3.1) The Lanczos algorithm is an iterative method such that [H()] 1 [M] x has to be evaluated for variousx’s. To solve the linear systems by forward and backward substitutions, the LDL factorization [124] is used [P()] T [H()][P()] = [L()][D()][L()] T ; (3.2) where [L()] is a lower triangular matrix, [D()] is a block-diagonal matrix, and [P()] is a permutation matrix. In the following of the text, is omitted for simplified notation [P] T [H][P] = [L][D][L] T : For a given shift, the number of eigenpairs in its neighborhood must be calculated, which is done by an inertia count. The number of eigenvalues that have < for the GEP Eq (2.2) is the number of negative eigenvalues of the matrix [H], which is called the inertia count. By [125], it is equivalent to the number of the negative eigenvalues of the block diagonal matrix [D]. Performing inertia count for the whole frequency bandB, we intend to obtain many frequency slices with a homogeneous number of eigenvalues in each slice. Using SIL, the GEP can be solved in parallel with each run with a different shift . Moreover, the computation cost of the orthogonalization step in the Lanczos iterative procedure grows superlinearly with the number of eigenmodes, and this could be prohibitive when there are numerous modes. This problem is addressed by partitioning the frequency band, and the orthogonality between shifts is satisfied automatically. 3.1.2 Domain decomposition To take advantage of localized connections between different structural scales, the structural domain is partitioned intoN sub-structures 1 ; 2 ;:::; N . Then the FE matrices can be written as [K] = " K SS K SB K BS K BB # ; [M] = " M SS M SB M BS M BB # ; (3.3) 45 where the boundary DOFs that connect the sub-structures is denoted as B, and the interior DOFs of the sub-structures is denoted as S. The matrices [M SS ] and [K SS ] are block diagonal, for instance [K SS ] = 2 6 6 6 6 6 4 K 11 0 ::: 0 0 K 22 . . . . . . . . . . . . . . . 0 0 ::: 0 K NN 3 7 7 7 7 7 5 : (3.4) The GEP involves solving the linear system [K]x =f for many times. The matrix [K] is first decomposed as " K SS K SB K BS K BB # = " I N S 0 K BS K 1 SS I N B #" K SS 0 0 S K BB #" I N S K 1 SS K SB 0 I N B # ; (3.5) where N S and N B are the number of sub-structures DOFs and internal boundary DOFs, respectively, and [S K BB ] = [K BB ] [K BS ][K SS ] 1 [K SB ]; (3.6) is the Schur complement [126]. From Eq (3.5), the inverse of [K] can be obtained as [126] [K] 1 = " K SS K SB K BS K BB # 1 = " I N S K 1 SS K SB 0 I N B #" K 1 SS 0 0 S K BB 1 #" I N S 0 K BS K 1 SS I N B # : (3.7) Then, by partitionsx = (x S ;x B ) T andf = (f S ;f B ) T , the solution of [K]x =f can be written as, x B = [S K BB ] 1 f B [K BS ][K SS ] 1 f S ; x S = [K SS ] 1 (f S [K SB ]x B ): (3.8) In Eq (3.8), the linear system [S K BB ]x 0 =f 0 needs to be solved, which is cheap when the number of boundary DOFs, B, is small (which is the case for localized connections). Eq (3.8) also involves solving the linear system of sub-structures (index S) twice. Note that the matrix, [K SS ], for the sub-structures, is block diagonal, and hence, the linear system of each sub-structure can be solved independently. 46 By using SIL to solve the GEP Eq (2.2) as described in section 3.1.1, the linear system [H]x =f has to be solved for many times. Similar to Eq (3.8), the solution of [H]x =f can be obtained as x B = [S H BB ] 1 f B [H BS ][H SS ] 1 f S ; x S = [H SS ] 1 (f S [H SB ]x B ); (3.9) where [S H BB ] = [H BB ] [H BS ][H SS ] 1 [H SB ] is the Schur complement associated with [H]. A block LDL factorization of [H] can be obtained as [P] T [H][P] = [L][D][L] T , in which [L] = " L S 0 P T B H BS H 1 SS P S L S L B # ; [D] = " D S 0 0 D B # ; [P] = " P S 0 0 P B # ; (3.10) where [L S ], [D S ], [P S ], [L B ], [D B ]and [P B ]areobtainedbyLDLfactorization [P S ] T [H SS ][P S ] = [L S ][D S ][L S ] T and [P B ] T [S H BB ][P B ] = [L B ][D B ][L B ] T . The inertia countI([H]), that is the negative eigenvalues of [H] is obtained as I([H]) =I([D]) =I([D B ]) +I([D S ]) =I([D B ]) + N X J=1 I([D J ]): (3.11) 3.1.3 Craig-Bampton sub-structuring technique The CB sub-structuring technique uses the following displacement transformation ' = ' S ' B ! = " S t SB 0 I N B # q S q B ! = [V ]q; (3.12) whereq = (q S ;q B ) T is the component generalized coordinates; and [V ] = " S t SB 0 I N B # ; (3.13) is the CB transformation matrix with [t SB ] =[K SS ] 1 [K SB ], and [ S ] is the component modes of the sub-structures that satisfies [K SS ][ S ] = [M SS ][ S ][ S ]. 47 The component generalized coordinatesq is the solution of the GEP that projects the GEP of Eq (2.2) to the space spanned by the columns of [V ], that is [V ] T [K][V ] q = [V ] T [M][V ] q: (3.14) Taking advantage of the domain decomposition as described in section 3.1.2, matrices [ S ] and [ S ] are block diagonal [ S ] = 2 6 6 6 6 6 4 1 0 ::: 0 0 2 ::: . . . . . . . . . . . . 0 0 ::: 0 N 3 7 7 7 7 7 5 ; [ S ] = 2 6 6 6 6 6 4 1 0 ::: 0 0 2 ::: . . . . . . . . . . . . 0 0 ::: 0 N 3 7 7 7 7 7 5 : (3.15) Hence, [ S ] and [ S ] in the GEP [K SS ][ S ] = [M SS ][ S ][ S ] can be obtained by solvingN independent GEP associated with the sub-structures as follows [K JJ ][ J ] = [M JJ ][ J ][ J ]; forJ = 1;:::;N: (3.16) In each sub-structure J, only n J modes out of N J DOFs (usually n J N J ) are kept by considering a suitable cutoff frequency. We can write the CB GEP Eq (3.14) as [K]q =[M]q where [K] = [V ] T [K][V ] = " S 0 0 S K BB # ; [M] = [V ] T [M][V ] = " I N S M T BS M BS M BB # ; (3.17) with [M BS ] = M BS +t T SB [M SS ] [ S ]and [M BB ] = [M BB ]+[M BS ][t SB ]+([M BS ][t SB ]) T +[t SB ] T [M SS ][t SB ]. The dimension of the CB GEP is =N B +n S =N B + P N J=1 n J with n S = P N J=1 n J . 3.2 Craig-Bamptonimplementationsuitableforhighmodaldensity 3.2.1 Craig-Bampton implementation for Shift-Invert Lanczos Incomplexstructureswheretheassociatedsub-structureshavehighmodaldensity, thenumberofcomponent modes of sub-structures can be large, leading to a high dimensional CB GEP. In such a case, the direct 48 approach to solve the CB GEP (Eq (3.14)) is prohibitive by the high computational cost. However, from section3.1.3, theblockstructurefor [K]and [M]isinheritedfromthedomaindecomposition, inwhich [M BS ] is sparse with dense blocks (N B is small with localized the connections between sub-structures). Hence we can solve the CB GEP by SIL as described in sections 3.1.1 and 3.1.2. Similar to Eq (3.9), the linear system [H]x =f with [H] = [K][M] can be obtained by solving x B = [S H BB ] 1 f B [H BS ][H SS ] 1 f S ; x S = [H SS ] 1 (f S [H SB ]x B ); (3.18) where [S H BB ] = [H BB ] [H BS ][H SS ] 1 [H SB ] is the Schur complement associated with [H]; and some other matrices are of the following forms [H BB ] = [K BB ][M BB ]; [H BS ] =[M BS ]; [H SS ] = [ S ][I N S ]; (3.19) which leads to the expression of the Schur complement as [S H BB ] = [K BB ][M BB ] 2 [M BS ] ( S [I N S ]) 1 [M BS ] T : (3.20) Again, similar to section 3.1.2, a block LDL factorization [P ] T [H][P ] = [L][D][L] T can be obtained by the LDL factorization [P S ] T [H SS ][P S ] = [L S ][D S ][L S ] T and [P B ] T [S H BB ][P B ] = [L B ][D B ][L B ] T , in which [L S ] = [I N S ], [P S ] = [I N S ] and [D S ] = [ S ][I N S ], as [L] = " I S 0 P T B M BS ( S I N S ) 1 L B # ; [D] = " S I N S 0 0 D B # ; [P ] = " I S 0 0 P B # : (3.21) The inertia countI([H]), that is the negative eigenvalues of [H] is obtained as I([H]) =I([D]) =I([D B ]) +I([ S ][I N S ]) =I([D B ]) + N X J=1 I([ J ][I N J ]): (3.22) 49 3.2.2 Craig-Bampton technique adapted to the fully-loaded spent nuclear fuel canister In this work, we are interested in the vibration analysis of the FLSNFC. The canister is a cylindrical structure with a thick bottom plate, top lid, and a stiff honeycomb basket placed inside, and the FE model can be seen in Figure 3.1. The fuel assemblies (FAs) that follow GE 14 design are inserted in each of the 68 cells of the basket structure. The detailed FE model of the FA is shown in Figure 2.18. The details of the GE 14 FA are described in section 2.4.2. Figure 3.1: FE model of frame sub-structure (canister + basket), represented without the canister thick top lid As can be seen, the FLSNFC is a multiscale structure with three scales, namely, the lower scale, which consists of the fuel rods, the middle scale of the 68 FAs, and the upper scale of the canister and the basket. The connections between the FAs and basket is localized, in which the number of boundary DOFs associated with one FA is N = 27 ( denotes the boundary DOFs of the FA), which leads to the total number of boundary DOFs N B = 1140 (note that several FAs share some boundary DOFs). By taking advantage of the localized connections, a domain decomposition of 68 FAs and the structural frame (the canister and basket) is used, and the CB sub-structuring technique is applied. The total number of sub-structures is N = 69. A fine mesh is used for the FLSNFC, and the FE model has N = N B +N F + 68N A = 1140 + 2; 209; 084 + 68 1; 935; 837 = 133; 847; 140 DOFs, where N F and N A are the number of DOFs for the frame sub-structure and one FA, respectively. Due to nearly one hundred flexible fuel rods, about half a million modes are considered in the vibration analysis for each FA, which results in a CB model of the FLSNFC with high modal density. SIL solves this CB model following the procedure described in section 3.2.2. Note that the FAs are identical and connected to the frame sub-structure in the same manner; hence 50 we can first consider only two sub-structures, namely, the frame sub-structure and the FA, and then assemble the sub-structures on their boundaries. For the frame sub-structure, denoting the interior DOFs by “F” and the boundary DOFs by “B”, the FE matrices can be written as [K F ] = " K FF K FB K BF K F BB # ; [M F ] = " M FF M FB M BF M F BB # : (3.23) The static constraint mode of the boundary for the frame sub-structure can be computed by [t FB ] = [K FF ] 1 [K FB ]. The component modes of the frame sub-structure are obtained by [K FF ][ F ] = [M][ F ][ F ]: (3.24) Thenumberofmodesin [ F ]isdenotedasn F andisdeterminedbyacertaincutofffrequency. Byintroducing the CB transformation matrix [V F ] = " F t FB 0 I N B # ; (3.25) the CB matrices for the frame sub-structure can be obtained as [K F ] = [V F ] T [K F ][V F ] = " F 0 0 K F BB # ; [M F ] = [V F ] T [M F ][V F ] = " I F M T BF M BF M F BB # ; (3.26) where [K F BB ] = [K F BB ] + [K BF ][t FB ], and [M BF ] = [M BF ] + [t FB ] T [M FF ] [ F ], and [M F BB ] = [M F BB ] + [M BF ][t FB ]+([M BF ][t FB ]) T +[t FB ] T [M FF ][t FB ]. For the CB model of the FLSNFC, the Schur complement [S H BB ] as in Eq (3.20) is a combination of contributions from all its sub-structures. The contribution from the frame sub-structure is [S H;F BB ] = [K F BB ][M F BB ] 2 [M BF ] ([ F ][I N F ]) 1 [M BF ] T ; (3.27) where N F is the number of DOFs of the frame sub-structure. 51 For the FA, denoting the interior FA DOFs by “A” and the boundary DOFs by “”, the FE matrices can be written as [K A ] = " K AA K A K A K A # ; [M A ] = " M AA M A M A M A # : (3.28) The static constraint mode of the boundary for the FA can be computed by [t A ] =[K AA ] 1 [K A ]. The component modes of the FA are obtained by [K AA ][ A ] = [M][ A ][ A ]: (3.29) The number of modes in [ A ] is denoted as n A and is determined by certain cutoff frequency. Introducing the CB transformation matrix [V A ] = " A t A 0 I # ; (3.30) the CB matrices for the FA can be obtained as [K A ] = [V A ] T [K A ][V A ] = " A 0 0 K A # ; [M A ] = [V A ] T [M A ][V A ] = " I A M T A M A M A # ; (3.31) where [K A ] = [K A ]+[K A ][t A ], [M A ] = [M A ] + [t A ] T [M AA ] [ A ], and [M A ] = [M A ]+[M A ][t A ]+ ([M A ][t A ]) T + [t A ] T [M AA ][t A ]. The contribution of the FA to the Schur complement [S H BB ] is [S H;A ] = [K A ][M A ] 2 [M A ] ([ A ][I N A ]) 1 [M A ] T ; (3.32) where N A is the number of DOFs in the frame sub-structure. Combining Eqs (3.27) (3.32), the Schur complement of for the CB model of the FLSNFC can be obtained [S H BB ] = N X J=1 [S H;J BB ]; (3.33) where [S H;1 BB ] = [S H;F BB ]; and [S H;J BB ] is given by [S H;A BB ] for J 2. Note that the indices of the boundary DOFs of each FA are different in the frame sub-structure. In this study, we are interested in making a risk assessment of the internal integrity (lower scales) based on the signals detected from the exterior of the 52 canister. Let “I” denote the interested excitation and observation DOFs located at the frame sub-structure. Consequently, the (N I N I ) interested FRF matrix [U II (!)] associated withN I excitation and observation DOFs (the excitation and observation pool as defined in section 2.2) is obtained, through Eq (2.6) and the definition of the compliance described in section 2.2, as [U II (!)] = [ I F ] ! 2 [I ne ] + 2i![][ ] + [] 1 [ I F ] T ; (3.34) where the columns' I F of [ I F ] is given by ' I F = [ I F ]q F + [t I FB ]' B : (3.35) Recall that q = (q S ;q B ) T is the component generalized coordinates associated with the CB model of the FLSNFC, and here,q B =' B . In Eqs (3.34) and (3.35), [ I F ] and [t I FB ] are respectively the submatrices of [ F ] and [t FB ] restricting the DOFs to those of interest. 3.2.3 ErrorquantificationoftheCraig-Bamptonmodelofthefully-loadedspent nuclear fuel canister The error induced by the CB approximation needs to be quantified. The approximated FRF in the interested excitation and observation DOFs may be compared to the FRF obtained by exact elastic modes. However, the exact elastic modes are not available, and the purpose of the CB model is to avoid computing them. To address the issue, a reduced set of n d n e modes (n e is the number of exact elastic modes necessary to represent the FRF in the interested frequency band accurately) is used to compute the FRFs for comparison purposes. For the exact model, these n d modes are computed using SIL solver with domain decomposition as described in section 3.1.2. For the CB model of the FLSNFC, these corresponding n d modes are also computed. Note that the correspondence of the modes is crucial, which is achieved by inertia count for both the exact model and the CB model. Otherwise, the error quantification is not associated with the CB approximation. Since the FRFs of interest are located on the exterior of the canister, the majority of the elastic modes resulting from the vibrations of the flexible fuel rods (see in Figure 2.19) have a negligible 53 contribution to the FRFs. For the robustness of the comparison, we are choosing the n d most dominant modes by the automatic mode selection procedure described in section 2.2. We will recall the definition of global FRF error. Let U ij (!) be one of the FRF from the reference FRF matrix in Eq (3.34), and ~ U ij (!) be the approximation of U ij (!), the error in decibel is defined in Eq (2.12) ij =(u ij ; ~ u ij ) = s 1 B Z B (u ij (!) ~ u ij (!)) 2 d!: By using theN I interested DOFs as the excitation and observation pool, the global error (using the criterion of Eq (2.14)) is given by = +k ; = 1 N 2 I N I X i=1 N I X j=1 ij ; 2 = 1 N 2 I 1 N I X i=1 N I X j=1 ( ij ) 2 ; (3.36) in which k controls the confidence level of . The n e approximated CB modes are first computed. The importance of one given mode is determined by the global error defined in Eq (3.36) through sole modal removal, in which the reference FRF is obtained by all n e CB modes while the approximated FRF is obtained by including all n e CB modes except that particular mode. Then n d most dominant modes are computed exactly using SIL with frequency shifts corresponding to the CB eigenvalues. The n d modes are chosen to be computationally manageable, and the FRFs obtained by these modes have resonance peaks along with the frequency bandB. 3.3 Further reduced-order model by using dominant modes of sub- structures Due to the numerous local modes, the dimension of the CB GEP of Eq (3.14) can still be large, making the solving of the GEP be computational intense. In section 3.2.3, we mentioned that due to the local vibrations of the flexible fuel rods, there are numerous modes of the FLSNFC that have a negligible contribution to the interested FRFs, and thus the removal of these modes will not affect the response. This filtering is carried out on the structure’s upper scale after the CB modes’ computation. A ROM is proposed hereinafter, in 54 which the strategy is to filter the local modes of the sub-structures, specifically the modes of the FAs, in the process of building the CB model of the FLSNFC, to reduce n S = P N J=1 n J . To still have an accurate approximation, the filtered modes should have a negligible contribution to the FRFs, in which the importance of the modes of the sub-structures is determined following the same procedure described in section 3.2.3. The importance of one mode in the sub-structure is quantified by the global FRF error induced by its sole removal. The error is regarding the CB model of the FLNSNFC, and we want to avoid solving the GEP of the CB model. The reference FRFs are obtained by the CB model of FLSNFC through a direct numerical approach, which is supposed to be efficient with sparse frequency sampling and a reduced set of interested excitation and observation DOFs. To solve the CB model by direct approaches, the hysteretic damping model [127] is used because the eigenpairs of the structure are not available (this is the purpose of solving GEP of the CB model), which leads to the following CB matrices [K] = " S 0 0 S K BB # ; [M] = " I N S M T BS M BS M BB # ; [D] = 2(!) ! [K]; (3.37) where [D] is the hysteretic damping matrix with (!) the frequency-dependent damping ratio and the term 2(!) is usually called the loss factor. The hysteretic damping model matches the modal damping at ! if (!) is given by (!) = [127], where ! and are the eigenfrequency and damping ratio of mode , respectively. For completeness, [M BS ] = M BS +t T SB [M SS ] [ S ] and [M BB ] = [M BB ] + [M BS ][t SB ] + ([M BS ][t SB ]) T + [t SB ] T [M SS ][t SB ]. The CB dynamic stiffness matrix [H(!)] =! 2 [M] + i![D] + [K] is frequency-dependent, and can be obtained as [H(!)] = " H SS (!) H T BS (!) H BS (!) H BB (!); # ; (3.38) where [H SS (!)] =! 2 [I N S ] + (1 + 2i(!)) [ S ]; [H BS (!)] =! 2 [M BS ]; [H BB (!)] =! 2 [M BB ] + (1 + 2i(!))[S K BB ]: 55 The displacement transformation of CB ' = [V ]q as in Eq (3.12) is translated to [U(!)] = [V ]q(!). The FRFs of the frame sub-structure are obtained as [U(!)] = U S (!) U B (!) ! = [V ]q = " S t SB 0 I N B # q S q B ! ; (3.39) yielding [U S (!)] = [ S ]q S + [t SB ]q B and [U F (!)] = [ F ]q F + [t FB ]q B . Similar to section 3.1.2, q F and q B can be obtained using the Schur complement. Let [S H BB (!)] be the Schur complement associated with [H(!)] such that [S H BB (!)] = [H BB (!)] [H BS (!)][H SS (!)] 1 [H BS (!)] T ; (3.40) then, under the external force ofF(!), q B (!) = [S H BB (!)] 1 f B (!) [H BF (!)][H FF (!)] 1 f F (!) ; q F (!) = [H FF (!)] 1 f F (!) [H BF (!)] T q B (!) ; (3.41) where f(!) = [V ] T F(!). Plugging the second row of Eq (3.41) into Eq (3.40), we can obtain [U F (!)] = [ F ][H FF (!)] 1 f F (!) + [t FB ] [ F ][H FF (!)] 1 [H BF ] T q B (!): (3.42) For the case of unit point loads associated with the FRF matrix [U II (!)], in which f F = [ I F ] T and f B = [t I FB ] T , one can obtain [U II (!)] by Eqs (3.41) and (3.42) as [U II (!)] = [L II (!)] + [L IB (!)][S H BB (!)] 1 [L IB (!)] T ; (3.43) where [L II (!)] = [ I F ][H FF (!)] 1 [ I F ] T ; [L IB (!)] = [t I FB ] [ I F ][H FF (!)] 1 [H BF ] T = [t I FB ] +! 2 [ I F ][H FF (!)] 1 [M BF ] T : 56 The Schur complement, in this case, can be obtained by Eqs (3.40) and (3.38) as [S H BB (!)] = [H BB (!)] [H BF (!)][H FF (!)] 1 [H BF (!)] T N=69 X J=2 [H BJ (!)][H JJ (!)] 1 [H BJ (!)] T = [H BB (!)]! 4 [M BF ][H FF (!)] 1 [M BF ] T + N=69 X J=2 [M BJ ][H JJ (!)] 1 [M BJ ] T ! = [H BB (!)]! 4 h s H;F BB (!) i + N=69 X J=2 h s H;J BB (!) i ! ; (3.44) where h s H;F BB (!) i = [M BF ][H FF (!)] 1 [M BF ] T ; and the nonzero terms of h s H;J BB (!) i is given by h s H;A i = [M A ][H AA (!)] 1 [M A ] T , for J = 2;:::;N. Let [ ~ U II (!)] be the approximation of [U II (!)] due to the removal of a mode from one given FA. Since [L II (!)] and [L IB (!)] is retained after approximation, [ ~ U II (!)] can be obtained as [ ~ U II (!)] = [L II (!)] + [L IB (!)][ ~ S H BB (!)] 1 [L IB (!)] T ; (3.45) where [ ~ S H BB (!)] is given by the following modification of [S H BB (!)] as [ ~ S H BB (!)] =[H BB (!)]! 4 [M BF ][H FF (!)] 1 [M BF ] T ! 4 N=69 X J=2 [M BJ ][H JJ (!)] 1 [M BJ ] T 1 h(!) v B v T B ! =[S H BB (!)] + ! 4 h(!) v B v T B ; (3.46) in whichh(!) =! 2 + (1 + 2i(!)) A is the component of [H AA (!)] with A the eigenvalue of the removed FA mode, and v B is the removed column of [M BA ] due to the removal of one FA mode. It can be seen from Eq (3.46) that [ ~ S H BB (!)] is the rank-1 modification of [S H BB (!)], and its inversion can be given by the Woodbury formula [128] as [ ~ S H BB !] 1 = [S H BB (!)] 1 [S H BB (!)] 1 v B h(!) ! 4 +v T B [S H BB (!)] 1 v B 1 v T B [S H BB (!)] 1 : (3.47) 57 Plugging Eq (3.47) into Eq (3.45) we can obtain [ ~ U II (!)] =[L II (!)] + [L IB (!)][S H BB (!)] 1 [L IB (!)] T [L IB (!)][S H BB (!)] 1 v B h(!) ! 4 +v T B [S H BB (!)] 1 v B 1 v T B [S H BB (!)] 1 [L IB (!)] T =[U II (!)] 1 a(!) w I (!)w T I (!); (3.48) where a(!) = h(!) ! 4 +v T B [S H BB (!)] 1 v B and w I (!) = [L IB (!)][S H BB (!)] 1 v B are defined for brevity expres- sion. The dimension N I and the frequency sampling can be adjusted to reduce the computational cost for Eqs (3.43) and (3.48). The FA modes can then be sorted by their importance measure, and truncation is made. A CB model of the FLSNFC with a smaller dimension is obtained due to the truncation, enabling an efficient computation of the associated GEP. In addition, the filtering of the numerous modes of the sub-structures will also lead to the filtering of system modes, which results in a reduced set ofn<n e modes that can accurately describe the global dynamics of the FLSNFC. 3.4 Application on the fully-loaded spent nuclear fuel canister In vibration analysis of the FLSNFC, the interested FRFs are within the frequency band ofB = 2]0;f u ], where f u = 1000 Hz. To take the spread contributions of the eigenmodes into consideration, the modes withinB c = 2]0;f c ] with f c = 1200 Hz are computed. The total number of system modes withinB c is n e = 458; 910. The objective is to obtain the FRF matrix [U II (!)] of Eq (3.34), which has the FRFs of all the combination of N I = 1581 DOFs located at the exterior of the canister bottom plate, see in Figure 3.2. These excitation and observation nodes are chosen for damage identification of the internal components. The FLSNFC is suspended at the top through two points, which are assumed to be fixed in all DOFs, see in Figure 3.3. The excitation and observation are applied at the exterior of the canister bottom plate, see in Figure 3.4. The FAs are connected to the frame sub-structure at the bottom through a slender nose piece and vertical and horizontal springs near the top. The vibrations of the FAs are thus transmitted to the bottom plate through the nose piece and to the basket through the springs. 58 Figure 3.2: Location of the nodes associated withN I = 1581 DOFs of interest on the exterior of the canister bottom plate Figure 3.3: Top of a longitudinal mid- section of the canister loaded with one fuel assembly. The structure is fixed at the top of the canister lid through 2 attachments. Figure 3.4: Bottom of a longitudi- nal mid-section of the canister loaded with one fuel assembly. The exter- nal loads (red) and dynamic response monitoring (blue) are restricted to the exterior of the canister bottom plate. 3.4.1 Reference Craig-Bampton model In order to have a model for accuracy certification, a reference CB model corresponding to the cutoff fre- quency at 20 kHz for the frame and FA sub-structures is proposed. Accordingly, the following preliminary computations are carried out first: • Component modes [ F ] (N F n F ) of Eq (3.24) with 20 kHz truncation, obtained in 60 hours with 180 compute nodes; 59 • Component modes [ A ] (N A n A ) of Eq (3.29) with 20 kHz truncation, obtained in 5 hours with 180 compute nodes; • Static constrained modes [t FB ] (N F n B ), obtained in 3.5 hours with 1 compute node; • Static constrained modes [t A ] ((N A n )), obtained in 15 minutes with 1 compute node; where N F = 2; 209; 084, N B = 1140, n F = 211; 057, N A = 1; 935; 837, N = 27, and n A = 46; 383. The compute nodes are of 64 GB of memory and 16 CPUs if not specified. The component modes [ F ] and [ A ] are computed with LS-DYNA [129] by spectrum slicing, while the static constrained modes [t FB ] and [t A ] are computed with MATLAB [130]. In the remainder of this chapter, the computations are carried out with MATLAB software. The reference CB model of dimension = N B +n F + 68n A = 3; 366; 241 can be solved for the first n e = 458; 910 modes below f c in 3.8 hours with 50 compute nodes. Then, the CB modes are removed one at a time, and the induced global error in Eq (3.36) is calculated with k = 4, in which, for efficient consideration, a reduced set of N I = 63 DOFs (see in Figure 3.5) is used as excitation and observation pool. The value k is chosen such that nearly 99% of the approximated FRFs in the excitation and observation pool have an error less than . The CB modes are sorted by their importance, quantified by , and the n d = 200 most dominant modes are selected for accuracy certification. The corresponding 200 modes are recomputed exactly without any approximation by SIL solver, in which, for each of the n d = 200 modes, an eigenvalue shifts = ( is the CB eigenvalue) is introduced, and the nearby elastic modes are computed. The inertia counts of both the CB model and the exact model ensure the mode correspondence of these 200 modes. Though the localized connection between sub-structures is utilized, the computation of the 200 exact modes is expensive, which takes 35 hours with 25 compute nodes of 128 GB memory and 24 CPUs. The FRF matrices [U II (!)] computed by using the 200 dominant modes from the CB model and the exact model are obtained, in which the FRFs associated with the N I = 1581 DOFs as shown in Figure 3.2. The associated error ij of the FRFs obtained by the CB model is computed, and the histogram and global error of which are shown in Figure 3.6. An example FRF of error approximated equals is shown in Figure 3.7, which shows that the error is very mild. Note that more than 98.9% of the FRFs have an error less than ; this implies that the CB model with 20 kHz component truncation is accurate and can serve as a 60 Figure 3.5: Location of the nodes associated with N I = 63 DOFs of interest on the exterior of the canister bottom plate reference model. Several approximations are introduced in the following to increase the efficiency further and are compared with the reference CB model. Figure 3.6: Histogram of the error ij for the 20 kHz Craig-Bampton model with respect to the exact modes, both restricted to the n d = 200 most dominant modes. Red line: global error = + 4 = 0:46 dB. 0 200 400 600 800 1000 Frequency (Hz) -140 -120 -100 -80 -60 Acceleration (dB) Figure 3.7: FRF comparison between the 20 kHz Craig-Bampton model (red) and the exact model (black), both restricted to the n d = 200 most dominant modes. Error: ij = 0:46 dB. 3.4.2 Two-level nested Craig-Bampton model Relaxed truncation of sub-structural modes. To have a CB model with relaxed truncation of sub- structural modes and quantify the associated errors, convergence analysis of the global error to the cutoff frequency of the sub-structures is carried out. The CB model with cutoff frequencies at 1.2 kHz, 2 kHz, 2.4 kHz, 2.8 kHz, 3.2 kHz, 3.6 kHz, 4 kHz, 5 kHz, 8kHz, and 12 kHz are computed, the associated global errors with k = 4 are plotted in Figure 3.8. It can be seen that the error decreases with the increase of cutoff frequency and will converge to zero; moreover, 20 kHz cutoff frequency is not necessary. The error for the case with a 4 kHz cutoff frequency is = 0:38 dB, which is smaller than the error of the reference CB model. 61 For the case of 4 kHz truncation, the histogram of the error ij for the FRFs of the N I = 1581 DOFs of interest is plotted in Figure 3.9, one example FRF comparison with the error that is representative of the global error = 0:38 dB is plotted in Figure 3.10. Note that more than 99% of the FRFs have less error than 0.38 dB; thus, the 4 kHz truncation CB model does not lose the accuracy significantly and is used for future applications. For this truncation, we have n F = 16; 623 andn A = 14; 064 which leads to a CB model with dimension of = 974; 115. The CB GEP of this model is solved for the first n e modes in about 118 minutes with 68 compute nodes. 2 3 4 5 6 7 8 9 10 11 12 Cutoff frequency (kHz) 0 2 4 6 8 Error (dB) Figure 3.8: Convergence of the Craig-Bampton model error = + 4 with respect to the cutoff frequency. Figure 3.9: Histogram of the 4 kHz Craig- Bampton model error ij . Global error = + 4 = 0:38 dB (red line). 0 200 400 600 800 1000 Frequency (Hz) -140 -120 -100 -80 -60 Acceleration (dB) Figure 3.10: FRF comparison between the 4 kHz Craig-Bampton model (red) and the 20 kHz reference Craig-Bampton (black). Error: ij = 0:38 dB. Inner Craig-Bampton model. Since the objective of the study is to detect internal damages from signals on the exterior of the canister, uncertainties are to be introduced to the FAs, and the eigenvalue problem Eq (3.29) of the FA has to be solved many times. To efficiently solve the FA model, a CB model is also introduced. The CB model for FA is referred to as “inner CB,” and the previous system-level CB model is called “outer CB”. Note that the connections of the subcomponents of the FA are also localized and can be 62 utilized for efficient calculation (see also in section 2.4.2 for details). The domain decomposition of the FA is shown in Figure 3.11, in which the left figure is in full length, the middle and right figures are the zoom-in figures of the top and bottom parts. It can be seen that the FA is partitioned into, from left to right and top to bottom, the channel, the two water rods, the upper tie plate, the fuel rods that partitioned into nine segments along with the height, the lower tie plate, and the extremity of the bottom nose piece. Similar Figure 3.11: Fuel assembly domain decomposition: all the sub-structures are represented as disconnected from the others. to the analysis of the different cutoff frequencies in the outer CB model, convergence analysis of the global error to the cutoff frequency of the subcomponents is carried out. It is discovered that an inner truncation of 6 kHz is sufficient to keep the same level of accuracy. The histogram of the error ij and one example of FRF comparison that has the same error as the global error ( = 0:41 dB) are presented in Figures 3.12 and 3.13. More than 99% of the FRFs have less error than 0.41 dB; thus, the 6 kHz truncation of the inner CB model does not lose the accuracy significantly. The eigenvalue problem of the inner CB is solved in about 6 minutes by one compute node. The combination of the outer and inner CB is referred to as the two-level nested CB model. It takes about 118 + 6 = 124 minutes to compute the GEP of the two-level nested CB model. 63 Figure 3.12: Histogram of the error ij of the two-level nested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner trunca- tion. Global error = + 4 = 0:41 dB (red line). 0 100 200 300 400 500 600 700 800 900 1000 Frequency (Hz) -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 Acceleration (dB) Figure 3.13: FRF comparison between the 20 kHz reference Craig-Bampton (black) and the two-level nested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner truncation (red). Error: ij = 0:41 dB. 3.4.3 Reduced-order model based on dominant fuel assembly modes Themethodologyproposedinsection3.3isapplied. Theimportanceofthen F = 14; 047FAmodescomputed by the inner CB model with 6 kHz truncation has to be quantified for all 68 FAs. The importance of one specific FA mode is measured by the global error induced by the sole removal of the mode. The algorithm in section 3.3 is used to obtain the FRF matrix [U II (!)] and its approximation [ ~ U II (!)] and the global error can be obtained by the Eq (3.36), in which N I = 63 DOFs of interest presented in Figure 3.5 is used as excitation and observation pool, and a coarse frequency sampling is considered (instead of using n ! = 3000, n ! = 1000 is used as the number of frequency point in the bandB). The importance of the FA modes below 4 kHz is sorted and presented in Figures 3.14 and 3.15. It can be seen from Figure 3.15 that the modal 0 1 2 3 4 5 6 7 8 9 10 Mode index ×10 5 10 -10 10 -5 10 0 Modal importance (dB) Figure3.14: Importancevalue(dB)ofallthefuel assembly modes (sorted). Figure 3.15: For each fuel assembly, importance value (dB) of its modes (sorted). importance is slightly different for different FAs, and this is because of their different locations in the canister. To determine how many FA modes are to be filtered and remain the same level of accuracy, convergence analysis of the global error to the percentage of the FA modes preserved is carried out; the result is shown 64 in Figure (3.16). It can be seen from Figure (3.16), with 25% of the most dominant modes retained for each FA, the induced global error = 0:45 dB (following the excitation and observation pool ofN I = 1581 DOFs) is of the same level as that of the two-level nested CB model. In this case (n A = 3512), the histogram of the error ij is plotted in Figure 3.17, and one example FRF comparison with the same error as the global error is shown in Figure 3.18. As can be seen, this ROM still has the same level of accuracy as previous 10 20 30 40 50 60 70 80 90 100 Dimension (%) 0 1 2 3 4 Error (dB) Figure 3.16: Convergence of the global error = + 4 with respect to the percentage of dominant fuel assembly modes preserved. Figure 3.17: Histogram of the error ij for the proposed reduced-order model with 75% removal of the fuel assembly modes below 4 kHz fre- quency. Global error = + 4 = 0:45 dB. 0 100 200 300 400 500 600 700 800 900 1000 Frequency (Hz) -110 -100 -90 -80 -70 -60 -50 -40 Acceleration (dB) Figure 3.18: FRF comparison between the 20 kHz Craig-Bampton reference (black) and the proposed reduced-order model with 75% removal ofthefuelassemblymodesbelow4kHzfrequency (red). Error: ij = 0:45 dB. approximations. With 75% of the FA modes are removed, the outer CB model has a smaller dimension of = 256; 579, and as a result, the associated GEP can be solved more efficiently, in about 21 minutes with 68 compute nodes (instead of 118 minutes). The time to determine the importance of the FA modes for all 68 FAs is 14 minutes with 68 compute nodes. The filtering of the FA modes leads to the filtering of the system mode in which n = 232; 586 modes are presented in the frequency band instead of n e = 458; 910 modes. To highlight the necessity of evaluating the contribution of the FA modes, Figure 3.19 plots the histogram of the error ij associated with removing the most dominant FA mode from the calculation. The global error 65 of which is = 7:43 dB. An example FRF comparison with the error of 7.43 dB is presented in Figure 3.20, which shows that a single removal of the sub-structural mode can induce a significant error in the system level FRF. Figure 3.19: Histogram of the error ij for the two-level nested Craig-Bampton model removal of the one most dominant fuel assembly mode. Global error = + 4 = 7:43 dB. 0 200 400 600 800 1000 Frequency (Hz) -150 -100 -50 Acceleration (dB) Figure 3.20: FRF comparison between the 20 kHz Craig-Bampton reference (black) and the two-level nested Craig-Bampton model with 4 kHz outer truncation and 6 kHz inner truncation from which the one most dominant fuel assembly mode is removed (red). Error: ij = 7:43 dB. To conclude this chapter, a two-level nested CB model with dominant sub-structural modes in the CB transformation matrix is proposed. A computational gain of four orders of magnitude is obtained by comparing the cost per core to the exact model at the expense of negligible errors that are exactly characterized. In addition, the inner CB of the FA enables one to have efficient reanalysis when there are localized damages (or uncertainties) in the FAs. 66 Chapter 4 Accelerated basis adaptation method in the framework of polynomial chaos expansion In this chapter, we proposed several strategies to accelerate the basis adaptation in homogeneous chaos spaces. The accelerated basis adaptation method can perform the forward UQ efficiently by reducing the dimension of the random input space. The basis adaptation in the framework of PCE is first recalled, and then the methods to accelerate the convergence of the classical basis adaptation method are proposed. The final part is the applications of the methodology to a borehole function and a space structure model. 4.1 Basis adaptation in polynomial chaos expansion 4.1.1 Polynomial chaos expansion Consider physics problems with input parameters!2R d specified as random variables and possibly jointly distributed. Also, consider a QoI, Q, that is typically low-dimensional and can be expressed in the form of Q =M(!). Numerical resolution of the physical modelM(!) will typically result in an enormous scale computational problem that renders the evaluation of Q as a function of ! prohibitively expensive. As a first step in the PCE-based methods, input parameters are mapped to independent random variables that characterize the independent causes of observed uncertainty in !. We thus introduce an operatorT such that! =T (). This operator can be represented, for instance, in the form of a Rosenblatt transform [131] 67 or other measure transport operators [132]. We can thus express the quantity of interest, Q, in the following form, Q() =M(!) =M(T ()): (4.1) The second step in the PCE is to express the above composite map as a function of , using multivariate approximations that take into account the weight presented by the probability density function of . One reason to rewrite the problem in terms of independent vector is that the associated probability measure is the product of one-dimensional measures, thus allowing a systematic approach to this approximation task. This thesis is interested only in where is a vector of independent standard Gaussian variables. Consider a probability space ( ;F;P), andH a d-dimensional Hilbert space which is a closed linear span of a random vector that has standard multivariate Gaussian distribution. LetF(H ) be the - field generated by 2 H and H :n: be the homogeneous chaos of exact order n, then the space L 2 = L 2 ( ;F(H );P) can be decomposed as L 2 = L 1 n=0 H :n: [133], which is called the chaos decomposition. In addition,H :n: isspannedbytheordernHermitepolynomialh n [133], wherefh n g 1 n=1 isthesetoforthogonal polynomials with respect to the standard Gaussian measure. Thenf := h =kh k 2 2 ;jj = n 0g is an orthonormal basis in L 2 , where the multi-index = ( 1 ;:::; d )2J d := (N 0 ) d andjj = P d i=1 i . We can then represent the QoI Q2L 2 by the PCE [14], Q() = X 2J d Q (); (4.2) where Q is the PC coefficient associated with the Hermite polynomial . LetJ d;p J d be the subset of indices of maximum length p (jjp). Restricting the PCE to an arbitrary subsetJ ofJ d , we obtain the following approximation to Q(), Q()Q J () = X 2J Q (): (4.3) It can be shown that the approximation Q J d;p ! Q as p!1 in L 2 norm. The number of PC terms in Eq (4.3) is P = d +p p , which grows factorially with respect to both d and p. 68 In the non-intrusive implementation, PCE coefficients in Eq (4.3) can be found independently by pro- jectingM(!) tof g, Q = hM; i h 2 i = hQ; i h 2 i =hQ; i = 1 (2) d=2 Z R d Q(s) (s)e 1 2 s T s ds; 2J; (4.4) in which the inner product defined in L 2 space ishu;vi =hu;vi L 2 =E[uv] for any u;v2L 2 , whereE[] is the expectation with respect to the standard Gaussian measure. By this orthogonal projection, the mean- squared error QQ J 2 2 is minimized. Eq (4.4) can be further approximated by numerical quadrature as Q = X ( (q) ;wq )2Q d Q (q) (q) w q ; 2J ; (4.5) where the summation is over pairs (q) ;w q of coordinates and weights that define ad-dimensional quadra- ture ruleQ d . The number of points inQ d is denoted as n d . Given the Gaussian measure associated with, d-dimensional Gauss-Hermite quadrature has been a common choice for evaluating the integral in Eq (4.5). The main computation cost relies on evaluating the physical modelM T (q) for n q times. With the full tensor product quadrature rules, the number of required quadrature points grows exponentially with dimensiond [56]. The sparse quadrature method proposed by Smolyak [59] offers an efficient way to temper the curse of dimensionality of the full quadrature rules. It utilizes a level parameter to control the accuracy of the integration, and then, with a specified level, the dimension of the function determines the number of required quadrature points. However, the number of quadrature points still grows superlinearly with the increase of dimension and, hence, remains prohibitive for functions with large dimensions or that are themselves too expensive to evaluate more than a modest number of times. 4.1.2 Basis adaptation and dimension reduction Basis adaptation. The main idea of basis adaptation is to rotate the Gaussian input variables properly to form a new basis such that the PCE in terms of the new basis has concentrated representation in the first several dimensions [6]. 69 Let matrixA2R dd be an isometry such thatAA T =I, and2R d be defined as =A; = 2 4 r :r 3 5 ; (4.6) where r denotes the vector of the first r components, and :r denotes the vector of the remaining (dr) components of . It is easy to see that also has a standard multivariate Gaussian distribution. The expansion of ^ Q in terms of can be written as Q()Q J () =Q J (A T ) : (4.7) This gives rise to two representations as follows, Q J () = X 2J Q (); (4.8) and Q J (A T ) = X 2J Q A () = X 2J Q A A (); (4.9) where A () = (A). If the indexing set is equal toJ d;p , then the sets of polynomialsf g andf A g over this indexing set span the same subspace of d-dimensional polynomials of total order p [133], and we have, X 2J d;p Q () = X 2J d;p Q A A () : (4.10) We then have, Q A = X 2J d;p Q h ; A i; Q = X 2J d;p Q A h A ; i; 2J d;p : (4.11) The analytical results in equation (4.11) will be used to assess and correct errors introduced by approxima- tions in the remainder of the paper. 70 Note that Eqs (4.6) to (4.11) hold for any isometry A, and different choices of A, for instance, the Gaussian and quadratic adaptation [6], result in different performances of how the QoI is adapted to the rotated basis. For the Gaussian adaptation case, a pilot PCE of order up to 1 in the original space (full- dimensional) is first built, then the first row of the rotation matrixA is defined such that 1 = X 2J d;1 Q () = d X i=1 Q ei i ; (4.12) where e i is a d-dimensional multi-index with 1 at the i-th location and zeros elsewhere. ThenfQ ei g d i=1 is the set of PCE coefficients associated withf i g d i=1 , that is, the first order PCE coefficients or Gaussian coefficients. By this particular way of construction, 1 captures the complete Gaussian components of Q. The remaining rows of A are constructed in two steps as follows. The entries offQ ei ; i = 1; ;dg are ordered by descending order of their absolute value asfq i ; i = 1; ;dg, and the indices ofq i in the original fQ ei g are stored in an array i . Each of the rows of A starting with the second, (i = 2 ;d), is set to zero except at the location corresponding to i1 , which is set to the value Q e i1 . Thus, the second row is zero except at the location of the most important variable, which is set to q 1 , and the third row is zero except at the location corresponding to the second most important variable, which is set to q 2 , and so on. In a second step, a Gram-Schmidt procedure is applied to matrix A, with its first row being normalized, to transform it into a rotation matrix [6]. In this manner, the first r rows of A span the vectors 1 and those (r 1) coordinates of with the largest contributions to the first order PCE. The procedure to build a rotation matrix is summarized in Algorithm 4.1. This construction provides the motivation for seeking a reduced representation of Q in terms of r . Dimensionreduction. DenoteA r thematrixobtainedbytakingthefirstr rowsofA, andA :r thematrix consisting of the remaining (dr) rows. A reduction is achieved by expandingQ in a lower dimensional basis f Ar g 2Jr;p as Q Ar ( r ) =Q A 0 @ 2 4 r 0 3 5 1 A = X 2J d;p Q A 0 @ 2 4 r 0 3 5 1 A = X 2Jr;p Q Ar ( r ); (4.13) 71 Algorithm 4.1: Construction of rotation matrix by Gaussian adaptation 1 Construct first order PCE Q() =Q 0 + P d i=1 Q ei i with where Q 0 the 0 order coefficient ande i the d-dimensional multi-index with 1 at i-th location and zeros elsewhere; 2 Compute Q 0 andfQ ei g i=1;:::;d (e.g. by Eq (4.4)); 3 Construct the first row of rotation matrixA2R dd such that Eq (4.12) is satisfied; 4 Rank the first order coefficients by absolute value in descending order and record their indices in the original coordinates asf j g d j=1 ; 5 for j = 2 to d do 6 Intialize the the j th row ofA to be zeros and let the j1 -th entry be qual to Q e j1 ; 7 Perform Gram-Schmidt procedure onA to make it an isometry. where r d is the dimension of the reduced space;J r;p J d;p is the set of multi-indices that only has nonzero entries associated with r ;A r is the submatrix ofA with only the first r rows ofA; the subscript A r on Q indicates that the reduction to r dimensions is achieved viaA r ; and Ar := ( r ) = (A r ). In the sequel, Eq (4.13) is called the “r-dimensional adapted PCE of Q”, or simply the “r-dimensional adaptation of Q”. Note that in Eq (4.13), the assumption :r = 0 is made, the error introduced by which is negligible if the restriction of the approximation to the polynomial span of r is of sufficient accuracy. The PCE coefficientsfQ Ar g can be obtained by projectingQ tof Ar g and further be approximated by the r-dimensional Gauss-Hermit quadrature as Q Ar =hQ; Ar i = X ( (q) r ;wq )2Qr Q r (q) r (q) r w r q ; 2J r;p ; (4.14) where the quadrature ruleQ r is now r-dimensional, and r (q) r =A 1 2 4 (q) r 0 3 5 =A T 2 4 (q) r 0 3 5 (4.15) is the transformation from r to. As already mentioned, we note that the inverse mapping from to in Eq (4.15) omits the contributions of :r which introduces errors to the values of the PCE coefficients, but these errors are small when most of the probabilistic information is captured by the first r components of . As previously mentioned, the cost of the non-intrusive implementation of PCE grows linearly with the 72 number of physical model evaluations, which is the number of quadrature points. Thus, if we assure that r d, PCE with basis adaptation would require significantly fewer quadrature points than the standard PCE. To quantify the error of expansion Eq (4.13), we first projectQ Ar to the basisf g 2J d;p in the original space as Q Ar ( r ) = X 2J d;p ~ Q (); (4.16) and ~ Q can be obtained (by Eqs (4.13) and (4.16)) as ~ Q = X 2Jr;p Q Ar h Ar ; i: (4.17) The “~” over the PCE coefficient indicates that the coefficient is the projection from one space to another and will be referred to as “projected coefficient”. By doing so, the PCE coefficientsfQ Ar g 2Jr;p in the rotated space are transformed to the PCE coefficientsf ~ Q g 2J d;p in the original space. Define the vectors of PCE coefficients ~ Q coeff :=f ~ Q g 2J d;p andQ coeff :=fQ g 2J d;p , then ~ Q coeff is an approximation ofQ coeff , and the accuracy is dependent on r andA. The relative 2-norm between these two PCE representations can be measured by D = Q coeff ~ Q coeff 2 kQ coeff k 2 : (4.18) The error measure in (4.18) quantifies the error of coefficients relative to ad-dimensional basis, and we hence refer to this error as the “error in the full-dimensional space”. A systematic way to carry out the basis adaptation approach starts from a 1d adaptation in (4.13) (with r = 1); the value of r is then incremented sequentially, with a stopping criterion related to the proximity between two successive approximations. In order to compare two successive adaptations, say, r- and (r + 1)- dimensional adaptation, we can either carry a comparison of the probability density functions (PDFs) of the QoI or examine the difference between the vectors of PCE coefficients. In the latter case, we can first transfer PCE coefficients of bothfQ Ar g 2Jr;p andfQ Ar+1 g 2Jr+1;p to the original space in the same manner as in Eq (4.16), and then compare the two resulting vectors. A more straightforward way that compares PCE 73 coefficients in the space off Ar+1 ;2J r+1;p g is also possible by noting the hierarchical nesting of the multivariate Hermite polynomials implying that polynomials of dimension (r + 1) include as a subset of the polynomials of dimension r. Thus, Q Ar can be first projected to the space off Ar+1 ;2J r+1;p g where the projected coefficients are denoted by ~ Q Ar . To complete the notation, we introduce the vector of PCE coefficients of the (r+1)-dimensional expansion,Q Ar+1 coeff :=fQ Ar+1 ;2J r+1;p g, and the vector of projected PCE coefficients of the r-dimensional expansion, ~ Q Ar coeff :=f ~ Q Ar ;2J r+1;p g, which allows us to express the “error in the reduced dimensional space” in the form, R = k ~ Q Ar coeff Q Ar+1 coeff k 2 kQ Ar+1 coeff k 2 : (4.19) The convergence in the coefficients means that the additional basis termsf Ar+1 ;2J r+1;p and = 2J r;p g have negligible contributions toQ. This error measure and PDFs of the QoI are our criteria to assess whether the basis adaptation process is converged. In the classical basis adaptation approach described above, only the first r random variables in are considered in the reduced space. The accuracy of the adapted expansion increases by increasing r and will recover the full dimension expansion with r =d. 4.2 Acceleration of the convergence of the basis adaptation method The convergence rate of the PCE to the dimension of the adaptation depends on isometry A. In the meantime, the efficiency of the basis adaptation method is directly related to the converged dimension of the adaptation; hence, a low converged dimension is preferable. This section proposes two methods that can accelerate the convergence rate of the classical Gaussian adaptation. The first one exploits the probabilistic information embedded in the pilot first order PCE in the original space obtained in the first step of basis adaptation. The second algorithm utilizes higher dimensional adaptations to update the rotation matrix A rather than leaving it constant as the dimension increases. 74 4.2.1 Correction of the mean and Gaussian coefficients To construct the rotation matrix in the classical Gaussian adaptation, a pilot PCE up to order 1 in the original space is first performed Q J d;1 () = X 2J d;1 Q () =Q e0 +Q e1 1 + +Q e d d ; (4.20) wheree 0 is ad-dimensional multi-index of zeros everywhere, and thus,Q e0 is the 0 th order coefficient, which represents the mean value of bothQ andQ J d;1 , and we refer to it as the mean coefficient of the PCE. Define Q e coeff :=fQ ei g d i=1 as the set of first order (or Gaussian) coefficients. From Eq (4.5), these PCE coefficients can be obtained as Q e0 = X ( (q) ;wq )2Q d Q (q) w q ; Q e1 = X ( (q) ;wq )2Q d Q (q) (q) 1 w q ; :::; Q e d = X ( (q) ;wq )2Q d Q (q) (q) d w q : (4.21) Then, the first row ofA, by the construction of Eq (4.12) and the subsequent Gram-Schmidt procedure, is given byQ e coeff =kQ e coeff k 2 . In the adapted reduced space, the PCE up to order 1 for r-d adaptation can be written as Q Ar ( r ) = X 2Jr;1 Q Ar ( r ) =Q Ar e0 +Q Ar e1 1 + +Q Ar er r ; (4.22) wheree i is the r-dimensional multi-index with 1 at the ith location and zeros elsewhere; Q Ar e0 is the mean value of Q Ar ; andfQ Ar ei g r i=1 is the set of Gaussian coefficients. From Eqs (4.14) and (4.15), the PCE coefficients in Eq (4.22) are given by Q Ar e0 = X ( (q) r ;wq )2Qr Q A T (q) r 0 T ! w q ; Q Ar e1 = X ( (q) r ;wq )2Qr Q A T (q) r 0 T ! (q) 1 w q ; :::; Q Ar er = X ( (q) r ;wq )2Qr Q A T (q) r 0 T ! (q) r w q : (4.23) 75 Note that Q in (4.23) and (4.21) is an operator evaluated by the physical modelM. As we mentioned in Eqs (4.14) and (4.15), the argument of Q in (4.23), = A T = A T [ r ; :r ] T , is approximated by A T [ r ; 0] T ; that is, the probabilistic contents embedded in the random vector :r is ignored. In the classical basis adaptation approach, it is claimed that the first row of the rotation matrix A captures the complete Gaussian components, which holds when the error induced on the Gaussian components by the approximation :r = 0 is negligible. To measure the error in the mean and Gaussian coefficients in Eq (4.23) compared to Eq (4.21), we can project Eq (4.22) to the basisf g 2J d;1 in the original space (in Eq (4.20)) as X 2Jr;1 Q Ar ( r ) = X 2J d;1 ~ Q (); (4.24) wheref ~ Q ;2J d;1 g are the projected PCE coefficients, which, by noting that r =A r , can be computed by substituting (4.20) and the basis of (4.22) into (4.24) as ~ Q e0 =Q Ar e0 ; ~ Q e1 = r X i=1 Q Ar ei A i1 ; ; ~ Q e d = r X i=1 Q Ar ei A id ; (4.25) whereA ij is the entry of the rotation matrixA ofi-th row andj-th column. We recall thatfA 11 ; ;A 1d g = Q e coeff =kQ e coeff k 2 , yielding that, for one-dimensional adaptation (r = 1), as expected, ~ Q e0 =Q Ar e0 ; ~ Q e1 = Q Ar e1 Q e1 kQ e coeff k 2 ; ::: ; ~ Q e d = Q Ar e1 Q e d kQ e coeff k 2 : (4.26) Define ~ Q e coeff = f ~ Q ei ;i = 1;:::;dg, the arithmetic errors of mean and Gaussian coefficients for the 1d adaptation are 0 D =Q e0 ~ Q e0 =Q e0 Q e 1 0 and 1 D =Q e coeff ~ Q e coeff = 1 Q Ar e1 kQ e coeff k 2 Q e coeff ; (4.27) where p D denotes the arithmetic error of PCE coefficients associated with polynomials of order p in the original space. From Eq (4.27) we can see that 1 D is proportional to the exact Gaussian coefficientsQ e coeff , 76 implying that larger arithmetic errors will appear in terms with larger Gaussian coefficients. This is alarming since a larger value of the Gaussian coefficient usually signifies greater sensitivity in the associated PCE term. Furthermore, larger errors in the sensitive PC terms will induce relatively more significant errors in the expansion. Eq (4.26) is special case of Eq (4.25) with r = 1. With the increase of r, more terms are considered in the projected Gaussian coefficients of the adaptation, which are approaching the exact ones gradually and will recover the coefficients with r =d. The Gaussian coefficients in the classical basis adaptation approach converge to the exact ones with increasing dimensions from the above analysis. However, since the exact mean and Gaussian coefficients of the PCE in the original space are the building blocks of the rotation matrixA and are obtained in the first step of the basis adaptation approach, we can use the exact mean and Gaussian coefficients to correct the associated coefficients in the rotated and reduced spaces. To do so, we can project the first order PCE in the original space to the basis of the reduced spacef Ar ;2J r;1 g (which is the reverse of Eq (4.24)) X 2J d;1 Q () = X 2Jr;1 ~ Q Ar ( r ); (4.28) wheref ~ Q Ar ;2J r;1 g are the projected PCE coefficients in the rotated and reduced space, and, can be obtained explicitly by substituting (4.22) and the basis of (4.20) into (4.28) as ~ Q Ar e0 =Q e0 ; ~ Q Ar e1 = d X i=1 Q ei A 1i ; ::: ; ~ Q Ar er = d X i=1 Q ei A ri : (4.29) Thenf ~ Q Ar e0 ; ~ Q Ar e1 ; ; ~ Q Ar er g are the exact mean and Gaussian coefficients in the rotated and reduced space and can be used to correct the values offQ Ar e0 ; Q Ar e1 ; ; Q Ar er g in the r-d adaptation (in (4.22)). Correctingthe meanandGaussiancoefficientsacceleratesthe convergenceof thebasisadaptationmethod since the zero, and first order terms are now guaranteed to be accurate. We now only need to consider the higher order terms in the convergence analysis, and the necessity to obtain converged PCE coefficients of zero and first order term is removed. The accelerating effects are more pronounced for the problems that can be well approximated by the zero and first order polynomials. 77 4.2.2 Updating the rotation matrix by higher dimensional adaptation With each increment of the reduced dimensionr in basis adaptation (4.13), the resulting PCE is more closely aligned with the full-dimensional PCE, and new information that becomes available at that stage should be pertinent for constructing a more adapted isometryA. The r-dimensional adaptation of PCE contains the sum of monomials in r random variables having a total degree between zero and p. This section presents an algorithm to update the initial rotation matrixA by incorporating additional information gained as the dimension, r, of the adaptation is sequentially increased. Denote the new rotation matrix updated from A asB. The idea is to buildB such that it provides a one-dimensional adaptation of order up to p that is optimally close to an r-dimensional (r 2) adaptation (of order up to p) by rotation matrix A. Here we are seeking a one-dimensional adaptation from B, and thus only its first row is of interest. We recall that the first row ofA was constructed to capture the relative sensitivities expressed by the first order coefficients and is the leading direction in the rotated space. Suppose that the r-d adaptation by rotation matrixA is known, and we seek the best one-dimensional basis adaptation associated withB, that is Q B1 ( 1 ) = X 2J1;p Q B1 ( 1 ) = X 2J1;p Q B1 B1 (); (4.30) where =B. We avoid additional physical model evaluations by leveraging the available r-d adaptation already computed using rotation A. The PCE coefficientsfQ B1 g 2J1;p in Eq (4.30) can be approximated by projecting the r-d adaptation Q Ar = P 2Jr;p Q Ar ( r ) to the one-dimensional basisf B1 g 2J1;p as Q B1 = X 2Jr;p Q Ar h Ar ; B1 i; 2J 1;p : (4.31) We note that the expansion in Eq (4.30) is not a first order expansion, and therefore these coefficients do not represent the first row of matrixB. Now, for any given rotationB, the one-dimensional adaptation up 78 to a specific order is built. To compare the PCE coefficientsfQ B1 g 2J1;p andfQ Ar g 2Jr;p , Eq (4.30) is projected to the space rotated byA (spanned byf Ar g 2Jr;p ) as X 2J1;p Q B1 ( 1 ) = X 2Jr;p ~ Q B1 ( r ); (4.32) then, by substituting Eq 4.31 into the above equation,f ~ Q B1 g 2Jr;p can be obtained as ~ Q B1 = X 2J1;p Q B1 h B1 ; Ar i = X 2J 1;p 0 @ X 2Jr;p ^ Q Ar h Ar ; B1 i 1 A h B1 ; Ar i; 2J r;p : (4.33) The coefficientsf ~ Q B1 g are actually the projected coefficients of the one-dimensional adaptation byB to the r-dimensional basisf Ar ;2J r;p g. Then, to capture the probabilistic information of the r-dimensional adaptation by rotationA optimally, rotationB is constructed such that the relative 2-norm error between PCE coefficients off ~ Q B1 ;2J r;p g andfQ Ar ;2J r;p g is minimized, which can be expressed as, B 1 = arg min B12R d kQ Ar coeff ~ Q B1 coeff (B 1 )k 2 kQ Ar coeff k 2 ; (4.34) where Q Ar coeff :=fQ Ar ;2J r;p g, and ~ Q B1 coeff (B 1 ) :=f ~ Q B1 ;2J r;p g is a function of B (as can be seen from Eq (4.33)). The Basin-hopping algorithm [134], which is particularly useful for global optimization of high-dimensional problems, is adopted in this study to obtain the global minimum of the objective function. The other rows (from 2 to d) ofB copy the same construction ofA, and thus, the idea of using the relative sensitivity of the original variables to construct these rows is kept. Although we only update the first row of the rotation, the change also updates the other rows by following the Gram-Schmidt procedure. We say thatB is built from the r-dimensional adaptation by rotationA when Eq (4.34) is used to constructB. In the optimization problem of Eq (4.34),Q Ar coeff is known, and ~ Q B1 coeff (B) is obtained using Eqs (4.30) to (4.33), which involves only projections between two PCEs. Since no additional evaluations of the physical model are needed, the cost of the optimization is negligible. Suppose the difference between the one- dimensional adaptation by B and the r-dimensional adaptation by A is small, then PCE Q B1 in (4.30) is approximately in the polynomial span of the first r rows ofA. 79 With the increase of dimension r in the adaptation byA, the first row ofB will gain more probabilistic information. Thus, the 1d adaptation based onB will have higher accuracy, enabling us to have accelerated convergence for the adaptations byB. We could also increase the dimension of the adapted representation obtained fromB, as needed for convergence. DefiningB (1) asA and renamingB asB (2) , this procedure can be iterated, using B (s) as the prior adaptation matrix and computing the best leading dimension of a new adaptation matrix, B (s+1) , s 1. The dimension of the reduced representation using B (s+1) will then be decided based on a low-dimensional convergence study. It is reasonable to expect that the above iterative process may exhibit saturation as s increases. We found that increasing the dimension of the reduced representation obtained from successive s values mitigates this problem. Specifically, when using an adaptation matrixB (s) , we typically compare the accuracy of an s and an (s + 1)-dimensional reduced models (i.e., we respectively retain s and (s + 1) rows inB (s) ). In more detail, we first perform 1d and 2d adaptations using A, and the procedure is stopped if the results are within tolerance. Otherwise, the updated rotation matrixB (2) is computed from the previous 2d adaptation, followed by one new 2d and one 3d adaptation. Again, the procedure is stopped if the results are within tolerance. Otherwise, another updated rotation matrix B (3) is constructed using the previously computed 3d adaptation (based on B (2) ). One 3d and one 4d adaptation are then performed using B (3) . We continue this procedure until the results are within tolerance. Taking advantage of the known exact mean and Gaussian coefficients, the second method is combined with the mean and Gaussian coefficients correction as described in section 4.2.1. That is, we continually update the mean and Gaussian coefficients of the adapted PCE. The flowchart in Fig 4.1 shows the complete procedure to perform adaptations using this sequential adaptation procedure, referred to as sequentially optimized adaptations (SOA). The convergence between two successive adaptations can be either based on PCE coefficients or PDFs of the QoI. In the former case, we compare the PCE coefficients of two successive basis adaptations (the coefficients of the lower dimensional adaptation are projected to higher dimensional space). At the same time, the convergence criterion of PDFs requires one to first construct the PDFs by, for example, the kernel density estimation (KDE) [135]. 80 start r = 2 B (1) =A (r 1)-d adaptation usingB (r1) r-d adaptation usingB (r1) Results in tol- erance? stop Construct new B (r) from r- d adaptation usingB (r1) Yes No r =r + 1 Figure 4.1: Workflow of the sequentially optimized adaptation method. 4.3 Applications In this section, two examples show the applicability of the methods mentioned in this chapter. The Uncer- tainty Quantification Toolkit (UQTk) is used for the stochastic analysis in the present work [136]. 4.3.1 Borehole function Consider the following borehole water flow function f(x) = 2T u (H u H l ) ln r rw 1 + 2LTu ln( r rw )r 2 w Kw + Tu T l ; (4.35) which models the water flow through a borehole drilled from the ground surface through two aquifers [137]. In Eq (4.35), r w = 0:1 (m) is the deterministic radius of the borehole, and the seven other parameters are random and modeled by distributions shown in Table 4.1. Uncertain Parameter Distribution Units Radius of influence (r) Lognormal( = 7:71; = 1:0056) m Transmissivity of upper aquifer (T u ) U[63070; 115600] m 2 /yr Potentiometric head of upper aquifer (H u ) U[990; 1110] m Transmissivity of lower aquifer (T l ) U[63:1; 116] m 2 /yr Potentiometric head of lower aquifer (H l ) U[700; 820] m Length of borehole (L) U[1120; 1680] m Hydraulic conductivity of borehole (K w ) U[9855; 12045] m/yr Table 4.1: Uncertain parameters and distributions for the borehole function. 81 The uncertainties of the input parameters are propagated to the QoI f(x), the water flow rate of the borehole (m 3 /yr). The dimension of the forward UQ problem is thus d = 7. The classical basis adaptation methodisfirstimplementedforcomparisonpurposes. TobuildtherotationmatrixA, afull-dimensionalpilot PCE, Eq (4.20), in the original space is carried out, and the mean and Gaussian coefficients are obtained by Eq (4.21), requiring level 1 Smolyak quadrature (15 physical model evaluations). Next, the Gaussian coefficients are used to build the rotation matrix A by Algorithm 4.1. Then, adaptations of Eq (4.13) are implemented, in which the adapted PCEs are up to order p = 3, and a level 3 Smolyak quadrature is used. The probability density functions (PDFs) of f(x) obtained by 1d to 5d adaptations are presented in Figure 4.2(a), and the comparison of PCE coefficients of 4d and 5d adaptations is depicted in Figure 4.2(b), in which the PCE coefficients of the 4d adaptation are projected to the 5-dimensional space. Note that the 20 40 60 80 100 120 140 Water flow rate (m 3 /yr) 0.000 0.005 0.010 0.015 0.020 0.025 PDF Full dimension PCE 1d classical adapt. 2d classical adapt. 3d classical adapt. 4d classical adapt. 5d classical adapt. (a) 10 20 30 40 50 PCE term 0 2 4 6 8 10 12 14 Coefficient of PCE term PC coeffs. of 4d adapt. PC coeffs. of 5d adapt. (b) Figure 4.2: (a) PDFs of classical basis adaptation of dimensions from 1 to 5; and (b) PCE coefficients comparison of 4d and 5d adaptations in the space spanned byf A5 ;2J 5;p g PDF of f(x) obtained by the full-dimensional PCE of order up to 3 in the original space is also plotted in Figure 4.2(a) to serve as the reference, and in Figure 4.2(b), the mean coefficients (that is the 0 order coefficients) are not presented since they are much greater than the other PCE coefficients and their presence will make the results hard to read. The mean coefficients are not included in any figures comparing the PCE coefficients of different expansions in this study. It can be seen from Figure 4.2(a), the PDF of 4d adaptation converges to the given reference PDF. However, the 5d adaptation is required to confirm the convergence since the reference is usually unavailable in practice. The relative 2-norm error of the PCE coefficients in 82 Figure 4.2(b) is 0.2%, which is small enough to suggest the convergence of the 4d adaptation. Using level 3 Smolyak quadrature, the number of physical model evaluations of 1d to 5d adaptations are 3, 73, 159, 289 and 471, respectively, which gives a total number of model evaluations as 3 + 73 + 159 + 289 + 471 = 995. Whileinmanyapplicationsthattheclassicalbasisadaptationmethodcanconvergeveryfast(inthefirst2 or3dimensions)[65,66], theboreholewaterflowproblemconvergesrelativelyslow. Themethodthatcorrects the mean and Gaussian coefficients as described in section 4.2.1 is used to accelerate the convergence of the classical basis adaptation method. The mean and Gaussian coefficients of any r-d adaptation are updated by the exact ones (obtained from the pilot PCE), see Eqs 4.28 and (4.29). The PDFs of f(x) obtained by the 1d to 4d adaptations with updated coefficients (UCs) are presented in Figure 4.3(a), and the comparison of PCE coefficients of 3d and 4d adaptations in the space spanned byf A4 g 2J4;p is depicted in Figure 4.3(b). In engineering, the tails of the PDFs are of high interest since they are more related to extreme events 20 40 60 80 100 120 140 Water flow rate (m 3 /yr) 0.000 0.005 0.010 0.015 0.020 0.025 PDF Full dimension PCE 1d adapt. with UCs 2d adapt. with UCs 3d adapt. with UCs 4d adapt. with UCs (a) 5 10 15 20 25 30 PCE term 0 2 4 6 8 10 12 14 Coefficient of PCE term PC coeffs. of 3d adapt. with UCs PC coeffs. of 4d adapt. with UCs (b) Figure 4.3: (a) PDFs of basis adaptation of dimensions up to 4 with updated coefficients (UCs); and (b) PCE coefficients comparison of 3d and 4d adaptations with UCs in the space spanned byf A4 ;2J 4;p g [38, 138]. Under such premise, Figure 4.3(a) suggests that the 3d adaptation with UCs converges to the reference PDF, while in practice, additional 4d adaptation with UCs is required to confirm the convergence. The PCE coefficients comparison in Figure 4.3(b) also indicates the convergence, where the relative 2-norm error of the PCE coefficients is 7.0%. Same as the classical adaptation, level 3 Smolyak quadrature is used, and the total number of physical model evaluations of 1d to 4d adaptations is 3 + 73 + 159 + 289 = 524, which implies that the computational cost by this method is 53% of the classical basis adaptation method. 83 The method that uses higher dimensional adaptations to update the rotation matrix as described in section 4.2.2 is applied to reduce the computational cost further. Following the procedure presented in flowchart 4.1, the 1d and 2d adaptations are first built based on the rotation matrix A with convergence check followed. The results are shown in Figure 4.4. This method is combined with the mean and Gaussian coefficients correction approach; thus, Figures 4.4(a) and 4.3(a) have the same PDFs for 1d and 2d adap- tations. Figure 4.4 shows that the 2d basis adaptation is not converged, and further processes are needed. The relative 2-norm error in Figure 4.4(b) is 10.5%. 20 40 60 80 100 120 140 Water flow rate (m 3 /yr) 0.000 0.005 0.010 0.015 0.020 0.025 PDF Full dimension PCE 1d adapt by B (1) 2d adapt by B (1) (a) 1 2 3 4 5 6 7 8 9 PCE term 2 0 2 4 6 8 10 12 14 Coefficient of PCE term PC coeffs. of 1d adapt. by B (1) PC coeffs. of 2d adapt. by B (1) (b) Figure 4.4: (a) PDFs of 1d and 2d adaptations byB (1) ; and (b) PCE coefficients comparison of 1d and 2d adaptations byB (1) in the space spanned byf A2 ;2J 2;p g Then the 2d adaptation information is utilized to obtain a new rotation matrix B (2) such that the L 2 difference is minimized between the 1d adaptation using this new matrix and the 2d adaptation achieved with the previous matrixB (1) . The comparison of absolute values of the first rows ofB (1) andB (2) is shown in Figure 4.5. The first row of these rotation matrices represents the identified leading direction. We can see in the figure thatB (2) actually has a much different leading direction, and this direction will also influence further directions constructed via a Gram-Schmidt procedure. However, the two most sensitive variables (from their PCE coefficients), variable 6 and 3, retain their dominance in the updated rotation. Using the new rotation matrixB (2) , one 2d and one 3d adaptation are carried out, and the PDFs off(x) and the comparison of PCE coefficients are shown in Figure 4.6. Figure 4.6(a) suggests the convergence of the 2d adaptation using B (2) , which is confirmed by its consistency with the 3d adaptation. The relative 84 1 2 3 4 5 6 7 Variable 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Normalized PCE coefficient First row of B (1) First row of B (2) Figure 4.5: Comparison of absolute values of the first rows of rotationsB (1) andB (2) . 20 40 60 80 100 120 140 Water flow rate (m 3 /yr) 0.000 0.005 0.010 0.015 0.020 0.025 PDF Full dimension PCE 2d adapt by B (2) 3d adapt by B (2) (a) 2 4 6 8 10 12 14 16 18 PCE term 2 0 2 4 6 8 10 12 14 Coefficient of PCE term PC coeffs. of 2d adapt. by B (2) PC coeffs. of 3d adapt. by B (2) (b) Figure 4.6: (a) PDFs of 2d and 3d adaptations byB (1) ; and (b) PCE coefficients comparison of 2d and 3d adaptations byB (2) in the space spanned byf A3 ;2J 3;p g 2-norm error of the PCE coefficients between these two adaptations in Figure 4.6(b) is 9%. Comparing with Figure 4.4(b), we note that even a tiny change in the PCE coefficients impacts the PDFs. Different basis terms in PCE have different contributions to the QoI, and the same error on a basis term with greater contribution will have a greater impact on the resulting expansion. As a result, the convergence of PCE coefficients usually implies the convergence of PDFs, while the contrary is not necessarily true. Although the relative 2-norm error in Figure 4.6(b) is only slightly smaller than in Figure 4.4(b), the improvement in some particular basis terms could significantly increase the accuracy of the expansion. For example, we see better convergence of the quadratic term of the first adapted variable in Figure 4.6(b) (the 4th term) than in Figure 4.4(b) (the 3rd term), which could be part of the reason for a converged PDF in Figure 4.6(a), since, 85 the first adapted variable by construction is the most important one, and we have reason to believe that its quadratic term will also be important. Thus, the PCE coefficient comparison can help us to identify the contributions of the PCE terms. In this method, one 1d, two 2d, and one 3d basis adaptation are performed to achieve a convergence, which requires 3 + 2 73 + 159 = 308 physical model evaluations. The total cost of the method is only 30% of the classical basis adaptation method. MC sampling is a popular candidate for high-dimensional problems, and it is applied to the borehole function with different numbers of MC samples; the full-dimensional PCE of order 3 will be served as a reference. The resulting PDFs are presented in Figure 4.7(a), and a zoom-in view on the right tail of the PDF is shown in Figure 4.7(b). From Figure 4.7(a), we see that the overall shapes of the PDFs obtained from the MC method with 50, 100, and 500 samples have significant discrepancies from the reference PDF. If we focus on the right half from the peak, the PDFs of the MC method with 1000 and 2000 samples also have noticeable differences from the reference. For the right tail region, see in Figure 4.7(b), with the increase of MC samples, the tails of the PDFs gradually approach the tail of the reference and exhibit a good match with 3000 samples. This should be compared to 308 physical model evaluations to achieve a similar accuracy using the present SOA approach (see Figure 4.6). 20 40 60 80 100 120 140 Water flow rate (m 3 /yr) 0.000 0.005 0.010 0.015 0.020 0.025 PDF Full dimension PCE MC 50 MC 100 MC 500 MC 1000 MC 2000 MC 3000 (a) 100 105 110 115 120 125 130 Water flow rate (m 3 /yr) 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 PDF Full dimension PCE MC 50 MC 100 MC 500 MC 1000 MC 2000 MC 3000 (b) Figure 4.7: (a) PDFs obtained from Monte Carlo method with different samples and (b) a zoomed in view in the right tail region; the full dimensional PCE of order 3 is served as reference model. 86 4.3.2 Space structure model In the second example, a major component of a space structure is considered. The computational model is shown in Figure 4.8. The upper and lower parts of the structure are connected by the three-point mounting pedestals (modeled by linear constraints). The upper part is open, and there is a large mass (approximately two orders of magnitude heavier than the rest of the model) at the center connected to the sidewalls by several rigid links. The essential parts of the structure are contained in the lower part, which is composed of a cylindrical outer shell (see in Figure 4.8(a)) and an inner shock absorption block that contains hollow spaces to hold the essential components. Figure 4.8(b) presents the structure without the outer shell and inner shock absorption block such that the essential components are exposed. The model is of free boundary condition and is subjected to an impulse-type shock loading on the center mass of the upper part in Z direction, with its time history shown in Figure 4.9(a). The QoI is the maximum acceleration in X direction along with the time history of a critical point located at one of the essential components. Figure 4.9(b) depicts the time history of acceleration in X direction of the critical point when the model is deterministic, of which the maximum X acceleration is about 1:9 10 5 (in/s 2 ). (a) (b) Figure 4.8: (a) Finite element model of the space structure; (b) Finite element model of the space structure without the outer shell and inner absorption block. The dynamic effects of the impulse loading are first transmitted to 3 major components of the upper part via the rigid links, then to the three-point pedestals, through which the effects of the loading are 87 0.00 0.01 0.02 0.03 0.04 0.05 0.06 Time (s) 2 1 0 1 2 Force (lbf) 1e5 (a) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Time (s) 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 Acceleration in X (in/s 2 ) 1e4 (b) Figure 4.9: (a) Time history of the impulse force and (b) time history of acceleration in X direction of the critical point transferred to the lower part. Suppose that the modulus of elasticity of 3 major components of the upper parts (nominal value of 1:6 10 7 psi), the stiffness of 3 linear joints (nominal of value 5:7 10 7 lbs/in), and modulus of elasticity of 18 components of the lower part (14 of which with nominal values of 1:0 10 7 psi, 2 of which with nominal values of 9:76 10 7 psi, 1 of which with a nominal value of 2:93 10 7 psi and 1 of which with a nominal value of 8829.0 psi) are log-normal random variables (RVs) with mean values of their nominal values and coefficients of variance 10%, 20%, and 20%, for these three groups, respectively. The 18 components of the lower parts are chosen such that they are connected and will affect the response at the critical point, and the choices of coefficients of variance are based on the variabilities of the materials. Thus, the dimension of the forward UQ problem is d = 24. For the PCE, it is assumed that order p = 3 is enough to capture the probabilistic information, and level 3 Smolyak quadrature is used for numerical integration. Under such circumstances, the full-dimension PCE that requires 24,449 evaluations of the physical model is too expensive to handle; hence, the classical and accelerated basis adaptation approaches are implemented. In this application, the physical model is evaluated by LS-DYNA software [129], which is coupled with UQTk [139] to perform the PC expansion and basis adaptation of the QoI. A full-dimensional pilot PCE in the original space of Eq (4.20)is first carried out, and the mean and Gaussian coefficients are obtained by Eq (4.21), in which a level 1 Smolyak quadrature (49 evaluations of the physical model) is used. Then, the rotation matrixA is built, and the classical basis adaptation method is 88 applied to obtain adaptations of the QoI. The PDFs of the QoI obtained by 1d to 9d adaptations of order up top = 3 are presented in Figure 4.10, in which the level 3 Smolyak quadrature is used. In addition, the PDFs obtained by the full-dimensional PCE of order up to 1 and 2 in the original space are also plotted in Figure 4.10, which serve as potential reference PDFs. The full-dimensional PCE of order up to 2 is obtained by 17000 17500 18000 18500 19000 19500 20000 20500 Maximum X acceleration 0.000 0.002 0.004 0.006 0.008 0.010 PDF 24d PCE of order 1 24d PCE of order 2 1d adapt. of order 3 2d adapt. of order 3 3d adapt. of order 3 4d adapt. of order 3 5d adapt. of order 3 6d adapt. of order 3 7d adapt. of order 3 8d adapt. of order 3 9d adapt. of order 3 Figure 4.10: PDFs of the maximum X acceleration of the critical point obtained by full dimensional PCE and classical basis adaptations. level 2 Smolyak quadrature that requires 1,297 evaluations of the physical model, which is the best reference model we can get, limited by the computing resources of 100 compute nodes (each has 63.4G of memory and 16 CPUs). As shown in Figure 4.10, regarding the full-dimensional PCEs, although the convergence of the order 2 PCE is not confirmed (which requires PCE of order 3), the order 2 and order 1 PCEs are in good consistency in both tails. Figure 4.10 also shows that different dimensional adaptations converge slowly to the dimension and are not converged even when the dimension goes to 9. Moreover, the PDFs of different dimensional adaptations are far from the full-dimensional PCEs in the original space. It might be interesting to check if the first order information is converged for the 9d adaptation. Figure 4.11(a) shows the PDF comparison of first order expansions between full-dimensional PCE in the original spaceand9dadaptationintherotatedspace. Incontrast,Figure4.11(b)presentsthecorrespondingGaussian coefficients comparison between these two expansions, in which the coefficients of the 9d adaptation are projected to the original space for comparison. Figure 4.11(a) indicates that the first order expansion of the 9d adaptation is not converged to the exact one, and this can also be seen from Figure 4.11(b), where the Gaussian coefficients of the 9d adaptation are very different from the exact one. To obtain the 9d 89 basis adaptation, the number of quadrature points required for level 3 Smolyak quadrature is 1879, which is already very expensive; the total number of quadrature points for adaptations from 1d to 9d is 4996. Hence, the classical basis adaptation approach is insufficient to achieve acceptable reduction due to the slow convergence. 16500 17000 17500 18000 18500 19000 19500 20000 20500 21000 Maximum X acceleration 0.0000 0.0005 0.0010 0.0015 0.0020 PDF 24d PCE of order 1 9d adapt. of order 1 (a) 5 10 15 20 PCE term 200 100 0 100 200 300 400 Coefficient of PCE term 24d PCE of order 1 9d adapt. of order 1 (b) Figure 4.11: (a) PDF comparison of first order expansions between full dimensional PCE and 9d adapta- tion; (b) and Gaussian coefficient comparison between full dimension PCE and 9d adaptation in the full dimensional space. The method described in section 4.2.1 that updates the mean and Gaussian coefficients by the exact ones is first applied to accelerate the convergence in the classical basis adaptation approach. The resulting PDFs of the QoI are shown in Figure 4.12(a), and Figure 4.12(b) presents the PCE coefficients comparison between 2d and 3d adaptations with UCs. From the PDF comparison, we can see that the adaptation of the QoI with UCs converges at dimension 2, which is confirmed by the 3d adaptation with UCs. Both 2d and 3d adaptations with UCs are close to the full-dimensional PCE of order up to 2, though there are some imperfections in both tails of the PDFs. The relative 2-norm error of the coefficients in Figure 4.12(b) is 3.8%, which also indicates the consistency of 2d and 3d adaptations with UCs. Instead of a slow convergence rate up to dimension 9 that fails to obtain a converged result in the classical approach, the basis adaptation with UCs is converged with dimension 3; moreover, the 3d adaptation is close enough to the best available reference model. By level 3 Smolyak quadrature, the total number of required physical model evaluations to obtain 1d, 2d, and 3d adaptations with UCs is 3 + 73 + 159 = 235. 90 16000 17000 18000 19000 20000 21000 Maximum X acceleration 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 PDF 24d PCE of order 1 24d PCE of order 2 1d adapt. of order 3 with UCs 2d adapt. of order 3 with UCs 3d adapt. of order 3 with UCs (a) 2 4 6 8 10 12 14 16 18 PCE term 0 200 400 600 Coefficient of PCE term PC coeffs. of 2d adapt. with UCs PC coeffs. of 3d adapt. with UCs (b) Figure 4.12: (a) PDFs of the maximum X acceleration of the critical point obtained by basis adaptation of dimensions up to 3 with updated coefficients (UCs); and (b) PCE coefficients comparison of 2d and 3d adaptations with UCs in the space spanned byf A3 ;2J 3;p g. The SOA method is then applied to this space structure model. Again, this method is coupled with the mean and Gaussian coefficients correction method. The 1d and 2d basis adaptation with UCs are first built using the rotation matrixA, and note that the initial rotation matrix (in flowchart Figure 4.1) B (1) =A. The PDFs by the 1d and 2d basis adaptation with UCs are shown in Figure 4.13(a), and the corresponding PCE coefficients comparison in space spanned byf A2 g 2J2;p is presented in Figure 4.13(b). It can be seen from Figure 4.13(a) that the 1d and 2d adaptations are not consistent with each other, and the comparison with the candidate reference models indicates that the 2d adaptation does not converge to the reference. The relative 2-norm error between the two sets of PCE coefficients in Figure 4.13(b) is 23.0%, indicating that the 1d and 2d adaptations are far from each other. Then the 2d adaptation with UCs is utilized to obtain a new rotation matrixB (2) by global optimization. The new rotation matrix aims to concentrate the probabilistic information of the 2d adaptation ofB (1) to the first dimension of a new adaptation. To see how different the updated rotation is from the previous rotation, the absolute values of the first rows of B (1) and B (2) are compared in Figure 4.14. We can see in this figure that the difference between these two rotations is not as significant as we saw on the borehole function example in Figure 4.5; however, there are some noticeable differences in some variables that are considered sensitive to the QoI (based on their coefficients), for example, variables 7, 14 and 18. By the new rotation matrixB (2) , the 2d and 3d adaptations with UCs are implemented, and the result PDFs of the QoI are shown 91 16000 17000 18000 19000 20000 21000 Maximum X acceleration 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 PDF 24d PCE of order 1 24d PCE of order 2 1d adapt. of order 3 by B (1) 2d adapt. of order 3 by B (1) (a) 1 2 3 4 5 6 7 8 9 PCE term 0 200 400 600 Coefficient of PCE term PC coeffs. of 1d adapt. by B (1) PC coeffs. of 2d adapt. by B (1) (b) Figure 4.13: (a) PDFs of the maximum X acceleration of the critical point obtained by 1d and 2d adaptations by B (1) ; and (b) PCE coefficients comparison of 1d and 2d adaptations by B (1) in the space spanned by f A2 ;2J 2;p g 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Variable 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Normalized PCE coefficient First row of B (1) First row of B (2) Figure 4.14: Comparison of absolute values of the first rows of rotationsB (1) andB (2) for the space model. in Figure 4.15(a), the corresponding PCE coefficients comparison in the space spanned byf A2 g 2J2;p are presented in Figure 4.15(b). One can see from Figure 4.15(a) that the 2d and 3d adaptations by B (2) are consistent with each other, and they also converge to the potential reference models. Though there is some discrepancy in the peak area of the PDFs between the 2d adaptation and the full-dimensional expansion of order up to 2, the convergence in the tail, especially the right tail, is nearly perfect. As the tails of the PDFs are related to the extreme events, which are more important in engineering, the 2d and 3d adaptations by B (3) can well approximate the maximum X acceleration at the critical point. The comparison of PC coefficient of 2d and 3d adaptations byB (3) in Figure 4.15 also indicates that they are consistent with each 92 16000 17000 18000 19000 20000 21000 Maximum X acceleration 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 PDF 24d PCE of order 1 24d PCE of order 2 2d adapt. of order 3 by B (2) 3d adapt. of order 3 by B (2) (a) 2 4 6 8 10 12 14 16 18 PCE term 0 200 400 600 Coefficient of PCE term PC coeffs. of 2d adapt. by B (2) PC coeffs. of 3d adapt. by B (2) (b) Figure 4.15: (a) PDFs of the maximum X acceleration of the critical point obtained by 2d and 3d adaptations by B (2) ; and (b) PCE coefficients comparison of 2d and 3d adaptations by B (2) in the space spanned by f A3 ;2J 3;p g other; the 2-norm error of the coefficients is 5.3%. In the whole process to obtain the convergence result, one 1d, two 2d, and one 3d basis adaptations are calculated, which, by using level 3 Smolyak quadrature, requires 3 + 2 73 + 159 = 308 times physical model evaluations. While the computational cost is similar to that of the mean and Gaussian coefficients correction method, the converged expansion by the method that updates the rotation matrix can better approximate the tails of the PDF. To conclude this chapter, two methods that can accelerate the convergence rate of the classical basis adaptation method in the framework of PCE are proposed. These two methods can be combined to achieve accelerated convergence. Particular application scenarios are stochastic analyses of complex systems where the classical basis adaptation approach has slow convergence. 93 Chapter 5 Stochastic reduction for the high-dimensional quantity of interest problems with parametric and nonparametric uncertainties The previous chapter proposed dimension reduction techniques adapted to scalar QoI or QoI with limited dimensions. This chapter will introduce the stochastic reduction technique adapted to the high-dimensional QoI problems with parametric and nonparametric uncertainties. The chapter is organized as follows: 1. The KLE is recalled, which will be used to construct an expansion of the high-dimensional QoI. 2. The integration of the KLE and adapted PCE will be introduced, with KLE the tool to replace the high-dimensionalQoIwithalimitednumberofKLtermsandtheadaptedPCEtoreducethedimension required to construct probabilistic models for individual KL terms. The probabilistic models of the KL terms can then be substituted back to the KLE of the QoI to construct a surrogate of the QoI. 3. The fundamentals of the nonparametric stochastic analysis are presented and incorporated into the parametric approach to consider the uncertainties induced by modeling errors. 4. The method introduced in this chapter is applied to the detailed FA and the FLSNFC. 5.1 Karhunen-Loève expansion Let X(t;#) be a random process defined on the probability space ( ;F;P), which is a collection of random variables, with #2 representing the random events and t2 T. The set T is the index or parameter 94 set of the random process that can be an interval. Also, let X (t) be the mean of X(t;#) and K(s;t) = Cov (X(s;#);X(t;#)) be the covariance function. Then the KLE of process X(t;#) can be written as X(t;#) = X (t) + 1 X i=1 p i Z i (#)e i (t); (5.1) where i and e i are eigenvalue and eigenvector pairs of the Covariance matrix K(s;t), that is i e i (s) = Z T K(s;t)e i (t)dt; s2T: (5.2) From Eq (5.1), the KL parameters (or KL coefficients)fZ i (#)g can be obtained as Z i (#) = 1 p i Z (X(t;#) X (t))e i (t)dt; 8i (5.3) which are uncorrelated random variables with zero mean and unit variance. The decay of eigenvalues of the covariance matrix depends on the correlation length of the random process. Based on the decay, one can truncate the KLE to m terms such that the eigenvalues of the KL terms that are greater than m have small values. Then the process can be approximated as X(t;#) X (t) + m X i=1 p i Z i (#)e i (t); with m 1: (5.4) The above equation is called the m-term KLE of X(t;#). Since the KL coefficientsfZ i g are uncorrelated, the variance of the process is the sum of the eigenvalues. Thus, a common choice for the truncation m would be that more than 95% of the variance in the random process is captured by the first m KL terms or modes. The fundamental formulation shows that KLE approximates a stochastic process as a finite linear com- bination of orthogonal functions. The importance of the KLE is that it yields the best orthonormal basis in the sense that it minimizes the total mean-squared error [140]. 95 5.2 Integration of Karhunen-Loève expansion and basis adaptation This section intends to integrate the KLE and basis adaptation method to the probabilistic modeling of complex problems with high-dimensional QoI instead of a scalar. The KLE is used to decrease the dimension of the QoI, while the basis adaptation provides dimension reduction of the forward probabilistic modeling. Suppose the QoI in our problems is a processH(!;), for example, one FRF of a system with!2R n the vector of discretized frequency points in the interested frequency bandB = 2]0;f u ] and with2R d the vector of standard Gaussian variables. Then by Eq (5.4),H(!;) can be approximated by H(!;) H (!) + m X i=1 p i Z i ()e i (!); (5.5) where m is the number of KL terms considered. To build expansion (5.5), the covariance matrix of K H (! 1 ;! 2 ) = Cov (H(! 1 ;);H(! 2 ;)) need either be a known function or evaluated numerically. For the FRF we are considering, the response at each frequency is a combination of various eigenmodes of the system, and the contributions of all the modes can be very different at that frequency. This, in turn, leads to the complexity of the correlation among the frequencies, andthus, anexponentialdecayingcorrelationkernel, achoiceformanyapplications, mightnotbeaccurateto capture the correlation. To account for that, we evaluate the covariance of the FRF by a low-level quadrature rule. Recall that in the basis adaptation method in chapter 4, in order to construct the rotation matrix of the QoI, a pilot first order PCE is performed in which only a handful of samples is required, and a specific application is to use the level 1 quadrature rule to compute the associated PCE coefficients. Inspired by that, we can also use the level 1 quadrature rule to approximate the covariance ofH(!;). However, before doing that, we will first use the level 1 quadrature to approximate the mean of theH(!;) as H (!) =E [H(!;)] = Z H(!;)p()d X ( (q) ;wq )2Q d H(!; (q) )w q ; (5.6) 96 where ( (q) ;w q ) are the coordinates and weights that define a d-dimensional level 1 quadrature ruleQ d . Then the covariance ofH(!;) can be approximated as K H (! 1 ;! 2 ) = Cov (H(! 1 ;);H(! 2 ;)) =E [(H(! 1 ;) H (! 1 )) (H(! 2 ;) H (! 2 ))] = Z (H(! 1 ;) H (! 1 )) (H(! 2 ;) H (! 2 ))p()d X ( (q) ;wq )2Q d H(! 1 ; (q) ) H (! 1 ) H(! 2 ; (q) ) H (! 2 ) w q ; (5.7) where H (!) is computed by Eq (5.6). SubstitutingK H to Eq (5.2), the eigenvalue and eigenvector pairs of covariance kernelK H are obtained asf i ;e i (!)g. Then, any FRFH(!;) can be projected to the basis off i e i (!)g and the associated KL parametersfZ i ()g can be computed as Z i () = 1 p i Z H(!;)e i (!)d!; for i = 1;:::;m: (5.8) That is, a mapping fromH(!;)2R n to Z i ()2R m for given is constructed, where n is the number of frequency point considered in the FRF. Note that, although (5.7) is an approximation of the covariance ofH(!;), the expansion (5.5) still has controllable accuracy with respect to m, and it is exact when m!1. However, the accuracy of the covariance approximation affects the eigenvalue decay, and the number of required KL terms will also be affected. Usually, a small number m of KL terms is enough to achieve good accuracy of the KLE. In such circumstances, we will have nm, thus reducing the dimension of the QoI to a large extent. The detailed steps to obtain the KL basis are listed in Procedure 5.1: Then, instead of building a probabilistic model of H(!;) that is n dimensions, a replacement m- dimensional probabilistic model of Z i () can be constructed in order to reduce the computational efforts. Moreover, sincefZ i ()g are uncorrelated variables, it is reasonable to consider constructing m probabilistic models for each KL parameter independently. Although the integration of KLE can reduce the dimension of QoI, the dimension of the input space remains unchanged since the argument2R d is the same for both H(!;) and Z i (). By constructing a probabilistic model independently for Z i () for i = 1;:::;m, the QoI 97 Procedure 5.1: Compute Karhunen-Loève basis of the FRF 1 Generate d-dimensional level 1 Gauss–Hermite quadrature ( (q) ;w q )2Q d ; 2 Map (q) to physical input parameters by (q) =T ( (q) ) withT the same operator as explained in chapter 4; 3 EvaluateM(!; (q) ) and letH(!; (q) ) =M(!; (q) ); 4 ComputeK H by substitutingH(!; (q) ) to (5.7); 5 Obtain the eigenvalue and eigenvector pairs off i ;e i (!)g following (5.2). in each problem is a scalar. Therefore, in the case whered is of high dimension, the basis adaptation method in chapter 4 can be used to reduce the dimension of the forward UQ problem, which in turn can accelerate the probabilistic modeling. The adapted PCE of Z i () can be written as Z i () = X 2J d;p Z i; () = X 2J d;p Z A i i; () Z i ( r ) = X 2Jr;p Z A i r i; ( r ); for i = 1;:::;m; (5.9) whereA i is the rotation matrix associated with Z i , which is different for different KL terms; =A i is the vector of rotated variables byA i ; r =A i r withA i r the firstr rows of the rotation is the vector of the first r variables of;fZ i; g,fZ A i i; g, andfZ A i r i; g are the PCE coefficients ind-dimensional space,d-dimensional space, andr-dimensional r space, respectively. The same level 1 quadrature that is used to construct an approximation of K H can also be used to build the rotation matricesfA i g offZ i ()g. More specifically, onceH(!; (q) ) at quadrature pointsf (q) g are evaluated, they can be mapped to Z i ( (q) ) by (5.8) for i = 1;:::;m. Then,Z i ( (q) ) can be substituted in the first order PCE of Z i to compute the first order PCE coefficients. Finally, to follow the other steps of Algorithm 4.1, one can obtain the rotation matrix A i . In Eq (5.9), the greatest polynomial order p of r can be greater than 1, and based on the complexity of the underlying physical problem, the order can be adjusted. The final dimension r of the basis adaptation can be determined by convergence analysis as discussed in chapter 4. Moreover, the accelerated basis adaptation method can also be adopted if necessary. 98 The next step is to compute the adapted PCE coefficientsfZ A i r i; g in (5.9). An r-dimensional Gauss- Hermite quadrature rule will be applied for this task, and the level, l, of the quadrature depends on the order of the adapted PCE and the computing resources. The detailed steps to compute the coefficients are listed in Procedure 5.2, where, for consistency, the step to build rotation matrices for different KL parameters from the d-dimensional level 1 quadrature is also included. This procedure considers the case where the Procedure 5.2: Construct surrogate model of Karhunen-Loève parameters. 1 for i = 1 to m do 2 Contruct rotation matrixA i from level 1 d-dimensional quadrature generated from Procedure 5.1 by following steps in Algorithm 4.1; 3 Generate r-dimensional level l Gauss-Hermite quadrature of ( (q) r ;w q )2Q r ; 4 for i = 1 to m do 5 Map (q) r to full dimensional variables by (q) = (A i r ) T (q) r ; 6 Map (q) to physical input parameters by (q) =T ( (q) ); 7 EvaluateM(!; (q) ) and letH(!; (q) ) =M(!; (q) ); 8 Compute KL parameters at the quadratures by (5.8) as Z i (q) = 1 p i Z H !; (q) e i (!)d!; 9 Substitute the Z i (q) into (5.9) and compute the PCE coefficients by quadratures as Z A i r i; = D Z i ; A i r E X ( (q) r ;wq )2Qr Z i (q) (q) r w q ; (5.10) with (q) = (A i r ) T (q) r . input space is of high dimension. Since different rotation matricesfA i g are involved in the computation, the surrogate modeling has to be done for each KL parameter separately. Note that, if the dimension is manageable and the basis adaptation is unnecessary, the procedure can be replaced by the normal process of constructing PCEs for every KL parameter. In the latter case, the same quadrature points can be used for all the KL parameters, and thus there is no need to construct the surrogate models of the KL parameters independently form times. Instead, once the physical model is evaluated at the quadrature points, resulting inH(!; (q) ), they can be substituted into (5.8) to obtain the KL parameters at the quadrature points, 99 Z i ( (q) ) for i = 1;:::;m. Then, Z i ( (q) ) are substituted into 5.9 to compute the PCE coefficients Z i; for i = 1;:::;m, and the surrogate PCE models offZ i g i=1;:::;m can be built. Suppose PCE coefficientsfZ i; g or adapted PCE coefficientsfZ A i r i; g in (5.9) are obtained, then the surrogate models forfZ i ()g are constructed. The last step for the probabilistic modeling of processH(!;) istosubstitutethePCEoradaptedPCEofZ i ()of (5.9)totheKLEofH(!;)of (5.5)toobtainasurrogate model ofH(!;). Then, given any arbitrary input , the associated FRFH(!;) can be obtained. The surrogate model ofH(!;) acts as an operator that maps toH, that isH : R d 7!R n , more specifically, H :7!H(!;). 5.3 Incorporation of nonparametric stochastic analysis The probabilistic modeling introduced, including the integration of KLE, is a parametric approach, which has broad applications and can capture the model parameter uncertainties for the underlying computational models. However, parametric approaches are, in general, unable to represent the uncertainties induced by modeling errors (or background noises). For example, the parametric approach is restricted by the admissible space of the randomized parameters, and the span of the matrices that take the parameters as arguments, for example,P(), might not cover all the possible values of these matrices. Also, the assumptions, for example, the connections are modeled by springs in the FLSNFC, and other model assumptions cannot be explained by parametric stochastic analysis. Hence, the nonparametric stochastic analysis will be introduced in this section, and the fundamentals of the method will be presented with the incorporation of parametric and nonparametric approaches followed. 5.3.1 Fundamentals of the nonparametric stochastic analysis Consider a linear dynamic system in the frequency domain introduced in chapter 2 ! 2 [M] +i! [D] + [K] U(!) =F(!); (5.11) 100 where [M], [D], and [K] are the FE mass, damping, and stiffness matrices, which, without loss of generality, are considered to be symmetric and positive definite;F(!) is the discretized external force at frequency !; i denotes the imaginary unit. This model is referred to as the mean computational model. Suppose that a ROM is constructed by using ROB [ r ] =f 1 ;:::; nr g, which is a set of eigenmodes of the mean computational model, the number of the modes isn r . For example, ROB [ r ] can be a collection of the important modes obtained by the CB substructure technique and automatic mode selection (introduced in chapter 3). Now, introducing the reduced matrices of mass, damping, and stiffness [M] = [ r ] T [M][ r ]; [D] = [ r ] T [D][ r ] and [K] = [ r ] T [K][ r ]; (5.12) These matrices are symmetric and positive definite for any fixed structure. Thus, Cholesky decomposition of [M], [D], and [K] can be performed and obtain [M] = [L M ] T [L M ]; [D] = [L D ] T [L D ] and [K] = [L K ] T [L K ]; (5.13) where [L M ], [L D ], and [L K ] are upper triangular matrices. Then the projected model is called a ROM of the original model. Instead of randomizing the input parameters of the system and propagating the uncertainties to the generalized mass, damping, and stiffness matrices, the nonparametric method intends to construct random uncertainty models for the generalized mass, damping, and stiffness matrices directly. Suppose there are uncertainties in the dynamic system, and the random mass, damping, and stiffness matrices associated with [M], [D], and [K] are [M], [D], and [K], respectively. Given that [M], [D], and [K] are positive definite, [M], [D], and [K] must also be positive definite. Assume the means of the random reduced matrices are given by the corresponding means of the reduced matrices, that is Ef[M]g = [M]; Ef[D]g = [D] and Ef[K]g = [K]: (5.14) 101 Then from Eqs (5.13) and (5.14), the random matrices can be written as [M] = [L M ] T [G M ][L M ]; [D] = [L D ] T [G D ][L D ] and [K] = [L K ] T [G K ][L K ]; (5.15) where [G M ], [G D ], and [G K ] are also positive definite, and thus invertible. The mean values of [G M ], [G D ], and [G K ] are Ef[G M ]g = [I nr ]; Ef[G D ]g = [I nr ]; and Ef[G K ]g = [I nr ]; (5.16) with [I nr ] the (n r n r ) unity matrix. Then, the associated [G] matrices can be randomized and carry the uncertainties instead of directly randomizing the mass, damping, and stiffness matrices. Furthermore, since the [G] matrices for mass, damping, and stiffness matrices have the same expected values, the identity matrices, the randomization of these matrices can be dealt with in a unified way. The next step is to build probabilistic models for positive definite matrices [G M ], [G D ], and [G K ]. Suppose that [C] is one of [M], [D], and [K], and [G] is one of [G M ], [G D ], and [G K ] with [C] = [L C ] T [G C ][L C ]. With the available information and by using the maximum entropy principle [74], Soize [75, 76] proved that the probability density function of [G] can be obtained as p [G] ([G]) =c G (det([G])) (nr +1) 1 2 C 2 2 C exp n r + 1 2 2 C tr[G] ; 8 [G] that is positive definite. (5.17) In Eq (5.17), “tr” denotes the trace of the matrices. The hyperparameter C allows one to control the level of uncertainties presented in [G] with C = 1 n E k [G] [I nr ]k 2 F 1=2 ; (5.18) wherekk 2 F = denotes the Frobenius norm of the underlying matrix. The value of C has to be chosen such that 0< C < r n r + 1 n r + 5 : (5.19) 102 In Eq (5.17), the positive constant c G is obtained as c G = (2) nr (nr1)=4 n r + 1 2 2 C nr (nr +1) 2 2 C 8 < : nr Y j=1 n r + 1 2 2 C + 1j 2 9 = ; 1 ; (5.20) where () is gamma function. The generator of realizations of any random matrix [G] that has distributions of (5.17) is given by reference [73, 75]. In the generator, the matrix [G] is written as [G] = [L G ] T [L G ] (5.21) with [L G ] an upper triangular matrix satisfying: 1. random variablesf[L G ] jj 0;jj 0 g are independent; 2. for j < j 0 , [L G ] jj 0 can be written as [L G ] jj 0 = % nr Z jj 0 with % the coupling coefficient (% = 1 by default), nr = C (n r + 1) 1=2 , andZ jj 0 a standard Gaussian variable; 3. for j = j 0 , [L G ] jj = nr p 2V j where V j is a Gamma random variable with shape parameter of nr +1 2 2 C + 1j 2 and scale parameter 1. This generator provides a procedure for the Monte Carlo simulation of random matrix [G] and thus the simulation of random matrix [C]. The same procedure can be applied to random matrices [M], [D], and [K], where the associated random matrices [G M ], [G D ], and [G K ] can be simulated by the generator defined by Eq (5.21). It is proved in [75] that random matrices [G M ], [G D ], and [G K ] are independent. The dispersion or the level of uncertainty of random matrices [M], [D], and [K] are controlled by hyperparameters M , D , and K with 0< M ; D ; K < r n r + 1 n r + 5 ; (5.22) 103 when n r is a large number, the upper bound is approximately 1. For the same dispersion parameter, the effects on problems with different dimensions n r can be different, and in such circumstances, the following parameter can be used to specify a dimension-independent scale [75], C = C r 2 n r + 1 ; (5.23) with 0< C < q 2 nr +5 . For the same parameter C , the dispersion for problems with different dimensions is approximately the same. Therefore, instead of using M , D and K , M , D , and K can be used to control the dispersion of random matrices [M], [D] and [K], and the advantage of which is that they can be used for problems with different dimensions to achieve the same level of uncertainty. 5.3.2 Incorporation of parametric and nonparametric stochastic analysis Although the nonparametric stochastic analysis can represent model parameter uncertainties and uncertain- ties induced by modeling errors simultaneously, to take advantage of both parametric and nonparametric stochastic analysis, these two families of approaches will be integrated to bring a unified probabilistic model for complex models. The parametric approaches are useful when a probabilistic model can be constructed for the considered set of input parameters. One of the biggest challenges of the parametric approaches is that they usually can only deal with low-dimensional problems, and with the increase of dimension, the methods have often suffered the curse of dimensionality. Thanks to KLE and basis adaptation of PCE, we can build parametric probabilistic models for problems with high dimensional inputs and high dimensional outputs, as described in chapter 4 and the first two sections of the current chapter. Hence, the parametric stochastic analysis will be used to capture the uncertainties embedded in the input parameters. On the other hand, the nonparametric stochastic analysis will consider the uncertainties induced by modeling errors, which are not explained by the parametric approach. In the nonparametric approach, the dispersion parameters are usually required to be treated as parameters when data is not available, to perform sensitivity analysis to have a robust stochastic solution as a function of the level of uncertainties, or to be estimated from data when data is available. The 104 advantage of using nonparametric only to capture the uncertainties in modeling errors is that the level of dispersion in the modeling errors is easier to quantify based on expert opinions and experiences, such as experimental experiences. Since the dynamics system in the present study is carried out by the eigenmode superposition method, the nonparametric approach to the eigenvalue analysis is first introduced, then this step will be combined with the parametric analysis. Suppose that the eigenmodes of an investigated system are represented by [ r ] =f 1 ;:::; nr g, which will be used as a ROB in the nonparametric approach, the associated eigenvalues aref 1 ;:::; nr g. For example, in the FLSNFC model, the two-level nested CB method can first be used to obtain CB eigenmodes, these modes can then be projected to obtain the system modes or FE modes, then the automatic mode selection on the resulting modes can be performed to select a handful of n r modes to construct ROB [ r ]. Assume that the modal damping is applied with fixed modal damping ratios, then only the mass and stiffness matrices will be randomized. Also, assume mass normalization is used to obtain the ROB [ r ], then the reduced mass and stiffness matrices can be written as [M] = [ r ] T [M][ r ] = [I nr ] and [K] = [ r ] T [K][ r ] = [ nr ]; (5.24) where [ nr ] is a diagonal matrix with its entriesf 1 ;:::; nr g. Then the randomized mass and stiffness matrices are [M] = [G M ] = [L M ( M )] T [L M ( M )] and [K] = [G K ] = [ nr ] T [L K ( K )] T [L K ( K )][ nr ]; (5.25) where [ nr ] = p [ nr ]. In the above equation, we write the random matrices [L M ] and [L K ] in terms of dispersion parameters M and K , respectively. The eigenvalue problem associated with [M] and [K] is [K][X] = [M][X][]: (5.26) 105 By solving the above problem, we can obtain the eigenvalues as diagonal entries of [] and the eigenvectors as the columns of [X]. Apply the projection [ r ] on [X], and get the FE eigenvector matrix as [] = [ r ][X]: (5.27) To combine the parametric and nonparametric approaches, the parametric stochastic model will need to be constructed, for instance, using the integration of KLE and basis adaptation (see in chapter 4 and the first two sections of the current chapter). The parametric random samples, FE eigenmodes [], will be generated to build the probabilistic model during the construction. For the parametric analysis without integrating the nonparametric analysis, samples of [] will be substituted back to the dynamic equation to obtain the samples of the QoI, for example, the FRF, and these samples will be used to construct the probabilistic model of the QoI. Now, to integrate the nonparametric approach to the parametric approach in order to account for the uncertainties presented in the modeling errors, the generalized mass and stiffness matrices by [] can be randomized as in Eq (5.25). However, one should also note that the direct application of the nonparametric approach to [M] = [] T [M][] and [K] = [] T [K][] could lead to an eigenvalue problem of (5.26) of high dimension when the dimension of [] is high. Thanks to the automatic mode selection as described in chapter 4, the important modes of the FE model could be first selected from [], and a reduced eigenvector matrix with a handful of n r most important modes [ r ] can be built. Then follow Eqs from (5.24) to (5.27), we can easily obtain MC samples of FE modes. Since the parametric approach already has a pool of random samples of [], we choose the number of nonparametric MC samples for each parametric sample as 1. For example, suppose the parametric samples aref[] i g N P i=1 where N P is the number of samples. Then for each parametric sample [] i , the random mass and stiffness matrices are randomized, and only one random sample is drawn, which leads to one sample of the associated randomized eigenvector [] i . Therefore, after applying the nonparametric approach to all parametric samplesf[] i g N P i=1 , the associated randomized samples are f[] i g N P i=1 . The way of integration ensures that the number of output random samples of the FE modes remains the same as the parametric analysis. The new random samples of [] can then be substituted back to the dynamic equation to obtain samples of the QoI, and following the same routine in parametric analysis, 106 the probabilistic model of the QoI can be built. Although the structure of the parametric approach is not changed, the uncertainties induced by modeling errors are considered by adding additional randomness on each parametric sample. Through the integration, the propagation of parameter uncertainties is captured by the well-developed parametric analysis, while the nonparametric approach will only account for the uncertainties from modeling errors. Expert opinions and experiences can determine the dispersion parameters of the nonparametric analysis, the error induced by which will be small and only affect the modeling errors. 5.4 Applications 5.4.1 The fuel assembly The integration of KLE and adapted PCE (or PCE) to the probabilistic modeling described in section 5.2 will be applied to a BWR FA, where the detailed description of the model can be found in chapter 2. The connections between fuel rods and spacer grids and between spacer grids and the FA channel are randomized. For both types of connections, the ones associated with the bottom 4 spacer grids are modeled by one random variable, and the ones associated with the top spacer grids and upper handle structure are modeled by another random variable, which makes the dimension of the UQ problem d = 4. All variables are assumed to haveBe(;) distributions with = 5:0 and = 1:5 and scale factors their associated design values, that is, 0.185 kN/mm and 1.85 kN/mm for the two variables of rod-grid connections and two variables of grid-channel connections, respectively. The choice of shape parameters and makes the coefficient of variance (COV) of the 4 random variables 20%. The mean values of these two sets of random variables are approximately 0.15 kN/mm and 1.5 kN/mm. The PDF of theBe(;) with a scale factor equal to 1.85 kN/mm is presented in Figure 5.1, and the PDFs of the other parameters have the same shape but with different scale factors. The quantity of interest in this problem is the FRF of the bottom tip node with unit excitation in Z direction and observation in Y direction. A free boundary condition is applied to the FA. The frequency band of interest isB =]0; 1000] Hz, where the frequency points are n = 691. 107 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 K (kN/mm) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 PDF beta pdf Figure 5.1: Probability density function of theBe(5:0; 1:5) distribution with scale factor 1.85 kN/mm. Following Procedures 5.1 and 5.2, the intention is to build a surrogate model for the FRF, where KLE first approximates the FRF with finite terms and PCE constructs the surrogate models of KLE parameters. Since the input space dimension is 4 and is manageable by conventional PCE, no adaptation of the KLE parameters is required. After the surrogate model via KLE and PCE is built, its accuracy will be evaluated by comparing it to the direct PCE on the FRF. The first step of Procedure 5.1 is to generate d-dimensional level 1 Gauss-Hermite quadrature points of f (q) g, the number of level 1 quadrature for the 4-dimensional problem is 9. The Gauss quadrature points are then mapped to the physical space variablesf (q) g by inverse transform sampling, on which the FRF is evaluated. Using these quadrature points, we can first present the mean FRF shown in Figure 5.2. Then, the 0 200 400 600 800 1000 Frequency (Hz) 120 100 80 60 40 Acceleration (dB) Figure 5.2: Mean FRF of the FA obtained from the level 1 quadrature. 108 covariance matrix of the FRF is approximated through Eq (5.7), and the eigendecomposition of which yields the eigenvalue and eigenvector pairsf i ;e i (!)g. The eigenvalues of KL term up to 30 are shown in Figure 5.3, from which we can see that the eigenvalue decays very fast, and it is almost 0 at the 9-th term. Thus, in this application, we will choose m = 8 KL terms to approximate the FRF, and the first 8 eigenvectors of the covariance are plotted in a single figure, see in Figure 5.4. It can be seen that there are more variations in the lower frequency region, while the middle- and high-frequency regions only have peaks in the mode shape locally. Also, the contributions to the frequency band are different from mode to mode. 0 5 10 15 20 25 30 KL term 0 500 1000 1500 2000 Eigenvalue Figure 5.3: Eigenvalues of the covariance matrix of the FRF of the FA approximated by level 1 quadrature. 0 200 400 600 800 1000 Frequency (Hz) 0.4 0.3 0.2 0.1 0.0 0.1 0.2 0.3 Mode shape Figure 5.4: The first 8 eigenvectors of the covari- ance matrix of the FRF of the FA approximated by level 1 quadrature. Once the KL basis of the FRF is constructed, the next step is to construct probabilistic models for the 8 KL parameters. When the input space is of high dimension, basis adaptation will be required for each KL parameter, described in Procedure 5.2. Since the dimension is 4, we will use the classical PCE for the task, and sparse quadrature will be used to construct PCEs of the KL parameters of order up to 3. The PCE coefficients,fZ i; g 2J d;p in the expansions in Eq (5.9) need to be evaluated using a new set of quadrature points. More specifically, a level 3 Smolyak sparse quadrature rule is used. The level 3 quadrature pointsf (q) g are first generated and mapped back to the physical parameters spacef (q) g, and the FRF onf (q) g is evaluated and projected to the 8-term KL basis as in Eq (5.8) to obtain the KL parameters at the quadrature points. Then, we can use the KL parameters at the quadrature points to compute the PCE coefficientsfZ i; g 2J d;p , for i = 1;:::;m. The PDFs of the PCE models of order up to 1, up to 2, and up to 3 of the 8 KL parameters are shown in Figure 5.5. We call the PCE model up to order p the “order p 109 2.50.0 2.5 kle coefficient #1 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #2 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #3 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #4 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #5 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #6 0.0 0.2 0.4 0.6 PDF 2.50.0 2.5 kle coefficient #7 0.0 0.2 0.4 PDF 2.50.0 2.5 kle coefficient #8 0.0 0.2 0.4 PDF Figure 5.5: Probability density functions of KL parameters obtained from PCEs of order up to 1 (in blue), up to 2 (in orange), and up to 3 (in green) for the FA. PCE”. In the figures, the order 1 PCEs are in blue, the order 2 PCEs are in orange, and the order 3 PCEs are in green color for all KL parameters. The order 1, order 2, and order 3 PCEs of the KLE parameters are obtained from a level 3 sparse quadrature rule. It can be seen in Figure 5.5 that order 2 and order 3 PCEs converge to each other for all 8 KLE parameters, and this means that order 2 or order 3 PCE is enough to capture the probabilistic behavior of the KL parameters. The convergence here also means that the PCE is accurate. Also, we see that some KL parameters have apparent non-Gaussian behavior by comparing the order 3 PCE to the order 1 PCE, for example, KLE #1, #4, #5, and #6. These non-Gaussian behaviors enable the KLE to represent non-Gaussian behaviors in the FRF. Another way to check the accuracy of the probabilistic model of the KL parameters is to compare the parameters obtained from the PCE models to the ones obtained from the MC simulation of the physical model. Therefore, we randomly generated 500 MC samples from the input space, the physical model is evaluated, and the FRF at the MC samples is obtained, then the FRF is projected to the 8-term KL basis to compute the corresponding KL parameters. Thus, the KL parameters from the MC simulation can be served as reference samples. Using the same input samples, we can also compute the KL parameters by the PCE models, and these parameters are compared to the associated ones by MC. The comparison is shown in 110 Figure 5.6 for all 8 KL parameters. Each red dot in the figures is an MC sample with the X-axis value the 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #1 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #2 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #3 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #4 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #5 2.5 0.0 2.5 Physical eval 5.0 2.5 0.0 2.5 PCE of KLE #6 2.5 0.0 2.5 Physical eval 2.5 0.0 2.5 PCE of KLE #7 2.5 0.0 2.5 Physical eval 0 5 PCE of KLE #8 Figure 5.6: Plot of 500 MC samples of the first 8 KL parameters for the FA with X-axis the values obtained from the physical model evaluations and Y-axis the values obtained from the PCE models; the blue lines in the figures are lines with unit slop. KL parameter obtained from the physical model and the Y-axis value the KL parameter obtained from the underlying PCE model. When the X-value is close to Y-value for a single dot, the KL parameters obtained from the physical and PCE models are close, which means the PCE model is accurate for that sample. Thus, an accurate PCE model of a KL parameter will have the majority of the red dots located around the line with unit slop, shown in blue in the figures. From Figure 5.6, the PCE models for KL parameters from #1 to #5 show good accuracy with most samples located around the unit 1 slope lines with narrow scattering. The MC samples of the remaining KL parameters are also located around the unit 1 slop lines but with greater scatters. The results here also indicate the accuracy of the PCE models of the KL parameters. Moreover, the lower KL terms that are supposed to be more critical with greater eigenvalues have higher accuracy than the higher KL terms. Therefore, we will use the order 3 PCE for all 8 KLE parameters as their probabilistic model. These probabilistic models of the KL parameters can then be substituted back to the KLE (5.5) to obtain a probabilistic surrogate model of the FRF. Note that the level 3 quadrature used to construct the PCE models of the KL parameters can also be used to construct a PCE model of the FRF directly. For 111 comparison purposes, we will construct an order 3 PCE of the FRF as a reference, referred to as “PCE of the FRF”. The previous surrogate model of the FRF, since it also involves PCEs of the KL parameters, will be referred to as “PCE of the FRF via KLE”. Note that the construction of KLE of the FRF is not necessary for practice when the dimension of the input space is small, and one could use sparse quadrature to build a PCE model of the QoI directly like in the reference model here. We take a detour to the KLE to show the capability of output dimension reduction, which will be beneficial in problems with high dimensional input space. From the PCE of the FRF via the KLE model, we can generate the PDFs of the FRF at each frequency point, and we plot the 68% and 95% intervals (CI) of the model for presentation purposes in Figure 5.7. The PDF will be compared to the direct PCE model to check the accuracy of the PCE of the FRF via 0 200 400 600 800 1000 Frequency (Hz) 140 120 100 80 60 40 Acceleration (dB) mean 95% CI 68% CI Figure 5.7: Confidence intervals of the probabilistic model of the FRF of the FA obtained from PCE via KLE; the green solid line presents the mean FRF. the KLE model. We are more interested in the peaks of the FRF, and thus the comparison of the FRF at these peaks is presented in Figure 5.8. The PDFs in blue are obtained from the PCE of the FRF directly and are served as references, while the PDFs in orange are obtained from the PCE of the FRF via KLE. The presenting frequencies contain all the major peaks along with the whole frequency band. For all the presenting frequencies, the PCE model of the FRF via KLE has good accuracy compared to the reference PCE model, in which the range is consistent, the shape is about the same, and the non-Gaussian behaviors are also captured. Then, for any given input parameters in the support of the distributions of the input parameters, we can use the surrogate model of the FRF to obtain the associated FRF. 112 125 100 75 FRF at 35.5 Hz 0.00 0.02 0.04 PDF 100 80 FRF at 93.5 Hz 0.00 0.02 0.04 0.06 PDF 75 70 65 FRF at 183.3 Hz 0.0 0.1 PDF 42 40 FRF at 235.8 Hz 0.0 0.5 1.0 PDF 65 60 55 FRF at 321.6 Hz 0.0 0.2 0.4 PDF 45.0 42.5 FRF at 368.1 Hz 0.0 0.5 1.0 PDF 45 40 FRF at 437.1 Hz 0.0 0.2 0.4 PDF 34 33 FRF at 641.6 Hz 0 1 2 PDF 35 30 FRF at 699.1 Hz 0.0 0.2 0.4 PDF 42.5 40.0 37.5 FRF at 800.2 Hz 0.0 0.2 0.4 0.6 PDF 60 40 FRF at 869.4 Hz 0.00 0.05 0.10 PDF 57.5 55.0 52.5 FRF at 910.3 Hz 0.0 0.2 PDF Figure 5.8: Comparison of the FRF of FA at the frequencies around peaks; the blue lines are obtained from PCE of the FRF and the orange lines are obtained from PCE of the FRF via KLE. 5.4.2 The fully-loaded spent nuclear fuel canister model The integration of the KLE and the adapted PCE is now applied to the FLSNFC model, with the nonpara- metric stochastic model incorporated to account for the uncertainties in the modeling errors. As described before, the FLSNFC model has three major structural levels, the canister-basket structure, the 68 FAs, and numerous fuel rods. There are three types of connections associated with the structural levels: theconnectionsbetweenfuelrodsandspacergridsintheFA,theconnectionsbetweenspacergridsand FA channel, and the connections between FA channel and basket structure. The three types of connections are referred to as rod-grid, grid-channel, and channel-basket connections, respectively, or simply type I, type II, and type III connections. For each of the 68 FAs, these three types of connections are randomized, and the same type of connections associated with one FA share the same random variable, and thus, there are 113 3 random variables for each FA, and there are d = 3 68 = 204 random variables in total for the FLSNFC model. The random variables are assumed to haveBe(;) distributions with = 5:0 and = 1:5 and scale factors their associated design values, that is 0.185 kN/mm, 1.85 kN/mm, and 100 kN/mm for the three types of connections, respectively. Same as in the FA model, the shape parameters and make the COV of the 204 random variables 20%. The mean values of the random variables associated with the three types of connections are approximately 0.15 kN/mm, 1.5 kN/mm, and 80 kN/mm. The quantity of interest in this problem is the FRF at the center of the exterior of the canister bottom plate with unit excitation in Z direction and observation in Y direction. The canister model is hung at the top. The frequency band of interest isB =]0; 1000] Hz, where the frequency points are n = 691. Nonparametric stochastic analysis. Before diving into the parametric stochastic analysis, the nonpara- metric approach for the FLSNFC model is first considered, which will be used to account for uncertainties in the modeling errors. The modeling errors or background noises in the FLSNFC include several assumptions, for example, the connections between different structural levels are assumed to be modeled by linear springs, and some unobserved modeling errors. For the nonparametric analysis, we first need to use the eigenmodes of the system to propose a ROB. Using the two-level nested CB model with dominant FA component modes, system modes are n e = 232; 586. However, not all these modes have significant contributions to the FRFs at the canister bottom plate. Therefore, an automatic mode selection of the system modes is performed to reduce the size of the later nonparametric analysis, which suggests that with n r = 14; 500 modes kept in the calculation, the accuracy of the FRFs at the canister bottom plate is comparable when all the system modes are kept. Thus, these n r modes build a ROB [ r ] for the nonparametric analysis. Then the reduced mass and stiffness matrices, [M] and [K], can be constructed using Eq (5.24), then the randomized mass and stiffness matrices, [M] and [K], can also be computed by Eq (5.25). The random samples of [M] and [K] can then be obtained following the process under Eq (5.21), then the eigenvalue problem for samples of [M] and [K] can be solved, and the samples of the system modes [] can be obtained by Eqs (5.26) and (5.27). The samples of these modes can then be used to construct probabilistic models of the interested FRF. This application will use the same dispersion parameters for [M] and [K] to let M = K = C . Using C = 0:8515 corresponds to let the dependent dispersion parameter C = 0:01, the probability density 114 function of the FRF can be computed by using 136 nonparametric random samples (which is verified to be enough to have converged results). The 95% CI of the FRF is plotted in Figure 5.9 on the left, while a zoom-in figure in the frequency band of [200, 800] Hz is shown on the right. The deterministic reference FRF with 20 kHz truncation (in black) and the FRF obtained by the most dominant 14; 500 system modes (in blue) are also presented. It can be seen that FRF obtained from the dominant modes is almost identical to the reference FRF. Moreover, with the increase of frequency, the level of dispersion on FRF increases. M = 0.8515, K = 0.8515, = 0.01 0 200 400 600 800 1000 Frequency (Hz) -150 -140 -130 -120 -110 -100 -90 -80 -70 Acceletation (dB) SROM 14500 Ref Dmt 14500 modes (a) in the interested frequency band M = 0.8515, K = 0.8515, = 0.01 200 300 400 500 600 700 800 Frequency (Hz) -150 -140 -130 -120 -110 -100 -90 -80 -70 Acceletation (dB) SROM 14500 Ref Dmt 14500 modes (b) in the frequency band of [200, 800] Hz Figure 5.9: The 95% confidence interval (pink shades) of the probabilistic model of the FRF obtained from nonparametric analysis; the reference deterministic FRF (in black) and the FRF obtained by the most dominant 14; 500 modes (in blue) are also presented. In this case, the dimension of the eigenvalue problem of [M] and [K] is n r = 14; 500, which is still not cheap for direct solvers. To further reduce the cost, we will truncate the dimension of this eigenvalue problem by restricting the number of modes considered in the ROB to reduce further the dimension of [ r ]. By gradually decreasing the dimension of the ROB, we found that by using 327 most dominant modes in the ROB, the resulting probabilistic model of the FRF is comparable to the one shown in Figure 5.9. The 95% CI of the FRF obtained by nonparametric analysis with the most 327 dominant modes randomized is plotted in Figure 5.10 on the left, while a zoom-in figure in the frequency band of [200, 800] Hz is shown on the right. The dimension independent dispersion parameter, C = 0:01, is the same as the previous case, which guarantees the same level of uncertainty in the randomized model. Comparing Figure 5.10 to Figure 5.9, the 95% CI is almost the same in the whole frequency band, which is clearer if we compare the zoom in figures on the right. This means randomizing the 327 most dominant modes of the system is comparable to 115 M = 0.12806, K = 0.12806, = 0.01 0 200 400 600 800 1000 Frequency (Hz) -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 Acceletation (dB) SROM 327 Ref Dmt 14500 modes (a) in the interested frequency band M = 0.12806, K = 0.12806, = 0.01 200 300 400 500 600 700 800 Frequency (Hz) -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 Acceletation (dB) SROM 327 Ref Dmt 14500 modes (b) in the frequency band of [200, 800] Hz Figure 5.10: The 95% confidence interval (pink shades) of probabilistic model of the FRF obtained from nonparametric analysis with the most 327 system modes randomized; the reference deterministic FRF (in black) and the FRF obtained by the most dominant 14; 500 modes (in blue) are also presented. randomizing the 14500 most dominant modes, and thus the nonparametric analysis can be restricted to the 327 most dominant modes. The computation effort of 327-dimensional eigenvalue problems is negligible. The next step is to determine the dispersion parameter C to develop a probabilistic model for model error uncertainty. By the effort of the two-level nested CB and the comparison to the reference model, it is confident to say that more than 99% of the interested FRFs at the canister bottom plate have errors less than 0.5 dB, which is a small error. Thus, the level of dispersion presented in Figures 5.9 and 5.10 is apparent too large since the spread of the 95% CI in these two models is much greater than the FRF errors of 0.5 dB. However, at the same time, the nonparametric approach also needs to account for some model assumptions; for example, the connections between different structural levels are assumed to be springs. Furthermore, the uncertainties brought by the model assumptions are assumed to be limited. Based on these considerations, the expert opinions, and the experience in the experiment study, we choose a dispersion parameter of C = 0:005, which corresponds to M = K = 0:064 for the eigenvalue problem of dimension 327. The 95% CI of the resulting probabilistic model of the FRF is shown in Figure 5.11. The level of dispersion in this figure is much less than in Figure 5.10, but it is noticeable dispersion, especially in the high-frequency range. This is because the high-frequency signals are usually easier to be subjected to errors. Theuncertaintiescapturedbythenonparametricapproachwillbeanadditivetotheparametricstochastic analysis. The additional cost associated with this approach is the automatic mode selection of the system 116 M = 0.064031, K = 0.064031, = 0.005 0 200 400 600 800 1000 Frequency (Hz) -150 -140 -130 -120 -110 -100 -90 -80 -70 Acceletation (dB) SROM 327 Ref Dmt 14500 modes Figure 5.11: The 95% confidence interval (pink shades) of probabilistic model of the FRF obtained from nonparametric analysis with the most 327 system modes randomized and dispersion parameter C = 0:005. modesandtheeigenvalueproblemofdimension327. However, thelatterisnegligible, andtheformerrequires 5 minutes by using 60 computer nodes. Parametric stochastic analysis. The parametric approach is mainly the combination of KLE and adapted PCE. The 204-dimensional level 1 Gauss-Hermite quadrature points off (q) g are first generated, the number of samples is 409. These samples in the Gaussian space are mapped to the physical space variables f (q) g, and the FRF is evaluated. Then, a nonparametric analysis is carried out, in which the automatic mode selection of the system modes is performed to incorporate the uncertainties in the model errors for each sample of the FRF, and an eigenvalue problem of the random mass and stiffness matrices, [M] and [K], with one single random sample, is solved. For each parametric sample of the FRF, the output after the nonparametric analysis is a single FRF with uncertainties from the model errors incorporated. Thus, the number of samples remains the same, and we can proceed with the parametric stochastic analysis. By the level 1 quadrature, the mean FRF of the model can be first obtained and shown in Figure 5.12. Compared to Figure 5.2, there are fewer peaks in Figure 5.12, because the canister-basket structure filters many local vibrations. Then, the covariance matrix of the FRF is approximated by Eq (5.7), the eigenvalue and eigenvector pairs,f i ;e i (!)g, of it, can be obtained. The eigenvalues up to 100 KL terms are shown in Figure5.13. ThefigureshowsthattheeigenvaluedecaysfastwiththeKLterm, andtomakeanaccuracy-cost balance, we will choosem = 20 KL terms to build high order probabilistic models. The first 20 eigenvectors 117 0 200 400 600 800 1000 Frequency (Hz) 180 160 140 120 100 80 Acceleration (dB) Figure 5.12: Mean FRF obtained from the level 1 quadrature. of the covariance are presented in Figure 5.14. Again, some regions have more variations than others, and the contributions to the mode shape along the frequency band are different from mode to mode. 0 20 40 60 80 100 KL term 0 2500 5000 7500 10000 12500 15000 17500 20000 Eigenvalue Figure5.13: Eigenvaluesofthecovariancematrix of the FRF approximated by level 1 quadrature. 0 200 400 600 800 1000 Frequency (Hz) 0.4 0.2 0.0 0.2 0.4 Mode shape Figure 5.14: The first 20 eigenvectors of the co- variance matrix of the FRF approximated by level 1 quadrature. OncetheKLbasisisbuilt,thenextstepistoconstructprobabilisticmodelsforthefirst20KLparameters. Since the input space dimension is d = 204 and the computational model still requires about 40 minutes with 60 computer nodes, it is impossible to construct a classical PCE model, and thus the basis adaptation will be used to reduce the dimension of the UQ problems. The basis adaptation is suitable for scalar QoI, and thus each of the KL parameters requires an independent adapted model. The detailed steps to construct adapted PCE of the KL parameters are described in Procedure 5.2. 118 In some detail, the samples of the FRF at the level 1 quadrature points are incorporated with the nonparametric approach to add the uncertainties induced by model errors. Then the resulting samples of the FRF are projected to the KL basis to obtain the KL parameters by Eq (5.8). Then for i = 1;:::;m, the first order PCE of the i-th KL parameter can be built, and the first order PCE coefficients can be used to construct a rotation matrix,A i , adapted to the KL parameter. Next, to construct anr-dimensional adapted PCE model of order up to 3, the r-dimensional level 2 sparse Smolyak quadrature points (restricted by computational resources) are generated for the adapted variablesf (q) r g. Then, for i = 1;:::;m,f (q) r g is mapped back to the full-dimensional variablesf (q) g through the associated rotation matrixA i , the physical model of the FRFH( (q) ;!) are evaluated and projected toZ i ( (q) ), the adapted PCE coefficients are then computed by Eq (5.10), and the adapted PCE model of thei-th KL parameter is accomplished. The classical basis adaptation method is first applied to the construction of adapted PCEs. The PDFs of the 2d and 3d adapted PCE models of the first 20 KL parameters are shown in Figure 5.15, where the 2d adapted PCEs are shown in orange, and the 3d adapted PCEs are shown in green color. For comparison purposes, the 204-dimensional order 1 PCEs of the KL parameters computed by the 204-dimensional level 1 quadrature are also presented in the figures in blue color. We can see that apparently, the 2d and 3d adaptations are not converged, so the increase of dimension is required. Moreover, the higher-order 2d and 3d adaptations differ significantly from the 204-dimensional order 1 PCE. Note that we start the adaptation from dimension 2 since the sparse quadrature starts from dimension 2. Also, the sparse Smolyak quadrature has a nested structure; that is, the r-dimensional quadrature points are included in the (r + 1)-dimensional quadrature points, and thus, the computational cost can be reduced by not repeating the evaluation of the physical model at lower dimensional quadrature points. For example, for a level 2 sparse quadrature rule, 2d quadrature contains 21 nodes, and 3d quadrature contains 37 nodes (the 21 2d quadrature points are included). Hence, for each KL parameter, the 2d and 3d adaptation models require 37 times evaluation of the physical model. Thus, the total number of model evaluations for the 20 KL terms is 37 20 = 740. We then apply the accelerated basis adaptation described in chapter 4 to build adapted PCE models for the KL parameters. More specifically, the mean and Gaussian coefficient correction is first applied, and the 119 2 0 2 kle coefficient #1 0 1 2 PDF 2 0 2 kle coefficient #2 0 2 4 2 0 2 kle coefficient #3 0 5 2 0 2 kle coefficient #4 0 5 2 0 2 kle coefficient #5 0 2 4 2 0 2 kle coefficient #6 0 2 4 PDF 2 0 2 kle coefficient #7 0.0 2.5 5.0 2 0 2 kle coefficient #8 0.0 0.5 1.0 2 0 2 kle coefficient #9 0 2 4 2 0 2 kle coefficient #10 0.0 2.5 5.0 2 0 2 kle coefficient #11 0 5 PDF 2 0 2 kle coefficient #12 0 2 2 0 2 kle coefficient #13 0 2 2 0 2 kle coefficient #14 0.0 2.5 5.0 2 0 2 kle coefficient #15 0 2 4 2 0 2 kle coefficient #16 0.0 2.5 5.0 PDF 2 0 2 kle coefficient #17 0 1 2 2 0 2 kle coefficient #18 0 1 2 0 2 kle coefficient #19 0.0 2.5 5.0 2 0 2 kle coefficient #20 0 5 Figure 5.15: Probability density functions of KL parameters obtained from classical 2d (in orange) and 3d (in green) adaptations, and 204-dimensional order 1 PCEs (in blue). PDFs of the resulting adapted PCEs are shown in Figure 5.16. Applying the accelerated basis adaptation method, the PDFs of the 2d and 3d adapted PCEs (the orange and green curves) converge for all the KL parameters, see Figure 5.16. This means that 2d adapted PCEs are enough to capture the probabilistic behavior of the KL parameters, although the input space is of dimension 204, which is verified by the convergence to the 3d adapted PCEs. Furthermore, by comparing to the order 1 PCE, some KL terms present non-Gaussian behaviors, for example, KL terms #1, #2, #3, #8, #9, #10, #13, #17 and #18. These non-Gaussian behaviors enable the KLE to represent non-Gaussian behaviors in the FRF. The PCE coefficients of the accelerated 3d adaptation models of the 20 KL parameters are shown in Figure 5.17, and the 0-th order PCE coefficients, which represent the mean values of the KL parameters, are not shown in the figures. The first three coefficients correspond to the first order coefficients of the 120 2.5 0.0 2.5 kle coefficient #1 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #2 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #3 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #4 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #5 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #6 0.0 0.2 0.4 PDF 2.5 0.0 2.5 kle coefficient #7 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #8 0.0 0.5 PDF 2.5 0.0 2.5 kle coefficient #9 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #10 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #11 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #12 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #13 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #14 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #15 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #16 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #17 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #18 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #19 0.00 0.25 0.50 PDF 2.5 0.0 2.5 kle coefficient #20 0.0 0.2 0.4 PDF Figure 5.16: Probability density functions of KL parameters obtained from accelerated 2d (in orange) and 3d (in green) adaptations, and 204-dimensional order 1 PCEs (in blue). first, second, and third adapted variables ( 1 , 2 , and 3 ). From Figure 5.17, for all the KL parameters, the first order coefficient associated with 1 is dominant over all coefficients, which means that the first order information is the most critical in the KL parameters and that the first adapted variable 1 is the most important variable. However, there are also noticeable high order contributions; for example, for most KL parameters, the second order term of 1 (the fourth coefficient) has noticeable contributions. Thus, the dominance of the first adapted variable proves from another perspective that low dimensional adaptations of the KL parameters are possible. The 3d adapted PCE model will be used for the first 20 KL terms. Although the other KL parameters have no high order probabilistic model, the first order PCE can still be constructed using the 204-dimensional level 1 quadrature. Then, a 100-term KLE surrogate model of the FRF can be constructed by substituting 121 5 10 15 kle coefficient #1 0.50 0.25 0.00 PCE coefficient 5 10 15 kle coefficient #2 0.5 0.0 5 10 15 kle coefficient #3 0.5 0.0 5 10 15 kle coefficient #4 0.0 0.5 5 10 15 kle coefficient #5 0.0 0.5 5 10 15 kle coefficient #6 0.5 0.0 PCE coefficient 5 10 15 kle coefficient #7 0.0 0.5 5 10 15 kle coefficient #8 0.5 0.0 5 10 15 kle coefficient #9 0.5 0.0 5 10 15 kle coefficient #10 0.0 0.5 5 10 15 kle coefficient #11 0.0 0.5 PCE coefficient 5 10 15 kle coefficient #12 0.0 0.5 5 10 15 kle coefficient #13 0.0 0.5 5 10 15 kle coefficient #14 0.0 0.5 5 10 15 kle coefficient #15 0.0 0.5 5 10 15 kle coefficient #16 0.0 0.5 PCE coefficient 5 10 15 kle coefficient #17 0.0 0.5 5 10 15 kle coefficient #18 0.5 0.0 5 10 15 kle coefficient #19 0.5 0.0 5 10 15 kle coefficient #20 0.5 0.0 Figure 5.17: PCE coefficients of the accelerated 3d adaptation models of the KL parameters. the probabilistic model of the KL parameters, with the first 20 KL parameters the order 3 adapted PCE and the remaining 80 KL parameters the order 1 PCE, to the KLE of the FRF (5.5). Again, we call this surrogate “PCE of the FRF via KLE” model, which, by construction, already considered the uncertainties of the model errors through additive nonparametric analysis on parametric samples. From this model, we can generate the PDFs of the FRF at each frequency point, and we plot the 68% and 95% intervals (CI) of the model for presentation purposes in Figure 5.18. The spreads of the CIs are more remarkable than in the FA case in Figure 5.7; this is reasonable because the input space is much greater, hence more sources of uncertainties. In addition, the uncertainties due to model errors are also included in the FLSNFC model. Unlike the FA model, where a reference order 3 PCE of the FRF can be computed by a level 3 sparse quadrature rule, it is almost impossible to construct an order 3 PCE through any sparse quadrature for the 122 0 200 400 600 800 1000 Frequency (Hz) 180 160 140 120 100 80 60 40 Acceleration (dB) mean 95% CI 68% CI Figure 5.18: Confidence intervals of the probabilistic model of the FRF obtained from PCE via KLE; the green solid line represents the mean FRF. FLSNFC model due to the high dimensional input space of d = 204. In order to make a comparison, the construction of adapted PCE models of the FRF at several frequencies along with the frequency band can be carried out. The frequencies of interest are 40.5 Hz (at around the highest major peak in low frequency), 130.4 Hz (at a low frequency that has a widespread CI), 160.6 Hz (at around a major peak in low to middle frequency), 460.5 Hz (at around a peak in middle frequency) and 669.7 Hz (at around a major peak in high frequency). These adapted PCE models are performed on the FRF directly, not through the KLE, thus can be served as reference models for these particular frequencies. Similar to the adaptations of the KL parameters, each of the FRF at one particular frequency is a QoI, and a basis adaptation needs to be constructed for the QoI. In order to do that, the 204-dimensional level 1 quadrature is first used to construct rotation matrices for the FRF at the selected frequencies. Then the r-dimensional level 2 sparse quadrature points are sampled for the adapted variablesf (q) r g. Then for each selected frequency,f (q) r g are mapped to full-dimensional variablesf (q) r g by the associated rotation matrix, the physical model of the FRFH( (q) ;!) are evaluated, in particular, the FRF at the selected frequency is evaluated at the quadrature points. Next, the adapted PCE coefficients for the FRF at selected frequencies can be obtained, the PCE models adapted to FRFs at these frequencies are built. Again, we first use the classical basis adaptation, and the PDFs obtained from the adapted models of the FRFs at the selected frequencies are shown in Figure 5.19. The PDFs obtained from the adaptations of the FRF directly are shown in solid lines with 2d adaptation in blue and 3d adaptation in orange, while the PDF obtained from the 3d adapted PCE of the FRF via KLE 123 90 85 80 75 70 FRF at frequency 40.53 Hz 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 PDF 2d adapt 3d adapt 3d adapt via KLE 100 80 60 FRF at frequency 160.64 Hz 0.00 0.05 0.10 0.15 0.20 PDF 2d adapt 3d adapt 3d adapt via KLE 105 100 95 90 85 80 FRF at frequency 460.52 Hz 0.0 0.5 1.0 1.5 2.0 PDF 2d adapt 3d adapt 3d adapt via KLE 90 80 70 60 FRF at frequency 669.73 Hz 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 PDF 2d adapt 3d adapt 3d adapt via KLE 140 120 100 FRF at frequency 130.36 Hz 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 PDF 2d adapt 3d adapt 3d adapt via KLE Figure 5.19: PDFs of the FRFs at selected frequencies obtained from classical 2d (in blue) and 3d (in orange) adaptations, and from 3d adaptation via KLE (in dashed green). is shown in dashed line in green color, and all the PCEs have ordered up to 3. Thus, we see that the 2d adaptation and 3d adaptation are not converged, and further increasing of adapted dimension is required by classical adaptation. Also, the PDFs obtained from surrogate 3d adapted PCE via KLE models are not close to the PDFs obtained from the direct adaptation of the FRF. The accelerated basis adaptation is then applied to the direct adaptation of the FRF models; more specifically, the mean and Gaussian coefficients are corrected. The resulting PDFs are shown in Figure 5.20. In this case, the 2d and 3d adaptations of the FRF are consistent with each other, meaning that 2d or 3d adaptation of the FRF is enough to capture the probabilistic behavior of the FRF at selected frequencies. In addition, the PDFs obtained from the surrogate 3d adapted PCE of the FRF via KLE are also close to the PDFs obtained from direct adaptations of the FRF. This verifies the accuracy of the surrogate model at these frequencies. Therefore, since the selected frequencies are representative of the FRF of the frequency band, we can expect that the surrogate model is accurate along with the interested frequency band. The reason to choose only a limited number of frequencies for verification is the restriction of the expensive computational 124 90 85 80 75 70 65 FRF at frequency 40.53 Hz 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 PDF 2d adapt 3d adapt 3d adapt via KLE 100 80 60 FRF at frequency 160.64 Hz 0.00 0.01 0.02 0.03 0.04 0.05 PDF 2d adapt 3d adapt 3d adapt via KLE 110 100 90 80 FRF at frequency 460.52 Hz 0.00 0.02 0.04 0.06 0.08 PDF 2d adapt 3d adapt 3d adapt via KLE 90 80 70 60 50 FRF at frequency 669.73 Hz 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 PDF 2d adapt 3d adapt 3d adapt via KLE 160 140 120 100 80 FRF at frequency 130.36 Hz 0.00 0.01 0.02 0.03 0.04 PDF 2d adapt 3d adapt 3d adapt via KLE Figure 5.20: PDFs of the FRFs at selected frequencies obtained from accelerated 2d (in blue) and 3d (in orange) adaptations, and from 3d adaptation via KLE (in dashed green). cost. For example, the 2d and 3d adaptation at one particular frequency requires evaluating the physical model 37 times (as explained before), with five selected frequencies, the number of total evaluations is already 37 5 = 185 times. In this application, all the techniques described in the current chapter are applied. The KLE of the FRF has successfully reduced the dimension of the output from n = 691 to m = 20. Thanks to the limited dimensionofQoI,theacceleratedbasisadaptationmethodcanbeappliedtobuildprobabilisticmodelsforthe KL parameters. The accuracy of the adapted models is verified by the convergence of the adapted dimension and by the examination of the adapted PCE coefficients. These probabilistic models of the KL parameters with high order information can then be substituted to the KLE of the FRF to construct a surrogate model of the FRF, which is the so-called adapted PCE of the FRF via KLE. The accuracy of the surrogate model is verified by two aspects, with the first one inherited from the accuracy of the probabilistic model of the KL parameters and the second one by comparison to reference adaptation models that are adapted to the FRF directly at selected frequencies. The physical model evaluations are embedded in constructing the first order 125 204-dimensional PCE and basis adaptations of the KL parameters, leading to 409 + 740 = 1149 physical model evaluations. With these evaluations, a probabilistic model with good accuracy for the problem with input dimension d = 204 and output dimension n = 691 is accomplished. To conclude the chapter, the integration of the KLE and the adapted PCE is introduced to deal with stochastic reduction for problems with high-dimensional QoI. The nonparametric stochastic analysis is in- corporated into the parametric approach to account for uncertainties induced by modeling errors. The applications on the detailed FA and the FLSNFC show that these methods have good applicability and accuracy. 126 Chapter 6 Efficient Bayesian inference for inverse analysis This chapter proposes an efficient Bayesian analysis for problems with high-dimensional input parameter space and high-dimensional output. The Bayesian method for inverse problems is first described, followed by an efficient algorithm adapted to the high-dimensional input and output. Then the classical MCMC and the block-update MCMC are introduced to increase the convergence of posterior sampling. Finally, these methods are applied to the FA and the FLSNFC to infer the internal damages given observed data of damaged models on the exterior surfaces of the models. 6.1 Bayesian method for inverse problems In mathematics, the Bayes’ formula states that p XjY (xjy) = p YjX (yjx)p X (x) p Y (y) (6.1) where X and Y are two different random variables; p XjY (xjy) is called the posterior probability of X given Y; p YjX (yjx) is usually referred to as the likelihood of X given fixed Y since p YjX (yjx) =L XjY (xjy) with L() the likelihood function;p X (x) is the marginal probability ofX and is usually called the prior probability of X; p Y (y) is the marginal probability of Y. The denominator in Eq (6.1) is usually not needed to explore the state space of the posterior, and hence p XjY (xjy)/p YjX (yjx)p X (x): (6.2) 127 In practice, X is usually used to denote the unknown random variable, and Y denotes another random variable where some observation (or data)y is given. Then the primary objective of the Bayes’ formula is to infer the posterior probability of some unknown quantities given some available data of related quantities. This is what an inverse problem is formulated in applications, especially when there are errors in the observed data, whichisinevitableinmanyapplications, andthustheBayes’formulawillbeusedtosolveourinterested inverse problems. Supposethataninverseproblemofthefollowingisinterested: calibratingthemodelM(!;)withrespect to parameters where!2R n is the set of interested frequencies, given an observation y as a solution of the model. Since the observation is frequently subjected to noise, the relationship between observation and physical model can be expressed as: y =M(!;) + (6.3) where denotes the noise. In Eq (6.3), assumey2R n ,2R n , and2R d . Then the probability ofy given has the form p(yj) =p (yM(!;)); (6.4) where p is the joint density of . This probability is the likelihood, the degree of belief in y given the parameters. The most common assumption for the observation error is that its components are assumed to be independent and identically distributed (i.i.d) random variables with i N (0; 2 ); for i = 1;:::;n: (6.5) Under this assumption, the likelihood can be written as, p(yj) = n Y i=1 p (y i M(!;) i ) = 1 (2) n=2 n exp [yM(!;)] T [yM(!;)] 2 2 : (6.6) Based on Bayes’ formula (6.2), the posterior probability p(jy) has p(jy)/p(yj)p (): (6.7) 128 Suppose the set of parameters has a probability density of p (), which is the prior probability of . Assume further that the input parameters i ’s are independent beta distributions with the same shape parameters and , which are the parameter distributions interested in the current study, and each with different scale parameter s i for i = 1;:::;d. Then p( i ) = i si 1 1 i si 1 ( +) ()()s i ; (6.8) where () denotes the gamma function. With this prior distribution and the observation errors being standard Gaussian, we can write the posterior probability as p(jy)/ 1 (2) n=2 n exp (yM(!;)) T (yM(!;)) 2 2 d Y i=1 8 > < > : i si 1 1 i si 1 ( +) ()()s i 9 > = > ; : (6.9) The samples of the high-dimensional posterior distribution can be generated by methods like Markov Chain Monte Carlo (MCMC) [18, 19, 20]; however, this process usually requires numerous evaluations of the physical model at various input parameters, which will quickly become prohibited for problems with the complex forward model. Another challenge of the posterior sampling is that exploring the high-dimensional parameter spaceR d that resides in is demanding since the acceptance rate of the posterior sample in the MCMC method will be small. In the sequel, these two major challenges will be addressed, in which, before introducing the MCMC method, an efficient formulation of Bayesian inference is first proposed in the next section. 6.2 Efficient formulation in Bayesian inference The efficient probabilistic model proposed in the previous chapter approximates the physical model with acceptable accuracy. With that model, each input of can be first mapped to a vector of standard Gaussian variables through i = 1 (F i ( i )) =T 1 i ( i ); for i = 1;:::;d; (6.10) 129 where F i () and are the cumulative distribution functions of i and the standard Gaussian variable, respectively. This mapping can be written in vector form as, =T 1 (): (6.11) Next, the surrogate model can be evaluated asH(!;). From chapter 5, we know that the probabilistic surrogateH(!;) builds a mapping from to the FRFH(!;). On the other hand, during the construction of the probabilistic model,M(!;)H(!;()) with() =T 1 () is evaluated for some samples of in the Galerkin projection in order to find the PCE coefficients. Thus, the operatorT 1 : 7! is one to one, and the operatorH : 7!H(!;) is one to one, which yields a one-to-one mapping from toH(!;T 1 ()), that isH :7!H(!;T 1 ()). Note thatH(!;T 1 ()) is a polynomial approximation of the physical model M(!;) and can be evaluated with negligible cost. The idea of efficient Bayesian inference method is to replaceM(!;) with H(!;T 1 ()) in the likelihood calculation of p(yj) =p (yM(!;)) as, p(yj) =p (yH(!;T 1 ())): (6.12) Then, for any physical parameter, the surrogate model of the FRF is evaluated to compute its associated likelihood of p(yj). Compared to the physical model computed by the two-level nested CB method, the surrogate model is more efficient by several computational orders. With the assumptions on observation error and prior probability, and Eq (6.12), we can get the posterior probability as p(jy)/ 1 (2) n=2 n exp [yH(!;T 1 ())] T [yH(!;T 1 ())] 2 2 d Y i=1 8 > < > : i si 1 1 i si 1 ( +) ()()s i 9 > = > ; : (6.13) In the above formula, both the likelihood and prior probabilities can be computed efficiently, and thus the MCMC method can be applied to generate samples of the posterior. 130 For the surrogate model of the FLSNFC, although the dimensions reduction techniques are applied to both the QoI and the input space, the reduced outputs (the limited number of KL parameters) are mapped back to the full-dimensional QoI (the FRF), and the reduced input parameters (the adapted variables) are mapped back to the full-dimensional input parameters through inverse transform sampling. Therefore, the surrogate model builds a map from the full-dimensional input parameter to the FRF, and the posterior distribution of all the input parameters can be inferred from this formula. The capability is crucial since the damages could happen anywhere inside the model, and these damages are all critical for decision-making. However, the requirement of maintaining the full-dimensional input space poses another significant challenge for posterior sampling. 6.3 Markov chain Monte Carlo simulation 6.3.1 Classical Metropolis-Hasting algorithm Generating random samples from high-dimensional joint distributions is not an easy task. MCMC methods [18,19,20]areproposedforthispurpose, whereaMarkovchainisconstructedsuchthatitwillasymptotically reach a stationary distribution that matches the desired distribution. In statistics, Metropolis-Hasting (MH) algorithm [18, 19, 20] is one of the most popular MCMC methods to produce random samples of some probability distribution in which direct sampling is difficult. Suppose that we are interested in generating samples from a probability distribution of p . In the MH algorithm, we first assume that a functionq(j) is easy to simulate, it is a transition probability transitioning from state to another state , and when a symmetric transition is applied, q( j) = q(j ). The transition probability uniquely defines the underlying Markov chain, and we call q( () j (i) ) the proposal distribution, which is the probability of transitioning from the current state (i) to a new state () . The proposal function is usually user-specified, and the Markov chain is defined once the proposal function is defined. Then, another important step of the MH algorithm is defining the acceptance distribution, 131 which is the probability of accepting the proposed state () . One common choice of acceptance probability A( () ; (i) ) in MH algorithm is A( () ; (i) ) = min ( 1; p () q (i) j () p (i) q () j (i) ) : (6.14) The definition ensures that a point located in a higher density region of p() has a higher acceptance probability. Withthesetwoessentialprobabilitiesdefined, theMHalgorithmissummarizedinAlgorithm6.1. In this algorithm, N MCMC is the number of steps in the Markov chain, u is a random sample drawn from a uniform distributionU(0; 1). Algorithm 6.1: Metropolis–Hastings (MH) algorithm to draw samples from a high-dimensional distribution 1 Initialize (0) ; 2 for i=0 to (N MCMC 1) do 3 Draw sample uU(0; 1) ; 4 Draw sample () q () j (i) ; 5 if u

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Algorithms for stochastic Galerkin projections: solvers, basis adaptation and multiscale modeling and reduction

PDF

Model, identification & analysis of complex stochastic systems: applications in stochastic partial differential equations and multiscale mechanics

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

PDF

Data-driven multi-fidelity modeling for physical systems

PDF

Comprehensive uncertainty quantification in composites manufacturing processes

PDF

A polynomial chaos formalism for uncertainty budget assessment

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Risk assessment, intrinsic interpolation and computationally efficient models for systems under uncertainty

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Impurity gas detection for spent nuclear fuel (SNF) canisters using ultrasonic sensing and deep learning

PDF

Multi-scale dynamic capture for high quality digital humans

PDF

Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms

PDF

Studies into computational intelligence approaches for the identification of complex nonlinear systems

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Plant substructuring and real-time simulation using model reduction

PDF

Efficient deep learning for inverse problems in scientific and medical imaging

PDF

Computational validation of stochastic programming models and applications

PDF

Modeling and dynamic analysis of coupled structure-moving subsystem problem

PDF

Multi-robot strategies for adaptive sampling with autonomous underwater vehicles

PDF

Structural nonlinear control strategies to provide life safety and serviceability

Asset Metadata

Creator Zeng, Xiaoshu (author)

Core Title Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Civil Engineering

Degree Conferral Date 2022-05

Publication Date 03/07/2022

Defense Date 12/15/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag basis adaptation in polynomial chaos expansion,dynamic reduction,large-scale model,multi-component systems.,OAI-PMH Harvest,probabilistic inverse analysis,stochastic reduction,substructuring

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ghanem, Roger Georges (committee chair), Gencturk, Bora (committee member), Johnson, Erik A. (committee member), Masri, Sami F. (committee member), Oberai, Assad A. (committee member)

Creator Email xiaoshuz@usc.edu,xiaoshuzen@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC110768086

Unique identifier UC110768086

Legacy Identifier etd-ZengXiaosh-10431

Document Type Dissertation

Format application/pdf (imt)

Rights Zeng, Xiaoshu

Type texts

Source 20220308-usctheses-batch-915 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu