Digital twin approach for intelligent operation planning and health management of mechanical systems

The digital twin paradigm aims to fuse information obtained from sensor data, physics models, and operational data for a mechanical component in use to make well-informed decisions regarding health management and operations of the component. In this work, we discuss a methodology for digital-twin-based operation planning in mechanical systems to enable: a) cost-effective maintenance scheduling, and b) resilient operations of the system. As properties of mechanical systems, as well as their operational parameters, loads and environment are stochastic in nature, our methodology includes probabilistic damage diagnosis, probabilistic damage prognosis, and system optimization under uncertainty. As an illustrative example, we consider the problem of fatigue crack growth in a metal component. We discuss a probabilistic, ultrasonic-guided-wave-based crack diagnosis framework that can handle both aleatory and epistemic uncertainties in the diagnosis process. We build a high-fidelity, finite element model to simulate the piezoelectric effect and ultrasonic guided wave propagation. We use test data obtained by conducting diagnostic experiments on the physical twin to calibrate the error in the diagnosis model. We perform Bayesian diagnosis of crack growth using the corrected diagnosis model, considering data corrupted by measurement noise, and fuse the information from multiple sensors. We build a finite-element-based high-fidelity model for crack growth under uniaxial cyclic loading, and calibrate a phenomenological (low-fidelity) fatigue crack growth model using the high-fidelity model output as well as data from fatigue loading experiments performed on the physical twin. We use the resulting multi-fidelity model in a probabilistic crack growth prognosis framework; thus achieving both accuracy and computational efficiency. Lastly, we utilize the damage diagnosis framework along with the damage prognosis model to optimize system operations under diagnostic and prognostic uncertainty. We perform simulation as well as laboratory experiments that show how the digital twin of the component of interest can be used for intelligent health management and operation planning for mechanical systems.


INTRODUCTION
The need for aerospace vehicles to operate for long periods of time without the opportunity for maintenance or repair makes strategies for extending the maintenance-free operation window important. These strategies include reconfiguration of the system or redesigning system operations to minimize the stress experienced by critical system components. The system reconfiguration strategy depends both on the health state and its diagnosis; and the prediction of how the damage will evolve given the task to be performed (damage prognosis). The digital twin paradigm (Glaessgen & Stargel, 2012;El Saddik, 2018;Lee, Bagheri, & Kao, 2015;Söderberg, Wärmefjord, Carlson, & Lindkvist, 2017;Yang, Shen, & Wang, 2018;Tao, Cheng, et al., 2018;Bruynseels, Santoni de Sio, & van den Hoven, 2018;Tao, Zhang, Liu, & Nee, 2019), which integrates the diagnostic and prognostic information with an optimization framework, is well-suited for redesigning aerospace system operations. In this article, we discuss the digitaltwin-based methodology and illustrate its utility by performing a simulation experiment. In the simulation experiment, both (physical and digital) twins are represented by computational models. As an example, we consider fatigue crack growth from a hole (a surrogate for initial damage) in a thin Aluminum plate under uni-axial, cyclic loading. The Aluminum plate is a surrogate for an aerostructural component, whereas the cyclic loading pattern applied to the plate represents the loading applied to the component during different vehicle maneuvers. Thus, the theory-critical problem considered in this work is a simplified version of designing damage-adaptive, resilience-enhancing aerospace vehicle maneuvers. We discuss probabilistic diagnosis using ultrasonic guided-wave pitch catch, probabilistic fatigue crack prognosis using linear-elastic-fracture-mechanics-based, subcritical fatigue crack growth models, and strategies for loadprofile optimization under uncertainty to restrict the fatigue crack growth below a threshold value while ensuring that the work done by the system is greater than a predefined threshold value.
Design of damage-adaptive system operations necessitate quantification of the current damage state. Both the severity of damage and the (aleatory and epistemic) uncertainty in the diagnosis need to be quantified. Aleatory uncertainty is the irreducible uncertainty associated with natural variability, whereas epistemic uncertainty arises due to lack of knowledge, and is reducible. We discuss damage diagnosis using the ultrasonic guided-wave pitch-catch method (Alleyne & Cawley, 1992;Chang & Mal, 1999;Michaels & Michaels, 2005;J.-B. Ihn & Chang, 2008;Michaels, 2008;Sbarufatti, Manson, & Worden, 2014;Yang et al., 2016;Janapati, Kopsaftopoulos, Li, Lee, & Chang, 2016;He, Ran, Liu, Yang, & Guan, 2017;Wang, He, Guan, Yang, & Zhang, 2018). The ultrasonic pitch-catch method utilizes a network of monolithic PZT (Pb(Zr − Ti)O 3 ) transducers bonded to the component of interest as actuators and sensors for the Lamb waves propagating in the plate. We use a physics-based damage index calculated using features of the time-varying electric potential measured by the sensors for damage diagnosis. We discuss a Bayesian framework to perform probabilistic diagnosis that enables treatment of aleatory and epistemic uncertainties involved in the process. We perform low-fidelity modeling of the governing physics, sensitivity analysis and dimension reduction, high-fidelity modeling of the governing physics to generate training data and build computationally inexpensive surrogate models to perform Bayesian estimation and information fusion for crack length prediction.
Given the damage state in a component of interest (along with an estimate of the diagnosis uncertainty), damage-adaptive system operations can be designed if the propagation of damage under various candidate system operations can be estimated. This can be achieved by means of probabilistic damage prognosis in aerospace structural components. Here, a hole at the center of the aluminum plate represents the initial damage, and the crack growth due to applied cyclic loads represent damage evolution. Thus, a fatigue crack growth model is needed to perform damage prognosis. We employ an analytical model based on the assumptions of linear elastic fracture mechanics (LEFM) with small-scale plasticity, where the sub-critical fatigue crack growth due to the applied cyclic loading is estimated by calculating the crack growth rate as a function of stress-intensity factors and other model parameters. We consider the uncertainty in the diagnosis of the initial crack size, as well as the aleatory uncertainty in the crack growth estimation. Once again, we use computationally expensive surrogate models of the underlying physics to ensure that sampling-based probabilistic analyses are performed in a computationally efficient manner.
If the damage severity (and the associated uncertainty) is known, and a well-calibrated model for probabilistic damage prognosis is available, then system reconfiguration to ensure extended maintenance-free operation can be performed. In this work, as an example, we consider the problem of (cyclic) load-profile optimization. We assume that a given aerospace vehicle maneuver is associated with a characteristic (cyclic) load level range, and we seek the optimal magnitude as well as the optimal duration of the load (intensity and duration of the maneuver) to ensure: a) maintenance-free operation period is extended beyond the specified duration (that is, the damage growth is below a specified threshold), b) minimum required work is performed. We assume that the component (plate) has to complete a fixed number of maintenance-free tasks (missions), and each task has a few sub-tasks (maneuvers). Each sub-task is characterized by the limits (minimum and maximum) on the tensile load, and duration for which the load acts. The load levels can be seen as metric of the intensity of a given maneuver, whereas the duration defines the time for which a specific maneuver is sustainable to ensure maintenance-free operations. The goal of the load-profile optimization is to minimize the probability of exceeding the critical crack size at the end of the last mission while satisfying other constraints. We discuss a hybrid scheme that uses minimization of expected damage growth for the first few missions, and reliability-based optimization for the remaining missions. We use a surrogate-based optimization framework (Gutmann, 2001) to perform the optimization.

Probabilistic damage diagnosis
In this section, the methodology used for performing damage diagnosis using Bayesian information fusion for Lamb wave pitch-catch is discussed. The aluminum plate and sensor network considered in this work are depicted in figure 1.

Global sensitivity analysis and dimension reduction
A low-fidelity model (Lanza di Scalea & Salamone, 2008) is employed to simulate the Lamb wave actuation-propagationsensing under plane stress conditions to perform global sensi-  (Crawley & Luis, 1987) to obtain actuated shear stress in the plate, c) the mode expansion formulation by Giurgiutiu (Giurgiutiu, 2005) coupled with Fourier transform of the shear stress in the plate to compute the Lamb strain, d) a modified shearlag model for calculating the sensor response (strain) to the Lamb strain in the plate, and e) the integrated sensor strain multiplied by a sensitivity term to arrive at the sensor output voltage. Using this model and the values of the model parameters, the response (root mean square of the sensed signal) corresponding to the S0 mode at a given frequency can be obtained. For a given root mean squared (RMS) value of the input voltage, analyze the variance of the RMS value of the output voltage is analyzed. All parameters do not contribute equally to the uncertainty in the output voltage, and Sobol' sensitivity indices can be employed to estimate the contribution of the parameters to the output. For a model of the form is the vector of stochastic model inputs, the first-order (individual effect) index is given by: where x −i denotes all the model inputs other than x i . The total effect index is given by where V is the variance of the quantity in the parenthesis. The first-order index quantifies the individual contribution to the overall uncertainty in the output variable due to the uncertainty in an input variable without considering its interactions with other input variables, whereas the total effects index considers interactions with other input variables while estimating the same. Here, the relative values of these indices are used to identify parameters that have significant contribution to the uncertainty of the output quantity. We retain the important parameters as random variables and treat other parameters as deterministic variables at their mean values in subsequent analysis.

Surrogate modeling
The governing physics for the problem of interest involves the piezoelectric effect (Gauss' law for electricity) and elastic wave propagation in isotropic, thin plates (balance of momentum). A sampling-based algorithm is to be used to perform Bayesian diagnosis and information fusion. Thus, an efficient computational model that can provide the quantity of interest for a given sample of the parameters is needed. This is achieved by training a surrogate model (or a response surface) using the physics model. A finite element model for the governing physics is built using a commercial finite element program (Abaqus (Abaqus 6.14 Documentation, 2014)).
The governing physics is simulated using this model, the given set of input parameters, and given crack lengths to obtain output voltage signals necessary to perform the Bayesian estimation and information fusion. The damage sensitive feature (as a scalar quantity, y) of the output voltage signals is obtained and the training data set that contains the damage sensitive features corresponding to the chosen parameter values is obtained. A GP surrogate model is built using the training data set. Given the training data (inputs X, corresponding outputs y) and calibrated (trained) hyper-parameters of the GP (Θ c ), the output corresponding to a new input x * is a Gaussian random variable:

Bayesian information fusion
The goal is to estimate the crack length using the damagesensitive feature of the sensed (voltage) signal (y data ) and Bayesian estimation. To this end, the uncertainty in the knowledge of the crack length is expressed by means of a probability distribution function. A prior distribution (p prior ) of the crack length is assumed based on intuition, experience, etc.; and the knowledge is updated using the data by computing the likelihood function (p likelihood ) as: p post (a|y data ) ∝ p likelihood (y data |a) * p prior (a). (4) The Gaussian process (GP) surrogate model and a model of the measurement noise are used to obtain the likelihood function. The posterior distribution of the crack length is com-puted using a Markov chain Monte Carlo method. Each actuator-sensor path provides a different value of the damage metric, and hence, some additional information about the state of the damage (crack length) (figure 1). If the estimation is performed for a path (say) A2-S2, then the updated posterior for the crack length can be used as a prior for the Bayesian estimation for the next path (say A2-S3). The result at the end of the second Bayesian update is the result of the fusion of information contained in the data obtained from the two actuator-sensor paths. In our simulation experiments, the fusion of information obtained from Lamb-wave pitch-catch is performed along seven different actuator-sensor paths, viz. A2S2, A2S3, A2S1, A3S2, A1S2, A3S1, and A1S3 (see figure 1).

Probabilistic damage prognosis
If the current damage state and the associated uncertainty are known, then a damage prognosis methodology is needed to estimate potential damage growth (and associated uncertainty) for different system operations. We are concerned with fatigue crack growth in an aluminum plate (dimensions 380 mm x 150 mm x 0.81 mm) with an initial flaw under uni-axial, cyclic loading, as shown in figure 1. The plate is loaded on a 75-mm-wide region at the center of the shorter edge, and subjected to cyclic loads. For probabilistic crack prognosis, mode I crack propagation is considered under the uni-axial, tension-tension loading. Many empirical formulas for fatigue crack growth analysis are proposed in the context of LEFM with small-scale plasticity, for example, Paris' law (Paris & Erdogan, 1963), modified Paris' law (Donahue, Clark, Atanmo, Kumble, & McEvily, 1972), Forman's equation (Forman, Kearney, & Engle, 1967), the NASGRO model, etc. These models predict crack growth rate (da/dN ) as a function of the stress intensity factor range (∆K) and other model parameters. In this work, the Forman's equation, which takes into the consideration the effects of the stress ratio (R) and the fracture toughness (K c ) on the crack growth rate is used. Thus, the crack growth rate is modeled as: where m = 3.17, C is the parameter that can be calibrated using data from (laboratory) experiments, and K c = 67M P a √ m ("Fracture processes of aerospace materials", 2012). Given a cyclic load profile, a cycle-by-cycle estimation of the fatigue crack growth is performed. Thus, an estimate of stress intensity factors (SIFs) at the crack tip for the minimum and maximum load during a cycle for various crack lengths is needed. SIF can be expressed as an explicit function of crack size and stress level for specimens with simple geometry and loading conditions. However, analytical expressions for SIF for the specimen and loading of our interest (partially-loaded-edges, finite-width plate with a hole in the center and crack on both sides of the hole) are not available. Hence, finite element model analysis is necessary to obtain the SIFs of interest. To this end, a two-dimensional finite element model is created in Abaqus (Abaqus 6.14 Documentation, 2014) (assuming plane stress conditions) to reflect the operational conditions: one edge (partially) fixed and the other edge loaded only at the center 75-mm-wide region, and compute stress intensity factors for different crack lengths and loads using the contour integration technique (figure 2).

Figure 2. Aluminum plate (finite element) model used to determine SIFs
Obtaining the SIF using a finite element model at each cycle in the cycle-by-cycle analysis is a computationally expensive task for sampling-based, probabilistic fatigue crack growth prognosis. To expedite the process of SIF computation, a GP surrogate model that accepts the load and the current crack length as inputs, and provides the SIF as the output is used. Thus, the training points for the GP model include a twodimensional input (load, crack length) and one-dimensional output (SIF). A series of finite element model runs are performed with different combinations of crack sizes and loads to obtain the training and testing data set. The GP surrogate model, and crack growth model (Forman model) are used for probabilistic, cycle-by-cycle crack growth prediction. Diagnosis uncertainty and the variability in Forman model parameter (C) are the uncertainty sources considered in this work. At each cycle, the minimum and maximum loads are known and the current crack size distribution is inherited from the crack growth analysis of the previous cycle (or from damage diagnosis). The loads and samples from the current crack size distribution are used as the input of the GP model to predict the range of SIFs (∆K). These ∆K values are then used in the Forman's equation to predict the crack growth rate (da/dN ) distribution for the current cycle. The crack growth rate distribution is used to obtain the crack size distribution at the end of the current cycle. In this manner, probabilistic damage prognosis is performed given the candidate operational parameters (load profile).
Intelligent operation planning needs to be performed to ensure that the component completes the required tasks. In order to evaluate system performance, a metric is needed to measure the performance of the component. Here the mechanical work done during the loading phase of the cyclic loading is chosen as the required performance metric. The displacements on the left edge of the finite element model is recorded, and approximate the nodal force from the thickness t, the applied traction T and the element size ∆ as: Given the displacement and the loading, the work done by an applied load (F ) can be computed. The work done is a function of the applied load and crack size. Another Gaussian process surrogate model is built that outputs the work done in a cycle, given the loading and the crack length during a particular cycle. This GP model to evaluate a constraint in the load-profile optimization problem.

Load-profile optimization under uncertainty
A methodology to design intelligent, condition-based operations of a mechanical system that enhance system resilience is needed. To this end, we discuss load-profile optimization aimed at delaying damage growth while ensuring the required amount of work is performed by the system. The key assumptions of the load-profile optimization problem are listed below: 1. Maintenance-free execution of four tasks (or missions) is to be performed by the component.
2. Each task is divided into three sub-tasks. A given subtask (maneuver) is characterized by the corresponding minimum and maximum load levels. We also set the minimum and maximum duration for each maneuver for each mission. The loads for all maneuvers lead to constant stress ratio, R = 0.5.
For the optimization problem, the (constant) loading amplitude and the number of loading cycles (for each maneuver) are the decision variables. They represent the intensity and duration for which the maneuver is performed. Three different optimization strategies are possible a) minimization of the expected final crack size after each mission, b) minimization of the probability of exceeding a (predefined) critical crack size (a crit ) for each mission (reliability-based design optimization), and c) a hybrid strategy where the first two missions aim to minimize the expected value of final damage, and the next two missions use reliability-based design optimization. The probability of failure is very low in the first few missions, and the sampling-based optimization process may yield inaccurate results for this case. A hybrid strategy is thus adopted in this work. The optimization problem for this case can be stated as: Missions 3,4 where x is the vector of decision variables, θ is the vector of damage prognosis model parameters, g(x) is the performance function that denotes the work done described in section 2.2, E[a f (x, θ)] is the expected value of the crack size after each mission, and E[g(x)] is the expected value of the non-linear performance function (it evaluates the total amount of work done during each mission), W min represents the amount of that needs to be done for each mission, x lb , x ub represents the lower and upper bounds for the deterministic design variables, respectively, and P f = P [a f (x, θ) ≥ a crit ] is the probability of failure.

RESULTS AND DISCUSSION
In this section, we discuss results of parameter space reduction and surrogate model training required for probabilistic diagnosis and prognosis, and illustrate the methodology by performing load profile optimization for simulation experiments.

Global sensitivity analysis
Based on the material properties given in table 1 the global sensitivity analysis is performed. The input variables are assumed to be uncorrelated Gaussian random variables. Sobol' indices are computed using the method proposed by Saltelli et al. (Saltelli et al., 2010), using two sets of 10 7 random samples of the input variables. The results of the sensitivity analysis are depicted in table 2. It can be seen that the piezoelectric coefficient of actuator and sensor have the largest contribution to the variability of the output. Hence, we build the computational physics model for a sample of these parameters to create the training data for the GP surrogate model, and all other parameters are assumed to be deterministic variables at their mean values.

Surrogate modeling for diagnosis
We build a three-dimensional finite element model for a 0.81 mm thick aluminum plate (figure 1) with PZT-5J transducers as actuators and sensors. We use three-dimensional continuum finite elements (C3D20) to model the plate and Table 1. Mean (µ) and standard deviation (σ) values for input variables for sensitivity analysis 2.8 × 10 9 0.01 × 10 −3 74 × 10 9 500 × 10 −12 σ 1.8 × 10 8 0.01 × 10 −4 3.7 × 10 9 5 × 10 −11 Parameter lPZT (m) tPZT (m) Ep (Pa) ρp (kg/m 3 ) µ 7 × 10 −3 0.5 × 10 −3 71.7 × 10 9 2810 σ 7 × 10 −4 0.5 × 10 −4 3.58 × 10 9 140.5 Parameter Es (Pa) νs d31s (m/V) e33s (F/m 3 ) µ 7 × 10 −3 0.5 × 10 −3 500 × 10 −12 1.86 × 10 − 8 σ 7 × 10 −4 0.5 × 10 −4 5 × 10 −11 1.86 × 10 − 9 the adhesive used to bond the piezoelectric transducers to the plate. We utilize three-dimensional piezoelectric finite elements (C3D20E) to model the transducers. We use a Hanning-modulated, three-cycle long sine pulse (frequency 200 kHz) to excite the actuator. We run the model for twenty five samples of the important input variables (pertinent dielectric coefficients), and for six different values of the length of the crack growing out of the hole in the center. Thus, we perform one hundred and fifty runs of the finite element model to generate the training data for the GP model. We record the output voltage time series for each of the training points, for all crack lengths. We then compute a shortterm-Fourier-transform-based energy metric (J.  for S0 wavelet at central frequency (f c ) of 200 kHz, 250 kHz and 300 kHz. We treat the S0 wave energy (S0 e ) as the scalar output, based on which, the diagnosis is to be performed. Thus we build the GP surrogate model: and the subscript (a or s) denotes actuator or sensor properties, respectively. We repeat the procedure for each actuatorsensor path. The goodness of fit for all crack lengths for paths A2-S2 is shown in figure 3. It can be seen that the surrogate is sufficiently accurate. 3.3. Surrogate modeling for prognosis and work performance computation We perform 144 different finite element model runs, which cover the combinations of load of 1000-8000 lbs with a 1000 lbs increment and the crack sizes of 5-90 mm with a 5 mm increment. The combinations cover the whole range of the operational conditions. We randomly select 130 points from the data set to train the GP model, and use the remaining 14 points to test the GP model. The mean absolute errors for 14 testing points is less than 2%. Thus, the number of training points appears to provide high, converged accuracy for the surrogate model. The trained GP model is used to predict the SIF for different crack sizes and loading cases. The forcedisplacement data recorded on the loaded boundary is used to build another GP surrogate that estimates the performance constraint.

Load profile optimization: simulation experiment
In the simulation experiment, the probabilistic diagnosis is performed using ultrasonic guided-wave pitch-catch and synthetic data to obtain an estimate of the crack size at the beginning of a mission. The value of crack size obtained from    Table 6. Optimal design variables for intelligent mission planning probabilistic diagnosis, and the associated uncertainty are passed on to the load profile optimizer to design the optimal loading profile for the mission. The optimal loading is applied to a numerical model of the component to simulate the crack growth. This process is repeated for all four missions. In this manner, the simulation experiment is used to illustrate the integration of probabilistic diagnosis, prognosis, and load profile optimization. For the computational model representing the physical twin, the true value of C for the specimen is assumed to be C * = 1.15 × 10 −08 m/cycle. For the digital twin, uncertainty in the knowledge of the parameter C is modeled by assuming C as a Gaussian random variable with the mean value of 1.15 × 10 −8 m/cycle, and a coefficient of variation of 0.1154. The bounds for the design variables and the minimum work done constraint (W min ) for each mission are given in Tables 4 and 5. A critical crack size a crit = 10.0 mm is used for the last two missions. We remark that in a real-world application of the methodology, the critical damage size can be decided using the relevant standard of practice (e.g. Department of Defense Standard of Practice MIL-STD-1531Dc1).
It is assumed that a crack size a 0 = 5 mm represents the initial damage in the component. For the probabilistic diagnosis, synthetic data generated from the computational model is fed to the Bayesian estimator to obtain the crack size and diagnosis uncertainty. The diagnosis results are obtained using Markov-chain Monte-Carlo method (Metropolis-Hastings algorithm (Hastings, 1970)) . In the GP surrogate model (equation 9), the two dielectric coefficients (d 15a , d 15s ) are assumed to be fixed at their mean values and the likelihood of a candidate crack length is computed. A Markov chain of 10 5 Monte Carlo samples is constructed and the initial 10000 samples are rejected (initial burn-in samples).
The results of load profile optimization are shown in Table 6. It can be seen in Table 6 that the optimizer successfully attained the minimum work performance constraint (E[g(x)] > W min ) for all missions. It can be seen from the results that for some of the maneuvers, the optimizer chose the lowest possible load value to complete the mission. This can be explained as follows. The work requirement constraint is approximately proportional to the square of the applied loading, whereas the crack growth is approximately propor-tional to the m-th order of applied loading (m ≈ 3). Thus, utilization of lower load for a large number of cycles may be sufficient to complete a mission (provided the time constraint allows for the required number of cycles).

CONCLUSIONS
We discussed a methodology for performing health management and operations related decisions for a mechanical system. This was achieved by designing operational parameters (loading intensity) for a simple mechanical component such that the damage growth in the component was minimized, while the component performed desired work. We considered three key aspects of intelligent decision making for a mechanical system: probabilistic damage diagnosis, probabilistic damage prognosis, and system optimization under uncertainty. We built digital twin of the component of interest using multi-physics, multi-fidelity models and considered aleatory uncertainty, as well as diagnosis uncertainty in our analysis. We explored a hybrid approach that combined crack growth minimization with reliability-based design optimization. We showed that the proposed methodology can be successfully used, with the help of a digital twin for the system of interest, to perform operational planning to achieve the desired system performance goal.