Vibration Analysis of Plate with Irregular Cracks by Differential Quadrature Finite Element Method

A universal method combining the differential quadrature finite element method with the virtual spring technique for analyzing the free vibration of thin plate with irregular cracks is proposed. Translational and rotational springs are introduced to restrain the vertical displacement and orientation of the plate.Themass matrix and stiffness matrix for each element are deduced involving the effects of the virtual springs. The connection relationships between elements can be modified by setting the stiffness of the virtual springs.The vibration of two rectangular plates with three irregular shaped cracks and different boundary conditions are presented. The results are compared with those obtained by ANSYS, where the good agreement between the results validates the accuracy and efficiency of the present method.


Introduction
Cracks occurring in the plate cause the changes of the natural frequencies and the dynamic response behavior under loads.The effects of the cracks on the vibration characteristic of the structures have been widely studied.Krawczuk [1] proposed a FEM model of a rectangular plate with a through crack by modifying the stiffness matrix of the plate element to simulate the effect of the crack.Yong-gang et al. [2] studied the nonlinear vibrations for moderate thickness rectangular plate with surface penetrating crack of different location and depth.Qu et al. [3] proposed an approach to analyze the vibration of a piezoelectric composite plate with multiple cracks, which may be used for damage detection.An extended finite element method is proposed by Bachene et al. [4] to study the vibration of cracked plates.Vibration of rectangular plate and square plate with through-edge and center cracks are studied.This method is later used to analyze the vibration of cracked functionally graded material plates by Natarajan et al. [5], where only straight cracks are studied.Huang and Leissa [6] used the Ritz method to analyze the vibration of rectangular plate with side cracks.A new set of functions added to the traditional admissible functions is proposed, which indeed accelerate the convergence of the numerical solutions.The effects of the location and orientation of the side crack are studied.Later the method is applied to the thick rectangular plate [7,8], where the internal cracks are taken into account.Hosseini-Hashemi et al. [9] firstly proposed a set of exact closed-form characteristic equations incorporating the shear deformation and rotary inertia to analyze the free vibration of moderately thick rectangular plates with multiple all-over part-through cracks.However, the boundary conditions of the two opposite edges are limited to be simply supported and the cracks are required to be perpendicular to the simply supported edges.Natarajan et al. [10] studied the vibration of functionally graded material plates with a through center crack using an 8-noded shear flexible element.However, only the effect of the lengths of cracks is studied without presenting further discussions on the cracks.Ismail and Cartmell [11] investigated the forced vibration of a plate with surface crack of variable angular orientation.Bose and Mohanty [12] analyzed the vibration of a rectangular thin isotropic plate with a part-through surface crack with different orientation and location.Torabi et al. [13] 2 Shock and Vibration applied the differential quadrature element method (DQEM) to analyze the free transverse vibration of a nonuniform Timoshenko beam with multiple cracks.The compatibility conditions at the damaged section are considered as a discontinuity in slope and vertical displacements.Heydari et al. [14] proposed a continuous model for flexural vibration of Timoshenko beams with vertical edge crack, where the effects of shear deformation and rotary inertia are considered.Broda et al. [15] studied the longitudinal vibration of beams with breathing cracks.Joshi et al. [16,17] proposed an analytical modeling and studied the vibration of a rectangular plate with internal crack.The cracks are limited to be perpendicular to each other and parallel to one of the edges and located at the center of the plate.Y. C. Chen developed a special boundary element method (BEM) to analyze the elastodynamic of anisotropic elastic plates with holes and cracks.Nguyen-Thoi et al. [18] used an extended cell-based smoothed discrete shear gap method to analyze the free vibration of cracked Mindlin plate, where triangular elements are used.The effects of location and orientation of the cracks on the vibration characteristic are studied.Azam et al. [19] proposed a FEM model to analyze the vibration of a square plate with side crack.Bhardwaj et al. [20] used the extended isogeometric analysis to simulate the cracked functionally graded material palates.
According to the literature available, most of the researches are focused on the straight cracks with various locations and orientations.Few literatures on more complex cracks can be found.Monfared and Bagheri [21] studied a functionally graded material palate with multiple interacting arbitrary shaped cracks, where circular cracks and parabolic cracks are considered.Raffo and Quintana [22] presented a general algorithm to obtain approximate analytical solutions for the vibrations of a rectangular plate with internal curved line hinge, which can be regarded as a kind of crack where only the translational displacement is equal on the two sides of the crack.
In order to propose a universal method for vibration analysis of plate with irregularly shaped through cracks, a differential quadrature finite element method (DQFEM) proposed by Xing et al. [23,24] is used.The DQFEM has been widely used to model plate structures with complex geometries [25][26][27][28][29][30], which has been proved to be an efficient method using less element numbers and obtaining high accurate results.This method makes it possible to model cracks with irregular shapes simply by disconnecting the adjacent elements so that a crack emerges between the elements.However, each time the locations of cracks change, the stiffness matrix is need to be modified to satisfy the new connection relationships.To overcome the disadvantages of the original DQFEM, the virtual spring technique is introduced.The virtual spring technique has been used widely in simulation of the elastic boundary conditions as well as the coupling conditions [31,32].Thus, the connection relationships between elements can be easily modified just by setting the stiffness of the virtual springs.
In this paper, the virtual spring technique is introduced into the DQFEM.Based on the DQFEM theory, the mass matrices and stiffness matrices involving the springs' stiffness of the elements are deduced.The vibration results of two rectangular plates with different shaped cracks are presented.Based on the DQFEM theory, the mass matrices and stiffness matrices involving the spring's stiffness of the elements are deduced.The convergence studies on the stiffness of the springs are carried out.The vibration results of two rectangular plates with different shaped cracks are presented.The results are compared with those obtained by ANSYS.

Theory and Formulations
2.1.The Differential Quadrature Rules.The differential quadrature (DQ) method approximates the th derivative of a onedimensional function (), which is derivable in the interval [, ], at a point  between the derivable interval by using a summation including the weighting coefficients and the function values at all grid points in the solution domain.According to the differential quadrature rules, the th-order derivative of function () can be written as where  is the grid numbers in the interval. ()  is the th-order weighting coefficient; thus, [ ()   ] is the th-order weighting coefficient matrix.The off-diagonal terms of the first-order weighting coefficient matrix are given by (2).The higher-order weighting coefficient matrix can be derived by the recurrence relationships given by  (1)   = According to the two-dimensional DQ rule, the partial derivatives of function (, ) at grid point (  ,   ) can be written as Shock and Vibration 3 where  ()  and  ()  are the weighting coefficients along  and  directions, respectively. and  are the numbers of the grid point in  and  directions in the solution domain.For simplicity, a single index notation is used to present the compact form of the two-dimensional DQ rule, which are given as follows [24]: where ,  = ( − 1) + , ( = 1, 2, . . ., ;  = 1, 2, . . ., ) and where A () and B () are both ( × ) × ( × ) matrix, representing the weighting coefficient matrices related to the derivative in  and  direction, respectively, in which Note that the two-dimensional DQ rules can be treated as one-dimensional DQ rule applied both in  and  direction.Thus,  ()  =  ()   , where  ()  represents the th-order weighting coefficient along  direction of grid point at th row and th column. ()  represents the th-order weighting coefficient along  direction of grid point at th row and th column.

Equations of Motion of the Rectangular Plate Element.
The concept of finite element in finite element method is adopted into the DQFEM [25].In order to deduce the equations of motion of the whole structure, we need to start from a single element.Thus, the equations of motion of a rectangular element or a square element are firstly given here.A rectangular plate element is shown in Figure 1, in which two sets of virtual springs are attached to the four sides of the element .The first set of virtual springs connected the plate element to the ground, which are combined by translational and rotational springs.The translational springs restrict the translations of the grid points along the normal direction of the element, while the rotational springs restrict the rotations of the points along the normal direction of each side of the element.The stiffness of the translational and rotational springs is denoted by    ,    , and    , where  means the eth element of the structure,  and  represent the translational and rotational springs, and  and  denote the normal directions of each side of an element. is the th side of a square element.For    , due to the fact that the only normal directions of side 1 and side 3 are in  direction, thus  = 1, 2 for the rotational springs, which is the same as    and can be observed in Figure 1.The other set of virtual springs connected the element to an adjacent element, which also consists of translational and rotational springs.The stiffness of these coupling springs is expressed by  ,  ,  ,  , and  ,  , where the superscript ,  represents the coupling springs connecting the eth element and eith element,  means the adjacent element is connected to the th side of element e, and c indicates the coupling springs.For simplicity, the first set of virtual springs is called element springs, while the other sets are called coupling springs.
By introducing the virtual springs into the DQFEM, the elastic boundary conditions and the elastic coupling relationships between elements can be easily achieved.In other words, by simply setting the stiffness of the coupling springs to be zero or an infinite large value, the connected or broken relationships between two adjacent elements can be easily achieved.For broken relationships, we can regard them as cracks in the structures.
Then we consider a rectangular Kirchhoff thin plate element divided into  grid points along the  direction and  grid points along the  direction with two sets of virtual springs attached to the element.The displacement function of the th element is defined as where   are the Lagrange polynomials and    =   (  ,   ) is the displacement at the Gauss Lobatto quadrature points [24].Then the potential energy and kinetic energy of eth element are given as where  is the flexural rigidity,  is Poisson's ratio,  is density, and ℎ is the thickness of the plate; for thin plate ℎ/ = 0.001, where  is the side length of the plate.The total potential energy of the translational element springs and the two types of rotational springs is given as The total potential energy in the coupling springs is given as where the displacements    and    are the adjacent grid points of two adjacent elements as shown in Figure 2(a).What should be noted is that the potential energy in the coupling springs belongs to the two adjacent elements; thus, for one element the potential energy should be half of that in (11).According to Hamilton's principle, one has where Ne is the total numbers of the element in the structure.
According to the DQ rules, the expressions can be transformed into matrix form.The integrals in the energy expressions are calculated by using the Gauss Lobatto quadrature [24].The potential and kinetic energy of the eth element is where F (2) = A (1) B (1) and the displacement vector is which can be illustrated clearly through Figure 2(b).Matrix C = diag( 1 ,  2 , . . .,   ) is a  ×  diagonal matrix where the diagonal elements are the weighting coefficients of Gauss Lobatto quadrature for one-dimensional problem.Thus, the matrix for two-dimensional problem is Element e2 The matrix form of the potential energy of element springs is where the stiffness matrices of the translational springs that belong to each side of the element are given as The stiffness matrices of the rotational springs are given as The above stiffness matrices are sparse matrices and the distributions of the matrix elements in the sparse matrix can be understood through the above expressions.The matrix form of the potential energy of the coupling springs is given as w where the sign "=" represents that the matrix or the vector consists of matrices or vectors of two adjacent elements, which can be expressed as The matrices of the stiffness of the coupling springs are given as  (17).In order to illustrate the matrix −K ,  , we can assume that the 3rd side of element  is connected to the 1st side of element ei; then the matrix −K

𝑒,𝑒3
3 can be expressed as One can observe that, in the coupling matrix K

𝑒,𝑒3
3 are in the same rows as those in matrix K  3 and in the same columns as those in matrix K 3 3 .
Thus, the matrix −K 3, 3 can be easily understood which is given as

Equations of Motion of the Whole
Structure.So far the energy expressions for a single element have been given and illustrated, the next step is to assemble the expressions of each element together, which is similar to the operation in FEM.According to (12), the potential and kinetic energy of the whole structure is where 2) , A (2) , . . ., A (2) ) ,  B (2) = diag (  ⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ B (2) , B (2) , . . ., B (2) ) , Shock and Vibration where  is the total number of elements.The potential energy of the element springs in the whole structure is given as where Due to the fact that the there are no coupling effects between element springs of two different elements, the global matrix of the element springs is diagonal matrix.However, the coupling effects should be considered in the assembling process of the coupling matrices.The global matrix of the coupling element is obviously a nondiagonal matrix, in which the distributions of the nondiagonal matrix are determined by the serial numbers defined in the meshing process.The potential energy of all the coupling springs is Consequently, the potential energy of the whole system including the structure and the virtual springs can be obtained by the summation of ( 24), (27), and ( 29) In order to satisfy the continuity requirements between multiple plate elements and simplify the boundary conditions imposing procedures, the displacement vector is modified as [24] w where   = / represents the rotation in  direction,   = / represents the rotation in  direction, and   =  2 / represents the torque at the corner grid point.Using the DQ rule, one can obtain the following relation: For more information about matrix Q, one can refer to [24].The global matrix of the transformation matrix Q is given as Substituting ( 33) into ( 30) and ( 25) and taking out the displacement vector, then the stiffness matrix and mass matrix of the whole system can be written as +  B (1) ( Finally, the vibration problem became eigenvalue problem which can be solved by solving the following equation: 2.4.Curvilinear Quadrilateral Plate Element.Since the regular DQ method cannot be applied directly to the irregular domain, the mapping technique should be used to deal with problems with irregular geometries.The mapping technique has been adopted in many literatures [24][25][26]28], so the detailed description of the mapping technique will not be given here.The mapping process in [24] is adopted in this paper, where the serendipity-family interpolation functions are used in the mapping of an irregular domain, which usually contains four sides, into a regular square domain defined in −1 ≤  ≤ 1, −1 ≤  ≤ 1.Since all the elements considered in this research are straight-edge elements, the 3rd-order serendipity functions are used, which means on each side of the irregular domain there are four uniformly distributed points.According to the literature results, the 3rd-order serendipity functions are enough in the mapping process to obtain accurate results.

Numerical Examples
In this section, several numerical examples will be presented to validate the accuracy of the method combined with DQFEM and the virtual springs technique.First of all, a convergence study of the numbers of grid points of a square plate element with free boundary conditions is presented in Figure 3, from which one can observe that the first six frequency parameters  =  2 (ℎ/) 1/2 show a good convergence property whenever the numbers are even or odd.Based on the convergence study, the numbers of the grid points are set as  =  = 12 for the following researches, which are sufficient to achieve accurate results.The material parameters of the structures are given as Young's modulus  = 2.111 Pa, density  = 7800 kg/m 3 , Poisson's ratio  = 0.3, and thickness ℎ = 0.001 m.

Convergence Analysis of the Virtual Springs.
Firstly, the convergence of stiffness of the element springs is studied, where a square plate element with virtual springs attached to each side while the coupling springs are not concerned is used.By changing the stiffness of the translational springs on four sides uniformly from 0 to an infinitely large number, firstly, the changing rule of the frequency parameters  =  2 (ℎ/) 1/2 of the free vibration of the square plate element can be observed in Figure 4.The frequency parameters remain constant as the stiffness of the spring changes into a relatively large number.This is when the free boundary conditions turn into the simply supported boundary conditions.Here we take   / = 111 as the infinitely large number for the translational springs.Next keep the   / = 111 and increasing the stiffness of the rotational springs on four sides uniformly from 0 to an infinitely large number, similar change rule of the frequency parameters can also be observed in Figure 4.And this is when the simply supported boundary conditions turn into the clamped boundary conditions.We take   / = 111 as the infinitely large number for the rotational springs.Therefore, the classical boundary conditions can be easily achieved by setting the stiffness of the virtual springs to be 0 or an infinitely large number and the elastic boundary conditions are achieved by setting the stiffness to be a required value.
Next the convergence studies of the stiffness of the coupling springs are carried out by the same way as those above.The two elements are clamped at the sides parallel to the coupling side, while the other sides are free.Increasing the stiffness of the translational springs on the coupling side uniformly while keeping the stiffness of rotational springs to be zero, the frequency parameters change and then remain constant when the stiffness changes into an infinitely large number according to Figure 5.The coupling condition between the two element plates now is like the hinge connection, where the translational displacement of one element is restricted by the other element, while the rotational displacements are not restrained.We can take   / = 111 as the infinitely large number for the translational springs.Keeping   / = 111 and changing the stiffness of the rotational springs on the coupling side uniformly, one can observe similar change rule of the frequency parameters.We can take   / = 111 as the infinitely  large number for the rotational springs.Thus, the elastically coupling conditions can be simulated as those for boundary conditions.

Results of Plate with Irregular Cracks.
In this section, two examples of a rectangular plate with several kinds of cracks are to be given.The geometry of the first example is shown in Figure 6, where length  = 2m and width  = 1m.The rectangular plate is uniformly divided into 8 elements, which are numbered as shown in Figure 6.Three types of cracks are presented here: T-shaped, L-shaped, and I-shaped cracks, which are shown in Figure 7, and the boundary conditions corresponding to each crack plate are CCCC, FCFC, and FCFF, respectively, where C and F denote the clamped and free boundary conditions, respectively.Take the 2nd element of the T-shaped crack plate as an example to illustrate the process of setting the stiffness of the virtual springs.According to the boundary conditions of the plate and the position of the crack, one can find that the 1st side is clamped, 2nd and 3rd sides are free because they are not coupled with other elements and on boundary conditions are applied, and 4th side is coupled with element 1.Thus, the stiffness of the element and coupling springs can be set as Similar operations can be applied for plate with different crack types and boundary conditions.Numerical results are presented in Table 1, where the results are compared with those obtained by ANSYS with Shell 181 element and 20000 elements.The good agreement between the results can be observed, which validates the accuracy and efficiency of the present method.The mode shapes are also presented in Figures 8 and 9.
In the second example, the rectangular plate with the same geometric parameters is divided uniformly into 8 right trapezoid elements shown in Figure 10.Three types of cracks including T-shaped, L-shaped, and I-shaped cracks are also considered, shown in Figure 11, where the boundary conditions corresponding to each crack type are the same as those in the first example.The results are presented in Table 2 and the mode shapes are shown in Figures 12 and 13.The results match well with those obtained by ANSYS.
Through the two examples above, the virtual springs introduced into the DQFEM simplified the procedures of implementation of the boundary conditions as well as the coupling conditions in the origin DQFEM.The virtual spring  technique makes it convenient to set the boundary and coupling conditions between elements.The research in this paper may serve as basic knowledge for further studies on complex structures with complicated boundary and coupling conditions in engineering problems.

Conclusions
The virtual spring technique is introduced into the DQFEM to simulate the elastic boundary and coupling conditions.An element plate is restrained by two sets of virtual springs, of which the first set is called element springs which only  The research in this paper may serve as basic knowledge for further studies on complex structures with complicated boundary and coupling conditions in engineering problems.

Figure 1 :
Figure 1: The coupling relationships between element e and four adjacent elements.

Figure 2 :
Figure 2: (a) The relative positions of the coupling grid points on the adjacent sides of two elements and (b) the distributions of the grid points.
the eth and eith element, respectively.The nondiagonal elements −K ,  and −K ,  represent the coupling effects between two elements.The elements distributions in the sparse matrix for diagonal elements K   and K   are the same as that in So the nonzero elements in the coupling matrices −K ,3 3 and −K 3, 3 are symmetrically distributed about the diagonal line of matrix K ,3 3 .For matrices K ,  and K ,  , the nonzero elements are distributed in the same way.

Figure 3 :
Figure 3: The first six frequency parameters of a square plate with different numbers of grid points.

Figure 4 :Figure 5 :
Figure 4: The changing rules of the frequency parameters as the stiffness of the element springs increasing.

Figure 6 :
Figure 6: A rectangular thin plate divided uniformly into 8 square elements and the serial numbers of each element and side.

Figure 7 :
Figure 7: Geometries of three types of cracks in a rectangular plate and the boundary conditions.

Figure 8 :
Figure 8: The first nine mode shapes of the CCCC plate with T-shaped crack.

Figure 9 :
Figure 9: The first nine mode shapes of the FCFC plate with L-shaped crack.

Figure 10 :Figure 11 :
Figure 10: A rectangular thin plate divided uniformly into 8 right trapezoid elements and the serial numbers of each element and side.

Figure 12 :Figure 13 :
Figure 12: The first nine mode shapes of the CCCC plate with irregular T-shaped crack.

Table 1 :
Natural frequencies of rectangular plates with several classical boundary conditions and T-shaped, L-shaped, and I-shaped cracks.
a Results from ANSYS.

Table 2 :
Natural frequencies of rectangular plates with several classical boundary conditions and irregular T-shaped, irregular L-shaped, and irregular I-shaped cracks.
a Results from ANSYS.