Parametric Model Order Reduction of Guided Ultrasonic Wave Propagation in Fiber Metal Laminates with Damage

Abstract

This paper focuses on parametric model order reduction (PMOR) of guided ultrasonic wave propagation and its interaction with damage in a fiber metal laminate (FML). Structural health monitoring in FML seeks to detect, localize and characterize the damage with high accuracy and minimal use of sensors. This can be achieved by the inverse problem analysis approach, which employs the signal measurement data recorded by the embedded sensors in the structure. The inverse analysis requires us to solve the forward simulation of the underlying system several thousand times. These simulations are often exorbitantly expensive and trigger the need for improving their computational efficiency. A PMOR approach hinged on the proper orthogonal decomposition method is presented in this paper. An adaptive parameter sampling technique is established with the aid of a surrogate model to efficiently update the reduced-order basis in a greedy fashion. A numerical experiment is conducted to illustrate the parametric training of the reduced-order model. The results show that the reduced-order solution based on the PMOR approach is accurately complying with that of the high fidelity solution.

1. Introduction

Fiber metal laminates (FMLs) are hybrid materials that combine the ductile properties of metal with the high-specific stiffness of fiber-reinforced plastics (FRP). They possess the potential to offer a significant weight reduction and excellent fatigue properties in structural applications. These advantages led to a surge of interest in such materials in several industries, including aerospace and automotive, in the past decade. The complex structure of FMLs and their extensive application activated unique challenges not only for post-manufacture certification, but also for in-service inspection. In aerospace engineering, impacts with ground support equipment are deemed to be the major cause of in-service damage to layered laminates [1]. Collision with birds, runway stones, and ballistic impacts are some of the other sources of impact. Common damage types in composites are matrix cracking, fiber-matrix debonding, delamination, and fiber fracture [2,3]. Of all the mentioned damages, delamination is the most dangerous type of impact damage because there is no visible indication of damage on the outer surface. However, it can lead to an extremely detrimental effect on the stiffness and strength of the structure causing extensive loss of life and property. Therefore, a reliable nondestructive evaluation (NDE) method is required to accurately detect and characterize such damage, preferably in real-time using a minimum number of sensors. This gave rise to a new stream of NDE known as structural health monitoring (SHM). Several SHM techniques based on vibration, acoustics, and instrumentation have been developed over the years to detect damage at an early stage in composite materials [4]. Guided ultrasonic waves (GUW) feature (a) the ability to inspect large structures over the entire cross-sectional area; (b) the lack of need for complicated and expensive devices; (c) great cost-effectiveness [5]. Vibration-based SHM with GUW is widely acknowledged as one of the prominent techniques for damage identification in composites, as their propagation behavior changes when interacting with a damage [6,7].

Numerical tools, typically finite element method (FEM), are utilized to efficiently investigate the wave propagation and assess the suitability of GUW for damage detection. The solution to the inverse setting that combines the forward model and measurement data describes the location and class of the damage in the structure. As the forward simulations are performed at high frequencies of GUW, fine spatial and temporal discretization of the model is required. This subsequently results in a large number of degrees of freedom to be solved. The computational effort rises drastically when such models are employed as forward simulators in inverse problem analysis and uncertainty quantification, where the forward model must be solved several times. Model order reduction (MOR) aims to diminish the computational effort while the accuracy of the result remains the same. In most cases, reduced-order models are created based on one parameter value and then used to represent the system at different parameters in the parametric domain. However, the application of such reduced-order models at off-design often leads to poor predictions and failure. Therefore, every time parameters alter, modeling, solving and reduction of the system should be repeated, which makes MOR futile. It is therefore imperative to develop methods that preserve the parametric dependence of the system during its order reduction. Parametric model order reduction (PMOR) methods precisely satisfy this objective. Once the model reduction is accomplished with a PMOR method, design optimization and parametric space exploration cycles can be swiftly achieved by conducting parametric simulations in the trained reduced space itself.

Ettefagh et al. [8] used model updating based on model reduction by the Guyan method for damage-detection deploying modal analysis in aluminum beams, which showed higher accuracy and faster convergence. Fairly reasonable results were obtained by the laminate element method simulation of GUW in notched plate-like isotropic structures proposed by Glushkov et al. [9]. Based on heat transfer behavior, Liang et al. [10] used a proper orthogonal decomposition (POD) based model reduction strategy to detect cracks in a rectangular plate made of aluminum and iron. Bova et al. [11] compared the solutions of POD and modal-singular value decomposition (SVD) methods based on an offline–online approach in composite structures parameterized by fiber orientation and lamina thickness. Reduced-order modeling using POD was performed on a flexible aluminum plate by Capellari et al. [12]. They further used hybrid-Kalman filtering for damage detection. However, there was no previous work that studied PMOR on FMLs.

Considering the existing studies regarding the model reduction of lamb wave propagation, this research focuses on generating a low-cost but accurate model that simulates GUW propagation and its interaction with damage in FML for different values of the system parameters. The proposed method in this paper is concerned with the order reduction of high-dimensional nonlinear ordinary differential equation systems. The rest of this paper is organized as follows. Section 2 describes the targeted high dimensional finite element model developed in COMSOL-Multiphysics^® FE software. In Section 3 and Section 4, the standard MOR method and the PMOR approach used in this project are presented. The results from the numerical experiments are reported and discussed in Section 5. Finally, Section 6 summarizes the findings of this benchmarking study and discusses the areas of future work.

2. Finite Element Analysis of Lamb Wave Propagation

The numerical modeling and analysis of anti-symmetric lamb wave mode ($A_0$) propagation as well as its interaction with the defect in FML were studied using the FEM on COMSOL-Multiphysics ^® software. A 16-layered two dimensional carbon fiber reinforced plastic (CFRP)-steel laminate model was considered for further analysis. The simulations were carried out for 120 kHz central frequency of excitation. The excitation was realized by a five-cycle Hanning window sinusoidal burst, that was applied on the top and bottom left node of the model as shown in Figure 1.

The transient force applied upon excitation is as depicted in Figure 2. The simulation was conducted for a total time of $2.083\times {10}^{-4}$ s with the excitation pulse longing for $4.167\times {10}^{-5}$ s. A defect was introduced in the structure by a local reduction of the Young’s modulus within the steel lamina. As a result, a reduced Young’s modulus was assigned to a limited number of elements.

The boundaries of FML were subjected to free boundary conditions. All the nodes on the left end of the laminate were subjected to symmetry boundary condition, while the center node on the right end of the laminate was fixed in the x-direction. The considered thicknesses of steel and CFRP laminas and therefore, of the numerical model are 0.12 mm and 0.125 mm, respectively. The length and thickness of the laminate are 300 mm and 1.98 mm, respectively, with a defect placed at 90 mm from the left end of the laminate. The damage is sized to be as thick as the steel lamina and 4 mm in length. A virtual sensor placed at 57.5 mm records the out-of-plane displacements during the simulation. It was assumed that plane strain conditions would be prevalent. Albeit COMSOL implicit solver was used, time step and mesh sizes were calculated using a Courant–Friedrichs–Lewy condition. Following this, the time step size evaluated and used for the numerical modeling was $4.17\times {10}^{-7}$ s. The two-dimensional quadrilateral plane strain elements with quadratic Lagrange shape functions are used to discretize both steel and CFRP composite layers as well as the defect introduced. To ensure a proper representation of the wave propagation the element size is selected in accordance with [13,14]. This leads to one element over the thickness of each layer and an element length of 0.5 mm in the direction of the wave propagation. Hence, the numerical model consists of 79,266 nodal degrees of freedom in total. Finally, the material data that was used to define the steel and CFRP composite material models are listed in

Table 1. Mechanical properties of steel.

Density	Young’s Modulus	Poisson’s Ratio
(kg/m${}^3$)	(GPa)
7900	180	$0.3$

and 

Table 2. Mechanical properties of CFRP.

	Young’s Modulus
	(GPa)
$E_{11}=127.28$	$E_{22}=9.24$	$E_{33}=9.24$
	Poisson’s Ratio
	(- )
${\nu }_{12}=0.30$	${\nu }_{13}=0.30$	${\nu }_{23}=0.37$
	Shear Modulus
	(GPa)
$G_{12}=4.83$	$G_{13}=4.83$	$G_{23}=3.37$

, respectively. The density of CFRP was 1580 kg/m${}^3$.

Figure 3 provides different plots of the simulated wave propagation. The displacement field of the $A_0$ mode over the cross-section at the time instant $t=5\times {10}^{-5}$ s is provided in Figure 3a. For a better visualization, the displacement is scaled by a factor of $1\times {10}^{13}$. In Figure 3b, the finite element simulated GUW signal measured by the virtual sensor embedded in the FML is plotted. The reflected $A_0$ signal from the defect has its peak based on the position and worseness of the defect. Throughout this paper, it was considered that only one defect was present in the FML and its height was the same as that of the steel lamina. It was also assumed that the defect crops up only in the steel lamina. The simulation time for a solve using this high-fidelity (HiFi) numerical model was 66.29 s.

3. Model Order Reduction by Proper Orthogonal Decomposition

Model order reduction (MOR) pursues reducing the computational time and complexity of extensive dynamical systems by approximating them to the lower dimensions without losing the response characteristics of the system. The advantage of dealing with the basis functions is that it reduces the model order significantly. There exist several techniques that reduce the dimension of the system. They are classified into two main groups viz. Krylov and singular value decomposition (SVD) methods [15]. The former methods rely on iterative approaches for finding approximations to large-scale dynamical systems, whereas, the latter methods are based on the notion of expanding the Hi-Fi solution into a sequence of orthonormal basis functions, describing the most significant features of spatial and temporal variation of the field variable. One of the most prominent techniques in SVD methods is the POD method. The existence of POD dates back to the work of Karl Pearson [16] in 1901. Since then, the POD method has been successfully applied to numerous fields, namely, signal analysis and pattern recognition [17,18], statistics, geophysical fluid dynamics [19,20,21] or meteorology. Based on the domain of its application POD is known under several names: principal component analysis (in statistics) and Karhunen–Loève expansion (in stochastics). The works by Liang and Wu et. al [22,23] provide detailed descriptions of POD theory and its above-mentioned methods. For nonlinear systems, POD combined with Galerkin projection is widely used to produce reduced-order models. Consider an undamped linear dynamic mechanical high-fidelity (Hi-Fi) model following the system of equations:

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{Ku}=\mathbf{f}$$

(1)

with $\mathbf{M}\left(\theta \right)\in {\mathbb{R}}^{N\times N}$ denoting the global mass matrix, $\mathbf{K}\left(\theta \right)\in {\mathbb{R}}^{N\times N}$ representing the global stiffness matrix, $\mathbf{f}\left(t\right)\in {\mathbb{R}}^{N\times m}$ denoting the force of excitation, and $\mathbf{u}(\theta ,t)\in {\mathbb{R}}^{N\times m}$ the displacement matrix. $\theta \in \mathcal{D}\subset {\mathbb{R}}^{N_p}$ is a vector of $N_p$ parameters defining the system from the parametric space $\mathcal{D}$ and $t\in [0,t_{max}]$ is the time variable. The total number of degrees of freedom involved in the system is given by N and the number of discretized time steps is indicated by m. The order of this Hi-Fi model can be reduced by projecting the basis functions from the high-dimensional space ${\mathbb{R}}^N$ to a lower dimensional space ${\mathbb{R}}^n$ using a projection matrix $\mathsf{\Phi }(\theta ,t)\in {\mathbb{R}}^{N\times n}$ where $n<<m$:

$$\mathbf{u}\approx \mathsf{\Phi }{\mathbf{u}}_{\mathbf{r}},\ddot{\mathbf{u}}\approx \mathsf{\Phi }{\ddot{\mathbf{u}}}_{\mathbf{r}}.$$

(2)

Applying the above set of equations into Equation (1) and projecting them onto the basis, results in:

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}^{\mathbf{T}}\mathbf{M}\mathsf{\Phi }{\ddot{\mathbf{u}}}_{\mathbf{r}}+{\mathsf{\Phi }}^{\mathbf{T}}\mathbf{K}\mathsf{\Phi }{\mathbf{u}}_{\mathbf{r}}& ={\mathsf{\Phi }}^{\mathbf{T}}{\mathbf{f}}_{\mathbf{h}}\hfill \\ \hfill {\mathbf{M}}_{\mathbf{r}}{\ddot{\mathbf{u}}}_{\mathbf{r}}+{\mathbf{K}}_{\mathbf{r}}{\mathbf{u}}_{\mathbf{r}}& ={\mathbf{f}}_{\mathbf{r}}\hfill \end{array}$$

(3)

where ${\mathbf{M}}_{\mathbf{r}}$, ${\mathbf{K}}_{\mathbf{r}}\in {\mathbb{R}}^{n\times n}$, and ${\mathbf{f}}_{\mathbf{r}}\in {\mathbb{R}}^{n\times m}$ are the reduced system matrices. The approximated load matrix evaluated using the reduced-order solution is represented by ${\mathbf{f}}_{\mathbf{h}}\in {\mathbb{R}}^{N\times m}$. The projection matrix $\mathsf{\Phi }$, which contains the main features of the system, can be extracted using the POD method. To start with, a sufficient number of observations from the Hi-Fi model was collected in a matrix called snapshot matrix. The high-dimensional model can be analytical expressions, a finely discretized finite difference or a finite element model representing the underlying system. In the current case, the snapshot matrix $\mathbf{S}(\theta ,t)\in {\mathbb{R}}^{N\times m}$ was extracted and is further decomposed by thin SVD as follows:

$$\mathbf{S}=[{\mathbf{u}}_1,{\mathbf{u}}_2,\dots ,{\mathbf{u}}_m]$$

(4)

$$\mathbf{S}=\mathbf{P}\Sigma {\mathbf{V}}^{\mathbf{T}}.$$

(5)

In (5), $\mathbf{P}(\theta ,t)=[{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_m]\in {\mathbb{R}}^{N\times m}$ is the left-singular matrix containing orthogonal basis vectors, which are called proper orthogonal modes (POMs) of the system, $\Sigma =diag({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_m)\in {\mathbb{R}}^{m\times m}$, with ${\sigma }_1\ge {\sigma }_2\ge \dots \ge {\sigma }_m>0$, denotes the diagonal matrix containing the singular values ${\left\{{\sigma }_k\right\}}_{k=1}^m$ and $\mathbf{V}\in {\mathbb{R}}^{m\times m}$ represents the right-singular matrix, which will not be of much use in this method of MOR. In general, the number of modes n required to construct the data is significantly less than the total number of modes m available. In order to decide the number of most influential mode shapes of the system, a relative energy measure E described as follows is considered:

$$E={\displaystyle \frac{{\sum }_{k=1}^n{\sigma }_k}{{\sum }_{k=1}^m{\sigma }_k}}.$$

(6)

The error $\epsilon$ from approximating the snapshots using POD basis can then be obtained by:

$$\epsilon =\frac{{\sum }_{k=n+1}^m{\sigma }_k}{{\sum }_{k=1}^m{\sigma }_k}.$$

(7)

Based on the preferred accuracy, one can select the number of POMs required to capture the dynamics of the system. The collection of POMs results in the projection matrix

$$\mathsf{\Phi }=[{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_n]\in {\mathbb{R}}^{N\times n}.$$

(8)

Once the projection matrix is obtained, the reduced system (3) can be solved for $u_r\mathrm{and}{\ddot{u}}_r$. Subsequently, the solution for the full order system can be evaluated using (2). The approximation of high-dimensional space of the system largely depends on the choice of extracting observations to ensemble them into the snapshot matrix. For a detailed explanation on the POD basis in general Hilbert space, the reader is directed to the work of Kunisch et al. [24].

4. Parametric Model Order Reduction

4.1. Overview

The reduced-order models produced by the method described in Section 3 usually lack robustness concerning parameter changes and hence must often be rebuilt for each parameter variation. In real-time operation, their construction needs to be fast such that the precomputed reduced model can be adapted to new sets of physical or modeling parameters. Most of the prominent PMOR methods require sampling the entire parametric domain and computing the Hi-Fi response at those sampled parameter sets. This avails the extraction of global POMs that accurately captures the behavior of the underlying system for any given parameter configuration. The accuracy of such reduced models depends on the parameters that are sampled from the domain. In POD-based PMOR, the parameter sampling is accomplished in a greedy fashion-an approach that takes a locally best solution hoping that it would lead to the global optimal solution [25,26,27,28,29,30,31,32,33,34,35,36,37]. It seeks to determine the configuration at which the reduced-order model yields the largest error, solves to obtain the Hi-Fi response for that configuration and subsequently updates the reduced-order model. Since the exact error associated with the reduced-order model cannot be computed without the Hi-Fi solution, an error estimate is used. Based on the type of underlying PDE several a posteriori error estimators [38,39,40,41,42], which are relevant to MOR, were developed in the past. Most of the estimators used in [26,27,28,29,30,31,32,33,36,37,38,39,40,41,43,44] focused on the affine parameter dependency of the Hi-Fi model, which resulted in an offline/online decomposition approach: expensive computations of lower-order matrices are done offline while the norm of the residual for any given parameter configuration was computed online with a minimal effort.

The POD-based global reduced-order models are well suited for approximating parametrized elliptic and parabolic PDE models. However, the PMOR of a wide range of hyperbolic problems with non-linearity and discontinuity still remain a challenge. Therefore, a strong investigation is going on in the MOR research community to reduce the Kolmogorov n-width of the solution manifold [45,46]. Projection-based MOR along with an active manifold was carried out by Boncoraglio et al. [32] to efficiently solve multidisciplinary design optimization problem, which blends both linear and nonlinear constraints in aerodynamics. The authors used a deep convolutional autoencoder to discover the pertinent active manifold for dimensionality reduction of a high-dimensional design parameter space (58 structural and shape parameters). Bui-Thanh et al. [36] achieved MOR for design optimization of a heat conduction fin by implementing an efficient adaptive algorithm to determine appropriate sample points over a large input parametric space (up to 21 design parameters). McBane et al. [43] proposed a component-wise reduced-order model based on [47,48] to optimize the topology of a lattice-type structure. They further went on to simplify the model to increase the speedup of the optimization process. A space-time MOR method built on least-squares Petrov–Galerkin projection was presented by Kim et al. [44] to solve linear dynamical systems. The approach was well demonstrated on 2D diffusion and 2D convection diffusion problems. Further contributions on PMOR span across the domains of contact in multibody nonlinear dynamics [49], nonlinear fluid–structure interaction problems [50], uncertainty quantification [51,52] and contact mechanics [53,54].

Paul-Dubois-Taine et al. [35] employed an alternative approach, built on the notion of optimization techniques, that samples the parameters adaptively and extracts the efficient global POMs. A surrogate model for the evaluated a posteriori error indicators was constructed, which enabled the localization of parametric domain where the probability of error is the largest. This facilitated an efficient training method to produce an accurate reduced-order model for the underlying Hi-Fi model with complex parameter dependencies. The authors have illustrated the effective functioning of the method on linear and nonlinear mechanical systems. Considering only the linear dynamical system and lower dimensional input parametric space, in this research work, the method presented by Paul-Dubois-Taine et al. was adopted to accomplish PMOR for GUW propagation in a defective FML.

4.2. An Adaptive POD-Greedy Approach

A standard POD-greedy method progresses by finding a parameter ${\theta }_i$, at every iteration i, that maximizes the norm of the error $\mathbf{e}\left(\theta \right)$ between the reduced-order solution and its underlying Hi-Fi solution defined as follows:

$${\parallel \mathbf{e}\left(\theta \right)\parallel }_2={\int }_0^{t_{max}}{\parallel \mathbf{u}(\theta ,t)-\mathsf{\Phi }{\mathbf{u}}_{\mathbf{r}}(\theta ,t)\parallel }_{2}dt.$$

(9)

The Hi-Fi model was then solved for parameter ${\theta }_i$ to extract the basis vector corresponding to ${\theta }_i$ and update the projection matrix $\mathsf{\Phi }$. As the Hi-Fi solution $\mathbf{u}(\theta ,t)$ in practice was unknown prior to solving the Hi-Fi model, an error indicator was utilized in lieu of the error norm. Besides several existing error estimators, the norm of the residual is a well-known alternative to the error estimates:

$${\parallel \mathbf{r}\left(\theta \right)\parallel }_2={\int }_0^{t_{max}}{\parallel \mathbf{f}\left(t\right)-{\mathbf{f}}_{\mathbf{h}}\left(t\right)\parallel }_{2}dt.$$

(10)

Throughout this paper, the a posteriori error indicator $\mathcal{J}\left(\theta \right)$ is referred to the norm of the residual ${\parallel \mathbf{r}\left(\theta \right)\parallel }_2$. In the classical POD-greedy approach, a finite set of candidate parameters $\zeta$ of cardinality $N_{\zeta }$ is searched iteratively to identify the parameter ${\theta }_i$ that yields the largest error norm. When the dimension $N_p$ of the parameter space $\mathcal{D}$ is large and if the number of randomly collected candidate parameter sets is small, it is likely that the target parameter configuration is not included. This challenge is dealt with by determining $N_{\zeta }$ set of candidate parameters at every iteration i in an adaptive manner following a greedy algorithm. It was illustrated by the works of [35,55] that the adaptive PMOR approach requires limited offline training time relative to that of the classical PMOR method. The objective of adaptive parameter sampling is to seek the optimal parameter ${\theta }_i$, in every iteration i, from the pool of error indicators evaluated over sets of candidate parameters of smaller cardinality. This procedure is initiated by selecting a parameter point from $\mathcal{D}$ and its associated reduced-order basis ${\mathsf{\Phi }}_1\in {\mathbb{R}}^{N\times 1}$ is computed. Next, the first set $q=0$ of candidate parameter points ${\zeta }_{i,0}\subset \mathcal{D}$ of smaller cardinality $N_{\zeta }^0<N_{\zeta }$ are randomly selected. For each of these points, the algorithm evaluates the reduced-order model and also their corresponding residual-based error indicators ${\left\{{\mathcal{J}}_j\left(\theta \right)\right\}}_{j=1}^{N_{\zeta }^0}$. These error indicators are then used to build a surrogate model ${\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)$ for the error estimator over the entire parametric domain $\mathcal{D}$. In this work, a multiple linear regression-based surrogate model is used. Subsequently, the designed surrogate model is employed to estimate the location of an additional set $q=1$ of candidate parameters ${\zeta }_{i,1}\subset \mathcal{D}$ with high probability to have largest error estimates. The cardinality of the newly added set is $N_{\zeta }^1<N_{\zeta }$. Once the surrogate model was constructed, the probability of candidate points neighboring the highest error indicator was evaluated by the following method proposed in [56]. This involves computing the maximum value ${\widehat{\mathcal{J}}}_{max}^{\left[q\right]}\left(\theta \right)$ of the surrogate model ${\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)$ over $\mathcal{D}$ and then selecting a series of targets ${\mathcal{T}}_j>{\widehat{\mathcal{J}}}_{max}^{\left[q\right]},j=1,\dots ,N_{\mathcal{T}}$. The target values were chosen similar to those used in [56]. Together with the mean-squared error $s^{\left[q\right]}\left(\theta \right)$ of the surrogate model ${\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)$, the associated probability ${\mathsf{\Xi }}_{{\mathcal{T}}_j}\left(\theta \right)$ for each of those target values is modeled assuming a Gaussian distribution for $\mathcal{J}\left(\theta \right)$ with mean ${\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)$ and standard deviation $s^{\left[q\right]}\left(\theta \right)$ as:

$${\mathsf{\Xi }}_{{\mathcal{T}}_j}\left(\theta \right)=\kappa \left(\frac{{\mathcal{T}}_j-{\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)}{s^{\left[q\right]}\left(\theta \right)}\right)$$

(11)

where $\kappa \left(·\right)$ represents the normal cumulative distribution function (CDF). The point ${\theta }^{\star }\in \mathcal{D}$ that maximizes ${\mathsf{\Xi }}_{{\mathcal{T}}_j}$ is then selected. The set of ${\left\{{\theta }_j^{\star }\right\}}_{j=1}^{N_{\mathcal{T}}}$ is clustered by means of K-means clustering. The optimal number $N_{clust}$ of cluster points are evaluated with the aid of the "evalclusters" function built-in MATLAB 2019b. As a result, the parameters corresponding to the cluster centers are added as the additional set of candidate parameters. The algorithm then determines the reduced-order model for the additional candidate points and estimates their error indicators ${\left\{{\mathcal{J}}_l\left(\theta \right)\right\}}_{l=1}^{N_{\zeta }^1}$. This process is then repeated until the maximum cardinality $N_{\zeta }$ is reached with $q=N_{add}$ sets of candidate parameters, i.e., $N_{\zeta }=N_{\zeta }^0+N_{\zeta }^1+\dots +N_{\zeta }^{N_{add}}$. The pool of error indicators:

$$\mathcal{J}\left(\theta \right)={\left\{{\mathcal{J}}_j\left(\theta \right)\right\}}_{j=1}^{N_{\zeta }^0}\cup {\left\{{\mathcal{J}}_l\left(\theta \right)\right\}}_{l=1}^{N_{\zeta }^1}\cup \dots \cup {\left\{{\mathcal{J}}_r\left(\theta \right)\right\}}_{r=1}^{N_{\zeta }^{N_{add}}}$$

(12)

is used to adapt a new surrogate model. The accuracy of the surrogate model increases with increasing iterations as more error indicator values will be available. The step-by-step procedure for determining the additional candidate points is described in Algorithm 1. Furthermore, the candidate parameters with different cardinality are concatenated to form ${\zeta }_i$:

$${\zeta }_i={\zeta }_{i,0}\cup {\zeta }_{i,1}\cup \dots \cup {\zeta }_{i,N_{add}}.$$

(13)

Algorithm 1 Identification of the additional points with maximum probability of having the highest error indicator value.

Input: Surrogate model $\widehat{\mathcal{J}}\left(\theta \right)$ using multiple regression technique and its associated mean-square error $s^{\left[q\right]}\left(\theta \right)$ Output:$N_{clust}$ number of points with high probability of having largest error indicator 1: Evaluate maximum ${\widehat{\mathcal{J}}}^{max}$ and minimum ${\widehat{\mathcal{J}}}^{min}$ values of the surrogate model over $\mathcal{D}$ 2: Pick a large set of trial points $\mathcal{S}\subset \mathcal{D}$ 3: Define $N_{\mathcal{T}}$ number of target values $\{{\alpha }_1,\dots ,{\alpha }_{N_{\mathcal{T}}}\}$ as used in [56] 4: for j = $1,\dots ,N_{\mathcal{T}}$ do 5:   Compute ${\mathcal{T}}_j={\widehat{\mathcal{J}}}^{max}+{\alpha }_j({\widehat{\mathcal{J}}}^{max}-{\widehat{\mathcal{J}}}^{min})$ 6:   Identify the points with maximum probability for each of the target values $${\theta }_j^{\star }=\underset{\theta \in \mathcal{S}}{arg\; max}\left({\mathsf{\Xi }}_{{\mathcal{T}}_j}\left(\theta \right)\right)=\underset{\theta \in \mathcal{S}}{arg\; max}\left(\kappa \left(\frac{{\mathcal{T}}_j-\widehat{\mathcal{J}}\left(\theta \right)}{s^{\left[q\right]}\left(\theta \right)}\right)\right)$$ 7: end for 8: Evaluate the optimal number $N_{clust}$ of clusters 9: Cluster the set {${\theta }_1^{\star }$,…,${\theta }_{N_{\mathcal{T}}}^{\star }$} by means of K-means clustering algorithm into $N_{clust}$ number of clusters 10: Compute their associated $N_{clust}$ cluster centers

Ultimately, the optimal parameter ${\theta }_i$ is identified as the maximizer of the error indicators collected over the set ${\zeta }_i$ as follows:

$${\theta }_i=\underset{\theta \in {\zeta }_i}{arg\; max}\mathcal{J}\left(\theta \right)$$

(14)

and is further used to update the parameter independent projection subspace $\mathsf{\Phi }$. It should be noted that ${\theta }_i$, identified by the above-mentioned process, need not be accurate as the parameter sampling is performed roughly around the maximum error indicator. Besides that, the relative error between the Hi-Fi model and reduced model at ${\theta }_i$ both before and after the reduced-order basis update, ${\mathbf{e}}^b\left(\theta \right)$ and ${\mathbf{e}}^a\left(\theta \right)$, respectively, were also evaluated. As a result, at every greedy iteration i, the following set of relative error and error indicators are available to construct an approximate error model:

$$\eta =\{({\mathcal{J}}_1^b,{\mathbf{e}}_1^b),({\mathcal{J}}_1^a,{\mathbf{e}}_1^a),\dots ,({\mathcal{J}}_i^b,{\mathbf{e}}_i^b),({\mathcal{J}}_i^a,{\mathbf{e}}_i^a)\}.$$

(15)

In this paper, a linear model was built to compute the relative error value for an estimated error indicator value. The error estimator becomes much more precise as the number of iterations i increases in the greedy algorithm. The training phase or the update of reduced-order basis in this approach continues until either of the conditions is satisfied: (a) the error estimate ${\mathbf{e}}_i$ equals its maximum threshold ${\mathbf{e}}_{max}$ and (b) the maximum number of basis updates are attained $i=I_{max}$. The complete adaptive POD-greedy approach is summarized in Algorithm 2.

Algorithm 2 Adaptive POD-Greedy algorithm.

Input: Maximum number of iterations $I_{max}$, initial cardinality of set of candidate parameters $N_{\zeta }^0$, parametric space $\mathcal{D}$, maximum error tolerance ${\mathbf{e}}_{max}$, Output: Global reduced-order basis $\mathsf{\Phi }$ 1: Randomly select an initial parameter point from the parameter domain ${\theta }_1\in \mathcal{D}$ 2: Compute the Hi-Fi solution ${\left\{\mathbf{u}({\theta }_1,t)\right\}}_{t=0}^{t_{max}}$ for that parameter 3: Generate reduced-order basis ${\mathsf{\Phi }}_1\in {\mathbb{R}}^{N\times 1}$ by POD using the Hi-Fi solution 4: Initiate an empty error pool $\eta =⌀$ 5: for$i=2,\dots ,I_{max}$do 6:       Randomly select an initial set of $N_{\zeta }^0$ candidates $\zeta =\{{\theta }^{\left(1\right)},\dots ,{\theta }^{\left(N_{\zeta }^0\right)}\}\subset \mathcal{D}$ 7:       for $j=1,\dots ,N_{\zeta }^0$ do 8:          Determine the reduced-order solution ${\left\{\mathbf{w}\left({\theta }^{\left(j\right)}\right)\right\}}_{t=0}^{t_{max}}$ for the parameter ${\theta }^{\left(j\right)}$    using reduced-order basis ${\mathsf{\Phi }}_{i-1}$ 9:          Estimate the corresponding residual-based a posteriori error indicator $\mathcal{J}\left({\theta }^{\left(j\right)}\right)$ 10:     end for 11:     Set q = 0 12:     while $N_{\zeta }^q<N_{\zeta }$ do 13:          $q=q+1$ 14:          Fit a surrogate model ${\widehat{\mathcal{J}}}^{\left[q\right]}\left(\theta \right)$ for the error indicators and compute its mean-    squared error $s^{\left[q\right]}\left(\theta \right)$ 15:          Determine ${\zeta }_{i,q}$ candidates points with maximum probability of having largest    error indicator value using Algorithm 1 16:          Form a set $\zeta =\zeta \cup {\zeta }_{i,q}=\{{\theta }^{\left(1\right)},{\theta }^{\left(2\right)},\dots ,{\theta }^{\left(N_{\zeta }^q\right)}\}$ 17:          for j = $N_{\zeta }^{q-1}+1,\dots ,N_{\zeta }^q$ do 18:             Compute the reduced-order solution ${\left\{\mathbf{w}\left({\theta }^{\left(j\right)}\right)\right\}}_{t=0}^{t_{max}}$ for the parameter $\theta \left(j\right)$        using reduced-order basis ${\mathsf{\Phi }}_{i-1}$ 19:             Evaluate their associated a posteriori error indicator $\mathcal{J}\left({\theta }^{\left(j\right)}\right)$ 20:          end for 21:       end while 22:       Find ${\theta }_i={arg\; max}_{\theta \in \zeta }\mathcal{J}\left(\theta \right)$ 23:       if $(i>2)$ & $({\mathbf{e}}_i>{\mathbf{e}}_{max})$ 24:          Set $\mathsf{\Phi }={\mathsf{\Phi }}_{i-1}$ 25:          break 26:       end if 27:       Evaluate the Hi-Fi solution ${\left\{\mathbf{u}({\theta }_i,t)\right\}}_{t=0}^{t_{max}}$ associated with ${\theta }_i$ 28:       Compute the exact relative error ${\mathbf{e}}_i^b$ for ${\left\{\mathbf{w}\left({\theta }_i\right)\right\}}_{t=0}^{t_{max}}$ using ${\mathsf{\Phi }}_{i-1}$ and its associated error indicator ${\mathcal{J}}_i^b$ 29:       Update the reduced-order basis ${\mathsf{\Phi }}_i$ by POD using ${\left\{\mathbf{u}({\theta }_i,t)\right\}}_{t=0}^{t_{max}}$ 30:       Compute the exact relative error ${\mathbf{e}}_i^a$ for the reduced-order solution ${\left\{\mathbf{w}\left({\theta }_i\right)\right\}}_{t=0}^{t_{max}}$ using updated ${\mathsf{\Phi }}_i$ and its associated error indicator ${\mathcal{J}}_i^a$ 31:       Form a set of error pool $\eta =\eta \cup ({\mathcal{J}}_i^b,{\mathbf{e}}_i^b)\cup ({\mathcal{J}}_i^a,{\mathbf{e}}_i^a)$ 32:       Establish an error model that estimates error for a given error indicator value using the error pool $\eta$ 33:       if $i=I_{max}$ then 34:          Set $\mathsf{\Phi }={\mathsf{\Phi }}_i$ 35:       end if 36: end for

5. Application

The presented adaptive POD-greedy algorithm is applied to train the reduced-order model of Lamb wave propagation in an FML structure with damage. The Hi-Fi model of the system under consideration was developed based on FEM using COMSOL-Multiphysics^® commercial software. As described in Section 2, the damage was modeled as the loss of stiffness in a localized area within the steel lamina. The training was carried out on a 3D parametric space $\mathcal{D}$, defined by Young’s modulus, position on the x-axis and length of the damage. The parameter is then $\theta =\{E_{damage},x_{damage},l_{damage}\}\in {\mathbb{R}}^3$. The parametric space is delimited by $E_{damage}\in [5\times {10}^6,50\times {10}^8]$ Pa, $x_{damage}\in [20,150]$ mm and $l_{damage}\in [2,15]$ mm. To avoid long simulation time and unnecessary boundary reflections of GUW, the position of the damage on the x-axis is limited to 150 mm. At every greedy iteration, the global reduced-order basis is enriched with a better approximating subspace using the POD method.

Since the proposed PMOR approach is model intrusive, it is required to solve the system of dynamic equations using the prominent Newmark’s integration method (the algorithm is presented in Appendix A. for the completeness of this paper). The comparison of the out-of-plane displacement of FML obtained by the Hi-Fi model and reduced-order model, produced using adaptive PMOR approach with Galerkin projection, is shown in Figure 4. It can be observed that almost accurate predictions are obtained with the reduced-order model. The signal reflected by the damage in the laminate was perfectly captured by the trained low-dimensional model. In addition, a reduced-order solution was obtained using the same global basis for several parameter configurations in $\mathcal{D}$ to test its robustness (see Figure 5). In all the cases, the predictions were very good and the damage signals were reconstructed flawlessly. The global reduced-order basis accumulated 340 modes to attain this accuracy. At every iteration i, an arbitrary cardinality of $N_{\zeta }^0=10$ initial set of candidate parameters were randomly sampled from the parametric space. The set was then enriched with additional parameters in an adaptive procedure until $N_{\zeta }=20$ was reached.

The evolution of residual-based a posteriori error indicator and relative error between the Hi-Fi model and the low-dimensional model are shown in Figure 6a,b, respectively. The values of relative error as well as error indicator decay with increasing number of modes implying the fine approximation of global projection matrix.

When the cardinality of initial parameter sets was increased to $N_{\zeta }^0=15$, it was required to enrich the reduced-order basis with 354 modes to achieve accurate prediction. This indicates that the adaptive sampling of parameters is definitely efficacious compared to random sampling. The linear error estimator built using the collected error pool at every iteration is shown in Figure 7. It is evident from Figure 7a,b that the error estimator is refined with increasing greedy iterations.

The numerical experiment was performed on a 4-core Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz processor with 16 GB RAM. The computational cost for computing Hi-Fi solution and reduced-order solution is summarized in

Table 3. Computation time for high-dimensional and reduced-order models.

Model	Training Time	Computational Time
High-dimensional	-	66.29 s
Reduced-order	17.6 h	1.96 s

. Using the adaptive POD-greedy approach, a speedup factor of 33.82 is achieved. This substantial decrease in the computational effort is very much appreciated in inverse problem analysis for the localization and characterization of the defect in the FML.

The global modes obtained for the trained parametric domain was then tested to build the reduced-order model for a parameter lying off the trained space. The damage was introduced in the CFRP lamina, which was two layers away in z-axis from the previously considered steel lamina. The low-dimensional model was produced for the parameter configuration $\theta =(3.55\mathrm{GPa},90\mathrm{mm},7\mathrm{mm})$. Figure 8 shows that the reduced-order solution was able to very well detect the damage present in the FML. The location of the damage signal produced by the reduced-order model complies precisely with that of the Hi-Fi model Moreover, there was only a subtle difference in the magnitude of the damage signal. By extending the training parametric space to the parameter set $\theta =\{E_{damage},x_{damage},z_{damage},l_{damage}\}\in {\mathbb{R}}^4$, the low-dimensional solution can be very accurately obtained using the adaptive POD-greedy approach. However, the training should be executed from the beginning when the parameter exploration space is varied.

The accuracy of POD-based global reduced-order model, which produces the response for any point in the parametric space, depends on the size of the extracted global reduced-order basis. In some cases, due to the extensively large parametric space obtaining an efficient global reduced-order model is quixotic (as in the case of [32]). To overcome this limitation of global reduced-order model approach, Y. Choi et al. [31] presented a novel methodology that accelerates the solution of design optimization problem. The approach uses a database of local parameterized reduced-order models constructed in offline and interpolates those reduced-order models online to produce a new reduced-order model for an unsampled point in the parametric space queried during the optimization process. The accuracy of the resulting reduced-order model depends on the database produced in the offline phase. Y. Choi et al. performed an efficient database construction based on a saturation assumption greedy procedure proposed by Hesthaven et al. [57]. According to this greedy procedure, a saturation constant that indicates the nature of an error estimate for a parameter is evaluated (see Definition 1 in [31]). Consequently, the computation of error estimates at some points are judiciously avoided and the overall computation time was considerably reduced.

Nevertheless, the approach of adaptive PMOR using a surrogate model employed in this research work was also capable of producing an efficient global reduced-order model for high dimensional parameter space problems. Binder et al. [55] also adopted it to speed up the computation of a convection-diffusion-reaction PDE with parameter space of dimension up to ${\mathbb{R}}^{10000\times 11}$ that arises in analyzing financial risks. Thus, in this work, the application of adaptive PMOR approach for GUW propagation in a defective FML in relatively smaller parametric space was successfully demonstrated. The resulting global reduced-order model significantly reduced the computation time without compromising on the accuracy.

6. Conclusions

In this paper, a parametric model reduction method was employed to produce reduced-order models for a high-dimensional linear dynamical structural system with a speedup factor of 33.82. A finite element method has been utilized to solve the high-dimensional system. The global reduced-order basis produced by the presented adaptive POD-greedy approach is robust for any parameter configuration from the considered parametric domain. An adaptive sampling technique using a multiple linear regression-based surrogate model was exploited to locate the parameters that are most likely to maximize the error indicator. The modes corresponding to these parameters were accumulated in a greedy fashion and the global reduced-order bases are enriched until the necessary accuracy is achieved. The method was tested and studied on a numerical experiment of guided ultrasonic wave propagation in a damaged carbon fiber reinforced epoxy-steel laminate. The reduced-order model generated using the presented approach was able to predict the solution and detect the damage which was even as small as 2 mm in length very accurately. Moreover, it was also capable of capturing a detailed response of the system for parameters that are even marginally away from the trained parameter space. In the future, this research will continue to use this expeditious low-cost model for the inverse analysis to (a) localize and characterize the damage in the fiber metal laminate and (b) quantify the uncertainties concerning the damage.

In this work, a POD-based projection method was used to produce the global reduced-order model. However, such a linear subspace based approach will often not be able to approximate the solutions of hyperbolic systems, which involve discontinuity and non-linearity, with small number of modes. One has to make use of the local reduced-order models interpolation [31] or nonlinear manifold based [32] approaches. Another interesting aspect of our future work within this research project will be to address these issues in the context pertinent to the guided ultrasonic wave propagation in damaged fiber metal laminates.