The exact solution of the eigenproblem for the parametric down conversion process in the Ke

H.G.GeorgiadisMechanics Division, National Technical University of Athens,1Konitsis Street, Zographou GR-15773,Greece e-mail:georgiad@central.ntua.grMem.ASME The Mode III Crack Problem in Microstructured Solids Governed by Dipolar Gradient Elasticity: Static and Dynamic AnalysisThis study aims at determining the elastic stress and displacementfields around a crack in a microstructured body under a remotely applied loading of the antiplane shear(mode III)type.The material microstructure is modeled through the Mindlin-Green-Rivlin dipo-lar gradient theory(or strain-gradient theory of grade two).A simple but yet rigorous version of this generalized continuum theory is taken here by considering an isotropic linear expression of the elastic strain-energy density in antiplane shearing that involves only two material constants(the shear modulus and the so-called gradient coefficient).In particular,the strain-energy density function,besides its dependence upon the standard strain terms,depends also on strain gradients.This expression derives from form II of Mindlin’s theory,a form that is appropriate for a gradient formulation with no couple-stress effects(in this case the strain-energy density function does not contain any rotation gradients).Here,both the formulation of the problem and the solution method are exact and lead to results for the near-tipfield showing significant departure from the predictions of the classical fracture mechanics.In view of these results,it seems that the conventional fracture mechanics is inadequate to analyze crack problems in microstructured materials. Indeed,the present results suggest that the stress distribution ahead of the tip exhibits a local maximum that is bounded.Therefore,this maximum value may serve as a measure of the critical stress level at which further advancement of the crack may occur.Also,in the vicinity of the crack tip,the crack-face displacement closes more smoothly as com-pared to the classical results.The latter can be explained physically since materials with microstructure behave in a more rigid way(having increased stiffness)as compared to materials without microstructure(i.e.,materials governed by classical continuum me-chanics).The new formulation of the crack problem required also new extended defini-tions for the J-integral and the energy release rate.It is shown that these quantities can be determined through the use of distribution(generalized function)theory.The boundary value problem was attacked by both the asymptotic Williams technique and the exact Wiener-Hopf technique.Both static and time-harmonic dynamic analyses are provided.͓DOI:10.1115/1.1574061͔1IntroductionThe present work is concerned with the exact determination of mode III crack-tipfields within the framework of the dipolar gra-dient elasticity͑or strain-gradient elasticity of grade two͒.This theory was introduced by Mindlin͓1͔,Green and Rivlin͓2͔,and Green͓3͔in an effort to model the mechanical response of mate-rials with microstructure.The theory begins with the very general concept of a continuum containing elements or particles͑called macromedia͒,which are in themselves deformable media.This behavior can easily be realized if such a macro-particle is viewed as a collection of smaller subparticles͑called micromedia͒.In this way,each particle of the continuum is endowed with an internal displacementfield,which is expanded as a power series in internal coordinate variables.Within the above context,the lowest-order theory͑dipolar or grade-two theory͒is the one obtained by retain-ing only thefirst͑linear͒term.Also,since these theories introduce dependence on strain and/or rotation gradients,the new material constants imply the presence of characteristic lengths in the ma-terial behavior,which allow the incorporation of size effects into stress analysis in a manner that the classical theory cannot afford. The Mindlin-Green-Rivlin theory and related ideas,after afirst development and some successful applications mainly on stress concentration problems during the sixties͑see,e.g.,Mindlin and Eshel͓4͔,Weitsman͓5͔,Day and Weitsman͓6͔,Cook and Weits-man͓7͔,Herrmann and Achenbach͓8͔,and Achenbach et al.͓9͔͒, have also recently been employed to analyze complex problems in materials with microstructure͑see,e.g.,Vardoulakis and Sulem ͓10͔,Fleck et al.͓11͔,Lakes͓12͔,Vardoulakis and Georgiadis ͓13͔,Wei and Huthinson͓14͔,Begley and Huthinson͓15͔,Exa-daktylos and Vardoulakis͓16͔,Huang et al.͓17͔,Zhang et al.͓18͔,Chen et al.͓19͔,Georgiadis and Vardoulakis͓20͔,Georgia-dis et al.͓21,22͔,Georgiadis and Velgaki͓23͔,and Amanatidou and Aravas͓24͔͒.More specifically,recent work by the author and co-workers͓13,20–23͔,on wave-propagation problems showed that the gradient approach predicts types of elastic waves that are not predicted by the classical theory͑SH and torsional surface waves in homogeneous materials͒and also predicts dispersion of high-frequency Rayleigh waves͑the classical elasticity fails to predict dispersion of these waves at any frequency͒.Notice that all these phenomena are observed in experiments and are also predicted by atomic-lattice analyses͑see,e.g.,Gazis et al.͓25͔͒.Contributed by the Applied Mechanics Division of T HE A MERICAN S OCIETY OFM ECHANICAL E NGINEERS for publication in the ASME J OURNAL OF A PPLIED M E-CHANICS.Manuscript received by the ASME Applied Mechanics Division,Apr.28,2002;final revision,Dec.19,2002.Associate Editor:B.M.Moran.Discussion onthe paper should be addressed to the Editor,Prof.Robert M.McMeeking,Depart-ment of Mechanical and Environmental Engineering University of California–SantaBarbara,Santa Barbara,CA93106-5070,and will be accepted until four months afterfinal publication of the paper itself in the ASME J OURNAL OF A PPLIED M ECHAN-ICS.Copyright©2003by ASMEJournal of Applied Mechanics JULY2003,Vol.70Õ517Thus,based on existing gradient-type results,one may conclude that the Mindlin-Green-Rivlin theory extends the range of appli-cability of continuum theories in an effort towards bridging the gap between classical͑monopolar or nongeneralized͒theories of continua and theories of atomic lattices.In the present work the concept adopted,following the afore-mentioned ideas,is to view the continuum as a periodic structure like that,e.g.,of crystal lattices,crystallites of a polycrystal or grains of a granular material.The material is composed wholly of unit cells͑micromedia͒having the form of cubes with edges of size2h.This size is therefore an intrinsic material length.We further assume͑and this is a rather standard assumption in studies applying the Mindlin-Green-Rivlin theory to practical problems͒that the continuum is homogeneous in the sense that the relative deformation͑i.e.,the difference between the macrodisplacementgradient and the microdeformation—cf.Mindlin͓1͔͒is zero andthe microdensity does not differ from the macrodensity.Then,weformulate the mode III crack problem by considering an isotropicand linear expression of the strain-energy density W.This expres-sion in antiplane shear and with respect to a Cartesian coordinatesystem Ox1x2x3reads Wϭ␮␧p3␧p3ϩ␮c(ץs␧p3)(ץs␧p3),where the summation convention is understood over the Latin indices,which take the values1and2only,(␧13,␧23)are the only iden-tically nonvanishing components of the linear strain tensor,␮is the shear modulus,c is the gradient coefficient͑a positive con-stant accounting for microstructural effects͒,andץs()ϵץ()/ץx s.The problem is two-dimensional and is stated in the plane(x1,x2).The above strain-energy density function is the simplest possible form of case II in Mindlin’s͓1͔theory and is appropriate for a gradient formulation with no couple-stress ef-fects,because W is completely independent upon rotation gradi-ents.Indeed,by referring to a strain-energy density function that depends upon strains and strain gradients in a three-dimensional body͑the Latin indices now span the range͑1,2,3͒͒,i.e.,a func-tion of the form Wϭ(1/2)c pqs j␧pq␧s jϩ(1/2)d pqs jlm␬pqs␬jlm with (c pqs j,d pqs jlm)being tensors of material constants and␬pqs ϭץp␧qsϵץp␧sq,and by defining the Cauchy͑in Mindlin’s nota-tion͒stress tensor as␶pqϭץW/ץ␧pq and the dipolar stress tensor ͑a third-rank tensor͒as m pqsϭץW/ץ(ץp␧qs),one may observe that the relations m pqsϭm p(qs)and m p[qs]ϭ0hold,where()and ͓͔as subscripts denote the symmetric and antisymmetric parts of a tensor,respectively.Accordingly,couple stresses do not appear within the present formulation by assuming dipolar͑internal͒forces with vanishing antisymmetric part͑more details on this are given in Section2below͒.A couple-stress,quasi-static solution of the mode-III crack problem was given earlier by Zhang et al.͓18͔. Note in passing that in the literature one mayfind mainly two types of approaches:In thefirst type͑couple-stress case͒the strain-energy density depends on rotation gradients and has no dependence upon strain gradients of the kind mentioned above ͑see,e.g.,͓11,17–19,23͔͒,whereas in the second type the strain-energy density depends on strain gradients and has no dependence upon rotation gradients͑see, e.g.,͓13,16,20–22͔͒.Exceptions from this trend exist of course͑see,e.g.,͓5–7͔͒and these works employ a more complicated formulation based on form III of Mindlin’s theory,͓1͔.Here,in addition to the quasi-static case,we also treat the time-harmonic dynamical case,which is pertinent to the problem ofstress-wave diffraction by a pre-existing crack in the body.In thelatter case,besides the standard inertia term in the equation ofmotion,a micro-inertia term is also taken into account͑in a con-sistent and rigorous manner by considering the proper kinetic-energy density͒and this leads to an explicit appearance of theintrinsic material length h.We emphasize that quasi-static ap-proaches cannot include explicitly the size of the material cell intheir governing equations.In these approaches,rather,a charac-teristic length appears in the governing equations only through the gradient coefficient c͑which has dimensions of͓length͔2)in the gradient theory without couple-stress effects or the ratio͑␩/␮͒͑which again has dimensions of͓length͔2)in the couple-stress theory without the effects of collinear dipolar forces,where␩is the couple-stress modulus and␮is the shear modulus of the ma-terial.Of course,one of the quantities c or͑␩/␮͒also appears within a dynamic analysis,which therefore may allow for an in-terrelation of the two different characteristic lengths͑the one in-troduced in the strain energy and the other introduced in the ki-netic energy—see relative works by Georgiadis et al.͓22͔and Georgiadis and Velgaki͓23͔͒.Indeed,by comparing the forms of dispersion curves of Rayleigh waves obtained by the dipolar ͑‘‘pure’’gradient and couple-stress͒approaches with the ones ob-tained by the atomic-lattice analysis of Gazis et al.͓25͔,it can be estimated that c is of the order of(0.1h)2,͓22͔,and␩is of the order of0.1␮h2,͓23͔.The mathematical analysis of the dynamical problem here pre-sents some novel features related to the Wiener-Hopf technique not encountered in dealing with the static case.The Wiener-Hopf technique is employed to obtain exact solutions in both cases,and also the Williams technique is employed for an asymptotic deter-mination of the near-tipfields.Also,since the gradient formula-tion exhibits a singular-perturbation character,the concept of a boundary layer is employed to accomplish the solution.On the other hand,the gradient formulation demands extended definitions of the J-integral and the energy release rate.It is further proved, by utilizing some theorems of distribution theory,that both energy quantities remain bounded despite the hypersingular behavior of the near-tip stressfield.Finally,physical aspects of the solution are discussed with particular reference to the closure of the crack faces and the nature of cohesive tractions.2Fundamentals of the Dipolar Gradient ElasticityA brief account of the Mindlin-Green-Rivlin theory,͓1–3͔,per-taining to the elastodynamics of homogeneous and isotropic ma-terials is given here.If a continuum with microstructure is viewed as a collection of subparticles͑micromedia͒having the form of unit cells͑cubes͒,the following expression of the kinetic-energy density͑kinetic energy per unit macrovolume͒is obtained with respect to a Cartesian coordinate system Ox1x2x3,͓1͔,Tϭ12␳u˙p u˙pϩ16␳h2͑ץp u˙q͒͑ץp u˙q͒,(1)where␳is the mass density,2h is the size of the cube edges,u p is the displacement vector,ץp()ϵץ()/ץx p,(˙)ϵץ()/ץt with t de-noting the time,and the Latin indices span the range͑1,2,3͒.We also notice that Georgiadis et al.͓22͔by using the concept of internal motions have obtained͑1͒in an alternative way to that by Mindlin͓1͔.In the RHS of Eq.͑1͒,the second term representing the effects of velocity gradients͑a term not encountered within classical continuum mechanics͒reflects the greater detail with which the dipolar theory describes the motion.Next,the following expression of the strain-energy density is postulated:Wϭ12c pqs j␧pq␧s jϩ12d pqs jlm␬pqs␬jlm,(2)where(c pqs j,d pqs jlm)are tensors of material constants,␧pq ϭ(1/2)(ץp u qϩץq u p)is the linear strain tensor,and␬pqsϭץp␧qs is the strain gradient.Notice that in the tensors c pqs j and d pqs jlm ͑which are of even rank͒the number of independent components can be reduced to yield isotropic constitutive relations.Such an isotropic behavior is considered here.Again,the form in͑2͒can be viewed as a more accurate description of the constitutive re-sponse than that provided by the classical elasticity,if one thinks of a series expansion for W containing higher-order strain gradi-ents.Also,one may expect that the additional term͑or terms͒will be significant in the vicinity of stress-concentration points where the strain undergoes very steep variations.Then,pertinent stress tensors can be defined by taking the variation of W518ÕVol.70,JULY2003Transactions of the ASME␶pq ϭץWץ␧pq,(3a )m pqs ϭץW ץ␬pqs ϵץWץ͑ץp ␧qs ͒,(3b )where ␶pq ϭ␶qp is the Cauchy ͑in Mindlin’s notation ͒stress tensor and m pqs ϭm psq is the dipolar ͑or double ͒stress tensor.The latter tensor follows from the notion of multipolar forces,which are antiparallel forces acting between the micro-media contained in the continuum with microstructure ͑see Fig.1͒.As explained by Green and Rivlin ͓2͔and Jaunzemis ͓26͔,the notion of multipolar forces arises rather naturally if one considers a series expansion for the mechanical power M containing higher-order velocity gra-dients,i.e.,M ϭF p u ˙p ϩF pq (ץp u ˙q )ϩF pqs (ץp ץq u ˙s )ϩ...,where F p are the usual forces ͑monopolar forces ͒within classical con-tinua and (F pq ,F pqs ,...)are the multipolar forces ͑dipolar or double forces,triple forces and so on ͒within generalized con-tinua.In this way,the resultant force on an ensemble of subpar-ticles can be viewed as being decomposed into external and inter-nal forces with the latter ones being self-equilibrating ͑see Fig.1͒.However,these self-equilibrating forces ͑which are multipolar forces ͒produce nonvanishing stresses,the multipolar stresses.Ex-amples of force systems of the dipolar collinear or noncollinear type are given,e.g.,in Jaunzemis ͓26͔and Fung ͓27͔.As for the notation of dipolar forces and stresses,the first index of the forces denotes the orientation of the lever arm between the forces and the second index the orientation of the pair of the forces;the same meaning is attached to the last two indices of the stresses,whereas the first index denotes the orientation of the normal to the surface on which the stress acts.The dipolar forces F pq have dimensions of ͓force ͔͓length ͔;their diagonal terms are double forces without moment and their off-diagonal terms are double forces with moment.The antisymmetric part F [pq ]ϭ(1/2)(x p F q Ϫx q F p )gives rise to couple stresses.Here,we do not consider couple-stress effects emphasizing that this is compat-ible with the particular choice of the form of W in ͑2͒,i.e.,a form dependent upon the strain gradient but completely independent upon the rotation gradient.Further,the equations of motion and the tractionboundary con-ditions along a smooth boundary can be obtained either from Hamilton’s principle ͑Mindlin ͓1͔͒or from the momentum balance laws and their application on a material tetrahedron ͑Georgiadis et al.͓22͔͒:ץp ͑␶pq Ϫץs m spq ͒ϭ␳u ¨q Ϫ␳h 23͑ץpp u¨q ͒,(4)n q ͑␶qs Ϫץp m pqs ͒ϪD q ͑n p m pqs ͒ϩ͑D l n l ͒n p n q m pqs ϩ␳h 23n r ͑ץr u ¨s͒ϭP s (n ),(5a )n q n r m qrs ϭR s (n ),(5b )where body forces are absent,D p ()ϭץp ()Ϫn p D (),D ()ϭn l ץl (),n s is the unit outward-directed vector normal to theboundary,P s(n )is the surface force per unit area ͑monopolar trac-tion ͒,and R s (n )is the surface double force per unit area ͑dipolar traction ͒.Finally,it is convenient for calculations to introduce another quantity,which is a kind of ‘‘balance stress’’͑see Eq.͑7͒below ͒,and is defined as␴pq ϭ␶pq ϩ␣pq ,(6)where ␣qs ϭ(␳h 2/3)(ץq u¨s )Ϫץp m pqs .With this definition,Eq.͑4͒takes the more familiar formץp ␴pq ϭ␳u ¨q .(7)Notice that ␴pq is not an objective quantity since it contains the acceleration terms (␳h 2/3)(ץq u ¨s ).These micro-inertia terms also are responsible for the asymmetry of ␴pq .This,however,does not pose any inconsistency but reflects the role of micro-inertia and the nonstandard nature of the theory.In the quasi-static case,where the acceleration terms are absent,␴pq is an objective tensor.On the other hand,the constitutive equations should definitely obey the principle of objectivity ͑cf.Eqs.͑9͒and ͑10͒below ͒.Now,the simplest possible form of constitutive relations is ob-tained by taking an isotropic version of the expression in ͑2͒in-volving only three material constants.This strain-energy density function readsW ϭ12␭␧pp ␧qq ϩ␮␧pq ␧pq ϩ12␭c ͑ץs ␧pp ͒͑ץs ␧qq ͒ϩ␮c ͑ץs ␧pq ͒͑ץs ␧pq ͒,(8)and leads to the constitutive relations␶pq ϭ␭␦pq ␧ss ϩ2␮␧pq ,(9)m spq ϭc ץs ͑␭␦pq ␧j j ϩ2␮␧pq ͒,(10)where ͑␭,␮͒are the standard Lame´’s constants,c is the gradient coefficient ͑material constant with dimensions of ͓length ͔2),and ␦pq is the Kronecker delta.Equations ͑9͒and ͑10͒written for a general three-dimensional state will be employed below only for an antiplane shear state.In summary,Eqs.͑4͒,͑5͒,͑9͒,and ͑10͒are the governing equa-tions for the isotropic dipolar-gradient elasticity with no couple bining ͑4͒,͑9͒,and ͑10͒leads to the field equation of the problem.Pertinent uniqueness theorems have been proved for various forms of the general theory ͑Mindlin and Eshel ͓4͔,Achenbach et al.͓9͔,and Ignaczak ͓28͔͒on the basis of positive definiteness of the strain-energy density.The latter restriction re-quires,in turn,the following inequalities for the material con-stants appearing in the theory employed here ͑Georgiadis et al.͓22͔͒:(3␭ϩ2␮)Ͼ0,␮Ͼ0,c Ͼ0.In addition,stability for the field equation in the general inertial case was proved in ͓22͔and to accomplish this the condition c Ͼ0is a necessary one ͑we notice incidentally that some heuristic gradient-like approaches not employing the rigorous Mindlin-Green-Rivlin theory appeared in the literature that take a negative c —their authors,unfortu-nately,do not realize that stability was lost in their field equation ͒.Finally,the analysis in ͓22͔provides the order-of-magnitude esti-mate (0.1h )2for the gradient coefficient c ,in terms of the intrin-sic material length h.Fig.1Monopolar …external …and dipolar …internal …forces act-ing on an ensemble of subparticles in a material with micro-structureJournal of Applied MechanicsJULY 2003,Vol.70Õ5193Formulation of the Quasi-Static Mode III Crack Problem,the J -Integral,and the Energy Release RateConsider a crack in a body with microstructure under a quasi-static antiplane shear state ͑see Fig.2͒.As will become clear in the next two sections,the semi-infinite crack model serves in a boundary layer type of analysis of any crack problem provided that the crack faces in the problem under consideration are trac-tion free.It is assumed that the mechanical behavior of the body is determined by the Eqs.͑4͒,͑5),(9),and ͑10͒of the previous section.An Oxyz Cartesian coordinate system coincident with the system Ox 1x 2x 3utilized previously is attached to that body,and an antiplane shear loading is taken in the direction of z -axis.Also,a pure antiplane shear state will be reached,if the body has the form of a thick slab in the z -direction.In such a case,the follow-ing two-dimensional field is generated:u x ϭu y ϭ0,(11a )u z ϵw 0,(11b )w ϵw ͑x ,y ͒,(11c )and Eqs.͑8)–(10͒take the formsW ϭ␮͑␧xz 2ϩ␧yz 2͒ϩ␮cͫͩץ␧xz ץx ͪ2ϩͩץ␧xzץyͪ2ϩͩץ␧yzץxͪ2ϩͩץ␧yz ץyͪ2ͬ,(12)␶xz ϭ␮ץw ץx ,(13a )␶yz ϭ␮ץw ץy,(13b )m xxz ϭ␮c ץ2wץx 2,(14a )m xyz ϭ␮cץ2wץx ץy ,(14b )m yxz ϭ␮c ץ2wץx ץy ,(14c )m yyz ϭ␮c ץ2wץy2.(14d )Further,͑4͒provides the equation of equilibriumץץx ͩ␶xz Ϫץm xxz ץx Ϫץm yxz ץy ͪϩץץy ͩ␶yz Ϫץm xyz ץx Ϫץm yyzץyͪϭ0,(15)which along with ͑13͒and ͑14͒leads to the following field equa-tion of the problem c ٌ4w Ϫٌ2w ϭ0,(16)where ٌ2ϭ(ץ2/ץx 2)ϩ(ץ2/ץy 2)and ٌ4ϭٌ2ٌ2.Finally,one may utilize ␴pq defined in ͑6͒for more economy in writing some equa-tions in the ensuing analysis.The antiplane shear components of this quantity are as follows:␴xz ϭ␮ͩץw ץx ͪϪ␮c ٌ2ͩץwץx ͪ,(17a )␴yz ϭ␮ͩץw ץy ͪϪ␮c ٌ2ͩץwץyͪ.(17b )Assume now that the cracked body is under a remotely applied loading that is also antisymmetric about the x -axis ͑crack plane ͒.Also,the crack faces are traction-free.Due to the antisymmetry of the problem,only the upper half of the cracked domain is consid-ered.Then,the following conditions can be written along the plane (ϪϱϽx Ͻϱ,y ϭ0):t yz ϵ␶yz Ϫץm xyz ץx Ϫץm yyz ץy Ϫץm yxzץxϭ0for ͑ϪϱϽx Ͻ0,y ϭ0͒,(18)m yyz ϭ0for ͑ϪϱϽx Ͻ0,y ϭ0͒,(19)w ϭ0for ͑0Ͻx Ͻϱ,y ϭ0͒,(20)ץ2wץy 2ϭ0for ͑0Ͻx Ͻϱ,y ϭ0͒,(21)where ͑18͒and ͑19͒directly follow from Eqs.͑5͒͑notice also that ͑18͒can be written as ␴yz Ϫ(ץm yxz /ץx )ϭ0by using the ␴pq quantity ͒,t yz is defined as the total monopolar stress,and ͑20͒together with ͑21͒always guarantee an antisymmetric displace-ment field w.r.t.the line of the crack prolongation.The definition of the stress t yz follows from ͑5a ͒.The problem described by ͑11)–(21͒will be considered by both the asymptotic Williams method and the exact Wiener-Hopf technique.Notice finally that no difficulty will arise by having zero boundary conditions along the crack faces since,eventually,the solution will be matched at regions where gradient effects are not dominant ͑i.e.,for x ӷc 1/2)with the K III field of the classical theory and in this way the remote loading will appear in the solution.Next,we present the new extended definitions of the J -integral and the energy release rate G .These definitions of the energy quantities are pertinent to the present framework of dipolar gradi-ent elasticity and to the aforementioned case of a crack in a quasi-static antiplane shear state.By following relative concepts from Rice ͓29,30͔,we first introduce the definitionJ ϭ͵⌫ͩWdy ϪP ¯z(n )ץw ץx d ⌫ϪR ¯z(n )D ͩץw ץxͪd ⌫ͪ,(22)where ⌫is a two-dimensional contour surrounding the crack tip͑see Fig.2͒,whereas the monopolar and dipolar tractions P ¯z (n )and R ¯z (n )on ⌫are given asP ¯z (n )ϭn q ͑␶qz Ϫץp m pqz ͒ϪD q ͑n p m pqz ͒ϩ͑D l n l ͒n p n q m pqz ,(23a )R ¯z (n )ϭn p n q m pqz .(23b )In the above expressions,n p with components (n x ,n y )is the unit outward-directed vector normal to ⌫,the differential operators D and D p were defined in Section 2,W is the strain-energy density function given by ͑12͒,and the indices (l ,p ,q )take the values x and y only.Of course,the above expressions for the tractions on ⌫are compatible with Eqs.͑5͒.Further,it can be proved that the inte-gral in ͑22͒is path independent by following Rice’s,͓29͔,proce-dure.Path independence is of great utility since it permits alter-nate choices of integration paths that may lead to adirectFig.2A crack under a remotely applied antiplane shear load-ing.The contour ⌫surrounding the crack tip serves for the definition of the J -integral.520ÕVol.70,JULY 2003Transactions of the ASMEevaluation of J .We should mention at this point that ͑22͒is quite novel within the present version of the gradient theory ͑i.e.,a form without couple stresses ͒,but expressions for J within the couple-stress theory were presented before by Atkinson and Leppington ͓31͔,Zhang et al.͓18͔,and Lubarda and Markenscoff ͓32͔.In particular,the latter work gives a systematic derivation of conser-vation integrals by the use of Noether’s theorem.Finally,we no-tice that the way the J -integral will be evaluated below is quite different than that by Zhang et al.͓18͔.Indeed,use of the theory of distributions in the present work leads to a very simple way to evaluate J ͑see Section 7below ͒.As for the energy release rate ͑ERR ͒now,we also modify the classical definition in order to take into account a higher-order term that is compatible with the present strain-gradient frameworkG ϭlim⌬x →0͵0⌬x ͫt yz ͑x ,y ϭ0͒•w ͑x ,y ϭ0͒ϩm yyz ͑x ,y ϭ0͒•ץw ͑x ,y ϭ0͒ץyͬdx⌬x,(24)where ⌬x is the small distance of a crack advancement.Of course,any meaningful crack-tip field given as solution to an associated mathematical problem,should result in a finite value for the energy quantities defined above.Despite the strong singu-larity of the stress field obtained in Sections 5and 6,the results of Section 7prove that J and G are indeed bounded.4Asymptotic Analysis by the Williams MethodAs is well known,Williams ͓33,34͔͑see also Barber ͓35͔͒de-veloped a method to explore the nature of the stress and displace-ment field near wedge corners and crack tips.This is accom-plished by attaching a set of (r ,␪)polar coordinates at the cornerpoint and by expanding the stress field as an asymptotic series in powers of r .By following this method here we are concerned,in a way,only with the field components in the sharp crack at very small values of r ,and hence we imagine looking at the tip region through a strong microscope so that situations like the ones,e.g.,on the left of Fig.3͑i.e.,a finite length crack,an edge crack or a crack in a strip ͒appear to us like the semi-infinite crack on the right of this figure.The magnification is so large that the other surfaces of the body,including the loaded remote boundaries,ap-pear enough far away for us to treat the body as an ‘‘infinite wedge’’with ‘‘loading at infinity.’’The field is,of course,a com-plicated function of (r ,␪)but near to the crack tip ͑i.e.,as r →0)we seek to expand it as a series of separated variable terms,each of which satisfies the traction-free boundary conditions on the crack faces.In view of the above,we consider the following separated form w (r ,␪)ϭr ␻ϩ1u (␪),where the displacement satisfies ͑16͒.Fur-ther,if only the dominant singular terms in ͑16͒are retained,the PDE of the problem becomes ٌ4w ϭ0,where ٌ4ϭٌ2ٌ2ϭ(ץ2/ץr 2ϩ1/r ץ/ץr ϩ1/r 2ץ2/ץ␪2)2.Also,in view of the defini-tions of stresses as combinations of derivatives of w and by re-taining again only the dominant singular terms,the boundary con-ditions t yz (x ,y ϭϮ0)ϭ0and m yyz (x ,y ϭϮ0)ϭ0will give at ␪ϭϮ␲ͩץ2ץr 2ϩ1r 2ץ2ץ␪2ϩ1r 2ͪץwץ␪ϭ0,(25a )ͩ1r ץץr ϩ1r 2ץ2ץ␪2ͪw ϭ0.(25b )In addition,the pertinent antisymmetric solution ͑i.e.,with odd behavior in ␪͒to the equation ٌ4w ϭ0has the following general form:w ϭr ␻ϩ1͑A 1sin ͓͑␻ϩ1͒␪͔ϩA 2sin ͓͑␻Ϫ1͒␪͔͒,(26)where ␻is ͑in general ͒a complex number and (A 1,A 2)are un-known constants.Now,͑25͒and ͑26͒provide the eigenvalue prob-lem͑␻ϩ1͒cos ͓͑␻ϩ1͒␲͔•A 1Ϫ3͑␻Ϫ1͒cos ͓͑␻Ϫ1͒␲͔•A 2ϭ0,(27a )͑␻ϩ1͒sin ͓͑␻ϩ1͒␲͔•A 1ϩ͑␻Ϫ3͒sin ͓͑␻Ϫ1͒␲͔•A 2ϭ0.(27b )For a nontrivial solution to exist,the determinant of the coeffi-cients of (A 1,A 2)in the above system should vanish and this gives the result:sin(2␻␲)ϭ0⇒␻ϭ0,1/2,1,3/2,2,....Next,by observing from ͑12͒that the strain-energy density W behaves at most as (ץ2w /ץr 2)or,by using the form w (r ,␪)ϭr ␻ϩ1u (␪),no worse than r ␻Ϫ1,we conclude that the maximum eigenvalue al-lowed by the integrability condition of the strain-energy density is ␻ϭ1/2.The above analysis suggests that the general asymptotic solu-tion is of the form w (r ,␪)ϭr 3/2u (␪),which by virtue of ͑26͒and ͑27b ͒becomesw ͑r ,␪͒ϭAr 3/2͓3sin ͑␪/2͒Ϫ5sin ͑3␪/2͔͒,(28)where A ϵϪA 1and the other constant in ͑26͒is given by ͑27b ͒as A 2ϭ3A 1/5.The constant A ͑amplitude of the field ͒is left un-specified by the Williams technique but still the nature of the near-tip field has been determined.Finally,the total monopolar stress has the following asymptotic behavior:t yz ͑x ,y ϭ0͒ϭO ͑x Ϫ3/2͒as x →ϩ0.(29)This asymptotic behavior will also be corroborated by the results of the exact analysis in the next section.5Exact Analysis by the Wiener-Hopf MethodAn exact solution to the problem described by ͑11͒–͑21͒will be obtained through two-sided Laplace transforms ͑see,e.g.,van der Pol and Bremmer ͓36͔and Carrier et al.͓37͔͒,theWiener-Fig.3William’s method:the near-tip fields of …i …a finite length crack,…ii …an edge crack,and …iii …a cracked strip correspond to the field generated in a body with a semi-infinite crackJournal of Applied MechanicsJULY 2003,Vol.70Õ521。

On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes qXiangxiong Zhang a ,Chi-Wang Shu b ,*a Department of Mathematics,Brown University,Providence,RI 02912,United StatesbDivision of Applied Mathematics,Brown University,Providence,RI 02912,United Statesa r t i c l e i n f o Article history:Received 4March 2010Received in revised form 6July 2010Accepted 13August 2010Available online 18August 2010Keywords:Hyperbolic conservation laws Discontinuous Galerkin method Positivity preserving High order accuracyCompressible Euler equations Gas dynamicsFinite volume schemeEssentially non-oscillatory scheme Weighted essentially non-oscillatory schemea b s t r a c tWe construct uniformly high order accurate discontinuous Galerkin (DG)schemes which preserve positivity of density and pressure for Euler equations of compressible gas dynam-ics.The same framework also applies to high order accurate finite volume (e.g.essentially non-oscillatory (ENO)or weighted ENO (WENO))schemes.Motivated by Perthame and Shu (1996)[20]and Zhang and Shu (2010)[26],a general framework,for arbitrary order of accuracy,is established to construct a positivity preserving limiter for the finite volume and DG methods with first order Euler forward time discretization solving one-dimen-sional compressible Euler equations.The limiter can be proven to maintain high order accuracy and is easy to implement.Strong stability preserving (SSP)high order time dis-cretizations will keep the positivity property.Following the idea in Zhang and Shu (2010)[26],we extend this framework to higher dimensions on rectangular meshes in a straightforward way.Numerical tests for the third order DG method are reported to dem-onstrate the effectiveness of the methods.Ó2010Elsevier Inc.All rights reserved.1.IntroductionIn this paper we are interested in constructing high order accurate schemes for solving hyperbolic conservation law sys-tems.For scalar conservation laws,the entropy solution is total variation diminishing (TVD),which is a desired property for numerical solutions.While traditional TVD schemes (e.g.[10])measure the total variation by that of the cell averages or grid values,leading to a necessary degeneracy of accuracy to first order at smooth extrema [18],genuinely high order TVD schemes can be constructed for one-dimensional scalar conservation laws [21,27]by measuring the total variation of the reconstruction polynomials.For multi-dimensional scalar conservation laws,it is difficult to enforce the TVD property for a high order scheme,however it is reasonable to insist on a strict maximum principle,which is satisfied by the entropy solu-tion.Genuinely high order accurate finite volume and discontinuous Galerkin (DG)schemes which satisfy a strict maximum principle have been constructed recently in [26].For hyperbolic conservation law systems,the entropy solutions in general satisfy neither the TVD property nor the max-imum principle.In this paper we are mainly interested in the Euler equations for the perfect gas,the one-dimensional ver-sion being given by0021-9991/$-see front matter Ó2010Elsevier Inc.All rights reserved.doi:10.1016/j.jcp.2010.08.016qResearch supported by AFOSR Grant FA9550-09-1-0126and NSF Grant DMS-0809086.*Corresponding author.Tel.:+14018632549;fax:+14018631355.E-mail addresses:zhangxx@ (X.Zhang),shu@ (C.-W.Shu).w t þf ðw Þx ¼0;t P 0;x 2R ;ð1:1Þw ¼qm E0B @1C A ;f ðw Þ¼m q u 2þp ðE þp ÞuB @1CAð1:2Þwithm ¼q u ;E ¼12q u 2þq e ;p ¼ðc À1Þq e ;where q is the density,u is the velocity,m is the momentum,E is the total energy,p is the pressure,e is the internal energy,and c >1is a constant (c =1.4for the air).The speed of sound is given by c ¼ffiffiffiffiffiffiffiffiffiffiffic p =q p and the three eigenvalues of the Jaco-bian f 0(w )are u Àc ,u and u +c .Physically,the density q and the pressure p should both be positive.We are interested in positivity-preserving high order schemes,which maintain the positivity of density and pressure at time level n +1,provided that they are positive at time level n .The techniques developed in this paper can be considered as generalizations of the maximum-principle-satisfying limiters in [26]and the positivity-preserving schemes in [20].We remark that failure of pre-serving positivity of density or pressure may cause blow-ups of the numerical algorithm,for example,for low density prob-lems in computing blast waves,and low pressure problems in the computing high Mach number astrophysical jets [9].We also remark that most commonly used high order numerical schemes for solving hyperbolic conservation law systems,including,among others,the Runge–Kutta discontinuous Galerkin (RKDG)method with a total variation bounded (TVB)lim-iter [2,4],the essentially non-oscillatory (ENO)finite volume and finite difference schemes [11,25],and the weighted ENO (WENO)finite volume and finite difference schemes [17,13],do not in general satisfy the positivity property for Euler equa-tions automatically.We now consider the Euler equations (1.1)in more detail.Let p ðw Þ¼ðc À1ÞE À1m 2qbe the pressure function.It can be easily verified that p is a concave function of w =(q ,m ,E )T if q P 0.For w 1=(q 1,m 1,E 1)T and w 2=(q 2,m 2,E 2)T ,Jensen’s inequality implies,for 06s 61,p s w 1þð1Às Þw 2ðÞP sp w 1ðÞþð1Às Þp w 2ðÞif q 1P 0;q 2P 0:ð1:3ÞDefine the set of admissible states byG ¼w ¼q m E 0B @1C Aq >0and p ¼ðc À1ÞE À12m 2q >08><>:9>=>;;then G is a convex set.If the density or pressure becomes negative,the system (1.1)will be non-hyperbolic and thus the ini-tial value problem will be ill-posed.We are interested in schemes for (1.1)producing the numerical solutions in the admissible set G .We start with a first order finite volume schemew n þ1j¼w n j Àk h w n j ;w n j þ1 Àh w n j À1;w n j h i;ð1:4Þwhere h (Á,Á)is a numerical flux,n refers to the time step and j to the spatial cell (we assume uniform mesh size only for sim-plicity),and k ¼D tD xis the ratio of time and space mesh sizes.w n j is the approximation to the cell average of the exact solution v (x ,t )in the cell I j ¼x j À1;x j þ1h i,or the point value of the exact solution v (x ,t )at x j ,at time level n .The scheme (1.4)and itsnumerical flux h (Á,Á)are called positivity preserving,if the numerical solution w n j being in the set G for all j implies the solu-tion w n þ1jbeing also in the set G .This is usually achieved under a standard CFL condition k kðj u j þc Þk 16a 0:ð1:5ÞExamples of positivity preserving fluxes include the Godunov flux [6],the Lax–Friedrichs flux [20],the Boltzmann type flux[19],and the Harten–Lax–van Leer flux [12].We now consider a general high order finite volume scheme,or the scheme satisfied by the cell averages of a DG method solving (1.1),which has the following formw n þ1j ¼w n j Àk h w Àj þ1;w þj þ12Àh w Àj À1;w þj À12h i;ð1:6Þwhere h is a positivity preserving flux under the CFL condition (1.5),w n j is the approximation to the cell average of the exact solution v (x ,t )in the cell I j ¼x j À12;x j þ12h iat time level n ,and w Àj þ12;w þj þ12are the high order approximations of the point values v x j þ1;t nwithin the cells I j and I j +1respectively.These values are either reconstructed from the cell averages w n j in a finitevolume method or read directly from the evolved polynomials in a DG method.We assume that there is a polynomial vectorX.Zhang,C.-W.Shu /Journal of Computational Physics 229(2010)8918–89348919q j (x )=(q j (x ),m j (x ),E j (x ))T (either reconstructed in a finite volume method or evolved in a DG method)with degree k ,wherek P 1,defined on I j such that w n j is the cell average of q j (x )on I j ,w þj À12¼q j x j À12and w Àj þ12¼q j x j þ12.A general framework to construct a high order positivity preserving finite volume scheme for the Euler equations wasintroduced in [20],in which a sufficient condition for the solution w n þ1jof (1.6)to be in the set G is that,all the nodal values w Æj þ12and w n þ1j Àa w Àj þ12þw þj À1 are in the set G under the CFL conditionk kðj u j þc Þk 16a a 0;where a 2(0,1]is a constant.Strong stability preserving (SSP)high order Runge–Kutta [25]and multi-step [24]time discret-ization will keep the positivity since G is convex.It is reasonable to require and easy to enforce the positivity of the point values w Æj þ1.However,it is more difficult to en-force the positivity of w n þ1j Àa w Àj þ12þw þj À12without destroying accuracy for an arbitrary high order scheme.We refer to Perthame and Shu [20]for more discussions on this point.In this paper,we provide a similar sufficient condition,whichis however much easier to enforce.We need the N -point Legendre Gauss–Lobatto quadrature rule on the intervalI j ¼x j À12;x j þ12h i,which is exact for the integral of polynomials of degree up to 2N À3.We would need to choose N such that2N À3P k .Denote these quadrature points on I j asS j ¼x j À12¼b x 1j ;b x 2j ;...;b x N À1j ;b x N j¼x j þ12n o:ð1:7ÞWe will prove that a sufficient condition for w n þ1j2G is simply q j b x a j2G for a =1,2,...,N ,under a suitable CFL condition.The same type of the linear scaling limiter used in [26]can enforce this sufficient condition without destroying accuracy.This limiter is also very easy to implement.Furthermore,we provide a straightforward extension of this result to arbitrary high order two-dimensional schemes on rectangular meshes.The main conclusion of this paper is,by adding a positivity preserving limiter which will be specified later to a high order accurate finite volume scheme or a discontinuous Galerkin scheme solving one or multi-dimensional Euler equations,with the time evolution by a SSP Runge–Kutta or multi-step method,we obtain a uniformly high order accurate scheme preserving the positivity in the sense that the density and pressure of the cell averages are always positive if they are positive initially.The paper is organized as follows:we first prove the positivity result for schemes in one space dimension in Section 2.In Section 3,we show a straightforward extension to two space dimensions on rectangular meshes.In Section 4,numerical tests for the third order DG method will be shown.Concluding remarks are given in Section 5.2.Positivity-preserving high order schemes in one dimension 2.1.A sufficient conditionWe consider the first order Euler forward time discretization (1.6);higher order time discretization will be discussed la-ter.Let b w a be the Legendre Gauss–Lobatto quadrature weights for the interval À12;12ÂÃsuch that P Na ¼1b w a ¼1,with2N À3P k .Motivated by the approach in [20,26],our first result isTheorem 2.1.For a finite volume scheme or the scheme satisfied by the cell averages of a DG method (1.6),if q j b x a j2G for all j and a ,then w n þ1j 2G under the CFL conditionk kðj u j þc Þk 16b w1a 0:ð2:1ÞProof.The exactness of the quadrature rule for polynomials of degree k impliesw n j¼1ZI jq j ðx Þdx ¼XN a ¼1b w a q j b x a j :By adding and subtracting h w þj À12;w Àj þ12,the scheme (1.6)becomes w n þ1j¼XN a ¼1b w a q j b x a jÀk h w Àj þ12;w þj þ1 Àh w þj À1;w Àj þ12þh w þj À1;w Àj þ12Àh w Àj À12;w þj À1h i ¼XN À1a ¼2b wa q jb x ajþb w N w Àj þ12Àk w N h w Àj þ12;w þj þ12 Àh w þj À12;w Àj þ12 h iþb w 1w þj À12Àk b w 1h w þj À12;w Àj þ12 Àh w Àj À12;w þj À12 h i ¼XN À1a ¼2b w a q j b x a jþb w N H N þb w 1H 1;8920X.Zhang,C.-W.Shu /Journal of Computational Physics 229(2010)8918–8934whereH1¼wþjÀ1Àkw1h wþjÀ1;wÀjþ12Àh wÀjÀ12;wþjÀ1h i;ð2:2ÞH N¼wÀjþ12Àkw Nh wÀjþ12;wþjþ1Àh wþjÀ1;wÀjþ12h i:ð2:3ÞNotice that(2.2)and(2.3)are both of the type(1.4),and b w1¼b w N,therefore H1and H2are in the set G under the CFLcondition(2.1).Now,it is easy to conclude that w nþ1jis in G,since it is a convex combination of elements in G.hRemark2.2.Here we only discuss Euler forward.Strong stability preserving high order Runge–Kutta[25]and multi-step [24]time discretization will keep the validity of Theorem2.1since G is convex.Remark2.3.From the proof of Theorem2.1,we can see that any type of quadrature rule will work as long as the quadrature points include the two cell ends and the quadrature is exact for polynomials of degree k.It would appear that there is a pos-sibility to achieve a larger CFL number if we canfind a better quadrature in the sense that b w1is larger.However,for k=2,3, we have verified that the Gauss–Lobatto quadrature is the best choice.Remark2.4.For the Lax–Friedrichsfluxhðu;vÞ¼1½fðuÞþfðvÞÀa0ðvÀuÞ ;where a0=k(j u j+c)k1,the CFL condition(1.5)was proven for k a0612in[20]by proving the numerical solution of thefirstorder Lax–Friedrichs scheme to be the cell average of the exact solution.Here,we prove that(1.5)for the Lax–Friedrichsflux can be relaxed to k a061.The Lax–Friedrichs scheme can be written asw nþ1j ¼w njÀk h w nj;w njþ1Àh w njÀ1;w njh i¼ð1Àk a0Þw njþk a0w njþ1À1f w njþ1!þk a0w njÀ1þ1f w njÀ1!:Assume w nj ;w njÀ1and w njþ1are in the set G,we want to show w nþ1j2G under the CFL k a061.Notice that w nþ1jis a convexcombination of the three vectors w nj ;w njþ1À1a0f w njþ1and w njÀ1þ1a0f w njÀ1,we only need to show w njÀ1þ1a0f w njÀ1andw njþ1À1a0f w njþ1are in the set G.It is easy to check that thefirst components of the both vectors are positive.The only non-trivial part is to check the positivity of the‘‘pressure”.For simplicity,we drop the subscripts and superscripts,i.e.,we provewÆ10fðwÞ2G if w2G.Let p¼1c EÀ1m2qand u=m/q.By a direct calculation,we havep wÆ1a0fðwÞ¼pð1Æua0Þq;ð1Æua0ÞmÆ1a0p;ð1Æua0ÞEÆua0pT"#¼1ÀpqcÀ12ða0ÆuÞ2!1Æua0pTherefore,p wÆ1a0fðwÞ>0()pqcÀ12ða0ÆuÞ<1()cpq<2ccÀ1ða0ÆuÞ2()ffiffiffiffiffiffiffic pqr<ffiffiffiffiffiffiffiffiffiffiffiffi2ccÀ1sða0ÆuÞ:Since c¼ffiffiffiffiffiffiffiffiffiffiffic p=qpand a0=k(j u j+c)k1,we have p wÆ1fðwÞ>0.The CFL condition(2.1)using the Lax–Friedrichsflux and the Gauss–Lobatto quadrature points for k=2,3,4,5are listed in Table2.1.We note that these conditions are comparable with and only slightly more restrictive than the standard CFL conditions for linear stability of discontinuous Galerkin methods[5].Table2.1The CFL condition(2.1)of Lax–Friedrichsflux for26k65and the Gauss–Lobatto quadrature points onÀ12;1 2ÂÃ.k CFL Quadrature points onÀ12;1 2ÂÃ2k a0616À12;0;12ÈÉ3k a0616À12;0;12ÈÉ4k a06112À12;À1ffiffiffiffi20p;1ffiffiffiffi20p;12n o5k a06112À12;À1ffiffiffiffi20p;1ffiffiffiffi20p;12n oX.Zhang,C.-W.Shu/Journal of Computational Physics229(2010)8918–893489212.2.A limiter to enforce the sufficient conditionGiven the vector of approximation polynomials q j (x )=(q j (x ),m j (x ),E j (x ))T ,either reconstructed for a finite volume scheme or evolved for a DG scheme,with its cell average w n j ¼q n j ;m n j ;E n jT 2G ,we would like to modify q j (x )into ~q j ðx Þsuch that it satisfiesAccuracy:For smooth solutions,the limiter does not destroy accuracyk ~qj ðx ÞÀq j ðx Þk ¼O ðD x k þ1Þ8x 2I j ;where kÁk denotes the Euclidean norm.Positivity:~q j b x a j2G for a =1,2,...,N . Conservativity:1xZI j~q j ðx Þdx ¼w n j :Definep n j¼ðc À1ÞE n j À12m n j 2=q n j .Then q n j >0and p nj >0for all j .Assume there exists a small number e >0such thatq n j P e and p nj P e for all j .For example,we can take e =10À13in the computation.The first step is to limit the density.Replace q j (x )byb q j ðx Þ¼h 1q j ðx ÞÀq n jþq n j ;ð2:4Þwhereh 1¼min q n j Àe q n jÀq min ;1();q min ¼min a q j b x a j :ð2:5ÞThen the cell average of b q j ðx Þover I j is still q n jand b q j b x a jP e for all a .The accuracy of b q j ðx Þcan be proven following the same lines as in [26].The second step is to enforce the positivity of the pressure.We need to introduce some notations.Letb q j ðx Þ¼ðb q j ðx Þ;m j ðx Þ;E j ðx ÞÞT and b q a j denote b q j b x a j.DefineG e¼w ¼q m E 1C A 0B @q P e and p ¼ðc À1ÞE À12m 2q P e 8><>:9>=>;;ð2:6Þ@G e ¼w ¼q m E 1C A 0B @q P eand p ¼ðc À1ÞE À1m 2q ¼e 8><>:9>=>;;ð2:7Þands aðt Þ¼ð1Àt Þw n jþt b q j b x a j;06t 61:ð2:8ÞG e is a convex set thanks to (1.3).o G e in (2.7)is a surface which contains part of the boundary of G e .s a (t )in (2.8)is the straight line passing through the two points w n j and b q j b x a j.If b q a j lies outside of G e ,namely p b q a j<e ,then there exists an intersection point of the straight line s a (t )and the surface @G e .Let s a e denote this intersection,then s a e ¼s a t a e ÀÁfor some t a e 2½0;1 satisfying p s a t a eÀÁÀÁ¼e .We will abuse the notation and let s a e ¼b q a j if p ðb q a j Þ2G e .So we haves ae ¼s a t a e ÀÁ;if p b q a j<e ;b q a j ;if p b q a jP e :8><>:ð2:9ÞWe consider the following new vector of polynomials~q j ðx Þ¼h 2b q j ðx ÞÀw n jþw n j ;ð2:10Þ8922X.Zhang,C.-W.Shu /Journal of Computational Physics 229(2010)8918–8934withh 2¼min a ¼1;2;...;N s a e Àw n jb q a j Àw n j:ð2:11ÞIt is easy to see that the cell average of ~q j ðx Þover I j is w n j .Next we would like to show the following lemma.Lemma 2.5.The ~q j ðx Þdefined in (2.10)and (2.11)satisfies ~q j b x a j2G e &G for all a .Proof.Notice that ~q j ðx Þis actually a convex combination of b q j ðx Þand w n j ,so the density of ~q j b x a jis no less than e .For the same reason,p ð~q j b x a j ÞP e if p b q j b x a jP e .If p b q j b x a j<e ,then p s a e ÀÁ¼e and ~qj b x aj¼h 2b q j b x a j Àw n j þw n j ¼h 2a e t a e b q j b x a j Àw n j þw n j h i þ1Àh 2a e w n j ¼h 2a e s a e þ1Àh 2a ew n j :Notice that s a e Àw n jb q a jÀw n j¼t a e .Therefore,(2.11)implies h 26t a e .So ~q j b x a j is a convex combination of s a e and w n j ,and thus p ~q j b x a j P e .h Finally,we need to show the limiter (2.10)and (2.11)does not destroy accuracy when q j (x )approximates a smooth solu-tion.Define d (z ,G )=min w 2G k z Àw k .Assume the exact solution v (x ,t n )is smooth and d (v (x ,t n ),G e )P M ,"x ,for some constant M >0.It suffices to show h 2=1+O (D x k +1).If h 2<1,then h 2¼s b e Àw n j =b q b j Àw n jfor some b where s b e is the intersection of thestraight line and the surface.Since w n j is a (k +1)th order approximation to the cell average of v (x ,t n ),we have d w n j ;G eP M þO ðD x k þ1ÞP M2if D x is smallenough.We can also assume the overshoot s b e Àb q b j ¼O ðD x k þ1Þsince b q b j is a (k +1)th order approximation to a point in G e .Thus,j 1Àh 2j ¼1Àh 2¼1Às b e Àw n j b q b j Àw n j ¼s b e Àb q b j b q bj Àw n j6s b e Àb q b jd w n j ;Ge ¼O ðD x k þ1Þ;where the third equality is due to the fact that b q b j ;s b e and w n j lie on the same line.Therefore,the limiting process (2.4),(2.5),(2.10)and (2.11)returns ~qj ðx Þsatisfying the accuracy,positivity and conservativity.2.3.Implementation for the DG methodAt time level n ,assuming the DG polynomial in cell I j is q j (x )=(q j (x ),m j (x ),E j (x ))T with degree k ,and the cell average of q j (x )is w n j ¼q n j ;m n j ;E nj T2G ,then the algorithm flowchart of our algorithm for the Euler forward is Set up a small number e ¼min j 10À13;q n j ;p w n jn o . In each cell,modify the density first:evaluate min a ¼1;...;N q j b x a j and get b q j ðx Þby (2.4)and (2.5),set b q j ðx Þ¼b q j ðx Þ;m j ðx Þ;E j ðx ÞÀÁT. Then modify the pressure:let b q a j denote b q j x a j ,for each a ,if p b q a j<e ,then solve the following quadratic equation for t a e ,p 1Àt a e ÀÁw nj þt a e b q j b x aj h i ¼e :ð2:12ÞIf p b q a jP e ,then set t a e ¼1.h 2in (2.11)is mathematically equivalent to h 2¼min a ¼1;...;N t a e .Get ~q j ðx Þby (2.10). Use ~qj ðx Þinstead of q j (x )in the DG scheme with Euler forward in time under the CFL condition (2.1).h For SSP high order time discretizations,we need to use the limiter in each stage for a Runge–Kutta method or in each stepfor a multi-step method.Remark 2.6.The implementation for a finite volume method is similar,but it will be a little bit more complicated for WENO since there are only nodal values but no polynomials in each cell after WENO reconstruction.One way to implement the limiter is to construct polynomials using the nodal values and cell averages,see [26]for details.We are also exploring other,simpler ways to implement this positivity preserving limiter for WENO finite volume schemes.These implementation details and numerical tests will be reported elsewhere.X.Zhang,C.-W.Shu /Journal of Computational Physics 229(2010)8918–89348923Remark 2.7.Theoretically,there is a complication regarding the CFL condition (2.1)for a Runge–Kutta time discretization.Consider the third order SSP Runge–Kutta method.To enforce (2.1)rigorously,we need to get an accurate estimation of k (j u j +c )k 1for all the three stages based only on the numerical solution at time level n ,which is highly nontrivial mathe-matically.In practice,we can simply multiply a factor,for example 2to 3,to the quantity k (j u j +c )k 1of w n ,as an estimation for all the stages.Although this is a rough estimation,it works well for us to choose a time step satisfying (2.1)in all the examples in Section 4.To be more efficient,we could implement this more stringent CFL condition only when a preliminary calculation to the next time step produces negative density or pressure.This complication does not exist if we use a SSP multi-step time discretization.3.Positivity-preserving high order schemes in two dimensions 3.1.A sufficient conditionIn this section we extend our result to finite volume or DG schemes of (k +1)th order accuracy on rectangular meshes solving two-dimensional Euler equationsw t þf ðw Þx þg ðw Þy ¼0;t P 0;ðx ;y Þ2R 2;ð3:1Þw ¼qm n EB B B @1C C C A;f ðw Þ¼mq u 2þp q u v ðE þp ÞuB B B @1C CC A;g ðw Þ¼nq u v q v 2þp ðE þp ÞvB B B @1C CC Að3:2Þwithm ¼q u ;n ¼q v ;E ¼12q u 2þ12q v 2þq e ;p ¼ðc À1Þq e ;where q is the density,u is the velocity in x direction,v is the velocity in y direction,m and n are the momenta,E is the totalenergy,p is the pressure,e is the internal energy.The speed of sound is given by c ¼ffiffiffiffiffiffiffiffiffiffiffic p =q p .The eigenvalues of the Jacobian f 0(w )are u Àc ,u ,u and u +c and the eigenvalues of the Jacobian g 0(w )are v Àc ,v ,v and v +c .The pressure function p is still concave with respect to w if q P 0and the set of admissible statesG ¼w ¼q m n E1C C C A 0B B B @q >0and p ¼ðc À1ÞE À12m 2q À12n 2q >08>>><>>>:9>>>=>>>;is still convex.For simplicity we assume we have a uniform rectangular mesh.At time level n ,we have a vector of approximation poly-nomials of degree k ,q ij (x ,y )=(q ij (x ,y ),m ij (x ,y ),n ij (x ,y ),E ij (x ,y ))T with the cell average w n ij ¼q n ij ;m n ij ;n n ij ;E nij Ton the (i ,j )cell x i À1;x i þ1h i Ây j À1;y j þ1h i .Letw þi À12;j ðy Þ;w Ài þ12;j ðy Þ;w þi ;j À12ðx Þ;w Ài ;j þ12ðx Þdenote the traces of q ij (x ,y )on the four edges respectively,see Fig.3.1.All of the traces are vectors of single variable polynomials of degree k.8924X.Zhang,C.-W.Shu /Journal of Computational Physics 229(2010)8918–8934。

Hindawi Publishing CorporationMathematical Problems in EngineeringVolume2011,Article ID745908,14pagesdoi:10.1155/2011/745908Research ArticleNumerical Investigations on SeveralStabilized Finite Element Methods forthe Stokes Eigenvalue ProblemPengzhan Huang,1Yinnian He,1,2and Xinlong Feng11College of Mathematics and System Sciences,Xinjiang University,Urumqi830046,China2Faculty of Science,Xi’an Jiaotong University,Xi’an710049,ChinaCorrespondence should be addressed to Yinnian He,heyn@Received27May2011;Accepted3July2011Academic Editor:Alexander P.SeyranianCopyright q2011Pengzhan Huang et al.This is an open access article distributed under the Creative Commons Attribution License,which permits unrestricted use,distribution,and reproduction in any medium,provided the original work is properly cited.Several stabilizedfinite element methods for the Stokes eigenvalue problem based on the lowest equal-orderfinite element pair are numerically investigated.They are penalty,regular, multiscale enrichment,and local Gauss integration parisons between them are carried out,which show that the local Gauss integration method has good stability,efficiency,and accuracy properties,and it is a favorite method among these methods for the Stokes eigenvalue problem.1.IntroductionIt is well known that numerical approximation of eigenvalue problems plays an important role in the analysis of the stability of nonlinear PDEs.Meanwhile,they are wildly used in many application areas:structural mechanics andfluid mechanics see 1 .Thus, development of an efficient and effective computational method for investigating these problems has practical significance and has drawn the attention of many peoples.At the time of writing,numerous works are devoted to these problems see 2–9 ,and the references cited therein .Yin et al. 10 derived a general procedure to produce an asymptotic expansion for eigenvalues of the Stokes problem on rectangular mesh.Chen and Lin 11 proposed the stream function-vorticity-pressure method for the eigenvalue problem.Rate of convergence estimates were derived for the approximation of eigenvalues and eigenvectors by Mercier et al. 12 .Xu and Zhou 13 proposed a two-grid discretization scheme for solving eigenvalue problems.Mixedfinite element methods are a natural choice for solvingfluid mechanics equations because these equations naturally appear in mixed form in terms of velocity and pressure 14,15 .In the analysis and practice of employing mixedfinite element methods in solving the Stokes equations,the inf-sup condition has played an important role because it ensures a stability and accuracy of the underlying numerical schemes.A pair offinite element spaces that are used to approximate the velocity and the pressure unknowns are said to be stable if they satisfy the inf-sup condition.Intuitively speaking,the inf-sup condition is something that enforces a certain correlation between twofinite element spaces so that they both have the required properties when employed for the Stokes equations.However,due to computational convenience and efficiency in practice,some mixedfinite element pairs which do not satisfy the inf-sup condition are also popular.Thus,much attention has been paid to the study of the stabilized method for the Stokes problem.Recent studies have focused on stabilization techniques,which include penalty 16–18 ,regular 19 ,multiscale enrichment 20 ,and local Gauss integration method 21,22 . There exist a lot of theoretical results for the stabilized mixedfinite element methods for the Stokes equations,and the comparisons between them are also given see 17,19,20,22–24 , and the references cited therein .In this paper,we mainly focus on the Stokes eigenvalue problem solved by these stabilizedfinite element methods based on the lowest equal-order finite element space pair.Moreover,we present the comparisons between these methods for the considered problem.The remainder of this paper is organized as follows.In the Section2,we introduce the studied Stokes eigenvalue problem,the notations,and some well-known results used throughout this paper.Then several stabilized mixedfinite element methods are reviewed, and their key stabilization techniques are recalled in parisons between these stabilized methods are performed numerically in Section4.Finally,we end with some short conclusions in Section5.2.PreliminariesLetΩbe a bounded,convex,and open subset of R2with a Lipschitz continuous boundary ∂Ω.We consider the Stokes eigenvalue problem as follows.Find u,p;λ ,such that−νΔu ∇p λu inΩ,div u 0inΩ,u 0on∂Ω,2.1where u u1 x ,u2 x represents the velocity vector,p p x the pressure,ν>0the viscosity,andλ∈R the eigenvalue.We will introduce the following Hilbert spaces:X H10 Ω 2,Y L2 Ω 2,M L20 Ωq∈L2 Ω :Ωq d x 0. 2.2The spaces L2 Ω m,m 1,2,are equipped with the L2-scalar product ·,· and L2-norm · L2 or · 0.The space X is endowed with the usual scalar product ∇u,∇v and the norm ∇u 0.Standard definitions are used for the Sobolev spaces W m,p Ω ,with the norm · m,p,m,p≥0. We will write H m Ω for W m,2 Ω and · m for · m,2.We define the continuous bilinear forms a ·,· and d ·,· on X×X and X×M, respectively,bya u,v ν ∇u,∇v ,∀u,v∈X,dv,qq,div v,∀v∈X,∀q∈M,2.3and a generalized bilinear form B ·,· ; ·,· on X×M × X×M byBu,p;v,qa u,v −dv,p−du,q,∀u,p,v,q∈X×M. 2.4With the above notations,the variational formulation of problem 2.1 reads as follows. Find u,p;λ ∈ X×M ×R with u 0 1,such that for all v,q ∈X×M,Bu,p;v,qλ u,v . 2.5Moreover,the bilinear form d ·,· satisfies the inf-sup condition 25 for all q∈M:sup v∈Xdv,q∇v 0≥β1q, 2.6whereβ1is a positive constant depending only onΩ.Therefore,the generalized bilinear form B ·,· ; ·,· satisfies the continuity property and inf-sup condition:Bu,p;v,q≤cν ∇u 0p∇v 0q,∀u,p,v,q∈X×M,supv,q ∈X×MBu,p;v,q∇v 0q≥β2ν ∇u 0p,∀u,p∈X×M,2.7where c andβ2are the positive constants depending only onΩ.3.Stabilized Mixed Finite Element MethodsFor h>0,we introduce thefinite-dimensional subspaces X h×M h⊂X×M,which are characterized by K h,a partitioning ofΩinto triangles K with the mesh size h,assumed to be uniformly regular in the usual sense.Then we defineX hu∈C0Ω2∩X:u|K∈P1 K 2,∀K∈K h,M hq∈C0Ω∩M:q|K∈P1 K ,∀K∈K h,3.1where P1 K represents the space of linear functions on K.Remark3.1 Nonconformingfinite element space .Denote the boundary edge byΓj ∂Ω∩∂K j and the interior boundary byΓjk Γkj ∂K j∩∂K k.Set the centers ofΓj andΓjk byζj and ζjk,respectively.The nonconformingfinite element space NC h for the velocity will be taken to beNC hv:v j v|Kj∈P1K j2,v jζjkv kζkj,vζj0,K j∈K h,∀j,k, 3.2where P1 K j is the set of all polynomials on K j of degree less than1.Note that the nonconformingfinite element space NC h is not a subspace of X.In this nonconforming case, the pair offinite element spaces is NC h×M h 26 ;that is,the conforming space is still used for pressure.Moreover,the discrete mixedfinite element formulation for the Stokes eigenvalue problem reads:find u h,p h;λh ∈ X h×M h ×R with u h 0 1,such that for all v,q ∈X h×M h,Bu h,p h;v,qλh u h,v . 3.3Next,we denote by U the array of the velocity and by P the array of the pressure.It is easy to see that 3.3 can be written in matrix form:A BB T0UPλhE000UP, 3.4where the matrices A,B,and E are deduced in the usual manner,using the bases for X h and M h,from the bilinear forms a ·,· ,d ·,· ,and ·,· ,respectively,and B T is the transpose of matrix B.Remark3.2.In the nonconforming case,we define the bilinear formsa u,v νK∈K h ∇u,∇v K,dv,pK∈K hdiv v,qK, 3.5and a generalized bilinear form B ·,· ; ·,· on X×M × X×M byBu,p;v,qa u,v −dv,p−du,q,∀u,p,v,q∈X×M. 3.6Accordingly,the discrete nonconforming formulation for the Stokes eigenvalue problem is:find u h,p h;λh ∈ NC h×M h ×R with u h 0 1,such that for all v,q ∈NC h×M h,Bu h,p h;v,qλh u h,v . 3.7Note that the lowest equal-order pair does not satisfy the discrete inf-sup condition:sup v h∈X h dv h,q h∇v h 0≥β3q h,or supv h∈NC hdv h,q h∇v h 0,h≥β3q h,∀q h∈M h, 3.8where the constantβ3>0is independent of h,and ∇v h 0,h is the discrete energy seminorm in the nonconforming case.In order to fulfill this condition,several ways have been used to stabilize the lowest equal-orderfinite element space pair.Method1 Penalty method .Find u h,p h;λh ∈ X h×M h ×R with u h 0 1,such that for all v,q ∈X h×M h,Bu h,p h;v,qενp h,qλh u h,v , 3.9whereε>0is a penalty parameter.The performance of this method obviously depends on the choice of the penalty parameterε.Then the matrix form of 3.9 can be expressed as⎡⎣A BB T ενD⎤⎦UPλhE000UP, 3.10where the matrix D is deduced in the usual manner,using the base for M h,from p h,q . Because the coefficient matrix of 3.10 is usually large and sparse,it is not easy to compute the numerical solution using a direct method.In general,one can use the Uzawa algorithm. By some simple calculation,we getE P i 1I−τB T A−λh E −1B−ενDE P i, 3.11where E PiP−P i andτis a positive constant.Method2 Regular method .Find u h,p h;λh ∈ X h×M h ×R with u h 0 1,such that for all v,q ∈X h×M h,Bu h,p h;v,q−δK∈K h∇p h−λh u h,∇qKλh u h,v , 3.12whereδ h2/ αν is a stabilization parameter andα>0.The regular method uses a simple way to stabilize the mixedfinite element approximation without a loss of accuracy.In fact,it can be treated the regular method as a special Douglas-Wang’s scheme 19 .Then the matrix form of this stabilized version can be expressed asA BB T−δD1UPλhE0F0UP, 3.13where additional blocks D 1and F correspond to the following respective terms:K ∈K h∇p h ,∇qK,−δK ∈K hu h ,∇qK.3.14As 3.11 ,we also haveE P i 1I −τB T−λh FA −λh E −1B δD 1 E P i .3.15Method 3 Multiscale enrichment method .Find u h ,p h ;λh ∈ X h ×M h ×R with u h 0 1,such that for all v,q ∈X h ×M h ,B u h ,p h ; v,q−δ1 K ∈K h∇p h −λh u h ,∇q K δ2e ∂K j ∩∂K kν∂n u h , ν∂n v e λh u h ,v ,3.16where δ1 h 2/ α1ν ,δ2 h/ α2ν are the stabilization parameters,α1,α2>0,n is the normal outward vector,∂n is normal derivative operator,and v denotes the jump of v across e .This stabilized method includes the usual Galerkin least squares stabilized terms on each finite element and positive jump terms at interelement boundaries.Moreover,a direct algebraic manipulation leads to the matrix formA BB Tδ2D 2−δ1D 1U PλhE 0F 0 UP,3.17where the matrix D 2is deduced in the usual manner,using the bases for X h ,from the terme ∂K j ∩∂K k ν∂n u h , ν∂n v e .As 3.11 ,we haveE P i 1I −τB Tδ2D 2−λh FA −λh E −1B δD 1E P i .3.18Method 4 Local Gauss integration method .Find u h ,p h ;λh ∈ X h ×M h ×R with u h 0 1,such that for all v,q ∈X h ×M h ,Bu h ,p h ; v,q −Gp h ,q λh u h ,v ,3.19where G p h ,q is defined byG p h ,qK ∈K hK,2p h q d x −K,1p h q d x,∀p h ,q ∈M h ,3.20whereK,i g x d x indicates a local Gauss integral over K that is exact for polynomials of degree i ,i 1,2.In particular,the trial function p h ∈M h must be projected to piecewiseconstant space when i 1.This stabilization technique is free of stabilization parameters and does not require any calculation of high-order derivatives,a specification of any mesh-dependent parameter or edge-based data structures.Then the corresponding matrix form of this stabilized method isAB B T −GU PλhE 000 UP, 3.21where the matrix G is deduced in the usual manner,using the bases for M h ,from the term G p h ,q .As 3.11 ,we getE P i 1I −τB TA −λh E −1B GE P i .3.22Remark 3.3.It is well known that if τis well chosen,then U i and P i converge,respectively,to U and P with a rate of convergence based on 3.11 , 3.15 , 3.18 ,and 3.22 .From these equations,we can find that Method 1converges faster than Method 4.Methods 2and 3,whose coe fficient matrices are not symmetric,may cost more time to converge.Remark 3.4.By using the regularity assumptions and well-established techniques for eigen-value approximation 1,8,10,12 ,the theoretical convergence rates should be of order O h 2 and O h for the velocity in the L 2-and H 1-norms,respectively,of order O h for the pressure in the L 2-norm,and of order O h 2 for the eigenvalue by using all these stabilized methods.4.Numerical ExperimentsIn this section we numerically compare the performance of the various stabilized mixed finite element methods discussed in the previous section.In all experiments,the algorithms are implemented using public domain finite element software 27 with some of our additional codes.Let the computation be carried out in the region Ω { x,y |0<x,y <1}.We considerthe Stokes eigenvalue problem in the case of the viscidity ν 1,and it will be numerically solved by the stabilized mixed methods on uniform mesh see Figure 1 .Here,we just consider the first eigenvalue of the Stokes eigenvalue problem for the sake of simplicity.The exact solution of this problem is unknown.Thus,we take the numerical solution by the standard Galerkin method P 2-P 1element computed on a very fine mesh 6742grid points as the “exact”solution for the purpose of comparison.Here,we take λ 52.3447as the first exact eigenvalue.As we know,the stabilized term of the regular and multiscale enrichment methodsmust be controlled by carefully designed stabilization parameters i.e.,δ,δ1,δ2 ,whose optimal values are often unknown.Hence,in Figures 2and 3we show the e ffect on the error of varying δ,δ1,and δ2on a fixed mesh h 1/64for the regular and multiscale enrichment approximations,respectively.An interesting thing can be observed that these two methods have an analogous convergence pattern with respect to the parameter δand δ1.Because they involve a similar stabilization term with respect to these two parameters,we note that the errors can become large with some values of the stabilization parameters.XY0.20.40.60.8100.20.40.60.81Figure 1:Uniform finite element partitioning of the unitsquare.10010−110−210−310−4L 2-e r r o r o f λh10−2100102104106108δFigure 2:E ffect of varying δfor the regular method.Results gotten from the penalty,regular,multiscale enrichment,local Gauss inte-gration,and nonconforming local Gauss integration methods are presented in Tables 1–5,respectively.Here,we choose ε 1.0e −5,α 8,α1 8,and α2 12.Because they can deal with the considered problem well.From Tables 1,2,3,4,and 5,we can see that these methods work well and keep the convergence rates just like the theoretical analysis except the multiscale enrichment method.Meanwhile,it can be seen that the penalty method requires the least CPU-time,which validates the analysis in Remark 3.3.As expected,we have an10010−110−210−310−2100102104106108δ1L 2-e r r o r o f λha10310210110010−110−210−310−410−2100102104106δ2L 2-e r r o r o f λhbFigure 3:E ffect of varying δ1 a and δ2 b for the multiscale enrichment method.Table 1:Results get from the penalty method with ε 1.0e −5.1/h CPU-time λh |λh −λ|/λRate 80.01660.26280.151269—160.09454.16880.0348472 2.1180240.23453.14260.015243 2.0393320.43752.79090.00852367 2.0205400.75152.62950.00544061 2.012048 1.15652.54330.00377371 2.006556 1.73452.48980.00277154 2.0023642.43752.45580.002122371.9986Table 2:Results get from the regular method with α 8.1/h CPU-time λh |λh −λ|/λRate 80.03156.72830.0837442—160.10953.48030.0216951 1.9486240.26652.85310.00971199 1.9822320.51552.63140.00547673 1.9913400.86152.52840.00350903 1.995048 1.31352.47230.00243822 1.996956 1.95352.43850.00179192 1.9979642.71952.41650.001372191.9986interesting observation that the error of the nonconforming local Gauss integration method is better than the conforming version,which is not surprising since the degrees of freedom of the nonconforming method are nearly three times than that of conforming one on uniform mesh see Figure 1 .Hence,it is natural that the nonconforming local Gauss integration method is more accurate and costs more CPU-time.Besides,to show the stability and e fficiency of these methods for the consideredproblem,we present the velocity streamlines and the pressure contours with h 1/64Table3:Results get from the multiscale enrichment method withα1 8andα2 12.1/h CPU-timeλh|λh−λ|/λRate 80.06367.82370.295713—160.28157.10620.0909634 1.7008 240.65654.73780.0457175 1.6968 32 1.26553.83040.0283839 1.6569 40 2.15653.38130.0198027 1.6133 48 3.35953.12290.0148665 1.5725 56 5.12552.95880.0117327 1.5357 647.18752.84720.00959935 1.5029 Table4:Results get from the local Gauss integration method.1/h CPU-timeλh|λh−λ|/λRate 80.03157.39510.096482—160.10953.62010.024366 1.9854 240.25152.91190.0108368 1.9983 320.46952.66380.00609553 2.0001 400.79752.54890.00390065 2.0006 48 1.23452.48650.00270843 2.0007 56 1.84452.44880.00198963 2.0008 64 2.54752.42440.00152315 2.0008Table5:Results get from the local Gauss integration method with the nonconforming element.1/h CPU-timeλh|λh−λ|/λRate 80.03150.21210.0407434—160.12551.73550.0116391 1.8076 240.34452.06190.00540181 1.8932 320.67252.18250.00309932 1.9311 40 1.15652.23970.00200547 1.9508 48 1.85952.27130.00140228 1.9624 56 2.79752.29050.00103505 1.9698 64 4.00152.30310.000795121 1.9749 in Figures5,6,7,8,and9.Meanwhile,we present the results by the standard Galerkin method P2-P1element computed on a veryfine mesh 6742grid points for the purpose of comparison in Figure4.From Figures5 a –9 a ,five resolved vortices are captured,which is consistent with that in Figure4 a .For the pressure,the penalty method is divergent from Figure5 b ,although it costs the least time.From Figures6 b -9 b ,the nonconforming local Gauss integration method shows the best numerical stability.5.ConclusionsIn this paper we have presented several stabilized mixedfinite element methods in solving the Stokes eigenvalue problem based on the lowest equal-orderfinite element space pair.By being compared numerically,we get some conclusions as follows.X Y0.20.40.60.8100.20.40.60.81 aY0.20.40.60.81X0.20.40.60.8154321−1−2−3−4−50.102798−0.108067 bFigure 4:Velocity streamlines a and pressure level lines b for the standard Galerkin method P 2-P 1element .X Y0.20.40.60.8100.20.40.60.81aY0.20.40.60.81X 00.20.40.60.814035302520151050−5−10−15−20−25−30−35−40bFigure 5:Velocity streamlines a and pressure level lines b for the penalty method with ε 1.0e −5.i The stability and e fficiency of the penalty method depend on the penalty parameter.The smaller this parameter,the more stable the method.However,if this parameter is too small,the condition number of the system matrix arising from this method will become too large to solve. ii The performance of the regular and multiscale enrichment method heavily dependson the choice of the stabilization parameters,which is a di fficult task in reality.Meanwhile,a poor choice of these stabilization parameters can also lead to serious deterioration in the convergence rates. iii The local integration method is free of stabilization parameters and showsnumerically the best performance among the methods considered for the given problem.X Y00.20.40.60.8100.20.40.60.81 a54321−1−2−3−4−50.107858−0.110432Y0.20.40.60.81X0.20.40.60.81bFigure 6:Velocity streamlines a and pressure level lines b for the regular method with α8.X Y0.20.40.60.8100.20.40.60.81 aY0.20.40.60.81X 00.20.40.60.810.1071354321−1−2−3−4−5−0.10957 bFigure 7:Velocity streamlines a and pressure level lines b for the multiscale enrichment method with α1 8and α2 12.iv From Tables 4and 5,we have an interesting observation that the value of λh byconforming method becomes small to converge to the exact solution and the one by nonconforming method becomes large to converge to the exact solution.We hold that it can obtain more accurate numerical solution by Lagrange interpolation between conforming and nonconforming results based on the same degrees of freedom of these two methods.It may get superconvergence result on this triangular mesh.AcknowledgmentsThe authors would like to thank the editor and reviewers for their valuable comments and suggestions which helped us improve the results of this paper.This work is in parts supported by the NSF of China no.10901131,no.10971166 ,the National High TechnologyX Y0.20.40.60.8100.20.40.60.81aY0.20.40.60.81X 00.20.40.60.8154321−1−2−3−4−50.113643−0.1104 bFigure 8:Velocity streamlines a and pressure level lines b for the local Gauss integration method.X Y0.20.40.60.8100.20.40.60.81 aY0.20.40.60.81X0.20.40.60.8154321−1−2−3−4−50.108278−0.106764 bFigure 9:Velocity streamlines a and pressure level lines b for the local Gauss integration method with the nonconforming element.Research and Development Program of China 863Program,no.2009AA01A135 ,the China Postdoctoral Science Foundation no.200801448,no.20070421155 ,and the Natural Science Foundation of Xinjiang Province no.2010211B04 .References1 I.Babu ˇs ka and J.E.Osborn,“Eigenvalue problems,”in Handbook of Numerical Analysis ,P .G.Ciarletand J.L.Lions,Eds.,vol.II,pp.641–787,Finite Element Method,North-Holland,Amsterdam,1991. 2 I.Babu ˇs ka and J.E.Osborn,“Finite element-Galerkin approximation of the eigenvalues andeigenvectors of selfadjoint problems,”Mathematics of Computation ,vol.52,no.186,pp.275–297,1989. 3 Q.Lin and H.Xie,“Asymptotic error expansion and Richardson extrapolation of eigenvalueapproximations for second order elliptic problems by the mixed finite element method,”Applied Numerical Mathematics ,vol.59,no.8,pp.1884–1893,2009.4 Q.Lin,“Fourth order eigenvalue approximation by extrapolation on domains with reentrantcorners,”Numerische Mathematik,vol.58,no.6,pp.631–640,1991.5 S.Jia,H.Xie,X.Yin,and S.Gao,“Approximation and eigenvalue extrapolation of Stokes eigenvalueproblem by nonconformingfinite element methods,”Applications of Mathematics,vol.54,no.1,pp.1–15,2009.6 H.Chen,S.Jia,and H.Xie,“Postprocessing and higher order convergence for the mixedfinite elementapproximations of the eigenvalue problem,”Applied Numerical Mathematics,vol.61,no.4,pp.615–629, 2011.7 H.Chen,S.Jia,and H.Xie,“Postprocessing and higher order convergence for the mixedfinite elementapproximations of the Stokes eigenvalue problems,”Applications of Mathematics,vol.54,no.3,pp.237–250,2009.8 C.Lovadina,M.Lyly,and R.Stenberg,“A posteriori estimates for the Stokes eigenvalue problem,”Numerical Methods for Partial Differential Equations,vol.25,no.1,pp.244–257,2009.9 K.A.Cliffe,E.Hall,and P.Houston,“Adaptive discontinuous Galerkin methods for eigenvalueproblems arising in incompressiblefluidflows,”SIAM Journal on Scientific Computing,vol.31,no.6, pp.4607–4632,2010.10 X.Yin,H.Xie,S.Jia,and S.Gao,“Asymptotic expansions and extrapolations of eigenvalues for theStokes problem by mixedfinite element methods,”Journal of Computational and Applied Mathematics, vol.215,no.1,pp.127–141,2008.11 W.Chen and Q.Lin,“Approximation of an eigenvalue problem associated with the Stokes problemby the stream function-vorticity-pressure method,”Applications of Mathematics,vol.51,no.1,pp.73–88,2006.12 B.Mercier,J.Osborn,J.Rappaz,and P.A.Raviart,“Eigenvalue approximation by mixed and hybridmethods,”Mathematics of Computation,vol.36,no.154,pp.427–453,1981.13 J.C.Xu and A.H.Zhou,“A two-grid discretization scheme for eigenvalue problems,”Mathematics ofComputation,vol.70,no.233,pp.17–25,2009.14 Z.X.Chen,Finite Element Methods and Their Applications,Springer,Heidelberg,Germany,2005.15 V.Girault and P.A.Raviart,Finite Element Method for Navier-Stokes Equations:Theory and Algorithms,Springer,Berlin,Germany,1987.16 B.Brefort,J.M.Ghidaglia,and R.Temam,“Attractors for the penalized Navier-Stokes equations,”SIAM Journal on Mathematical Analysis,vol.19,no.1,pp.1–21,2001.17 Y.N.He,“Optimal error estimate of the penaltyfinite element method for the time-dependent Navier-Stokes equations,”Mathematics of Computation,vol.74,no.251,pp.1201–1216,2005.18 J.Li,L.Q.Mei,and Y.N.He,“A pressure-Poisson stabilizedfinite element method for thenon-stationary Stokes equations to circumvent the inf-sup condition,”Applied Mathematics and Computation,vol.182,no.1,pp.24–35,2006.19 J.Douglas Jr.and J.Wang,“An absolutely stabilizedfinite element method for the Stokes problem,”Mathematics of Computation,vol.52,no.186,pp.495–508,1989.20 R.Araya,G.R.Barrenechea,and F.Valentin,“Stabilizedfinite element methods based on multiscaleenrichment for the stokes problem,”SIAM Journal on Numerical Analysis,vol.44,no.1,pp.322–348, 2006.21 J.Li,Y.N.He,and Z.X.Chen,“A new stabilizedfinite element method for the transient Navier-Stokesequations,”Computer Methods in Applied Mechanics and Engineering,vol.197,no.1–4,pp.22–35,2007.22 J.Li and Y.N.He,“A stabilizedfinite element method based on two local Gauss integrations for theStokes equations,”Journal of Computational and Applied Mathematics,vol.214,no.1,pp.58–65,2008.23 X.L.Feng,I.Kim,H.Nam,and D.Sheen,“Locally stabilized P1-nonconforming quadrilateral andhexahedralfinite element methods for the Stokes equations,”Journal of Computational and Applied Mathematics.In press.24 J.Li,Y.N.He,and Z.X.Chen,“Performance of several stabilizedfinite element methods for theStokes equations based on the lowest equal-order pairs,”Computing,vol.86,no.1,pp.37–51,2009.25 F.Brezzi and M.Fortin,Mixed and Hybrid Finite Element Methods,vol.15,Springer,New York,NY,USA,1991,Springer Series in Computational Mathematics.26 J.Li and Z.X.Chen,“A new local stabilized nonconformingfinite element method for the Stokesequations,”Computing,vol.82,no.2-3,pp.157–170,2008.27 FreeFem ,version2.23,/.。