HIGH-ORDER WENO SIMULATIONS OF THREE-DIMENSIONAL RESHOCKED RICHTMYER-MESHKOV INSTABILITY TO LAT

合集下载

efficient implementation of weighted ENO scheme

JOURNAL OF COMPUTATIONAL PHYSICS126,202–228(1996)ARTICLE NO.0130Efﬁcient Implementation of Weighted ENO Schemes*G UANG-S HAN J IANG†AND C HI-W ANG S HUDi v ision of Applied Mathematics,Brown Uni v ersity,Pro v idence,Rhode Island02912Received August14,1995;revised January3,1996or perhaps with a forcing term g(u,x,t)on the right-handside.Here uϭ(u1,...,u m),fϭ(f1,...,f d),xϭ(x1,...,x d) In this paper,we further analyze,test,modify,and improve thehigh order WENO(weighted essentially non-oscillatory)ﬁnite differ-and tϾ0.ence schemes of Liu,Osher,and Chan.It was shown by Liu et al.WENO schemes are based on ENO(essentially non-that WENO schemes constructed from the r th order(in L1norm)oscillatory)schemes,which wereﬁrst introduced by ENO schemes are(rϩ1)th order accurate.We propose a new wayHarten,Osher,Engquist,and Chakravarthy[5]in the form of measuring the smoothness of a numerical solution,emulatingthe idea of minimizing the total variation of the approximation,of cell averages.The key idea of ENO schemes is to use which results in aﬁfth-order WENO scheme for the case rϭ3,the‘‘smoothest’’stencil among several candidates to ap-instead of the fourth-order with the original smoothness measure-proximate theﬂuxes at cell boundaries to a high order ment by Liu et al.Thisﬁfth-order WENO scheme is as fast as theaccuracy and at the same time to avoid spurious oscillations fourth-order WENO scheme of Liu et al.and both schemes arenear shocks.The cell-averaged version of ENO schemes about twice as fast as the fourth-order ENO schemes on vectorsupercomputers and as fast on serial and parallel computers.For involves a procedure of reconstructing point values from Euler systems of gas dynamics,we suggest computing the weights cell averages and could become complicated and costly for from pressure and entropy instead of the characteristic values to multi-dimensional ter,Shu and Osher[14,15] simplify the costly characteristic procedure.The resulting WENOdeveloped theﬂux version of ENO schemes which do not schemes are about twice as fast as the WENO schemes using thecharacteristic decompositions to compute weights and work well require such a reconstruction procedure.We will formulate for problems which do not contain strong shocks or strong reﬂected the WENO schemes based on thisﬂux version of ENO waves.We also prove that,for conservation laws with smooth solu-schemes.The WENO schemes of Liu et al.[9]are basedtions,all WENO schemes are convergent.Many numerical tests,on the cell-averaged version of ENO schemes.including the1D steady state nozzleﬂow problem and2D shockFor applications involving shocks,second-order schemes entropy wave interaction problem,are presented to demonstratethe remarkable capability of the WENO schemes,especially the are usually adequate if only relatively simple structures WENO scheme using the new smoothness measurement in resolv-are present in the smooth part of the solution(e.g.,the ing complicated shock andﬂow structures.We have also applied shock tube problem).However,if a problem contains rich Yang’s artiﬁcial compression method to the WENO schemes tostructures as well as shocks(e.g.,the shock entropy wave sharpen contact discontinuities.ᮊ1996Academic Press,Inc.interaction problem in Example4,Section8.3),high ordershock capturing schemes(order of at least three)are more1.INTRODUCTION efﬁcient than low order schemes in terms of CPU time andmemory requirements.In this paper,we further analyze,test,modify,and im-ENO schemes are uniformly high order accurate rightprove the WENO(weighted essentially non-oscillatory)up to the shock and are very robust to use.However,theyﬁnite difference schemes of Liu,Osher,and Chan[9]for also have certain drawbacks.One problem is with the freelythe approximation of hyperbolic conservation laws of adaptive stencil,which could change even by a round-offthe type perturbation near zeroes of the solution and its derivatives.Also,this free adaptation of stencils is not necessary in u tϩdiv f(u)ϭ0,(1.1)regions where the solution is smooth.Another problem isthat ENO schemes are not cost effective on vector super-computers such as CRAY C-90because the stencil-choos-*Research supported by ARO Grant DAAH04-94-G-0205,NSFGrants ECS-9214488and DMS-9500814,NASA Langley Grant NAG-ing step involves heavy usage of logical statements,which 1-1145and Contract NAS1-19480while the second author was in resi-perform poorly on such machines.Theﬁrst problem could dence at ICASE,NASA Langley Research Center,Hampton,VA23681-reduce the accuracy of ENO schemes for certain functions 0001,and AFOSR Grant95-1-0074.[12];however,this can be remedied by embedding certain †Current address:Department of Mathematics,UCLA,Los Angeles,CA90024.parameters(e.g.,threshold and biasing factor)into the2020021-9991/96$18.00Copyright©1996by Academic Press,Inc.All rights of reproduction in any form reserved.WEIGHTED ENO SCHEMES203 stencil choosing step so that the preferred linearly stable the total variation of the approximations.This new mea-surement gives the optimalﬁfth-order accurate WENO stencil is used in regions away from discontinuities.See[1,3,13].scheme when rϭ3(the smoothness measurement in[9]gives a fourth-order accurate WENO scheme for rϭ3). The WENO scheme of Liu,Osher,and Chan[9]is an-other way to overcome these drawbacks while keeping the Although the WENO schemes are faster than ENOschemes on vector supercomputers,they are only as fast robustness and high order accuracy of ENO schemes.Theidea is the following:instead of approximating the numeri-as ENO schemes on serial computers.In Section4,wepresent a simpler way of computing the weights for the calﬂux using only one of the candidate stencils,one usesa con v ex combination of all the candidate stencils.Each approximation of Euler systems of gas dynamics.The sim-pliﬁcation is aimed at reducing theﬂoating point opera-of the candidate stencils is assigned a weight which deter-mines the contribution of this stencil to theﬁnal approxi-tions in the costly but necessary characteristic procedureand is motivated by the following observation:the only mation of the numericalﬂux.The weights can be deﬁnedin such a way that in smooth regions it approaches certain nonlinearity of a WENO scheme is in the computation ofthe weights.We suggest using pressure and entropy to optimal weights to achieve a higher order of accuracy(anr th-order ENO scheme leads to a(2rϪ1)th-order WENO compute the weights,instead of the local characteristicquantities.In this way one can exploit the linearity of the scheme in the optical case),while in regions near disconti-nuities,the stencils which contain the discontinuities are rest of the scheme.The resulting WENO schemes(rϭ3)is about twice as fast as the original WENO scheme which assigned a nearly zero weight.Thus essentially non-oscilla-tory property is achieved by emulating ENO schemes uses local characteristic quantities to compute the weights(see Section7).The same idea can also be applied to the around discontinuities and a higher order of accuracy isobtained by emulating upstream central schemes with the original ENO ly,we can use the undivideddifferences of pressure and entropy to replace the local optimal weights away from the discontinuities.WENOschemes completely remove the logical statements that characteristic quantities to choose the ENO stencil.Thishas been tested numerically but the results are not included appear in the ENO stencil choosing step.As a result,theWENO schemes run at least twice as fast as ENO schemes in this paper since the main topic here is the WENOschemes.(see Section7)on vector machines(e.g.,CRAY C-90)andare not sensitive to round-off errors that arise in actual WENO schemes have the same smearing at contact dis-continuities as ENO schemes.There are mainly two tech-computation.Atkins[1]also has a version of ENO schemesusing a different weighted average of stencils.niques for sharpening the contact discontinuities for ENOschemes.One is Harten’s subcell resolution[4]and the Another advantage of WENO schemes is that itsﬂux issmoother than that of the ENO schemes.This smoothness other is Yang’s artiﬁcial compression(slope modiﬁcation)[20].Both were introduced in the cell average context. enables us to prove convergence of WENO schemes forsmooth solutions using Strang’s technique[18];see Section Later,Shu and Osher[15]translated them into the pointvalue framework.In one-dimensional problems,the sub-6.According to our numerical tests,this smoothness alsohelps the steady state calculations,see Example4in Sec-cell resolution technique works slightly better than theartiﬁcial compression method.However,for two or higher tion8.2.In[9],the order of accuracy shown in the error tables dimensional problems,the latter is found to be more effec-tive and handy to use[15].We will highlight the key proce-(Table1–5in[9])seemed to suggest that the WENOschemes of Liu et al.are more accurate than what the dures of applying the artiﬁcial compression method to theWENO schemes in Section5.truncation error analysis indicated.In Section2,we carryout a more detailed error analysis for the WENO schemes In Section8,we test the WENO schemes(both theWENO schemes of Liu et al.and the modiﬁed WENO andﬁnd that this‘‘super-convergence’’is indeed superﬁ-cial:the‘‘higher’’order is caused by larger error on the schemes)on several1D and2D model problems and com-pare them with ENO schemes to examine their capability coarser grids instead of smaller error on theﬁner grids.Our error analysis also suggests that the WENO schemes in resolving shock and complicatedﬂow structures.We conclude this paper by a brief summary in Section can be made more accurate than those in[9].Since the weight on a candidate stencil has to vary ac-9.The time discretization of WENO schemes will be imple-mented by a class of high order TVD Runge–Kutta-type cording to the relative smoothness of this stencil to theother candidate stencils,the way of evaluating the smooth-methods developed by Shu and Osher[14].To solve theordinary differential equationness of a stencil is crucial in the deﬁnition of the weight.In Section3,we introduce a new way of measuring thesmoothness of the numerical solution which is based uponminimizing the L2norm of the derivatives of the recon-dudt ϭL(u),(1.2)struction polynomials,emulating the idea of minimizing204JIANG AND SHUTABLE I where L (u )is a discretization of the spatial operator,the third-order TVD Runge–Kutta is simplyCoefﬁcients a r k ,lrk l ϭ0l ϭ1l ϭ2u (1)ϭu n ϩ⌬tL (u n )u (2)ϭ u n ϩ u (1)ϩ ⌬tL (u (1))(1.3)20Ϫ1/23/211/21/2u n ϩ1ϭ u n ϩ u (2)ϩ ⌬tL (u (2)).301/3Ϫ7/611/61Ϫ1/65/61/3Another useful,although not TVD,fourth-order Runge–21/35/6Ϫ1/6Kutta scheme isu (1)ϭu n ϩ ⌬tL (u n )f (u )ϭf ϩ(u )ϩf Ϫ(u ),(2.4)u(2)ϭu nϩ ⌬tL (u (1))(1.4)where df (u )ϩ/du Ն0and df (u )Ϫ/du Յ0.For example,u (3)ϭu n ϩ⌬tL (u (2))one can deﬁneu n ϩ1ϭ (Ϫu n ϩu (1)ϩ2u (2)ϩu (3))ϩ ⌬tL (u (3)).f Ϯ(u )ϭ (f (u )ϮͰu ),(2.5)This fourth-order Runge–Kutta scheme can be made TVD by an increase of operation counts [14].We will mainly where Ͱϭmax ͉f Ј(u )͉and the maximum is taken over the use these two Runge–Kutta schemes in our numerical tests whole relevant range of u .This is the global Lax–Friedrichs in Section 8.The third-order TVD Runge–Kutta scheme (LF)ﬂux splitting.For other ﬂux splitting,especially the will be referred to as ‘‘RK-3’’while the fourth-order (non-Roe ﬂux splitting with entropy ﬁx (RF);see [15]for details.TVD)Runge–Kutta scheme will be referred to as ‘‘RK-4.’’Let f ˆϩj ϩ1/2and f Ϫj ϩ1/2be,resp.the numerical ﬂuxes obtained from the positive and negative parts of f (u ),we then have2.THE WENO SCHEMES OF LIU,OSHER,AND CHANf ˆj ϩ1/2ϭf ˆϩj ϩ1/2ϩf ˆϪj ϩ1/2.(2.6)In this section,we use the ﬂux version of ENO schemes Here we will only describe how f ˆϩj ϩ1/2is computed in [9]as our basis to formulate WENO schemes of Liu et al.and on the basis of ﬂux version of ENO schemes.For simplicity,analyze their accuracy in a different way from that used we will drop the ‘‘ϩ’’sign in the superscript.The formulas in [9].We use one-dimensional scalar conservation laws for the negative part of the split ﬂux are symmetric (with (i.e.,d ϭm ϭ1in (1.1))as an example:respect to x j ϩ1/2)and will not be shown.As we know,the r th-order (in L 1sense)ENO scheme u t ϩf (u )x ϭ0.(2.1)chooses one ‘‘smoothest’’stencil from r candidate stencils and uses only the chosen stencil to approximate the ﬂux Let us discretize the space into uniform intervals of sizeh j ϩ1/2.Let’s denote the r candidate stencils by S k ,k ϭ0,⌬x and denote x j ϭj ⌬x .Various quantities at x j will be 1,...,r Ϫ1,whereidentiﬁed by the subscript j .The spatial operator of the WENO schemes,which approximates Ϫf (u )x at x j ,will S k ϭ(x j ϩk Ϫr ϩ1,x j ϩk Ϫr ϩ2,...,x j ϩk ).take the conservative formIf the stencil S k happens to be chosen as the ENO interpola-L ϭϪ1⌬x(f ˆj ϩ1/2Ϫf ˆj Ϫ1/2),(2.2)tion stencil,then the r th-order ENO approximation of h j ϩ1/2iswhere the numerical ﬂux f ˆj ϩ1/2approximates h j ϩ1/2ϭf ˆj ϩ1/2ϭq r k (f j ϩk Ϫr ϩ1,...,f j ϩk ),(2.7)h (x j ϩ1/2)to a high order with h (x )implicitly deﬁned by [15]wheref (u (x ))ϭ1⌬x͵x ϩ⌬x /2x Ϫ⌬x /2h (␰)d ␰.(2.3)q r k (g 0,...,g r Ϫ1)ϭ͸r Ϫ1l ϭ0ark ,l g l .(2.8)Here a r k ,l ,0Յk ,l Յr Ϫ1,are constant coefﬁcients.For We can actually assume f Ј(u )Ն0for all u in the range of our interest.For a general ﬂux,i.e.,f Ј(u )Ն͞0,one can later use,we provide these coefﬁcients for r ϭ2,3,in Table I.split it into two parts either globally or locally,WEIGHTED ENO SCHEMES205TABLE II To just use the one smoothest stencil among the r candi-dates for the approximation of h j ϩ1/2,is very desirable nearOptimal Weights C r kdiscontinuities because it prohibits the usage of informa-C r k k ϭ0k ϭ1k ϭ2tion on discontinuous stencils.However,it is not so desir-able in smooth regions because all the candidate stencils r ϭ21/32/3—carry equally smooth information and thus can be used r ϭ31/106/103/10together to give a higher order (higher than r ,the order of the base ENO scheme)approximation to the ﬂux h j ϩ1/2.In fact,one could use all the r candidate stencils,which all together contain (2r Ϫ1)grid values of f to give f ˆj ϩ1/2ϭq 2r Ϫ1r Ϫ1(f j Ϫr ϩ1,...,f j ϩr Ϫ1)(2.12)a (2r Ϫ1)th-order approximation of h j ϩ1/2:ϩ͸r Ϫ1k ϭ0(ͶkϪC r k )q rk (f j ϩk Ϫr ϩ1,...,f j ϩk ).f ˆj ϩ1/2ϭq 2r Ϫ1r Ϫ1(f j Ϫr ϩ1,...,f j ϩr Ϫ1)(2.9)Recalling (2.9),we see that,the ﬁrst term on the right-which is just the numerical ﬂux of a (2r Ϫ1)th-order up-hand side of the above equation is a (2r Ϫ1)th-orderstream central scheme.As we know,high order upstreamapproximation of h j ϩ1/2.Since ͚r Ϫ1k ϭ0C r k ϭ1,if we require central schemes (in space),combined with high order ͚r Ϫ1k ϭ0Ͷk ϭ1,the last summation term can be written asRunge–Kutta methods (in time),are stable and dissipative under appropriate CFL numbers and thus are convergent,according to Strang’s convergence theory [18]when the ͸r Ϫ1k ϭ0(ͶkϪC r k)(q r k(fj ϩk Ϫr ϩ1,...,f j ϩk )Ϫh j ϩ1/2).(2.13)solution of (1.1)is smooth (see Section 6).The above facts suggest that one could use the (2r Ϫ1)th-order upstream central scheme in smooth regions and only use the r th-Each term in the last summation can be made O (h 2r Ϫ1)iforder ENO scheme near discontinuities.As in (2.7),each of the stencils can render an approxima-Ͷk ϭC r k ϩO (hr Ϫ1)(2.14)tion of h j ϩ1/2.If the stencil is smooth,this approximation is r th-order accurate;otherwise,it is less accurate or even for k ϭ0,1,...,r Ϫ1.Here,h ϭ⌬x .Thus C r k will bearnot accurate at all if the stencil contains a discontinuity.the name of optimal weight.One could assign a weight Ͷk to each candidate stencil S k ,The question now is how to deﬁne the weight such that k ϭ0,1,...,r Ϫ1,and use these weights to combine the (2.14)is satisﬁed in smooth regions while essentially non-r different approximations to obtain the ﬁnal approxima-oscillatory property is achieved.In [9],the weight Ͷk for tion of h j ϩ1/2asstencil S k is deﬁned byͶk ϭͰk0r Ϫ1,(2.15)fˆj ϩ1/2ϭ͸r Ϫ1k ϭ0Ͷk q r k (f j ϩk Ϫr ϩ1,...,f j ϩk ),(2.10)wherewhere q r k (f j ϩk Ϫr ϩ1,...,f j ϩk )is deﬁned in (2.8).To achieveessentially non-oscillatory property,one then requires the weight to adapt to the relative smoothness of f on each Ͱk ϭC r k(␧ϩIS k )p,k ϭ0,1,...,r Ϫ1.(2.16)candidate stencil such that any discontinuous stencil is ef-fectively assigned a zero weight.In smooth regions,one Here ␧is a positive real number number which is intro-can adjust the weight distribution such that the resultingduced to avoid the denominator becoming zero (in our approximation of the ﬂux fˆj ϩ1/2is as close as possible to later tests,we will take ␧ϭ10Ϫ6.Our numerical tests that given in (2.9).indicate that the result is not sensitive to the choice of ␧,Simple algebra gives the coefﬁcients C r k such thatas long as it is in the range of 10Ϫ5to 10Ϫ7);the power p will be discussed in a moment;IS k in (2.16)is a smoothness measurement of the ﬂux function on the k th candidate q 2r Ϫ1r Ϫ1(f j Ϫr ϩ1,...,f j ϩr Ϫ1)ϭ͸r Ϫ1k ϭ0C r k q rk(fj ϩk Ϫr ϩ1,...,f j ϩk )(2.11)stencil.It is easy to see that ͚r Ϫ1k ϭ0Ͷk ϭ1.To satisfy (2.14),it sufﬁces to have (through a Taylor expansion analysis)and ͚r Ϫ1k ϭ0C r k ϭ1for all r Ն2.For r ϭ2,3,these coefﬁcients are given in Table II.Comparing (2.11)with (2.10),we getIS k ϭD (1ϩO (h r Ϫ1))(2.17)206JIANG AND SHUfor k ϭ0,1,...,r Ϫ1,where D is some non-zero quantity IS 1ϭ ((f Јh Ϫ f Љh 2)2ϩ(f Јh ϩ f Љh 2)2)independent of k .As we know,an ENO scheme chooses the ‘‘smoothest’’ϩ(f Љh 2)2ϩO (h 5)(2.24)ENO stencil by comparing a hierarchy of undivided differ-IS 2ϭ ((f Јh ϩ f Љh 2)2ϩ(f Јh ϩ f Љh 2)2)ences.This is because these undivided differences can be used to measure the smoothness of the numerical ﬂux on ϩ(f Љh 2)2ϩO (h 5).(2.25)a stencil.In [9],IS k is deﬁned asWe can see that the second-order terms are different from stencil to stencil.Thus (2.22)is no longer valid at critical IS k ϭ͸r Ϫ1l ϭ1͸r Ϫli ϭ1(f [j ϩk ϩi Ϫr ,l ])2r Ϫl,(2.18)points of f (u (x ))which implies that the WENO scheme of Liu et al.for r ϭ3is only third-order accurate at these points.In fact,the weights computed from the smoothness where f [и,и]is the l th undivided difference:measurement (2.18)diverge far away from the optimal weights near critical points (see Fig.1in the next section)f [j ,0]ϭf jon coarse grids (10to 80grid points per wave).But on f [j ,l ]ϭf [j ϩ1,l Ϫ1]Ϫf [j ,l Ϫ1],k ϭ1,...,r Ϫ1.ﬁne grids,since the smoothness measurements IS k for all k are relatively smaller than the non-zero constant ␧in For example,when r ϭ2,we have(2.16),the weights become close to the optimal weights.Therefore the ‘‘super-convergence’’phenomena appeared IS k ϭ(f [j ϩk Ϫ1,1])2,k ϭ0,1.(2.19)in Tables 1–5in [9]are caused by large error commitment on coarse grids and less error commitment on ﬁner grids When r ϭ3,(2.18)giveswhen using the errors of the ﬁfth-order central scheme as reference (see Tables III and IV).IS k ϭ ((f [j ϩk Ϫ2,1])2ϩ(f [j ϩk Ϫ1,1])2)(2.20)At discontinuities,it is typical that one or more of the r candidate stencils reside in smooth regions of the numerical ϩ(f [j ϩk Ϫ2,2])2,k ϭ0,1,2.solution while other stencils contain the discontinuities.The size of the discontinuities is always O (1)and does not In smooth regions,Taylor expansion analysis of (2.18)change when the grid is reﬁned.So we have for a smooth givesstencil S k ,IS k ϭ(f Јh )2(1ϩO (h )),k ϭ0,1,...,r Ϫ1,(2.21)IS k ϭO (h 2p )(2.26)where f Јϭf Ј(u j ).Note the O (h )term is not O (h r Ϫ1)that we and for a non-smooth stencil S l ,would want to have (see (2.17)).Thus in smooth monotone regions,i.e.,f Ј϶0,we haveIS l ϭO (1).(2.27)Ͷk ϭC r k ϩO (h ),k ϭ0,1,...,r Ϫ1.(2.22)From the deﬁnition of the weights (2.15),we can see that,for this non-smooth stencil S l ,the corresponding weight Recalling (2.12),we see that the WENO schemes with theͶl satisﬁessmoothness measurement given by (2.18)is (r ϩ1)th-order accurate in smooth monotone regions of f (u (x )).This re-Ͷl ϭO (h 2p ).(2.28)sult was proven in [9]using a different approach.For r ϭ2,this is optimal in the sense that the third-order upstream Therefore for small h and any positive integer power p ,central scheme is approximated in most smooth regions.the weight assigned to the non-smooth stencil vanishes as However,this is not optimal for r ϭ3,for which this h Ǟ0.Note,if there is more than one smooth stencils in measurement can only give fourth-order accuracy while the r candidates,from the deﬁnition of the weights in (2.15),the optimal upstream central scheme is ﬁfth-order accu-we expect each of the smooth stencils will get a weight rate.We will introduce a new measurement in the next which is O (1).In this case,the weights do not exactly section which will result in an optimal order accurate resemble the ‘‘ENO digital weights.’’However,if a stencil WENO scheme for the r ϭ3case.is smooth,the information that it contains is useful and When r ϭ3,Taylor expansion of (2.20)givesshould be utilized.In fact,in our extensive numerical ex-periments,we ﬁnd the WENO schemes in [9]work very IS 0ϭ ((f Јh Ϫ f Љh 2)2ϩ(f Јh Ϫ f Љh 2)2)well at shocks.We also ﬁnd that p ϭ2is adequate to obtain essentially non-oscillatory approximations at leastϩ(f Љh 2)2ϩO (h 5)(2.23)WEIGHTED ENO SCHEMES207for r ϭ2,3,although it is suggested in [9]that p should IS 1ϭ (f j Ϫ1Ϫ2f j ϩf j ϩ1)2ϩ (f j Ϫ1Ϫf j ϩ1)2(3.3)be taken as r ,the order of the base ENO schemes.We will use p ϭ2for all our numerical tests.Notice that,IS 2ϭ(f j Ϫ2f j ϩ1ϩf j ϩ2)2ϩ (3f j Ϫ4f j ϩ1ϩf j ϩ2)2.(3.4)in discussing accuracy near discontinuities,we are simply concerned with spatial approximation error.The error due In smooth regions.Taylor expansion of (3.2)–(3.4)gives,to time evolution is not considered.respectively,In summary,WENO schemes of Liu et al.deﬁned by (2.10),(2.15),and (2.18)have the following properties:IS 0ϭ (f Љh 2)2ϩ (2f Јh Ϫ f ٞh 3)2ϩO (h 6)(3.5)1.They involve no logical statements which appear in IS 1ϭ(f Љh 2)2ϩ (2f Јh ϩ f ٞh 3)2ϩO (h 6)(3.6)the base ENO schemes.IS 2ϭ (f Љh 2)2ϩ (2f Јh Ϫ f ٞh 3)2ϩO (h 6),(3.7)2.The WENO scheme based on the r th-order ENO scheme is (r ϩ1)th-order accurate in smooth monotone where f ٞϭf ٞ(u j ).If f Ј϶0,thenregions,although this is still not as good as the optimal order (2r Ϫ1).IS k ϭ(f Јh )2(1ϩO (h 2)),k ϭ0,1,2,(3.8)3.They achieve essentially non-oscillatory property by emulating ENO schemes at discontinuities.which means the weights (see (2.15))resulting from thismeasurement satisfy (2.17)for r ϭ3;thus we obtain a 4.They are smooth in the sense that the numerical ﬂux ﬁfth-order (the optimal order for r ϭ3)accurate f ˆj ϩ1/2is a smooth function of all its arguments (for a general WENO scheme.ﬂux,this is also true if a smooth ﬂux splitting method is Moreover,this measurement is also more accurate at used,e.g.,global Lax–Friedrichs ﬂux splitting).critical points of f (u (x )).When f Јϭ0,we have3.A NEW SMOOTHNESS MEASUREMENTIS k ϭ (f Љh 2)2(1ϩO (h 2)),k ϭ0,1,2,(3.9)In this section,we present a new way of measuring thesmoothness of the numerical solution on a stencil which which implies that the weights resulting from the measure-can be used to replace (2.18)to form a new weight.As we ment (3.1)is also ﬁfth-order accurate at critical points.know,on each stencil S k ,we can construct a (r Ϫ1)th-To illustrate the different behavior of the two measure-order interpolation polynomial,which if evaluated at x ϭments (i.e.,(2.18)and (3.1))for r ϭ3in smooth monotone x j ϩ1/2,renders the approximation of h j ϩ1/2given in (2.7).regions,near critical points or near discontinuities,we com-Since total variation is a good measurement for smooth-pute the weights Ͷ0,Ͷ1,and Ͷ2for the functionness,it would be desirable to minimize the total variation for the approximation.Consideration for a smooth ﬂux and for the role of higher order variations leads us to the f j ϭͭsin 2ȏx j if 0Յx j Յ0.5,1Ϫsin 2ȏx j if 0.5Ͻx j Յ1.(3.10)following measurement for smoothness:let the interpola-tion polynomial on stencil S k be q k (x );we deﬁneat all half grid points x j ϩ1/2,where x j ϭj ⌬x ,x j ϩ1/2ϭx j ϩ⌬x /2,and ⌬x ϭ .We display the weights Ͷ0and Ͷ1in IS k ϭ͸r Ϫ1l ϭ1͵x j ϩ1/2x j Ϫ1/2h 2l Ϫ1(q (l )k )2dx ,(3.1)Fig.1(Ͷ2ϭ1ϪͶ0ϪͶ1is omitted in the picture).Notethe optimal weight for Ͷ0is C 30ϭ0.1and for Ͷ1is C 31ϭ0.6.We can see that the weights computed with (2.18)where q (l )kis the l th-derivative of q k (x ).The right-hand side (referred to as the original measurement in Fig.1)are of (3.1)is just a sum of the L 2norms of all the derivatives far less optimal than those with the new measurement,of the interpolation polynomial q k (x )over the interval especially around the critical points x ϭ , .However,near (x j Ϫ1/2,x j ϩ1/2).The term h 2l Ϫ1is to remove h -dependent the discontinuity x ϭ ,the two measurements behave factors in the derivatives of the polynomials.This is similar similarly;the discontinuous stencil always gets an almost-to,but smoother than,the total variation measurement zero weight.Moreover,for the grid point immediately left based on the L 1norm.It also renders a more accurate to the discontinuity,Ͷ0Ȃ and Ͷ1Ȃ ,which means,when WENO scheme for the case r ϭ3,when used with (2.15)only one of the three stencils is non-smooth,the other two and (2.16).stencils get O (1)weights.Unfortunately,these weights do When r ϭ2,(3.1)gives the same measurement as (2.18).not approximate a fourth-order scheme at this point.A However,they become different for r Ն3.For r ϭ3,similar situation happens to the point just to the right of (3.1)givesthe discontinuity.For simplicity of notation,we use WENO-X-3to standIS 0ϭ (f j Ϫ2Ϫ2f j Ϫ1ϩf j )2ϩ (f j Ϫ2Ϫ4f j Ϫ1ϩ3f j )2(3.2)208JIANG AND SHUFIG.1.A comparison of the two smoothness measurements.for the third-order WENO scheme (i.e.,r ϭ2,for which u 0(x )ϭsin 4(ȏx ).Again we see that WENO-RF-4is more accurate than WENO-RF-5on the coarsest grid (N ϭ20)the original and new smoothness measurement coincide)but becomes less accurate than WENO-RF-5on ﬁner grids.(where X ϭLF,Roe,RF refer,respectively,to the global The order of accuracy for WENO settles down later than Lax–Friedrichs ﬂux splitting,Roe’s ﬂux splitting,and Roe’s in the previous example.Notice that this is the example ﬂux splitting with entropy ﬁx).The accuracy of this scheme for which ENO schemes lose their accuracy [12].has been tested in [9].We will use WENO-X-4to represent the fourth-order WENO scheme of Liu et al.(i.e.,r ϭ3 4.A SIMPLE WAY FOR COMPUTING WEIGHTS FORwith the original smoothness measurement of Liu et al .)EULER SYSTEMSand WENO-X-5to stand for the ﬁfth-order WENO scheme resulting from the new smoothness measurement.In later For system (1.1)with d Ͼ1,the derivatives d f i /dx i ,i ϭsections,we will also use ENO-X-Y to denote conventional 1,...,d ,are approximated dimension by dimension;forENO schemes of ‘‘Y’’th order with ‘‘X’’ﬂux splitting.We caution the reader that the orders here are in L 1sense.So TABLE IIIENO-RF-4in our notation refers to the same scheme as ENO-RF-3in [15].Accuracy on u t ϩu x ϭ0with u 0(x )ϭsin(ȏx )In the following we test the accuracy of WENO schemes Method N L ȍerror L ȍorder L 1error L 1order on the linear equation:WENO-RF-410 1.31e-2—7.93e-3—u t ϩu x ϭ0,Ϫ1Յx Յ1,(3.11)20 3.00e-3 2.13 1.32e-3 2.5940 4.27e-4 2.81 1.56e-4 3.08u (x ,0)ϭu 0(x ),periodic.(3.12)80 5.17e-5 3.05 1.13e-5 3.79160 4.99e-6 3.37 6.88e-7 4.04320 3.44e-7 3.86 2.74e-8 4.65In Table III,we show the errors of the two schemes at t ϭ1for the initial condition u 0(x )ϭsin(ȏx )and compare WENO-RF-510 2.98e-2— 1.60e-2—them with the errors of the ﬁfth-order upstream central 20 1.45e-3 4.367.41e-4 4.4340 4.58e-5 4.99 2.22e-5 5.06scheme (referred to as CENTRAL-5in the following ta-80 1.48e-6 4.95 6.91e-7 5.01bles).We can see that WENO-RF-4is more accurate than 160 4.41e-8 5.07 2.17e-8 4.99WENO-RF-5on the coarsest grid (N ϭ10)but becomes 320 1.35e-9 5.03 6.79e-10 5.00less accurate than WENO-RF-5on the ﬁner grids.More-over,WENO-RF-5gives the expected order of accuracy CENTRAL-510 4.98e-3— 3.07e-3—20 1.60e-4 4.969.92e-5 4.95starting at about 40grid points.In this example and the 40 5.03e-6 4.99 3.14e-6 4.98one for Table IV,we have adjusted the time step to ⌬t ȁ80 1.57e-7 5.009.90e-8 4.99(⌬x )5/4so that the fourth-order Runge–Kutta in time is 160 4.91e-9 5.00 3.11e-9 4.99effectively ﬁfth-order.3201.53e-105.009.73e-115.00In Table IV,we show errors for the initial condition。

Comparison of Fifth-OrderWENO Scheme and (2)

put.Phys.doi:10.4208/cicp.110212.021112a Vol.14,No.3,pp.599-620 September2013Comparison of Fifth-Order WENO Scheme andFinite Volume WENO-Gas-Kinetic Scheme forInviscid and Viscous Flow SimulationJun Luo,Lijun Xuan and Kun Xu∗Mathematics Department,Hong Kong University of Science and Technology,Clear Water Bay,Kowloon,Hong Kong.Received11February2012;Accepted(in revised version)2November2012Available online25January2013Abstract.The development of high-order schemes has been mostly concentrated onthe limiters and high-order reconstruction techniques.In this paper,the effect of theﬂux functions on the performance of high-order schemes will be studied.Based on thesame WENO reconstruction,two schemes with differentﬂux functions,i.e.,theﬁfth-order WENO method and the WENO-Gas-kinetic scheme(WENO-GKS),will be com-pared.Theﬁfth-orderﬁnite difference WENO-SW scheme is a characteristic variablereconstruction based method which uses the Steger-Warmingﬂux splitting for invis-cid terms,the sixth-order central difference for viscous terms,and three stages Runge-Kutta time stepping for the time integration.On the other hand,theﬁnite volumeWENO-GKS is a conservative variable reconstruction based method with the sameWENO reconstruction.But,it evaluates a time dependent gas distribution functionalong a cell interface,and updates theﬂow variables inside each control volume byintegrating theﬂux function along the boundary of the control volume in both spaceand time.In order to validate the robustness and accuracy of the schemes,both meth-ods are tested under a wide range ofﬂow conditions:vortex propagation,Mach3step problem,and the cavityﬂow at Reynolds number3200.Our study shows thatboth WENO-SW and WENO-GKS yield quantitatively similar results and agree witheach other very well provided a sufﬁcient grid resolution is used.With the reduc-tion of mesh points,the WENO-GKS behaves to have less numerical dissipation andpresent more accurate solutions than those from the WENO-SW in all test cases.Forthe Navier-Stokes equations,since the WENO-GKS couples inviscid and viscous termsin a singleﬂux evaluation,and the WENO-SW uses an operator splitting technique,itappears that the WENO-SW is more sensitive to the WENO reconstruction and bound-ary treatment.In terms of efﬁciency,theﬁnite volume WENO-GKS is about4timesslower than theﬁnite difference WENO-SW in two dimensional simulations.The cur-rent study clearly shows that besides high-order reconstruction,an accurate gas evolu-tion model orﬂux function in a high-order scheme is also important in the capturing of600J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620 physical solutions.In a physicalﬂow,the transport,stress deformation,heat conduc-tion,and viscous heating are all coupled in a single gas evolution process.Therefore,it is preferred to develop such a scheme with multi-dimensionality,and uniﬁed treat-ment of inviscid and dissipative terms.A high-order scheme does prefer a high-ordergas evolution model.Even with the rapid advances of high-order reconstruction tech-niques,theﬁrst-order dynamics of the Riemann solution becomes the bottleneck forthe further development of high-order schemes.In order to avoid the weakness of thelow orderﬂux function,the development of high-order schemes relies heavily on theweak solution of the original governing equations for the update of additional degreeof freedom,such as the non-conservative gradients ofﬂow variables,which cannot bephysically valid in discontinuous regions.PACS:02.60Cb,47.11.Df,47.45.AbKey words:WENO scheme,gas-kinetic scheme,Euler equations,Navier-Stokes equations,high-order methods.J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620601 smoothness-dependent weight to each stencil and results in a(2r−1)th order accuracy. The WENO reconstruction is more accurate,efﬁcient and stable than the ENO reconstruc-tion.So far,the WENO schemes have found wide applications.The general impression is that the WENO scheme is not sensitive to theﬂuxes used,such as Lax-Friedrichs or Steger-Warming.The full accuracy of the scheme mainly depends on the order of the re-construction.One of the purpose of the current paper is to test the effect ofﬂux functions on the performance of high-order schemes.It turns out that besides the high-order ini-tial reconstruction theﬂux modeling also plays an essential role to capture accurateﬂow evolution,and to reduce the sensitive dependence of the solution on the initial recon-struction,especially in a barely resolvedﬂow region.Allﬂux functions are equivalent in a well-resolvedﬂow region,because theﬂux consistency plays a dominant role here.But their performance may deviate from each other in not well-resolved cases,such as3or4 points inside a boundary layer.Currently,in order to increase the accuracy of the WENO scheme,many attempts have been tried to develop hybrid schemes,where the WENO is used in the discontinuous region and high-order compact scheme is used in the smooth region[13,18].The aim of the hybrid scheme is to basically develop a method which could make a smooth transition between the upwind and central difference method,be-cause the fundamental basis of the1st-order Riemann solver for wave decomposition and theﬂux splitting technique contradicts with theﬂow physics in the smooth region,where the high-order spatial and temporal evolution are fully coupled.In the smooth region,the traditional central difference approximation with Cauchy-Kowalevskaya technique is far more appropriate to describeﬂow evolution than upwind concept.The Riemann solu-tion is a low order dynamic model,which is needed to model discontinuousﬂow in order to introduce enough dissipation.And this amount of dissipation is closely related to the initial jump at the cell interface[26].Without modifying the1st-order Riemann dynamic model,the development of high-order schemes becomes a game of reconstruction,since the interface jump is the only freedom many high-order schemes are able to control.Un-fortunately,there is no principle to design such an optimal and universal interface jump, and this kind of research will be endless.In order to get out of this dilemma,the use of a high-order dynamic evolution model,which could make a smooth transition between the upwind and central difference scheme,is necessary.For a second order scheme for the inviscidﬂow,we have such a high-order dynamic model,which is the generalized Riemann solver and the gas-kinetic scheme[1,11].In the past decades,a gas-kinetic scheme(GKS)based on the kinetic equation has been developed for the modeling of gas evolution process starting from a discontinu-ity[12,16,26].Theoretically,the GKS does not target to solve accurately the gas kinetic BGK model[2],but uses the kinetic equation to do the modeling around a cell interface. In the GKS evolution,the whole process from particle free transport to the Navier-Stokes (NS)or Euler solution formation has been recovered[9,28].Theﬂow regime of the gas evolution depends on the ratio of time step to the particle collision time.In the smooth ﬂow region,based on the Chapman-Enskog expansion,a time dependent gas distribu-602J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620 tion function corresponding to the NS or Euler solutions can be obtained accurately by the GKS.This limiting formulation is basically the central difference method.In the dis-continuity region,where the physical solution cannot be well resolved by the numerical cell size,theoretically it is not necessary to know the precise”macroscopic”governing equations here,because it is not needed to incorporate a precise amount of numerical dissipation.But,in such a region the gas evolution model should follow a path which is consistent with the physical one,such as keeping a non-equilibrium gas distribution function at a cell interface to cope with the dissipative mechanism of a physical shock layer[11,29].This limiting case is theﬂux vector splitting upwind method.The advan-tage of the GKS is that theﬂux function makes a smooth transition from a upwind(kinetic scale)to a central difference(hydrodynamic scale)method.Recently,based on the WENO reconstruction and high-order gas-kinetic solution [12],aﬁnite volume high-order WENO-gas-kinetic scheme(WENO-GKS)has been de-signed and tested for multi-dimensionalﬂows[15].Different from the traditional WENO scheme,the WENO technique is only used in the data reconstruction of the conservative ﬂow variables at the beginning of each time step.Starting from the high-order recon-struction,a space and time dependent gas distribution function is obtained along the tangential direction of a cell interface,from which the numericalﬂux is evaluated and used in aﬁnite volume scheme.In GKS,there is no need to use the Runge-Kutta time stepping and Gaussian point integration along a cell interface for the update ofﬂow variables inside each control volume.Among all WENO schemes we tested,theﬁnite difference WENO method with Steger-Warmingﬂux splitting gives the best results.The detailed formulation of this WENO scheme is presented in Section3.In order to show the effect of aﬂux function on the per-formance of high-order schemes,the numerical solutions of theﬁnite difference WENO-Steger-Warming scheme(WENO-SW)will be compared with that of theﬁnite volume WENO-GKS.The test cases are carefully chosen in order to test the accuracy,the shock-capturing ability,the robustness,and the stability of the schemes.In order to eliminate the numerical error due to the complicated geometry,all tests are conducted in rectangu-lar meshes.The5th-order WENO reconstruction is used on the characteristic variables in the WENO-SW,and the same WENO reconstruction is used on the conservative vari-ables in the WENO-GKS.The reason we use conservative variables for the WENO-GKS is that the scheme is not sensitive to the variable used in the reconstruction,because the evolution of the whole reconstructed curves will be followed dynamically in theﬂux evaluation.2The5th order WENO reconstructionThe5th order WENO reconstruction on a uniform rectangular mesh is a standard recon-struction method[20].Assume that Q is the variable to be reconstructed.¯Q i is the cell averaged value inJ.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620603 the ith cell.Q l i and Q r i are the two pointwise values reconstructed at the left and right interfaces of the ith cell.The5th WENO reconstruction is deﬁned as,Q r i=2∑s=0w s q(s)i,Q l i=2∑s=0˜w s˜q(s)i,whereq(0)i=16¯Qi+1−16¯Qi−1+53¯Qi+1,q(2)i=16¯Qi−1+116¯Qi−73¯Qi+2,˜q(1)i=16¯Qi−16¯Qi−2+53¯Qi,w s=αs(ǫ+βs)2,˜w s=˜αs(ǫ+βs)2,s=0,1,2,β0=134(3¯Q i−4¯Q i+1+¯Q i+2)2,β1=134(¯Q i−1−¯Q i+1)2,β2=134(¯Q i−2−4¯Q i−1+3¯Q i)2,d0=˜d2=35,d2=˜d0=1604J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620 The BGK equation in2-D isf t+ u·∇f=g−f2(u2+v2+ξ2) T,dξ=dξ1dξ2···dξK,and K is the number of degrees of internal freedom,i.e.,K=(4−2γ)/(γ−1)for2-Dﬂow andγis the speciﬁc heat ratio.Since the mass,momentum, and energy are conserved during particle collisions,f and g satisfy the conservation con-straint, (g−f)ψαd u d v dξ=0,α=1,2,3,4,(3.3)at any point in space and time.The integral solution of(3.1)isf( x,t, u,ξ)=12g l,r0((a l,r1)2+d l,r11)x2+1J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620605Here g l0and g r0are two Maxwellian distribution functions which correspond to the left and right macroscopic variables respectively.For example,g l0corresponding to W l= (ρl,(ρl U l),(ρl V l),(ρl E l))T has the formg l0=ρl λl2eλl((u−U l)2+(v−V l)2+ξ2),(3.7)whereλl equals to m/2kT l,m is the molecular mass,k is the Boltzmann constant,and T l is the temperature.For the modeling of the local equilibrium distribution function g,we can use the Taylor expansion of the equilibrium state and getg( x,t, u,ξ)=¯g+¯g¯a1x+¯g¯a2y+¯g¯At+12¯g(¯a22+¯d22)y2+¯g(¯a1¯a2+¯d12)xy+1τ t0g(−ut′,y−vt′,t′, u,ξ)e−(t−t′)/τdt′=C1¯g+C2¯g¯a1u+C1¯g¯a2y+C2¯g¯a2v+C3¯g¯A+12C1¯g(¯a22+¯d22)y2+C2¯g(¯a22+¯d22)vy+12C5¯g(¯A2+¯B)+C6¯g(¯A¯a1+¯b1)u+C3¯g(¯A¯a2+¯b2)y+C6¯g(¯A¯a2+¯b2)v,(3.10) ande−t/τf0(−ut,y−vt, u,ξ)= C7f l0(−ut,y−vt, u,ξ),u>0,C7f r0(−ut,y−vt, u,ξ),u<0,(3.11) whereC1=1−e−t/τn,C2=(t+τ)e−t/τn−τ,C3=t−τ+τe−t/τn,(3.12a) C4=(−t2−2τt)e−t/τn,C5=t2−2τt,C6=−τt(1+e−t/τn),C7=e−t/τn.(3.12b) In any numerical simulation,numerical dissipation is necessary in order to cope with the limited cell resolution.In order to add numerical dissipation but not change theﬂow606J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620 property,in(3.12),we distinguish two particle collision times.One is the physical one (τ)which corresponds to the real collision time in the Chapman-Enskog expansion for the Navier-Stokes solution,and another one is the numerical one(τn)which takes into account the effect of initial pressure jump at the cell interface.As considered in[15],the BGK-NS solver uses√τ=µ/¯p,τn=µ/¯p+β∆x√¯λ+β∆xJ.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-620607 becomesW n+1ij =W n ij+12∆x i−1∆x∆y t n+1t n12∆y j[F i−1/2(t,y)−F i+1/2(t,y)]dydt,(3.16)where F j−1/2(t,x),F j+1/2(t,x),F i−1/2(t,y),F i+1/2(t,y)are theﬂuxes along the four cell in-terfaces respectively.Because the time-dependentﬂuxes can be explicitly integrated along a cell interface in the GKS,theﬁnal scheme presents a high-orderﬂow transport through the interface within a time step without using Runge-Kutta time stepping and ﬂux construction at Gauss points.3.2WENO-Steger-Warming schemeThe WENO-SW solves the hydrodynamic equationsW t+F(W)x+G(W)y=F v(W,W x,W y)x+G v(W,W x,W y)y,(3.17) where W is the conservative variables,F and G are the inviscidﬂuxes,and F v and G v are the viscousﬂuxes.For aﬁnite-difference scheme,we need to construct both inviscid and viscousﬂuxes at the cell interface.The following WENO scheme is the5th orderﬁnite-difference WENO-Steger-Warming scheme,where the Steger-Warming splitting[23]is used to obtain the inviscidﬂuxes at the cell interface.3.2.1Inviscidﬂux reconstructionThe x-directionﬂux F can be decomposed asF=RΛLW,(3.18) whereΛ=diag[λ1,λ2,···,λn]is a diagonal matrix with eigenvalues of∂F/∂W,R is the right eigenvector matrix and L is the left one,n is the number of the equations.Then,F can be split according toF=F++F−.(3.19) In the Steger-Warming splitting,F±=RΛ±LW,Λ=Λ++Λ−.(3.20) andΛ±=diag[λ±1,λ±2,···,λ±n]withλ±i=λi± 2(i=1,···,n),(3.21)whereεis a small constant.In our simulation,ε=1.0e−3.For the construction of the numerical x-directionﬂuxˆF i+1/2,j at the cell interface (x i+1/2,y j)of the cell(i,j),we use the following steps.608J.Luo,L.Xuan and K.Xu/put.Phys.,14(2013),pp.599-6201.Split theﬂuxes F i+l,j(l=−2,···,3)in the surrounding cells to get F±i+l,j(l=−2,···,3).e the variable W i+1/2,j=(W i,j+W i+1,j)/2to calculate the left and right eigenvector matrixL i+1/2,j and R i+1/2,j at the cell interface.3.Change the conservativeﬂuxes F±i+l,j(l=−2,...,3)to the characteristicﬂuxes˜F±i+l,j(l=−2, (3)by˜F±i+l,j=L i+l/2,j F±i+l,j(l=−2,···,3).e the5th order WENO reconstruction(see Section2)on the characteristicﬂuxes˜F±i+l,j(l=−2,···,3)to get the characteristicﬂux˜F±i+1/2,j,and calculate theﬁnal characteristicﬂux˜Fi+1/2,j =˜F+i+1/2,j+˜F−i+1/2,j.5.Construct theﬂuxes for the conservative variables at the cell interface byˆFi+1/2,j=R i+1/2,j˜F i+1/2,j.In a similar way,we can get the y-directional numericalﬂuxˆG i,j+1/2at the cell inter-face(x i,y j+1/2).3.2.2Viscousﬂux reconstructionThe viscousﬂuxes F v and G v are related to the derivatives ofﬂow variables.For the NS equations,we need to calculate T x,U x,U y,V x and V y in order to get the x-direction viscousﬂux F v,where T is the temperature,U is the x-direction velocity,and V is the y-direction velocity.So,the basic step is to construct the derivatives of a variable q.Firstly,in each cell(i,j),we use a6th-order central difference to calculate the deriva-tive,i.e., ∂q60∆x.(3.22)Then,with all the derivatives,we can calculate the x-directional and y-directional viscous ﬂuxes F v i,j and G v i,j in each cell.Finally,the numerical viscousﬂuxes at the cell interfaces can be calculated by the6th order central interpolationˆF vi+1/2,j =37(F v i,j+F v i+1,j)−8(F v i−1,j+F v i+2,j)+F v i−2,j+F v i+3,j60.(3.23b)3.2.3Time evolutionTheﬁnal numericalﬂuxes at the cell interface are denoted byˆF i+1/2,j=ˆF i+1/2,j−ˆF v i+1/2,j,(3.24a)ˆG i,j+1/2=ˆG i,j+1/2−ˆG v i,j+1/2.(3.24b)With the above numericalﬂuxes,we can get the increment∆W i,j=ˆF i+1/2,j−ˆF i−1/2,j∆y.(3.25)A3rd-order TVD Runge-Kutta method is used toﬁnish the time evolution,W(1)i,j=W n i,j−∆t∆W i,j(W n),(3.26a)W(2)i,j =34[W(1)i,j−∆t∆W i,j(W(1))],(3.26b)W(n+1)i,j =13[W(2)i,j−∆t∆W i,j(W(2))].(3.26c)4Numerical examplesFor all tests in this paper,the time step is determined by the CFL condition with CFL number0.5.4.1Mach3step problemThe Mach3step problem wasﬁrst proposed by Woodward and Colella in[25].The com-putational domain is[0,3]×[0,1].A step with height0.2is located at x=0.6.The upstream velocity is(U,V)=(3,0).The adiabatic slip Euler boundary condition is implemented at all boundaries.As explained in[25],the corner of the step is the center of a rarefaction fan and itis a singular point of theﬂow.Theﬂow will be much affected by the numerical errorsgenerated just in the neighborhood of this singular point.Almost for all WENO schemes,in order to get better results,a special treatment introduced in[25]has been applied.Or,a reﬁned mesh is used around the corner.In order to avoid confusion and compare theresults of the WENO-GKS and the WENO-SW truthfully,there is no special treatmentanywhere in the current study.The parameterǫ=10−6in the5th order WENO recon-struction is used.In the WENO-GKS,α=0andβ=100for numerical collision time in (3.14).Figs.1,2and3present the numerical solutions from both WENO-SW and WENO-GKS with different number of mesh points,i.e.,120×40,60×20,and30×10.All these results are consistent with theﬂow structures in[25].Based on the theseﬁgures,the WENO-SW gives a large tilted Mach stem above the up surface of the step,and the WENO-GKS presents a much straighter and shorter one.After a few shock reﬂections, the shock front of the WENO-SW gets smeared and disappeared in the case with30×10 mesh points case.But,for the WENO-GKS,the shock reﬂection can be observed clearly in the coarse mesh case.This illustrates that WENO-SW has large numerical dissipation in comparison with WENO-GKS.In order to understand why the WENO-SW has a tiltedFigure1:Mach3the WENO-GKS5.(upper one)and the WENO-SW(lower one).In eachﬁgure,there are50contours from0.5toFigure2:Mach3the WENO-GKS (upper one)and the WENO-SW(lower one).In eachﬁgure,there are50contours from0.5to5.Mach stem,the WENO-GKS is tested for the NS solution with Reynolds number1000, but with the Euler boundary condition.Fig.4shows the comparison of WENO-SW for the inviscidﬂow and WENO-GKS for the viscousﬂow.It shows that the solution given by the WENO-SW is more close to a viscous solution.This indicates again that a high level of numerical dissipation exists in the WENO-SW.Figure3:Mach3the WENO-GKS (upper one)and the WENO-SW(lower one).In eachﬁgure,there are50contours from0.5to5.Figure4:Mach3Upper one:NS solution with Reynolds number1000by the WENO-GKS.Lower one:Euler solution by the WENO-SW.In each ﬁgure,there are50contours from0.5to5.4.2Isentropic periodic vortex propagationThis is a test for the accuracy of the Euler solutions(see[13,20]).The initial condition is given byκ(U(x,x,0),V(x,y,0))=(1,1)+e1−r2,S(x,y,0)=1,8γπ2Figure5:Isentropic periodic vortex propagation:density distribution at t=10(1period)by the WENO-GKSFigure6:Isentropic periodic vortex propagation:density distribution at t=100(10periods)by the WENO-GKS (WGKS)and the WENO-SW(WSW).where the temperature T and the entropy S are related to the densityρand the pressurep byT=pργ,and(¯x,¯y)=(x−5,y−5),r2=¯x2+¯y2,and the vortex strengthκ=5.The computational domain is[0,10]×[0,10].The periodic boundary condition is used in both directions.The numerical results with80×80and40×40cells at t=10(1period)and t=100 (10periods)are shown in Figs.5and6for both WENO-SW and WENO-GKS.In order to show the relationship between the parameterǫin the WENO reconstruction and the accuracy of the numerical solution,we test this case with two differentǫ,which are EPS1 (ǫ=10−6)and EPS2(ǫ=10−2).Theoretically,a large value ofǫpresents a reconstruction with more equally weights for different stencils.For a smoothﬂow,a smallǫmakes thereconstruction to reduce to a3rd order interpretation.A largerǫbalances the weights of different stencils and make the reconstruction to be close to the5th order interpolation. As shown in Figs.5and6,for such a smoothﬂow,in general,the reconstruction with EPS2will give more accurate and less dissipative results than that from EPS1for both WENO-SW and WENO-GKS.At the same time,for both meshes the results from WENO-GKS are more accurate than these from WENO-SW,especially at a mesh size40×40and t=100,see Fig.6.At time t=10and40×40mesh points,see Fig.5,EPS1introduces numerical dissipation and presents undershoot for both WENO-SW and WENO-GKS. However,with EPS2,the WENO reconstruction introduces overshoot for both schemes. But,in both cases due to the time accurate gas evolution model in the gas-kinetic scheme, the WENO-GKS could increase the undershoot and decrease the overshoot in compari-son with WENO-SW.In other words,the WENO-GKS has a better capacity to drive the solution in the correct direction due to the participation of the subcellﬂow distribution in the construction of the dynamicalﬂux function.More speciﬁcally,the gas-kineticﬂux function does not only depend on theﬂow variables at the cell interface,but also takes into account the whole curve evolution around the interface.Based on the above test,we can clearly realize that even though the gas-kinetic scheme is not solving the inviscid Euler equations directly,it gives accurate inviscid solution. In general,the numerical dissipation of the gas-kinetic scheme is less than the schemes based on the Riemann solvers,since the GKS can make a smooth transition from the up-wind to the central difference.This is one of the reason for the reduction of numerical dissipation in GKS in the smoothﬂow region.4.3CavityﬂowThe cavityﬂow at low Mach number is a standard test case for validating incompress-ible or low speed NSﬂow solvers.This is also a good test case for the shock capturing scheme in validating its capacity in capturing the Navier-Stokes solutions,especially with the non-linear limiters involved.For a directional splittingﬁnite volume scheme,it will be difﬁcult to present an accurate viscous solution with strong vortex structure.Fortu-nately,theﬁnite difference WENO-SW is an intrinsic multi-dimensional scheme due to its simultaneous evaluation of∂x and∂y terms of the NS equations at a grid point.At the same time,theﬁnite volume WENO-GKS has a multi-dimensionalﬂux function as well, where both x-and y-directionﬂow derivatives contributes to theﬂux in the interface normal direction.Theﬂow is bounded by a unit square and is driven by a uniform translation of the top boundary.In the simulation,the diatomic gasγ=1.4is considered.The boundary is isothermal and nonslip,and theﬂow condition is Ma=0.3and Re=3200,where Ma is the Mach number and Re is the Reynolds number.Since the benchmark solution is from incompressible NS equations,in order to avoid kinematic dissipation[6,27],most simu-lations in the past are based on the numerical methods for the incompressible equations or the artiﬁcial compressibility ones,where a continuous initial reconstruction across acell interface is assumed.However,here we are going to use the same shock capturing WENO reconstruction in the cavity simulation.This is challenge for any shock captur-ing NSﬂow solver,because the cell interface discontinuity may generate large numerical dissipation.In order to get the best results for the WENO-SW,we tested many boundary condi-tion treatments andﬁnally choose the following reconstruction.The special treatment is the following.The temperature and velocities are given directly since the boundary is isothermal and nonslip.Other data at the boundary are extrapolated from theﬂow.For the WENO-SW,the three layer interfaces close to the cavity wall but inside the compu-tational domain are specially treated.For theﬁrst layer cell interface,a3rd-order extrap-olation is used to reconstruct the conservative variableﬂuxes.For the second interface, a4th-order interpolation is used to reconstruct the conservative variableﬂuxes.For the third interface,a5th-order upwind interpolation is used to reconstruct the characteristic ﬂuxes.We use a4th-order interpolation to reconstruct the viscousﬂuxes at all of the three interfaces.For the WENO-GKS,the boundary treatment is relatively simple.A5th-order extrapolation for the conservative variables is used at the boundary.Firstly,we compare the efﬁciency in the cavity simulation for both schemes.The results are shown in Table1.As shown in the table,the WENO-GKS is about4times slower than the WENO-SW.This shows that the speed of WENO-GKS is similar to aﬁnite volume WENO scheme.As shown in[21],the computational cost of theﬁnite volume WENO scheme is indeed at least4times more expensive than theﬁnite difference one in a two-dimensional simulation.Table1:Average computational time for one time-step.WENO-GKS33×33cells 4.3276e-0037.4369e-002For the cavityﬂow simulations,we use three sets of meshes,which are101×101, 65×65,and33×33for both schemes.Fig.7shows the distributions of streamline for both schemes with a mesh of65×65points.Except around the up left corner,the streamlines from both WENO-SW and WENO-GKS are close to each other.The results of U-velocity along the vertical symmetric line at x=0.5and V-velocity along the horizontal symmetric line at y=0.5with different mesh sizes will be presented in detail.Fig.8shows U and V velocity distributions for a mesh with101×101points. The benchmark solution is from[5,24].As shown in thisﬁgure,the results from both WENO-SW and WENO-GKS are consistent with the reference solutions.The reconstruc-tion with EPS1presenting a better shock capturing property turns out to be more dissipa-tive than those with EPS2.This is clearly observed in the V-velocity for the WENO-SW. The WENO-GKS solutions are less sensitive to the values ofǫ.At such a reﬁned mesh, it is hard to distinguish the two solutions from WENO-GKS.The overall solutions from WENO-GKS are closer to the benchmark ones than the WENO-SW results.。

Schools Wise up to Ways of High中英翻译

Unit 6 Text1 Schools Wise up to Ways of High-techCheatsGeorge Mason University instructor Anne Marchant calls them “patchwork plagiarists”—the students who copy and paste together passages from articles they’ve found on the Internet, then turn in t he work as their own. She catches at least one such student every semester in their computer-science class, she said. She even discovered such plagiarism in her computer ethics course. “Certainly, cheating is persavise,” Marchant said. “It’s usually deadly obvious. The introduction will be written in broken English; then it will have this flawlessly written, almost doctoral-quality body, then a conclusion that goes back to broken English.”乔治梅森大学讲师安妮·马钱特所谓的“拼凑剽窃者”----是指学生从因特网上找到文章，经过复制、粘贴，然后把它当作自己的作业上交，她每个学期在她的计算机课上至少抓住一名这样的学生，她说，甚至在计算机道德教育课上也发现这样的“剽窃者”。

Advanced Mathematical Modeling Techniques

Advanced Mathematical ModelingTechniquesIn the realm of scientific inquiry and problem-solving, the application of advanced mathematical modeling techniques stands as a beacon of innovation and precision. From predicting the behavior of complex systems to optimizing processes in various fields, these techniques serve as invaluable tools for researchers, engineers, and decision-makers alike. In this discourse, we delve into the intricacies of advanced mathematical modeling techniques, exploring their principles, applications, and significance in modern society.At the core of advanced mathematical modeling lies the fusion of mathematical theory with computational algorithms, enabling the representation and analysis of intricate real-world phenomena. One of the fundamental techniques embraced in this domain is differential equations, serving as the mathematical language for describing change and dynamical systems. Whether in physics, engineering, biology, or economics, differential equations offer a powerful framework for understanding the evolution of variables over time. From classical ordinary differential equations (ODEs) to their more complex counterparts, such as partial differential equations (PDEs), researchers leverage these tools to unravel the dynamics of phenomena ranging from population growth to fluid flow.Beyond differential equations, advanced mathematical modeling encompasses a plethora of techniques tailored to specific applications. Among these, optimization theory emerges as a cornerstone, providing methodologies to identify optimal solutions amidst a multitude of possible choices. Whether in logistics, finance, or engineering design, optimization techniques enable the efficient allocation of resources, the maximization of profits, or the minimization of costs. From linear programming to nonlinear optimization and evolutionary algorithms, these methods empower decision-makers to navigate complex decision landscapes and achieve desired outcomes.Furthermore, stochastic processes constitute another vital aspect of advanced mathematical modeling, accounting for randomness and uncertainty in real-world systems. From Markov chains to stochastic differential equations, these techniques capture the probabilistic nature of phenomena, offering insights into risk assessment, financial modeling, and dynamic systems subjected to random fluctuations. By integrating probabilistic elements into mathematical models, researchers gain a deeper understanding of uncertainty's impact on outcomes, facilitating informed decision-making and risk management strategies.The advent of computational power has revolutionized the landscape of advanced mathematical modeling, enabling the simulation and analysis of increasingly complex systems. Numerical methods play a pivotal role in this paradigm, providing algorithms for approximating solutions to mathematical problems that defy analytical treatment. Finite element methods, finite difference methods, and Monte Carlo simulations are but a few examples of numerical techniques employed to tackle problems spanning from structural analysis to option pricing. Through iterative computation and algorithmic refinement, these methods empower researchers to explore phenomena with unprecedented depth and accuracy.Moreover, the interdisciplinary nature of advanced mathematical modeling fosters synergies across diverse fields, catalyzing innovation and breakthroughs. Machine learning and data-driven modeling, for instance, have emerged as formidable allies in deciphering complex patterns and extracting insights from vast datasets. Whether in predictive modeling, pattern recognition, or decision support systems, machine learning algorithms leverage statistical techniques to uncover hidden structures and relationships, driving advancements in fields as diverse as healthcare, finance, and autonomous systems.The application domains of advanced mathematical modeling techniques are as diverse as they are far-reaching. In the realm of healthcare, mathematical models underpin epidemiological studies, aiding in the understanding and mitigation of infectious diseases. From compartmental models like the SIR model to agent-based simulations, these tools inform public health policies and intervention strategies, guiding efforts to combat pandemics and safeguard populations.In the domain of climate science, mathematical models serve as indispensable tools for understanding Earth's complex climate system and projecting future trends. Coupling atmospheric, oceanic, and cryospheric models, researchers simulate the dynamics of climate variables, offering insights into phenomena such as global warming, sea-level rise, and extreme weather events. By integrating observational data and physical principles, these models enhance our understanding of climate dynamics, informing mitigation and adaptation strategies to address the challenges of climate change.Furthermore, in the realm of finance, mathematical modeling techniques underpin the pricing of financial instruments, the management of investment portfolios, and the assessment of risk. From option pricing models rooted in stochastic calculus to portfolio optimization techniques grounded in optimization theory, these tools empower financial institutions to make informed decisions in a volatile and uncertain market environment. By quantifying risk and return profiles, mathematical models facilitate the allocation of capital, the hedging of riskexposures, and the management of investment strategies, thereby contributing to financial stability and resilience.In conclusion, advanced mathematical modeling techniques represent a cornerstone of modern science and engineering, providing powerful tools for understanding, predicting, and optimizing complex systems. From differential equations to optimization theory, from stochastic processes to machine learning, these techniques enable researchers and practitioners to tackle a myriad of challenges across diverse domains. As computational capabilities continue to advance and interdisciplinary collaborations flourish, the potential for innovation and discovery in the realm of mathematical modeling knows no bounds. By harnessing the power of mathematics, computation, and data, we embark on a journey of exploration and insight, unraveling the mysteries of the universe and shaping the world of tomorrow.。

Mathematical modeling of drying of pretreated

Mathematical modeling of drying of pretreated and untreated pumpkinT.Y.Tunde-Akintunde &G.O.OgunlakinRevised:13February 2011/Accepted:26April 2011#Association of Food Scientists &Technologists (India)2011Abstract In this study,drying characteristics of pretreated and untreated pumpkin were examined in a hot-air dryer at air temperatures within a range of 40–80°C and a constant air velocity of 1.5m/s.The drying was observed to be in the falling-rate drying period and thus liquid diffusion is the main mechanism of moisture movement from the internal regions to the product surface.The experimental drying data for the pumpkin fruits were used to fit Exponential,General exponential,Logarithmic,Page,Midilli-Kucuk and Parabolic model and the statistical validity of models tested were determined by non-linear regression analysis.The Parabolic model had the highest R 2and lowest χ2and RMSE values.This indicates that the Parabolic model is appropriate to describe the dehydration behavior for the pumpkin.Keywords Pumpkin fruits .Hot air drying .Effective diffusivity .Mathematical modellingNomenclature a Drying constant b Drying constant c Drying constant DR Drying rate (g water/g dry matter*h)k Drying constant,1/min M e Equilibrium moisture content(kg water/kg dry matter)M i Initial moisture content (kg water/kg dry matter)M R Dimensionless moisture ratioMR exp,i Experimental dimensionless moisture ratio MR pre,i Predicted dimensionless moisture ratio M t Moisture content at any time of drying (kg water/kg dry matter)M t +dtMoisture content at t +dt (kg water/kg dry matter)N Number of observationsn Drying constant,positive integer R 2Coefficient of determination t Time (min)W Amount of evaporated water (g)W 0Initial weight of sample (g)W 1Sample dry matter mass (g)z Number of constants χ2Reduced chi-squareIntroductionPumpkin (cucurbita mixta )is a fruit rich in Vitamin A,potassium,fiber and carbohydrates.It is a versatile fruit that can be used for either animal feed or for human consumption as a snack,or made into soups,pies and bread.The high moisture content of the fruit makes it susceptible to deterioration after harvest.The most common form of preservation being done locally is drying.It is a means of improving storability by increasing shelf-life of the food product.Dried products can be stored for months or even years without appreciable loss of nutrients.Drying also assist in reducing post harvest losses of fruits and vegetables especially which can be as high as 70%(Tunde-Akintunde and Akintunde 1996).Sun drying is the common method of drying in the tropical region.However the process is weather depen-T.Y .Tunde-Akintunde (*):G.O.Ogunlakin Department of Food Science and Engineering,Ladoke Akintola University of Technology,PMB 4000,Ogbomoso,Oyo State,Nigeria e-mail:toyositunde@J Food Sci TechnolDOI 10.1007/s13197-011-0392-2dent Also,the drying of food products with the use of sun-drying usually takes a long time thus resulting in products of low quality.There is a need for suitable alternatives in order to improve product quality.Hot air dryers give far more hygienic products and provide uniform and rapid drying which is more suitable for the food drying processes(Kingsly et al.2007a;Doymaz 2004a).Many types of hot-air dryers are being used for drying agricultural products.However the design of such driers for high-moisture foods constitutes a very complex problem owing to the characteristics of vegetable tissues.One of the most important factor that needs to be considered in the design of driers is the proper prediction of drying rate for shrinking particles and hence the appropriate drying time in the dryer needs to be investigated.This can be achieved by determination of the drying characteristics of the food material.Therefore the thin layer drying studies which normally form the basis of understanding the drying characteristics of food materials has to be carried out for each food material.Studies have been carried out on thin layer drying of some food materials(Gaston et al.2004), fruits(Doymaz2004a,c;Simal et al.2005),leaves and grasses(Demir et al.2004).Though there has been some literature on drying of pumpkin(Alibas2007;Doymaz2007b;Sacilik2007),the selection of appropriate model to describe the drying process for the variety common in the country and within the experimental conditions considered in this study is yet to be done.Thin layer drying models used in the analysis of drying characteristics are usually theoretical,semi-theoretical or purely empirical.A number of semi-theoretical drying models have been widely used by various researchers(Sharma and Prasad2004;Simal et al.2005;Sogi et al.2003;Togrul and Pehlivan2004).A number of pre-treatments can be applied depending on the food to be dried,its end use,and availability(Doymaz 2010).Pretreatment of food materials which includes; blanching,osmotic dehydration,soaking in ascorbic acid before or on drying have been investigated to prevent the loss of colour by inactivating enzymes and relaxing tissue structure.This improves the effect of drying by reducing the drying time and gives the eventual dried products of good nutritional quality(Kingsly et al.2007b;Doymaz 2010).Various commercially used pretreatments include potassium and sodium hydroxide,potassium meta bisul-phate,potassium carbonate,methyl and ethyl ester emul-sions,ascorbic and citric acids,(Kingsly et al.2007a,b; Doymaz2004a,b;El-Beltagy et al.2007).However non-chemical forms of pretreatment are generally preferred especially among small-scale processors in the tropics. Blanching as a pre-treatment is used to arrest some physiological processes before drying vegetables and fruits.It is a heat pre-treatment that inactivates the enzymes responsible for commercially unacceptable darkening and off-flavours.Blanching of fruits and vegetables is generally carried out by heating them with steam or hot water(Tembo et al.2008).However other forms of blanching i.e.oil-water blanching have been used by Akanbi et al.(2006)in a previous study for pretreating chilli.These forms of blanching were observed to have an effect on the drying rate and quality of the dried chili. However oil-water blanching pretreatments have not been studied for pumpkin.The aim of this study was:(a)to study the effect of the different blanching pre-treatments on the drying times and rate,and(b)to fit the experimental data to seven mathematical models available in the literature.Materials and methodsExperimental procedureFresh pumpkin fruit were purchased from Arada,a local market in Ogbomoso,Nigeria.The initial average moisture content of fresh pumpkin samples was deter-mined by oven drying method(AOAC1990),and it was found to be91.7%(wet basis)or10.90(g water/g dry matter).The samples were washed and peeled after which the seeds were removed.The pumpkin was then sliced into pieces of5mm×5mm dimensions.The blanching pre-treatments for inactivation of enzymes are as indicated below while untreated samples(UT)were used as the control.i)Samples submerged in boiling water for3minutes andcooling immediately in tap water-(WB)ii)Samples steamed over boiling water in a water bath (WBH14/F2,England)for3minutes and cooled immediately in tap water-(SB)iii)Samples dipped for3minutes in a homogenized mixture of oil and water of ratio1:20(v/v)with0.1g of butylated hydroxyl anisole(BHA)heated to95°C-(O/W B)In all the pretreatments the excess water was removed by blotting the pretreated pumpkin samples with tissue paper.The moisture contents after pretreatment were12.69(g water/g dry matter)for water and oil-water blanching and11.5(g water/g dry matter)for steam blanching.Drying procedureThe drying experiments of pumpkin samples were carried out in a hot-air dryer(Gallenkamp,UK)having three tiers of trays. Perforated trays having an area of approximately0.2m2wereJ Food Sci Technolplaced on each tier and the trays were filled with a single layer of the pumpkin samples.The air passes from the heating unit and is heated to the desired temperature and channeled to the drying chamber.The hot air passes from across the surface and perforated bottom of the drying material and the direction of air flow was parallel to the samples.The samples utilised for each experimental condition weighed 200±2g.Drying of the pumpkin was carried out at drying temperatures of 40to 800C with 20°C increment,and a constant air velocity of 1.5m/s for all circumstances.The dryer was adjusted to be selected temperature for about half an hour before the start of experiment to achieve the steady state conditions.Weight loss of samples was measured at various time intervals,ranging from 30min at the beginning of the drying to 120min during the last stages of the drying process by means of a digital balance (PH Mettler)with an accuracy of ±0.01g.The samples were taken out of the dryer and weighed during each of the time intervals.Drying was stopped when constant weight was reached with three consecutive readings.The experiments were repeated twice and the average of the moisture ratio at each value was used for drawing the drying curves.Model nameModelReferencesExponential modelMR ¼exp Àkt ðÞEl-Beltagy et al.(2007)Generalized exponential model MR ¼Aexp Àkt ðÞShittu and Raji (2008)Logarithmic model MR ¼aexp ðÀkt Þþc Akpinar and Bicer (2008)(Page ’s model)MR ¼exp Àkt n ðÞSingh et al.(2008)Midilli –Kucuk model MR ¼aexp Àkt n ðÞþbt Midilli and Kucuk (2003)Parabolic modelMR ¼a þbt þct 2Sharma and Prasad (2004),Doymaz (2010)Table 1Mathematical models fitted to pretreated pumpkin drying curves100200300400Drying Time (min)WB SBO/W B UT0.10.20.30.40.50.60.70.80.910100200300Drying Time (min)M o i s t u r e R a t i o00.10.20.30.40.50.60.70.80.91M o i s t u r e R a t i oWB SB O/W B UT00.10.20.30.40.50.60.70.80.91M o i s t u r e R a t i o100200300Drying Time (min)WB SB O/W B UTbacFig.1Drying curve of pump-kin slices dried at (a )40,(b )60and (c )80°C (WB-water blanched,SB-steam blanched,O/W B-oil water blanched,UT-untreated).Each observation is a mean of two replicate experiments (n =2)J Food Sci TechnolMathematical modelingThe moisture content at any time of drying (kg water/kg dry matter),M t was calculated as follows:M t ¼W o ÀW ðÞÀW 1W 1ð1ÞWhere W 0is the initial weight of sample,W is the amount of evaporated water and W 1is the sample dry matter mass.The reduction of moisture ratio with drying time was used to analyse the experimental drying data.The moisture ratio,M R,was calculated as follows:M R ¼M t ÀM eM i ÀM e¼exp ÀKt ðÞð2ÞWhere M t ,M i and M e are moisture content at any time of drying (kg water/kg dry matter),initial moisture content (kg water/kg dry matter)and equilibrium moisture content (kg water/kg dry matter),respectively.The equilibrium moisture contents (EMCs)were deter-mined by drying until no further change in weight wasobserved for the pumpkin samples in each treatment and drying condition (Hii et al.2009)The drying constant,K ,is determined by plotting experimental drying data for each pretreatment and drying temperature in terms of Ln M R against time (t)where M R is moisture ratio.The slope,k,is obtained from the straight line graphs above.The drying rate for the pumpkin slices was calculated as follows:DR ¼M t þdt ÀM tdtð3ÞWhere DR is drying rate,M t +dt is moisture content at t +dt (kg water/kg dry matter),t is time (min).The moisture ratio curves obtained were fitted with six semi-theoretical thin layer-drying models,Exponential,Generalized exponential,Page,Logarithmic,Midilli -Kucuk and Parabolic models (Table 1)in order to describe the drying characteristics of pretreated pumpkin.These linear forms of these models were fitted in the experimental data using regression technique.To evaluate the models,a nonlinear regression procedure was per-formed for six models using SPSS (Statistical Package for00.020.040.060.080.10.12Drying Time (h)D r y i n g r a t e (g w a t e r /g D M .h )WB SB O/W B UT00.020.040.060.080.10.120.140.160.180246Drying time (h)D r y i n g r a t e (g w a t e r /g D M .h )WB SB O/W B UT00.010.020.030.040.050.060246Drying Time (h)D r y i n g R a t e (g w a t e r /g D M .h )WB SB O/W B UTbacFig.2Drying rate curve of pretreated pumpkin dried at (a )40,(b )60and (c )80°C (WB-water blanched,SB-steam blanched,O/W B-oil water blanched,UT-untreated).Each observation is a mean of two replicate experiments (n =2)J Food Sci Technolsocial scientists)11.5.1software package.The correlation coefficient (R 2)was one of the primary criteria to select the best equation to account for variation in the drying curves of dried samples.In addition to R 2,other statistical parameters such as reduced mean square of the deviation (χ2)and root mean square error (RMSE)were used to determine the quality of the fit.The higher the value of R 2and the lower the values of χ2and RMSE were chosen as the criteria for goodness of fit.(Togrul and Pehlivan 2002;Demir et al.2004;Doymaz 2004b ;Goyal et al.2006).The above parameters can be calculated as follows:#2¼PN i ¼1MR ðexp ;i ÞÀMR ðpred ;i ÞÀÁN Àz2ð4ÞRMSE ¼1N XN i ¼1MR pred ;i ÀMR exp ;i ÀÁ2"#12ð5ÞWhere MR exp,i and MR pre,i are experimental and pre-dicted dimensionless moisture ratios,respectively;N is number of observations;z is number of constantsResults and discussion Drying characteristicsThe characteristics drying curves showing the changes in moisture ratio of pretreated pumpkin with time at drying temperatures of 40,60and 80°C are given in Fig.1.Figure 2show the changes in drying rate as a function of drying time at the same temperatures.It is apparent that moisture ratio decreases continuously with drying time.According to the results in Fig.1,the drying air temperature and pre-treatments had a significant effect on the moisture ratio of the pumpkin samples as expected.This is in agreement with the observations of Doymaz (2010)for apple slices.The figures show that the drying time decreased with increase in drying air temperature.Similar results were also reported for food products by earlierTable 2Drying constant (K)values for pretreatment methods and drying temperaturesDrying Temperature (°C)Pretreatment method Steam Blanching Water Blanching Oil-Water Blanching Untreated 400.45290.44580.42260.4029600.62230.61630.60150.5637800.78840.78710.78370.7803Model nameCoefficient of determination (R 2)Reducedchi-square (χ)Root mean square error (RMSE)40°CGeneralized Exponential Model 0.95180.0058240.071384Exponential Model 0.96550.0051960.062424Logarithmic Model 0.96840.0056210.059272Page Model0.973750.0032790.048397Midilli –Kucuk model 0.863980.0277520.109059Parabolic model 0.99630.0004730.0164460°CGeneralized Exponential Model 0.98740.0013820.035049Exponential Model 0.98820.0013460.032351Logarithmic Model 0.99420.0007850.022878Page Model0.980390.001110.028847Midilli –Kucuk model 0.818570.0145140.085189Parabolic model 0.99410.0004480.01672680°CGeneralized Exponential Model 0.95640.0038620.058588Exponential Model 0.96220.0038610.054796Logarithmic Model 0.97250.0107950.082138Page Model0.95950.0086160.080386Midilli –Kucuk model 0.937020.022370.105759Parabolic model0.98650.0032920.04685Table 3Curve fitting criteria forthe various mathematical models and parameters for pumpkin pre-treated with water blanching and dried at temperatures of 40,60and 80°CJ Food Sci TechnolModel nameCoefficient of determination (R 2)Reducedchi-square (χ)Root mean square error (RMSE)40°CGeneralized Exponential Model 0.95930.0048430.065095Exponential Model 0.96890.0042740.056615Logarithmic Model 0.968960.0050650.056262Page Model0.95990.0045870.057237Midilli –Kucuk model 0.818470.0342380.121134Parabolic model 0.99490.0006490.01925260°CGeneralized Exponential Model 0.987370.0013540.034688Exponential Model 0.98760.0014480.033563Logarithmic Model 0.990810.0008640.023994Page Model0.97490.0012870.031073Midilli –Kucuk model 0.982270.0169220.091982Parabolic model 0.99410.0005980.01933680°CGeneralized Exponential Model 0.94360.0054840.069818Exponential Model 0.951870.0108440.060609Logarithmic Model 0.958190.0052010.082325Page Model0.922370.0065670.070179Midilli –Kucuk model 0.955020.0210930.102696Parabolic model0.979590.0047230.058882Table 4Curve fitting criteria for the various mathematical models and parameters for pumpkin pre-treated with steam blanching and dried at temperatures of 40,60and 80°CModel nameCoefficient of determination (R 2)Reducedchi-square (χ)Root mean square error (RMSE)40°CGeneralized Exponential Model 0.94540.0067290.076734Exponential Model 0.96030.0059530.066819Logarithmic Model 0.961750.0069530.065922Page Model0.96450.0045960.057295Midilli –Kucuk model 0.838010.0328130.118587Parabolic model 0.99560.0008220.02167260°CGeneralized Exponential Model 0.976650.0017910.039903Exponential Model 0.982710.002020.039642Logarithmic Model 0.986950.0016620.033287Page Model0.974690.0012580.03072Midilli –Kucuk model 0.88840.0083740.064708Parabolic model 0.99280.0004530.01682280°CGeneralized Exponential Model 0.947850.0047940.065279Exponential Model 0.95540.00480.061101Logarithmic Model 0.96250.0120260.086696Page Model0.932230.0106440.089346Midilli –Kucuk model 0.94940.0200940.100234Parabolic model0.983570.0043540.053874Table 5Curve fitting criteria for the various mathematical models and parameters for pumpkin pre-treated with oil-water blanching and dried at temperatures of 40,60and 80°CJ Food Sci TechnolModel nameCoefficient of determination (R 2)Reducedchi-square (χ)Root mean square error (RMSE)40°CGeneralized Exponential Model 0.948760.0058960.071829Exponential Model 0.95850.0056790.065264Logarithmic Model 0.959110.0067350.06488Page Model0.94120.0064750.068007Midilli –Kucuk model 0.78980.0399250.130808Parabolic model 0.985950.0018690.03267860°CGeneralized Exponential Model 0.97350.0024360.046534Exponential Model 0.97490.0028840.047363Logarithmic Model 0.981070.0011850.028035Page Model0.95840.0021470.040127Midilli –Kucuk model 0.76690.0163070.090296Parabolic model 0.990390.0011790.0272280°CGeneralized Exponential Model 0.949980.0045840.063836Exponential Model 0.95720.0047680.060896Logarithmic Model 0.959850.0050460.057999Page Model0.920160.0077380.069542Midilli –Kucuk model 0.95070.0179590.09476Parabolic model0.97750.0042850.056692Table 6Curve fitting criteria for the various mathematical models and parameters for untreated pumpkin and dried at temper-atures of 40,60and 80°CDrying Time (h)M o i s t u r e R a t i oDrying Time (h)M o i s t u r e R a t i oDrying Time (h)M o i s t u r e R a t i obacFig.3Comparison of experi-mental and predicted moisture ratio values using Parabolic model for pretreated pumpkin dried at (a )40,(b )60and (c )80°C (WBE-water blanched experimental,SBE-steam blanched experimental,O/W BE-oil water blanchedexperimental,UTE-untreated experimental,WBP-waterblanched predicted,SBP-steam blanched predicted,O/W BP-oil water blanched predicted,UTP-untreated predicted).Each observation is a mean of two replicate experiments (n =2)J Food Sci Technolresearchers(Sacilik and Elicin2006;Lee and Kim2009; Kumar et al.2010).The drying time required to lower the moisture ratio of water blanched samples to0.034when using an air temperature of40°C(5h)was approximately twice that required at a drying air temperature of80°C (2.5h).This same trend occurred for both untreated and other pretreatment methods.The drying time required to reach moisture ratio of0.018for drying temperature of60°C for samples pretreated with steam blanching was3h while the corresponding values for water blanched,oil-water blanched and control samples were3.5,3.9and4h respectively.The difference in drying times of pre-treated samples with steam blanching was16.7%,30%and33.3%shorter than water blanched,oil-water blanched and control samples,re-spectively.Similar trends were observed at drying temper-atures of40and80°C.All the pretreated samples had higher drying curves than the untreated samples generally(Fig.1).This is an indication of the fact that various forms of blanching pretreatments increase the drying rate for pumpkin samples. This is similar to the observations of Goyal et al.(2008), Doymaz(2007a),Kingsly et al.(2007a),Doymaz(2004b) for apples,tomato,peach slices and mulberry fruits.The difference in drying is more pronounced at the initial stages of drying when the major quantity of water is evaporated while at the latter stages of drying the difference in the amount of water evaporated is not as pronounced as the early drying stages.The difference in the drying of the different pretreatments experienced in the initial stages becomes less pronounced with increase in drying temper-ature.This may be because at higher temperatures the driving force due to diffusion from the internal regions to the surface is higher thus overcoming hindrances to drying more effectively resulting in more uniform drying.The K values for pretreated samples were higher than that of untreated samples for all the drying temperatures(Table2). This confirms the fact that the various forms of blanching pretreatments increased the rate at which drying took place.Analysis of the drying rate curves(Fig.2)showed no constant-rate period indicating that drying occurred during the falling-rate period.Therefore,it can be considered a diffusion-controlled process in which the rate of moisture removal is limited by diffusion of moisture from inside to surface of the product.This is similar to the results reported for various agricultural products such Amasya red apples (Doymaz2010),peach(Kingsly et al.2007a),yam(Sobukola et al.2008),and tumeric(Singh et al.2010).Fitting of drying curveSPSS statistical software package for non-linear regression analysis was used to fit moisture ratio against drying time to determine the constants of the four selected drying models.The R2,χ2and RMSE used to determine the goodness of fit of the models are shown in Tables3,4,5 and6.The Parabolic model gave the highest R2value which varied from0.9963to0.9775for experimental conditions considered in this study.The values ofχ2and RMSE for the Parabolic model which varied from0.000448 to0.004723and0.01644to0.058882respectively were the lowest for all the models considered.From the tables,it is obvious that the Parabolic model therefore represents the drying characteristics of pretreated pumpkin(for individual drying runs)better than the other models(Generalized exponential,Exponential,Logarithmic,Page or Midilli–Kucuk)considered in this study.The comparison between experimental moisture ratios and predicted moisture ratios obtained from the Parabolic model at drying air temperature of40°C,60°C and80°C for pumpkin samples are shown in Fig.3.The suitability of the Parabolic model for describing the pumpkin drying behaviour is further shown by a good conformity between experimental and predicted moisture ratios as seen in Fig.3.This is similar to the observations of Doymaz(2010)for drying of red apple slices at55,65and75°C.ConclusionThe effect of temperature and pre-treatments on thin layer drying of pumpkin in a hot-air dryer was investigated. Increase in drying temperature from40to80°C decreased the drying time from5hours to4hours for all the samples considered.The pretreated samples dried faster than the untreated samples.Samples pretreated with steam blanch-ing had shorter drying times(hence higher drying rates) compared to water blanched,oil-water blanched and control samples The entire drying process occurred in falling rate period and constant rate period was not observed.The suitability of four thin-layer equations to describe the drying behaviour of pumpkin was investigated.The model that had the best fit with highest values of R2and lowest values ofχ2,MBE and RMSE was the Parabolic model. Thus this model was selected as being suitable to describe the pumpkin drying process for the experimental conditions considered.ReferencesAkanbi CT,Adeyemi RS,Ojo A(2006)Drying characteristics and sorption isotherm of tomato slices.J Food Eng73:141–146 Akpinar AK,Bicer Y(2008)Mathematical modelling of thin layer drying process of long green pepper in solar dryer and under open sun.Energy Conver Manag49:1367–1375Aliba I(2007)Microwave,air and combined microwave–air-drying parameters of pumpkin slices.LWT40:1445–1451J Food Sci TechnolAOAC(1990)Official methods of analysis,15th edn.Association of Official Analytical Chemists,ArlingtonDemir V,Gunhan T,Yagcioglu AK,Degirmencioglu A(2004) Mathematical modeling and the determination of some quality parameters of air-dried bay leaves.Biosys Eng88(3):325–335 Doymaz I(2004a)Drying kinetics of white mulberry.J Food Eng61(3):341–346Doymaz I(2004b)Pretreatment effect on sun drying of mulberry fruit (Morus alba L.).J Food Eng65(2):205–209Doymaz I(2004c)Convective air drying characteristics of thin layer carrots.J Food Eng61(3):359–364Doymaz I(2007a)Air-drying characteristics of tomatoes.J Food Eng 78:1291–1297Doymaz I(2007b)The kinetics of forced convective air-drying of pumpkin slices.J Food Eng79:243–248Doymaz I(2010)Effect of citric acid and blanching pre-treatments on drying and rehydration of Amasya red apples.Food Bioprod Proc 88(2–3):124–132El-Beltagy A,Gamea GR,Amer Essa AH(2007)Solar drying characteristics of strawberry.J Food Eng78:456–464Gaston AL,Abalone RM,Giner SA,Bruce DM(2004)Effect of modelling assumptions on the effective water diffusivity in wheat.Biosys Eng88(2):175–185Goyal RK,Kingsly ARP,Manikantan MR,Ilyas SM(2006)Thin layer drying kinetics of raw mango slices.Biosys Eng95(1):43–49 Goyal RK,Mujjeb O,Bhargava VK(2008)Mathematical modeling of thin layer drying kinetics of apple in tunnel dryer.Int J Food Eng 4(8):Article8Hii CL,Law CL,Cloke M,Suzannah S(2009)Thin layer drying kinetics of cocoa and dried product quality.Biosys Eng102:153–161 Kingsly RP,Goyal RK,Manikantan MR,Ilyas SM(2007a)Effects of pretreatments and drying air temperature on drying behaviour of peach slice.Int J Food Sci Technol42:65–69Kingsly ARP,Singh R,Goyal RK,Singh DB(2007b)Thin-layer drying behaviour of organically produced tomato.Am J Food Technol 2:71–78Kumar R,Jain S,Garg MK(2010)Drying behaviour of rapeseed under thin layer conditions.J Food Sci Technol47(3):335–338 Lee JH,Kim HJ(2009)Vacuum drying kinetics of Asian white radish (Raphanus sativus L.)slices.LWT-Food Sci Technol42:180–186Midilli A,Kucuk H(2003)Mathematical modeling of thin layer drying of pistachio by using solar energy.Energy Conver Manag 44(7):1111–1122Sacilik K(2007)Effect of drying methods on thin-layer drying characteristics of hull-less seed pumpkin(Cucurbita pepo L.).J Food Eng79:23–30Sacilik K,Elicin AK(2006)The thin layer drying characteristics of organic apple slices.J Food Eng73:281–289Sharma GP,Prasad S(2004)Effective moisture diffusivity of garlic cloves undergoing microwave convective drying.J Food Eng65(4):609–617Shittu TA,Raji AO(2008)Thin layer drying of African Breadfruit (Treculia africana)seeds:modeling and rehydration capacity.Food Bioprocess Technol:1–8.doi:10.1007/s11947-008-0161-zSimal S,Femenia A,Garau MC,Rossello C(2005)Use of exponential,page and diffusion models to simulate the drying kinetics of kiwi fruit.J Food Eng66(3):323–328Singh S,Sharma R,Bawa AS,Saxena DC(2008)Drying and rehydration characteristics of water chestnut(Trapa natans)as a function of drying air temperature.J Food Eng87:213–221Singh G,Arora S,Kumar S(2010)Effect of mechanical drying air conditions on quality of turmeric powder.J Food Sci Technol47(3):347–350Sobukola OP,Dairo OU,Odunewu A V(2008)Convective hot air drying of blanched yam slices.Int J Food Sci Technol43:1233–1238 Sogi DS,Shivhare US,Garg SK,Bawa SA(2003)Water sorption isotherms and drying characteristics of tomato seeds.Biosys Eng 84(3):297–301Tembo L,Chiteka ZA,Kadzere I,Akinnifesi FK,Tagwira F(2008) Blanching and drying period affect moisture loss and vitamin C content in Ziziphus mauritiana(Lamk.).Afric J Biotech7:3100–3106Togrul IT,Pehlivan D(2002)Mathematical modeling of solar drying of apricots in thin layers.J Food Eng55:209–216Togrul IT,Pehlivan D(2004)Modeling of thin layer drying kinetics of some fruits under open air sun drying process.J Food Eng65(3):413–425Tunde-Akintunde TY,Akintunde BO(1996)Post-harvest losses of food crops:sources and solutions.Proceedings of the Annual Conference of the Nigerian Society of Agricultural Engineers, Ile-Ife,Nigeria from November19–22,1996.V ol18:258–261J Food Sci Technol。

高英unit 3

Buildings in Old Shanghai
Skyscrapers in Shanghai
Crowded r-3) Globalization is a reality but it is not something completely new. What is new is the speed and scope of changes. Definition: a worldwide reformation of cultures, a tectonic shift of habits and dreams; a wild assortment of changes in politics, business, health, entertainment
Culture manifestation
The first wave--- Agrarian and agricultural culture replaced huntergatherer culture
Centralization concentration
The second wave
Sesame Street uses combinations of animation and live actors to stimulate young children's minds, improve their letter and word recognition, basic arithmetic, geometric forms, classification, simple problem solving, and socialization by showing children or people in their everyday lives. Since the show's inception, other instructional goals have been basic life skills, such as how to cross the street safely, proper hygiene, healthy eating habits, and social skills.

simulations in the high-temperature phase

J.Phys.A:Math.Gen.31(1998)10045–10052.Printed in the UK PII:S0305-4470(98)94980-0 Glauber dynamics of the SK model:theory and simulations in the high-temperature phaseGrzegorz SzamelDepartment of Chemistry,Colorado State University,Fort Collins,CO80523,USAReceived9June1998Abstract.Glauber dynamics of the Sherrington–Kirkpatrick spin-glass model in the high-temperature phase is studied using theory and computer simulations.The theoretical approachfollows the spirit of Mori’s continuous fraction expansion.Its predictions agree very well withthe computer simulation results.1.IntroductionThe dynamical properties of the Sherrington–Kirkpatrick(SK)spin-glass model[1]have been a subject of continuous interest in recent years[2,3].Most of the theoretical studies considered Langevin dynamics of the soft-spin version of the SK model[4–6].The soft-spin version,while showing very interesting dynamical properties[6],lacks the original motivation of the SK model:neither its statics nor its dynamics is exactly solvable.The results are obtained perturbatively with respect to the four-spin coupling constant u.To recover the Ising limit one has to let u approach inﬁnity.In practice,this procedure allows one to analyse the long-time asymptotic behaviour of the spin correlations.It is not well suited to study the time dependence for all times(even for T T c)or to calculate the so-called absolute frequency scale.Glauber dynamics of the hard-spin version of the SK model has been analysed in the original SK paper[1]and in other early spin-glass publications[7].These works used an ad hoc method to introduce a dynamical version of the Onsager reactionﬁeld term into the equations of motion.A more systematic approach was started by Sommers[8].However, his method was criticized by L usakowski[9]and its validity is uncertain.Recently a novel approach to Glauber dynamics of spin glasses has been proposed by Coolen,Sherrington et al[10,11].The simple version of their theory[10]describes very well the order parameter ﬂow direction above the de Almeida–Thouless(AT)[12]line but misses the slowing down which sets in when the former line is approached from above.The more advanced version [11]agrees well with the simulation data for short times but it remains to be seen whether it predicts divergent relaxation times at and below the AT line.In another interesting recent study Nishimori and Yamana[13]analysed Glauber dynamics of the SK model via a high-temperature expansion.Their method has an obvious virtue of a systematic expansion in a small parameter.However,it cannot be applied near or below the spin-glass transition.The original motivation for the work presented here was to improve upon the simple CS theory[10].Brieﬂy,we have noticed that the simple CS theory is quite similar in spirit to the Enskog kinetic theory of hard sphereﬂuids[14].We expected that by using kinetic theory techniques we might be able to incorporate divergence of the relaxation times 0305-4470/98/5010045+08$19.50c 1998IOP Publishing Ltd1004510046G Szamelinto the CS theory.It turned out that a more efﬁcient way to proceed was to depart from one of the ingredients of the CS theory:while CS derive equations of motion for sample-averaged quantities we postpone sample averaging until after the equations of motion are solved.It should be emphasized that in all other respects our zeroth-order approximation follows the spirit of the Enskog theory(and that of the simple CS theory).In theﬁrst-order approximation we include terms that would be analogous to ring and repeated ring terms in the kinetic theory.Because we are dealing with sample-dependent quantities the structure of our theory is more complicated than that of the CS theory.Therefore we restrict ourselves to studying time-dependent spin correlations in equilibrium.Elsewhere we have presented a derivation that parallels closely the kinetic theory analysis[15].This approach seems feasible only in the high-temperature regime.In order to analyse the low-temperature dynamics we have reformulated the theory in the spirit of Mori’s[16]continuous fraction expansion.Here we present a derivation of the basic equations and an application to the high-temperature phase.In the following paper[17]we analyse dynamics in the low-temperature phase.2.Theory2.1.DeﬁnitionsThe hard-spin SK model consists of Ising spinsσi=±1interacting via inﬁnite-range exchange coupling constants J ij,H=−i<jJ ijσiσj.(1) The coupling constants J ij are quenched random variables distributed according to thesymmetric distribution P(J ij)∼exp(−J2ij /(2J2/N)).The time-dependent spin–spin correlations can be written in the following form:δσi(t)δσj eq= δσi exp( t)δσj eq.(2)Hereδσi is theﬂuctuation of the value of the i th local spin,δσi=σi− σi eq and ··· eq denotes the equilibrium ensemble average with the Boltzmann distribution,··· eq≡σ1,σ2,......P eq P eq∼exp(−βH).(3) Finally, is the evolution operator=−i(1−S i)w i(4) with S i being the spin–ﬂip operator,S iσi=−σi,and w i being the transition ratew i=(1−σi tanh(βh i))/2(5) where h i is a local magneticﬁeld acting on the i th spinh i=j=iJ ijσj.(6)Note that in equation(2)the equilibrium probability distribution P eq stands to the right of the quantity being averaged and the evolution operator exp( t)acts on everything to its right(including the probability distribution).Glauber dynamics of the SK model100472.2.Zeroth-order approximation To calculate the time-dependent spin–spin correlation function we follow the spirit of Mori’s[16]continuous fraction expansion.In the zeroth-order approximation we restrict ourselves to the subspace spanned by the spin ﬂuctuations.We deﬁne a projection operator on this subspace P 0f (σ1,σ2,...)=ijδσi A ij δσj f (σ1,σ2,...) eq .(7)Here the matrix A is the inverse matrix of spin correlations that is deﬁned by the following equation: jA ij δσj δσk eq =δik .(8)Note that A is equal to the Hessian of the TAP [18]free energy.At and above T c A is known explicitly [19,2]:A ij =−βJ ij +δij (1+(βJ )2).(9)We start by writing down a formal expression for the time derivative of the spin–spin correlations∂t δσi (t)δσj eq = δσi exp ( t)δσj eq .(10)Next,we insert into it the identity operator written as a sum of the projection P 0and the orthogonal projection Q 0=I −P 0δσi exp ( t)δσj eq = δσi (P 0+Q 0)exp ( t)δσj eq .(11)In the zeroth-order approximation we neglect the contribution involving Q 0.Then,using the deﬁnition of the projection operator P 0we get ∂t δσi (t)δσj eq =klδσi δσk eq A kl δσl (t)δσj eq .(12)Finally,using the explicit form of the evolution operator we rewrite the above result in the following way:∂t δσi (t)δσj eq =−(1− σi tanh (βh i ) eq ) lA il δσl (t)δσj eq .(13)This equation is equivalent to a disorder-dependent version of the local equilibrium approximation of Kawasaki [20].It is identical to the evolution equation obtained from the zeroth-order kinetic theory [15].One sees immediately that the evolution equation (12)with the Hessian (9)is almost identical to that derived in the original SK paper [1].It should be emphasized,however,that a dynamical version of the Onsager reaction ﬁeld term has been naturally included in the matrix A .Straightforward calculation shows that the sample-averaged solution of (12)has the same form as the sample-averaged solution of the SK equation if the timescale of SK is rescaled by a factorτ=1/(1− tanh (βh)σ eq ).(14)Here the overline denotes the ‘spatial average’:f =(1/N) i f i .Explicitly,we ﬁnd [ δσi (t)δσi (0) eq ]=2π 1−1d x (1−x 2)1/21+(βJ )2−2βJ x exp −(1+(βJ )2−2βJ x)t τ .(15)10048G SzamelHere [...]denotes sample averaging over the probability distribution P (J ij ).Note that in order to get the above result we have used the usual approximation that the eigenvectors of the Hessian are uncorrelated.As discussed by SK [1]the formula (15)leads to the algebraic t −1/2decay of the spin correlations at the transition temperature T c =J .Note that the prefactor of the t −1/2term is different from that predicted by SK.The factor (14)can be calculated using the probability distribution of the local ﬁelds calculated in [21].Explicitly,we getτ−1=1−[ tanh 2(βh) eq ]=1−1J √2π ∞−∞d h tanh 2(βh)cosh (βh)exp (−β2J 2/2−h 2/2J 2).(16)Note that here we have implicitly used the usual approximation of replacing the spatial average by the sample average.We will compare the theoretical prediction (15)with the simulations in section 4.Here we only note that at T =T c numerical integration gives τ≈2.2242which roughly corresponds to the discrepancy between the original SK result (equation (15)with τ=1or equation (5.27)of [1])and the simulation data they presented (ﬁgure 13(c )of [1]).Finally,it should be noticed here that the zeroth-order approximation is exact at short times.More precisely,equation (13)reproduces exactly the ﬁrst time-derivative of the spin correlations at t =0.This is of course analogous to the simple CS theory being exact at short times.2.3.First-order approximationIn order to go beyond the zeroth-order theory we ﬁrst use the standard projection operator manipulations [22]in order to rewrite the exact evolution equation (11)in the following way: t 0 j(δij δ(t −t )+M irr ij (t −t ) 1−σj tanh (βh j ) −1eq )∂t δσj (t )δσk eq =− 1−σi tanh (βh i ) eql A il δσl (t)δσj eq .(17)Here M irr is the irreducible [22]memory matrixM irr ij (t)= D i exp (Q 0 irr Q 0t)D j eq (18)with D i deﬁned asD i =−Q 0(σi −tanh (βh i ))(19)and irr is the irreducible part of the evolution operator, irr = − ij δσi eq [ δσi δσj eq ]−1 δσj .(20)D i can be interpreted as the part of the time derivative of the i th spin ﬂuctuation that is orthogonal to all the spin ﬂuctuations.To this end we have to deﬁne a conjugate evolution operator ∗∗=− iw i (1−S i )(21)Glauber dynamics of the SK model10049 that acts on the spin variables,and then we check that the time derivative ofδσi is given by the following expression:∂tδσi= ∗δσi=−(σi−tanh(βh i)).(22)It should be emphasized here that the formal manipulations leading to the evolution equation(17)have a simple physical interpretation[22]:the irreducible memory matrix renormalizes the inverse rate of spinﬂips.In the zeroth-order approximation we neglected the renormalization of the inverse rate completely.In the next step(ﬁrst-order approximation)we include it approximately.To this end weﬁrst deﬁne a projection operator on the subspace spanned by D iP1f(σ1,σ2,...)=D i C ij D j f(σ1,σ2,...) eq.(23)ijHere C is the inverse matrix of the correlations of the projected time derivativesC ijD j D k eq=δik.(24)jNext,we write down a formal expression for the time derivative of the memory function and insert into it the identity operator written as a sum of the projection P1and the orthogonal projection Q1=I−P1∂t M irr ij(t)= D i Q0 irr Q0exp(Q0 irr Q0t)D j eq= D i Q0 irr Q0(P1+Q1)exp(Q0 irr Q0t)D j eq.(25) In theﬁrst-order approximation we neglect the contribution involving Q1∂t M irr ij(t)≈ D i Q0 irr Q0P1exp(Q0 irr Q0t)D j eq.(26)Then,using the deﬁnition of the projection operator P1,we getD i Q0 irr Q0D k eq C kl M irr kj(t).(27)∂t M irr ij(t)=klAt this point we note that in the high-temperature limit the off-diagonal equilibrium correlations between irreducible quantities(i.e.quantities that have spin correlations subtracted out)are of higher order in N−1/2than the spin–spin correlations.This can be easily shown by a high-temperature expansion.Moreover,it can also be shown that the off-diagonal part of the memory matrix evolution operator D i Q0 irr Q0D j eq,i=j,can be neglected in the thermodynamic limit.Hence, only the diagonal partD i Q0 irr Q0D i eq=−(βJ)2 1−tanh2(βh) eq( tanh4(βh i) eq− tanh2(βh i) 2eq)(28) is relevant.Finally we see that in theﬁrst-order approximation only the diagonal elements of the memory matrix contribute and,moreover,that their time evolution is given by∂t M irr ii(t)= D i Q0 irr Q0D i eq C ii M irr ii(t)(29) where= 1−tanh2(βh i) eq(1−A ii 1−tanh2(βh i) eq).(30) M irr ii(t=0)=C−1iiEquations(17)and(29)with(28)and(30)constitute ourﬁrst-order approximation. In the high-temperature regime it is actually possible to go beyond theﬁrst-order approximation.The structure of the theory remains the same:at any higher order level a new(diagonal)irreducible memory matrix appears.10050G SzamelIt should be noted here that one should not expect any divergence of the memory function relaxation time at the transition temperature: D i Q0 irr Q0D i eq stays positive and C ii does not approach0as T→T c.To get explicit results from theﬁrst-order approximation we have to solve a system of equations(17)and(29).In the high-temperature phase we pre-average equation(29)to get ∂t M irr(t)= D Q0 irr Q0D eq C M irr(t)(31) with(t=0)= 1−tanh2(βh) eq(1−(1+(βJ)2) 1−tanh2(βh) eq)(32) and then we use the pre-averaged irreducible memory function in equation(17).The main argument for the pre-averaging is analytical simplicity.However,we expect that in the high-temperature phase pre-averaging should be justiﬁable.Equations(17)and(31)can be easily solved.Weﬁnd that theﬁrst-order approximation results in a minute correction of the zeroth-order one.This fact could have been anticipated: numerical integration using the localﬁeld distribution of[21]shows that at the transition temperature M irr(t=0)≈0.045.puter simulationsWe performed a series of Glauber dynamics simulations of the SK model at temperatures 2T c,1.5T c,1.25T c,and T c using the algorithm of Mackenzie and Young[23].Above the transition temperature we used equilibration times varying between100Monte Carlo steps per spin(MCS)at2T c and200MCS at1.25T c.At the transition temperature we used very long equilibration time of10000MCS.We used sample sizes of N=3000:at each temperature we simulated100samples.After equilibrating the system we collected the data for the time-dependent spin–spin correlation function[ σi(t)σi(0) eq]=(1/N)iσi(t)σi(0) eq(33) and for its time derivative,i.e.for the following correlation function:d d t [ σi(t)σi(0) eq]=(1/N)i(tanh(βh i(t))−σi(t))σi(0) eq.(34)We monitored the derivative independently having in mind the low-temperature phase:in contrast to the correlation function the derivative decays to zero even in the low-temperature phase.4.Results and discussionInﬁgure1we compare time-dependent spin correlations calculated from theﬁrst-order theory with the results of the computer simulations.On the scale of theﬁgure the predictions of the zeroth-order theory are indistinguishable from those of theﬁrst order.Inﬁgure2we compare the time derivative of the spin correlations calculated from the ﬁrst-order theory with the results of the computer simulations.Again,on the scale of the ﬁgure the predictions of the zeroth-order theory are indistinguishable from those of theﬁrst order.Glauber dynamics of the SK model10051Figure1.Time-dependent spin–spin correlation function as a function of time.Symbols:Glauber dynamics simulation data;squares:T=2T c;circles:T=1.5T c;triangles:T=1.25T c;crosses:T=T c.Solid lines:predictions of theﬁrst-order theory.Figure2.Time derivative of the spin–spin correlation function.Symbols:Glauber dynamicssimulation data;squares:T=2T c;circles:T=1.5T c;triangles:T=1.25T c;crosses:T=T c.Solid lines:predictions of theﬁrst-order theory.Qualitatively and quantitatively there is perfect agreement between the theory and simulations.It should be emphasized here that it is theﬁrst time such an agreement has been achieved for the time-dependent spin correlations,even in the high-temperature phase.It should also be realized that in spite of this quantitative agreement theﬁrst-order theory is not exact.There is an inﬁnite number of higher-order irreducible memory functions and we do not have any argument that would allow us to neglect them.In this regard it seems that,even in the high-temperature phase,dynamics of the SK model is much more complicated than its equilibrium properties.Finally,we would like to comment on the relation between our theory and the simple10052G SzamelCS approach[10].CS,in effect,use a local in time closure hypothesis in the evolution equations for the sample-averaged quantities.If we were to derive an equation of motion for the sample averaged spin–spin correlation function from the zeroth-order approximation (equation(13)),our closure approximation would be non-local in time.Hence,we do not disagree with the main CS assumption that(certain)dynamic quantities are self-averaging and therefore it should be possible to derive equations of motion that would involve sample-averaged quantities only.On the other hand,we suspect that these equations of motion might be non-local in time.One should note at this point that the equations of motion for the time-dependent correlations that are derived in the soft-spin approach are,in fact,non-local in time.AcknowledgmentsI would like to thank Marshall Fixman for stimulating discussions.This work was partially supported by NSF grant no CHE-9624596.References[1]Sherrington D and Kirkpatrick S1975Phys.Rev.Lett.351972Kirkpatrick S and Sherrington D1978Phys.Rev.B174384[2]Binder K and Young A P1986Rev.Mod.Phys.58801[3]Fisher K H and Hertz J A1991Spin Glasses(Cambridge:Cambridge University Press)[4]Sompolinsky H and Zippelius A1981Phys.Rev.Lett.47359Sompolinsky H and Zippelius A1982Phys.Rev.B256860[5]Sompolinsky H1981Phys.Rev.Lett.47935[6]See e.g.Cugliandolo L F and Kurchan J1994J.Phys.A:Math.Gen.275749[7]Kinzel W and Fisher K H1977Solid State Commun.23687Fisher K H1983Solid State Commun.46309[8]Sommers H J1987Phys.Rev.Lett.581268[9] L usakowski A1991Phys.Rev.Lett.662543[10]Coolen A C C and Sherrington D1993Phys.Rev.Lett.713886Coolen A C C and Sherrington D1994J.Phys.A:Math.Gen.277687[11]Laughton S N,Coolen A C C and Sherrington D1996J.Phys.A:Math.Gen.29763[12]de Almeida J R L and Thouless D J1978J.Phys.A:Math.Gen.11983[13]Nishimori H and Yamana M1996J.Phys.Soc.Japan653See also Yamana M,Nishimori H,Kadowaki T and Sherrington D1997J.Phys.Soc.Japan661962 [14]For a formulation of the Enskog kinetic theory that resembles the simple CS theory see Resibois P1978J.Stat.Phys.19593General discussion of the Enskog theory can be found in R´e sibois P and de Leener M1977Classical Kinetic Theory of Fluids(New York:Wiley)[15]Szamel G1997J.Phys.A:Math.Gen.305727[16]Mori H1965Prog.Theor.Phys.34399[17]Szamel G1998J.Phys.A:Math.Gen.3110053[18]Thouless D J,Anderson P W and Palmer R G1977Phil.Mag.35593[19]Bray A J and Moore M A1979J.Phys.C:Solid State Phys.12L441[20]Kawasaki K1966Phys.Rev.145224Kawasaki K1966Phys.Rev.148375[21]Thomsen M,Thorpe M F,Choy T C,Sherrington D and Sommers H J1986Phys.Rev.B331931[22]Kawasaki K1995Physica A21561[23]Mackenzie N D and Young A P1983J.Phys.C:Solid State Phys.165321。

simultaneous three-equation model

simultaneous three-equation model1. 引言1.1 概述在经济学和社会科学研究中，同时考虑多个相关方程的模型被广泛应用。

这些模型被称为"simultaneous three-equation models"，因为它们包含三个相互关联的方程。

这些方程通常代表着不同的经济变量之间的关系，通过分析这些关系可以更好地理解经济系统的运作机制。

1.2 文章结构本文将对simultaneous three-equation model进行详细的探讨和分析。

首先，我们将介绍该模型的理论背景，包括概述、模型假设以及应用领域。

接下来，我们将阐述模型的建立过程和参数估计方法，并介绍模型选择准则和诊断检验。

然后，我们将展示实证分析结果，并对其进行解释和讨论。

最后，我们将总结主要研究发现并展望存在的问题与改进研究方向。

1.3 目的本文旨在提供一个全面而深入的介绍simultaneous three-equation model，并通过实证分析来展示其在经济学研究中的应用价值。

我们希望通过对该模型的探讨和分析，能够增加读者对经济系统运作机制的理解，并为相关领域的研究和政策制定提供有益的参考。

同时，我们也将指出该模型存在的问题，并展望未来改进研究的方向。

2. 理论背景2.1 三方程模型概述在经济学和社会科学领域，三方程模型是一种常用的统计分析方法。

该模型通过建立三个联立方程来描述多个变量之间的关系和相互作用。

这些方程通常代表着不同变量之间的因果关系，可以揭示出导致某一现象或结果的因素。

三方程模型通常包括三个主要成分：自变量、因变量和错误项。

自变量是研究者选定的与因变量有关的解释变量，它们对因变量产生直接或间接影响。

因变量是研究者希望解释或预测的目标变量。

错误项则代表了未被观察到或遗漏的其他影响因素，即未能被自变量完全解释并捕捉到的部分。

2.2 模型假设为了建立一个有效的三方程模型，需要满足一些基本假设。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。