Space-time Discontinuous Galerkin Method Based on a New Generalized Flux Vector Splitting Method for Multi-dimensional Nonlinear Hyperbolic Systems

P.A.Trapperand P.Z.Bar-Yoseph

1 Introduction

The Discontinuous Galerkin(DG)method stems from investigations of numerical solution of the linear neutron transport equation, first Reed and Hill(1973)and subsequently by Lesaint and Raviart(1974).Since then the method has been widely developed and analyzed,and used extensively in different fields ranging from computational fluid dynamics and acoustics to electro magnetics and elasticity.For an extensive overview of DG methods see Cockburn,Karniadakis,and Shu(2000),Hesthaven and Warburton(2008),Di Pietro and Ern(2012),and Feng,Karakashian,and Xing(2014),and the references therein.

DG methods involve discontinuous approximations over finite elements with weakly enforced connectivity.Consequently,these methods can easily handle irregular meshes,complex geometries and polynomial approximations of different degrees in different elements,making them suitable forhp-adaptivity.Moreover,they are stable,locally conservative,high-order accurate and highly parallelizable.The other advantage of these methods lies in their ability to accurately capture discontinuities,sharp gradients and shocks in the solution,making them attractive for high frequency response of the system.

In this paper we focus on a particular family of DG methods,known as explicit space-time DG,which treat time as an additional element dimension and assume the unknown fields to be discontinuous in time.These methods have another exclusive advantage in that they easily allow for unstructured meshes in the space-time domain(i.e.,different time steps may be used in different elements).This is because the time step is no longer governed by the smallest elements in the mesh via a CFL condition,thus reducing computational cost.This approach was introduced by Bar-Yoseph(1989),who expanded the basic idea of Lesaint and Raviart(1974)to multi-dimensional nonlinear and quasi-linear hyperbolic systems of equations with shock fronts.Flux vector splitting with an alternating sweep in the forward and backward space directions was used.In this algorithm the discontinuities of the split fluxes are weighted along all boundaries,resulting in a physically meaningful up winding effect.

Bar-Yoseph and Elata(1990)further developed this notion to provide an answer to the efficiency problem by moving the nodes to Gauss points,thus cutting down the number of operations needed.They were also able to reconstruct the exact solution of some problems by using titled elements,thus offering ana posteriorierror study.Bar-Yoseph,Elata,and Israeli(1993)offered a qualitative and quantitative presentation of the stability,dissipation and dispersion of this method.Later,Aharoni and Bar-Yoseph(1992)developed a new approach for the integration of governing nonlinear ODE’s in time.Zrahia and Bar-Yoseph(1994a)further generalized this using the time spectral element method,which is a high-order method with high numerical efficiency and a high degree of accuracy that has been subsequently successfully used[Ben-Tal,Bar-Yoseph,and Flashner(1995,1996);Bar-Yoseph,Fisher,and Gottlieb(1996a,b);Bar-Yoseph(1998);Weill,Shitzer,and Bar-Yoseph(1993);Zrahia and Bar-Yoseph,(1994b);Bar-Yoseph,Moses,Zrahia,and Yarin(1995);Naveh,Bar-Yoseph,and Halevi(1999)].

The flux vector splitting technique employed in Bar-Yoseph(1989)and Bar-Yoseph and Elata(1990)was developed by Steger and Warming(1981).This technique worked for systems in which the flux vectors are homogeneous functions of degree one of variables,e.g.,Euler equations.Yet not all hyperbolic systems satisfy this property.Among those that do not are shallow water equations and equations of elasticity,for example.In this study we enhance the method first proposed by Bar-Yoseph(1989),Bar-Yoseph and Elata(1990),and Bar-Yoseph,Elata,and Is-raeli(1993)with a new generalized technique for splitting the flux vector that is not limited to its specific properties.specifically,the homogeneity property of the flux vector is no longer required.This technique is based on the flux’s characteristic decomposition,thus extending the scope of applicability of Bar-Yoseph(1989),Bar-Yoseph and Elata(1990),and Bar-Yoseph,Elata,and Israeli(1993)to a wider range of problems,and particularly to equations of elasticity.At this point we should emphasize that this method has a particular advantage for problems in solid mechanics that involve large deformations in soft materials because,due to its Eulerian nature,this method uses a fixed mesh,and as opposed to its Lagrangian counterpart,no element distortions will occur.

This paper is organized as follows.The mathematical formulation is described in the second section.We present our computational model in the third section,together with a new generalized technique for flux vector splitting.The computational results of various problems based on models of vibrating string and vibrating rod are provided and discussed in the fourth section.The critical time required for solution to reach a breakdown in nonlinear problems is estimated analytically,based on the work of Lax(1964).Finally,the last section offers some conclusions.

2 Formulation

Let ??RN+1be an open space-time region with piecewise smooth boundary Γ.Let(x0,x1,...,xN)be the set of Cartesian coordinates of pointxin ?:x0denotes the temporal coordinate t,and(x1,x2,...,xN)are the spatial co-ordinates,where N is the number of space dimensions.denote the canonical basis vectors of RN+1and letbe the inward unit vector normal to Γ(summation convention on repeated indices operates unless specifically stated otherwise).For simplicity,we assume that

and Γ is an N-dimensional hypersurface,admitting the following decomposition:

and ? is the empty set.

We consider the following system ofmfirst-order non-linear differential equations:

where

Eq.(3)is a system of balance laws derived forubelonging to an open space ofRm.The flux vectors,fffi∈Rm,i=0(1)N,are non-linear functions ofu.The vector ggg ∈Rmis a source vector.

System(3)can be also rewritten in a quasi-linear form

whereAAAi(uuu)=?ufffi,i=0(1)Nare the correspondingm×mJacobian matrices.We assume AAAi( uuu)are defined such that(5)is a first-order hyperbolic system(i.e.it has real eigenvalues).

The present Initial Boundary Value Problem(IBVP)consists of finding a functionu,which satisfies(3)or(5)subject to the initial condition

together with boundary conditions of the form

Here,uuu0andbbbjare given functions andBBBjare given matrices.

Consider a hypersurfaceS(x)which divides the region ? into two subregions ?+and ??.Let ψ be a tensor-valued function which is continuous in ?+and ??,and has definite limits ψ+and ψ?asxapproaches a point on the hypersurfaceSfrom paths entirely within the regions ?+and ??,respectively.The surface is called a singular surface with respect to ψ if

A singular surface is said to be a wave front if and only if it coincides with the one induced by physics.Here,ψ+is the region ahead of the"in fl ow"direction,while ψ?is the region behind it.Ifuis discontinuous across a space-time hypersurfaceS,the integrated balance laws imply that the jump inuacrossSsatisfies in whichnis the space-time normal to the singular hypersurface.

Let L2(?),L2(Γi)denote the Hilbert spaces of real-valued functions,square integrable in the Lebesgue sense on ?,Γirespectively.The inner products formdimensional vector functions can be defined as

A weak form of the balance laws including jump terms can be written as follows:

3 Computational model

Let the given domain ? be replaced by a collection ?hof hyperbrick elements ?h,satisfying certain regularity conditions[Hughes(1987)].Suppose that the temporal and spatial domains are discretized by a uniform mesh of elements,i.e.,h0=htand hi=hx,i=1(1)N,wherehtandhxare the mesh parameters representing the element size in the time and space directions,respectively.Thus,the element aspect ratio can be defined asr=ht/hx.

We introduce the following space-time finite element space of admissible functions:

whereQkldenotes the space of polynomials on ?eof degreekin space andlin time,i.e.Vhis the space of piecewise polynomials with no continuity requirement across inter-element boundaries.

Letnnne=eei,i=0(1)Nbe the unit inward normal vector to the element boundary Γe.The element boundaries are defined by

where

here+and?are called the positive and negative element boundaries.

3.1 Flux Vector Splitting

Although techniques for multi-dimensional splitting of flux vectors have been more extensively utilized in finite difference computations,they are equally applicable to finite element methods.In order to obtain a numerically stable explicit scheme for solution of(3)or(5),it is useful to split the flux vector according to the direction of propagation of information in the space-time domain,e.g.according to the sign of the wave propagation speed in the differential equation.The theory of multidimensional splitting of flux vectors was developed by Steger and Warming(1981)and later used in Bar-Yoseph(1989),and Bar-Yoseph and Elata(1990).It involves an attempt to systematically stabilize finite difference schemes employed for inviscid gas dynamic equations.The approach was based on the homogeneity property of the Euler equations in combination with specific equations of state,which allowed splitting the Jacobian matrixAAAinto+AAAand?AAAwith respect to positive and negative eigenvalues and acquiring the positive and negative fluxes as a multiplication of the correspondent matrices with the variablesuuu.This is also true for linear systems with constant coefficients.During the last decades,various flux-splitting techniques have proposed[Toro(2009)],mostly developed for the Euler equations of gas dynamics and usually relying on the above homogeneity property.However,not all hyperbolic systems satisfy this property.Among those that do not are shallow water equations and equations of elasticity,for example.

In this chapter we introduce a generalized technique for flux vector splitting that is based on characteristic decomposition offand that no longer requires the above homogeneity property.Let,∈R1,j=1(1)m,i=1(1)N,be themeigenvalues of the(m×m)Jacobian matrix AAAi,called characteristic speeds,and∈Rm,j=1(1)m,i=1(1)Nthemcorresponding linearly independent right eigenvectors that express the corresponding characteristic directions.Physically,eigenvalues represent speeds of propagation of information.Speeds will be measured as positive in the direction of increasing x and as negative otherwise.

The total flux vector fffimay be decomposed with respect to characteristic directions in the following way:

where

are the flux components associated with the positive and negative direction of xirespectively.It should be noted that in the general case of nonlinear problems,the eigenvalues may change their signs from point to point in the space-time domain.This implies that flux splitting must be consequently performed at each iteration,so the positive and negative sets of fluxes will change their compound accordingly.Further we formulate the positive and negative fluxes in a quasi-linear form,which will be useful later in this paper

3.2 Discontinuous Galerkin

Two different types of styles can be used for equations and mathematical expressions.They are:in-line style,and display style.The discontinuous Galerkin finite element method of(5)is obtained by posing the following formulation on a finite dimensional subspaceVhof the space of admissible functions.specifically we seekuuuh∈VVVhsuch that

whereh(.,.)andhh.,.i denote the discrete inner products,and[[ψ]]=ψin?ψout,ψin,ψoutare the values of ψ at the element boundary on the inside and outside of the element,respectively,anddenotes the element flux.

Since the approximate solutionuhis uniquely determined by(19),it is possible to computeuhsuccessively on each element ?e∈?h,e=1(1)NE,starting at the’in fl ow’boundary Γ?where the initial boundary conditions are weakly imposed.Here the initial and boundary conditions,equations(6)and(7),are replaced by weak conditions on the flux vector components across Γ?.This means that we have to impose a set of conditions that are combinations of the physical variables instead of the physical variables themselves.We use the standard discontinuous Galerkin method in which the weighting functions vector is the same within the element domain and on the element boundary.

Alternately,the discontinuous Galerkin method can be generalized by the Petrov-Galerkin method in which not only the jump discontinuity and the residual terms may be weighted by different test functions,but also the test functions and the base functions are different.Recently,Han and Atluri(2014a,b)presented an approach,which blends the(Meshless Local Petrov Galerkin)MLPG Methods of Atluri(1998,2004)and the energy conservation laws of Noether(1918)and Eshelby(1951,1975),and showed that it converges much faster and leads to much better accuracies than the classical FEM based on the global weak forms of the Newtonian Momentum Balance Laws.

3.3 Computational aspects

In order to explain how the method works,we review the technique of constructing the element coefficient matrix and the right-hand-side vector,using the proposed scheme;further details may be found in Bar-Yoseph and Elata(1990).

The base functions are discontinuous both in space and time.We use the Gauss-Lagrange interpolation,where the base functions,are defined at the ngv points of the Gauss-Legendre quadrature of the master element,i.e.,()=δij,i,j+1(1)ngv,where,,...,are the integration point coordinates in each element.

We emphasize at this point that the approximation space of the Gauss-Lagrange interpolation is exactly the same as for the standard interpolation in which the nodes are located at the corners and midside points of the element(i.e.,the discretization error is the same),but the Gauss-Lagrange interpolation for DG is more computationally efficient[Bar-Yoseph and Elata(1990)].

Bilinear and biquadratic Lagrangian base functions are used in two dimensions over the master element.Similarly,in three dimensions trilinear and triquadratic Lagrangian bases are used over a brick.

We apply the Gaussian quadrature formula to define the discreteL2(?e)inner product as

Similarly,for the surface integral we define the discreteL2()inner product,as applied on the discontinuities in the split flux vectors,as

where±are the positive and negativemultiplicativematrices ofuwith respect toare the coordinates of the surface Gauss points,are the quadrature weights on,ngs is the number of nodal points per element face,andngsis the number of integration points on an element face.

Substituting(20-24)into equation(20)yields a matrix equation in the following form:

whereKKKVis the volume matrix;KKKSis the surface matrix;qrepresents the contribution of adjacent elements to the flux jump conditions;dis the unknown element values ofuh;and,the element surfaces,are singular surfaces with respect touh.As mentioned in section 3.1,in the general case of quasi-linear equations(5)the sign of every characteristic speed can be a function of(x,u).Therefore,it would be more natural to solve(25)by a semi-iterative technique with an alternating sweep in the forward and backwardx-directions.

The algorithm is composed of two sequential steps.In the first step,using the solution of the previous iteration for imposing weak conditions along the element boundaries,the solution within each element is directly defined by a LU decomposition or an iterative solver.In the second step,the element flux vectors are iteratively updated.Since the solution procedure follows an element-by-element iteration,the problems of bandwidth and front width associated with direct solvers do not arise.Consequently,storage and computation are not influenced by element or node numbering.

The following examples demonstrate the application of this procedure to linear and quasi-linear problems.

4 Numerical examples

In this section we provide examples of different vibrating string and elastic rod models in which discontinuities appear.All problems are solved by bilinear(BL)and biquadratic(BQ)discontinuous finite elements,since LQ and QL elements were found not to be effective[Bar-Yoseph(1989)].Based on the work of Lax(1964),we analytically estimate the critical time required for solution of nonlinear problems to reach a breakdown and compare this to the one obtained from numerical results.In addition,we calculate and analyze the rate of convergence.

4.1 Nonlinear string

In the following examples we consider the standing vibrations of a finite,continuous,and nonlinear string[Zabusky(1962)],stretched along thex-axis from 0 toL,fixed at its end points.The reference mass density of the string denoted by ρ0,T0is the tension of the string and φ(x,t)is the deflection of the string in the transversal direction.

The governing equation is given by

where

is the square effective wave speed[Zabusky(1962)],the nonlinearity of interest,is a reference wave speed,and ε and γ are real positive scalars.The initial and boundary conditions are as follows

The second-order equation(26)may be transformed to the following system of first-order equations

where

In general these variables do not have to have a definite physical interpretation but rather some combination of some physical variables.In the case of a vibrating string,choosing them in this way gives them a physical meaning of velocity and slope of the string,respectively.

The initial and boundary conditions become

The Jacobian matrix of(29)is

with the eigenvalues

and the corresponding eigenvectors

The total flux of the system can be obtained by integration

It should be emphasized that the similar flux obtained using the Steger and Warming(1981)technique,is essentially linearized.

To recapture the displacement φ fromu2we may use the Newton-Cotes numerical integration for every time levelj,in a trapezoidal sense for bilinear elements

and in a Simpson sense for biquadratic elements

where φi,jandui2,jare the values at the pointsi,jof the space-time domain andhxis the element size in the space direction,thus obtaining an exact integration of u2.

4.1.1 Linear spring

As a first example we consider the simplest case of linear infinite string(γ=0),so(27)becomes

The analytical solution for this IBVP is given by:

The expressions for the total flux(35),the right eigenvectors,and the positive and negative flux vectors were obtained using the technique described in section 3.1,and the multiplicative matrices(18)are presented in the Appendix.

In this particular example the Steger and Warming(1981)splitting technique yields identical results becausefis linear.For this example we use the following parameters,L=1,c0=1.

In this paper we concentrate only on a posteriori error analysis.In the example problems,the L2error norm is considered

whereuis the exact solution(or reference solution at dense mesh)anduhis the approximated solution.

The error as a function of the element aspect ratiorfor different elements is shown in Fig.3a(here htis increased while hxis held fixed;hx=20 for BL and BQ elements).

Figure 1:Linear spring:(a)displacement,(b)velocity and(c)slope for hx=1/20,r=1.

Figure 2:Linear spring:slope at t=3.0 for hx=1/5,r=1.

Figure 3:Linear spring:solution accuracy and convergence in L2norm as function of(a)r,and(b)hx.

Above a certain value ofr,the error is dominated by the temporal error,and the optimal rate of convergenceh1+1for BL andh1+2for BQ is recovered.Moreover the scheme is unconditionally stable(independent of Courant No.=|c|r).We can use BL or BQ elements with Courant No.=1 without affecting the solution accuracy.The rate of convergence is depicted in Fig.3b(here the element aspect ratioris held fixed=1/2).Again,the optimal rate of convergence is obtained.

4.1.2 Nonlinear spring withγ=1

As a second example we consider a nonlinear case with γ=1,so that(27)becomes:

The expressions for the total flux(35),the right eigenvectors,and the positive and negative flux vectors are obtained with the technique described in section 3.1,and themultiplicativematrices(18)are presented in the Appendix.

In this example we use the following parameters,ε=0.2,L=1,c0=1.

The critical time required for a solution to reach a breakdown in non-linear problems can be estimated analytically based on the work of Lax(1964)

Based on(42)the breakdown is expected to occur in this example atTcr～=4L2?(c0επ2)=2.0264.

Figure 4:γ=1:(a)displacement,(b)velocity and(c)slope for hx=1/20,r=1.

Figure 5:γ=1:(a)velocity and(b)slope at t=3.0 for hx=1/32,r=1.

Figure 6:γ=1:solution accuracy and convergence in L2norm as function of hx.

Fig.5 shows that biquadratic elements produce a steeper shock front and have less dissipation than bilinear ones.

Fig.6 shows the rate of convergence in a regular region.Since it is impossible to say anything about the rate of convergence in the vicinity of the shock front[Bar-Yoseph(1989)],in most of the examples the regions away from the shock front(at every regular point)are considered.The examples show that the full rate of convergence is achieved.In the BL case the convergence is even faster than expected:≈2.3.

4.1.3 Nonlinear spring withγ=2

As a third example we consider a nonlinear case with γ=2,so that(27)becomes:

The expressions for the total flux(35),the right eigenvectors,and the positive and negative flux vectors are obtained with the technique described in section 3.1,and themultiplicativematrices(18)are presented in the Appendix.

For this example we used the following parameters,ε=0.2,L=1,c0=1.

Based on(42),in this example the breakdown is expected to occur atTcr～=4L2?(c0επ2)=1.0132.

In this example the breakdown occurs faster than in the previous example.Therefore,we only discuss the solution up untilt=1.0.

Figure 7:γ=2:(a)displacement,(b)velocity and(c)slope for hx=1/20,r=1.

Figure 8:γ=2:(a)velocity and(b)slope at t=1.0 for hx=1/20,r=1.

Figure 9:γ=2:(a)velocity and(b)slope at t=3.0 for hx=1/40,r=1.

This example again shows that BQ elements produce a steeper shock front but also more wiggles.BL elements produce more dissipation,which damps the waves with the high frequencies that cause those wiggles.Fig.8-9 show that as we refine the mesh,the solution becomes more accurate,producing a steeper shock front and fewer oscillations.

4.2 Nonlinear elastic rod

4.2.1 Nonlinear elastic rod with clamped-free boundary conditions

Next we consider the longitudinal vibration of the clamped-free finite elastic rod of lengthLexposed to initial deflection φ0at the free edge.The reference mass density of the rod is denoted by ρ0,E0is the Young modulus,and φ(x,t)is the deflection of the rod in the longitudinal direction.The governing equation is the same equation(26),as in the previous example,but with a different effective wave speed

Now the governing equation is transformed,as in section 4.3,to the system of first-order equations as(29)and(30),whereu1andu2are the velocity in the longitudinal direction and the axial strain,respectively.

The initial and boundary conditions become

and(44)becomes

The expressions for the total flux,the right eigenvectors,and the positive and negative flux vectors are obtained with the technique described in section 3.1,and themultiplicativematrices(18)are presented in the Appendix.

For this example we use the following parameters,L=1,c0=1,φ0=0.2.

Figure 10:Nonlinear clamped-free elastic rod:(a)displacement,(b)velocity and(c)slope for hx=1/20,r=1.

Figure 11:Nonlinear clamped-free elastic rod:(a)velocity and(b)strain at t=3.0 for hx=1/20,r=1.

Figure 12:Nonlinear clamped-free elastic rod:strain in the middle of the rod as a function of time(0≤t≤10),hx=1/20,r=1.

Fig.12 shows that in this example the BQ elements are absolutely superior to the BL elements:they yield steeper solution profiles with less dissipation and a smaller amount of overshoot(away from boundaries)than in the BL case.The overshoot is greater for the BQ case only untilt=0.5.

Fig.13 demonstrates the rate of convergence of the problem in the whole spacetime domain.Because there are shock fronts in the solutions for both the linear and the nonlinear cases,we should not expect the full rate of convergence.In all these cases the rate of convergence is≈1.

In all these examples,10-15 iterations per time level were required for convergence to meet π/10?4error criteria.

Though the discontinuities in this example are present from the beginning,no breakdown occurs,as can be seen in Fig.10-11.This is compatible with the Lax(1964)prediction,in which for φ,xx(0)=0 the breakdown is not expected to occur,Tcr=∞.

Figure 13:Nonlinear clamped-free elastic rod:solution accuracy and convergence in L2norm as function of hx(0≤t≤3).

4.2.2 Nonlinear elastic rod with clamped-clamped boundary conditions

We still want to examine the emerging breakdown for the rod problem,so we change the initial and boundary conditions as in the vibrating string problem in the previous examples(31).In this example we use the following parameters,ε=0.5,L=1,c0=1,φ0=0.2.

Figure 14:Nonlinear clamped-clamped elastic rod:(a)displacement,(b)velocity and(c)strain for hx=1/20,r=1.

Based on(42)the breakdown is expected to occur atTcr～=2(1+ε)L2/(0.2c0επ2)=3.0396.From all these examples we can see that when ε=0 no breakdown will occur.Again it is evident that BQ elements produce a steeper shock front and clearly a greater amount of overshoot.Lower order elements produce more numerical dissipation and therefore can yield solution profiles with a smaller amount of overshoot compared to those of higher order elements.

Fig.15 also indicates a phase difference between the BL and BQ cases:The average location of the shock front for the BL case is to the right of the one for the BQ case.As we have seen in the previous example,BL elements produce more dispersion errors than BQ elements.Fig.16 shows the rate of convergence in a regular region.Again we can see that the full rate of convergence is achieved.

Figure 15:Nonlinear clamped-clamped elastic rod:(a)velocity and(b)strain at t=3.0 for hx=1/20,r=1.

Figure 16:Nonlinear clamped-clamped elastic rod:solution accuracy and convergence in L2norm as function of hx(0≤t≤1).

5 Concluding remarks

The main conclusions of this work are as follows:

1.The space-time discontinuous finite element algorithm enhanced with a new generalized technique for flux vector splitting has been proved to be an efficient algorithm in numerical approximation for IBVP,especially for those solutions that develop shock layers.Moreover,it has also been shown to be valid for problems in elasticity.

2.The proposed technique for the splitting of flux vectors based on characteristic decomposition is valid for any kind of flux.Its computational efficiency for other types of problems should be further investigated.

3.Based on a posteriori error analysis we found that for smooth solutions,the optimal rate of convergencehn+1was recovered inL2norm for both linear and nonlinear hyperbolic systems, while in some cases of bilinear elements even a slight super convergence was observed.For regions with discontinuities in the solution,the convergence is slower than optimal.

4.All the profiles have a typical overshoot,although the wave front is quite crisp.Lower order elements have a greater degree of numerical dissipation and therefore can yield solution profiles with a smaller amount of overshoot compared to those of higher order elements,though higher order elements produce a steeper shock front.

5.Numerical results demonstrate that the present method may be effective in suppressing spurious oscillations with short wavelength.It is concluded that a proper combination of anh?prefinement strategy can be a viable alternative to schemes equipped with artificial stabilizers.The exact jump conditions can be reconstructed by tilting the element boundaries in the direction of a shock wave front[Bar-Yoseph and Elata(1990)].

6.The breakdown time in all relevant examples was found to be in agreement with theoretical prediction[Lax(1964)].

7.This method has a particular advantage for problems in solid mechanics that involve large deformations in soft materials,because this method,being essentially Eulerian,uses the fixed mesh,and contrary to its Lagrangian counterpart,no element distortions will occur.

Aharoni,D.;Bar-Yoseph,P.(1992):Mixed finite element formulation in the time domain for solution of dynamic problems.Computational Mechanics,vol.9,pp.359–374.

Atluri,S.N.;Zhu,T.(1998):A new meshless local Petrov-Galerkin(MLPG)approach in computational mechanics.Comput.Mech.,vol.22,pp.117-127.

Atluri,S.N.(2004):The Meshless Local Petrov-Galerkin(MLPG)Method for Domain&Boundary Discretizations.Tech Science Press,665 pages.

Bar-Yoseph,P.(1989):Space-time discontinuous finite element approximations for multi-dimensional non-linear hyperbolic systems.Computational Mechanics,vol.5,pp.145-160.

Bar-Yoseph,P.Z.(1998):Novel spectral and finite element method for unsteady heat transfer problems,in:de Vahl Davis,G.,Arinc F.(Eds.),Advances in Computational Heat Transfer.Begell House:New York.

Bar-Yoseph,P.;Elata,D.(1990):An efficientL2Galerkin finite element method for multi-dimensional nonlinear hyperbolic systems.International Journal for Numerical Methods in Engineering,vol.29,pp.1229-1245.

Bar-Yoseph,P.;Elata,D.;Israeli,M.(1993):On the generalizedL2Galerkinfinite element method for linear hyperbolic equations.International Journal for Numerical Methods in Engineering,vol.36,pp.679-694.

Bar-Yoseph,P.Z.;Fisher,D.;Gottlieb,O.(1996a):Spectral element method for nonlinear temporal dynamical systems.Computational Mechanics,vol.18,pp.302–313.

Bar-Yoseph,P.Z.;Fisher,D.;Gottlieb,O.(1996b):Spectral element method for nonlinear spatio-temporal dynamics of an Euler Bernoulli beam.Computational Mechanics,vol.19,pp.136–151.

Bar-Yoseph,P.;Moses,E.;Zrahia,U.;Yarin,A.L.(1995):Space-time spectral element methods for one-dimensional nonlinear advection-diffusion problems.Journal of Computational Physics,vol.119,pp.62-74.

Ben-Tal,A.;Bar-Yoseph,P.Z.;Flashner,H.(1995):Optimal maneuver of a flexible arm by space-time finite element method.AIAA,Journal of Guidance,Control,and Dynamics,vol.18,pp.1459–1462.

Ben-Tal,A.;Bar-Yoseph,P.Z.;Flashner,H.(1996):Space-time spectral element method for optimal slewing of a flexible beam.International Journal for Numerical Methods in Engineering,vol.39,pp.3101–3121.

Cockburn,B.;Karniadakis,G.;Shu,C.W.(2000):Discontinuous Galerkin Methods:Theory,Computation and Applications.Lecture Notes in Computational Science and Engineering 11,Springer Verlag:Berlin.

Di Pietro,D.A.;Ern,A.(2012):Mathematical Aspects of Discontinuous Galerkin Methods.Springer Verlag:Berlin,Heidelberg.

Eshelby,J.D.(1951):The Force on an Elastic Singularity.Phil.Trans.R.Soc.Lond.A,vol.244,pp.87-112.

Eshelby,J.D.(1975):The elastic energy-momentum tensor.Journal of Elasticity,vol.5,nos.3-4,pp.321-335.

Feng,X.;Karakashian,O.;Xing,Y.(2014):Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations:2012 John H.Barrett Memorial Lectures,Springer.

Han,Z.D.;Atluri,S.N.(2014a):Eshelby Stress Tensor T:a Variety of Conservation Laws for T in Finite Deformation Anisotropic Hyperelastic Solid&Defect Mechanics,and the MLPG-Eshelby Method in Computational Finite Deformation Solid Mechanics-Part I.CMES:Computer Modeling in Engineering&Sciences,vol.97,no.1,pp.1-34.

Han,Z.D.;Atluri,S.N.(2014b):On the(Meshless Local Petrov-Galerkin)MLPG-Eshelby Method in Computational Finite Deformation Solid Mechanics-Part II.CMES:Computer Modeling in Engineering&Sciences,vol.97,no.3,pp.119-237.

Hesthaven,J.S.;Warburton,T.(2008):Nodal Discontinuous Galerkin Methods:Algorithms,Analysis,and Applications.Springer Verlag:New York Inc.

Hughes,T.J.R.(1987):The Finite Element method:Linear Static and Dynamic Finite Element Analysis.Prentice-Hall:Englewood Cliffs,NJ.

Lax,P.D.(1964):Development of singularities of solutions of nonlinear hyperbolic partial differential equations.Journal of Mathematical Physics,vol.5,pp.611-613.

Lesaint,P.;Raviart,P.(1974):On a finite element method for solving the neutron transport equation,in:deBoor,C.(Ed.),Mathematical Aspects of finite Elements in Partial Differential Equations.Academic Press:New York,pp.89–145.

Naveh,Y.;Bar-Yoseph,P.Z.;Halevi,Y.(1999):Nonlinear modeling and control of a unicycle.Dynamics and Control,vol.9,pp.279-296.

Noether,E.(1918):Invariante Variationsprobleme.G?ttinger Nachrichten,Mathematisch-physikalische Klasse,vol.2,p.235.(English translation by M.A.Tavel,Transport Theory and Statistical Physics,vol.1,p.183,1971).

Reed,W.H.;Hill,T.R.(1973):Triangular mesh methods for the neutron transport equation.Los Alamos scientific Laboratory Report LA-UR-73-479,Los Alamos,NM.

Steger,J.L.;Warming,R.F.(1981):Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods.Journal of Computational Physics,vol.40,pp.263-293.

Toro,E.F.(2009):Riemann Solvers and Numerical Methods for Fluid Dynamics:A practical introduction.Third edition.Springer Verlag.

Weill,A.;Shitzer,A.;Bar-Yoseph,P.(1993):Finite element analysis of the temperature field around two adjacent cryo-probes.ASME Trans.Journal of Biomechanical Engineering,vol.115,pp.374-379.

Zabusky,N.J.(1962):Exact Solution for the Vibrations of a Nonlinear Continuous Model String.Journal of Mathematical Physics,vol.3,pp.1028-1039.

Zrahia,U.;Bar-Yoseph,P.(1994a):Space time spectral element method for solution of second-order hyperbolic equations.Computer Methods in Applied Mechanics and Engineering,vol.116,pp.135–146.

Zrahia,U.;Bar-Yoseph,P.(1994b):Alternative designs towards thermal optimization of coated valves using space-time finite elements.International Journal of Heat and Fluid Flow,vol.4,pp.189-206.

Appendix

In this appendix,for each example,we present the expressions for the total flux(35),right eigenvectors,the positive and the negative flux vectors,obtained with the technique described in the section 3.1,and themultiplicativematrices(18).

4.1.1.

4.1.2.

4.1.3.

4.2.