The strong formulation finite element method : stability and accuracy

The Strong Formulation Finite Element Method (SFEM) is a numerical solution technique for solving arbitrarily shaped structural systems. This method uses a hybrid scheme given by the Differential Quadrature Method (DQM) and the Finite Element Method (FEM). The SFEM takes the best from DQM and FEM giving a highly accurate strong formulation based technique with the adaptability of finite elements. The present study investigates the stability and accuracy of SFEM when applied to 1D and 2D structural components, such as rods, beams, membranes and plates using analytical and semi-analytical well-known solutions. The numerical results show that the present approach can be very accurate using a small number of grid points and elements, when it is compared to standard FEM.


INTRODUCTION
his manuscript illustrates a general view of a class of methods which approximates functions at points, called collocation or sampling points.The complete review of the present methodology can be found in the previous works by the authors [1][2][3][4][5][6] and in the book by the authors [7].The name which comprehends all this approaches is termed Strong Formulation Finite Element Method and it has its basis in the Differential Quadrature (DQ) Method (DQM) [8][9][10][11][12] as far as the numerical solution is concerned and in the Finite Element concept [13] as far as element connection is considered.The DQ method is part of a more general family of methods called Spectral Methods (SMs) [14][15][16] with the advantage of considering a free collocation once certain basis functions are set.However, DQM is not always stable increasing the number of collocation points.For this reason Generalized Differential Quadrature (GDQ) method was introduced by Shu [17].Shu followed the former studies of Quan and Chang [10,11], who developed a mathematical procedure for the definition of the weighting coefficients in exact form, thus the accuracy and stability is not lost increasing the number of grid points.Other papers that used the mapping technique for strong form based formulations can be found in [18,19].It should be mentioned that, the authors in their previous papers gave explicit formulae for the corner point conditions [1][2][3][4][5][6] that are an extension of the conditions given by other researchers.The following sections expose the mathematical basics of the DQ and GDQ methods.Only some details are given for the sake of conciseness, nevertheless for further details the interested reader can read the references [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34].In the numerical applications, the stability, reliability and accuracy of the Strong Formulation Finite Element Method (SFEM) and the Weak Formulation Finite Element Method (WFEM) are shown through several figures.

DIFFERENTIAL AND INTEGRAL QUADRATURE
enerally an unknown sufficiently smooth function   f x can be approximated by a set of basis functions , where N is the total number of collocation points in a closed definition interval.A polynomial set uniformly converges to the unknown function when the number of grid points tends to infinity and the unknown function is smooth in a closed interval.Hence, the approximate solution of a function   These polynomials constitute a linear vector space N V and are used for the approximation (1) [7,17].A list of several basis functions is given in Tab. 1 and it should be used as a reference for the present manuscript.It should be mentioned that Tab. 1 presents the polynomials in their definition interval indicated by r .The functional approximation is performed using the Cartesian coordinates, as indicated by Eq. ( 1).The approximation can be done if the domain is discretized in N discrete points, such as k x for , , , Tab. 2 reports the most common grid point distributions that can be found in literature.It is assumed that every grid is defined in the interval . In addition, some cases need extra points for the enforcement of the boundary conditions, known also as  -point technique [9, 12]).Hence, a grid distribution whereas considering  -points it takes the following aspect ..., 1 The stretching formulation [7,17] can be defined using the previous definitions ( 2)-( 3), a grid distribution 2) and a non-zero constant  , as This technique takes the points in the domain and stretches them into the domain . The parameter  cannot take any value because for some entries the distribution (see for more details [140]).Finally, the collocation points in dimensionless form , thus, a general coordinate transformation [7] can be used as where is a dimensionless discretization.Using the so-called differentiation matrix procedure [7,10,11,17], the approximation (1) for the one-dimensional case can be written in matrix form as where , , , is the vector of the unknown function values, λ is the vector of the unknown coefficients j  and the components of the coefficient matrix A are given by where  is the domain length, the transformation (5) corresponds to , for 1, 2,..., Expression (1) can be derived for the general n -th order polynomial and the derivative is transferred to the functions   j x


, because the unknown coefficients j  do not depend on the variable x Expression ( 8) can be written in matrix form as follows where is the vector which contains the derivative of the function   f x at all the discrete points i x .In order to perform the derivation, the matrix A should be invertible.This property depends on the basis functions  and on the location of the grid points.The unknown parameters vector λ can be evaluated from Eq. ( 6) by simply inverting the matrix A , such as Subsequently, substituting expression (10) into (9), one obtains The differentiation matrix   n D is found as a matrix product between the inverse of the matrix A containing the values of the basis functions x  in all the grid points and the derivatives of the same functions     All the steps presented above are valid for any derivative order n .In general, Eq. ( 11) takes the form In conclusion, a generic n -th order derivative can be expressed by the following relation Although Eq. ( 12) is valid for any basis function and grid point distribution, the coefficient matrix A , since it is like the Vandermonde matrix, can become ill-conditioned when the number of grid points N is large.It is important to note that, when the Lagrange polynomials   j l x , Lagrange trigonometric polynomials   j S x or the Sinc function   j Sinc x are chosen as a basis of the linear vector space N V , the coefficient matrix A becomes an identity matrix This is due to the fact that, the three previous basis functions have the following properties In all these three cases Eq. ( 12) becomes The important consequence that derives from relations ( 14) and ( 15) is that the matrix ) and thus the ill-conditioning drawback does not occur when the Lagrange polynomials   j l x , Lagrange trigonometric polynomials   j S x or the Sinc function   j Sinc x are chosen.Further details about the weighting coefficients of the previously reported methodologies can be found in [7].a) b) Figure 1: a) Static analysis for single element structural components varying the number of points N for a CC Rod, SS EB beam with  -points (1/2) and a CC Tim beam using PDQ basis functions and Che-Gau-Lob grid distribution.b) Static analysis for structural components varying the number of elements e n with 7 N  for a CC Rod, SS EB beam with  -points (1/2) and a CC Tim beam using PDQ basis functions and Che-Gau-Lob grid distribution.a) b) Figure 2: a) Static analysis for structural components using Weak Formulation (WFEM) varying the number of elements e n with 7 N  for a CC Rod and a CC Tim beam using PDQ basis functions and Leg-Gau-Lob grid distribution.b) Dynamic discrete spectra for single element structural components with 151 N  for a CC Rod, a SS EB beam with  -points (1/2) and a CC Tim beam using PDQ basis functions and Che-Gau-Lob grid distribution.

Uniform (Unif)
Chebyshev-Gauss-Lobatto (Che-Gau-Lob) Quadratic (Quad) (only for M odd) Extended Chebyshev (Extd Cheb) or (Cheb I) 1, 2,..., , 1,1 Jacobi-Gauss (Jac-Gau) Ding et al. [37] distribution  However, when a domain decomposition method is taken into account, it could not be necessary to use a polynomial of high degree inside each element, because a good approximation might be captured with a number of points less than 13 N  .It has been demonstrated in [8,9] that the Vandermonde matrix becomes ill-conditioned when 13 N  .Thus, if a limited number of points is considered the matrix inversion ( 12) can be used for the evaluation of the weighting coefficients.It is also noted that considering a limited number of points any polynomial basis (Tab. 1) within any distribution (Tab.2) can be taken into account for the evaluations of the weighting coefficients using the presented general scheme (12).On the contrary, when a single element is taken into account the weighting coefficients in exact form must be used to avoid ill-conditioning problems.The present work also illustrates some applications related to the Weak Formulation Finite Element Method (WFEM), where an integration procedure based on the weighting coefficients of the GDQ method is developed.The present technique has been described in the works [7,17].This approach has been termed the Generalized Integral Quadrature (GIQ) method.In general, the numerical integration of a function   f x over a domain   , a b can be defined as where k w are the weighting coefficients.The integral of   f x over a given domain can be generally approximated by a linear combination of all the functional values in the whole domain as

 
  The limits i x , j x of Eq. ( 18) can be changed.When  18) becomes a conventional integral (17).
The GIQ weighting coefficients can be computed as The weighting coefficients for the integral are evaluated by inverting the matrix   , which depends on the matrix of the weighting coefficients of the first order derivative [7,17] and they can be calculated by the following relations It is important to note that c is an arbitrary constant and can be set equal to  5), the following relation can be written It is noted that 1N

APPLICATIONS AND RESULTS
n the present section the static and dynamic behaviors of 1D and 2D structural components, summarized in Tab. .In all the static calculations the uniform applied load for both 1D and 2D problems is taken equal to 100 Pa q  . The convergence, stability and reliability of the SFEM varying the basis functions, collocation, number of grid points and number of elements are discussed in the following.The spacing of floating point numbers, due to the machine in use, is 52 16 . The static and dynamic convergence and stability of Rods, Euler-Bernoulli (EB) and Timoshenko (Tim) beams are shown.For the present studies the basis functions, the grid point distributions, the number of points in each domain and the number of elements for the domain decomposition can be selected.The EB beam has the additional parameter given by the  -points location.The authors defined a new way of inserting the  -points other than the classic one [9,12].This new approach consists in adding the  -points as functions of the first and third (last and two-before the last) points.The  -points location is indicated in brackets "( )".As far as the classic  -points technique is concerned, the extra points locations are directly reported as ( 510  ), whereas for the second approach the distance between the first and the third points is indicated.For instance, if (1/2) is used the  -points are the mid-points of the first and the third ones.In the following applications the GDQ method based on PDQ basis functions is considered.The reader can find the complete nomenclature of the grid point locations in Tab. 2. The one dimensional tests are divided into static and dynamic analysis.For the static benchmarks, the L2-norm of the absolute error of the displacement field is taken into account, since comparing a displacement of a single point of the domain can rise numerical instabilities.For the dynamic cases, the discrete spectra and the relative error on the first natural frequency will be shown in the following.Fig. 1a shows the L2-norm of the absolute error of the displacement for several structural elements composed of a single finite element varying the number of points N .The error given by the EB beam model is higher than the one exhibited by the Rod and Tim beam, because extra points have to be added in order to enforce the boundary conditions.The Che-Gau-Lob grid distribution is chosen within PDQ basis functions.For the EB beam model  -points are taken as (1/2), showing the best accuracy with respect to the classic value ( 5 10  ).The present method has been demonstrated to be stable even for 151 N  .Fig. 1b studies the accuracy of the domain decomposition by varying the number of elements with 7 N  in each element.The same basis functions and grid distribution of the previous case was used.Increasing the number of grid points or the elements, the same trends are obtained when Fig. 2a and 1b are compared.In the following the free vibration problems of the Rod, EB and Tim beams are reported in Fig. 2b, 3 ) is investigated.Secondly the same structural elements are divided into 100 e n  of low degree ( 7 N  ).The dynamic discrete spectra give a very close picture of the complete dynamic behavior given the total number of degrees of freedom.The percentage error with respect to the exact solution is plotted against the mode number.Globally in Fig. 2b, the accuracy is around 60% (for the Rod and the EB beam) and 50% (for Tim beam), when a single element is taken into account.However, it is uncommon to use a single element with many points in it.Thus, the SFEM approach is investigated in Fig. 3a where 100 e n  and 7 N  .The  -points for the EB beam are taken as (1/2) and the global accuracy results to be good until 30% of the modes, whereas the other Rod and Tim beam models reach the 60% and 40% respectively.These plots report some discrete jumps through the whole spectrum that were not investigated in the present paper.Hence, these aspects will be deeply exhibited in a following work.Finally the dynamic discrete spectra of the WFEM are represented in Fig. 3b for the Rod and the Tim beam.It can be noted that the curves start to detach from the horizontal line with respect to Fig. 3a where SFEM approach is presented.In fact, Rod and Tim beam are accurate until the 20% of the modes.Fig. 4a reports the convergence rate of the SFEM in the double logarithmic plot of the relative error of the first frequency versus the number of elements.
In the study the degree of the approximating polynomials is 6 ( 7 N  ) and the number of element varies from 1 to 100.PDQ basis functions and Leg-Gau grid distribution are taken into account.It is remarked that the Rod and the Tim beam have a rate of convergence higher than the EB beam due to the fact that fourth order differential equations need extra points for the imposition of the boundary conditions, hence 2 grid points (per boundary) are lost for the imposition of the boundary conditions.Furthermore, the maximum accuracy that can be reached by the Rod is around 15 10  and the Tim beam is around 12 10  , on the contrary the EB beam reaches 9 10  .For the sake of comparison, the interested reader can refer to the works [36,37] where the free vibrations of Rods and EB beams are investigated within the Isogeometric Analysis (IGA).For instance Fig. 10-12 of reference [37] show the convergence rates of a Rod with different polynomial orders.As expected, the convergence rate changes from 1/4 to 1/8 increasing the polynomial order.Nevertheless for all three cases the author shows only four points for the Rod convergence.In fact Fig. 10 and  e n  .The author is assuming that the minimum convergence is reached, hence no more elements are needed.Nevertheless, the machine error noise is never shown throughout the paper.Ultimately, it is noted that the results proposed by [37] are in perfect agreement with the ones presented by the authors in Fig. 4a.The only exception is given by EB beams with  -points, because a lower trend 1/6 occurs together with a lower accuracy 9 10  .The interested reader can find more details about the so-called round-off error in the book by Boyd [15].In conclusion, Fig. 4b, 5a, 5b exhibit a summary of the Rod structure behavior using four different numerical approaches: SFEM, WFEM, SEM (Spectral Element Method) and FEM.The first solves the strong form and uses 1  C compatibility conditions, the second solves the weak form and uses 1  C compatibility conditions, the third and the fourth solve the weak form with 0 C continuity, but the former with a general polynomial degree (Lagrange polynomials) and the latter with well-known linear and quadratic functions.Fig. 4b shows that increasing the polynomial order inside each element the error increases.In fact the SFEM error with 7 N  is higher than the FEM one with 2 N  .Nevertheless, the maximum error is always below the machine working precision (black dashed line).A comparison in terms of dynamic discrete spectra in shown in Fig. 5a.It can be noticed that SFEM offers the best accuracy in terms of percentage of accurate modes and the maximum errors are always more limited than all the other techniques.Fig. 5b shows several convergence rates with respect to the first natural frequency.Obviously, increasing the polynomial degree the convergence rate passes from 1/2 to 1/10.Linear FEM ( 2 N  ) has the same trend (1/2) of the ones obtained by SFEM and WFEM with 3 N  and 1 C continuity.The strong formulation requires a higher derivation order with respect to a variational or weak formulation.Considering the second order Rod problem, for a given basis function of order n the strong form derives the approximate solution two times, whereas the weak form only one.Therefore, the strong form solution becomes of order 2 n  and the weak form is 1 n  .Moreover, since the 1  C continuity is considered in the strong form the boundary points are excluded from the computation due to the kinematic condensation.

Axially-loaded beam, Bar or Rod problem Static analysis
Dynamic analysis

Static analysis
Dynamic analysis Boundary conditions Simply Supported-Simply Supported (SS) Simply Supported-Simply Supported (SS)

Reissner-Mindlin (RM) thick plate problem Static analysis
Dynamic analysis Boundary conditions: Completely Simply-Supported (SSSS) ,0 0, ,0 0, 0, 0, 0 , 0, , 0, , , 0    In conclusion, in order to have the same accuracy trend SFEM or WFEM with 1 C compatibility conditions should have a higher number of grid points with respect to SFEM or SEM with 0 C compatibility conditions.Remark that the first point of the curves in Fig. 5b refer to a single element mesh.The comparison between SFEM and the other weak methods is meaningless for this case because, excluding the boundary conditions, the SFEM has an 1 n  order with respect to the others due to the strong formulation (second order derivatives with respect to first order ones).It should be also noted that not only SFEM and WFEM but also SEM has the increasing error machine effect when the maximum accuracy is reached (around 12 10  ).Fig. 6-8 show some results related to the static and dynamic stability and accuracy of Membranes, Kirchhoff-Love (KL) and Reissner-Mindlin (RM) plates.Fig. 6a is analogous to Fig. 1a where the number of grid points N varies from 5 to 31 using PDQ basis functions and Che-Gau-Lob grid distribution.The KL plate model is considered with  -points ( 5 10  ).Note that when 15 N  the error is stable increasing the number of points (round-off plateau).Fig. 6b collects the static convergence behavior of a single domain decomposed into several regular elements (no mapping is involved).The number of grid points is fixed 7 N  and the number of elements is variable 1, 4, 9, 16, 25, 36, 49, 64.For the present cases the PDQ basis functions and Leg-Gau grid distribution are considered.As expected the error decreases increasing the number of elements.The dynamic convergence of several structural components is shown in Fig. 7 and 8. Fig. 7a shows the dynamic discrete spectra of the Membrane, the KL and the RM plates using PDQ basis functions and Che-Gau-Lob grid distribution with 31 N  using a single element.The Membrane and KL plate calculate accurate modes until the 25% of total modes, whereas the RM plate gets 10% of accurate modes due to the comparison with a semi-analytical solution.An analogous plot using PDQ basis functions and Leg-Gau grid distribution with 7 N  and 64 e n  is shown in Fig. 7b for the Membrane and the RM plate.Due to the element decomposition the accuracy decreases both for the Membrane and the RM plate.However, the error is limited within the 5% for the 40% and 50% of the modes of the Membrane and the RM plate, respectively.In conclusion a comparison in terms of dynamic discrete spectra for structured and distorted meshes is presented in Fig. 8 for the Membrane and RM plate, respectively.Different number of grid points is considered, in order to get the same number of degrees of freedom changing the number of elements.As expected the modes are more accurate when the mesh is structured, whereas a lower accuracy is observed when irregular meshes are considered.

CONCLUSIONS AND REMARKS
he present paper references the DQM and other methods that, according to the authors' knowledge, represent similar techniques.The paper aims to present a general view on strong formulation methods.The main novelty of this manuscript is given by the presentation in several forms of the stability and accuracy of the SFEM technique when applied to simple 1D and 2D models and compared to exact solutions (or semi-analytical ones as far as RM plate is concerned).This work should help researchers, of the computational mechanics community, to understand the advantages and disadvantages of a strong formulation approach and, in particular, to the ones who are keen on weak (or T variational) formulation based approaches.Although it seems that DQM has been widely developed for several engineering problems, some aspects are still in a developing stage (e.g.boundary conditions in a domain decomposition approach using 1 C conditions).In a future paper the authors want to expose in detail the boundary conditions implementation from the mathematical point of view, giving all the formulae in discrete form.Moreover, the numerical results will be compared to classical FEM and SEM solutions.
a) Dynamic discrete spectra for various structural components with 7 N  and 100 e n  for a CC Rod, a SS EB beam with -points (1/2) and a CC Tim beam using PDQ basis functions and Leg-Gau grid distribution.b) Dynamic discrete spectra for several structural components using Weak Formulation (WFEM) with 7 N  and 100 e n  for a CC Rod and a CC Tim beam using PDQ basis functions and Leg-Gau-Lob grid distribution.
a) Relative error of the first frequency for various structural components with 7 N  for a CC Rod, a SS EB beam with  - points (1/2) and a CC Tim beam using PDQ basis functions and Leg-Gau grid distribution.b) Effect of the grid point number inside each element for the static analysis of a CC Rod using SFEM with PDQ basis functions and Leg-Gau grid distribution, WFEM and SEM with PDQ basis functions and Leg-Gau-Lob grid distribution and FEM with linear and quadratic shape functions.
Rod using SFEM with PDQ basis functions and Leg-Gau grid distribution, WFEM and SEM with PDQ basis functions and Leg-Gau-Lob grid distribution and FEM with linear and quadratic shape functions.b) Relative error of the first frequency varying the number of elements e n for a CC Rod using SFEM with PDQ basis functions and Leg-Gau grid distribution, WFEM and SEM with PDQ basis functions and Leg-Gau-Lob grid distribution and FEM with linear and quadratic shape functions.
a) Static analysis for single element structural components varying the number of points N for a Membrane, a KL plate with  -points ( 5 10  ) and a RM plate using PDQ basis functions and Che-Gau-Lob grid distribution.b) Static analysis for structural components varying the number of elements e n with 7 N  for a Membrane and a RM plate using PDQ basis functions and Leg-Gau grid distribution.a) b) Figure 7: a) Dynamic discrete spectra for single element structural components with 31 N  for a Membrane, a KL plate with  - points ( 5 10  ) and a RM plate using PDQ basis functions and Che-Gau-Lob grid distribution.b) Dynamic discrete spectra for various structural components with 7 N  and 64 e n  for a Membrane and a RM plate using PDQ basis functions and Leg-Gau grid distribution.
obtain stable and accurate results.It must be pointed out that in the applications, the coordinate transformation from an interval defined, especially when different basis functions are used for the definitions of the weighting coefficients.Thus, the transformation of the GIQ weighting coefficients 1Nk w allows to switch from the interval   , c d to a generic one   , a b .Recalling Eq. (

kw
are the weighting coefficients in the shifted interval   , c d and the 1N k w are the ones in the physical interval   , a b .
, 4a.In the computations d N represents the total number of degrees of freedom I and m n indicates the mode number which is plotted in the discrete spectra.It is remarked that plate.Fig.2band 3a illustrate the dynamic discrete spectra by following the same analyses of Fig.1a and 1b.Firstly a single element of high degree ( 151 N  a) Dynamic discrete spectra using different meshes (structured and distorted) of a Membrane using PDQ basis functions and Leg-Gau grid distribution.b) Dynamic discrete spectra using different meshes (structured and distorted) of a RM plate using PDQ basis functions and Leg-Gau grid distribution.

Table 2 :
List of several grid distributions used in structural mechanics applications.

Table 3 :
List of one dimensional static and dynamic exact solutions.

Table 4 :
List of two dimensional static and dynamic exact solutions.