On the modeling of beam reinforced thin plates using the spectral element method

Modeling beam reinforced thin plates at mid and high frequencies through the most commonly used methods such as finite and boundary element methods frequently leads to unsatisfactory results, since the accuracy of these methods depends on the relation between the dimensions of the elements in which the structure was discretized and the wavelength. Due to this characteristic, the modeling using these techniques will require that the size of the elements becomes smaller as the frequency increases, while its number needs to be increased. For structures that are usual in some areas, like the aerospace industry, this will be possible only with an unreasonable computational effort, which is responsible for restricting the use of these methods practically to low-frequency applications. Semi-analytical methods such as the spectral element method do not need mesh refinement at higher frequencies, but they were very limited in the geometries and boundary conditions that can be treated. This paper presents a spectral element for rectangular thin plates reinforced symmetrically along the sides with Euler beams, which can be used to model plates with arbitrary boundary conditions. The method was verified by comparing its results with those obtained from a Finite Element model.


Introduction
Reinforced thin plates are widely used nowadays, for instance in the aerospace industry, where reinforced panels are assembled to form aircraft fuselage and launch vehicle fairing structures.The attempt to obtain reliable models for dynamic simulations of this kind of mechanical structure is mostly done by using the Finite Element Method (FEM) [1] and the Boundary Element Method (BEM) [2].A main drawback of this approach when trying to simulate the structural behavior at mid and high frequencies is that it will demand an increase on the number of elements associated with a reduction in their dimensions as the frequency range to be analyzed increases, what leads to the need of more powerful computational resources.Even for mid-frequency ranges the models obtained with these techniques become already prohibitively large.This characteristic of element-based methods limits their use mostly to low frequency applications.
An alternative to element-based methods is the Statistical Energy Method (SEA) [3], a probabilistic technique capable to predict the averaged energy levels.The good accuracy provided by this approach for high frequencies is fair, but it is unable to give results at different locations on the domain and, in the mid-frequency range, some of the assumptions that are necessary for its application do not hold, such as high modal density.
At the same time, other methodologies based on an indirect Trefftz approach [4,5] have been recently developed.These methods, which can be classified as wave-based methods, do not require mesh discretization to model a domain with constant geometric and physical properties, since the pressure or displacement fields are described by wave functions that exactly satisfy the differential equation of the problem.The solutions, obtained as an infinite series truncated accordingly to the desired precision, are able to describe an infinite number of modes and are obtained by determining the unknown contribution of the wave factors, what is done by introducing the boundary conditions of the problem.The matrices produced are smaller than the ones from FEM and BEM, and in spite of the fact that they are fully populated and frequency dependent, it has been proven that the methods are computationally more efficient for the analysis of steady-state vibroacoustic problems.A comprehensive overview of the methodologies used in the wave based methods is presented in reference [6].
Among the wave methods, it is worth mentioning the Image Method [7,8], the Superposition Method [9] and the Spectral Element Method (SEM), also known as the Continuous Element Method, which is the subject of this work.In the field of spectral methods.Kulla [10] presented a high precision finite element method, which was able to model beams and plates with arbitrary boundary conditions.The same approach was used by Kevorkian and Pascal [11] and Casimir et al. [12] who called it the continuous element method.Lee and Lee [13] applied the spectral element method to model Levy type plates and Doyle [14] gave a Fourier approach to it.Arruda et al. [15] extended the work of Lee and Doyle, developing a spectral element for reinforced Levy type plates and Campos and Arruda [16] presented a spectral element capable to model plates with any boundary conditions.
In this paper, the spectral element method is extended to treat thin plates reinforced with Euler beams and with any boundary condition.The dynamic stiffness matrix for a spectral reinforced plate element is developed and the problem is solved both for the homogeneous and forced cases.A numerical example is developed to demonstrate the accuracy of the method and the results are compared with those obtained with FEM.

The Kirchhoff plate equation solution
In this section, the Fourier series solution for the differential equation of a Kirchhoff plate is used to obtain the matrices of displacements and forces, dependents on x, y and the circular frequency ω, at the boundaries of the plate.From these rectangular matrices it is possible, by a new series expansion, to obtain square matrices dependent only on ω, which combined will result in the dynamic stiffness matrix of the problem.
Starting from the governing equation in the frequency domain of the Kirchhoff plate, which can be expressed as with where w is the transversal displacement, P 0 is the intensity of a punctual load applied to position (x 0 , y 0 ), ρ is the material density, h the plate thickness, ν the Poisson coefficient and E the Young modulus.Assuming now that there is no loading in the domain of the plate, i.e., all excitations will be introduced as boundary conditions, the problem can be addressed by the homogeneous form of Eq. ( 1), which has a solution of the form w(x, y; ω) = C e p x e q y (4) with C, p and q to be determined.By introducing this solution on Eq. (3), it is obtained its characteristic equation There are infinite values of p and q that satisfy Eq. ( 5).Considering a rectangular plate with dimensions L x =2a and L y =2b and assuming that the solution presented in Eq. ( 4) can be expanded in the x direction as an exponential Fourier series, it will be obtained Introducing the expression for p m into Eq.( 5), two expressions for q m will be obtained and, therefore, a given m will yield eight basis solutions for Eq. ( 1), which grouped into a set can be expressed as where the terms C 1 m , C 2 m , . . ., C 8 m , are unknown constants to be determined.
Applying the same approach to obtain the solution in the y direction, another set of eight basis solutions is obtained as Therefore, the general solution for the homogenous differential Eq. ( 1) is whose explicit expression for a given n 1, with the Fourier series in trigonometric form and with new constants C 1 m , C 2 m , . . ., C 16 m in which are contained the constants presented in Eqs ( 8) and ( 9) and the imaginary number i that multiplies the sine functions, is and, for cosine functions, Defining Equation (11) becomes For n = 0, ψ 0 = e i k y e −i k y e k y e −k y e i k x e −i k x e k x e −k x T (18) Assuming now n = 0, 1, 2, . .., it can be defined and a general expression for the displacements, in vector form, will be given by

Spectral stiffness matrix for an element with beam reinforcement
In order to develop a spectral element reinforced with Euler beams, the contribution of the beam will be taken into account by considering it ideally attached to the plate [16], i.e., the displacements and rotations of the beam and the edge of the plate where it is attached are exactly the same at all points along the beam.In this case, the moments and the effective shear force at the boundaries of the plate, shown in Fig. 1 for the x direction, will equal those in the beam.
The beam will be subjected to both flexural and torsional motion and the equations of motion, for the case presented in Fig. 1, have the form where for the beam w(a, y) is the transverse displacement, θ(a, y) is the torsion angle for a beam, EI b and GJ b are the flexural and torsional rigidities, respectively, and ρA and ρI p are the mass and the polar mass moment of inertia per unit length, respectively.The thin plate shear forces and moments can be expressed as functions of the displacement, using the well known relations From Eqs (23) and ( 24), a modified boundary condition can be established to take into account the effective shear force and the moment including the beam effects for the reinforcement placed at x = −a or x = a, as and for the reinforcement placed at y = −b or y = b as The slope of the thin plate will be given by the usual expressions Introducing Eq. ( 22) in Eqs (26) to (35) will make it possible to derive a spectral plate element with beam reinforcement by writing, in a matrix form, the displacements and forces at the boundaries along the x and y directions.This can be done simply by introducing in the expressions of displacements and forces the values of x or y of the respective boundary and, after that, arranging them as shown in Eqs (36) and (37), which will yield the boundary displacement and force vectors d and f , given by which can be written as a product of rectangular matrices of displacement and force coefficients, and, by the vector of unknown constants c, Eq. ( 21), as The matrices D and F have the form presented in Eq. (39).

D(x, y)
and symbolic computation should be used to determine the values of each matrix element.Equations ( 40) and ( 41) show examples of the expressions obtained in this process for the case of just one reinforcement beam placed along y = b.
In order to make the matrices D and F square, they and the vectors d and f will be expanded in a trigonometric Fourier series truncated at an adequate number of terms m equal to n, and the sine and cosine coefficients of each expanded term will be placed in different rows.This procedure results in the square matrices Ds and Fs associated with the expanded vectors ds and fs , and has the additional advantage of making it possible to obtain a new relation between these terms without the dependence on x and y, which is done by assuming that where D, F, d, and f, all independent of x and y, are the coefficient matrices and vectors, respectively, for displacements and forces at the boundary of the plate.The matrix T, which depends on x and y, is a diagonal matrix which contains the trigonometric terms of the Fourier series expansion and has the form and T 0 is an 8 × 8 identity matrix.Introducing Eq. (42) in Eq. (38) will produce two expressions independent of x and y given by Eliminating c from the expressions in Eq. ( 45), the spectral stiffness matrix S will be obtained as For plates with all edges with homogeneous force boundary conditions or with all edges clamped, the natural frequencies and mode shapes can be obtained by setting respectively the vectors f or d equal to zero in Eq. ( 47) and solving the resulting eigenproblem.

Case of inhomogeneous boundary conditions
If forces or displacements are imposed to the plate at the boundaries, the displacement field can still be obtained with the homogeneous formulation by introducing the boundary displacements and forces in vectors d and f, respectively.
The boundary conditions can be imposed by expanding them in the same Fourier series used in obtaining the dynamic stiffness matrix and inserting the coefficients obtained in the corresponding positions of vectors d or f.For the particular case of a punctual unitary force -considered as a boundary condition -applied to an edge along the y direction and at an arbitrary position y = y 0 , which can be expressed as P = P 0 δ (y-y 0 ), the expansion in a Fourier series will produce coefficients given by the expressions Introducing these coefficients in f and imposing the boundary conditions on Eq. (47), in the same way as it is done in the finite element method, the value of d can be obtained which, introduced in Eq. (50), will completely determine the displacement field as

Numerical validation
In real cases, the reinforced structures to be modeled have periodicity in the reinforcement, which can be in one or two directions.In order to save computational time it should be taken advantage of this periodicity by modeling basic cases of reinforced plates which can be assembled in the same way as it is done for FEM elements.In the following numerical example, the reinforcement is placed only in one edge of the plate, which makes it suitable for modeling panels with reinforcements just in one direction.If a panel with reinforcements in two directions was to be modeled, it could be done by assembling a basic case with two reinforcing beams concurring at one corner.
In this section a numerical implementation of the reinforced SEM is presented to validate its formulation.A rectangular plate with free-free-free-free boundary conditions was modeled with SEM and the results were compared with those obtained with a FEM model.The plate and the reinforcement beam have the dimensions shown in Fig. 2 and the following properties: E = 69 GPa, ρ = 2700 kg/m 3 and ν = 0.3.
As the purpose of this example is only to compare the modeling methods, it is assumed that for both methods the area moment associated with the torsional rigidity, J b , can be approximated by where I b is the beam flexural rigidity.Expressions for J b can be found, for instance, in reference [17].
The polar moment of area of the beam, I p , is given by A unitary punctual harmonic load P = δ(y-y0) was applied at x 0 = −0.15 and y 0 = 0.00 and the frequency response functions, measured at the same point, were obtained for both methods.Twenty five terms of the Fourier expansion were used for the SEM and 8100 four-node quadrilateral plate elements and 90 two-node beam elements were used for the first finite element model, FEM1, and 32400 four-node quadrilateral plate elements and 180 two-node beam elements for the second finite element model, FEM2.Convergence was assumed to be reached when visually the FRFs stop changing with the refinement of the mesh and the increase in the number of Fourier terms for FEM and SEM models, respectively.The results obtained were compared in the frequency ranges 0-1 kHz, 7-8 kHz and 14-15 kHz, as shown in Figs 3, 4, and 5.The agreement between then in the two lower frequency range is quite good and in the higher frequency range the FEM results approach the SEM solution as the mesh is refined.In terms of the computational cost, SEM presented a better performance, even with the FEM code being a commercial package and the SEM code being neither optimized nor compiled.

Conclusions
A spectral element for thin plates which can be used to model beam reinforced plates with arbitrary boundary conditions was developed, and a detailed description on how to obtain all the terms needed to implement it was presented.The results obtained using this reinforced spectral plate element proved to be appealing and its accuracy makes it a potential tool for structural analysis of aerospace panels in mid and high frequency ranges, with a smaller model and a computational cost very low compared to FEM.The use of symbolic computation is essential in obtaining the expressions contained in the matrices of displacements and forces.Further development of this method is necessary to allow its application for any kind of load in the domain and for the modeling of domains other than rectangular.

Fig. 1 .
Fig. 1.(a) Reinforced plate; (b) shear forces and moments acting on the plate and on the reinforcing beam.