Crack Identification of Cantilever Plates Based on a Kriging Surrogate Model

Cantilever plates appear in a wide range of structural systems in civil, ocean, and aeronautical engineering, which may be subjected to known and unknown cyclic loads. Consequently, fatigue cracks may be initiated. A crack can lead to reduced performance and shorter lifetime, and, even worse, may induce catastrophic failure of structures, which should be considered seriously in both theoretical and experimental research.

A significant number of methods for health monitoring and crack identification have been proposed during recent years [1,2]. In the early stages of the crack detection development, the frequency contour plot method was one of the most popular techniques to identify a single crack in a beam by using the first few natural frequencies. Gudmundson [3] discussed the effect of geometrical imperfection on the eigenvalue by means of perturbation analysis. Liang et al. [4] proposed a massless rotational spring model for a crack based on the Euler–Bernoulli beam theory, in which the location and size of a crack could be identified via finding the intersection point of a few frequency contour lines. Ostachowicz and Krawczuk [5] studied the forced vibration of a beam and the effects of the crack locations and sizes on the vibration behavior of the structure. Lele and Maiti [6] extended the frequency contour plot method in beams based on the Timoshenko beam theory. However, the frequency contour plot method suffers from the drawback that the curves of the frequency contour plot might not intersect because of inaccuracies in the model with respect to the measured results on practical structures.

Many efforts have been devoted to crack identification in beams. However, there are few works on crack identification for plates. In the work by Krawczuk et al. [7], special spectral elements were used to analyze waves in cracked plates. Horibe and Watanabe [8] applied a genetic algorithm to determine the crack location using finite element models. Hadjileontiadis and Douka [9] presented an effective method for detecting cracks in plate structures based on fractal dimension analysis. The location and length of the crack were identified as the abrupt changes in the spatial variation of the 2D fractal dimension signal for plates with cracks of varying lengths at different locations. In the paper by Lam and Yin [10], cracks were modeled as a set of linear springs with varying stiffness and a two-stage method for multicrack detection was proposed using statistical analysis. One shortcoming of these methods is that they are only suitable for cracks that are parallel to the plate edge, which of course is a limitation in real world situations. In order to overcome this shortcoming, Bagheri et al. [11] investigated the use of curvelet transform to identify the location of damage. More recently, Moore et al. [12] used a Markov–Chain Monte Carlo implementation of Bayes' rule to estimate the crack parameters such as location, orientation, and size. Although this method gave a good prediction with acceptable accuracy, it required a large amount of samples and hence computation time.

In essence, crack identification is an inverse problem. Optimization methods are commonly employed to estimate the most suitable crack parameters by minimizing the objective function related to the output discrepancies such as natural frequencies, model shapes, and dynamical responses. It is desirable to construct a simple relationship between the variables and the objective or constraint functions to avoid the repeated analysis of computationally expensive FE models during the optimization process. Recently, surrogate models based on neural networks (NNs) [13,14] and the polynomial response surface (PRS) method [15,16] were investigated by many researchers. Atalla and Inman [13] employed a neural network to estimate the updated model parameters. Real and imaginary components of the frequency response are integrated over selected frequency intervals and used to provide information about the system to the network. Lee et al. [14] used mode shape differences or mode shape ratios between before and after damage as the input to the NN to reduce the effect of the modeling errors in the baseline FE model, from which the training patterns are to be generated. While the NNs-based method is successful to some extent, the required number of training samples would exponentially increase, resulting in considerable computational efforts.

In the thesis by Cundy [17], a four-step process based on PRS was developed for simple physical systems and damage identification was performed successfully, given the limited amount of “training” data used. Faravelli and Casciati [18] and Casciati [19] employed acceleration time histories collected under different loading conditions for identifying the presence of damage and locating cracks by a comparison of the sum of the squared errors (SSE) histograms. Moreover, a derivative of SSE, defined as ADM, was proposed and used for detecting distributed cracks in an actual masonry [20]. Huang et al. [21] stated that PRS method not only tended to show severe oscillations, but also required too many support points, which might result in limited accuracy when the response data to be modeled had multiple local extrema.

Compared with NNs and PRS, Kriging surrogate models might offer a good alternative. It provides explicit functions to represent the relationships between the inputs and outputs of a linear or nonlinear system for response estimation and parameter identification. The use of Kriging models often requires only a small number of samples and can reduce significantly the computing time in obtaining the optimal parameters. Therefore, Kriging models have drawn much attention and have already been increasingly used in industrial design [22]. However, to date, little research has been carried out in applying Kriging models to crack identification in plates. Bilicz et al. [23] provides a methodology for the characterization of a 3D defect embedded in a conductive nonmagnetic plate from the measurement of the impedance variations of an air-cored pancake coil at eddy current frequencies based on Kriging interpolation using the Expected Improvement algorithm. However, the boundary of the defect model must be parallel or perpendicular to the plate, which somewhat diminishes its utility. Other adaptive metamodels (e.g., radial basis functions, RBF) for crack characterization also can be found in Refs. [24,25].

The objective of this work is to apply a Kriging model to single arbitrary crack identification for plate-type structures. Firstly a number of FE models with varying tip locations are constructed in Ansys code to solve dynamical responses. To speed up the spectral analysis, PEM is used to guarantee accuracy and reduce computational cost. The initial Kriging model is then constructed by the samples of crack parameters and their corresponding RMS values at specified points. To estimate the actual crack parameters, the SPSO algorithm is employed, which is demonstrated to be effective, accurate, and robust for searching for a global optimal solution. After that, the current optimal solution will be used in FE analysis and inserted into the initial sample set to update the initial Kriging model, until the surrogate model is sufficiently accurate and the optimization process converges. Finally, the effectiveness and efficiency of the presented method is demonstrated by numerical studies on a cantilever plate with an internal crack in the presence of noise. The identified results show that the proposed method holds considerable promise for practical engineering problems.

A model of an elastic, rectangular plate with an internal crack at an arbitrary position is shown in Fig. 1. The plane dimensions of the plate are a and b in the x and y directions, respectively. x_i and y_i (i = 1,2) represent the true coordinates of the two crack tips in a Cartesian coordinate system. For simplicity, the nondimensional locations are used in the form of ξ_i = x_i/a and η_i = y_i/b (i = 1,2). Hence, the values of crack parameter are limited in the interval (0,1).

In previous works [9,10], a crack was modeled as massless line springs with varying stiffness along the crack, which divided the plate into two segments. Then a hypothetical displacement function and the governing differential equations of the plate segments was established to estimate the crack parameters such as location and size. However, this process is likely to be awkward when the crack slants to the plate edge because of the difficulties in describing the displacement field using the displacement function. As such, finite element techniques may be a promising choice for the construction of a plate model with an arbitrary crack. In the current study, the finite element models are established with a mesh of isoparametric elements that contain a singularity. As shown in Fig. 2, the elements away from the crack tips are eight-node quadrilateral, and the elements surrounding the tips are triangular with three nodes. It should be mentioned that it is assumed that the crack remains always open without changing the stiffness and mass properties of the plate during the motion, and the contacting nonlinearity and the crack propagation process are beyond the scope of our paper.

In order to assess the accuracy of the cracked plate model, a comparative study on a completely free square plate with an internal crack is carried out. Table 1 displays the first five nondimensional frequency parameters $ωa2ρh/D$ calculated by using the finite element method and Ritz method [26] in six different cases. ω is natural frequencies of the plate, h is the thickness of the plate, ρ is the mass density of the plate; D = Eh³/12(1−ν²) is the flexural rigidity, E is the modulus of elasticity, and ν is Poisson's ratio.

Natural frequencies as an overall global characteristic of a mechanical system are not very effective for crack identification in plates. Location-based characteristics should be considered to locate the crack tips. For this reason, vibration amplitudes and impulse responses are widely used for damage detection. However, practical engineering problems are very different from simulated cases due to the presence of modeling errors and measurement noise. As observed by the investigators of Refs. [27,28], these outputs are somewhat sensitive to environmental changes and noise interference, which may lead to unreliable results of crack localization. Considering this, some researchers used stochastic response analysis for damage detection and system identification [29–31]. Liberatore and Carman [32] suggested that the RMS values of the power spectral density (PSD) could be adopted to identify damage in beams and give satisfactory predictions even for noisy cases. Inspired by their work, the RMS values at specified points are used in this work due to its strong capability in describing local features and against noise.

After the cracked plate models with varying tip locations are established, excitation and measurement points are located as shown in Fig. 3. The plate is excited by a simulated white-noise force spectrum with a constant amplitude. With four measurement points and four excitations points, 16 acceleration spectra are calculated for each damaged case.

The area under the acceleration PSD curve is the mean square (MS) value. Its square root (RMS) will be used to identify the crack parameters.

Varying tip locations can result in different FE meshes. To calculate the corresponding RMS values based on FE models, remeshing should be made at every iterative step of optimization. This is awkward for complex problems that need a great deal of iterative steps. Therefore it is desirable to construct a simple relationship between crack parameters and the associated RMS values before optimization for avoiding the expensive FE analysis.

where $rmsl*$ and $rmsl(FE)$ are the lth component of the RMS vector of surrogate model and the true value calculated via FE analysis, respectively; $rms(Average)$ is the mean of all true values.

where rms^* is a vector of the predicted root mean square values at a new set of crack parameters x^* related to the locations of crack tips, and $rms(Target)$ represents the given ones which may be measured on the real cracked plate. lb and ub denote the lower and upper bounds. w is the weighting imposed on the components of rms^*. In this study, all its element are set to 1.

where c₁ and c₂ are two positive constants called acceleration coefficients, P_i and P_g are local and global best locations, respectively, and r₁ and r₂ are random numbers in the interval (0,1). If P_g = P_i, S_i(t +1) = S_i(t). As such, the particle at the global best position will stop evolution. To improve the global searching ability, this algorithm randomly generates an extra particle labeled j with position S_j to continue evolution in the search domain. It means that at least one particle is generated in the searching domain randomly to improve the global searching ability and final searching quality of PSO algorithm.

When the optimal solution from the initial Kriging model is obtained, the FE analysis should be carried out by using these parameters to assess the accuracy of the initial Kriging model. Owing to the small number of samples, the initial surrogate model should be updated. The optimal solution and its corresponding RMS values by FE analysis will be then inserted into the sample database. The Kriging model will be rebuilt by this point-adding process until it is sufficiently accurate according to SMC and EISE criterion. The proposed method can be organized as a surrogate model updating and optimizing estimation procedure as illustrated in Fig. 4. The identification of crack tip locations consists of the following steps:

1. Step 1: Generate initial sampling parameters X (crack tip locations) and run the FE analysis program to obtain the corresponding output RMS.

2. Step 2: Construct the initial Kriging surrogate model with samples and output RMS values obtained in Step 1.

3. Step 3: Find the optimal parameters $xk*$ by minimizing discrepancies between the targeted values $rms(Target)$ and the calculated rms^* by means of SPSO algorithm based on the initial Kriging surrogate model, and set the iterative index k = 1.

4. Step 4: Check criteria: If the current RMS values predicted by Ansys and the one based on Kriging surrogate model satisfy the given criteria, then stop updating surrogate model and predict crack parameters with this surrogate model; else, go to Step 5.

5. Step 5: Add $xk*$ and $rms(xk*)$ into the initial sampling database X and RMS generated in Step 1, respectively, then update the current Kriging model and update the iterative index: k = k + 1.

6. Step 6: Loop to Step 3 and repeat the process till the criteria are satisfied, and output the optimal crack parameters.

To demonstrate the effectiveness and efficiency of the presented method, a numerical study of a cantilever plate with an internal crack is carried out. The plate under consideration has a square area of 1.4 m × 1.4 m and thickness of 0.01 m. The material properties of the steel plate are: Young's modulus E = 210 GPa, Poisson's ratio ν = 0.29, and density ρ = 7850 Kg/m³. All the DOFs on the cantilever side are fixed in the FE model. Fifty samples are firstly generated via Latin Hypercube design [38] and prepared for the construction of the initial Kriging model by varying the locations of crack tips. Forced vibration tests are simulated assuming a random white-noise force spectrum with constant amplitude 0.004 kgf²/Hz acting at four specified points, respectively. Consequently, sixteen RMS values can be extracted and expressed as a 1 × 16 response vector for each case, so that q = 16. The concerned frequency range selected to compute the RMS is chosen between 1 Hz and 180 Hz. The Kriging model is constructed with $f(xi)$ in Eq. (7) taken to be a constant one for any sample.

Once the initial Kriging surrogate model has been constructed, the SPSO algorithm is used to estimate the optimal crack parameters based on the surrogate model by minimizing the objective function with optimization parameters given in Table 2. Meanwhile, the accuracy of the current surrogate model is checked according to the given criteria (SMC > 0.996 and EISE < 0.005) to decide whether the model should be updated.

where d denotes the normally distributed random values between −1 and 1; the noise level $nM=0.03$ is used. $yωM$ and $yωS$ represent measured and simulated harmonic responses, respectively.

Six cases are considered to assess the performance of the proposed method, as shown in Fig. 5. These cases cover typical crack situations. The actual and identified crack parameters are reported in Table 3. Very good agreement between the actual and the identified locations of crack tips based on the Kriging surrogate model is obtained. It is worth noting that the proposed method provides a correct identification of crack tip locations even for noisy cases, and it can be seen that the Kriging surrogate model seems to be insensitive to random noise.

Figure 6 illustrates the convergence of the objective function with the Kriging model updating during crack identification. The horizontal axis is the number of updating steps and the vertical axis is the value of the objective function. The relative error of the identified nondimensional crack parameters are within ±0.62%. The cost values seem to be reduced in general with the iteration of surrogate model updating, which illustrates the process of surrogate model updating is essentially a continuous improvement in describing the changing trends of the surrogate model around the current optimal solution. It is worth noting that after certain steps, some cost values are increasing (e.g., k = 5 in case 1, k = 3,4,7 in case 3, and so on), which indicates another minimum may exist around this optimal solution. That is to say, the current solution, which will be selected as an additional sampling point, may not be the global minimum. Consequently, it is stated that the accuracy of the surrogate model and the global searching ability of the algorithm are equally important in the whole crack identification process.

The values of SMC and EISE at the last iterative step for each case are shown in Fig. 7 and it can be found that the final Kriging model, after several updating steps, is accurate enough for the estimation of crack parameters. Figure 8 compares, for each of the considered cases, the RMS values predicted by the final Kriging model and corresponding FE model.

To investigate the effect of sampling size on the precision of the identified results, all the cases are recalculated by using (75, 100, 125, 150, 175, 200) samples given in the Appendix, respectively. The identified crack parameters by using 100 and 200 samples are listed in Table 4. In addition, the required surrogate model updating steps are compared in Fig. 9. It seems that the initial sampling sizes do not directly affect the precision of the identified crack parameters, but do change the required updating steps. From Fig. 9, it can be seen that the optimal solution can be found with k = 1 when 200 samples are initially used for the construction of the surrogate model so that there is no need to carry out the updating process in Cases 1, 3, and 6, while the initial Kriging model cannot provide the “best” solution for other cases. In this sense, the surrogate model is not perfect for describing the relation between crack parameters and RMS values around the optimal solution. As such, more samples should be used and a process of surrogate model updating should be then considered.

All the results indicate that the Kriging model is a powerful tool to give a quantitative estimation of how the locations of crack tips affect RMS values. The actual crack parameters can be determined by this method.

An effective method based on a Kriging surrogate model for crack identification in a plate is presented in the paper. The results clearly show its usefulness and accuracy. Some appealing features of the presented method include: (1) the Kriging model is applied to provide a simple relation between crack parameters and corresponding RMS values to avoid the remeshing problem at every iterative step of optimization, (2) the values of RMS at specified points of the plate are used as output responses to construct the surrogate model, which is impervious to random noise and has seldom been used in crack identification, and (3), PEM, instead of traditional spectral analysis, is used to produce accurate results and reduce computational cost.

It should be mentioned that only simulated noise is considered in this work. There may be some undesired discrepancies (also known as “model error”) between FE model predictions and experimental results due to inevitable uncertainties in material properties and boundary conditions. FE model updating may be a better way to solve this problem. Therefore, the model updating strategy should be considered before crack identification for use in real applications. In addition, only one crack was considered in this paper. More cracks with more complex shapes also may appear in practical problems. As cracks can be complex, more crack parameters should be considered to describe the characterization of intercrossing between two or more arbitrary cracks. This issue will be discussed in the future.