A Hybrid Finite Element-Fourier Spectral Method for Vibration Analysis of Structures with Elastic Boundary Conditions

A novel hybrid method, which simultaneously possesses the efficiency of Fourier spectral method (FSM) and the applicability of the finite element method (FEM), is presented for the vibration analysis of structures with elastic boundary conditions. The FSM, as one type of analytical approaches with excellent convergence and accuracy, is mainly limited to problems with relatively regular geometry. The purpose of the current study is to extend the FSM to problems with irregular geometry via the FEM and attempt to take full advantage of the FSM and the conventional FEM for structural vibration problems. The computational domain of general shape is divided into several subdomains firstly, some ofwhich are represented by the FSMwhile the rest by the FEM.Then, fictitious springs are introduced for connecting these subdomains. Sufficient details are given to describe the development of such a hybrid method. Numerical examples of a one-dimensional Euler-Bernoulli beam and a two-dimensional rectangular plate show that the present method has good accuracy and efficiency. Further, one irregular-shaped plate which consists of one rectangular plate and one semi-circular plate also demonstrates the capability of the present method applied to irregular structures.


Introduction
The needs for engineers to accurately predict performance of dynamic systems and to produce optimal designs as well as fast progress in hardware performance and constant decrease in the price of computers all lead to increasing popularity and application of the FEM in engineering.Though the FEM is applicable to problems of any geometric and material properties and boundary conditions, it can be very expensive and time consuming to solve large scale finite element models using a desktop computer, especially when the frequency involved is relatively high.When a parametric study is required in the analysis of some problems, such as for the rotor shaft-oil film bearing-flexible foundation system including several different components, remeshing and reanalysis are needed each time [1][2][3].Additionally, the FEM model would become impractically large if wave propagation characteristics of a structure at higher frequencies are sought [4].
Besides simplifying a model to reduce computational workload, many attempts have been made within the last decades in order to replace or improve the conventional FEM with other methods [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].Some of these are able to solve mid-frequency vibration problems with some degrees of success.Apparently, techniques for solving mid-frequency vibration problems of general boundary conditions at high accuracy are still needed [20].How to obtain an accurate model with fewer degrees of freedom is one of the main issues in structural dynamic analysis.Recently, Fourier spectral method (FSM) which was put forward by Li [21] has been widely applied to various vibration analyses of regular structures such as straight beams [21], multispan beams [22], and in-plane [23] and out-of-plane [24,25] vibration of rectangular plates.At the early stage of this method, the generalized Fourier series are first substituted into the governing differential equations with relevant boundary equations which are derived from the force equilibrium relationship, and then the Galerkin discretization procedure 2 Mathematical Problems in Engineering is incorporated to transform the overall equations into a standard matrix eigenvalue problem.And recently, one more general formulation based on an energy principle where all the unknown coefficients are regarded as generalized independent coordinates [26] has been introduced, using FSM, to vibration analysis of, for example, single rectangular plates with arbitrary nonuniform elastic edge restraints [27], and more complex structures such as elastically coupled rectangular plates [28] and circular cylindrical shells [29].FSM is even applied to vibroacoustic analysis of a threedimensional acoustic cavity with both the displacement of the plate and the sound pressure constructed in the form of Fourier series with supplemental terms [30].As a novel Fourier expansion method, the convergence speed of the FSM is improved substantially compared with the traditional Fourier series method.It is a suitable systematic method for solving a wide range of engineering and scientific problems.However, compared with local methods such as the FEM, as an analytical method which uses global basis functions, the FSM also shows intrinsic difficulties in satisfying boundary conditions for irregular-shaped domains like other global methods [16] such as spectral methods or spectral finite element methods [13,14].
For many problems such as ships and aeroplanes, the domain considered is geometrically complex and the application of the available spectral techniques to such structures would become difficult or even impossible.The main purpose of this paper is to combine the advantages of the FSM and FEM to develop a new hybrid method for problems with partly regular-shaped boundaries.Firstly, the structure is divided into two zones: one with regular boundaries described by the FSM and the other with irregular boundaries discretised by the FEM.At their interface, the coupling between these two zones is implemented by introducing several fictitious springs on the common edges.These two subsystems are coupled together using Hamilton's equation.After these two methods are combined, it can be found that the computational efficiency is enhanced compared with the standard FEM, and limitations that the FSM can hardly be applied to complex structures are also overcome.Sections 2 and 3 give, respectively, the theoretical formulations and numerical examples of one-dimensional (beam) and twodimensional (plate) structures based on the new method.

Governing Equation of Beam Vibration. Consider an
Euler-Bernoulli beam elastically supported at its both ends.The equations of motion as well as the corresponding boundary conditions can be derived based on the energy principle, in which the total potential energy  and total kinetic energy  are written as where  FSM is the beam length,  is the flexural displacement,  is the angular frequency, and , , and  are, respectively, the flexural rigidity, the mass density, and the cross-sectional area;  0 and  1 are the linear spring constants, and  0 and  1 are the rotational spring constants at  = 0 and  =  FSM , respectively.If FEM is used, the equation for free vibration can be written as in which Λ = { 1 ,  1 ,  2 ,  2 , . . .,   ,   }  and   and   are, respectively, the deflection and rotation at the th node.For a beam with elastic boundary conditions, it is not easy to determine its eigenfunctions which can be used as displacement functions.In this study, the displacement function for FSM will be sought in the form of a simple series expansion as where   = / FSM .When cosine functions are used, the subsequent mathematical derivation becomes much easier.However, the elastic boundary conditions cannot be satisfied by such a series and discontinuities in the first and third derivatives occur at the boundaries.Four supplemental functions   are introduced to remove the potential discontinuities at both ends; namely, Actually, the choice of these auxiliary functions is not unique; however, the appropriate construction of these functions can simplify the subsequent mathematical formulation significantly.These functions can be chosen as long as the following conditions are satisfied:   1 (0) =   2 ( FSM ) =   3 (0) =   4 ( FSM ) = 1 and all the rest first and third derivatives of these functions at these two end points are equal to zero.Detailed derivations and explanations can be found in [21,27].In order to find all the unknown coefficients in (3), which actually determines the structural dynamic behaviour, the Hamilton equation is employed here; namely, where  denotes the total potential energy and  the total kinetic energy.Substituting (3) into the Hamilton equation and minimizing it with respect to each unknown Fourier coefficient and   will lead to the following equation in matrix form: where Ψ = { 0 ,  1 , . . .,   ,  1 ,  2 ,  3 ,  4 }  and K FSM and M FSM are the system stiffness and mass matrices which are explicitly defined in Appendix A in supplementary material available online at http://dx.doi.org/10.1155/2014/834834.In (6), it is assumed that the series expansion (3) is truncated to  =  terms in practical calculation; namely,  = 0, 1, 2, . . ., .

Combination of the FEM and the FSM.
After briefly reviewing the main principles of FEM and FSM, development of a hybrid method from these two will be described.As illustrated in Figure 1, a beam structure is divided into two zones: the left is handled by the FEM and the right by the FSM, which are, respectively, named as FEM zone and FSM zone for convenience.To combine these two methods, two types of fictitious springs with linear and rotational spring coefficients of   and   are introduced at the interface where two zones meet.On the FEM side, it is assumed that the beam is divided into (-1) elements, namely, with  nodes.The displacement and rotation at the interface are  FEM  and  FEM  .On the other hand, the displacement and rotation at the interface from the FSM are The potential energy stored in these two springs can now be calculated as When the spring constants of the fictitious springs involved are assigned very big values, the transverse displacements and rotations of the two subdomains at the interface will then be effectively equal.
Substituting ( 8) and ( 9) into ( 5), the coupled matrix K couple can then be obtained.The overall equation of motion can be derived as where the elements of the coupled matrices are given in Appendix B.

Numerical Results and Discussion
. In this section, several numerical examples are presented to show the accuracy of this hybrid method.Firstly, a straight beam with the following parameters is used in the calculation: beam length  =  FEM +  FSM = 1 m , in which  FEM and  FSM represent the length of the FEM and FSM zones on the beam structure; Young's modulus  = 2.1 × 10 11 Pa; the second moment of inertia of the cross-section  =  4 /4 with  = 0.01 m .The Fourier series are truncated to  = 20 for the FSM, while the element number for the FEM is set as 50.In order to study the influence of the relative proportion of these two subdomains on the prediction results, the length of the FSM zone is varied from 0 to 1 with the increment of 0.2 while the overall length of the beam is kept unchanged as  = 1 m .The results are given in Table 1 for a beam with clamped-clamped boundary condition, which can be readily achieved by setting all the end restraint stiffness into a very large number (the translational springs are  0 = 5 × 10 7 ,  1 = 5 × 10 7 while the rotation springs at the boundary are  0 = 5 × 10 7 ,  1 = 5 × 10 7 and the coupled springs are   = 5 × 10 7 ,   = 5 × 10 7 ) for FSM zone side and fixing the end node on FEM zone side.As can be observed, the relative proportion between these two computational zones almost has no effect on the accuracy of the predicted frequencies.Such good results indicate that the proposed hybrid method is effective.Table 2 compares the efficiency and accuracy of the FEM, the FSM, and the hybrid method.All these three methods are implemented on the same computational platform in MATLAB.A relative high frequency range up to the 20th mode is computed by these three approaches.Firstly, FEM and FSM are used to make prediction and it can be seen that FSM is much more efficient than FEM, since its dimension of final system matrix is much smaller than that of FEM.For the proposed hybrid method, it can be observed again that the prediction is accurate enough for capturing the modal characteristics of this free-free beam.
In the first three lines of Table 2, we see that in these three cases the accuracy class of the 20th modal frequency calculated is equal.The degrees of freedom of the matrix used for calculation are 85, 600, and 345 for FSM, FEM, and current hybrid method, respectively.Although the CPU times are not a good measure of an algorithm's performance as they do not provide much insight about an algorithm's behaviour and they are still presented here as references which are 0.016(s), 9.235(s), and 7.66(s), respectively, for the first three cases.Compared with FEM, since the FSM with rapid convergence and high efficiency is incorporated into this structure, the overall degrees of freedom will be reduced considerably, say, from 600 to 345 for the first three cases.The CPU time consumed is also reduced compared with pure FEM.The remaining cases in this table are all the same: the current method would reduce DOF significantly compared with pure FEM while the accuracy class of the frequency calculated is still kept in the same level.

Governing Equation of Rectangular
Plates.Similar to the framework for vibration analysis of beam structures, vibrating plate structures with elastic edge supports can also be dealt with by the energy principle.For a rectangular plate with width  and length , the total potential energy  and total kinetic energy  are given by 2 ) d d 2 ) d d where  is the flexural displacement field, and , ], , and ℎ are, respectively, the flexural rigidity, the Poisson ratio, the mass density, and the plate thickness. 0 and   ( 0 and   ) are the constants of linearly distributed springs and  0 and   ( 0 and   ) the rotational spring constants at  = 0 and  =  ( = 0 and  = ), respectively.In the FEM description of a plate structure, each node has a lateral displacement and two rotations.The equation for vibration analysis of a plate structure can also be obtained in the form of (2).
In this study, the displacement function of a rectangular plate for the FSM will be sought in the form of series expansions as follows [24]: where   = /,   = /, and    (or    ) represent a set of closed-form supplementary functions, defined as The form of    is much the same as    except that  should be replaced by  while  should be replaced by .A procedure similar to that for the beam in Section 2 can be established to determine all the coefficients in (12), in which the series is truncated to  =  and  = .

Combination of the FEM and the FSM.
Compared with the hybrid beam model, the coefficients which couple these two methods are a little different.As illustrated in Figure 2, a plate structure is divided into two zones: the left is handled by the FEM and the right by the FSM.Geometrically, in the coupled beam model, the two different beam zones are connected point to point, while for the coupled plate the two zones are connected line to line.As the rectangular plate is one type of structure which possesses symmetric characteristics, the edge on the left side of the FSM plate is considered here.For this paper, the model for FSM is based on classical Poisson-Kirchhoff theory, which requires C 1 continuity in the framework of the finite element method [1].Certain types of elements based on Poisson-Kirchhoff theory have compatibility problems as the rotational displacements would become discontinuous between adjacent elements on the edge.Hence, for the sake of generality, the transverse springs between these two zones are set to be linearly distributed while the rotation springs are concentrated to nodes.For the FSM zone, the displacement and rotation of the edge at  = 0 can be expressed as For the   element bordering the interface, the four nodes of the element are defined as   ,   ,   , and   .
And then the potential energy stored in linear spring can now be calculated as where   and   are the coordinates of the nodes    and    , respectively, which are located on the interface line.
The type of element used here is a 4-node incompatible plate element, in which the normal derivative at the edge is discontinuous.As a result, the rotation displacement at the edge cannot be obtained directly from derivatives of (12).For overcoming difficulty, one procedure is introduced here taking care of the uniformly distributed rotational spring along the interface line: they are replaced by a certain number of discrete point-to-point rotational springs between each pair of the nodes having the same spatial coordinates on either side of the interface line.Suppose that the number of the pairs of such nodes is   .For these discrete rotation springs, the potential energy stored in them can be expressed as where    represents the -coordinate of the   nodes on the interface line and     is the rotations about the -axis at the   nodes.Substituting (17) into (18) and (19) and then setting the derivatives of the resultant equations (with respect to each unknown Fourier coefficient and unknown displacement   and   ) to zero, the coupled matrix K couple can be obtained.In this way, the overall system equation can be written as in which As the detailed coupled matrices can also be easily obtained according to the same procedures the beam coupled matrices are obtained and details of ( 20) are not given here.

Numerical Results and Discussion
. In this section, the accuracy and effectiveness of the proposed hybrid method will be demonstrated for predicting the modal behaviour of a 2-D plate structure with various boundary conditions and/or different shapes.As the first example, a rectangular plate with simply supported boundary condition at all its edges is considered.For the FSM, a simply supported edge can be treated as a special case while the stiffness coefficients for the translational springs become infinitely large and the rotational spring becomes zero.The plate dimension is 5 m wide and 10 m long with thickness ℎ = 0.008 m; the material parameters used for the plate are  = 2.1 × 10 11 Pa and  = 7800 Kg/m 3 .As mentioned above, in numerical calculations the series expansion in (13) has to be truncated to  =  = 12.The plate is divided by a straight line at the centre of the length into two square zones with all the width and length equal to 5 m.For the FEM model, 20×40 elements are chosen so that higher (such as the 8th, 9th, and 10th) frequencies given in Table 3 would then become accurate enough for comparison.Table 3 shows the first ten natural frequencies of a rectangular plate simply supported at all its edges (the translational springs are  0 = 0,   ,  0 ,   = 5 × 10 7 while the rotation springs at the boundary are  0 = 0,   ,  0 ,   = 5 × 10 7 and the coupled springs  l = 5 × 10 7 ,   = 5 × 10 7 ).Several mode shapes are also given in Figure 3.It is possible that the matrix K becomes ill-conditioned if the coupled springs,  l and   , are set to very large numbers.In Figure 4, the first 6 mode frequencies of the rectangular plate are given with the coupled springs,  l and   , varied from 0 to 10 10 .It can be clearly seen from this figure that the influence of the coupled springs on the frequencies calculated is minimal when the springs changed within a wide range from 10 6 to 10 10 .
In order to verify the accuracy and stability of the proposed hybrid method thoroughly, frequencies of higher modes are provided in Table 4.It can be seen clearly from this table that the hybrid method is slightly more accurate than the FEM compared with analytical results.However, because the FEM is known to be incapable of producing good results at very high frequencies and the hybrid method uses the FEM, it is expected that the hybrid method will produce more accurate results of high frequencies than the FEM but still cannot produce very accurate results of very high frequencies.
Another example of a rectangular plate (identical to the previous FSM rectangular plate with its width and length equal to 5 m) plus a matching semicircular plate (with its diameter equal to the width of the rectangular plate) is also solved and the results are presented in Table 5 with their corresponding mode shapes plotted in Figure 5.For the semicircular plate, the finite element matrices for bending are directly extracted from the ANSYS software based on the shell 63 element.As shown in Figure 5, elements used in this example are irregular-or even triangular-shaped.And the results show the generality of the present method when using finite elements with irregular shapes.
It can also be seen from this example that the present method can be applied to structures with geometrically irregular boundary conditions.In these two 2-D examples, the validity and accuracy of the presented hybrid method are also confirmed by comparison with FEM results.

Conclusions
In this paper, a hybrid FEM-FSM modelling method for the vibration analysis of structures with elastic boundary conditions is proposed.The key point of this method is to divide the whole structure into FEM and FSM zones, and then fictitious springs are introduced to couple the two zones at their interface within the framework of energy principle.Comparative studies of the hybrid method with the FSM and FEM are carried out for modal characteristics, which confirm both the validity of the approach and its efficiency.The number of the degrees of freedom can also be significantly reduced compared with the conventional FEM.Although the concept of such a hybrid method is demonstrated only for a 1-D structure (Euler-Bernoulli beam) and two 2-D structures (plate) in this paper, the proposed approach can easily be extended to more complicated 3-D structures such as box-type structures, as well as structures with holes or notches.

Nomenclature
: Total potential energy : Total kinetic energy : Flexural rigidity : Young's modulus : Second moment of inertia of the cross-section ]: Poisson's ratio

𝑤:
Flexural displacement : Rotational displacement  0 ,  1 : Linear spring constants at  = 0 and  = LFSM, respectively  0 ,  1 : Rotational spring constants at  = 0 and  = LFSM, respectively ,  FSM ,  FEM : Total beam length, beam length of the FSM zone, and FEM zone, respectively : A n g u l a r f r e q u e n c y : Mass Total number of nodes on the interface line   : F r e q u e n c y p a r a m e t e r  0 ,   : Translational types of stiffness at  = 0 and  = , respectively   ,   : Translational types of stiffness at  = 0 and  = , respectively  0 ,   : Rotational types of stiffness at  = 0 and  = , respectively  0 ,   : Rotational types of stiffness at  = 0 and  = , respectively  l : Fictitious linear spring constants for plate   : Fictitious rotational spring constants for plate K: S t i ff n e s s m a t r i x M: M a s sm a t r i x .

Figure 1 :
Figure 1: Schematic of elastically restrained beam structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.

Figure 2 :
Figure 2: Schematic of elastically restrained plate structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.

Figure 3 :
Figure 3: Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangular.The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) tenth mode shapes.

Figure 4 :
Figure 4: Influence of coupled springs on the former 6 mode frequencies of the rectangular plate.

Figure 5 :
Figure 5: Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangularcircular plate.The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) twelfth mode shapes.

Table 1 :
Effect of the division ratio of the two computational domains on the modal frequencies of the beams structure with clamped-clamped boundary condition.Mode order  = /(  √/EI) 1/2 ,  = 20,   = 50

Table 2 :
Comparison of computing time of the hybrid method with pure FEM and FEM with different truncated Fourier terms and/or FEM element numbers, in which a completely free-free beam structure is used.

Table 3 :
Comparison of the first ten natural frequencies of the rectangular plate.: a FEA results are calculated by using ANSYS with 20 × 40 elements.b Results with  =  = 12 at the FSM zone and 20 × 20 elements at the FEM zone. Note

Table 4 :
Comparison of higher natural frequencies of the rectangular plate.
Note: a FEA results are calculated by using ANSYS with 20 × 40 elements.b Results with  =  = 12 at the FSM zone and 20 × 20 elements at the FEM zone.

Table 5 :
Comparison of the first ten natural frequencies of the rectangular-half circular plate. = 12 at the FSM zone and 22 elements at the edge of the coupled line and the element matrix of the semicircular plate were derived from the ANSYS software.
Note: a FEA results are calculated by using ANSYS with 773 elements.b Results with  = Double Fourier series coefficients for plate   ,   : Single Fourier series coefficients for plate   : N u m b e ro f e l e m e n t s   ,   ,   ,   : Four nodes number of element     ,   : Coordinates of the nodes   and   in direction   : Index of nodes on the interface line   :