Quasi-Explicit Calibration of Gatheral's SVI model

地震勘探发展史第一篇：地震勘探发展史地震勘探发展史利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理勘探方法叫作地震勘探。

地震勘探是钻探前勘测石油与天然气资源的重要手段。

地震勘探起始于19世纪中叶1845年，R.马利特曾用人工激发的地震波来测量弹性波在地壳中的传播速度。

1913年前后R.费森登发明反射法地震勘探。

1921年，J.C.卡彻将反射法地震勘探投入实际应用。

1930年，通过反射法地震勘探工作,在该地区发现了3个油田。

从此，反射法进入了工业应用的阶段。

20世纪早期德国L.明特罗普发现折射法地震勘探。

20世纪30年代，苏联Г。

А。

甘布尔采夫等吸收了反射法的记录技术，对折射法作了相应的改进。

20世纪50～60年代，反射法的光点照相记录方式被模拟磁带记录方式所代替，从而可选用不同因素进行多次回放，提高了记录质量。

20世纪70年代，模拟磁带记录又为数字磁带记录所取代，形成了以高速数字计算机为基础的数字记录、多次覆盖技术、地震数据处理技术相互结合的完整技术系统，大大提高了记录精度和解决地质问题的能力。

从20世纪70年代初期开始，采用地震勘探方法研究岩性和岩石孔隙所含流体成分。

我国的地震勘探发展1955年，我国煤炭工业上开始采用地震勘探技术，并在华东组建了全国第一支地震勘探队伍。

1971年，由煤炭科学研究总院西安分院、渭南煤矿专用设备厂研制成功MD-1型半导体磁带记录地震仪，这是我国第一套自行设计制造的煤田地震勘探仪器，并在国内煤田地震队中推广应用。

1979年我国打破了西方国家的技术封锁，成功研制出MDS-1型数字地震仪，对数字地震勘探起到了很大的推动作用。

1984~1985年，随着对外改革开放政策的实施，我国煤田地震勘探队伍开始从国外引进21套以DFS-V和SN338为主的数字地震仪，同时引进了以IBM-4381为主机的地震数据处理系统。

1978年，中国煤田地质总局在伊敏河矿区开展煤田三维地震勘探技术前提性研究。

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。

1、microscopic world 微观世界2、macroscopic world 宏观世界3、quantum theory 量子[理]论4、quantum mechanics 量子力学5、wave mechanics 波动力学6、matrix mechanics 矩阵力学7、Planck constant 普朗克常数8、wave-particle duality 波粒二象性9、state 态10、state function 态函数11、state vector 态矢量12、superposition principle of state 态叠加原理13、orthogonal states 正交态14、antisymmetrical state 正交定理15、stationary state 对称态16、antisymmetrical state 反对称态17、stationary state 定态18、ground state 基态19、excited state 受激态20、binding state 束缚态21、unbound state 非束缚态22、degenerate state 简并态23、degenerate system 简并系24、non-deenerate state 非简并态25、non-degenerate system 非简并系26、de Broglie wave 德布罗意波27、wave function 波函数28、time-dependent wave function 含时波函数29、wave packet 波包30、probability 几率31、probability amplitude 几率幅32、probability density 几率密度33、quantum ensemble 量子系综34、wave equation 波动方程35、Schrodinger equation 薛定谔方程36、Potential well 势阱37、Potential barrien 势垒38、potential barrier penetration 势垒贯穿39、tunnel effect 隧道效应40、linear harmonic oscillator线性谐振子41、zero proint energy 零点能42、central field 辏力场43、Coulomb field 库仑场44、δ-function δ-函数45、operator 算符46、commuting operators 对易算符47、anticommuting operators 反对易算符48、complex conjugate operator 复共轭算符49、Hermitian conjugate operator 厄米共轭算符50、Hermitian operator 厄米算符51、momentum operator 动量算符52、energy operator 能量算符53、Hamiltonian operator 哈密顿算符54、angular momentum operator 角动量算符55、spin operator 自旋算符56、eigen value 本征值57、secular equation 久期方程58、observable 可观察量59、orthogonality 正交性60、completeness 完全性61、closure property 封闭性62、normalization 归一化63、orthonormalized functions 正交归一化函数64、quantum number 量子数65、principal quantum number 主量子数66、radial quantum number 径向量子数67、angular quantum number 角量子数68、magnetic quantum number 磁量子数69、uncertainty relation 测不准关系70、principle of complementarity 并协原理71、quantum Poisson bracket 量子泊松括号72、representation 表象73、coordinate representation 坐标表象74、momentum representation 动量表象75、energy representation 能量表象76、Schrodinger representation 薛定谔表象77、Heisenberg representation 海森伯表象78、interaction representation 相互作用表象79、occupation number representation 粒子数表象80、Dirac symbol 狄拉克符号81、ket vector 右矢量82、bra vector 左矢量83、basis vector 基矢量84、basis ket 基右矢85、basis bra 基左矢86、orthogonal kets 正交右矢87、orthogonal bras 正交左矢88、symmetrical kets 对称右矢89、antisymmetrical kets 反对称右矢90、Hilbert space 希耳伯空间91、perturbation theory 微扰理论92、stationary perturbation theory 定态微扰论93、time-dependent perturbation theory 含时微扰论94、Wentzel-Kramers-Brillouin method W. K. B.近似法95、elastic scattering 弹性散射96、inelastic scattering 非弹性散射97、scattering cross-section 散射截面98、partial wave method 分波法99、Born approximation 玻恩近似法100、centre-of-mass coordinates 质心坐标系101、laboratory coordinates 实验室坐标系102、transition 跃迁103、dipole transition 偶极子跃迁104、selection rule 选择定则105、spin 自旋106、electron spin 电子自旋107、spin quantum number 自旋量子数108、spin wave function 自旋波函数109、coupling 耦合110、vector-coupling coefficient 矢量耦合系数111、many-partic le system 多子体系112、exchange forece 交换力113、exchange energy 交换能114、Heitler-London approximation 海特勒-伦敦近似法115、Hartree-Fock equation 哈特里-福克方程116、self-consistent field 自洽场117、Thomas-Fermi equation 托马斯-费米方程118、second quantization 二次量子化119、identical particles全同粒子120、Pauli matrices 泡利矩阵121、Pauli equation 泡利方程122、Pauli’s exclusion principle泡利不相容原理123、Relativistic wave equation 相对论性波动方程124、Klein-Gordon equation 克莱因-戈登方程125、Dirac equation 狄拉克方程126、Dirac hole theory 狄拉克空穴理论127、negative energy state 负能态128、negative probability 负几率129、microscopic causality 微观因果性本征矢量eigenvector本征态eigenstate本征值eigenvalue本征值方程eigenvalue equation本征子空间eigensubspace (可以理解为本征矢空间)变分法variatinial method标量scalar算符operator表象representation表象变换transformation of representation表象理论theory of representation波函数wave function波恩近似Born approximation玻色子boson费米子fermion不确定关系uncertainty relation狄拉克方程Dirac equation狄拉克记号Dirac symbol定态stationary state定态微扰法time-independent perturbation定态薛定谔方程time-independent Schro(此处上面有两点)dinger equation 动量表象momentum representation角动量表象angular mommentum representation占有数表象occupation number representation坐标（位置）表象position representation角动量算符angular mommentum operator角动量耦合coupling of angular mommentum对称性symmetry对易关系commutator厄米算符hermitian operator厄米多项式Hermite polynomial分量component光的发射emission of light光的吸收absorption of light受激发射excited emission自发发射spontaneous emission轨道角动量orbital angular momentum自旋角动量spin angular momentum轨道磁矩orbital magnetic moment归一化normalization哈密顿hamiltonion黑体辐射black body radiation康普顿散射Compton scattering基矢basis vector基态ground state基右矢basis ket ‘右矢’ket基左矢basis bra简并度degenerancy精细结构fine structure径向方程radial equation久期方程secular equation量子化quantization矩阵matrix模module模方square of module内积inner product逆算符inverse operator欧拉角Eular angles泡利矩阵Pauli matrix平均值expectation value （期望值）泡利不相容原理Pauli exclusion principle氢原子hydrogen atom球鞋函数spherical harmonics全同粒子identical partic les塞曼效应Zeeman effect上升下降算符raising and lowering operator 消灭算符destruction operator产生算符creation operator矢量空间vector space守恒定律conservation law守恒量conservation quantity投影projection投影算符projection operator微扰法pertubation method希尔伯特空间Hilbert space线性算符linear operator线性无关linear independence谐振子harmonic oscillator选择定则selection rule幺正变换unitary transformation幺正算符unitary operator宇称parity跃迁transition运动方程equation of motion正交归一性orthonormalization正交性orthogonality转动rotation自旋磁矩spin magnetic monent(以上是量子力学中的主要英语词汇，有些未涉及到的可以自由组合。

高分子物理常见名词Θ溶剂(Θ solvent)：链段-溶剂相互吸引刚好抵消链段间空间排斥的溶剂，形成高分子溶液时观察不到远程作用，该溶剂中的高分子链的行为同无扰链Θ温度(Θ temperature)：溶剂表现出Θ溶剂性质的温度Argon理论(Argon theory)：一种银纹扩展过程的模型，描述了分子链被伸展将聚合物材料空化的过程Avrami方程(Avrami equation)：描述物质结晶转化率与时间关系的方程：Kt-α，α为转化率，K与n称Avrami常数(Avrami constants) =-1n)exp(Bingham流体(Bingham liquid)：此类流体具有一个屈服应力σy，应力低于σy时不产生形变，当应力大于σy时才发生流动，应力高于σy的部分与应变速率呈线性关系Boltzmann叠加原理(Blotzmann superposition principle)：Boltzmann提出的粘弹性原理：认为样品在不同时刻对应力或应变的响应各自独立并可线性叠加Bravais晶格(Bravais lattice)：结构单元在空间的排列方式Burger's模型(Burger's model)：由一个Maxwell模型和一个Kelvin模型串联构成的粘弹性模型Cauchy应变(Cauchy strain)：拉伸引起的相对于样品初始长度的形变分数，又称工程应变Charpy冲击测试(Charpy impact test)：样品以简支梁形式放置的冲击强度测试，测量样品单位截面积的冲击能Considère构图(Considère construction)：以真应力对工程应作图以判定细颈稳定性的方法Eyring模型(Eyring model)：一种描述材料形变过程的分子模型，认为形变是结构单元越过能垒的跳跃式运动Flory-Huggins参数(Flory-Huggins interaction parameter)：描述聚合物链段与溶剂分子间相互作用的参数，常用χ表示，物理意义为一个溶质分子被放入溶剂中作用能变化与动能之比2.11.2Flory构图(Flory construction)：保持固定拉伸比所需的力f对实验温度作图得到，由截距确定内能对拉伸力的贡献，由斜率确定熵对拉伸力的贡献Flory特征比(characteristic ratio)：无扰链均方末端距与自由连接链均方末端距的比值Griffith理论(Griffith theory)：一种描述材料断裂机理的理论，认为断裂是吸收外界能量产生新表面的过程Hencky应变(Hencky strain)：拉伸引起的相对于样品形变分数积分，又称真应变Hermans取向因子(Hermans orientation factor)：描述结构单元取向程度的参数，是结构单元与参考方向夹角余弦均方值的函数Hoffman-Weeks作图法(Hoffman-Weeks plot)：一种确定平衡熔点的方法。

Numerical simulation of quench propagation at early phaseby high-order methodsShaolin Maoa,*,Cesar A.Luongo a ,David A.KoprivabaCenter for Advanced Power Systems,and National High Magnetic Field Laboratory,Florida State University,Tallahassee,FL 32310,United StatesbDepartment of Mathematics,Florida State University,Tallahassee,FL 32304,United StatesAbstractIn this paper a high-order discontinuous Galerkin (DG)spectral element method (SEM)is introduced to deal with thermo-hydraulic quench simulation in superconducting magnets,specifically in the case of cable-in-conduit conductors (CICC).We will focus on 1D quench propagation at early phase which is important to understand the mechanism of quench propagation and protect the supercon-ducting magnets.A second-order polynomial interpolation to the helium heat capacity is introduced in this study,and boundary con-ditions are imposed along the characteristics to overcome the typical difficulties encountered in numerical studies.Explicit fourth order Runge–Kutta time integration is used to match the high accuracy in space.Numerical results demonstrate the advantages in accuracy of this method when used to simulate early quench development.Ó2006Published by Elsevier Ltd.Keywords:Quench propagation;Superconducting magnets;Discontinuous Galerkin;Spectral element1.IntroductionIn large-scale superconducting magnet systems,thermal stability is a key issue and quench propagation is always an important consideration.Superconductors are designed to operate at very high current density so that when an external perturbation is strong enough,the superconduc-tor will go from the superconducting state to the normal state (a resistive conductor).This transition constitutes a quench.When the magnet quenches,the normal zone evolves and expands with time.The temperature gradient in the normal zone spans a range from liquid helium tem-perature in the undisturbed conductor,to a maximum tem-perature at the center of the normal front (where the original perturbation took place)[1].We focus on cable-in-conduit conductors (CICC),which are commonly used in large-scale magnets.As the name implies,they consist of a superconducting cable inside a metal pipe that contains the liquid helium coolant.Almost invariably the liquid helium inside the conduit is supercrit-ical.The high Reynolds number,low Mach number fluid flow in the CICC is governed by the unsteady convec-tion–diffusion equations.One difficulty in simulating the flow is the complicated heat transfer coupling between the fluid,the conductors and the conduit wall.Another dif-ficulty comes from the highly non-linear physical properties of solid materials and liquid helium [2].Many numerical methods have been considered to sim-ulate quench in superconducting magnets.Examples include finite element,finite volume and finite difference methods [3–10].However,most of them are first or sec-ond-order accurate in space and time.For long time wave traveling problems,large dissipation and dispersion errors seriously degrade the solutions.The main goal of our research is to seek an algorithm with both high accuracy and efficiency to solve the quench0011-2275/$-see front matter Ó2006Published by Elsevier Ltd.doi:10.1016/j.cryogenics.2006.01.006*Corresponding author.Present address:Department of Mechanical Engineering and National Institute for Global Environmental Change,Southcentral Regional Center,Tulane University,New Orleans,LA 70118,United States./locate/cryogenicsCryogenics 46(2006)589–596propagation problem.In this paper,we focus on a1D dis-continuous Galerkin(DG)spectral element method(SEM) with explicit Runge–Kutta time integration[10,11].The SEM is a weighted-residual technique for the solution of partial differential equations that combines the geometric flexibility of low-orderfinite element methods with the rapid convergence rate of spectral methods[12].Roe’s approximate Riemann solver with modifications for real gas effects is used to treat the numericalflux generated in the weak form of the governing system of equations[13].2.Modeling of quench propagation in CICCA typical CICC has a large ratio of length-to-diameter, up to105,so the compressible heliumflow can be simplified by considering it only along the longitudinal direction of the channel[1].Fig.1is a schematic of the physical system under study.The cable of the CICC can be assumed to be at uniform temperature across the conductor.The helium exchanges heat between the conductors and the conduit wall.Friction also plays an important role in the helium flow.The basic governing equations consist of the continu-ity of mass,momentum and energy for the helium,and the coupled heat balance equations for the conductors and the conduit.The system of equations is written in the conserva-tive form:oo tQþoo xF¼S;ð1Þq st C sto T sto t¼oo xk sto T sto xþs3;ð2Þq jk C jko T jko t¼s4;ð3ÞwhereQ¼qmeB@1C A;F¼mpþm2=qðeþpÞm=qB@1C A;S¼s1s2B@1C A;m¼q v;e¼qðeþv2=2Þ¼q E;S1¼Àf q v j v j2d h;s2¼h st P stA HeðT stÀTÞþh jk P jkA HeðT jkÀTÞ;s3¼Àh st P ststðT stÀTÞÀh st–jk P ststðT stÀT jkÞþJðx;t;T stÞþa; s4¼Àh jk P jkA jkðT jkÀTÞÀh st–jk P st–jkA jkðT jkÀT stÞ;NomenclatureA Cu area of copper in the conductor,m2A st area of the conductor,m2A He area of helium,m2c specific heat of helium,J/kg Kc st specific heat of conductor,J/kg Kc jk specific heat of conduit,J/kg Kd h hydraulic diameter,me specific internal energy,m2/s2f Darcy friction factor for CICCh st heat transfer coefficient(helium–conductor),W/K m2h jk heat transfer coefficient(helium–conduit),W/K m2h st–jk heat transfer coefficient(conductor–conduit),W/K m2I operating current,A J Joule heating,W/m3k st thermal conductivity of conductor,W/K m m q vp helium pressure,PaP st wet perimeter for conductor,mP jk wet perimeter for conduit,mP st–jk wet perimeter for conductor–conduit,mT liquid helium temperature,KT st conductor temperature,KT jk conduit temperature,Kv helium velocity,m/sa external heat sources,W/m3q density of helium,kg/m3q st density of helium,kg/m3q jk density of the conductor,kg/m3q e resistivity of copper,Xm Fig.1.Schematic of one-dimensional CICC,all strands are combined as acomponent,and there is a thermal coupling between each two componentsof conduit(jacket),strands,and liquid helium(fluid).590S.Mao et al./Cryogenics46(2006)589–596and J (x ,t ,T st )=q e I 2/(A Cu A st )is the Joule heating in the conductors.The equation of state of helium is added to close the system:p ¼p ðq ;T Þ:ð4ÞEqs.(1)–(4),along with well-defined boundary conditionsand initial conditions,define the governing system of equa-tions for the quasi-1D quench propagation problem (see Nomenclature for a definition of symbols).The left hand side of Eq.(1)is the Euler equations,so it is reasonable to analyze helium flow by characteristic meth-ods.Riemann variables are required to evaluate the bound-ary fluxes in the DG spectral element computation.To solve the Riemann integral with global second-order accu-racy,an approximating polynomial based on a carefully chosen equilibrium state (p 0,T 0)is used at each time step [14,15]given by 1q c¼1q cþa 1ðp Àp 0Þþa 2ðp Àp 0Þ2þa 3ðT ÀT 0Þþa 4ðT ÀT 0Þ2þa 5ðp Àp 0ÞðT ÀT 0Þð5Þand illustrated in Fig.2.It is very important to note that the boundary condi-tions should be imposed carefully to keep the numerical computation stable.The governing Eq.(1)are hyperbolic,so only conditions imposed along the characteristic will not reflect spurious waves into the computation domain.At the later phase of the quench propagation (after 0.01–0.1s),more disturbances (Joule heating and friction forces)are added to the system and parabolic behavior dominates the physical evolution of quench.Therefore,it is reason-able to replace non-reflecting boundary condition by resorting to some empirical process [7].3.Discontinuous Galerkin spectral element methods Spectral element methods can be regarded as high-order finite element methods.They incorporate the flexibility of the finite element simulation and rapid convergence of spectral discretization.Spectral methods have the advan-tage of small phase errors for long time evolution of quench propagation.Spectral element methods can also be parallelized.Due to the large non-linear source terms in the governing equations (1)–(3),there is a prohibiting restriction on explicit time integration which also make DG-spectral element methods a reasonable choice to get high accuracy and high efficiency by using parallel comput-ing techniques.We discretize the system of Eqs.(1)–(3)by a discontin-uous Galerkin spectral element method (DG-SEM)[11].Advantages of the approach include ease of changing the approximation order,exponentially small dispersion and dissipation errors [16],mesh flexibility,and efficient parallel implementation [17]at any approximation order.In 1D,the region under consideration is divided into non-overlapping line segments,which correspond to the elements.The space coordinate within each element is then mapped onto the interval n 2[À1,1]by an affine transformation:x ¼x k À1þ1þn 2ðx k Àx k À1Þ¼x k À1þ1þn 2d k:ð6ÞWithin each element the solution and the fluxes are approx-imated by N th-order polynomials defined at the Legendre–Gauss quadrature points [12].That is F ¼X N i ¼0F i ‘i ðn Þ;Q ¼X N i ¼0Q i ‘i ðn Þ;S ¼X N i ¼0S i ‘i ðn Þ;where‘i ðn Þ¼Y N j ¼1;j ¼i n Àn jn iÀn j :After some algebraic manipulation we have the finalapproximation of 1D quench propagation equations in the collocation (nodal)form:d Q i d t þF Ãð1Þ‘i ð1Þw i ÀF ÃðÀ1Þ‘i ðÀ1Þw i ÀX jF j ð‘0i ;‘j ÞNw i "#¼S i ;ð7Þwhere the w i are the Gauss quadrature weights.A DG spectral element discretization can also be applied to the heat balance equations by adding an equation to rep-resent the diffusive flux,viz.,h Ào ðc T st Þ¼0;Fig.2.Interpolation errors to heat capacity q c of supercritical helium by using Eq.(5)to fit the data from HEPAK software package [14].The Riemann integral is not sensitive to the second-order fitting curve methods used to deal with supercritical helium.S.Mao et al./Cryogenics 46(2006)589–596591where c¼k stq st c st is the equivalent diffusion coefficient of theconductor.The spectral discretization of the heat balance is given byh i þ"ðc T stÞjn¼1‘ið1Þw iÀðc T stÞjn¼À1‘iðÀ1Þw iþXjðc T stÞjð‘0i;‘jÞNw i#¼0:ð8ÞRoe’s approximate Riemann solver,modified for the real gas effects,is chosen to compute the numericalflux at the interface[12].The use of the Riemann solver at the element faces makes the imposition of boundary conditions simple, since one needs only to specify the external state as the in-put for the Riemann solver[13].The semi-discrete approximations(7)and(8)are inte-grated by the Runge–Kutta method.The coefficients for the fourth-order Runge–Kutta method can be found in[18].For unsteady problems,the time step is restricted bya Courant–Friedrichs–Lewy(CFL)condition for the approximation of the Euler equations,and by the sec-ond-order derivative terms in the heat balance equations for the conductor and conduit.4.Numerical results and discussionsThe method is tested on two problems,both for a short length coil of NbTi superconductor.Thefirst has an exter-nal heat pulse added to the conduit to put the system over the stability margin.The numerical results will be com-pared to those obtained by the commercialfinite element code,GANDALF[19].The second case is an experiment by Ando et al.[20]will be used to benchmark our numer-ical simulation.Thefirst benchmark problem that we solve was reported by Arp[3],who solved three cases:One case with an initial perturbation that gives an unquenchedflow(case(a)),and two cases in which the initial heat pulses are large enough to lead to quench(cases(b)and(c)).In these problems,the initial helium temperature and pressure were4.5K and 3·105Pa.The operating current and magneticfield are 104A and6.0T.We simulate theflow in a20m CICC cable,with heat pulses up to60J/kg of conductor applied to a1m central section of the dummy coil.Two different external heat pulses were chosen to test the DG-spectral element method.Input data are the same except that explicit time integration was used for DG-SEM and implicit methods were used in GANDALF.To Fig.3.Evolution and distribution of quenching pressure in a short CICCsample[3]using(a)DG-SEM and(b)GANDALF code,the external heatpulse was8·103W/m in a1m zone for0.001s.Fig.4.Conductor temperature distribution at different time using(a)DG-SEM and(b)GANDALF code,the external heat pulse was8·103W/min a1m zone for0.001s.592S.Mao et al./Cryogenics46(2006)589–596compare the accuracy,the minimum time step is set to the same value in DG-SEM and GANDALF.The numerical results in Figs.3–5focus on the early phase of quench propagation.Clearly,a better resolution in large gradient region and strong stability are shown in the simulation by DG-SEM.When the external heat pulse is4·104W/m or higher,GANDALF becomes unstable,while the DG-SEM remains stable for this high heatflux.DG-SEM has the advantage of numerical conservation by using Riemann Solvers to deal with theflux at each element faces.The second benchmarking is an experiment from Ando et al.[20],in which a short NbTi superconductor coil was used.A quench was initiated with a fast heat pulse by a heater that was added for0.1ms in a4mm zone of a 26m dummy coil.The operating current was varied between1.5kA and2.0kA.The boundary conditions were an assumed subsonic velocity at the exit,and a symmetry condition at the middle point.The initial helium tempera-ture and pressure were given by 4.2K and1MPa.A Gaussian heat pulse(for numerical resolution reasons) was imposed in the middle to initiate the quench simulation.The numerical results by DG-SEM and the GANDALF code were compared to the same experimental data[20]. The GANDALF code uses linearfinite elements,for a sec-ond-order spatial approximation,and either afirst or sec-ond-order approximation in time[7].The DG-SEM used a fourth-order approximate within each element.The same mesh size was used in both simulations so the accuracy can be compared for the same number of degrees of freedom. The GANDALF code can use adaptive time stepping to improve the efficiency of the computation.Our DG-SEM only uses explicit fourth-order Runge–Kutta time integra-tion without time step control.The main reason for us to use explicit methods is that the quench propagation simu-lation is an unsteady problem with large non-linear rge time steps can speed up the simulation at the price of the accuracy of the solution.The evolution and distribution of the helium density, conductor temperature,helium quench pressure,and Fig.5.Distribution of the induced velocity of helium at different time[3]using(a)DG-SEM and(b)GANDALF code,the external heat pulse was8·103W/m in a1m zone for0.001s.Fig.6.Evolution and distribution of helium density for I=1.5kA inAndo’s case by(a)DG-SEM and(b)FEM code GANDALF.S.Mao et al./Cryogenics46(2006)589–596593helium induced-flow velocities are shown in Figs.6–9.The results shown are for an operating current I =1.5kA.The numerical results by DG-SEM exhibit a very good agreement with GANDALF code.The maximum quench pressure profile (Fig.7a and b)shows that the difference in maximum quench pressure is about 3–5%between the DG-SEM and GANDALF.Considering the different procedures in which the helium data and their derivatives are applied in two numerical methods,the numerical results by DG-SEM appear to be sufficiently accurate for engi-neering applications and in good agreement with GANDALF.The GANDALF code cannot impose the non-reflecting boundary conditions,which were used in DG-SEM simulation.The governing system of equations for helium flow is hyperbolic,and we claim that a pressure expansion procedure at the two ends of the CICC tube is much better approximation to the actual physics,than a fixed pressure assumed in GANDALF.The predicted exit velocity and density of helium is about the same for both codes.The numerical results by using GANDALF have some obvious oscillations in the large gradient regions duringquench propagation (Fig.7b)but a good resolution was obtained by using DG-SEM (Fig.7a).Discontinuous Galerkin spectral element methods are based on the conser-vative form of the governing equations and the Roe’s Riemann solver was used to keep conservation of the numerical flux across the interfaces of every element.In the GANDALF code the central discretization in space was applied and artificial damping terms are required to make it stable.For quench evolution,the choice of artifi-cial damping is a non-trivial issue.If the numerical damp-ing is too large,the numerical system may not be conservative.Nevertheless,the simulation of conductor temperature and normal zone is still accurate.All these potential problems are avoided with DG-SEM.Finally,the numerical results are benchmarked by com-paring with experimental data for the normal zone speed in Figs.10and 11.Good agreement (5%difference)is found for the normal zone evolution between the numerical simulation and the experimental data for different operat-ing currents.At the lower operating current situation (I =1.5kA)both DG-SEM and GANDALF approxi-mated the experimental results well,but DG-SEM shows better accuracy in the high operating currentcaseFig.7.Evolution and distribution of quench pressure for I =1.5kA in Ando’s case by (a)DG-SEM and (b)FEM codeGANDALF.Fig.8.Evolution and distribution of induced-flow velocity in CICC for I =1.5kA in Ando’s case by (a)DG-SEM and (b)FEM code GANDALF.594S.Mao et al./Cryogenics 46(2006)589–596(I =1.8kA),especially with increased number of mesh elements.In quench propagation simulations the computational time is always an important issue.Because two computa-tions of Roe’s approximate Riemann solver are needed in each element when using the discontinuous Galerkin spec-tral element methods (DG-SEM),the computation time by using DG-SEM is longer than that by using GANDALF code (21.5h CPU time for DG-SEM vs.18h for GAN-DALF when simulating a 3s quench propagation on our SUN UNIX workstation).However,DG-SEM shows more robust results than GANDALF in terms of keeping numerical stability for large perturbations in quench prop-agation.The efficiency of simulation can be improved by resorting to parallel computing techniques,which is an area of potential future work.5.ConclusionA discontinuous Galerkin spectral element method was successfully used to track a special moving front problem,namely quench propagation in CICC.Special effort was made to deal with the real gas/fluid properties (supercritical fluid)of helium.A simple approach was proposed in this study to approximate the Riemann integral,which is the critical step of applying Roe’s approximate Riemann solver in DG-SEM.A second-order curve-fitting was used to compute helium properties.The DG-SEM was used to overcome numerical difficulties encountered by traditional methods,such as finite difference and finite element,in which artificial damping is required.The DG-SEM is a robust method that can obtain high accuracy (resolution)in large gradient regions and demonstrate numerical stabil-ity for large external disturbances (strong heat pulses).Differences are seen between the DG-SEM results and GANDALF code and experimental data.These are largely attributed to the limitation of preconditioning of helium properties (equation of state),boundary conditions,and forcing terms in the governing equations.More important,however,is how to improve the efficiency.ParallelizationofFig.9.Evolution and distribution of conductor temperature for I =1.5kA in Ando’s case by (a)DG-SEM and (b)FEM codeGANDALF.parison of evolution of normal zone by numerical results and experimental data (I =1.5kA).parison of evolution of normal zone by numerical results and the experimental data (I =1.8kA).S.Mao et al./Cryogenics 46(2006)589–596595the discontinuous Galerkin spectral element methods for quench propagation simulations is the main consideration of our research work in the future.AcknowledgementThe research was partially supported by the Center for Advanced Power Systems,Florida State University,with funding from the Office of Naval Research(ONR). References[1]Dresner L.Stability of superconductors.New York:Plenum;1995.[2]Van Sciver SW.Helium cryogenics.New York:Plenum;1986.[3]Arp V.Stability and thermal quenches in force-cooled superconduc-ting cables.In:Proceedings of superconducting MHD magnet design conference.Cambridge:MIT;1980.p.142–57.[4]Marinucci C.A numerical model for the analysis of stability andquench characteristics of forced-flow cooled superconductors.Cryo-genics1983;23:579–86.[5]Wong R.Program CICC,flow and heat transfer in cable-in-conduitconductors equations and verification.LLNL report,May22,1989 [unpublished].[6]Luongo CA,Loyd RJ,Chen FK,Peck SD.Thermal-hydraulicsimulation of helium expulsion from a cable-in-conduit conductor.IEEE Trans Mag1989;25(2):1589–95.[7]Bottura L.A numerical model for the simulation of quench in theITER magnets.J Comput Phys1996;125:26–41.[8]Shajii A,Freidberg JP.Quench in superconducting magnets.I.Modeland numerical implementation.J Appl Phys1994;76(5):3149–58. [9]Koizumi N,Takahashi Y,Tsuji H.Numerical model using an implicitfinite difference algorithm for stability simulation of a cable-in-conduit superconductor.Cryogenics1996;36(9):649–59.[10]Mao S,Luongo CA,Kopriva DA.Discontinuous Galerkin spectralelement simulation of quench propagation in superconducting magnets.IEEE Trans Appl Supercond2005;15(2):1675–8.[11]Kopriva DA,WoodruffSL,Hussaini MY.Discontinuous spectralelement approximation of Maxwell’s equations.In:Cockburn B, Karniadakis GE,Shu CW,editors.Lecture notes in computational science and engineering,vol.11.Springer-Verlag;2000.p.355–61.[12]Canuto C,Hussaini MY,Quarteroni A,Zang T.Spectral methods influid dynamics.New York:Springer-Verlag;1987.[13]Roe PL.Approximate Riemann solvers,parameter vectors,anddifference-schemes.J Comput Phys1981;43(2):357–72.[14]Arp V,McCarty RD,Friend DG.Thermophysical properties ofHelium-4from0.8to2000K with pressure to2000MPa.NIST Technical Note1334,1998[revised].[15]Arp V.Private communication,April2004.[16]Stanescu D,Kopriva DA,Hussaini MY.Dispersion analysis fordiscontinuous spectral element methods.J Sci Comput2001;15: 149–71.[17]Hesthaven JS,Warburton T.Nodal high-order methods on unstruc-tured grids.J Comput Phys2002;181:186–221.[18]Press WH,Flannery BP,Teukolsky SA,Vetterling WT.Numericalrecipes.Cambridge:Cambridge University Press;1989.[19]GANDALF:A computer code for quench analysis of dualflowCICC’s,V2.1,CryoSoft,1999.[20]Ando T et al.Propagation velocity of the normal zone in a cable-in-conduit conductor.Adv Cryogen Eng1988;35A:701–8.596S.Mao et al./Cryogenics46(2006)589–596。