Crack detection in a pipe by adaptive subspace iteration algorithm and least square support vector regression

2.1. Stable completion

The stable completion [30] is a method to construct interpolating multiwavelets with specific defined properties in the multiresolution finite element space. Comparing to the lifting scheme [19-20], the stable completion is more flexible and efficient to construct new multiwavelets with high vanishing moments, which is important for the scale decoupling of stiffness matrix and mass matrix in the multiscale computation. The stable completion assembles a new compact multiwavelet as a linear combination of several primitive wavelets and interior scaling functions:

where $g_{j,m,m^{\text{'}}}$ and $s_{j,k,m}$ are the stable completion coefficients with $m^{\text{'}}\in M\left(j\right)$. The stable completion provides us much freedom to design the wavelets by choosing appropriate stable completion coefficients. The stable completion coefficients can be computed by computing a basis for the null space of the local interaction matrix $R_j$ in the scale-orthogonal relationship equations:

where, ${\varphi }_{j,k^{*}}$ are all the scaling functions on a given domain ${\mathrm{\Omega }}_j$, ${\varphi }_{j,k}$ are the subset of scaling functions that are all interior in the domain ${\mathrm{\Omega }}_j$ and $a(•,•)$ is the operator in the differential equations, engineering problems, etc. According to the scale-decoupling conditions in Eq. (2), we can get the principle of construction of the scale-decoupling wavelet bases according to the inheritance of vanishing moments property [23].

2.2. Rayleigh-Timoshenko beam-type multiwavelets

According to Rayleigh-Timoshenko beam theory, the equation of motion is of the form [31]:

where $w$ denotes the transverse displacements, $\theta$ the rotation of cross-section, $\rho$ the mass density per unit length, $A$ the area of cross section, $E$ the modulus, $I$ the second moment of area, $G$ the shear modulus, $K_s$ is the shear correction factor. The operators corresponding to the stiffness matrix in the element free vibration have two types of multiwavelets in the form:

where $i\in J$. The operators for mass matrix are:

According to the operator of Rayleigh-Timoshenko beam problems, we can construct multiwavelets based on the stable completion method, which ensure the scaling functions and wavelet functions are orthogonal with respect to the operator (5)-(9). Figs. 1(a) and 1(b) show quadratic and cubic Rayleigh Timoshenko beam-type multiwavelets on the support (–1, 1), respectively.

Fig. 1Rayleigh Timoshenko beam-type multiwavelets: a) Quadratic, b) Cubic

a)

b)

2.3. Rayleigh-Euler beam-type multiwavelets

Considering rotary inertia of beam, the governing equation of transverse vibration of Rayleigh-Euler beam is (Reddy 2004 [31]):

where $\rho$, $A$, $E$ and $I$ have the same meaning as Rayleigh-Timoshenko beam. The operators for Rayleigh-Euler beam element can be derived as:

Fig. 2(a) shows a cubic Euler-Bernoulli beam multiwavelet with three vanishing moments on the support (0, 3). For quintic Hermite scaling functions, we can construct new multiwavelets on the support of two elements with two vanishing moments as shown in Fig. 2(b). To satisfy the operator-orthogonality in Eq. (12), new multiwavelets must be constructed with two more vanishing moments.

Fig. 2Rayleigh Euler-Bernoulli beam-type multiwavelets: a) Cubic, b) Quintic

a)

b)

3.1. Pipe crack model

Fig. 3 shows a transverse crack of depth $h$ is considered on a pipe and the crack section. The total length of the pipe is $L_t$, the lengths of pipe support are $L_a$ and $L_b$. The inner and outer diameter of the pipe are $D_a$ and $D_b$, resectively. The normalized crack location and depth are $\beta =l/L$ and $\alpha =h/D_b$, resectively. The rotational spring model is used to model the cracked pipe. To overcome the difficulty in calculating local flexibility in the model, we consider a cracked pipe as a combination of a series of thin annuli. The stress intensity factor of a thin annulus can be calculated by analytical method presented by Tada et al. [32]. The stress intensity factor can be resolved to the desired accuracy when the number of divided thin annulus is large enough.

Considering the crack as a non-dimensional rotational linear spring [7-8], the cracked stiffness matrix can be written as:

where $c_m$ is dimensionless compliance matrix calculated by local flexibility according to the strain energy theory in Reference [33]. The cracked stiffness submatrix $K_c$ can be assembled into the global stiffness matrix and the cracked structural model based on the finite element multiwavelets is constructed. The solution of the eigenvalue problem for pipe structure can proceed as that of uncracked pipe structure.

Fig. 3Simply supported pipe with a crack

3.2. Adaptive subspace iteration algorithm (ASIA)

The eigenvalue equation for pipe structure:

can be discretized by interpolating multiwavelets at level $j$ and the generalized eigenvalue problems can be obtained as:

where $K_{j+1}$ and $M_{j+1}$ are multilevel stiffness and mass matrices at level $j+1$, ${\lambda }_{j+1}$ and ${\omega }_{j+1}^{}$ are the eigenvectors and eigenvalues at level $j+1$. The multilevel stiffness matrices $K_{j+1}$ and $M_{j+1}$ of eigenvalue equation can be denoted in the two-level form:

and:

Since the reference or exact solution may not be found in general engineering eigenvalue problems, the eigenvalue solution at level $j+1$ is used to replace the exact eigenvalue solution and the error estimator is given:

Given an initial domain ${\mathrm{\Omega }}_j$, an initial scale $j=\text{0}$, a global tolerance $\epsilon$ and an iteration tolerance $\tau$, an adaptive subspace iteration algorithm is presented to solve structural eigenvalue problems.

The multiscale method for eigenvalue problems is programmed as follows:

(1) Select initial vector $X_j^0$ at level $j$, let $Y_j=M_j{X}_j^0$. The initial eigenvector can be chosen randomly, but not orthogonal to ${\lambda }_j$, i.e. satisfy the equation ${\left(X_j^0\right)}^T{M}_j{\lambda }_j\ne 0$.

(2) Solve the equation $K_j{X}_j^1=Y_j$. The scale-decoupling property of multilevel stiffness matrices ensures that it is unnecessary to factorize the stiffness matrices at each level.

(3) Construct new stiffness and mass matrices from:

and let ${Y^{\text{'}}}_j=M_j{X}_j^1$, then:

(4) Compute the approximate value of ${\mathrm{\Lambda }}_j^{*}$ from:

where:

(5) Compute the iteration error estimator ${\eta }_{j,i}$ in the form:

and test whether it satisfies ${\eta }_{j,i}<\tau$. If it does not, go to step (6), otherwise, go to step (7).

(6) Regenerate $Y_j={Y^{\text{'}}}_j{\lambda }_j^{*}$., go to step (2).

(7) Let ${\omega }_j^{*}={\left({\mathrm{\Lambda }}_j^{*}\right)}^{-1/2}$, compute the global tolerance ${\xi }_j$ and test whether ${\xi }_j<\epsilon$. If it does not, go to step (8), otherwise, go to step (9).

(8) Calculate the local error estimators ${\delta }_{j,n}$ in each local domains ${\mathrm{\Omega }}_{j,n}$ in the following form:

and let ${\delta }_j^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{a}\mathrm{x}\left({\delta }_{j,n}\right)$. If ${\delta }_{j,n}\ge \theta {\delta }_j^{\mathrm{m}\mathrm{a}\mathrm{x}}$, the wavelet bases into the local domains and form new stiffness and mass matrices at the level $j$, go to step (1).

(9) Output the eigenvalue ${\omega }_j={\omega }_j^{*}$ and the eigenvectors are computed in the form:

3.3. LSSVR algorithm

According to the structural risk minimization principle, the regression problems can be expressed as a constrained optimization problem in the form [28-29]:

where $e_i$ is the error. The solution of LS-SVM regressor will be obtained after we construct the Lagrangian function:

and a system of linear equations can be derived based on KKT condition:

where $P={[p_1,\cdots ,p_N]}^T$, $1_N={[1,\cdots ,1]}^T$ and $u={[u_1,\cdots ,u_N]}^T\text{.}$ Here, $I_N$ is an $N$×$N$ identity matrix, and $\mho \in {\mathbb{R}}^{N\times N}$ is the kernel matrix defined by:

The kernel function $\mathrm{\Xi }\left(•,•\right)$ typically has several options, such as the linear, polynomial, radial basis function (RBF) kernel, etc. We choose the RBF as the kernel of LS-SVR.

3.4. Crack detection

According to linear fracture mechanics theory, the additional localized flexibility in crack vicinity can be represented by a lumped parameter element [16]. With the crack propagation, the local flexibility and dynamic characteristics of the pipe structure will change accordingly, which can be used to identify crack location and size. The determination of natural frequencies $f$ for a given crack location and depth can be obtained by the inputs of normalized crack location $\beta$ and depth $\alpha$. The relationship between natural frequencies $f_r,$$r=\text{1, 2, 3,…}$ and crack parameter can be denoted as:

where $F_r$ is the function relation. The identification inverse problem of pipe crack, determining the normalised crack location $\beta$ and depth $\alpha$, can be written as:

Usually, the crack location and size of the pipe structure can be identified by the first two frequencies. Because natural frequencies change a little when the crack appears at the node of corresponding order, the first three frequencies are always used as input parameters of inverse problem to identify structural crack.

The pipe crack can be accurately detected by the combination of multiwavelet-based element method and LSSVR algorithm. The identification process can be summarized as following.

(1) The cracked structure model was built using multiwavelet-based elements and the database of pipe structure employing first three natural frequencies of the pipe structure with different location and size are solved by ASIA.

(2) Construction of prediction model of pipe crack. The training sample set for the LSSVR model is built as $\left\{\begin{array}{cc}{f}_i& \alpha \end{array}\right\}$, $\left\{\begin{array}{cc}{f}_i& \beta \end{array}\right\}$ by use of the detection data and the regression parameter was solved by LSSVR algorithm.

(3) Extracting the first three frequencies of the pipe structure form vibration experiments based on FFT method.

(4) The first three measured natural frequencies are employed as input parameters and identifying the crack location and depth by the trained LSSVR.

4.1. Numerical investigation

To verify the validity of multiwavelet-based Rayleigh Euler and Rayleigh Timoshenko beam elements, a typical pipe with different length under different boundary condition are illustrated. The first three frequencies of uncracked and cracked pipe models are solved by ASIA. To predict the crack location and depth, multiwavelet-based cracked pipe models are used to detect the crack in a pipe under a given prediction error estimator.

Example 1. Given a simple supported pipe, the pipe geometries and the material properties are as follows: pipe length $L=$ 1 m, $La=$ 0.005 m, $Lb=$ 0.005 m, inner diameter $D_a=$0.05 m, outer diameter $D_b=$0.06 m, Young’s modulus $E=\text{1×}{\text{10}}^{\text{11}}{\text{N/m}}^{\text{2}}$, $G=\text{6×}{\text{10}}^{\text{10}}{\text{N/m}}^{\text{2}}$, material density $\rho =\text{7900}{\text{kg/m}}^{\text{3}}$, Poisson’s ratio $\mu =\text{0.3}$, shear correction coefficient $K_s=$5/6.

Fig. 4LSSVM computation

The first three frequencies of the pipe crack model are obtained by using ASIA. Fig. 4 shows the LS-SVM results of the crack detection model in the region of the training points by using a Matlab/C toolbox for least squares support vector machines named LS-SVMlab [34]. In order to obtain an LS-SVM model (with the RBF kernel), two extra tuning parameters are given in Fig. 4, where $\gamma$ is the regularization parameter, determining the trade-off between the training error minimization and smoothness of the estimated function and ${\sigma }^2$ is the kernel function parameter. Let the error estimator of first three frequencies be 0.001 %, Table 1 illustrates the first three natural frequencies using ASIA and those of closed-form solutions in Reference [13] under various crack location and crack depth. It shows that the numerical solutions are in very good agreement with the closed-form solutions. Table 2 show the frequency error estimators in different spaces using quadratic Lagrange multiwavelets under $\beta =$0.5 and $\alpha =$0.5. Fig. 5 shows the error estimators of first three natural frequencies for the quadratic Lagrange multiwavelet solution. It can be seen that the solution of natural frequencies is convergent as the scale becomes larger and the convergence rate of the solution of first natural frequency is highest in that of three natural frequencies.

To detect the crack in the pipe, the exact first three frequencies are used to be inputs of cracked pipe models based on the cubic Hermite multiwavelets shown in Fig. 2. Table 3 shows the prediction results of crack parameters based on the cubic Hermite multiwavelets while the error estimator of actual parameters is set to be 0.001 %. It can be concluded that the prediction results of crack parameters can efficiently approximate to high precision by ASIA. It can be found that the relative errors of the predicted normalized crack location and size are not more than 0.01 %.

Fig. 5Error estimators of first three natural frequencies for the quadratic Lagrange multiwavelet solution

4.2. Experiment investigation

In the vibration experiment of the cracked pipe, a Polytec Doppler laser vibrometer OFV-505/5000 and Sony EX system were used to measure the vibration of one point in the pipe. Force-hammer is used as current pulse exciting source. To reduce the reflection of the laser beam and spectral noise, retro-reflective tapes were put on the measurement point in the pipe. For each cracked pipe, the high-metrical frequencies can be obtained by using Fast Fourier transform (FFT) and spectrum zoom technique of Matlab R2010b. Figs. 6 and 7 show the experiment setup and cracked pipe, respectively. The pipe geometries and the material properties are as follows: $Lt=$ 1.450 m, $La=$ 0.005 m, $Lb=$ 0.005 m, inner diameter $D_a=$0.05 m, outer diameter $D_b=$0.06 m, Young’s modulus $E=\text{2.06×}{\text{10}}^{\text{11}}{\text{N/m}}^{\text{2}}\text{,}$ material density $\rho =\text{7390}{\text{kg/m}}^{\text{3}}\text{,}$ Poisson’s ratio $\mu =\text{0.3}$. The crack of each pipe is generated by wire-cut machining with a molybdenum wire of diameter 0.18 mm.

Fig. 6Experimental setup

Fig. 7Cracked pipe

The quintic Hermite multiwavelets are used to construct the cracked pipe model, which is solved by ASIA. For the detection of pipe crack, the first three experimental measured frequencies are emplolyed as inputs of the prediction problems of pipe crack. For the prediction of crack location and depth, LSSVR algorithm is applied. For each cracked pipe, the high-metrical frequencies can be obtained by using the standard FFT program. A typical record including the impulse response signal and the power spectrum for $\beta =\text{0.480}$ and $\alpha =\text{0.308}$ is shown in Fig. 8. Table 4 shows the comparison of actual normalized crack parameters $\beta$ and $\alpha$ and the predicted crack parameters ${\beta }^{*}$ and ${\alpha }^{*}$. For the given cases, the relative errors of ${\beta }^{*}$ are not more than 5.31 % while the relative errors of ${\alpha }^{*}$ arrive at 11.82 %. Therefore, the proposed model-based crack detection method by the combination of multiwavelet-based elements and LSSVR algorithm is considered to be valid for actual application to detect fatigue cracks in a pipe.

Fig. 8FFT analysis on cracked pipe with β=0.480 and α=0.308