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

Authors: Zhongyuan Xu; Wei Chen
Journal: Shock and Vibration
 (2017)

---

## Abstract

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. The mass 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.


---

## Body

## 1. 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] 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–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.


## 2. Theory and Formulations

### 2.1. The Differential Quadrature Rules

The differential quadrature (DQ) method approximates thesth derivative of a one-dimensional function f(x), which is derivable in the interval [a,b], at a point x 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 sth-order derivative of function f(x) can be written as(1)dsfxdxsi=∑j=1NAijsfj,where N is the grid numbers in the interval. Aij(s) is the sth-order weighting coefficient; thus, [Aij(s)] is the sth-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(2)Aij1=∏k=1,k≠iNxi-xkxi-xj∏k=1,k≠jNxj-xk(3)Aijs=sAiis-1Aij1-Aijs-1xi-xjFori≠j,2≤s≤N-1Aiis=-∑j=1,j≠iNAijs1≤s≤N-1.According to the two-dimensional DQ rule, the partial derivatives of functionf(x,y) at grid point (xi,yj) can be written as(4)∂rf∂xrij=∑m=1MAimrfmj,∂sf∂ysij=∑n=1NBjnsfin,∂r+sf∂xr∂ysij=∑m=1MAimr∑n=1NBjnsfmn,where Aim(r) and Bjn(s) are the weighting coefficients along x and y directions, respectively. M and N are the numbers of the grid point in x and y 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]:(5)∂rf∂xrk=∑m=1M×NA-kmrf-m,∂sf∂ysk=∑m=1M×NB-kmsf-m,∂r+sf∂xr∂ysk=∑m=1M×NF-kmrf-m,where k, m=(j-1)M+i, (i=1,2,…,M;j=1,2,…,N) and(6)f-=f11⋯fM1f12⋯fM2⋯f1N⋯fMNTA-r=Ar0⋯00Ar⋯0⋮⋮⋱⋮00⋯Ar,B-s=B11sB12s⋯B1NsB21sB22s⋯B2Ns⋮⋮⋱⋮BN1sBN2s⋯BNNs,where A-(r) and B-(s) are both (M×N)×(M×N) matrix, representing the weighting coefficient matrices related to the derivative in x and y direction, respectively, in which(7)Ar=A11rA12r⋯A1MrA21rA22r⋯A2Mr⋮⋮⋱⋮AM1rAM2r⋯AMMr,Bijs=Bijs0⋯00Bijs⋯0⋮⋮⋱⋮00⋯Bijs.Note that the two-dimensional DQ rules can be treated as one-dimensional DQ rule applied both inx and y direction. Thus, Aijs=Bijs, where Aijs represents the sth-order weighting coefficient along x direction of grid point at jth row and ith column. Bijs represents the sth-order weighting coefficient along y direction of grid point at jth row and ith column.


### 2.2. 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 e. 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 ktie, kryie, and krxie, where e means theeth element of the structure, t and r represent the translational and rotational springs, and x and y denote the normal directions of each side of an element. i is the ith side of a square element. For kryie, due to the fact that the only normal directions of side 1 and side 3 are in y direction, thus i=1,2 for the rotational springs, which is the same as krxie 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 kctie,ei, kcryie,ei, and kcrxie,ei, where the superscript e,ei represents the coupling springs connecting theeth element andeith element, ei means the adjacent element is connected to the ith side of elemente, andc indicates the coupling springs. For simplicity, the first set of virtual springs is called element springs, while the other sets are called coupling springs.Figure 1
The coupling relationships between element e and four adjacent elements.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 intoM grid points along the x direction and N grid points along the y direction with two sets of virtual springs attached to the element. The displacement function of the eth element is defined as (8)w-ex,y=∑i=1M∑j=1Mlixljxwije,where li are the Lagrange polynomials and wije=we(xi,yj) is the displacement at the Gauss Lobatto quadrature points [24]. Then the potential energy and kinetic energy ofeth element are given as(9)Upe=D2∫∫Se​∂2w-e∂x22+∂2w-e∂y22+2υ∂2w-e∂x2∂2w-e∂y2+21-υ∂2w-e∂x∂y2dxdyTe=12∫∫Se​ρhw-˙e2dxdy,where D is the flexural rigidity, υ is Poisson’s ratio, ρ is density, and h is the thickness of the plate; for thin plate h/a=0.001, where a 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(10)Vse=12∑i=14∫0liektiew-ie2dlie+12∑i=12∫0liekrxie∂w-ie∂x2dxie+12∑i=12∫0liekryie∂w-ie∂y2dyie.The total potential energy in the coupling springs is given as(11)Vce,ei=12∑i=14∫0liekctie,eiw-iei-w-ie2dlie+12∑i=12∫0liekcrxie,ei∂w-iei∂x-∂w-ie∂x2dxie+12∑i=12∫0liekcryie,ei∂w-iei∂y-∂w-ie∂y2dyie,where the displacements w-iei and w-ie 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(12)∑e=1NeUpe+Vse+12Vce,ei=∑e=1NeTe,whereNe is the total numbers of the element in the structure. According to the DQ rules, the expressions can be transformed into matrix form.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.
(a)
(b)The integrals in the energy expressions are calculated by using the Gauss Lobatto quadrature [24]. The potential and kinetic energy of theeth element is (13)Upe=D2w-eTA-2TC-A-2+B-2TC-B-2+2υA-2TC-B-2+21-υF-2TC-F-2w-eTe=12ρhω2w-eTC-w-e,where F-(2)=A-(1)B-(1) and the displacement vector is (14)w-e=w11e⋯wM1ew12e⋯wM2e⋯w1Ne⋯wMNeT,which can be illustrated clearly through Figure 2(b). Matrix C=diag(C1,C2,…,CM) is a M×M 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 (15)C-=diagC1C1,…,CMC1,C1C2,…,CMC2,…,CMCN.The matrix form of the potential energy of element springs is(16)Vse=12w-eT∑i=14C-K-tie+∑i=12A-1TC-K-rxieA-1+∑i=12B-1TC-K-ryieB-1w-e,where the stiffness matrices of the translational springs that belong to each side of the element are given as(17)K-t1e=diagkt1e,kt1e,…,kt1e︷M,0,…,0︷N-1MK-t3e=diag0,…,0︷N-1M,kt3e,kt3e,…,kt3e︷MK-t2e=diagKt2e,Kt2e,…,Kt2e︷N,Kt2e=diag0,…,0︷M-1,kt2eK-t4e=diagKt4e,Kt4e,…,Kt4e︷N,Kt4e=diagkt4e,0,…,0︷M-1.The stiffness matrices of the rotational springs are given as(18)K-rx1e=diagKrx1e,Krx1e,…,Krx1e︷N,Krx1e=diagkrx1e,0,…,0︷M-1K-rx2e=diagKrx2e,Krx2e,…,Krx2e︷N,Krx2e=diag0,…,0︷M-1,krx2eK-ry1e=diagkry1e,kry1e,…,kry1e︷M,0,…,0︷N-1MK-ry2e=diag0,…,0︷N-1M,kry2e,kry2e,…,kry2e︷M.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(19)Vce,ei=12∑i=14w--e,eiTC--K--ctie,eiw--e,ei+∑i=12w--e,eiTA--1TC--K--crxie,eiA--1w--e,ei+∑i=12w--e,eiTB--1TC--K--cryie,eiB--1w--e,ei,where the sign “=” represents that the matrix or the vector consists of matrices or vectors of two adjacent elements, which can be expressed as(20)w--e,ei=w-ew-eiT,C--=diagC-,C-,A--=diagA-,A-,B--=diagB-,B-.The matrices of the stiffness of the coupling springs are given as(21)K--ctie,ei=K-ctie-K-ctie,ei-K-ctiei,eK-ctiei,K--crxie,ei=K-crxie-K-crxie,ei-K-crxiei,eK-crxiei,K--cryie,ei=K-cryie-K-cryie,ei-K-cryiei,eK-cryiei.For matrixK--ctie,ei, the diagonal elements K-ctie and K-ctiei are only related to theeth andeith element, respectively. The nondiagonal elements -K-ctie,ei and -K-ctiei,e represent the coupling effects between two elements. The elements distributions in the sparse matrix for diagonal elements K-ctie and K-ctiei are the same as that in (17). In order to illustrate the matrix -K-ctie,ei, we can assume that the 3rd side of element e is connected to the 1st side of elementei; then the matrix -K-ct3e,e3 can be expressed as(22)-K-ct3e,e3=-00⋯0⋮⋮⋱⋮00⋯0kct1e,e10⋯00kct1e,e1⋯0⋮⋮⋱⋮00⋯kct1e,e1︸M0⋯0⋮⋱⋮0⋯00⋯00⋯0⋮⋱⋮0⋯0︸N-1M.One can observe that, in the coupling matrixK--ct3e,e3, the nonzero elements in -K-ct3e,e3 are in the same rows as those in matrix K-ct3e and in the same columns as those in matrix K-ct3e3. Thus, the matrix -K-ct3e3,e can be easily understood which is given as (23)-K-ct3e3,e=-0⋯00⋯0⋮⋱⋮0⋯00⋯0⋮⋱⋮0⋯0︷N-1Mkct1e,e10⋯00kct1e,e1⋯0⋮⋮⋱⋮00⋯kct1e,e100⋯0⋮⋮⋱⋮00⋯0︷M.So the nonzero elements in the coupling matrices-K-ct3e,e3 and -K-ct3e3,e are symmetrically distributed about the diagonal line of matrix K--ct3e,e3. For matrices K--crxie,ei and K--cryie,ei, the nonzero elements are distributed in the same way.


### 2.3. 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 (24)Up=D2w⌢TA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2w⌢(25)T=12ρhω2w⌢TC⌢w⌢,where (26)w⌢=w-1w-2⋯w-NeTA⌢2=diagA-2,A-2,…,A-2︷Ne,B⌢2=diagB-2,B-2,…,B-2︷Ne,C⌢=diagC-,C-,…,C-︷Ne,where Ne is the total number of elements. The potential energy of the element springs in the whole structure is given as(27)Vs=12w⌢T∑i=14C⌢K⌢ti+∑i=12A⌢1TC⌢K⌢rxiA⌢1+∑i=12B⌢1TC⌢K⌢ryiB⌢1w⌢,where (28)K^ti=diagK-ti1,K-ti2,…,K-tiNe.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(29)Vct=14w⌢T∑i=14C⌢K⌢cti+∑i=12A⌢1TC⌢K⌢crxiA⌢1+∑i=12B⌢1TC⌢K⌢cryiB⌢1w⌢.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) (30)U=12w⌢TDA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2+∑i=14C⌢K⌢ti+12K⌢cti+∑i=12C⌢A⌢1TK⌢rxi+12K⌢crxiA⌢1+B⌢1TK⌢ryi+12K⌢cryiB⌢1w⌢.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](31)we=wme,wmye,wmxe,wmxyei=1,M;j=1,N,wme,wmyei=3,…,M-2;j=1,N,wme,wmxei=1,M;j=3,…,N-2,wmei=3,…,M-2;j=3,…,N-2T,where wmx=∂w/∂x represents the rotation in x direction, wmy=∂w/∂y represents the rotation in y direction, and wmxy=∂2w/∂x∂y represents the torque at the corner grid point. Using the DQ rule, one can obtain the following relation:(32)we=Q-w-e.For more information about matrix Q, one can refer to [24]. The global matrix of the transformation matrix Q is given as(33)w=Q⌢w-(34)Q⌢=diagQ-,Q-,…,Q-︷Ne.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(35)K=12Q⌢-TDA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2+∑i=14C⌢K⌢ti+12K⌢cti+∑i=12C⌢A⌢1TK⌢rxi+12K⌢crxiA⌢1+B⌢1TK⌢ryi+12K⌢cryiB⌢1Q⌢-1M=12ρhQ⌢-TC⌢Q⌢-1.Finally, the vibration problem became eigenvalue problem which can be solved by solving the following equation: (36)K-ω2M=0.


### 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–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.


## 2.1. The Differential Quadrature Rules

The differential quadrature (DQ) method approximates thesth derivative of a one-dimensional function f(x), which is derivable in the interval [a,b], at a point x 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 sth-order derivative of function f(x) can be written as(1)dsfxdxsi=∑j=1NAijsfj,where N is the grid numbers in the interval. Aij(s) is the sth-order weighting coefficient; thus, [Aij(s)] is the sth-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(2)Aij1=∏k=1,k≠iNxi-xkxi-xj∏k=1,k≠jNxj-xk(3)Aijs=sAiis-1Aij1-Aijs-1xi-xjFori≠j,2≤s≤N-1Aiis=-∑j=1,j≠iNAijs1≤s≤N-1.According to the two-dimensional DQ rule, the partial derivatives of functionf(x,y) at grid point (xi,yj) can be written as(4)∂rf∂xrij=∑m=1MAimrfmj,∂sf∂ysij=∑n=1NBjnsfin,∂r+sf∂xr∂ysij=∑m=1MAimr∑n=1NBjnsfmn,where Aim(r) and Bjn(s) are the weighting coefficients along x and y directions, respectively. M and N are the numbers of the grid point in x and y 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]:(5)∂rf∂xrk=∑m=1M×NA-kmrf-m,∂sf∂ysk=∑m=1M×NB-kmsf-m,∂r+sf∂xr∂ysk=∑m=1M×NF-kmrf-m,where k, m=(j-1)M+i, (i=1,2,…,M;j=1,2,…,N) and(6)f-=f11⋯fM1f12⋯fM2⋯f1N⋯fMNTA-r=Ar0⋯00Ar⋯0⋮⋮⋱⋮00⋯Ar,B-s=B11sB12s⋯B1NsB21sB22s⋯B2Ns⋮⋮⋱⋮BN1sBN2s⋯BNNs,where A-(r) and B-(s) are both (M×N)×(M×N) matrix, representing the weighting coefficient matrices related to the derivative in x and y direction, respectively, in which(7)Ar=A11rA12r⋯A1MrA21rA22r⋯A2Mr⋮⋮⋱⋮AM1rAM2r⋯AMMr,Bijs=Bijs0⋯00Bijs⋯0⋮⋮⋱⋮00⋯Bijs.Note that the two-dimensional DQ rules can be treated as one-dimensional DQ rule applied both inx and y direction. Thus, Aijs=Bijs, where Aijs represents the sth-order weighting coefficient along x direction of grid point at jth row and ith column. Bijs represents the sth-order weighting coefficient along y direction of grid point at jth row and ith column.


## 2.2. 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 e. 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 ktie, kryie, and krxie, where e means theeth element of the structure, t and r represent the translational and rotational springs, and x and y denote the normal directions of each side of an element. i is the ith side of a square element. For kryie, due to the fact that the only normal directions of side 1 and side 3 are in y direction, thus i=1,2 for the rotational springs, which is the same as krxie 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 kctie,ei, kcryie,ei, and kcrxie,ei, where the superscript e,ei represents the coupling springs connecting theeth element andeith element, ei means the adjacent element is connected to the ith side of elemente, andc indicates the coupling springs. For simplicity, the first set of virtual springs is called element springs, while the other sets are called coupling springs.Figure 1
The coupling relationships between element e and four adjacent elements.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 intoM grid points along the x direction and N grid points along the y direction with two sets of virtual springs attached to the element. The displacement function of the eth element is defined as (8)w-ex,y=∑i=1M∑j=1Mlixljxwije,where li are the Lagrange polynomials and wije=we(xi,yj) is the displacement at the Gauss Lobatto quadrature points [24]. Then the potential energy and kinetic energy ofeth element are given as(9)Upe=D2∫∫Se​∂2w-e∂x22+∂2w-e∂y22+2υ∂2w-e∂x2∂2w-e∂y2+21-υ∂2w-e∂x∂y2dxdyTe=12∫∫Se​ρhw-˙e2dxdy,where D is the flexural rigidity, υ is Poisson’s ratio, ρ is density, and h is the thickness of the plate; for thin plate h/a=0.001, where a 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(10)Vse=12∑i=14∫0liektiew-ie2dlie+12∑i=12∫0liekrxie∂w-ie∂x2dxie+12∑i=12∫0liekryie∂w-ie∂y2dyie.The total potential energy in the coupling springs is given as(11)Vce,ei=12∑i=14∫0liekctie,eiw-iei-w-ie2dlie+12∑i=12∫0liekcrxie,ei∂w-iei∂x-∂w-ie∂x2dxie+12∑i=12∫0liekcryie,ei∂w-iei∂y-∂w-ie∂y2dyie,where the displacements w-iei and w-ie 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(12)∑e=1NeUpe+Vse+12Vce,ei=∑e=1NeTe,whereNe is the total numbers of the element in the structure. According to the DQ rules, the expressions can be transformed into matrix form.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.
(a)
(b)The integrals in the energy expressions are calculated by using the Gauss Lobatto quadrature [24]. The potential and kinetic energy of theeth element is (13)Upe=D2w-eTA-2TC-A-2+B-2TC-B-2+2υA-2TC-B-2+21-υF-2TC-F-2w-eTe=12ρhω2w-eTC-w-e,where F-(2)=A-(1)B-(1) and the displacement vector is (14)w-e=w11e⋯wM1ew12e⋯wM2e⋯w1Ne⋯wMNeT,which can be illustrated clearly through Figure 2(b). Matrix C=diag(C1,C2,…,CM) is a M×M 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 (15)C-=diagC1C1,…,CMC1,C1C2,…,CMC2,…,CMCN.The matrix form of the potential energy of element springs is(16)Vse=12w-eT∑i=14C-K-tie+∑i=12A-1TC-K-rxieA-1+∑i=12B-1TC-K-ryieB-1w-e,where the stiffness matrices of the translational springs that belong to each side of the element are given as(17)K-t1e=diagkt1e,kt1e,…,kt1e︷M,0,…,0︷N-1MK-t3e=diag0,…,0︷N-1M,kt3e,kt3e,…,kt3e︷MK-t2e=diagKt2e,Kt2e,…,Kt2e︷N,Kt2e=diag0,…,0︷M-1,kt2eK-t4e=diagKt4e,Kt4e,…,Kt4e︷N,Kt4e=diagkt4e,0,…,0︷M-1.The stiffness matrices of the rotational springs are given as(18)K-rx1e=diagKrx1e,Krx1e,…,Krx1e︷N,Krx1e=diagkrx1e,0,…,0︷M-1K-rx2e=diagKrx2e,Krx2e,…,Krx2e︷N,Krx2e=diag0,…,0︷M-1,krx2eK-ry1e=diagkry1e,kry1e,…,kry1e︷M,0,…,0︷N-1MK-ry2e=diag0,…,0︷N-1M,kry2e,kry2e,…,kry2e︷M.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(19)Vce,ei=12∑i=14w--e,eiTC--K--ctie,eiw--e,ei+∑i=12w--e,eiTA--1TC--K--crxie,eiA--1w--e,ei+∑i=12w--e,eiTB--1TC--K--cryie,eiB--1w--e,ei,where the sign “=” represents that the matrix or the vector consists of matrices or vectors of two adjacent elements, which can be expressed as(20)w--e,ei=w-ew-eiT,C--=diagC-,C-,A--=diagA-,A-,B--=diagB-,B-.The matrices of the stiffness of the coupling springs are given as(21)K--ctie,ei=K-ctie-K-ctie,ei-K-ctiei,eK-ctiei,K--crxie,ei=K-crxie-K-crxie,ei-K-crxiei,eK-crxiei,K--cryie,ei=K-cryie-K-cryie,ei-K-cryiei,eK-cryiei.For matrixK--ctie,ei, the diagonal elements K-ctie and K-ctiei are only related to theeth andeith element, respectively. The nondiagonal elements -K-ctie,ei and -K-ctiei,e represent the coupling effects between two elements. The elements distributions in the sparse matrix for diagonal elements K-ctie and K-ctiei are the same as that in (17). In order to illustrate the matrix -K-ctie,ei, we can assume that the 3rd side of element e is connected to the 1st side of elementei; then the matrix -K-ct3e,e3 can be expressed as(22)-K-ct3e,e3=-00⋯0⋮⋮⋱⋮00⋯0kct1e,e10⋯00kct1e,e1⋯0⋮⋮⋱⋮00⋯kct1e,e1︸M0⋯0⋮⋱⋮0⋯00⋯00⋯0⋮⋱⋮0⋯0︸N-1M.One can observe that, in the coupling matrixK--ct3e,e3, the nonzero elements in -K-ct3e,e3 are in the same rows as those in matrix K-ct3e and in the same columns as those in matrix K-ct3e3. Thus, the matrix -K-ct3e3,e can be easily understood which is given as (23)-K-ct3e3,e=-0⋯00⋯0⋮⋱⋮0⋯00⋯0⋮⋱⋮0⋯0︷N-1Mkct1e,e10⋯00kct1e,e1⋯0⋮⋮⋱⋮00⋯kct1e,e100⋯0⋮⋮⋱⋮00⋯0︷M.So the nonzero elements in the coupling matrices-K-ct3e,e3 and -K-ct3e3,e are symmetrically distributed about the diagonal line of matrix K--ct3e,e3. For matrices K--crxie,ei and K--cryie,ei, the nonzero elements are distributed in the same way.


## 2.3. 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 (24)Up=D2w⌢TA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2w⌢(25)T=12ρhω2w⌢TC⌢w⌢,where (26)w⌢=w-1w-2⋯w-NeTA⌢2=diagA-2,A-2,…,A-2︷Ne,B⌢2=diagB-2,B-2,…,B-2︷Ne,C⌢=diagC-,C-,…,C-︷Ne,where Ne is the total number of elements. The potential energy of the element springs in the whole structure is given as(27)Vs=12w⌢T∑i=14C⌢K⌢ti+∑i=12A⌢1TC⌢K⌢rxiA⌢1+∑i=12B⌢1TC⌢K⌢ryiB⌢1w⌢,where (28)K^ti=diagK-ti1,K-ti2,…,K-tiNe.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(29)Vct=14w⌢T∑i=14C⌢K⌢cti+∑i=12A⌢1TC⌢K⌢crxiA⌢1+∑i=12B⌢1TC⌢K⌢cryiB⌢1w⌢.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) (30)U=12w⌢TDA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2+∑i=14C⌢K⌢ti+12K⌢cti+∑i=12C⌢A⌢1TK⌢rxi+12K⌢crxiA⌢1+B⌢1TK⌢ryi+12K⌢cryiB⌢1w⌢.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](31)we=wme,wmye,wmxe,wmxyei=1,M;j=1,N,wme,wmyei=3,…,M-2;j=1,N,wme,wmxei=1,M;j=3,…,N-2,wmei=3,…,M-2;j=3,…,N-2T,where wmx=∂w/∂x represents the rotation in x direction, wmy=∂w/∂y represents the rotation in y direction, and wmxy=∂2w/∂x∂y represents the torque at the corner grid point. Using the DQ rule, one can obtain the following relation:(32)we=Q-w-e.For more information about matrix Q, one can refer to [24]. The global matrix of the transformation matrix Q is given as(33)w=Q⌢w-(34)Q⌢=diagQ-,Q-,…,Q-︷Ne.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(35)K=12Q⌢-TDA⌢2TC⌢A⌢2+B⌢2TC⌢B⌢2+2υA⌢2TC⌢B⌢2+21-υF⌢2TC⌢F⌢2+∑i=14C⌢K⌢ti+12K⌢cti+∑i=12C⌢A⌢1TK⌢rxi+12K⌢crxiA⌢1+B⌢1TK⌢ryi+12K⌢cryiB⌢1Q⌢-1M=12ρhQ⌢-TC⌢Q⌢-1.Finally, the vibration problem became eigenvalue problem which can be solved by solving the following equation: (36)K-ω2M=0.


## 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–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.


## 3. 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 Figure3, from which one can observe that the first six frequency parameters λ=ωa2(ρh/D)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 M=N=12 for the following researches, which are sufficient to achieve accurate results. The material parameters of the structures are given as Young’s modulus E=2.1e11 Pa, density ρ=7800 kg/m3, Poisson’s ratio υ=0.3, and thickness h=0.001 m.Figure 3
The first six frequency parameters of a square plate with different numbers of grid points.


### 3.1. 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λ=ωa2(ρh/D)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 kt/D=1e11 as the infinitely large number for the translational springs. Next keep the kt/D=1e11 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 kr/D=1e11 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.Figure 4
The changing rules of the frequency parameters as the stiffness of the element springs increasing.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 Figure5. 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 kct/D=1e11 as the infinitely large number for the translational springs. Keeping kct/D=1e11 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 kcr/D=1e11 as the infinitely large number for the rotational springs. Thus, the elastically coupling conditions can be simulated as those for boundary conditions.Figure 5
The changing rules of the frequency parameters as the stiffness of the coupling springs increasing.


### 3.2. 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 Figure6, where length a=2 m and width b=1 m. 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(37)kt12=kry12=1e11D,kt22=kt32=kt42=kry22=krx12=krx22=0kct42,24=kcrx12,21=1e11D,kct12,21=kct22,22=kct32,23=kcrx22,22=kcry12,21=kcry22,22=0.Similar operations can be applied for plate with different crack types and boundary conditions.Figure 6
A rectangular thin plate divided uniformly into 8 square elements and the serial numbers of each element and side.Figure 7
Geometries of three types of cracks in a rectangular plate and the boundary conditions.Numerical results are presented in Table1, 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.Table 1
Natural frequencies of rectangular plates with several classical boundary conditions and T-shaped, L-shaped, and I-shaped cracks.
Types
Mode numbers
1
2
3
4
5
6
7
8
9
T-shapedCCCC
4.6809
5.2512
6.5963
8.7289
10.3130
14.5327
17.5577
17.9553
19.8488
4.6640a
5.2355
6.5828
8.7118
10.3050
14.5200
17.5550
17.9550
19.8507
L-shapedFCFC
0.4842
1.5996
2.4210
3.0750
3.8014
4.8670
5.8165
6.8436
7.7874
0.4892a
1.5991
2.4196
3.0661
3.7987
4.8616
5.8049
6.8421
7.7821
I-shapedCFFF
0.2382
0.2675
1.3373
1.5138
1.8848
2.1912
3.7463
4.2576
5.7725
0.2146a
0.2930
1.3338
1.5195
1.8827
2.2125
3.7466
4.2565
5.7554
aResults from ANSYS.Figure 8
The first nine mode shapes of the CCCC plate with T-shaped crack.Figure 9
The first nine mode shapes of the FCFC plate with L-shaped crack.In the second example, the rectangular plate with the same geometric parameters is divided uniformly into 8 right trapezoid elements shown in Figure10. 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.Table 2
Natural frequencies of rectangular plates with several classical boundary conditions and irregular T-shaped, irregular L-shaped, and irregular I-shaped cracks.
Types
Mode numbers
1
2
3
4
5
6
7
8
9
T-shapedCCCC
4.2295
6.1142
6.9490
9.2834
10.5213
13.6350
17.1236
18.5392
19.1972
4.2847a
5.8950
6.7348
9.4992
10.5300
13.5580
17.2230
18.6560
19.2360
L-shapedFCFC
0.4945
1.5181
2.2444
3.1847
3.5957
4.4550
5.2570
6.7095
7.6263
0.4871a
1.5894
2.2282
3.0957
3.5863
4.5244
5.1729
6.7097
7.5596
I-shapedCFFF
0.2392
0.2519
1.3157
1.4969
1.7301
2.0860
3.6811
4.2399
5.1381
0.2124a
0.2889
1.3141
1.4914
1.7564
2.0674
3.6928
4.2317
5.1639
aResults from ANSYS.Figure 10
A rectangular thin plate divided uniformly into 8 right trapezoid elements and the serial numbers of each element and side.Figure 11
Geometries of three types of irregular cracks in a rectangular plate and the boundary conditions.Figure 12
The first nine mode shapes of the CCCC plate with irregular T-shaped crack.Figure 13
The first nine mode shapes of the FCFF plate with irregular I-shaped crack.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.


## 3.1. 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λ=ωa2(ρh/D)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 kt/D=1e11 as the infinitely large number for the translational springs. Next keep the kt/D=1e11 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 kr/D=1e11 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.Figure 4
The changing rules of the frequency parameters as the stiffness of the element springs increasing.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 Figure5. 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 kct/D=1e11 as the infinitely large number for the translational springs. Keeping kct/D=1e11 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 kcr/D=1e11 as the infinitely large number for the rotational springs. Thus, the elastically coupling conditions can be simulated as those for boundary conditions.Figure 5
The changing rules of the frequency parameters as the stiffness of the coupling springs increasing.


## 3.2. 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 Figure6, where length a=2 m and width b=1 m. 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(37)kt12=kry12=1e11D,kt22=kt32=kt42=kry22=krx12=krx22=0kct42,24=kcrx12,21=1e11D,kct12,21=kct22,22=kct32,23=kcrx22,22=kcry12,21=kcry22,22=0.Similar operations can be applied for plate with different crack types and boundary conditions.Figure 6
A rectangular thin plate divided uniformly into 8 square elements and the serial numbers of each element and side.Figure 7
Geometries of three types of cracks in a rectangular plate and the boundary conditions.Numerical results are presented in Table1, 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.Table 1
Natural frequencies of rectangular plates with several classical boundary conditions and T-shaped, L-shaped, and I-shaped cracks.
Types
Mode numbers
1
2
3
4
5
6
7
8
9
T-shapedCCCC
4.6809
5.2512
6.5963
8.7289
10.3130
14.5327
17.5577
17.9553
19.8488
4.6640a
5.2355
6.5828
8.7118
10.3050
14.5200
17.5550
17.9550
19.8507
L-shapedFCFC
0.4842
1.5996
2.4210
3.0750
3.8014
4.8670
5.8165
6.8436
7.7874
0.4892a
1.5991
2.4196
3.0661
3.7987
4.8616
5.8049
6.8421
7.7821
I-shapedCFFF
0.2382
0.2675
1.3373
1.5138
1.8848
2.1912
3.7463
4.2576
5.7725
0.2146a
0.2930
1.3338
1.5195
1.8827
2.2125
3.7466
4.2565
5.7554
aResults from ANSYS.Figure 8
The first nine mode shapes of the CCCC plate with T-shaped crack.Figure 9
The first nine mode shapes of the FCFC plate with L-shaped crack.In the second example, the rectangular plate with the same geometric parameters is divided uniformly into 8 right trapezoid elements shown in Figure10. 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.Table 2
Natural frequencies of rectangular plates with several classical boundary conditions and irregular T-shaped, irregular L-shaped, and irregular I-shaped cracks.
Types
Mode numbers
1
2
3
4
5
6
7
8
9
T-shapedCCCC
4.2295
6.1142
6.9490
9.2834
10.5213
13.6350
17.1236
18.5392
19.1972
4.2847a
5.8950
6.7348
9.4992
10.5300
13.5580
17.2230
18.6560
19.2360
L-shapedFCFC
0.4945
1.5181
2.2444
3.1847
3.5957
4.4550
5.2570
6.7095
7.6263
0.4871a
1.5894
2.2282
3.0957
3.5863
4.5244
5.1729
6.7097
7.5596
I-shapedCFFF
0.2392
0.2519
1.3157
1.4969
1.7301
2.0860
3.6811
4.2399
5.1381
0.2124a
0.2889
1.3141
1.4914
1.7564
2.0674
3.6928
4.2317
5.1639
aResults from ANSYS.Figure 10
A rectangular thin plate divided uniformly into 8 right trapezoid elements and the serial numbers of each element and side.Figure 11
Geometries of three types of irregular cracks in a rectangular plate and the boundary conditions.Figure 12
The first nine mode shapes of the CCCC plate with irregular T-shaped crack.Figure 13
The first nine mode shapes of the FCFF plate with irregular I-shaped crack.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.


## 4. 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 restrict the displacement of single element while the second set is called coupling springs which restrict the relative displacement between two adjacent elements. Based on the DQFEM theory, the mass matrices and stiffness matrices involving the spring’s stiffness of the elements are deduced. The vibration results of two rectangular plates with different shaped cracks are presented. Through the examples, the procedures of setting the stiffness values of the springs of an element according to the boundary conditions and the crack’s positions are well illustrated. Three types of cracks with different shapes are easily achieved just by modifying the stiffness values of the corresponding springs. The results are compared with those obtained by ANSYS and the good agreement between the results can be observed, which validates the accuracy and efficiency of the present method. 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.



---

*Source: 2073453-2017-12-07.xml*