Cmp Mod Fluid Dynm
Cmp Mod Fluid Dynm MPO 662
Popular in Course
Popular in Meteorology Phys Oceanography (Mpo)
verified elite notetaker
verified elite notetaker
BLAW 3310 - 001
verified elite notetaker
verified elite notetaker
verified elite notetaker
This 38 page Class Notes was uploaded by Kip Turner on Thursday September 17, 2015. The Class Notes belongs to MPO 662 at University of Miami taught by Staff in Fall. Since its upload, it has received 57 views. For similar materials see /class/205735/mpo-662-university-of-miami in Meteorology Phys Oceanography (Mpo) at University of Miami.
Reviews for Cmp Mod Fluid Dynm
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 09/17/15
JOURNAL OF COMPUTATIONAL PYSICS 54 325 362 1984 A Fully Multidimensional Positive Definite Advection Transport Algorithm with Small Implicit Diffusion PIOTR K SMOLARKIEWICZ National Center for Atmospheric Research Boulder Colorado 80307 Received July 8 1983 The idea of the simple positive de nite advection scheme presented previously in Monthly Weather Review 11 1983 479 is improved for an optional multidimensional case and is presented in a generalized format The accuracy of the scheme is discussed and a review of existing options is presented and illustrated through numerical tests 1 INTRODUCTION In numerical modeling of physical phenomena it is often necessary to solve the advective transport equation for positive de nite scalar functions Numerical schemes of second or higherorder accuracy can produce negative values in the solution due to the dispersive ripples Lower order schemes such as the donor cell or Lax Friedrichs or higher order schemes with zerothorder diffusion added produce no ripples but suffer from excessive implicit diffusion In the last ten years a possible resolution of this dilemma has been developed in the form of hybrid schemes in which the advective fluxes are given as a weighted average of a rst order positive de nite scheme s uxes and a higher order scheme s fluxes The difference in deter mination of the weights in the calculation of the average advective fluxes has led to different hybrid schemes Two main hybrid type schemes have been developed One the so called uxcorrected transport FCT method was originated by Boris and Book 1 3 and generalized by Zalesak 14 the other was developed by Harten and Zwas 6 9 in the form of the self adjusting hybrid schemes SAHS method Both methods were constructed to deal effectively with shocks and contact discontinuities Solutions of the advection transport equation obtained by using FCT or SAHS maintain positive de niteness of the initial condition and as can be seen from presented tests Zalesak 14 Harten 6 be very accurate Unfortunately application of these methods to the modeling of complex multidimensional hydrodynamical systems like atmospheric phenomena is rather limited due to the excessive computer time required Furthermore in many hydrodynamical systems On leave from the Institute of Geophysics Warsaw University Poland The National Center for Atmospheric Research is sponsored by the National Science Foundation 325 0021999184 300 Copyright 1984 by Academic Press Inc All rights of reproduction in any form reserved 326 PIOTR K SMOLARKIEWICZ especially in the presence of turbulent diffusion dealing with a shock or discontinuity is not as important as maintaining the positive de niteness of the evaluated scalar quantity Compromising between computer time ef ciency and accuracy Clark and Hall 45 applied the idea of hybridization to the numerical evaluation of the cloud s water transport equation developing possibly the simplest known hybrid scheme The numerical diffusion in this scheme is larger than in the FCT but the time consumption is about half of that in the FCT Recently developed by Harten 7 8 the totalvariationdiminishing TVD schemes like FCT or SAHS have been designed to treat shocks and discontinuities The explicit version of the second order accurate TVD scheme 8 is computationally ef cient enough to be considered for application in modeling the phenomena mentioned above in Section 4 of this paper one version of the TVD scheme will be discussed in detail All FCT SAHS and TVD methods have been constructed to maintain monotonicity of the initial condition which implies the maintenance of the positive de niteness of the transported quantity This paper presents another solution Using an iterative approach one can construct from the basis of the upstream scheme a class of nonlinear multidimen sional positive de nite advection transport algorithms with small implicit diffusion In comparison to the FCT SAHS or TVD schemes the general iterative principle of the algorithm seems to be simple and can be easily developed from the Taylor series expansion applied to the upstream scheme 13 In the rst iteration an upstream scheme is used in its classical sense while each following corrective iteration reapplies the upstream scheme but with a specially defined antidiffusive velocity eld The number of iterations is optional and each additional iteration increases the solution s accuracy Such a procedure has been applied successfully in Smolarkiewicz 13 to onedimensional advection problems and to a multidimensional problem when the scheme has been used in a time splitting form In a multidimensional case when the scheme was applied to a combined form the difference between the combined and timesplitting forms of the advection schemes was discussed in detail in Smolarkiewicz 12 only the rst corrective iteration gave improved results The second corrective iteration unrealistically deformed the solution and in some cases even resulted in instability of the scheme In 13 it was suggested that this error was caused by an effect due to the crossderivative terms as was also discussed in detail in 12 Introducing cross terms to the antidiffusive velocities eliminates this problem and results in a fully multidimensional uni ed algorithm that strictly maintains the positive de niteness of the transported quantity The proposed algorithm represents an open family of schemes of varying levels of accuracy complication and computational ef ciency The simplest version of the scheme is secondorder accurate in both space and time As was shown in 13 it gives results of a quality comparable with that obtained from more complicated hybrid schemes while considerably reducting computational time The most accurate and complicated version presented in this paper is third order accurate in time and fifthorder accurate in space In principle it is possible to construct the algorithm with an optional order of accuracy For smooth initial conditions the schemes also preserve monotonicity POSITIVE DEFINITE ADVECTION SCHEME 327 but for the shocktype initial conditions they may produce locally in the closest neighbourhood of the shock a small ampli cation This ampli cation is signi cantly smaller than the amplitude of oscillation produced by linear schemes of a high order of accuracy and is of minor importance in a broad category of problems However the example of the modi cation of the algorithm for the case of a shocktype initial condition will be discussed as one of the possible options of the scheme In Section 2 of this paper the basic form of the algorithm will be developed in general optionally dimensional form In Section3 such features of the scheme as consistency stability and accuracy will be discussed Section4 will present the results of some tests and compare the scheme with a version of the TVD scheme because the rst version of the algorithm was compared with the FCT SAHS and Clark and Hall schemes in 13 this comparison will not be repeated here In Section 5 some of the possible options of the scheme will be discussed 2 DEVELOPMENT OF THE ALGORITHM The equation to be solved is the continuity equation describing transport of the nondiffusive scalar quantity in Mdimensional space the proposed algorithm was placed in a general M optional format for compactness but the author is able to prove the stability of the scheme for M g 3 only which does not necessarily mean that the scheme is unstable in a case of M gt 3 811 M 8 E 1 8t 18x wu O where WE ultx39x quot is the nondiffusive scalar quantity assumed to be non negative u E u t x x is the Ith velocity component I l M t x E x x are the time and spaceindependent variables To describe compactly the numerical equations that will be used later it is convenient to introduce a few symbols 14quot is a numerical approximation of the solutin of Eq 1 de ned in points tquot xi where tquot n At xi 1quot AX39 i2 AX2 i AX n 0 NT 1quot 0 NX and AX is the constant spatial increment in the Ith direction in this paper the indices described by capital letters eg I J always indicate vector components and the indices described by lowercase letters eg i j indicate the position on a grid space e E 0 0 0 1 0 0 is a unity vector in the Ith direction 328 PIOTR K SMOLARKIEWICZ quotui Hlm is the I th velocity component in the nth time step de ned on a staggered grid the I th component is staggered AX in the I th direction Later to simplify the fomulas the index n beside velocity symbols will be omitted where it is possible to do so without introducing any misunderstandings With the above symbols the well known rst orderaccuracy upstream scheme for numerical evaluation of Eq 1 can be written in the form M Win1 Win E F M39a Wine5 uil2e FIWiLe Viquot Lu2 2 11 where F is the advective flux in I th direction evaluated in the same staggered points as the Ith velocity component and de ned as follows mi u u lul wi u lul w 3 A suf cient stability condition for the scheme in Eq 2 can be written in the form M V 2 laiI2gtel Cum2mg g lt1 4 t xiuzel 1 cf Eqs 3 140 of Roache 11 where ai12eE ui12eAtAXI 5 Note that originaly u was not de ned in xHUzm for 1 9E J and that this de nition is not necessary for numerical evaluation of Eq 1 For the stability condition it always can be done at least by the appropriate arithmetical averaging see eg Eq 14 The constant lt l was placed in Eq 4 to point out that in the case of a divergent flow eld the traditional restriction can be insuf cient and as follows from experience and simple analysis of the behavior of the upstream scheme around the point in which each velocity component changes from u to u should be replaced by 9 lt Under the condition given in 4 the scheme in 2 is positive de nite which means lt w9gt0 v W20 6 I lquotxi The properties in Eqs 4 and 6 of the scheme in Eq 2 as well as its low computer time consumption are very useful for the application of Eq 2 to the numerical evaluation of Eq 1 Unfortunately the scheme in Eq 2 is of rstorder accuracy in both time and space and has strong implicit diffusion The rate of the implicit diffusion may be easily estimated for the case of a uniform flow POSITIVE DEFINITE ADVECTION SCHEME 329 V I u const Expanding g4 1 whet ML in a secondorder Taylor sum about the point tquot xi scheme 2 may be written as all n M a I n M a 05 AX At at i 121 8x1 u i 1 8x iu i u ax I J 121 05 Atu u ax i 7 lat From Eq 7 it can be seen that when At and AX a O for all I Eq 7 approaches Eq 1 but during a realistic computational process the scheme in Eq 2 with nite At and AX approximates more accurately with secondorder accuracy an advective transport equation with additional diffusive terms rather than the original Eq 1 On the other hand these implicit diffusion terms are important for the stability of the scheme and must not be explicitly subtracted from the scheme An intuitively obvious approach is to make the advection step using Eq 2 and then reverse the effect of the diffusion equation 8 M a l I 1 2 8394 M 1 J 8W 797 ng 05u IAX Atu 1W J 05Atuu W 8 1 in the next corrective step The diffusion process and the equation that describes it are irreversible But this does not mean that the effect of the diffusion process cannot be reversed in time To do this one can design a process that will return to the previous state or at least to a state close to it It is enough to notice that introduction of the arti cial diffusive velocities uf allows one to write Eq 8 in the form 8w M 8 aT Z 37013 W 9 1 where 1 8w 05 I AX At I 2 iui u l W ax M 1 6 us J 05Atu uJVa JZIJ if Wgt0 10 1 0 if w0 Now de ning antidiffusive velocities E ufj 11 330 PIOTR K SMOLARKIEWICZ the reversal in time of the diffusion equation 8 may be simulated by the advection equation 9 with the antidiffusive velocities 17 instead of u Based on this concept the corrective step is suggested in the form M W w Z FRIIiquot17iH12F wiaviii um 12 I l where Ili is given by the righthand side rhs of Eq 2 and Vii IZII12eE luil12eiAXl Atuiu22 X We We 8 AX ie i l J 05 At uil2euil2e l ZEN Wite1e1 with Vitee lfliaie Wit ee Wiie witeI e Wife 6 AXJ 13 where J J J J J uil2e 025uie12e Home uic l2cJ ui 12eJ 14 and e is a small value eg 10 to ensure 12quot 0 when W eI w 0 or wit e19 wit wield wife 0 The corrective step is again the upstream scheme and also contains implicit diffusion that again can be corrected by the next corrective step The number of corrective steps is optional so the final algorithm may be written in the form M k Wiwk Wiwkil 2 FI39ikl Wiflf l Lil2 k FIWijltek l Wimk l Lu2n 15 where k IORD W 5011 w evaluated from Eqs 13 and 14 and l w E Mquot y quot D 5 W uleZME quotaha2 Note that when IORD l the algorithm in Eq 15 is exactly the same as the classical upstream scheme Theoretically IORD may be any value but as can be concluded from the performed tests using IORD gt 4 only negligibly improves the accuracy of the solution while increasing the computational costs of the scheme In 13 the second term in Eq 13 which compensates the effect of the cross spatial derivatives in Eq 7 was not evaluated In a multidimensional case for IORD gt 2 this neglect resulted in deformation of the solution or even in instability of the scheme POSITIVE DEFINITE ADVECTION SCHEME 331 3 CONSISTENCY STABILITY AND ACCURACY OF THE ALGORITHM The proposed algorithm is constructed from the well known conditionally stable and consistent upstream scheme To show the consistency of the entire scheme in Eq 15 it is enough to show that when At AX gt 0 for all I the corrective iterations k gt 1 in Eq 15 do not affect the solution of the rst iteration k l in Eq 15 From Eqs 13 11 and 10 it is easy to see that when At AXi gt 0 then wk diamm 0 which implies ugwmiao for all k The latter implies that Eq 15 for k gt 1 converges to Q W t0 16 which means that the scheme is consistent To show the stability of the scheme it is enough to show that stability of the rst iteration implies stability of all subsequent iterations An optional grid point and time step may be taken into consideration to prove stability for any arbitrary grid point and time step The notation is simpli ed for this analysis by omitting the indices that indicate the position on the grid space The application of Eq 13 results in quota 5 la l G39Y a a 17 xx Mt where 07 is related to d as in Eq 5 and 6 and 6 are the Aww ratios in Eql3 Because all w values have been obtained from the positive de nite upstream scheme VJ l lt 1 Because of Eq 4 Va ltlt 1 which implies 0 g a 1 2 lt025 in Eq 17 Using this information it can be shown that Id39ltta39I a 2la39z f la l 18 iii Because of Eq 4 M 2 a C 19 l J The sum of Eq 18 over all I and the use of Eq 19 lead to M id ltCCZ zla 2 20 11 M I l 5815429 332 PIOTR K SMOLARKIEWICZ From the Schwartz s inequality M E la392gtCZM 21 I l which allows us to write Eq 20 in the form M 2 WI ltC C2M 3M 22 11 The rhs of the inequality in 22 is less than or equal to C for M g 3 so nally M M V 21 di12el lt 2 law 23 tquotxiu2eJMlt3 1 1 which results in M k1 i4 k I aham lt l aha2 24 quotXilt12eJMlt 3 1 Additionally the maximum values of the rhs of Eq 22 are and C for M l 2 and 3 respectively This is an important conclusion because even if the original velocity eld is nondivergent the antidiffusive velocity eld is divergent So to be positive that the whole scheme is stable it is safer to ensure that te rhs of Eq 24 is less than or equal to Finally we have proved that for M l 2 the stability of the rst iteration implies the stability of the whole scheme and for M 3 this is true if the constant in Eq 4 is less than or equal to The performed tests eg M 3 and gO95 suggest that the algorithm is always stable if the original upstream scheme is stable but this conclusion is still unproved for M 2 3 Because each iteration of the algorithm is in the form of the upstream scheme ful lling the stability criteria ensures that the scheme is positive de nite The scheme is nonlinear even in the case of a uniform velocity eld According to Harten 8 p 6 stability and consistency are the only necessary conditions for its convergence From Eqs 12 and 13 it is easy to see that when At AX gt 0 for all I the scheme converges to an upstream scheme The latter point as well as the results of the performed tests suggests that convergence of the scheme is at least as well ensured as for the upstream scheme The algorithm s order of accuracy can be determined by estimating the truncation error It can be shown and con rmed by the results of tests that the scheme is at least secondorderaccurate in space When time dependence of the velocity eld can be ignored the scheme is also second orderaccurate in time Let s assume for the sake of simplicity that the considered scheme has IORD 2 and that the speci ed velocity eld is time independent and nondivergent Although for development of Eq 7 a uniform velocity eld has been assumed it is easy to show that to obtain Eq 7 it is enough to assume a nondivergent velocity eld in Section 5 of this paper it will be t 02AXAtAt 29 334 PIOTR K SMOLARKlEWICZ Because the third term on the rhs of 29 is of the leading order AX 2 At Eq 29 may be nally written in the form tquot At 3 Wm aquot J E W dt 02AXAt At 30 n which means that 911quot is a solution of Eq 1 that is second order accurate in time and space The same conclusion could be reached in a more heuristic manner It is enough to note that the corrective step compensates the rstorder small diffusive term with the accuracy to the rst order of the rst order small term This means that the residual of the compensation is a secondorder small term It is still unclear whether increasing the number of iterations increases the order of accuracy of the scheme or only decreases the amplitude of the existing 02AX At At error in 30 From Eqs 13 and 15 it can be concluded that a two iteration scheme IORD 2 results in a three point scheme while a three iteration scheme YiLBC TRUNC ERRER 2 400 I 1 l l l 00 10 20 30 40 50 60 70 NUMBER F DTDX HRLMINGS FIG 1 The dependence of the measure of the truncation error on the number of halvings of the grid and time increments The curves Sl S8 represent schemes with 0RD 1 8 respectively POSITIVE DEFINITE ADVECTION SCHEME 335 IORD 3 results in a vepoint scheme Generally an IORDiteration scheme determines gradients of uxes on the basis of the information from 2 IORD l points To illustrate the behavior of the truncation error the following test has been performed In a uniform one dimensional velocity eld the initial condition in the form of the Gausstype function has been assumed After an arbitrarily chosen fixed time period T NT At the average error per time and per space step between the numerical and analytical solutions has been evaluated as follows 12 TRERNTNX a wTx WTYNT NX 31 where wTx39 w i are the analytical and numerical solutions respectively at the point T xi Dividing successively At AX by 2 and doubling simultaneously NT NX the sequence of TRERNTNX has been obtained The dependence of the log2TRER on the number of halvings of the original At AX is presented in Fig l for different IORD s and the chosen Courant number 05 The curves S1 52 S4 S6 YILQGZITRUNC ERRQR 25i0 400 l l I l l 00 10 20 30 40 50 60 70 NUMBER BF DTDX HRLVINGS FIG 2 As in Fig 1 but for a leapfrogtrapezoidal scheme secondorder accurate in time and second L2 to eight L8 order accurate in space 336 PIOTR K SMOLARKIEWICZ and S8 are equivalent to IORD 1 2 4 6 and 8 respectively In Fig 2 the same dependence is presented but for the well known leapfrog trapezoidal scheme Zalesak 14 Appendix which is secondorder accurate in time and second fourth sixth and eighth orderaccurate in space The curves L2 L4 L6 and L8 represent the different versions of the leapfrog scheme Note that the quantity TRER de ned in Eq 31 does not represent the truncation error according to its de nition 10 p 20 but rather some measure of that error So the order of the accuracy of the various schemes cannot be estimated directly from the plots presented in Figs and 2 Generally the truncation error may be written in the form truncation error Atquot39Cl C2AXquot quot 32 where p px are the leading orders of the truncation in time and space respectively and C1 C2 are coef cients that are generally solution and Courant number YZLBGZ 1 TRUNC ERRBR 50 l I l I I 41 J J 1 1 l 00 10 20 30 40 50 60 70 NUMBER 0F DTDX HFILVINGS 4010 FIG 3 As in Fig 1 but for an IORD 2 scheme Curves C1 C2 and C3 are for CFL 025 05 and 075 respectively POSITIVE DEFINITE ADVECTION SCHEME 337 dependent From 32 it can be seen that when At 2 quot Ato AX 2 quot AXO and n is big then log2C1C2 if pxpt lo tr error n g2 p39 log2 C1 if px gt p 33 The relation in 33 explains the shape of the curves in Figs 1 and 2 Note that for a big n all lines L2 L8 have the same slope and the lines L4 L8 are identical Generally it can be concluded that the slopes of the lines in Figs 1 and 2 illustrate the order of truncation error in time while the shift between the lines for the same type of scheme is related to the order of the truncation in space Comparing S2 with L2 and S4 S8 with L4 L8 in Figs 1 and 2 it can be seen that S2 has the same slope as L2 and lies entirely below L2 while S4 S8 have greater slopes than L4 L8 and lie entirely below them To make a final conclusion it is necessary to present the dependence of the TRER quantity on the Courant number In Figs 3 and 4 the two YILEJG iTRUNC ERRDR 2 5 0 I l f f I I 7 3 J 4100 T F 15 0 Tc E C3 200 x C 1 Cix r 250 a t 300 e 35 o t 400 I L l l l l 00 10 20 30 40 50 60 70 NUMBER BF DTDX HRLVINGS FIG 4 As in Fig 3 but for a leapfrogtrapezoidal scheme secondorderaccurate in both time and space 338 PIOTR K SMOLARKIEWICZ iteration scheme and second order in the time and space leapfrog scheme are compared for Courant numbers Cl 025 C2 05 and C3 075 From Figs 3 and 4 it can be seen that for all Courant numbers the lines for the two iteration scheme have the same slopes as the lines for the leapfrog scheme It also can be seen that the curves in Fig 3 lie entirely below the curves in Fig 4 These features suggest that the twoiteration scheme is secondorder accurate in time and space In Figs 5 and 6 the same comparison is presented for a four iteration scheme and a second order accuracy in time and fourthorder accuracy in space leapfrog scheme Comparing Figs 1 2 and 5 it can be concluded that the greater slope for S4 SS than for L4 L8 does not indicate a difference in order of accuracy in time between say S4 and L4 but rather some particular dependence of the truncation error on the Courant number Comparing Figs 5 and 6 it can be noticed that the values of TRER for the fouriteration scheme are smaller than those for the 24 leapfrog scheme but have a different dependence on the Courant number Note that C1 in Fig 5 lies above C but below C3 in Fig 6 This suggests that the fouriteration scheme is fourth YILBBZITRUNC ERRBR 750 T7 1 i l l 10 0 DP N 1 F ai 715 0 7K 1 ecx i t c 20 0 j 09 2 C g 250 e lt4 7 Ce l l 4 l g gt 10 20 30 40 50 00 70 NUMBER BF DTDX HRLVINGS g g 390 39o and 3 FIG 5 As in Fig 3 but for an IORD 4 scheme POSITIVE DEFINITE ADVECTION SCHEME 339 YILUO TRUNC ERRBR 2 50 I I I I I I E 10o gt3 clf03 150 V C9 03 C 2oro Ce C 3 C a x 09 2so Cu 300 quot 350 39 400 I 1 l 1 1 00 10 20 30 40 50 6 0 70 NUMBER F DTDX HRLVINGS F 10 6 As in Fig 3 but for a leapfrogtrapezoidal scheme secondorder accurate in time and fourth orderaccurate in space orderaccurate in space and generally that increasing the number of iterations increases the order of accuracy in space It also can be concluded that the truncation error reaches a minimum for Courant number 05 4 RESULTS OF TESTS Two and three dimensional solidbody rotation tests were performed to illustrate the behavior of the proposed algorithm The onedimensional version of the scheme was tested in 13 The twodimensional case used 101 points in each direction with AX 1 AX 2 AX1 At0l and a constant angular velocity of w 01 The velocity components are u1 wx2 x02 and u2 wxl x01 where x01 x02 50 AX 50 AX One full rotation around the point 39xo x02 was equivalent to 628 time steps In this circumstance the constant in the stability criterion 4 was 099 A cone was used as the initial condition with a base radius of 340 PIOTR K SMOLARKIEWICZ 15 AX and a maximum value of 40 at the point xm1xm2 75 AX 50 AX Fig 7 For all tests performed the same boundary conditions were used The rst partial spatial derivative in the normal direction was assumed to vanish at the out ow boundary vanishing of the second derivative does not ensure the positive de niteness of the scheme The undisturbed value of the eld was assumed at the in ow boundary The solution after six full rotations 3768 time steps is presented in Figs 8 9 and 10 for the scheme in Eq 15 with IORD 2 3 and 4 respectively The maximal values of the presented solutions are 216 317 and 325 respectively and the energy error ER2 is 052 020 and 014 where BRZ is defined as ER2 E 1 w2x x2 I alx1 dx2 1 outflowwz dtl 0 ltgt 71 x w2x x20dx 1ng 34 FIG 7 Initial condition for the twoedimensional tests The scale values in left39front and right back corners are 2 and 4 respectively The scale values are the same in all gures POSITIVE DEFINITE ADVECTION SCHEME 341 FIG 8 Solution for an IORD 2 scheme after six full rotations 3768 time steps where a signi es the whole domain of x x2 All versions of the scheme that were considered are in a conservative form and the error of the conservation of u was to a roundoff error level of 10 Furthermore the minimum value of the solution obtained was 0 It can be seen by comparing Figs 8 9 and 10 and taking into consideration that when a simple upstream scheme IORD 1 is used the solution almost vanishes 13 Fig 2 that signi cant improvements of the algorithm s accuracy are obtained when IORD changes between 1 and 4 For IORD gt 4 the differences between two successive IORD schemes are unnoticeable on the gures eg for IORD 6 the maximum value of the tested solution is 327 and ER2 is 012 When the results presented in Figs 8 9 and 10 are compared with the equivalent results obtained from the timesplitting form of the onedimensional version of the scheme 13 Figs 13 14 and 15 it can be concluded that the combined form of time differencing gives slightly better results than timesplitting The results shown in Figs 8 9 and 10 have been compared with the solutions obtained from the second 342 PIOTR K SMOLARKIEWICZ FIG 9 As in Fig 8 but for an IORD 3 scheme orderaccurate explicit version of the TVD scheme 8 following Theorem 43 and Eqs 35 with 170 310 316b with 620 41 45 48 with g 0 Figure 11 shows the solution after six rotations for the timesplitting form of the chosen version of the TVD scheme Solution maximum and ER2 values are 163 and 052 respectively In Fig 12a the same solution is presented but for the scheme applied in a combined timedifferencing form Solution maximum and ER2 are 152 and 060 Comparing Fig 12a with Fig 11 it can be seen that the combined form of the scheme results in a deformation of the solution elongation normal to the direction of motion this is especially clear in Fig 12b where the solution is shown a quarter of a rotation earlier than in Fig 12a This deformation is typical for all multidimensional higherorder schemes that do not take into account the cross spatial partial derivatives 12 13 It can also be concluded that the chosen version of the TVD scheme is even more diffusive than the IORD 2 version of the proposed algorithm Computer time consumption for both schemes is practically the POSITIVE DEFINITE ADVECTION SCHEME 343 FIG 10 As in Fig 8 but for an IORD 4 scheme same ie about three times that required for the simple upstream scheme For versions of the scheme with IORD 3 and 4 computer time consumption is about ve and seven times that required for the upstream scheme Compare this with the FCT scheme 13 which requires about eight times the upstream amount In the threedimensional case the number of grid points chosen was 41 in each direction with AX AX 2 AX 3 AX 25 At 02 and a constant angular velocity 9 02 w2 w2 12 where w 01 The velocity components are u1 Q3x2 x02 22x3 xo3 u2 23x1 x01 Qlx3 x03 and u3 Qzxl x01 Q x2 x02 where x039 x02 x0 20 AX 20 AX 20 AX The initial condition was a sphere with the radius 7 AX and linearly variable density from 0 on the edge to a maximum value of 4 in the center 20 AX 7 AX6 2 20 AX 7AX6 2 20 AX l4 AX63912 In Fig 13 the values greater than or equal to 05 are plotted The sphere is rotating around the diagonal axis of the grid space that passes through the corner shown in Fig 13 One full revolution around the 344 PIOTR K SMOLARKIEWICZ FIG 11 As in Fig 8 but for the time splitting version of the TVD scheme point xo39x02x03 is equivalent to 314 time steps The maximum value of the Courant number in the stability condition 4 is 095 In Figs 14 and 15 the values greater than or equal to 05 of the solutions after ve revolutions are shown for the four iteration scheme 15 and timesplitting TVD scheme respectively The maximum solution values are 167 and 055 and the ERZ s are 063 and 088 respec tively for both schemes The results for the combined form of the TVD scheme are not presented because the maximum solution value is 050 ER2 091 From Fig 14 it can be seen that the solution is deformed ie elongated toward the center of the rotation Because similar dif culties in a twodimensional advective transport problem have been solved by introducing the secondorder cross terms explicitly to the scheme it was suspected that taking into account the thirdorder cross terms would solve deformation trouble in a three dimensional case Repeating the whole procedure discussed in Section 2 but applying the third order Taylor series expansion POSITIVE DEFINITE ADVECTION SCHEME FIG 12a As in Fig 8 but for the combined version of the TVD scheme instead of the secondorder one leads to a corrected form of the antidi usive velocity 1 l I Ha2 1 ilt12e1qu3 Hume 35 where 1 uilHlZk F nan2n qu uIIZe3 432 ll5 w with 1 X i 2 sum of the above terms EMA1 2 M 2 1 27 2 At Ha2m quotHa2 J31 Jil 346 PIOTR K SMOLARKIEWICZ wite e Wite17 Vite Vi cl sum of the above terms 6 AX AX M if M I l J L 3 1 ui12euil2c quotHume Jl L l thl L J l gtIlt gt1 lVie1eL Wi eeL WieJ eL Wi eJ cL Wi e eLe X WlteleL eV eeJ eu Wig el eL sum of the above terms a AX AXL 39 Note that the structure of the scheme remains the same as in Eq 15 and only new terms are included in the antidiffusive velocities Extending the results of Section 3 it can be concluded that this newly obtained version of the algorithm is thirdorder FIG 12b As in Fig 12a but shown a quarter of a rotation earlier POSITIVE DEFINITE ADVECTION SCHEME 347 FIG 13 Initial condition for the threedimensional tests All points in which the value of a function is greater than or equal to 05 are plotted FIG 14 Solution for the IORD 4 scheme after ve rotations The values greater than or equal to 05 are plotted 58154210 348 PIOTR K SMOLARKIEWICZ FIG 15 As in Fig 14 but the time splitting version of the TVD scheme 5439 FIG 16 As in Fig 14 but for a modi ed thirdorder accurate in time version of the IORD4 scheme POSITIVE DEFINITE ADVECTION SCHEME 349 accurate in time when the original velocity eld is time independent and third or higherorderaccurate in space depending upon the number of iterations The application of four iterations to the new version of the algorithm produces the results shown in Fig 16 after ve rotations The maximum eld value and ER2 are 269 and 054 respectively It can be seen that the deformation that occurred with the basic version of scheme has almost totally vanished Also the implicit diffusion of the scheme is signi cantly decreased On the other hand the timesplitting form of the original scheme gives results similar to those shown in Fig 14 maximum value 173 ERZ 062 Finally it can be concluded that including all thirdorder terms explicitly in the scheme is essential to eliminate the occurring deformation and improve the solution s accuracy The results obtained suggest there is not enough evidence to form a stronger conclusion that isotropic distribution of the truncation error can be reached only for those advection schemes that are at least the same order of accuracy in both time and space as the dimensionality of the problem The results obtained in 12 the cross terms have been applied explicitly to different types of schemes suggest that this formulated hypothesis is generally applicable 5 SOME POSSIBLE OPTIONS OF THE SCHEME In the previous sections the basic form of the algorithm has been presented and discussed in detail The different assumptions on the truncation of the Taylor series expansion in the rst step of the scheme development were shown to lead to the new versions of the algorithm Using the same logical procedures as in the previous sections one can construct many new schemes Some of the existing options will be described below 51 Divergent Flow Field To obtain Eq 7 a uniform velocity eld had been assumed When u E u x Eq 7 takes the form av 71 M al I M aul J 7 quot2157 054 297 W i Jl n 37 where the terms that are not written explicitly are the same as in 7 Note that 37 proves that 7 and 13 are valid not only for a uniform velocity eld but for any velocity eld that is time independent and nondivergent Using 37 results in a new corrected form of 13 1 1 1 1 1 Ha2 quot025 At ull2eul32e ui 12eAX M I J J 025 At Ha2 Z HO2 uii2e e 11 J I ii inn uil l2e eAXJ 38 350 PIOTR K SMOLARKIEWICZ where the terms omitted are the same as in 13 When the original velocity eld is nondivergent the new terms in 38 are insigni cant They vanish in the second iteration of the scheme while in higher iterations they are always an order smaller in AX than the original terms It has been found that application of Eq 38 to the tests presented in Section 4 practically does not affect the previously obtained results On the other hand when the velocity is strongly convergent application of Eq 38 to the problem of the evolution of the droplet size distribution due to the evaporation con densation process improves the results William Hall personal communication 52 Time Dependent Velocity Field When the time dependence of the velocity eld is taken into account Eq 7 should be written in the form 39 FIG 17 Initial condition for the shock type initial condition problem tests Scale values in corners are the same as for Fig 7 POSITIVE DEFINITE ADVECTION SCHEME 351 where the terms not written explicitly are the same as in 7 The antidiffusive velocity 13 is written in the form aul n ugwm 05 At 40 at il2e where the omitted terms are the same as in 13 It is still unclear how to properly approximate the velocity time derivative in 40 in the case of a nonlinear equation such as the momentum equation In applications in which time dependence of the velocity eld is given a priori Eq 40 may be used in the second iteration of the scheme Tests performed where the angular velocity from the previous section was allowed to be time dependent wcosa I suggest that the effect of Eq 40 is negligible This problem will be investigated more systematically in the future FIG 18 Solution for the IORD 4 scheme after one revolution 352 PIOTR K SMOLARKIEWICZ 53 Schemes with a Higher Order of Accuracy in Time In Section 4 the development of the scheme that is third orderaccurate in time was presented One can similarly construct schemes with an optional order of accuracy in time and space Increasing the order of accuracy of the scheme results in additional terms in Eq 7 and dramatically complicates the form of the antidif fusive velocity 13 It has been found that the application of the versions of the scheme that are third and fourthorder accurate in time in both a timesplitting and a combined form has little effect on the results of the two dimensional tests from the previous section It has been found that the third or fourthorder accurate in time version of the scheme may have an important application to the moving shock problem 54 The Shock Type Initial Condition Problem To illustrate the behavior of the scheme for a shocktype initial condition it is convenient to repeat the test proposed by Zalesak 14 In the case of the two FIG 19 As in Fig 18 but for the IORD 2 scheme posmvra DEFINITE ADVECTION SCHEME 353 IA n l I V 0 1 quot FIG 20 As in Fig 18 but for the TVD scheme dimensional test from the previous section the cone is replaced by a cylinder with the same base radius amplitude and position as the cone with a deep groove as shown in Fig 17 In Figs 18 19 and 20 the solutions after one revolution as in 14 are shown for IORD4 IORD2 and a combined TVD scheme The maximum values and ER2 are respectively 476 and 028 382 and 046 and 334 and 050 It can be seen that while the TVD scheme nearly maintains the at shape of the cylinder s top the proposed algorithm especially for IORD 4 creates an artificial maxima This effect is an obvious consequence of the scheme structure and will always occur in the closest neighbourhood of the shock or contact discon tinuity In the category of problems in which this local ampli cation is unacceptable it is possible to eliminate this effect for example by turning off the corrective procedure in the closest neighbourhood of the shock As an example the semi empirical switch is presented 354 prom K SMOLARKIEWICZ Swil mm min 32e5 3 12 413 lV Vigtel Wig l yl e 1 0531 sin ESLZ 41b I1Zt 8 Wfwitelg where the constant E w w e 8 l S was found experimen tally to be 10 22E8203 10 The sw1tch lS constructed only after the rst iteration of the scheme and it remains constant in all consequent iterations Note that SW is u turns to zero when sually equal to unity and auxquot ax il2c w liq 3 umquot N I n FIG 21 As in Fig 18 but for the time splitting version of the IORD 4 scheme with the switch 41 applied POSITIVE DEFINITE ADVECTION SCHEME 355 and 17 H 2e39 0 Modi ed antidiffusive velocity takes the form lquot k ui12 uii2eEq13 X SWi12e k 2 1am IORD 1 42 The Optimal value of the constant ES in 41b was found to be 03 x 10 2 for all cases tested Larger values of the ES produce more diffusive solutions while smaller ones do not eliminate the ampli cation effect totally The solutions for the IORD 4 scheme applied in a timesplitting and a combined form are shown in Figs 21 and 22 respectively The slight difference between these two solutions is caused by the switch rather than the scheme itself The switch should not be universally applied because it affects the solutions for nonshocktype initial conditions In Fig 23 the solution after six rotations of the cone is shown for the IORD 4 scheme with the switch 41a b applied The solution obtained is less diffusive than for the TVD scheme FIG 22 As in Fig 21 but for the combined scheme 53154241 356 PIOTR K SMOLARKIEWICZ Figs 123 b but signi cantly more diffusive than for the original IORD 4 scheme Fig 10 The universal character of the constant ES is still unknown and in this point the switch presented needs to be improved In its current form it should be interpreted only as an example indicating one of the possibilities of the shock problem solution In the special case of a shock occurring on the lee side of the initial condition a monotonicity preserving solution is provided by the thirdorderaccurate in time version of the algorithm 36 In Figs 24a and 24b the results of the one dimensional tests u AX 1 At 05 are shown respectively for the basic l3 and modi ed 36 versions of the scheme For comparison the solutions for the TVD scheme are also plotted dashed lines It can be seen that application of the antidiffusive velocities given by Eq 36 totally eliminates the ampli cation effect produced by the basic version of the algorithm The solution obtained is less diffusive than the one FIG 23 Solution after six rotations for the scheme shown in Fig 22 but applied to the initial condition in Fig 7 POSITIVE DEFINITE ADVECTION SCHEME 357 produced by the TVD scheme On the other hand when the shock occurs on the wind side of the initial condition the ampli cation effect produced by the basic version of the scheme Fig 25a transforms to oscillations for the modi ed version of the algorithm Fig 25b This result suggests a different construction of the switch than previously discussed For example if in the rst term of the onedimensional version of Eq 36 uAX2 is replaced by the sign xi3quot wi l ul AX2 the result will be a scheme with properties similar to those of the tested version of the TVD scheme 55 Practical Simpli cation of the Scheme In cases when computational ef ciency is essential for a problem one may nd the IORD 2 scheme is the only practical one to apply Although this simplest version of the algorithm produces results of the quality comparable with that of the TVD scheme it is significantly less accurate than the versions with IORD gt 2 In 13 a simple compromise between accuracy and computational ef ciency was presented FUNC oo l 00 200 400 600 800 1000 FIG 24a Different stages of the shock heavy solid line propagation 40 80 120 time steps Solid lines IORD 4 scheme dashed lines TVD scheme 358 PIOTR K SMOLARKIEWICZ FUNC l 00 200 400 600 800 1000 FIG 24b As in Fig 24a but for IORD 4 third order accurate in time scheme Because usually the antidiffusive velocities are smaller than the ones allowed by the stability criteria it is possible to increase them by multiplying by some factor Sc 17i 12 il2elEq 13 X 50 43 It was found experimentally that even a small increase of Se over unity eg Sc 106 signi cantly improves the solution for the IORD 2 scheme In Fig 26 the solution of the two dimensional test from Section 4 is shown after six rotations for the IORD 2 scheme with Eq 43 applied and Sc 106 cf Fig 8 The maximum solution value and ER2 are respectively 317 and 031 It was found that the optimal value of the coefficient Sc depends upon the CFL number and dimensionality of the problem The character of this dependence is still unknown In some problems it may be worthwhile to tune Sc using some simple tests The scheme can also be simpli ed in a different way In the presented basic version of the algorithm the general form of the antidiffusive velocity is 22 17w 14 It was found that the form 171701 u can be used as well This new form of the 20 r I r r I h Ia Ifr I I I I I I 4 s I 0 i u L 1 Z 3 LI I o0 r o I 4 I gI I I I I II I I I40 I I I L 00 200 400 000 800 1000 X Flo 25a As in Fig 24a but for a shock occurring on the wind side of the initial condition 2I0 quot I v I n I i I 7 i I I v I I I 1 0 LJ 39 0 2 D LL 0o I 10 I 4 I g I I I I I II l I I g I I I I 00 200 400 000 300 1000 X FIG 25b As in Fig 24b but for an initial condition as in Fig 25a 360 PIOTR K SMOLARKIEWICZ antidiffusive velocity results in slightly larger truncation errors but gives nal results close to those produced by the original algorithm Using this information the antidiffusive velocity for the IORD 2 scheme may be modi ed that is modi ed 17 quot17 W 443 or modi ed I W 44b In Fig 27 the solution of the cone rotation test after six revolutions is shown for the IORD2 scheme with 44a applied The maximum solution value and ER2 are 302 and 023 respectively cf Figs 8 and 26 The disadvantage of this method is a more restrictive stability criterion of g g 05 in 4 FIG 26 As in Fig 8 but with Sc 106 in Eq 43 applied POSITIVE DEFINITE ADVECTION SCHEME 361 FIG 27 As in Fig 8 but with Eq 44a applied 6 CONCLUSIONS 1 Using an iterative approach based upon the upstream scheme a class of fully multidimensional nonlinear computationally ef cient positive de nite advective transport algorithms has been constructed The simplest version of the developed schemes is secondorderaccurate in both time and space The most accurate version tested in this paper is thirdorder accurate in time and fth orderaccurate in space In principle it is possible to construct an algorithm with an optional order of accuracy in time and space 2 Results presented in this paper suggest that there exists a general connection between the distribution of the truncation error of any advection scheme and the dimensionality of the problem It is suggested that to ensure that the numerical 362 PIOTR K SMOLARKIEWICZ solution of the advection equation be free of strong arti cial deformations it is necessary to use schemes of at least the same order of accuracy in both time and space as the dimensionality of the problem 3 The procedures were discussed whereby the basic form of the algorithm can be modi ed depending on the user s requirements Among others such options as a generalization of the scheme for a divergent flow eld case or the construction of a monotonicitypreserving scheme have been mentioned ACKNOWLEDGMENTS The author is grateful to Dr Terry L Clark of NCAR for his suport in this project as well as for his helpful discussion and criticism of the manuscript The author also wishes to thank Dr William D Hall for his criticism of the manuscript and Mary Davis for typing the manuscript I wish to acknowledge the support received from the Advanced Study Program of NCAR while doing this research REFERENCES L BOOK J P BORIS AND K HAIN J Comput Phys 18 1975 248 BORIS AND D L BOOK J Comput Phys 11 1973 38 BORIS AND D L BOOK J Comput Phys 20 1976 397 CLARK J Atmospheric Sci 36 1979 2191 CLARK AND W D HALL J Atmospheric Sci 36 1979 470 A ARTEN Math Camp 32 1978 363 A HARTEN High Resolution Scheme for Hyperbolic Conservation Laws New York Univ Report DOEERO3077 167 March 1982 A HARTEN On a Class of High Resolution Total VariationStable FiniteDifference Schemes New York Univ Report DOEERO3077 l76 October 1982 9 A HARTEN AND G ZWAs J Comput Phys 9 1972 568 10 R D RICHTMYER AND K W MORTON Difference Methods for InitialValue Problems p 405 Interscience Wiley New York 1967 11 P J ROACHE Computational Fluid Dynamics p 446 Hermosa Publ Albuquerque 1976 12 P K SMOLARKIEWICZ Monthly Weather Rev 110 1982 1968 13 P K SMOLARKIEWICZ Monthly Weather Rev 111 1983 479 14 S T ZALESAK J Comput Phys 31 1979 335 l 2 P 3 P 4 L 5 L 6 H 7 9
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'