Identification of structural parameters including crack using one dimensional PZT patch model

Abstract

This article presents a new concept of using the one-dimensional piezo-electric patch on beam model for structural identification (SI). A hybrid element constituted of one-dimensional beam element and a PZT sensor is used with reduced material properties. Accuracy of this element is first verified against a corresponding 3D finite element model. Then SI is carried out as an inverse problem whereby parameters are identified by minimizing the deviation between the predicted and measured voltage response of the patch, when subjected to impulse excitation. A non-classical optimization algorithm Particle Swarm Optimization is used to minimize this objective function. Identified parameters involve stiffness, damping as well as the depth and location of crack in a beam. The validity of the proposed approach is proved by numerical studies on a beam, nine member frame and crack depth and location identification using various patch lengths. The signals are polluted with 5% Gaussian noise to simulate experimental noise. The results show there is a significant improvement in identification accuracy compared to other methods. The proposed method is also successfully verified experimentally.

Keywords:

• Structural identification
• voltage matching
• PZT patches
• inverse problem
• Particle Swarm Optimization

1. Introduction

Structural identification (SI) is a class of inverse problem which deals with how to estimate structural parameters such as stiffness and damping from the measurement of input excitation and the corresponding output response. SI takes part important role in model updating and structural health monitoring. SI algorithms are generally categorized into the frequency, modal and time domain, parametric and nonparametric models, and deterministic and stochastic approaches, and classical and non-classical methods.[1,2]

Chakraborthy and Mukhopadhyay [3] estimated simultaneously elastic parameters and cross-sectional dimensions of a stiffened plate from experimental modal data and Finite Element (FE) predictions, using model updating technique. The method used mode shape and frequency error function between simulated ‘experimental’ and FE model values were minimized by Levenberg Marquardt algorithm. Kianian et al. [4] investigated damage detection of offshore jacket platforms exploiting a method which applies frequency response function data-set in some frequency points and also the proper methodology applied to the selection of frequency points to improve the damage detection efficiency. Yang and Wang [5] introduced a new damage detection method based on the new conception of natural frequency vector and natural frequency vector assurance criterion and verified by both simulative and experimental examples. Wang et al. [6] suggested a two-stage identification algorithm utilizing the changes in natural frequencies and static displacements, for identifying the structural damages. Liu et al. [7] identified the presence of structural damage using multi-objective optimization, keeping variations in natural frequency and mode shapes as individual objective functions. However, all these methods based on frequency domain may not be suitable for the large complex system which requires high frequency resolution, local damage assessment or non-linear parameter determination.[2] Moreover, damping properties of structures cannot be estimated and also need some transformation of dynamic response data which requires more computational effort.[8]

In this respect, time domain algorithm SI schemes are a better alternative and it has an advantage that structural response can be directly used without any transformation.[9] Time domain algorithms are generally classified as classical (least square (LS), Extended Kalman filter and maximum likelihood method etc) or non-classical methods (neural network and evolutionary algorithm).[10] Goa and Lu [11] investigated a method for the time-domain analysis of noise-contaminated vibration signals for global structural damage diagnosis. The proposed method is based on conjunction with time-domain analysis approach based on Auto-Regressive with exogenous input (ARX) and Kalman filter technique. Kang et al. [8] presented a time domain algorithm to estimate the stiffness and damping parameters of the structure using measured acceleration. The proposed SI algorithm estimates structural parameters through the minimization of an error function defined by the time integral of the least-squared error between the measured and the calculated accelerations. Koh et al. [12] applied substructural and progressive identification method in time domain using interface measurements. It will reduce the number of unknown structural parameters to be identified, reducing the computational effort. In general the classical methods work well, but they do face problems in getting stuck in local optima, sensitivity to noise, and requirement to give an initial value by the user. These are avoided when using non-classical optimization methods.[12]

A non-classical method is usually based on heuristic concepts such as Genetic Algorithm (GA) based on Evolutionary principle or Particle Swarm Optimization (PSO) on behavioural principle. Varghese and Shankar [13] used PSO for SI based on the combined objective function of power flow balance and acceleration matching. It is based on weighted aggregation approach and identified the crack location and depth with better accuracy than the conventional acceleration matching methods. Nandakumar and Shankar [9] proposed damped transfer matrix and state vectors for structural parameter estimation using PSO. Koh et al. [10] proposed hybrid of GA and compatible local search operator for the identification of structural system with 52 unknown parameters. Zhenghao Ding et al [14]identified multiple cracks using improved artificial bee colony algorithm (I-ABC) based on objective function which consists of limited number of measured natural frequencies. Hu et al. [15] identified cracks using homotopy continuation algorithm based on using several measured acceleration responses. Lia et al. [16] used a derivative of eigen-parameters-based model updating approach for the damage identification of multiple cracks.

Cracks in a structural member reduce the local stiffness which affects vibrational response. Gounaris et al. [17] identified the crack depth and location of a Timoshenko beam FE of six degrees of freedom per node, using a coupled response measurement. However, it is suitable for only single edge non-propagating crack present in the beam and it requires many nomograms for enhancing the accuracy and utility of the identification procedure. Krawczuk et al. [18] developed single edge cracked beam FE which is based on Elasto-plastic fracture mechanics theory. The crack is transverse open edge and is found that crack tip plasticity increases the predictive capability of the model and also inertia matrices are same as for non-cracked element. Viola et al. [19] formulated shape functions of a cracked Timenshenko beam element and used them for identification of the magnitude and location of the crack based on experimentally measured frequencies and mode shapes. Viola et al. [20] investigated the effect of cracks in the stiffness matrix, neglecting the crack closure effect by assuming an open crack.

Sumant and Maiti [21] predicted crack location and depth using a experimentally measured voltage drop between healthy and damaged beam with pair of PZT patches subjected to static loading. Here, the identification is only significant when the distance between a crack and a patch along the edge is approximately less than the depth of the beam. Two stage identification is proposed by Fukunaga et al. [22] using a limited number of PZT sensor in time and frequency domain. It identified crack depths up to 10%, but it strongly depends on reliable modal data of undamaged structure.

The high reliability, sensitivity and electromechanical coupling property of PZT has gained significant attention for potential application as sensors for monitoring and as actuators for controlling the response of structures. Shen [23] proposed a one-dimensional mathematical model for determining the mechanical responses of beams with PZT as actuator and as a sensor. Bendary and Raid [24] formulated one-dimensional integrated beam element using hermit cubic and Lagrangian interpolation function and carried out static and dynamic analysis. Zemcik and Sadilek [25] developed one-dimensional hybrid PZT element based on the Euler Bernoulli beam using bilinear Lagrangian interpolation polynomial for electric potential and carried out modal analysis. Later, the same element is used for frequency response analysis and the result is experimentally verified by Sadilek and Zemcik [26]. Sulbhewar and Raveendranath [27] formulated one-dimensional Euler Bernoulli beam with PZT structure using coupled field polynomial which is independent of material configuration of the piezoelectric beam cross section.

In most of the literature presented accelerometer signal is used as inverse quantity for SI. However, compared to accelerometer, PZT patches have the advantages of low cost, light weight, wide dynamic range. The most notable feature of PZT patches is their miniaturized appearance and ability to be implanted in civil structure for in situ health monitoring of structure. Moreover, because of their light weight the PZT patches are not likely to alter the dynamic properties and also output voltage is less contaminated with noise signal. To the best of the author’s knowledge, this is the first attempt to use voltage output from one-dimensional model PZT patch as an inverse problem for the SI with crack. The theory is validated with numerical and experimental examples in the section given below.

2. Constitutive equations of PZT

Piezo Electric Materials are used to transform the mechanical displacement into an electrical field (voltage potential) in which case the piezoelectric material acts as a sensor(direct effect), and its converse effect acts as an actuator. The constitutive equations for the transversely isotropic piezo electric medium which define the interaction between the stress ($\mathit{\sigma }$), strain ($\mathit{\epsilon }$), electric displacement (D) and electric field (E) in the form [28](1) $$\begin{array}{cc}\hfill {\mathit{\sigma }}_j& =C_{jk}{\mathit{\epsilon }}_k-e_{jm}{E}_m\hfill \end{array}$$(1) (2) $$\begin{array}{cc}\hfill D_l& =e_{lj}{\mathit{\epsilon }}_k+{\mathit{ϵ}}_{lm}{E}_m\hfill \end{array}$$(2)

where $C_{jk}$, ${\mathit{ϵ}}_{lm}$ and $e_{lj}$$(j,k=1,\dots ,6$ and $l,m=1,\dots ,3)$ are elastic, dielectric and piezo electric coupling coefficient, respectively. The transversely poled piezoelectric material is bonded/embedded in host structure. The complete coupled three-dimensional constitutive equation of a piezo electric material with principal material axis (x, y and z) can be written as follows [27](3) $$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathit{\sigma }}_x\\ {\mathit{\sigma }}_y\\ {\mathit{\sigma }}_z\\ {\mathit{\tau }}_{yz}\\ {\mathit{\tau }}_{xz}\\ {\mathit{\tau }}_{xy}\\ D_x\\ D_y\\ D_z\end{array}\right]=\left[\begin{array}{ccccccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0& 0& 0& -e_{31}\\ C_{12}& C_{22}& C_{23}& 0& 0& 0& 0& 0& -e_{32}\\ C_{13}& C_{23}& C_{33}& 0& 0& 0& 0& 0& -e_{33}\\ 0& 0& 0& C_{44}& 0& 0& 0& -e_{24}& 0\\ 0& 0& 0& 0& C_{55}& 0& -e_{15}& 0& 0\\ 0& 0& 0& 0& 0& C_{66}& 0& 0& 0\\ 0& 0& 0& 0& e_{15}& 0& {\mathit{ϵ}}_{11}& 0& e_{32}\\ 0& 0& 0& e_{24}& 0& 0& 0& {\mathit{ϵ}}_{22}& 0\\ e_{31}& e_{32}& e_{33}& 0& 0& 0& 0& 0& {\mathit{ϵ}}_{33}\\ \end{array}\right]\left[\begin{array}{c}{\mathit{\epsilon }}_x\\ {\mathit{\epsilon }}_y\\ {\mathit{\epsilon }}_z\\ {\mathit{\gamma }}_{yz}\\ {\mathit{\gamma }}_{xz}\\ {\mathit{\gamma }}_{xy}\\ E_x\\ E_y\\ E_z\end{array}\right]\end{array}$$(3)

where $\mathit{\tau }$ and $\mathit{\gamma }$ are shear stress and shear strain, respectively

Figure 1. Simple three dimensional beam.

3. One dimensional beam geometry

The simple three-dimensional beam with axis is shown in Figure 1. The Euler–Bernoulli beam theory is applied for one-dimensional beam with piezo electric patch, which neglect the shear effect. The model is assumed to be plane stress and width in the y direction has zero stress. Therefore, it is possible to set ${\mathit{\sigma }}_y={\mathit{\sigma }}_z={\mathit{\tau }}_{xy}={\mathit{\tau }}_{xz}={\mathit{\tau }}_{yz}={\mathit{\gamma }}_{xy}={\mathit{\gamma }}_{xz}={\mathit{\gamma }}_{yz}=0$ while ${\mathit{\epsilon }}_y\ne 0$; ${\mathit{\epsilon }}_z\ne 0$ [27]. The polarization axis z is aligned with the thickness direction of the beam, thus only $D_z$ is taken and thus for electric field $E_x=E_y=0$. Applying these condition in Equation (3(3) $$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathit{\sigma }}_x\\ {\mathit{\sigma }}_y\\ {\mathit{\sigma }}_z\\ {\mathit{\tau }}_{yz}\\ {\mathit{\tau }}_{xz}\\ {\mathit{\tau }}_{xy}\\ D_x\\ D_y\\ D_z\end{array}\right]=\left[\begin{array}{ccccccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0& 0& 0& -e_{31}\\ C_{12}& C_{22}& C_{23}& 0& 0& 0& 0& 0& -e_{32}\\ C_{13}& C_{23}& C_{33}& 0& 0& 0& 0& 0& -e_{33}\\ 0& 0& 0& C_{44}& 0& 0& 0& -e_{24}& 0\\ 0& 0& 0& 0& C_{55}& 0& -e_{15}& 0& 0\\ 0& 0& 0& 0& 0& C_{66}& 0& 0& 0\\ 0& 0& 0& 0& e_{15}& 0& {\mathit{ϵ}}_{11}& 0& e_{32}\\ 0& 0& 0& e_{24}& 0& 0& 0& {\mathit{ϵ}}_{22}& 0\\ e_{31}& e_{32}& e_{33}& 0& 0& 0& 0& 0& {\mathit{ϵ}}_{33}\\ \end{array}\right]\left[\begin{array}{c}{\mathit{\epsilon }}_x\\ {\mathit{\epsilon }}_y\\ {\mathit{\epsilon }}_z\\ {\mathit{\gamma }}_{yz}\\ {\mathit{\gamma }}_{xz}\\ {\mathit{\gamma }}_{xy}\\ E_x\\ E_y\\ E_z\end{array}\right]\end{array}$$(3) ) and it is reduced in the form(4) $$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathit{\sigma }}_x\\ D_z\end{array}\right]=\left[\begin{array}{cc}\widehat{C}& -\widehat{e}\\ \widehat{e}& \widehat{\mathit{ϵ}}\\ \end{array}\right]\left[\begin{array}{c}{\mathit{\epsilon }}_x\\ E_z\end{array}\right]\end{array}$$(4)

where $\widehat{C}=Q_{11}-\frac{{(Q_{12})}^2}{Q_{22}}$ and $Q_{ij}=C_{ij}-(\frac{C_{i3}{C}_{j3}}{C_{33}})(i,j=1,2)$; $\widehat{e}=\overline{e_{31}}-\overline{e_{32}}(\frac{Q_{12}}{Q_{22}})$ and $\overline{e_{3i}}=e_{3i}-e_{33}(\frac{C_{13}}{C_{33}})(i=1,2)$; $\widehat{\mathit{ϵ}}={\overline{\mathit{ϵ}}}_{33}+\frac{{({\overline{e}}_{32})}^2}{Q_{22}}$ and ${\overline{\mathit{ϵ}}}_{33}={\mathit{ϵ}}_{33}+\frac{{(e_{32})}^2}{C_{33}}$.[29] Here, $\widehat{C}$, $\widehat{\mathit{ϵ}}$ and $\widehat{e}$ are reduced elastic, dielectric and piezo electric coupling coefficient, respectively, and these calculated values are shown in Table 1. This reduced property is used for further numerical study of one-dimensional beam with PZT patch in MATLAB®.

Table 1. Material properties of the beam.

Material	Properties
Aluminium [30]	$\mathrm{E}=71$ GPa, $\mathit{\nu }=0.3,\mathit{\rho }=2210{\mathrm{kg}\mathrm{m}}^{-3}$
PZT 5H [29]	$C_{11}=C_{22}=126$ GPa ,$C_{12}=79.5$ Gpa
	$C_{13}=C_{23}=84.1$ GPa ,$C_{33}=117$ Gpa
	$C_{44}=C_{55}=23$ GPa ,$C_{66}=23.25$ Gpa
	$e_{31}=e_{32}=-6.5\mathrm{C}{\mathrm{m}}^{-2},e_{33}=23.3\mathrm{C}{\mathrm{m}}^{-2},e_{15}=e_{24}=17\mathrm{C}{\mathrm{m}}^{-2}$
	${\mathit{ϵ}}_{11}={\mathit{ϵ}}_{22}=1.503\times {10}^{-8}\mathrm{F}{\mathrm{m}}^{-1},{\mathit{ϵ}}_{33}=1.3\times {10}^{-8}\mathrm{F}{\mathrm{m}}^{-1}$
	$\mathit{\rho }=7500\mathrm{kg}{\mathrm{m}}^{-3}$
Reduced Properties	$\widehat{C}=60.013$ GPa, $\widehat{e}=-16.4921\mathrm{C}{\mathrm{m}}^{-2}$
	$\widehat{\mathit{ϵ}}=2.5885\times {10}^{-8}\mathrm{F}{\mathrm{m}}^{-1}$

Figure 2. One-dimensional beam element: PZT and supporting structure sharing common nodes.

4. Analytical and FEM formulation

One-dimensional beam formulation has been presented in Refs. [26,27], but not used for time domain dynamic response with damping and inverse problem. The FE formulation of the structure is briefly explained here. The beam element is based on the Euler–Bernoulli theory and element has two nodes. The two independent polynomials are used for interpolation of mechanical and electrical field variable. First hermit cubic polynomial is used for the interpolation of mechanical quantities of vertical displacement (w) and rotation ($\mathit{\theta }$) as shown in Figure 2. Let the vertical displacement (w) is approximated across the length as(5) $$\begin{array}{c}\hfill w(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3=[N_w]\{w^e\}\end{array}$$(5)

where $N_w$ is the shape interpolation function, and structural nodal degree of freedom per element is arranged as(6) $$\begin{array}{c}\hfill w^e={[w_1,{\mathit{\theta }}_1,w_2,{\mathit{\theta }}_2]}^T\end{array}$$(6)

Thus Equation (5(5) $$\begin{array}{c}\hfill w(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3=[N_w]\{w^e\}\end{array}$$(5) ) can be written as(7) $$\begin{array}{c}\hfill w(x)=N_1{w}_1+N_2{\mathit{\theta }}_1+N_3{w}_2+N_4{\mathit{\theta }}_2=\sum _{j=1}^4{N}_j{w}_j\end{array}$$(7)

where(8) $$\begin{array}{cc}\hfill N_1& =1-3x^2/{l_e}^2+2x^3/{l_e}^3,N_2=x-2x^2/{l_e}^2+x^3/{l_e}^3\hfill \\ \hfill N_3& =3x^2/{l_e}^2-2x^3/{l_e}^3,N_4=x^3/{l_e}^2-x^2/l_e\hfill \end{array}$$(8)

The axial displacement can be written as(9) $$\begin{array}{c}\hfill u(x,z)=-z\mathit{\theta }(x)=-z\frac{\mathit{\partial }w}{\mathit{\partial }x}\end{array}$$(9)

Thus bending strain is(10) $$\begin{array}{c}\hfill \mathit{\epsilon }(x,z)=\frac{\mathit{\partial }u}{\mathit{\partial }x}=-z\frac{{\mathit{\partial }}^{2}w}{\mathit{\partial }{x}^2}=[B_w]\{w^e\}\end{array}$$(10)

where $[B_w]$ is the strain–displacement matrix consisting of derivatives of shape functions.

The electric potential is $\mathit{\varphi }(x,z)$ considered as a function of the thickness and the length of the beam. Hence, let langrangian bilinear function is estimated for the interpolation as(11) $$\begin{array}{c}\hfill \mathit{\varphi }(x,z)=a_4+a_{5}x+a_{6}z+a_{7}xz=[{\mathrm{\Phi }}_{\mathit{\varphi }}]\{{\mathit{\varphi }}^e\}\end{array}$$(11)

As shown in Figure 2, $\mathit{\varphi }$ and $\mathit{\psi }$ are lower and upper electrical potential of PZT surface. The electrical nodal degree of freedom per element can be ordered as(12) $$\begin{array}{c}\hfill \{{\mathit{\varphi }}^e\}={[{\mathit{\varphi }}_1,{\mathit{\psi }}_1,{\mathit{\varphi }}_2,{\mathit{\psi }}_2]}^T\end{array}$$(12)

Thus Equation (11(11) $$\begin{array}{c}\hfill \mathit{\varphi }(x,z)=a_4+a_{5}x+a_{6}z+a_{7}xz=[{\mathrm{\Phi }}_{\mathit{\varphi }}]\{{\mathit{\varphi }}^e\}\end{array}$$(11) ) can be written as(13) $$\begin{array}{c}\hfill \mathit{\varphi }(x,z)={\mathrm{\Phi }}_1{\mathit{\varphi }}_1+{\mathrm{\Phi }}_2{\mathit{\psi }}_1+{\mathrm{\Phi }}_3{\mathit{\varphi }}_2+{\mathrm{\Phi }}_4{\mathit{\psi }}_2=\sum _{j=1}^4{\mathrm{\Phi }}_j{\mathit{\varphi }}_j\end{array}$$(13)

where(14) $$\begin{array}{cc}\hfill {\mathrm{\Phi }}_1& =(l_e-x)(h+2h_p-2z)/2h_p{l}_e,{\mathrm{\Phi }}_2=-((h-2z)(l_e-x))/2hp{l}_e\hfill \\ \hfill {\mathrm{\Phi }}_3& =x(h+2h_p-2z)/2h_p{l}_e,{\mathrm{\Phi }}_4=-x(h-2z)/2h_{p}le\hfill \end{array}$$(14)

The electric field E(x, z) and electric potential $(\mathit{\varphi })$ can be related as follows(15) $$\begin{array}{c}\hfill E(x,z)=\frac{\mathit{\partial }\mathit{\varphi }(x,z)}{\mathit{\partial }z}\end{array}$$(15)

Thus it can be written as follows(16) $$\begin{array}{c}\hfill E(x,z)=[B_{\mathit{\varphi }}]\{{\mathit{\varphi }}^e\}\end{array}$$(16)

where $[B_{\mathit{\varphi }}]$ is the Electrical field-potential matrix consisting of derivatives of shape functions.

Here, homogeneous electrical boundary condition is imposed on the bottom surface of PZT patch to eliminate rigid body modes i.e. Lower surface is grounded with $\mathit{\varphi }=0$ V, while the upper surface was left open.

5. Evaluation of elemental matrices

The dynamic potential energy density G for the piezo electric medium can be written as [31](17) $$\begin{array}{c}\hfill G=\frac{1}{2}{\mathit{\epsilon }}^{T}c\mathit{\epsilon }-\frac{1}{2}{E}^T\mathit{ϵ}E-{\mathit{\epsilon }}^{T}eE\end{array}$$(17)

It includes contribution from mechanical and electrostatic energy and kinetic energy density is(18) $$\begin{array}{c}\hfill K=\frac{1}{2}\mathit{\rho }{\dot{w}}^2\end{array}$$(18)

Then apply variational principle to potential, kinetic energy and external forces (mechanical and electrical loading). It must be satisfied for any arbitrary variation of the displacements and electrical potentials and thus equation of motion of elemental matrices can be represented as [32](19) $$\begin{array}{c}\hfill \left[\begin{array}{cc}[M_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\ddot{\{w^e\}}\\ \ddot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[K_{ww}^e]& [K_{w\mathit{\varphi }}^e]\\ {[K_{w\mathit{\varphi }}^e]}^T& [K_{\mathit{\varphi }\mathit{\varphi }}^e]\end{array}\right]\left\{\begin{array}{c}\{w^e\}\\ \{{\mathit{\varphi }}^e\}\end{array}\right\}=\left\{\begin{array}{c}\{F^e\}\\ \{Q^e\}\end{array}\right\}\end{array}$$(19)

The Equation (19(19) $$\begin{array}{c}\hfill \left[\begin{array}{cc}[M_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\ddot{\{w^e\}}\\ \ddot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[K_{ww}^e]& [K_{w\mathit{\varphi }}^e]\\ {[K_{w\mathit{\varphi }}^e]}^T& [K_{\mathit{\varphi }\mathit{\varphi }}^e]\end{array}\right]\left\{\begin{array}{c}\{w^e\}\\ \{{\mathit{\varphi }}^e\}\end{array}\right\}=\left\{\begin{array}{c}\{F^e\}\\ \{Q^e\}\end{array}\right\}\end{array}$$(19) ) is the elemental equilibrium in the discretized form, where $[M_{ww}^e]$ is the mass matrix, $[K_{ww}^e]$ is the stiffness matrix corresponding to the mechanical degree of freedom, $[K_{w\mathit{\varphi }}^e]$ is the stiffness matrix due to electromechanical coupling, $[K_{\mathit{\varphi }\mathit{\varphi }}^e]$ is the stiffness matrix due to the electrical degrees of freedom alone, $\{F^e\}$ is the mechanical load vector and $\{Q^e\}$ is the electrical charge load vector. The size of each element matrices are $[4\times 4]$ and are given by(20) $$\begin{array}{c}\hfill [M_{ww}^e]={\int }_V{[N_w]}^T\mathit{\rho }[N_w]\mathrm{d}V\\ \hfill [K_{ww}^e]={\int }_V{[B_w]}^T\widehat{C}[B_w]\mathrm{d}V\\ \hfill [K_{w\mathit{\varphi }}^e]={\int }_V{[B_w]}^T\widehat{e}[B_{\mathit{\varphi }}]\mathrm{d}V\\ \hfill [K_{\mathit{\varphi }\mathit{\varphi }}^e]={\int }_V{[B_{\mathit{\varphi }}]}^T\widehat{\mathit{ϵ}}[B_{\mathit{\varphi }}]\mathrm{d}V\end{array}$$(20)

In general, all structures are lightly damped. Thus, adding an artificial linear viscous damping to Equation (19(19) $$\begin{array}{c}\hfill \left[\begin{array}{cc}[M_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\ddot{\{w^e\}}\\ \ddot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[K_{ww}^e]& [K_{w\mathit{\varphi }}^e]\\ {[K_{w\mathit{\varphi }}^e]}^T& [K_{\mathit{\varphi }\mathit{\varphi }}^e]\end{array}\right]\left\{\begin{array}{c}\{w^e\}\\ \{{\mathit{\varphi }}^e\}\end{array}\right\}=\left\{\begin{array}{c}\{F^e\}\\ \{Q^e\}\end{array}\right\}\end{array}$$(19) ) [32](21) $$\begin{array}{c}\hfill \left[\begin{array}{cc}[M_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\ddot{\{w^e\}}\\ \ddot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[C_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\dot{\{w^e\}}\\ \dot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[K_{ww}^e]& [K_{w\mathit{\varphi }}^e]\\ {[K_{w\mathit{\varphi }}^e]}^T& [K_{\mathit{\varphi }\mathit{\varphi }}^e]\end{array}\right]\left\{\begin{array}{c}\{w^e\}\\ \{{\mathit{\varphi }}^e\}\end{array}\right\}=\left\{\begin{array}{c}\{F^e\}\\ \{Q^e\}\end{array}\right\}\end{array}$$(21)

where the damping matrix $[C_{ww}^e]$ is defined as a proportional damping i.e.(22) $$\begin{array}{c}\hfill [C_{ww}^e]=a[M_{ww}^e]+b[K_{ww}^e]\end{array}$$(22)

where a and b are Rayleigh’s damping coefficients.

The matrices of equations are then assembled to obtain the global dynamic system equation. In order to solve this, the global matrix Equation (21(21) $$\begin{array}{c}\hfill \left[\begin{array}{cc}[M_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\ddot{\{w^e\}}\\ \ddot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[C_{ww}^e]& [0]\\ [0]& [0]\end{array}\right]\left\{\begin{array}{c}\dot{\{w^e\}}\\ \dot{\{{\mathit{\varphi }}^e\}}\end{array}\right\}+\left[\begin{array}{cc}[K_{ww}^e]& [K_{w\mathit{\varphi }}^e]\\ {[K_{w\mathit{\varphi }}^e]}^T& [K_{\mathit{\varphi }\mathit{\varphi }}^e]\end{array}\right]\left\{\begin{array}{c}\{w^e\}\\ \{{\mathit{\varphi }}^e\}\end{array}\right\}=\left\{\begin{array}{c}\{F^e\}\\ \{Q^e\}\end{array}\right\}\end{array}$$(21) ) can be expanded and written as:(23) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+[C_{ww}]\dot{\{w\}}+[K_{ww}]\{w\}+[K_{w\mathit{\varphi }}]\{\mathit{\varphi }\}=\{F\}\\ \hfill {[K_{w\mathit{\varphi }}]}^T\{w\}+[K_{\mathit{\varphi }\mathit{\varphi }}]\{\mathit{\varphi }\}=\{Q\}\end{array}$$(23)

The second part of Equation (23(23) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+[C_{ww}]\dot{\{w\}}+[K_{ww}]\{w\}+[K_{w\mathit{\varphi }}]\{\mathit{\varphi }\}=\{F\}\\ \hfill {[K_{w\mathit{\varphi }}]}^T\{w\}+[K_{\mathit{\varphi }\mathit{\varphi }}]\{\mathit{\varphi }\}=\{Q\}\end{array}$$(23) ) is simplified as(24) $$\begin{array}{c}\hfill \{\mathit{\varphi }\}={[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}\{Q\}-{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}{[K_{w\mathit{\varphi }}]}^T\{w\}\end{array}$$(24)

The above equation is substituted in the first part of the (23(23) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+[C_{ww}]\dot{\{w\}}+[K_{ww}]\{w\}+[K_{w\mathit{\varphi }}]\{\mathit{\varphi }\}=\{F\}\\ \hfill {[K_{w\mathit{\varphi }}]}^T\{w\}+[K_{\mathit{\varphi }\mathit{\varphi }}]\{\mathit{\varphi }\}=\{Q\}\end{array}$$(23) ) and is written as(25) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+[C_{ww}]\dot{\{w\}}+{[K_{ww}]}^{\ast }\{w\}={\{F\}}^{\ast }\end{array}$$(25)

where(26) $$\begin{array}{cc}\hfill {[K_{ww}]}^{\ast }& =[K_{ww}]-[K_{w\mathit{\varphi }}][K_{\mathit{\varphi }\mathit{\varphi }}]{[K_{w\mathit{\varphi }}]}^T,\hfill \\ \hfill {\{F\}}^{\ast }& =\{F\}-[K_{w\mathit{\varphi }}]{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}\{Q\}\hfill \end{array}$$(26)

For the sensor problem, displacement histories are obtained by solving Equation (25(25) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+[C_{ww}]\dot{\{w\}}+{[K_{ww}]}^{\ast }\{w\}={\{F\}}^{\ast }\end{array}$$(25) ) and substituted in (24(24) $$\begin{array}{c}\hfill \{\mathit{\varphi }\}={[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}\{Q\}-{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}{[K_{w\mathit{\varphi }}]}^T\{w\}\end{array}$$(24) ), and voltage vector across the sensor patch is obtained. In order to perform modal analysis, the above system is statically condensed and equation is modified to(27) $$\begin{array}{c}\hfill [M_{ww}]\ddot{\{w\}}+([K_{ww}]-[K_{w\mathit{\varphi }}]{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}{[K_{w\mathit{\varphi }}]}^T)w=0\end{array}$$(27)

and the problem is reduced to the eigen value of the matrix(28) $$\begin{array}{c}\hfill A={[M_{ww}]}^{-1}([K_{ww}]-[K_{w\mathit{\varphi }}]{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}{[K_{w\mathit{\varphi }}]}^T)\end{array}$$(28)

6. Numerical validation of one dimensional FE model

A computer code for the above-described hybrid beam element in MATLAB® environment has been developed to study the dynamic analysis of PZT sensor bonded to the top surface of cantilever beam as shown in Figure 3. The dimension of the rectangular beam is taken from Sandesh and Shankar [30] and PZT sensor dimensions are ($50\times 25.4\times 1{\mathrm{mm}}^3$). The material properties of host structure (Aluminium) and PZT sensor (reduced properties) are used as shown in Table 1. As no commercial one-dimensional PZT element is available, the present formulation is validated using ANSYS® 3D simulation.

For Ansys® 3D simulation a mesh of SOLID 186 element is used for the host structure and a mesh of SOLID 226 element is used for sensor(PZT 5H) patch.

Modal Analysis is carried out first for the given model and first three natural frequencies are tabulated in Table 2. As shown from table, natural frequency from present formulation is in good agreement with ANSYS® 3D simulation. The convergence test has been done for Ansys model using smaller elements.

Figure 3. Cantilever beam with PZT patch.

Table 2. Natural frequency of cantilver beam structure.

		Natural frequency in Hz
Formulation	First	Second	Third
Matlab	23.27	145.74	407.05
Ansys	23.24	145.12	402.72

Figure 4. Comparison of (a) Voltage and (b) Tip Displacement time histories.

In time domain analysis, numerically simulated response histories of all DOFs of the model are calculated in terms of displacement, velocity, acceleration and voltage using Newmark-Beta method. An impulsive force of 2 N in time period of 0.01 s in time step of 0.001 is applied in the free end of a cantilever beam. The effect of the damping was accounted by Rayleigh’s damping with modal damping ratio of 1% at its first two modes. The tip deflection time histories at point B in vertical direction and maximum voltage potential (near the clamped end) time history at PZT patch A are plotted in Figure 4. As seen from the plot, the present MATLAB® code gives precise predictions as given by the ANSYS® simulation. Comparable validation is also done in the frequency domain and which are shown in Figure 5. From this plot, two well-defined peaks are identified as first two modal frequencies. The result compares well with ANSYS® 3D simulation.

7. PSO algorithm

A heuristic optimization technique referred to as PSO is used here which mimics the social behaviour of swarms. It was first proposed by James and Ebehart [33]. It imitates the social behaviour of a swarm of birds. Each bird will tend to follow the general swarm direction in search of the target (food), but it has a component of its own intelligence and memory (i.e. local search) which influences its action. Each bird is visualized as a ‘particle’ which approaches the target (i.e. the global optima) with a ‘velocity’. The number of particles (i.e. population) and their initial random positions are specified. As the particles progress to the global optima through many generations, their current position is updated using two parameters $G_{\text{best}}$ which represent the historically best co-ordinate of all the particles in the population and $P_{\text{best},i}$ the historically best co-ordinate of the ith particle. The equations giving the velocity v and position x for the ith particle in the $k+1$ generation are given by,(29) $$\begin{array}{cc}\hfill v_i(k+1)& =\mathit{\phi }(k)v_i(k)+{\mathit{\alpha }}_1[{\mathit{\gamma }}_{1i}(P_{\text{best},i}-x_i(k))]+{\mathit{\alpha }}_2[{\mathit{\gamma }}_{2i}[(G_{\text{best}}-x_i(k))]\hfill \end{array}$$(29) (30) $$\begin{array}{cc}\hfill x_i(k+1)& =x_i(k)+v_i(k+1)\hfill \end{array}$$(30)

where i is the particle index, k the discrete time index, v the velocity of the ith particle, x the position of the ith particle/present solution. Here ${\mathit{\gamma }}_1$ and ${\mathit{\gamma }}_2$ represent two random numbers between 0 and 1, $\mathit{\phi }$ is an inertia term uniformly decreasing from 0.9 to 0.4 with generations and ${\mathit{\alpha }}_1$ and ${\mathit{\alpha }}_2$ are two acceleration constants set to two.[34] Several studies have pointed out the superiority of PSO algorithm over the more conventional heuristic algorithms such as GA for inverse problem applications.[35,36]

Figure 5. Comparison of (a) Voltage and (b) Tip Displacement in frequency domain.

8. Structural parameter identification using PZT patch

The time domain-based approach is used and voltage history is used as the main response quantity in the identification procedure. A few experimentally measured voltage potential responses ${\mathit{\varphi }}^m$ are measured from PZT patches. The estimated voltage potential ${\mathit{\varphi }}^e$ is obtained from mathematical model using Equation (24(24) $$\begin{array}{c}\hfill \{\mathit{\varphi }\}={[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}\{Q\}-{[K_{\mathit{\varphi }\mathit{\varphi }}]}^{-1}{[K_{w\mathit{\varphi }}]}^T\{w\}\end{array}$$(24) ). For exact identification, ${\mathit{\varphi }}^e$ has to match with experimentally measured responses ${\mathit{\varphi }}^m$. In this method, experimental responses are simulated from a known numerical model and polluted with Gaussian noise of zero mean and a certain standard deviation for practicality. Using PSO algorithm the following fitness (objective) function is minimized, which is the sum of squares of deviations between the measured and estimated voltage.(31) $$\begin{array}{c}\hfill f=\frac{{\sum }_{i=1}^M{\sum }_{j=1}^L{\left|{\mathit{\varphi }}^m(i,j)-{\mathit{\varphi }}^e(i,j)\right|}^2}{ML}\end{array}$$(31)

The superscripts m and e denote measured and estimated responses for fitness evaluation, M is the number of measurement sensors used and L is the number of time steps. Ideally, it must be minimized to zero, but usually it approaches a value close to zero.

9. Numerical results and discussion

In order to validate the proposed method, various numerical simulations are carried out. The purpose of these numerical simulations is to study the effect of patch size on the accuracy of identification. Three cases studied here are (a) stiffness identification of a beam (b) a more complex structure i.e. sub structure of nine member frame and (c) crack identification (magnitude and location) in a beam structure. Experimental measured responses are numerically simulated from a fully defined model in MATLAB® using Newmark’s constant acceleration scheme. The structure is excited by a impulse force and responses are measured at PZT patches. The mass of the structure is assumed to be known a priori and the stiffness parameters are the unknowns to be identified for the first two cases. These are estimated using the inverse formulation with a single objective approach. In order to simulate the effect of noise in experiments, Gaussian random noise level of 5%(standard deviation) and zero mean is added to all the measured signals. The noise is added with simulated responses in such a way that,(32) $$\begin{array}{c}\hfill \text{Noise}=g\times {\text{RMS}}_{\text{signal}}\times NL\end{array}$$(32)

where g is the standard Gaussian variable, ${\text{RMS}}_{\text{signal}}$ is the Root Mean Square (RMS) value of the numerically simulated signal and NL is Noise level.

In this numerical study, three different PZT patch lengths are investigated and are PZT:5, PZT:10 and PZT:20%. They, respectively, represent the length of the patch expressed as percentage of beam length. The width of PZT is same as host structure and a constant thickness of 1 mm is used for this study.

9.1. Fixed-fixed beam structure

The proposed method is first applied to fixed–fixed beam as shown in Figure 6. Model geometry and properties of host structure (Aluminium) are taken from [30]. In this study, two PZT patches are bonded on either side of the structure and the beam is divided into five Euler–Bernoulli beam elements. The Flexural rigidity (EI) of each element is 38.48 $\text{N}{\text{m}}^2$. The first and second natural frequency of modes of vibration are 148.46 and 415.06 Hz respectively. Rayleigh damping with the modal damping ratio of 3% is used for the first two modes of vibration. The beam is excited by providing an impulse force of 2 N at node 4 over a time of 0.01 s in the time step of 0.001.

Figure 6. Fixed-fixed beam with PZT patches.

The unknown parameters such as stiffness and damping are obtained by PSO in the range of $\pm 50\%$ of actual value. In this PSO parameters are set to 100 particles (swarm size) and 200 generations. First, in this study we compare the identification of structural parameters by voltage matching with different lengths of PZT. The identified results are shown in Tables 3 and 4, respectively. From tables, it can be seen that for noise free case Mean Absolute Error (MAE) of PZT-10% i.e. $(50\times 25.4\times 1{\mathrm{mm}}^3$ PZT) is 0.089% and whereas PZT-5% i.e. $(25\times 25.4\times 1{\mathrm{mm}}^3$ PZT) is 0.136%. The percentage increase of error compared with these two cases is only 34.5%. Similarly for 5% noise case, it is 0.111 and 0.167%, respectively, and percentage change of error of these two cases is 33% only. The percentage of error in identifying damping constants $\mathit{\alpha }$ and $\mathit{\beta }$ of PZT-10% are 0.032% and 0.002% for noise free case and is 0.109 and 0.242% for 5% noise, respectively. Similarly for PZT-5% the results are 0.044 and0.442% for no noise and 0.109 and 0.563 for 5% noisy case, respectively.

The objective function convergence study in identifying stiffness and damping coefficient for beam structure without and 5% noise contamination is presented in Figure 7. The Figure shows convergence of voltage matching ($\mathit{\varphi }$) using different lengths of PZT sensor and acceleration matching (A) with accelerometer sensor (placed at node no:3). The method of acceleration matching is discussed in studies on Ref. [10,12,13]. It is clear that as length of PZT increases, convergence is faster and thus identification error can be reduced. In Figure 8(a) and (b) shown, objective function of PZT:20% has minimized better than PZT:5 and PZT:10% and is most accurate in identification. However, it is not required to use such large patch length. PZT:5% error in identification is only 0.167% which is sufficient for practical use. As compared with pure acceleration matching, proposed method even with PZT:5% shows better convergence and thus identifies the parameters with better accuracy.

A numerical study is carried out using 20 elements for the same configuration. The MAE of structural parameter identification using 20 elements with PZT:5% is increased to 5.35 and 8.29%, respectively, for noise free and noisy signals. The reason could be the difficulty in convergence of the optimization algorithm because of the number of unknown variables (element stiffness) which is now increased with more elements. Thus it appears that for a uniform beam, five elements are sufficient for SI.

9.1.1. Comparison of results with other time domain methods

The same problem is compared with acceleration matching using the accelerometer sensor without PZT patches. The results of structural parameter identification are shown in Table 5. The Mean Absolute Error for identifying stiffness for noise free and 5% noise are 0.343 and 0.841%, respectively, for acceleration matching method. The same problem of the cantilever was solved by sub structuring [30] using LS technique with accelerometer. There, MAE of identified stiffness parameters was 4.5% at 3% noise level and damping parameters are very high of 13.72 and 21.24%, respectively, and average computational time was 300 s but identification error by voltage matching with PZT is much lower than the above methods and here, mean computational time is only 60 s. This shows that using PZT as SI shows improved results when compared with other time domain methods.

Table 3. Identified structural parameters for only voltage matching: PZT:10%.

		Identified EI		Identified EI	% of Error
	Exact EI	$\mathrm{N}.m^2$	% of Error	$\text{N}{\text{m}}^2$
Element	${\mathrm{N}\mathrm{m}}^2$	Noise free		5% Noise
1	38.619	38.583	0.093	38.569	0.132
2	38.479	38.539	-0.156	38.508	-0.074
3	38.479	38.469	0.0254	38.531	-0.133
4	38.479	38.442	0.098	38.427	0.137
5	38.619	38.648	-0.075	38.588	0.081
Max absolute error (%)			0.156	0.137
MAE (%)			0.089	0.111
	Exact	Identified		Identified
Damping $\mathit{\alpha }$	41.124	41.138	-0.032	41.169	-0.109
constants $\mathit{\beta }$	$1.702\times {10}^{-5}$	$1.703\times {10}^{-5}$	0-.002	$1.698\times {10}^{-5}$	0.242

Table 4. Identified structural parameters for only voltage matching:PZT:5%.

		Identified EI		Identified EI	% of Error
	Exact EI	$\mathrm{N}{m}^2$	% of Error	$\text{N}{\text{m}}^2$
Element	${\mathrm{N}\mathrm{m}}^2$	Noise free		5% Noise
1	38.619	38.637	-0.046	38.529	0.231
2	38.479	38.428	0.134	38.514	$-0.09$
3	38.479	38.465	0.038	38.442	0.098
4	38.479	38.567	-0.226	38.538	0.152
5	38.619	38.528	0.236	38.517	0.263
Max absolute error (%)			0.236	0.263
MAE (%)			0.136	0.167
	Exact	Identified		Identified
Damping $\mathit{\alpha }$	41.124	41.106	0.044	41.396	-0.656
constants $\mathit{\beta }$	$1.702\times {10}^{-5}$	$1.69\times {10}^{-5}$	-0.442	$1.693\times {10}^{-5}$	0.563

Table 5. Identified structural parameters for only acceleration matching.

		Identified EI		Identified EI	% of Error
	Exact EI	${\mathrm{N}\mathrm{m}}^2$	% of error	$\text{N}{\text{m}}^2$
Element	${\mathrm{N}\mathrm{m}}^2$	Noise free		5% Noise
1	38.479	38.742	$-$0.678	38.024	1.197
2	38.479	38.471	0.021	38.8804	$-$1.0319
3	38.479	38.376	0.268	38.8704	$-$1.0065
4	38.479	38.656	$-$0.4579	38.1744	0.7997
5	38.479	38.590	$-$0.2879	38.544	$-$0.1671
Max absolute error (%)			0.678	1.197
MAE (%)			0.343	0.841
	Exact	Identified		Identified
Damping $\mathit{\alpha }$	41.124	40.954	0.415	41.496	$-$0.894
constants $\mathit{\beta }$	$1.702\times {10}^{-5}$	$1.710\times {10}^{-5}$	$-0.377$	$1.713\times {10}^{-5}$	$-$0.639

Figure 7. Convergence plot for (a) Noise-free and (b) 5% Noise case.

9.2. Sub-structure of a frame structure

In this example, a steel frame structure consisting of nine members is fixed at two supports as shown in Figure 8 as per Nandakumar and Shankar [9]. Each member has flexural rigidity (EI) of 43.2 N.m${}^2$ and cross section of $12\times 6$mm. The fundamental natural frequency of the system is 12 Hz. Here, two members of structure as shown in Figure 8(a) and (b) were identified by the different PZT patch length as attempted in the previous example. Each frame is divided into five one-dimensional FEs. The structure is excited by an impulse response of 10 N at the midpoint of the member 6. Rayleigh damping with the modal damping ratio of 5% is used. The identified results of frame No:4 and frame No:5 using PZT:10% i.e. $(100\times 12\times 1{\mathrm{mm}}^3$ PZT) and PZT:5% i.e. $(50\times 12\times 1{\mathrm{mm}}^3$ PZT) is shown in Tables 6 and 7, respectively. The MAE of member No:4, member No:5 without noise using PZT:10 and PZT:5% are 0.043, 0.053 and 0.109, 0.129%, respectively. The percentage increase of identified error compared with these two PZT patch lengths is 23.25 and 18.34%, respectively, for frame No:4 and frame No:5. Similarly for 5% noise case, these are 0.378, 0.436 (frame No:4) and 0.765, 0.959% (frame No:5), respectively, for PZT:10 and PZT:5%. The percentage increase of error compared with these two cases is 15.34 and 25.35%, respectively. The maximum percentage error of damping constant with 5% noise of frame No:5 are 4.23 and 5.56%, respectively, for patch length PZT:10 and PZT:5%. The frame No:4, it is 0.473 and 0.589 respectively, for PZT:10 and PZT:5%.

Convergence of fitness(objective) functions of different PZT length for the SI of frame N:4 (Noise free and 5% Noise) are presented in Figure 9. From the plot, it can be seen that PZT:20% has been minimized very well. However, it is not necessary to use in practice.

9.2.1. Comparison of results with other time domain methods

It may be noted that Nandakumar and Shankar [9] have solved the same problem using transfer matrix in time domain method and identified the stiffness parameters with MAE of 0.47 and 1.01% in measurement with 3% noise at all measured response. Also damping constant was identified with maximum percentage error of 7.14% at the same noise level and mean computational time was 50 s There, the error is slightly more than the proposed approach and also implementation of current PZT method is much simpler than that of combined accelerometer and strain gauge for transfer matrices. The parameters of the same frame were identified by using the complete structural identification algorithm proposed by Koh et al [12] with same input conditions and by same PSO parameters. The mean absolute error was 7.66% with complete measurement at all nodes without noise and required a computational time of 58 min for convergence. It is also compared with Yun and Babng [37] using artificial neural network in frequency domain for structural identification of a very similar frame structure. The substructure of frame was identified with 5% noise with an average error of 8% and also it requires more computational time for training the data. It is shown that structural identification with PZT gives improved results and here, mean computational time is only 90 s.

10. Crack location and depth identification of a structure

The FE formulation of cracked model is taken from Viola et al. [18,19] and briefly presented here. However, it was not used for inverse problem, or combined with One-Dimensional PZT patch. The cracked beam element and corresponding FE model are shown in Figures 10 and 11, respectively. Let $l_e$ be the length of the element, $l_1$ be the location of the crack from its left end node, a be the crack depth measured from the top of the crack section. The cross-sectional dimensions of the beam are $b\times h$.

Table 6. Percentage of error in identified parameters of frame structure: PZT:10%.

		% of error
		Frame no: 4		Frame no: 5
	Exact EI	Noise	With 5%	Noise	With 5%
Element	N m${}^2$	free	Noise	free	Noise
1	43.265	0.029	0.004	$-0.157$	$-0.4663$
2	43.2	$-0.024$	0.081	0.068	$-0.623$
3	43.2	0.161	$-0.573$	$-0.080$	1.272
4	43.2	$-3\mathrm{e}-3$	0.221	0.144	$-1.071$
5	43.265	$3.21\mathrm{e}-5$	$-1.01$	$-0.096$	$-0.332$
Mean absolute error(%)		0.043	0.378	0.109	0.765
Damping $\mathit{\alpha }$	3.944	$-0.019$	0.473	$-0.032$	$-4.23$
constants $\mathit{\beta }$	$6.32\mathrm{e}-4$	0.008	$-0.206$	0.025	2.793

Table 7. Percentage of error in identified parameters of frame structure: PZT:5%.

		% of error
		Frame no: 4		Frame no: 5
	Exact EI	Noise	With 5%	Noise	With 5%
Element	N m^2	free	Noise	free	Noise
1	43.265	0.037	$-0.445$	$-0.179$	$-0.489$
2	43.2	0.116	$-0.773$	$-0.163$	$-0.668$
3	43.2	$-0.075$	0.166	$-0.08$	1.273
4	43.2	$-0.0295$	0.229	0.121	$-1.746$
5	43.265	$-0.009$	0.564	$-0.098$	0.621
Mean absolute error(%)		0.053	0.436	0.129	0.959
Damping $\mathit{\alpha }$	3.944	$-0.03$	0.589	$-0.04$	$-5.56$
constants $\mathit{\beta }$	6.32e-4	0.007	$-0.316$	0.03	3.493

Figure 8. Nine member frame structure.

Figure 9. Convergence plot for (a) noise-free and (b) 5% noise case: frame member 4.

Figure 10. Cracked beam element.

Figure 11. Equivalent model of cracked beam element.

It is derived based on elasto-plastic fracture mechanics and is as Euler–Bernoulli beam element having three segments.[18] They are two intact segments and crack with negligible length and mass. The element has two nodes and with two degrees of freedom (transverse displacement and rotation) in each node. Since it is discontinued at the crack plane, two different polynomials are used for the interpolation of the field variable as follows [18](33) $$\begin{array}{cc}\hfill v_1(x)& =a_1+a_{2}x+a_3{x}^2+a_4{x}^3;0\le x\le l_1\hfill \\ \hfill {\mathit{\theta }}_1(x)& =\frac{\mathrm{d}{v}_1(x)}{\mathrm{d}x}=a_2+2a_{3}x+3a_4{x}^2;0\le x\le l_1\hfill \\ \hfill v_2(x)& =a_5+a_{6}x+a_7{x}^2+a_8{x}^3;l_1\le x\le l_e\hfill \\ \hfill {\mathit{\theta }}_2(x)& =\frac{\mathrm{d}{v}_2(x)}{\mathrm{d}x}=a_6+2a_{7}x+3a_8{x}^2;l_1\le x\le l_e\hfill \\ \hfill \end{array}$$(33)

where $a_1$ to $a_8$ are the constants, subscript 1 and 2 represents left and right side segment of the element, respectively. The four boundary conditions, at the cracked section where the elastic hinge is located are as follows:(34) $$\begin{array}{cc}\hfill v_1(0)& =Y_1;{\mathit{\theta }}_1(0)={\mathrm{\Theta }}_1\hfill \\ \hfill v_2(l_e)& =Y_2;{\mathit{\theta }}_2(l_e)={\mathrm{\Theta }}_2\hfill \\ \hfill v_1(l_1)& =v_2(l_1);v_2^{″}(l_1)=\frac{1}{K_c}({\mathit{\theta }}_2(l_1)-{\mathit{\theta }}_1(l_1));\hfill \\ \hfill v_1^{″}(l_1)& =v_2^{″}(l_1);v_1^{″\prime }(l_1)=v_2^{″\prime }(l_1)\hfill \\ \hfill \end{array}$$(34)

where $K_c=EIc$ is the flexibility coefficient at the crack location, EI is the flexural stiffness and c be the torsional flexibility of the crack, $Y_1,{\mathrm{\Theta }}_1,Y_2$ and ${\mathrm{\Theta }}_2$ are the displacement and rotation at the nodes 1 and 2, respectively. The flexibility coefficient at the crack location can be also expressed as given in [19](35) $$\begin{array}{c}\hfill K_c=hf(\mathit{\xi })/2\end{array}$$(35)

where $f(\mathit{\xi })$ is a correction function which takes into account the body and crack geometry and depends on the dimensionless crack ratio $\mathit{\xi }=a/h$ which can be represented in as follows [19](36) $$\begin{array}{cc}\hfill f(\mathit{\xi })& ={\mathit{\xi }}^2(12-19.5\mathit{\xi }+70.1{\mathit{\xi }}^2-97.6{\mathit{\xi }}^3+142{\mathit{\xi }}^4-138{\mathit{\xi }}^5\hfill \\ \hfill & +128{\mathit{\xi }}^6-132{\mathit{\xi }}^7+379{\mathit{\xi }}^8-417{\mathit{\xi }}^9\hfill \\ \hfill & +131{\mathit{\xi }}^{10}+313{\mathit{\xi }}^{12}-357{\mathit{\xi }}^{13}+102{\mathit{\xi }}^{14})0\le \mathit{\xi }<0.5\hfill \\ \hfill f(\mathit{\xi })& =\frac{1.32}{(1-{\mathit{\xi }}^2)}-1.780.5\le \mathit{\xi }\le 1\hfill \end{array}$$(36)

Using Equation (34(34) $$\begin{array}{cc}\hfill v_1(0)& =Y_1;{\mathit{\theta }}_1(0)={\mathrm{\Theta }}_1\hfill \\ \hfill v_2(l_e)& =Y_2;{\mathit{\theta }}_2(l_e)={\mathrm{\Theta }}_2\hfill \\ \hfill v_1(l_1)& =v_2(l_1);v_2^{″}(l_1)=\frac{1}{K_c}({\mathit{\theta }}_2(l_1)-{\mathit{\theta }}_1(l_1));\hfill \\ \hfill v_1^{″}(l_1)& =v_2^{″}(l_1);v_1^{″\prime }(l_1)=v_2^{″\prime }(l_1)\hfill \\ \hfill \end{array}$$(34) ) and nodal values, the reduced shape function matrix [N(x)] of size 4 $\times$ 4 relates the nodal DOF with the field variables as follows [18]:(37) $$\begin{array}{c}\hfill \left\{\begin{array}{c}{v}_1(x)\\ {\mathit{\theta }}_1(x)\\ v_2(x)\\ {\mathit{\theta }}_2(x)\end{array}\right\}=\left[\begin{array}{cccc}{N}_{11}(x)& N_{12}(x)& N_{13}(x)& N_{14}(x)\\ N_{21}(x)& N_{22}(x)& N_{23}(x)& N_{24}(x)\\ N_{31}(x)& N_{32}(x)& N_{33}(x)& N_{34}(x)\\ N_{41}(x)& N_{42}(x)& N_{43}(x)& N_{44}(x)\\ \end{array}\right]\left\{\begin{array}{c}{Y}_1\\ {\mathrm{\Theta }}_1\\ Y_2\\ {\mathrm{\Theta }}_2\end{array}\right\}\end{array}$$(37)

The displacement in the cracked beam element can be written in terms of shape function from Equation (7(7) $$\begin{array}{c}\hfill w(x)=N_1{w}_1+N_2{\mathit{\theta }}_1+N_3{w}_2+N_4{\mathit{\theta }}_2=\sum _{j=1}^4{N}_j{w}_j\end{array}$$(7) ),(38) $$\begin{array}{cc}\hfill v_1(x)& =N_{11}(x)Y_1+N_{12}(x){\mathrm{\Theta }}_1+N_{13}(x)Y_2+N_{14}(x){\mathrm{\Theta }}_2\hfill \\ \hfill v_2(x)& =N_{31}(x)Y_1+N_{32}(x){\mathrm{\Theta }}_1+N_{33}(x)Y_2+N_{34}(x){\mathrm{\Theta }}_2\hfill \\ \hfill \end{array}$$(38)

From FE procedure, the stiffness and mass matrices of the cracked element can be obtained as follows(39) $$\begin{array}{cc}\hfill [k_{cr}^e]& =EI\left({\int }_0^{l_1}{\left\{\frac{{\mathrm{d}}^2[N_1(x)]}{\mathrm{d}{x}^2}\right\}}^T\left\{\frac{{\mathrm{d}}^2[N_1(x)]}{\mathrm{d}{x}^2}\right\}\mathrm{d}x+{\int }_{l_1}^{l_e}{\left\{\frac{{\mathrm{d}}^2[N_3(x)]}{\mathrm{d}{x}^2}\right\}}^T\left\{\frac{{\mathrm{d}}^2[N_3(x)]}{\mathrm{d}{x}^2}\right\}\mathrm{d}x\right)\hfill \end{array}$$(39) (40) $$\begin{array}{cc}\hfill [m_{cr}^e]& =\mathit{\rho }A\left({\int }_0^{l_1}{[N_1(x)]}^T[N_1(x)]\mathrm{d}x+{\int }_{l_1}^{l_e}{[N_3(x)]}^T[N_3(x)]\mathrm{d}x\right)\hfill \end{array}$$(40)

where $N_1(x),N_2(x),N_3(x)$ and $N_4(x)$ are the respective rows of shape function matrix N(x). The value of $K_c$ is zero, structure is considered as an intact element. The crack in the beam element has more dependence on change of stiffness, but it does not affect mass distribution appreciably.

10.1. Crack identification (depth and location) using PZT Patch

In this case, damage detection by the proposed method of voltage matching is applied to determine the magnitude (Depth) and location simultaneously for a single crack in the beam structure. The same cracked steel beam was studied in Ref. [13,19] is 520 mm long, 50 mm width and 8 mm thickness with Young’s modulus (E) 206 GPa and density of 7850 kg/m${}^3$.

The fixed–fixed beam is divided into five Euler FEs with an open edge crack and is at a distance of 275.6 mm from the left fixed end as shown in Figure 12. The absolute normalized crack location measured from the left end of the beam is $\mathit{\lambda }(l/L)=0.53$. In the FE model of the fixed–fixed beam, the crack is located in element 3. In this study, four different cases of normalized crack depth ($\mathit{\xi }=a/h$) are considered viz; 0.03, 0.05,0.1 and 0.5, respectively. The normalized crack location is measured from the left node of the element 3 is ${\mathit{\lambda }}_e(l_1/l_e)=0.6$.

In this study, as the first case two PZT patches with size PZT:5% i.e. $(26\times 50\times 1{\mathrm{mm}}^3$ PZT) is bonded on either side of the structure as shown in Figure 12. The first and second natural frequency of the modes of vibration are 156.18 and 436.08 Hz, respectively. Rayleigh damping with the modal damping ratio of 3% is used for the first two modes of vibration. First, the displacement, velocity and acceleration time history data are calculated for each nodal point using Newmark’s method with constant time step of 0.002 s in 0.01 s using the impulse force of 5 N in the vertical (upward) direction at node 4. Using the displacement response history, voltage responses are measured through two PZT patches. The numerically calculated voltages are polluted by artificially adding Gaussian white noise with zero mean and standard deviation of 5% to simulate experimental errors.

Here, Voltage difference between healthy and damaged states of structure for a small crack depth $\mathit{\xi }=0.03$ and $\mathit{\xi }=0.05$ is observed to be very small. The Figure 13 shows the undamaged and damaged responses available at the measurement PZT patch sensor A for normalized crack depths $\mathit{\xi }=0.1$ and $\mathit{\xi }=0.5$ cases, respectively. The shift in damaged response for a large crack is greater than a smaller crack.

Figure 12. FE Model of cracked beam with PZT patch.

It is assumed that the mass and flexural rigidity (EI) are known. Crack magnitude and location are set as the unknown variables for each element and thus there are 10 unknown variables in this problem. Thus the optimization variables to be identified are normalized crack depth $\mathit{\xi }=a/h$ and location ${\mathit{\lambda }}_e=l_1/l_e$ in five elements. The experimentally measured voltage response of PZT patches is required for the fitness evaluation. The mean square error (MSE) between measured and predicted voltage response at PZT patches are minimized by PSO. The lower and upper bound for PSO optimization to identify crack magnitude are set at zero and one. Similarly, the lower and upper bound to identify normalized crack location are set as zero and $l_e$, respectively. In this PSO parameters are to 100 particles( swarm size) and 500 generations as per previous literature.

The identified crack parameters using voltage matching with PZT:5% patches in various normalized crack depths of 0.03, 0.05, 0.1 and 0.5 with constant normalized absolute crack location ($\mathit{\lambda }$) 0.53 as shown in Table 8. The MAE in normalized crack depth is 1.4875 and 4.8%, respectively, for noise free and 5% noise case. Similarly for the normalized absolute crack location, it is 0.1575 and 1.6325%, respectively. It can be seen that identified parameters in noisy environment show more error as compared to noise free case.

Figure 13. Response at measurement point (PZT A) with single crack for (a) $\mathit{\xi }=0.1$ and (b) $\mathit{\xi }=0.5$ case.

Table 8. Identified values of normalized damage magnitude and location using PZT:5%.

	Normalized		Normalized
	Crack septh $\mathit{\xi }$		Crack location $\mathit{\lambda }$
Normalized	(Identified)		(Identified)
Crack depth $\mathit{\xi }$	Noise	With 5%	Noise	With 5%
(Exact)	free (% error)	Noise (% error)	free (% error)	Noise (% error)
0.03	0.031($-3.33$)	0.0327($-9$)	0.532($-0.377$)	0.5421($-2.2$)
0.05	0.0489(2.2)	0.053($-6$)	0.53075($-0.14$)	0.5405($-1.98$)
0.1	0.0996(0.4)	0.104($-4$)	0.5295(0.094)	0.5375($-1.41$)
0.5	0.4999(0.02)	0.499(0.2)	0.5249(0.019)	0.535($-0.94$)
Abs.Normalized	MAE (%)		MAE (%)
Crack location ($\mathit{\lambda }$)	1.4875	4.8	0.1575	1.6325
Exact	Maximum Absolute Error (%)		Maximum absolute error (%)
0.53	3.33	9	0.377	2.2

Figure 14(a) and (b) depict the element wise detection of crack depth and location for the case of $\mathit{\xi }=0.03$. From the plot, it can be seen that all elements have very small magnitudes of normalized crack depth and location except the third element which is the correct cracked element. It is seen from Table 8 that absolute error in noise is 3.33 and 0.377% for crack depth and location. But in the noisy case, it is 9 and 2.2%, respectively.

The Figure 15(a) and (b) shows the element wise details of normalized damage magnitude and location of crack of $\mathit{\xi }=0.1$ with same crack location. In the noise free case errors are very small such as 0.4 and 0.094%, respectively, for crack magnitude and location. For Noisy case, it is in the range of 4 and $-1.41$%, respectively.

Next, the improvement of crack (depth and location) identification of smallest depth i.e. $\mathit{\xi }=0.03$ at the same location using PZT:10% i.e. $(52\times 50\times 1{\mathrm{mm}}^3$ PZT) is identified. The identified results are shown in Table 9. Here the absolute error in noise free case is 2.67 and 0.32% for crack depth and location, respectively. The percentage decrease of error compared to PZT:5% is 19.8% and 15.11%, respectively. But for the noisy case, respective absolute errors are 6 and 2.075% and percentage decrease of error 33.33% and 5.68%, respectively.

Figure 14. Normalized damage magnitude and location of crack with $\mathit{\xi }=0.03$.

Figure 15. Normalized damage magnitude and location of crack with $\mathit{\xi }=0.1$.

Figure 16(a) and (b) shows the convergence of fitness (objective) functions PZT:5 and PZT:10% for $\mathit{\xi }=0.03$ (noise-free and noisy case), respectively. The better minimization of PZT:10% when compared with PZT:5% is seen from the plot. Hence, the error of both magnitude and location is reduced.

10.1.1. Comparison with other damage detection methods

The results of Viola et al. [19] are compared with results obtained from the proposed method. In Ref. [19], a single crack with depth $\mathit{\xi }=0.5$ and absolute crack location $\mathit{\lambda }=0.76$ was identified using a frequency domain method from experimentally obtained measurements of unspecified noise content. The identified structural parameters are 0.486 and 0.756 and the percentage of the absolute error is 2.8 and 0.53%, respectively, for crack depth and location. The very same problem was also solved by Varghese and Shankar [13] using the combined method of power flow and acceleration matching using substructure. The percentage of error in identifying crack depth and location was 0.92 and 0.5%, respectively, without noise and 1.72 and 2.8%, respectively, with 5% Noise. There, mean computational time was about 3000 s whereas with the proposed method it is 140 s Sumant and Maiti [21] predicted the magnitude and location of $\mathit{\xi }=0.03$ crack with errors of 9.33 and 4.95%, respectively, using static response of pair of PZT patch. But, in that case distance between sensor patch and crack should be less than crack depth (sensor is very near to crack) and also require undamaged structure details. The proposed method using voltage matching by PZT:5% identified the normalized crack depth ($\mathit{\xi }=0.5$) and location ($\mathit{\lambda }=0.53$) was 0.02 and 0.019%, respectively, for noise free case and 0.2 and 0.94%, respectively, for noisy condition. And $\mathit{\xi }=0.03$, it was 9 and 2.2%, respectively, for depth and location in noisy condition. The results show that PZT-based damage identification using voltage is a promising method for crack depth and location identification.

Table 9. Identified values of normalized damage magnitude and location using PZT:10%.

	Normalized		Normalized
	Crack depth $\mathit{\xi }$		Crack location $\mathit{\lambda }$
Normalized	(Identified)		(Identified)
Crack depth $\mathit{\xi }$	Noise	With 5%	Noise	With 5%
(Exact)	free (% error)	Noise (% error)	free (% error)	Noise (% error)
0.03	0.0308($-2.67$)	0.0318($-6$)	0.5317($-0.320$)	0.541($-2.075$)

Figure 16. Convergence plot for $\mathit{\xi }=0.03$ (a) Noise-free and (b) 5% Noise case.

The fixed–fixed beam described in the previous case is studied with two cracks and assumed that the cracks in different elements do not interact with each other. The cracks are assumed to be located at 143 and 221 mm from left end of beam, respectively. The normalized crack depth are 0.05 and 0 .5 and the corresponding normalized locations($\mathit{\lambda }$) from left end of fixed beam are 0.275, 0.425, respectively, for two cracks. Here, the crack lies in the second, third element and normalized crack location $({\mathit{\lambda }}_e)$ from left end of the corresponding element are 0.75 and 0.25. The crack depth and location are identified accurately as listed in Table 10. The accuracy of identification is as good as single crack identification.

11. Experimental work:fixed–fixed beam with PZT Structure

In this preliminary work, a brief experiment is carried out to identify stiffness and damping parameters of a fixed–fixed beam. A fixed–fixed beam made of acrylic material with dimensions $(450\times 25\times 12)$ mm is used for the experimental study. Here two PZT patch of dimensions $(25\times 25\times 1)$ mm each are bonded at the fixed ends of the structure as shown in Figure 17. The modulus of Elasticity (E) was estimated to be 3.9 GPa from a simple bending test and the density was measured to be 1190 kg/m${}^3$. The actual flexural rigidity (EI) of the beam is 14.04 N.m${}^2$. The damping ratio($\mathit{\zeta }$) was calculated from a simple free vibration decay test using logarithmic decrement and estimated as 8%. The natural frequencies for the first two modes of the structure are calculated from the frequency domain as 110 and 306 Hz. Assuming Rayleigh’s proportional damping model, exact values of damping constants $\mathit{\alpha }$ and $\mathit{\beta }$ are calculated as 81.64 and 6.12$\times {10}^{-5}$, respectively, for two modes. The beam was divided into five elements as shown in Figure 6.

Table 10. Double crack identification (magnitude and location) using PZT:5%.

		Normalized		Normalized
		Crack depth $\mathit{\xi }$		Crack location ${\mathit{\lambda }}_e$
Normalized	Normalized	(Identified)		(Identified)
Crack depth $\mathit{\xi }$	Crack location $\mathit{\lambda }$	Noise-free	5% Noise	Noise-free	5% Noise
(Exact)	(Exact)	(%Error)	(%Error)	(%Error)	(%Error)
0.05	0.275	0.0521($-4.2$)	0.0536($-7.2$)	0.2715(1.27)	0.2803($-1.93$)
0.5	0.425	0.4956(0.88)	0.5057($-1.14$)	0.4266($-0.38$)	0.4288($-0.89$)

Table 11. Experimentally identified parameters of fixed–fixed beam.

	Exact	Identified
Parameter	$N{m}^2$	$N{m}^2$	% of Error
$E{I}_1$	14.1656	15.24	$-7.58$
$E{I}_2$	14.04	13.35	4.91
$E{I}_3$	14.04	13.56	3.42
$E{I}_4$	14.04	14.96	$-6.55$
$E{I}_5$	14.1656	13.46	4.98
Mean Absolute Error (MAE)(%)			5.49
$\mathit{\alpha }$	81.64	90.5	$-10.8$
$\mathit{\beta }$	$6.12\times {10}^{-5}$	$6.65\times {10}^{-5}$	$-8.6$

Figure 17. Experimental set-up of fixed–fixed beam with PZT Patch.

Figure 18. Comparison of experimental and simulated voltage response at a PZT patch.

The beam is excited by a sinusoidal force of $2.8sin(2\mathit{\pi }\times 60.7t)$ N at the middle of the structure by a LDS permanent magnet 20 N modal shaker with a maximum displacement of 5 mm with an operating frequency range of 5 Hz$-13$ kHz. The applied force is measured by using a KISTLER force transducer and is acquired with a sampling frequency of 1000 Hz using DEWE 43 DAQ system. The dynamic voltage response is measured through two piezo patches and it is sampled at the rate of 1000 Hz. From the acquired data, a portion of data length 5 s is considered for parameter identification. The unknown structural parameters are identified by minimizing the MSE between the measured and predicted response using Equation (31(31) $$\begin{array}{c}\hfill f=\frac{{\sum }_{i=1}^M{\sum }_{j=1}^L{\left|{\mathit{\varphi }}^m(i,j)-{\mathit{\varphi }}^e(i,j)\right|}^2}{ML}\end{array}$$(31) ). The PSO parameters are set to 100 particles (swarm size) and 500 generations. The search range for structural parameters is set as $\pm 50\%$ of actual value. The identified parameters are shown in Table 11. The structural stiffness parameters are identified with mean absolute error 5.49% and damping parameters are identified with maximum absolute error of 10.8%. The mean CPU time required for convergence is 125 s.

Figure 18 shows the matching between the experimental voltage of the left patch and numerically predicted voltage at successful identification.

12. Conclusions

In this article, a novel structural parameter identification based on voltage matching by using PZT patch sensor is presented. Unlike conventional accelerometer sensor, PZT patch is low cost, negligible weight,wide band width and dynamic measurement of distributed nature. Mathematical model of the PZT sensor with host structure is developed on the basis of one-dimensional geometry and reduced material property. The structural identification is done by the minimization of MSE of the deviation between the measured and predicted voltage responses with unknown structural parameters as optimization variables using PSO. The effect of different PZT length is investigated with numerical example and convergence of the fitness functions of those case are studied. The result shows that identification using PZT:5% (smallest patch length under study) is even sufficient for accurate structural parameter identification and it is successfully verified experimentally. The Mean Absolute Error (PZT:5%) are very low in the range of 0.167–0.959% (noisy condition) for structural parameter identification as compared to other existing methods. It is also used for identifying the crack depth and location accurately with an MAE of 4.8 and 1.63%, respectively, for noisy signal. The PZT method was capable of detecting crack depth more accurately than previous time and frequency domain methods.

Notes

No potential conflict of interest was reported by the authors.