Carlos A.Felippa

Department of Aerospace Engineering Sciences and Aerospace Mechanics Research Center,University of Colorado at Boulder,Boulder,CO 80309-0429,USA

ABSTRACT This paper considers Lagrangian finite elements for structural dynamics constructed with cubic displacement shape functions.The method of templates is used to investigate the construction of accurate mass-stiffness pairs.This method introduces free parameters that can be adjusted to customize elements according to accuracy and rank-sufficiency criteria.One-and two-dimensional Lagrangian cubic elements with only translational degrees of freedom(DOF)carry two additional nodes on each side,herein called side nodes or SN.Although usually placed at the third-points,the SN location may be adjusted within geometric limits.The adjustment effect is studied in detail using symbolic computations for a bar element.The best SN location is taken to be that producing accurate approximation to the lowest natural frequencies of the continuum model.Optimality is investigated through Fourier analysis of the propagation of plane waves over a regular infinite lattice of bar elements.Focus is placed on the acoustic branch of the frequency-vs.-wavenumber dispersion diagram.It is found that dispersion results using the fully integrated consistent mass matrix (CMM) are independent of the SN location whereas its lowfrequency accuracy order is O(κ8),where κ is the dimensionless wave number.For the diagonally lumped mass matrix(DLMM)constructed through the HRZ scheme,two optimal SN locations are identified,both away from third-points and of accuracy order O(κ8).That with the smallest error coefficient corresponds to the Lobatto 4-point integration rule.A special linear combination of CMM and DLMM with nodes at the Lobatto points yields an accuracy of O(κ10)without any increase in the computational effort over CMM.The effect of reduced integration(RI)on both mass and stiffness matrices is also studied.It is shown that singular mass matrices can be constructed with 2-and 3-point RI rules that display the same optimal accuracy of the exactly integrated case,at the cost of introducing spurious modes.The optimal SN location in two-dimensional,bicubic,isoparametric plane stress quadrilateral elements is briefly investigated by numerical experiments.The frequency accuracy of flexural modes is found to be fairly insensitive to that position,whereas for bar-like modes it agrees with the one-dimensional results.

KEYWORDS Structural dynamics;Lagrangian elements;finite elements;cubic shape functions;bar;plane stress;mass;stiffness;vibration;wave propagation;Fourier analysis;dispersion;templates

1 Introduction

The accompanying contribution in this issue [1] narrates the author’s effort in understanding and draining the “variational swamp” that emerged during 1965-1985.In those two decades the Finite Element Method (FEM) was liberated from the Rayleigh-Ritz connection,but fragmented developments had resulted in a dizzying array of variational tools.The effort resulted in the unification published in [2],which made a wide range of classical principles in structural mechanics (e.g.,Total Potential Energy,Hellinger-Reissner,Veubeke-Hu-Washizu,etc.) particular cases of a single parametrized functional.Setting parameters to specific values produced well known instances,plus an infinity of new ones.

A practical side effect:element derivations could be carried out by splitting the formulation and using two or more functionals.The split could be geometric:element boundary plus interior,as in the by-now large class of hybrids and incompatible models [3–5].Or additive:basic plus higher order response,as in Bergan’s Free Formulation [6].Or splitting done in the matrix computation process,as in elements produced by selective integration [7].

Investigation of common features emerging in those combinations led to the method of templates,which provides tools for the present paper.Templates are algebraic forms for the mass,material and geometric stiffness matrices of an individual element.Distinguishing feature:templates contain free parameters.These may be either simply scaling coefficients,or carry immediate physical meaning.In this paper two parameters are used:

· A geometric parameterγthat adjusts the side node location.

· A scaling parameterμthat defines the combination of consistent and diagonally-lumped mass matrices.

A short history of templates,adapted from a recent survey paper [8] is given in Appendix A.The present paper is a specialized continuation of a 2015 survey paper [9] that studies massstiffness templates for structural dynamics in a more general context,following up on earlier work [10,11].

2 Lagrangian Cubic Elements

Lagrangian structural elements constructed with cubic displacement shape functions possess only translational degrees of freedom (DOF).They carry additional nodes in addition to corners,Fig.1.These will be denoted as “side nodes” (SN) in the sequel,although for one-dimensional(1D) elements such as that shown in Fig.1a,they are also internal nodes.SN are usually placed at the thirdpoints,but their location may be symmetrically adjusted,within geometric limits,using a template parameter.For multidimensional Lagrangian elements such as those pictured in Figs.1c–1e,the SN location determines placement of internal and face nodes.

The SN adjustment effect is studied here in detail for the simplest case:the four-node,elastic,prismatic,cubic bar element of Fig.1a.TheMathematicacomputer algebra system is used for symbolic computations.Key results are then verified with numerical calculations.The stiffness and mass matrices are constructed through isoparametric cubic shape functions parametrized as per SN location.

The approximation properties are determined by Fourier analysis of plane waves propagating over an infinite lattice of identical elements.This provides dispersion diagrams that connect natural frequencies to dimensionless wavenumbers.Those diagrams include one acoustical branch and two optical branches.Of these the acoustic branch is of paramount interest as regards approximation of continuum natural frequencies.On the other hand,optical branches are physically spurious and a byproduct of the discretization process (They affect primarily the so-called numerical pollution).

The result of changing the SN location on acoustic branch accuracy is further studied for several mass-stiffness combinations differing on Gauss integration rules.The effect of mass lumping,as well as taking a linear combination of the consistent and diagonaly-lumped matrices,is considered.

The paper concludes with a numerical investigation of the performance of 2D plane stress Lagrangian elements in a beam vibration problem.Their inherent complexity makes symbolic analysis unfeasable on present computer algebra systems,a limitation explained in Appendix A.

3 The Four-Node Bar(Bar4)Element

Consider a four-node prismatic bar element of lengthl,constant mass densityρ,elastic modulusE,and cross section areaA,which can only move along its longitudinal axis,Fig.2.This model is acronymed Bar4 for brevity’s sake in the sequel.Nodes are locally numbered as shown.The axial coordinate isxwith origin at element center.The axial displacement isu=u(x).The isoparametric (isoP) natural coordinateξvaries fromξ=-1 toξ=1 from end Node 1 to end Node 2,respectively.The position of the side nodes (SN) is specified by the template parameterγshown in the figure,so that their isoP coordinates areξ=±γ.The SN arrangement is constrained to be symmetric respect to the centerξ=0.If the sign ofγis reversed,the nodes just flip,whence the allowable position range is taken to beγ∈(0,1),which excludesγ=0 andγ=1.

3.1 Bar4 Shape Functions

The isoparametric formulation based on theγ-parametrized shape functions is used.Forγ=1/3,(1) are the well known cubic shape functions for equidistant nodes [12].Ifγ=≈0.4472135955 they become the 4-point Lobatto shape functions.The shape function matrix is

The overall element isoparametric definition can be concisely written [13]

Here both x and u are column 4-vectors,with x=[-1 1-γ γ]T.Expanding:x=Nex=ℓ(-N1e+N2e-γ N3e+γ N4e)=ξ ℓ,whence the Jacobian determinant is

This is constant over the element and independent ofγ,which greatly simplifies integration.

3.2 Bar4 Strain-Displacement Equations

The natural derivatives of (1) are

The axial strain is given by

Here Beis the 1×4 strain-displacement matrix with entriesBei=J-1dNei/dξfori=1,...4.

4 Bar4 Mass and Stiffness Matrices

The element consistent mass matrix (CMM),parametrized in terms ofγ,is defined as

Here Neis the 1×4 shape function matrix (3),which depends onγ.For a prismatic homogeneous element,ρandAmay be taken out of the integral,as well asJ=l/2.The expression (7) can be evaluated using Gauss integration rules withpMpoints,in which casepMis appended as second argument of MeC.For constantρ,AandJ:

in whichwkandξkare weights and abcissas,respectively,of thepM-point Gauss rule.Since the entries of Neare cubic polynomials inξ,exactness is attained ifpM≥4.

The element stiffness matrix,parametrized in terms ofγ,is given by

in which Beis theγ-dependent strain-displacement matrix (6).For a prismatic homogeneous element,EandAmay be taken out of the integral,as well asJ=l/2.The expression (9) can be evaluated by a Gauss rule withpKpoints,in which casepKis appended as second argument of Ke:

The entries of Beare quadratic polynomials inξ,whence exactness is attained ifpK≥3.

4.1 Fully Integrated Mass Matrices

The fully integrated (FI) consistent mass matrix (CMM) is obtained by evaluating the integral(7) either analytically,or through the Gauss rule (8) withpM=4.The resulting matrix is rank sufficient.This exact CMM is denoted by

HereMCijandχCijdenote physical and dimensionless entries,respectively.The entry dependence onγis omitted from (11) to reduce clutter.Their expressions may be found in Table 1.

The associated diagonally-lumped mass matrix (DLMM) is denoted by

The diagonal entries are obtained via the HRZ scheme [14] as the row sumsχL1=χC11+χC12+χC13+χC14,χL2=χC21+χC22+χC23+χC24,χL3=χC31+χC32+χC33+χC34,andχL4=χC41+χC42+χC43+χC44.The entry dependence onγis omitted from (12) to reduce clutter.Their expressions are listed in Table 1.

Table 1:Bar4 consistent and diagonally-lumped mass matrices

Notes:FI:Fully integrated,RI:Reduced integration,GRp:p-point Gauss rule.Physical mass entries:MCij=MCji=(ρAℓ)χCij and MLii=(ρAℓ)χLi.

4.2 Fully Integrated Stiffness Matrix

The fully integrated (FI) stiffness matrix is obtained by evaluating the integral (9) either analytically,or through the Gauss rule (10) withpK=3.It will be denoted as

HereKijandψijdenote physical and dimensionless entries,respectively.To reduce clutter their dependence onγis omitted.Their expressions may be found in Table 2.This matrix has the full rank of 3.Its only zero-energy mode is the rigid body translation alongx.

Table 2:Bar4 stiffness matrices

4.3 Reduced Integration Mass and Stiffness Matrices

There has been recent work on nonstandard methods for constructing mass matrices,e.g.,[15–17].In particular,singular element mass matrices have attracted attention for direct time integration (DTI) computations,as a scheme to reduce numerical pollution caused by spurious optical branches [18,19] This is important in problems exhibiting rapid transients,as well as in multibody dynamics.Three procedures for generating singular mass matrices are described in the survey paper [9]:

· Spectrally truncated templates

· Projection onto a null subspace

· Reduced integration (RI)

Of these,only RI will be considered here because it links up smoothly with the fully integrated(FI) case.RI mass matrices of interest are obtained using Gauss rules withpM=2 orpM=3 points.The resulting CMM are denoted byMeC(γ,pM)and likewise for the associated DLMM constructed through the HRZ scheme.

Terminology:an eigenvector corresponding to a zero eigenvalue is called anull eigenvector.If the zero eigenvalue is multiple,the set of all null eigenvectors forms anull eigenspace.Spurious modesare null eigenvectors devoid of physical meaning.

The CMM withpM=3 has rank deficiency 1,which pertains to a spurious cubic mode that vanishes at all three Gauss points.The CMM withpM=2 has rank deficiency 2,which pertains to quadratic and cubic modes that vanish at both Gauss points.Their entries may be found in Table 1.(ThepM=1 version is of no interest.) Associated diagonally lumped matrices do not change.

A RI stiffness matrix is produced by a Gauss rule with 1 ≤pK≤2 points,denoted by K(γ,pK).Of the two possibilities,only Ke(γ,2) computed through the 2-point rule,is of some interest.This matrix has a rank deficiency of 2.The null subspace includes the physical rigid body mode,(uniform longitudinal motion) as well as a cubic mode whose strains vanish at the 2 Gauss points.The matrix is expressed in the form (9).Its entries may be found in Table 1.

4.4 Tradeoffs of Reduced Integration

The alleged benefits of RI in reducing DTI pollution may be negated by the presence of spurious modes.These can be benign or harmful.Three cases are considered below.

The 3-point integrated CMM of an individual element,denoted byMeC(γ,3),has an antisymmetric cubic spurious mode (the null eigenvector ofMeC).This mode,denoted byuA(ξ),vanishes at the three Gauss pointsξ=0 andξ±(SubscriptAstands for antisymmetric).Normalizing so thatuA=±1 atξ=±1 we have

This is independent ofγ.The mode is pictured in Fig.3a.If this mass matrix is used to model a free-free bar,the mode propagates over elements in the repeating pattern illustrated in Fig.3b.If the motion of any corner node (or of a side node ifγ/=) is set to zero,the mode disappears.Consequently it can be characterized asbenign,since it stays unique no matter how many elements are assembled,and can be effectively suppressed by just one homogeneous boundary condition.

The 2-point RI mass matrixMeC(γ,2)of an individual element is rank deficient by two.The spurious modes are a linear combination of quadratic and cubic functions that vanish at the two Gauss pointsξ=±1/It is convenient to select the following normalized null eigenvectors as basis of the null eigenspace:in which subscriptsSandAstand for symmetric and antisymmetric,respectively.These are pictured in Fig.3c.Their shapes are independent ofγ.To remove the null eigenspace,both ends must be fixed.For multiple elements the modes propagate as illustrated in Fig.3d.This plot shows three null eigenvectors for a fixed-fixed bar discretized with 4 elements;note that they vanish at the 8 Gauss points.If this problem is discretized withNeelements,the assembled M hasNe-1 zero eigenvalues.Thus this particular CMM may be characterized asdangerousfor anyγ,since spurious modes are hard to control with usual boundary conditions.

The 2-point RI stiffness matrix Ke(γ,2) has rank deficiency of 2.The null eigenspace is spanned by

HereuR(ξ) is the rigid body translational mode,which is physical.On the other hand the antisymmetric modeuA(ξ) is spurious:its derivative is zero at the Gauss points ±so associated strains vanish there.The spurious mode of (16) has been normalized so that its value is±1 at those points,as shown in Fig.3e.Its shape is independent ofγ.SinceuS(ξ) vanishes at the end nodes,it cannot be controlled by homogeneous boundary conditions there:it will “sprout”over each element of a mesh assembly.This is illustrated in Fig.3f for a fixed-fixed bar discretized with 4 elements:the assembled stiffness has three zero eigenvalues with the eigenvectors shown.If this problem is discretized withNeelements,the assembled K hasNe-1 zero eigenvalues.It follows that this particular stiffness may be characterized asdangerousfor anyγ.

The preceding analysis applies to the mass and stiffness matrix in isolation.When combined for dynamic analysis,more complex spurious mode effects may occur.

4.5 A New Result:Gamma Invariance

The fully integrated CMM and its RI cousins display a remarkable property:the computed natural frequenciesdo not depend on SN location.More precisely,they are independent ofγif 0＜γ＜1.This property will be generically calledγ-invariance.It can be illustrated through the frequency characteristic function (FCF) of an individual element:

in whichωdenotes circular frequency.Expanded FCF expressions are given in Table 3 forpM=4,3,2 andpK=3,2.To reduce clutterE,A,ρandlare set to 1.It can be observed that the numerators are polynomials inω2that do not depend onγ,while their denominators do not vanish ifγ∈(0,1).The natural frequenciesωi,which are roots ofF=0,are tabulated in the last column of that Table.(Note that if the reduced stiffness withpK=2 is used,one extra root becomes zero,whereas if a reduced mass matrix is used,some roots move to ∞because theω-polynomial order drops).As regards vibration modes,which are the eigenvectorsφiof the frequency eigensystem Keφi=ω2iMeCφi,their expressions as functions ofξover the element are also independent ofγif their values at the SN are adjusted accordingly.

The property is somewhat surprising since the eigenvalues of MeCand Ke in isolationdo depend onγ.Symbolic computations indicate thatγ-invariance remains valid for arbitrary assemblies of Bar4 elements that use CMM.But it does not hold for other mass matrices such as DLMM.

Conclusion:If a CMM is used,whether FI or RI,the choice ofγmakes no difference in frequency accuracy.Is there an optimal value in another sense? From a conditioning standpoint one would like to maximize the FCF denominators.As shown in Table 3,these have the formC γ2(1-γ2)4,in whichCis a numerical coefficient independent ofγ.The maximizing value in the range (0,1) isγ=which is an interior Lobatto point of the 4-point Lobatto integration rule [20].This value recurs in several optimization scenarios studied later.

Table 3:Frequency characteristic function (FCF) of a CMM-endowed Bar4 element

4.6 Noteworthy Instances

Some mass-stiffness instances worth noting are collected here.If the interior nodes are placed at the thirdpoints,γ=1/3 and the FI mass matrices are

whereas the FI stiffness matrix becomes

For this particular configuration,the DLMM displays the well known three-eighths integration formula,which is also the 4-point Newton-Cotes quadrature rule [20].

If the interior nodes are placed at the Lobatto points,γ=and the FI mass matrices become

whereas the FI stiffness matrix becomes

The diagonal entries and node locations of the DLMM in (20) correspond to the abcissas and weights of the 4-point Lobatto integration rule [20].By analogy,the foregoing matrices will be denoted as theLobatto DLMMandLobatto stiffness matrix,respectively.

Forγ=which (as shown later) gives a near-optimal DLMM,the FI mass matrices are

in whichχ13=andχ14=

5 Plane Wave Propagation over Bar4 Lattice

We consider the propagation of plane waves over an infinite lattice of identical Bar4 elements obtained by space discretization of a continuum bar model.Fig.4 is a schematics of the continuum-to-lattice semi-discretization process.The formulation concludes with the extraction of a two-element patch about a typical end noden.The necessary terminology is collected in Table 4.

Table 4:Nomenclature for harmonic plane wave propagation over Bar4 lattice

5.1 Plane Wave Equations

Plane wave propagation over a regular spring-mass lattice is governed by the semidiscrete,linear equation of motion (EOM):

Here u is the vector of axial nodal displacements while M and K are infinite,tridiagonal,Toeplitz mass and stiffness matrices.This EOM can be solved by Fourier methods.Fig.4b displays two characteristic lengths:λandl.The element length-to-wavelength ratio is calledϒ=l/λ.The floor function of its inverse:Neλ=λ/ℓis the number of elements per wavelengths.Those ratios characterize the fineness of the discretization,as pictured in Fig.4b.

Within constraints noted later the lattice can propagate real,travelling,harmonic plane waves of wavelengthλand group velocityc,as depicted in Figs.4b,4c.The wavenumber isk=2π/λand the circular frequencyω=2π/T=2π c/λ=k c.The range of physical wavelengths that the lattice may transport is illustrated in Fig.5.

To study plane wave solutions of (23) it is sufficient to extract a repeating two-element patch,a process schematized in Fig.4d.The seven patch nodes are renumbered as shown in Fig.6,and thexaxis origin relocated to the interface noden.The natural patch coordinate=x/ℓis defined as shown in the figure.A harmonic plane wave of amplitudeBis described by the function

In the second form the dimensionless wavenumberκ,circular frequencyω,and timeτare defined asκ=k l=2πl/λ,ω=ω l/c0,andτ=tc0/l,respectively.Herec0=is the elastic continuum bar group wave velocity,which for that model is the same as the phase velocity.(In physical acousticsc0is thesound speedof the material.) In both formsBhas dimension of displacement.

Using the stiffness and mass matrix notation of (11)–(13) the patch equations can be written

in which

Here theχijandψijare the dimensionless coefficients introduced in (11) and (13).(The extra subscriptCis suppressed since this form applies to all admissible mass and stiffness matrices.)From the patch EOM (25) take the middle 3 equations (Rows 3,4 and 5 of MPand KP),which repeat in the infinite lattice:

5.2 Fourier Analysis

The axial wave motion is represented using three basis functions:uB1(ξ,τ),uB2(ξ,τ) anduB3(ξ,τ) with the same structure as the last term in (24)

CoefficientsBiare linked toBby

in whichgiare functions of patch position,chosen so they are linearly independent and none is identically zero.For future use theBiare grouped into the column 3-vector B=B1B2BT3,and thegiinto the 3×3 diagonal matrix G=diag[g1,g2,g3],whence B=GB.

For convenient implementation we pick thegias sketched in Fig.6,namelyB1vanishes at all nodes but 1,4 and 7,B2vanishes at all nodes but 2 and 5,andB3vanishes at all nodes but 3 and 6.There is an infinite number of suchgifunctions.However,their actual form is irrelevant because the eigenvalues of the characteristic system (36) derived below do not depend upon theAccordingly,the nodal patch displacement vector is set as follows

Substituting (30) into (27) yields the wave propagation condition

If C B=for B/=0,C must be singular.This furnishes thecharacteristic equationdet(C)=0,from which factors depending onτonly may be dropped.The dimensionless frequencyΩappears only in even powers;more preciselyΩ2,Ω4andΩ6.It is convenient to define=Ω2as new variable.Then det(C)=0 is a (complicated) cubic equation inthat provides three solutions denoted as=Ω21,=Ω22,and=Ω23.One of those vanishes atκ=0 and is renamedΩafor acoustic branch.The remaining two roots correspond to optical branches and are relabeledΩ2o1andΩ2o2,which are ordered soΩ2o1≤Ω2o2atκ=0 (these are the so called cutoff frequencies).When expressed as functions ofκ:these represent the equations of the acoustic and optical branches.The combined plot of these branches asΩ(κ) curves is called thedimensionless dispersion diagramor DDD.This is done for the Brillouin zoneκ∈[0,2π].

6 Low Frequency Performance

The dimensionless series expansions of the acoustic branch aboutκ=0 is written

whereas those of the optical branches are

The performance of a mass-stiffness (MS) pair in structural dynamics is controlled by the accuracy of the acoustic branch for low frequencies,i.e.,small wavenumbers.Since the acoustic branch for the continuum isΩa=κ,orΩ2a=κ2,the accuracy is controlled by the first nonzero term afterκ2in theΩ2aseries (33).The will be called theprincipal error term,or PET,in the sequel.For example ifA4=0 andA6=-24,the PET isA6κ6/6!=-κ6/30.

6.1 Performance of CMM and DLMM

Table 5 shows the coefficients of (33) for fully integrated (pM=4) CMM mass matrices,as well as two RI instances obtained withpM=2 and 3 Gauss points.The stiffness matrix is either the exactly integrated one withpK=3 or the RI instance withpK=2.As can be observed the expansion does not depend onγ,confirming theγ-invariance analysis of Section 4.5.Thus the accuracy order displayed in the last column cannot be improved by moving nodes.

Table 5:Acoustic branch series for CMM

Table 6 shows the coefficients of (33) for DLMM mass matrices.These matrices are identical forpM=2,3,4,as shown in Table 1.Consequently only two rows are needed in Table 6 for ML(γ,pM) paired with the FI and RI stiffness matrices.The series coefficients now depend onγ—more precisely its even powers.The first row show that the pairing of Me(γ,pM) with Ke(γ,3)has generally accuracy orderO(κ6).This may be elevated toO(κ8) by solving-2A6=1-15γ2+50γ4=0.This yields two solutions:γ2B1=1/5 andγ2B2=1/10,in which subscriptBstands for“best.” Taking the positive square roots delivers the DLMM-optimal node locations

Table 6:Acoustic branch series For DLMM

The corresponding PETs are-(2/15)κ8/8! and-109/240κ8/8!,respectively.The first coefficient is about 3 times smaller (actually 109/32 ≈3.406 smaller).Thus the internal node locationγB1=which corresponds to the 4-point Lobatto rule,is optimal for DLMM together with Ke(γ,3).Note that ifγ=1/3,the PET is+(4/81)κ6/6!,which forκ=1 is about 151 times bigger than optimal.As regards the RI stiffness matrix Ke(γ,2),the error orderO(κ4) may not be raised by solvingA4=4γ2/(1-γ2)=0 becauseγcannot vanish.

Dimensionless dispersion diagrams for various instances of CMM and DLMM are shown in Figs.7 and 8,respectively.

6.2 Performance of CMM-DLMM Combinations

The truncation error-(2/15)κ8/8! for the optimal DLMM withγ=andpM=2,3,4 is smaller than that of the exactly integrated CMM,which is (2/5)κ8/8!,and of opposite sign.Consequently a linear combination of those matrices,henceforth denoted by LCMM,looks promising as far as elevating the error order toO(κ10).The combination weights areμand 1-μ,whereμis another template parameter.This follows the notation of [9].

Table 7 lists coefficients (computed byMathematica) of the acoustic branch series (33) for the two-parameter mass matrix combinations

Table 7:Acoustic branch series for DLMM

paired with Ke(γ,pK) forpK=3,2.The catalogued coefficients extend up toA10forpK=3,and toA6forpK=2,in accordance to the expected principal error terms (PET).

If the stiffness matrix is fully integrated (pK=3),a PET ofO(κ10) can be obtained by settingA6=A8=0 and solving forγ2andμ.This produces a system of two equations,cubic inγ2and linear inμ.Of the three solutions,one gives negativeγ2,and another one yields an indefinite LCMM.The remaining one is

On taking the positive square root,we get the optimal locationsγ=forpM=3,4 andγ=forpM=2.The latter:γ≈0.338062 is very close to 1/3.All three MS pairs have exactly the same PET:-(2/15)κ10/10!,as may be verified.The optimal LCMM associated with the 3-and 4-point Gauss rules actually coalesce:

The optimal LCMM associated with the 2-point rule is different

As regards the RI stiffness withpK=2,order elevation toO(κ8) requiresA4=A6,a system that lacks solutions.Thus onlyO(κ6) may be attempted by settingA4=0,which leads to a quadratic equation inγ2.The solutiongives imaginaryγforμ∈(0,1).

7 Reducing Numerical Pollution

In structural dynamics,optical branches have no physical meaning.They are a byproduct of the superlinear spatial discretization.Their presence is harmless in vibration analysis or in transient response computations carried out via mode superposition,as they are automatically filtered by the eigensolver.In direct time integration (DTI),however,thay can be responsible for“numerical pollution” that can render the computed solution worthless.In fact,this is the main reason forpreferring low-order elementsin that case,since those do not carry along that harmful extra baggage.

If one is determined to use elements such as Bar4 in DTI,pollution may be reduced by moving or flattening optical branches so that acoustoptical gaps are widened,because noise that falls into those gaps does not propagate.Several techniques for doing that are discussed in detail in [9] for the 3-node bar.Here we simply consider the effect ofreduced integrationon the position of optical branches.For a quick assessment we look at the so-calledcutoff frequencies,which are the values taken by the optical branches at zero wavelength:κ=0.Those are denoted byΩc1=Ωo1|κ→0andΩc2=Ωo2|κ→0,ordered so thatΩc1≤Ωc2.

Table 8 gives cutoff frequency expressions for parametrized LCMM-stiffness pairs,including fully integrated (FI) as well as several reduced integration (RI) versions.Table 9 lists numerical values for selected mass-stiffness pair instances,again showing the effect of the integration rule.

An effective way to cut pollution is to move an optical branch to infinity.As Table 8 shows,this can be done by using a RI stiffness paired with any mass matrix.To get rid of the other branch one may setγandμso that its cutoff frequency denominator vanishes.For exampleμ=1 andγ=(which places the SN at the 2-point Gauss locations) gets rid of both optical branches.Inevitably this reduces low frequency accuracy,which is onlyO(κ4) as shown in Section 6.2.In that case (μ=1 andγ=) the acoustic branch PET is down toκ4/12.This is worse that the best-fit instance of the much simpler 2-node bar element,which does not possess an optical branch.

Table 8:Cutoff frequencies of LCMM-stiffness pairs for several integration rules

Table 9:Cutoff frequencies of mass-stiffness instances for several integration rules

Notes:Cutoff frequencies are values of optical branch frequenciesΩo1andΩo2atκ=0.Ifμ=0(CMM)cutoff frequencies do not depend onγ.Ifμ/=0 there areγvalues(e.g.,γ=1/forμ=1)at which a finite cutoff frequency“takes off”to ∞.

To maintain low-frequency accuracy while alleviating pollution it would be necessary to consider more general mass templates,as done in [9] for the 3-node bar element.RI is insufficiently powerful to accomplish that goal.Such study has not been carried out.

8 Bar4 Test:Vibrations of a Fixed-Free Bar Member

Natural frequency predictions using several Bar4 mass instances are compared for predicting the first 3 natural frequencies of longitudinal vibrations of the fixed-free elastic bar member pictured in Fig.9.The member is prismatic,with constantE=1,A=1,andρ=1.The total member length is taken asL=π/2 for convenience.With those numerical properties the continuum eigenfrequencies are

The member is divided intoNeidentical elements,withNe=1,2,...8.Fig.9a illustrates the caseNe=4.Five template instances,identified in Table 10,are tested.Both mass and stiffness are fully integrated.Numerical results obtained for the first three frequencies are collected there.The higher convergence of BLCD is obvious.For example,2 elements giveω2correct to 4 digits while both CMM and DLMM give only 2.CMM,SDMM and LDMM overestimate frequencies whereas TDMM and BLCD underestimates them.

Table 10:Bar4 results for vibrations of a fixed-free bar member

The results of Table 10 are graphically reformatted in Fig.10,as log-log plots that display accuracyvs.element count.The horizontal axis shows number of elementsNein log2scale.The vertical axis displays correct digits of computed frequency,computed as

HereΔωiis the frequency error of computed valuesωFEMiwith respect to the continuum frequencies (40).Note that accuracy-log plots such as those in Fig.10 are unable to show whether the convergence is from above or below,because of the taking of absolute values in (41).That visualization deficiency should be kept in mind should frequency bounding be important.The plots clearly show at a glance that,for the sameNe,BLCD roughly doubles the number of correct digits provided by the other instances.It also illustrates that CMM and LDMM gave the same results,while SDMM is the worst performer.

9 Two-Dimensional Elements

This section extends some of the previous ideas to two-dimensional plane stress elements.A comprehensive symbolic Fourier analysis of higher order 2D elements of arbitrary geometry is presently beyond the capability of computer algebra systems.In the interest of brevity,the study is limited to a simple vibration benchmark,from which some preliminary conclusions are drawn.

9.1 The(n-1)/n Conjecture

For bar elements withn≤4 nodes,the optimal combination of the consistent and HRZlumped mass matrices (in the sense of best fit to the acoustic branch) has been found to be

For the linear element (n=2) this result was first presented in [21] and corroborated in [22].For the quadratic element (n=3) see [9,11].For the present cubic element (n=4) see Section 6.1.Whether (42) holds forn＞4 has not been investigated.Confirmation would have limited practical value,however,since shape functions beyond cubic are rarely used in dynamics.It does suggest,however,that lumped mass matrices perform relatively better as the order rises.More interesting is the question:Can (42) be extended to two dimensional elements? The numerical study that follows suggest that the extension is only partly correct.

9.2 Cantilever Vibration Benchmark

The benchmark 2D problem concerns the lowest fundamental frequencies of the slender 8:1 cantilever beam shown in Fig.11a.The beam has a narrow rectangular cross section.It is fabricated of isotropic material withE=2880,ν=0,andρ=1.The dimensions shown in the Figure are:beam heightH=1,beam spanL=8,and cross section widthh=1.

Figure 1:Cubic elements with side nodes:(a) 4-node bar;(b) 12-node serendipity quadrilateral;(c) 16-node Lagrangian quadrilateral;(d) 32-node serendipity brick;(e) 64-node Lagrangian brick.Non-corner nodes drawn at the thirdpoints for convenience.For (d,e) only visible nodes are shown

Figure 2:Configuration of the four-node bar (Bar4) element,showing local node numbering as well as geometric,constitutive and fabrication properties.Side Nodes 3–4 are symmetrically placed about the center,at distances ±γ l/2

Figure 3:Spurious modes in RI mass and stiffness matrices and their interelement propagation features:(a) Normalized null eigenvector of MeC(γ,3) over individual element;(b) Computed null eigenvector over a 4-element mass discretization of free-free bar modeled with MeC(γ,3);(c)Normalized null eigenvectors (symmetric and antisymmetric) of MeC(γ,2) over individual element;(d) Computed null eigenvectors over a 4-element mass discretization of fixed-fixed bar modeled with MeC(γ,2);(e) Normalized spurious null eigenvector of Ke(γ,2) over individual element;(f) Computed null eigenvector over a 4-element mass discretization of fixed-fixed bar modeled with Ke(γ,2).Eigenvectors displayed in (b),(d),(f) computed using Mathematica Eigensystem function;those in (d),(f) are not orthonormalized

Figure 4:Schematics of plane wave fourier analysis:(a) Continuum bar model;(b) FEM discretization as Bar4 lattice;(c) Propagation of plane wave over lattice;(d) Extraction of twoelement,seven-node patch

Figure 5:Plane waves of various wavelengths that a Bar4 lattice may transport.Yellow and red circles mark corner and side nodes,respectively

Figure 6:Two-element Bar4 patch that shows:renumbered nodes,x-axis origin relocated to interface node 4,natural patch coordinate =x/ℓ,and the amplitude coefficients Bi of the plane wave expansion (29) that are nonzero at the indicated nodes

Figure 7:Dimensionless dispersion diagrams (DDD) of several Bar4 combinations of full and reduced-integration consistent-mass and stiffness matrices:(a) MeC(γ,4) and Ke(γ,3);(b) MeC(γ,3) and Ke(γ,3);(c) MeC(γ,4) and Ke(γ,2);(d) MeC(γ,3) and Ke(γ,2).The value of γ ∈(0,1)is irrelevant

Figure 8:Dimensionless dispersion diagrams (DDD) of optimal Bar4 DLMM and LCMM instances paired with the fully integrated stiffness matrix.(a) MeL(1/,4) and Ke(γ,3);(b)MeL(1/,4) and Ke(γ,3);(c) MeLC(1/,4,3/4) and Ke(γ,3);(d) Detail of (c) Near k=Ω=p (Nyquist frequency and wavelength)

Figure 9:Fixed-free homogeneous prismatic elastic bar member used in vibration test for several Bar4 mass matrix instances.The figure shows a 2-element discretization with 7 nodes.Side nodes 2,3,5,and 6 are drawn at the lobatto points

Figure 10:Performance of selected Bar4 template instances in predicting the first three natural frequencies ωi,i=1,2,3 of the fixed-free prismatic homogeneous bar shown in Fig.9a.This is a graphical,log-log representation of the results of Table 10.Horizontal axis shows number of elements while vertical axis displays correct digits of computed frequency.See (41) for the definition of d

Figure 11:Cantilever vibration benchmark:(a) Structure;(b) A 16×2 FEM mesh idealization(nodes not pictured)

Figure 12:First six vibration modes for problem of Fig.11 obtained with a 32×4 Quad4 FEM discretization.Displacements grossly exaggerated for visibility

Figure 13:Accuracy of computed Quad4 frequencies for four meshes (8×1 through 64×8) as affected by parameter μ of the consistent-lumped mass combination MLC=(1-μ)MC+μML

Figure 14:Accuracy of computed Quad9 frequencies for four meshes (8×1 through 64×8),as affected by parameter μ of the consistent-lumped mass combination MLC=(1-μ)MC+μML

Figure 15:Accuracy of computed Quad16 frequencies for three meshes (8×1 through 32×4)with μ=3/4,as affected by parameter γ

Regular meshes of square elements are considered.The numberNxof elements along the span varies from 8 through 64 while the numberNythrough the height varies from 1 through 8.A 16×2 mesh is shown in Fig.11b.The root clamping condition is imposed by settinguxanduyto zero at all root nodes.

The first six lowest natural frequencies are listed in Table 11.They comprise four flexural(bending) modes with circular frequenciesω1,ω2,ω4andω5,as well as two bar-like (axial,longitudinal) modes:ω3andω6.It should be noted that flexural frequencies given by the Bernoulli-Euler and Timoshenko beam models differ markedly from those listed there because of 1D behavioral approximations as well as neglect of root clamping details.The associated mode shapes are pictured in Fig.12.

Table 11:Benchmark cantilever frequencies

9.3 Plane Stress Finite Element Modeling

The benchmark cantilever problem was modeled with regular FEM plane stress meshes ranging from 8×1 through 64×8.Three standard Lagrangian isoparametric models were tested:

· Quad4:Taig bilinear quadrilateral,4 nodes,8 DOF.

· Quad9:biquadratic quadrilateral,9 nodes,18 DOF.

· Quad16:bicubic quadrilateral,16 nodes,32 DOF,pictured in Fig.1c.

Frequencies for the first six lowest frequency modes (those depicted in Fig.12) were compared to the benchmark frequencies in Table 11.The relative accuracy is measured by

in whichandωiare the computed and benchmark frequencies,respectively,andiranges from 1 to 6.The valuedindicates the number of significant digits obtained.As noted in Section 8 this measure cancels out the error sign.However,in the present benchmark the convergence is often (but not always) monotonic.

9.4 Plane Stress Finite Element Results

For the Quad4 models,the mass matrix is taken as the LCMM

in which MCis computed with full integration and MLis its HRZ-lumped form.Six values ofμare tested:0 (consistent mass),1/4,1/2,2/3,3/4,and 1 (lumped mass).The results are presented in Fig.13.For the two bar-like modes (#3 and #6) the superior convergence ofμ=1/2 is obvious,in agreement with (42).For the four flexural modes,the convergence rate is roughly similar for allμtested,withμ=1 (lumped mass) consistently providing the best result andμ=0 (consistent mass) the worst.The advantage is insignificant for the fundamental mode but becomes more visible as the frequency increases.

For Quad9 the LC combination (44) is again selected with full integration.The same six values ofμused for Quad4 are tested.The results are presented in Fig.14.For the two bar-like modes the superior convergence ofμ=2/3 is obvious,in agreement with (42).For the flexural modes,the effect ofμis barely discernible.

On finally passing to Quad16,(44) is also assumed with full integration.The flexural mode results were found to be barely sensitive toμin the [0,1] range.To save space only the results forμ=3/4 are shown in Fig.15.An additional parameter now appears:γ,which localizes the 8 side and 4 interior nodes.Four values ofγwere tested:and 1/2.As before,the performances for bar-like and flexural modes are quite different.The combinationμ=3/4 andγ=(nodes at Lobatto points) was found to be optimal in Section 6.2.The results for the bar-like modes 3 and 6 verify this in 2D.For flexural modes the optimal value is unclear,althoughγ=appears to perform slightly better,followed byγ=1/3.

Three comments are in order for the Quad16 results.First,the flexural-mode convergence lines in Fig.15 appear more “spread out” than those in Figs.13 and 14.The reason is that thedscale is compressed on account of the high accuracy provided by the coarsest 8×1 mesh,which yields 5 to 3 exact digits for the first 4 flexural modes (Note that Quad4 gives less than one).Second,the jagged appearance of some plots is due to nonmonotonicity (onlyμ=0 would guarantee monotonic convergence).Third,theωiaccuracy of the eigensolver levels out at approximately 12 digits,which explains the incipient “flattening” of the bar-like modes.

10 Conclusions

The main results of this paper are summarized next.

· Presents the first detailed study,using plane wave packets as test functions,of the effect of adjusting side node locations of Lagrangian cubic finite elements on low-frequency accuracy.

· The combined effect of reduced integration of the mass and/or stiffness matrix is considered.

· For cubic bar elements,acoustic branch superconvergence is obtained by placing the side nodes at Lobatto points,and selecting a special combination of the mass and stiffness matrix.One remarkable result:6 digits inω1withone element,Table 10.

· For Lagrangian cubic plane stress elements,superconvergence is retained for bar-like vibration modes.This behavior continues in regular 2D meshes for bar-like modes,and pressumably in 3D.But it disappears for modes that involve rotation,flexure and shear.This is not surprising since axial plane wave packets (P-waves) are plainly insufficient;S-waves should be also included as test functions in 2D and 3D Fourier analysis.

Whether superconvergent stiffness-mass combinations can be constructed for general 2D and 3D motions using the general theory of mass-stiffness templates [9] is an open question.Challenges associated with symbolic analysis in multiple dimensions are discussed in Appendix E of that paper.

Funding Statement:This paper expands on work conducted during the 2005–2006 summer academic recesses while the author was a visitor at CIMNE (Centro Internacional de Métodos Numéricos en Ingenieria) at Barcelona,Spain.The visits were partly supported by fellowships awarded by the Spanish Ministerio de Educación y Cultura during May-June of those years,and partly by the National Science Foundation under grant High-Fidelity Simulations for Heterogeneous Civil and Mechanical Systems,CMS-0219422.Special thanks are due to Professor Loc Vu-Quoc for organizing the Special CMES Issue in celebration of Karl Pister’s 95th birthday.

Conflicts of Interest:The author declares that he has no conflict of interest regarding the present study.

Appendix A

The template approach to finite element construction gradually emerged from a series of serendipitous discoveries.The starting point was work of the author with Pål Bergan in 1982 while he was on a sabbatical at Stanford.The work focused on the development of high performance 3-node triangle elements for plane stress [24] and Kirchhoff plate bending [25] using the Free Formulation (FF) of Bergan and Nygård [6,26,27].The elements were eventually combined to form a shell element and implemented into the corotational nonlinear program FENRIS used by Der Norske Veritas for analyzing marine structures under extreme conditions.The plane stress element was the first 3-node membrane triangle with drilling DOF that was rank sufficient and free of aspect ratio locking.

As the name implies,FF allows substantial freedom in element development.In particular conformity can be ignored;all that matters is satisfying the Individual Patch Test that can be tested on a single element.Ensuing developments showed that free parameters naturally appeared in the stiffness matrix—as well as the mass and geometric stiffness matrices when the formulation was extended to dynamic and stability analysis,respectively.Which values should be given?Obviously the decision cannot be left to the discretion of program users,who might have no idea of the background theory.They were assigned using various optimality criteria.Free parameters kept appearing in FF successors,such as the Assumed Natural Deviatoric Strain (ANDES)method [28–30].

From special cases a general approach eventually took shape.It was named the template method in a 1994 journal paper [2].A fuller historical account is provided in a tutorial chapter [31].The general concept of template as parametrized forms of FEM matrix equations is discussed in [10,11,32–34].The main advantage of templates is flexibility.This brings the attendent ability to customize elements without worries about complying with restrictions such as interelement conformity.The approach can generate all possible convergent elements;the main question left is which one to pick up.This rosy picture,however,is clouded by practical difficulties:

· Element development computations must be done symbolically,carrying along material,fabrication and geometric information as variables.Such manipulations lie beyond human powers,and require a computer algebra system (CAS).

· Symbolic computational work effort grows exponentially in the number of free parameters as well as matrix order.In addition,symbolic eigenvalue analysis in stability and dynamic problems rapidly becomes difficult or impossible as the chacteristic polynomial order increases.

With current CAS the treatment of 1D elements is possible,as well as that of 2D elements of simple geometry.But 2D and 3D elements of arbitrary geometry (for instance,doubly curved shells) lie outside the scope of current computers and CAS power.