（自然地理学专业论文）质点跟踪法及其应用.pdf

摘要 用数值模式对河口海岸地区的物质输运进行计算时，平流项的数值格式必须 要能正确处理物质浓度锋面，以避免产生过多的数值耗散或频散。最基本的中央 差格式、迎风格式和l a x 一7 9 e n d r o f f 格式不是耗散过大就是频散太厉害。欧拉一 拉格朗日方法虽然能消除频散，且耗散比较小，但仍然存在一定的耗散，所以有 必要进一步改进。此外，欧拉拉格朗日方法的水平插值过程如果在。面内进行， 那么。面的弯曲还会使得插值依据不合理。所以本文首先对欧拉一拉格朗日方法 作了改进，把其中的插值过程改在绝对坐标下进行，使得水平插值过程在同一高 度的水平面内进行，从而克服了插值不合理的问题。然后，本文又设计了质点跟 踪法，该方法通过在网格内设置一些质点并对质点进行跟踪来计算平流项。因为 没有像欧拉一拉格朗目方法那样反复插值，所以质点跟踪法极大地降低了数值耗 散。质点跟踪法的低耗散性在本文的一维试验和二维试验中得到了体现，而改进 后的欧拉一拉格朗日方法也在二维试验中优于改进前的欧拉一拉格朗日方法。最 后，本文在三维e c o m 模型下把质点跟踪法应用到了对长江口污染物扩散的数值 模拟中。在该模拟中，本文通过对c o d 浓度的计算，又把质点跟踪法和改进后 的欧拉一拉格朗日方法作了比较。结果显示质点跟踪法不仅耗散更低，而且守恒 性也更好。同时，通过对污染物扩散的数值模拟，本文发现长江口污染物的扩散 范围能够到达祟明东滩、深水航道和南汇边滩，而没冒沙附近区域的c o d 浓度 更是尤其地高，所以污染物扩散对没冒沙地区的水质也形成了一定的威胁。 关键词：平流项：数值格式；质点跟踪法；污染物扩散；c o d a b s t r a c t t h en u m e r i c a ls c h e m eo ft h ea d v e c t i o nt e m ls h o u l db ed e s i g n e dc o r r e c t l yf o r t h e f r o n t so fs h a r pc o n c e n t r a t i o ng r a d i e n t st op r e v e n te x c e s s i v en u m e r i c a ld i f f u s i o no r d i s p e r s i o nw h e nc a l c u l a t i n gi n a s st r a n s p o r ti nt h ea r e a so fe s t u a r ya n dc o a s t s o m e b a s i cs c h e m e s ，s u c ha sc e n t r a l d i f f e r e n c i n g ，u p w i n da n dl a x w e n & o f f , a r ee i t h e r h i g h l y d i f f u s i v eo re x t r e m e l yd i s p e r s i v e e u l e r i a n - l a g a l l g i a nm e t h o d ( e l m ) c a n p r e v e n t t h ed i s p e r s i o nw h i l el o w e rt h ed i f f u s i o n ，b u ti t sd i f f u s i o ni ss t i l lc o n s i d e r a b l e ， s oi ti sn e c e s s a r yt ob ef u r t h e ri m p r o v e d i na d d i t i o n ，i ft h eh o r i z o n t a li n t e r p o l a t i o n p r o c e s sw a st a k e ni na6p l a i n ，t h e nt h eb e n d e d 巧p l a i nw o u l d m a k et h ei n t e r p o l a t i o n s o u r c e su n r e a s o n a b l e s oi nt h i s p a p e r ，f t r s t ，t h ee l mi si m p r o v e db yt a k i n gt h e i n t e r p o l a t i o np r o c e s su n d e ra b s o l u t ec o o r d i n a t e s ，s ot h a tt h eh o r i z o n t a li n t e r p o l a t i o n c o u l db et a k e ni nah o r i z o n t a l p l a i no ft h e s a m ea l t i t u d e ，t h u st h e p r o b l e mo f u v x e a s o n a b l ei n t e r p o l a t i o ns o u r c e si sp r e v e n t e d t h e n ，t h ep a r t i c l e t r a c i n gm e t h o di s d e s i g n e d ，w h i c hc a l c u l a t e st h ea d v e c t i o nt e r mb ys e t t i n gs o m ep a r t i c l e si nt h e 鲥d c e l l sa n d t r a c i n g t h e m w i t h o u t r e p e a t e di n t e r p o l a t i o n s s u c ha s i n e l m ，t h e p a r t i c l e t r a c i n gm e t h o d r e d u c e st h en u m e r i c a ld i f f u s i o nq u i t eal o t t h ei o wd i f f u s i o n o ft h e p a r t i c l e - - t r a c i n gm e t h o d i ss h o w ni nao n e - - d i m e n s i o n a le x p e r i m e n ta n da t w o - d i m e n s i o n a lo n e ，w h i l et h ei m p r o v e de l m p e r f o r m sb e t t e rt h a nt h eo r i g i n a le l m i nt h et w o d i m e n s i o n a l e x p e r i m e n t f i n a l l y , u n d e rt h e t h r e e d i m e n s i o n a le c o m m o d e l ，t h ep a r t i c l e t r a c i n gm e t h o di sa p p l i e d t ot h en u m e r i c a lm o d e l i n go ft h e c h a n g i i a n ge s t u a r ys e w a g ed i f f u s i o n i nt h i sm o d e l i n g ，t h ep a r t i c l e 一订a c m gm e t h o d i s a g a i nc o m p a r e dw i t ht h ei m p r o v e de l mb yc a l c u l a f t n gt h ec o n c e n t r a t i o no fc o d t h er e s u l ts h o w st h a tt h ep a r t i c l e - t r a c i n gm e t h o di sn o to n l yl e s sd i f f u s i v e ，b u ta l s o m o r ec o n s e r v a t i v e m e a n w h i l e ，b ym o d e l i n gt h es e w a g ed i f f u s i o n ，i ti sf o u n dt h a tt h e r a n g eo f t h ec h a n n i a n ge s t u a r ys e w a g ed i f f u s i o na r e ac a nr e a c ht h ee a s tc h o n g m i n g b a n k ，t h ed e e p w a t e rc h a n n e l ，a n dt h en a n h u ib a n k ，w h i l et h ec o dc o n c e n t r a t i o n n e a rm e i m a o s h ai se s p e c i a l l yh i 曲s ot h a tt h es e w a g ed i f f u s i o nc o n s t i t u t e sat h r e a tt o t h ew a t e r q u a l i t yo f m e i m a o s h a a r e a k e y w o r d s ：a d v e c t i o nt e r m ；n u m e r i c a l s c h e m e ；p a r t i c l e - t r a c i n gm e t h o d ；s e w a g e d i f f u s i o n ；c o d 陈晒睿硕士学位论文答辩委员会成员名单 姓名职称单位备注 沈焕庭教授华东师范大学主席 丁平兴教授华东师范大学 端义宏研究员上海市气象局 上海市环境科 林卫青研究员 学研究院 交通部科学研 戚定满研究员究院河口海岸 研究中心 学位论文独创性声明 本人所呈交的学位论文是我在导师的指导下进行的研究工作及 取得的研究成果。据我所知，除文中已经注明引用的内容外，本论文 不包含其他个人已经发表或撰写过的研究成果。对本文的研究做出重 要贡献的个人和集体，均已在文中作了明确说明并表示谢意。 作者签名：堡堕垒日期：里一! ! 茎曼璺印 学位论文授权使用声明 本人完全了解华东师范大学有关保留、使用学位论文的规定，学 校有权保留学位论文并向国家主管部门或其指定机构送交论文的电 子版和纸质版。有权将学位论文用于非赢利目的的少量复制并允许论 文进入学校图书馆被查阅。有权将学位论文的内容编入有关数据库进 行检索。有权将学位论文的标题和摘要汇编出版。保密的学位论文在 解密后适用本规定， 学位论文作者签名：豫蛹苷 日期：2 口。壁6 甩1 2 日 导师签名： 魄槲跳6 , j 文 第一章关于物质输运方程中平流项的数值格式 1 i 研究意义 本文所要介绍的质点跟踪法是一种在数值模式中求解物质输运方程中平流 项的方法。河口海岸地区水体中的物质输运问题当前正日益受到人们的关注。与 河口地区人民群众的生存环境密切相关的河口盐水入侵、泥沙淤积和冲郦、污染 物扩散等问题都和物质输运有关。而各项重大工程又将会给河口地区物质输运带 来什么样的影响，也是人们所迫切希望了解的。作为一种有效的手段，用数值模 式对河口地区物质输运问题进行模拟计算正变得越来越重要。 在数值模式中，处理物质输运的过程就是数值求解物质输运方程。求解物质 输运方程时，先把物质输运方程分解为平流项和扩散项等，再先后对不同的项分 别求解。在对平流项的求解过程中，根据数值格式的不同，一般总是或多或少会 存在一些计算耗散或频散i ”，给模拟结果带来误差。比如时间前差，空间中央差 的欧拉差分格式会导致耗散和频散；迎风差分格式会产生较大的耗散： l a x w e n d r o f f 格式虽然减少了耗散，但是其频散却相当严重。正因为这样，人 们又不断改进发展了许多格式，如t v d 格式、p p m 格式、欧拉一拉格朗日方法等， 它们用不同的方法减小了数值误差。比如t v d 格式和p p m 格式都是数值耗散较小 的格式，而欧拉一拉格朗目方法能在完全消除数值频散的同时降低数值耗散。但 是这些格式仍然存在一定的频散或者耗散，在某些情况下，比如在锋面处，这些 格式仍然显得不够精确。此外，很多平流项数值格式不仅会在水平方向产生一定 的数值耗散，减弱锋面的强度，而且在。坐标下更有可能在水深变化剧烈的地方 造成强烈的垂向数值耗散，从而使物质浓度在垂向不合理地变得均匀。所以，数 值耗散和数值频散是求解平流项时的一个令人头疼的问题，其所带来的误差有时 甚至会使人改变对模拟结果的判断。为了克服数值耗散或数值频散，防止数值耗 散或数值频散给模拟带来误差，人们仍然在不断探索更好的平流项数值格式。本 文所设计的质点跟踪法，正是一种为了减少数值耗散，提高平流项的计算精度而 作的平流项数值格式。 1 2 研究现状 在各种处理物质输运方程中平流项的数值格式中，有很多格式都是通过考虑 网格之间的物质变量通量来计算各网格的物质变量值的，所以可以把它们统称为 通量格式。在通量格式中，由于网格的物质变量值的变化是用通量来计算的，所 以物质守恒性好。也有很多格式不是通过通量来计算的，比如欧拉一拉格朗日方 法，它是采用逆向跟踪水质点并插值来计算网格的物质变量值的。虽然非通量格 式不一定能保证物质守恒，但是它们也有别的优势，比如能更方便地消除频散。 本文所设计的质点跟踪法就不是通量格式。 一般人们往往从最简单的维情形开始研究数值格式。在一维情形下，关于 物质变量a ( x ，f ) 的平流方程在速度u 为常数时为： o t ：一u o a 出 ( 1 1 ) 对于通量格式，把( i 1 ) 式写为通量形式看起来更方便： 如，o t = o ( u 趴，6 ) ： ( 1 2 ) 欧拉格式、迎风格式和l a x w e n d r o f f 格式是通量格式中最为基本的3 种。 欧拉格式是时间前差、空间中央差的有限差分格式。欧拉格式里的通量f ， 件j 为： f 。= “( 口? + a 三 ) 2 ( 1 3 ) j + = 欧拉格式对方程( 1 2 ) 的解为： 口? “= d ? 一u a t ( a i + j 一口三1 ) 2 & c ( 1 4 ) 该解是无条件不稳定的。 l a x w e n d r o f f 格式的通量，为： 件i fl = u ( a ? + 口：1 ) 2 + i t 2 a t ( a ? 一口三1 ) 2 a x ( 1 - 5 ) ”j l a x - w e n d r o f f 格式对方程( 2 ) 的解为： n ? “= a ? 一u a t ( a ? + 1 一口三1 ) 2 a x + 2 f2 ( 日三l + 口三1 2 a ；) 2 a x 2 ( 1 6 ) 当训址缸时( 这是c o u r a n t - f r i e d r i c h s l e w y 即c f l 条件) ，l a x w e n d r o f f 格 式是稳定的。 迎风格式的通量f ，为： 什j f 】= “( d ? + 口：1 ) 2 + i “i ( 口? 一口三1 ) 2 ( 1 7 ) “i 迎风格式对方程( 2 ) 的解为： 臼? “= 口? + “r ( d 三l a ；) a x ，当“0 时； ( 1 8 ) d ? “= 口? + u a t ( a , 一口：1 ) a x ，当“ 0 时； ( 1 9 ) 在满足c f l 条件的情况下，迎风格式是稳定的。 以上3 种格式一般认为是最基本的三种通量格式。后来通量格式有了很大的 发展，演变出很多种方法，而一般可以像y a n ga n dp r z e k w a s 【2 那样把这些方法 中的绝大多数分为代数格式和几何格式两大类。代数格式基本上都是混合格式， 其通量由一个低阶耗散基本格式( 比如一阶迎风格式) 的通量和一个高阶频散基 本格式( 比如l a x w e n d r o f f 格式) 的通量组合而成。各种代数格式之间的不同在 于它们计算两个基本格式的通量所占权重的方法不同。几何格式则是在每个网格 单元内构造一个关于平流物质变量的变化函数，通过保持其单调性来保证不产生 新的极值，交界面处的通量可以用某种上游方法计算。最低阶的几何格式为每个 网格单元内没有变化，此肘相当于一阶迎风格式。可以在每个网格单元内构造线 性变化既而种名口qp p m 的几何格式则在网格单元内构造了抛物线变化i d j a m e s 5 1 曾测试了一种叫做t v d 的代数格式和几何格式p p m 在陆架浅海的表现， 其结果显示这两种格式在该试验中都无频散、低耗散，而p p m 的数值耗散更小。 本文由于需要把t v d 格式和其它格式作比较，所以在此介绍一下本文所使用的 t v d 格式： t v d 格式是一种代数格式，即它的通量由个低阶格式的通量和一个高阶格 式的通量组合而成。这里用的t v d 格式是由一阶迎风格式和l a x 一- w e n d r o f f 格式 组合而成( 5 1 ，其总通量为： f ，= f 吁+ ( r ) ( f 。，一f 吁) ( 1 1 0 ) + j+ i “i1 + j 其中f ”和f “分别是由迎风格式和l a x w e n d r o f f 格式所提供的通量。函数 痧( ，) 是根据物质浓度值口的变化确定的限制函数闭： r = f 掣f 甲 ( 1 1 1 ) h j 哪h j 其中： f 彳“= f - f 咋 ( i 1 2 ) 什j件i什i fj = ( h1 1 ) 一a t u 2 l x ) ( d 一口，) ，2 ( 1 1 3 ) lh ij什i 以及： f 2 s f g n ( u ，。! ) ( 1 1 4 ) 2 如果，掣= 0 ，则取，= 0 ，此时等式( 1 1 0 ) 对任何有限的( r ) 都有： 什j f 。= ，吁= f 。 什ih j件j 有很多限制函数可供t v d 格式调用1 x , z a 】。其中包括： m i n i m o d e ：( ，) = m a x ( o ，m i n ( 1 ，) )( 1 1 5 ) m o n o t o n i c ：( ，) = ( ，+ ，j ) “1 + ，)( 1 ，】6 ) m u s c l ：声( ，) = m a x ( 0 , m i n ( 2 ，2 r ，( 1 + r ) 1 2 ) i )( 1 1 7 ) s u p e r b e e ：( ，) = m a x ( o ：。m i n ( 2 r ，1 ) ，r a i n ( r , 2 ) )( 1 1 8 ) 这些限制函数具有不同的保持锋面的能力。本文的例子里采用s u p e r b e e 限制函 数，它被认为是耗散最小的限制函数之一f 2 ，- 甜。 欧拉一拉格朗目方法则是与通量格式不同的另一类平流项数值格式弘。1 ”。这 类方法在计算物质输运方程中平流项时，首先采用逆向跟踪法确定水质点的历史 位置，再用拉格朗目插值的方法确定该点的平流物质变量的值。以一维平流方程 ( 1 1 ) 为例，该式可以写为拉格朗日形式： d a d t ：0 ( 1 ，1 9 ) 即水质点的物质浓度值在只考虑平流的情况下不随时间而改变。所以如果n + l 时刻位于网格点q 的水质点在n 时刻位于尸点，那么只考虑平流有： 口“( q ) = a ”( 尸) ( 1 2 0 ) 所以只要确定尸点在1 9 时刻的坐标位置，就可以用拉格朗目插值法求出a “( p ) 从 而得到a n * l ( 1 2 ) 的值。在二维或三维情形中的处理方法也同样如此。确定p 点在n 时刻的坐标位置，可以用多步逆向跟踪法或者预估修正法等方法确定；插值 ( p ) 则根据插值方法的不同有好几种方法。一般来说，一阶拉格朗日插值耗散 太大，二阶插值往往比较适用。z h u 8 】等曾介绍了种用一阶与二阶相结合的插 值方法插值日“( p ) 的欧拉一拉格朗日方法。欧拉一拉格朗日方法的一大优点是安全， 它一方面较少受到c f l 条件的限制，另一方面通过对插值结果的控制可以完全 避免出现新的极值。 1 3 本文的工作 本文的工作主要是对物质输运方程中平流项的数值格式作改进。人们之所以 仍然在大量地研究平流项数值格式，一个很重要的原因就是这些格式总还有不足 之处。比如像t v d 格式或者p p m 格式那样的通量格式，它们可能在一般情况下 能够做到无频散、低耗散，但是在c f l 条件不满足的情况下就难免要出现频散， 这使得网格局部加密或者时间步长的加大受到了更多的限制。又比如欧拉一拉格 朗日方法，它虽然能够有效防止频敞的发生，同时耗散也比较小，但是在物质守 恒性方面却不如通量格式。欧拉拉格朗目方法的水平插值过程如果在。面内进 行，那么。面的弯曲还会使得插值依据不合理。所以本文首先对欧拉拉格朗目 方法作了改进，把其中的插值过程改在绝对坐标下进行，以克服插值不合理问题。 而无论是t v d 格式还是欧拉一拉格朗日方法，它们都还有一定的耗散，在模拟精 度要求较高的场台仍显不足。所以本文又设计了质点跟踪法，以进一步减小数值 耗散。通过本文所作的几个试验，可以看到质点跟踪法不仅能消除频散，而且耗 散比其它各种格式更是明显减小。最后，本文将质点跟踪法试用于对长江口污染 物扩散的计算中，并通过和欧拉- 拉格朗日方法的比较，进一步说明其优越性。 第二章平流项数值格式的改进 本文对平流项数值格式的改进分为两个部分；第一部分为对欧拉一拉格朗曰 方法的改进，该改进使得欧拉拉格朗日方法的插值过程更加合理；第二部分设 计了质点跟踪法，该方法相对于欧拉拉格朗日方法，又极大地降低了数值耗散。 2 1 对欧拉一拉格朗日方法的改进 用欧拉一拉格朗日方法对平流项进行计算时，其水平插值过程如果是在。面 内进行的，那么在水深变化剧烈的区域可能会导致插值不当。在这些区域，0 一 面相邻网格由于。一面的弯曲，往往在垂向距离悬殊。这时，由于采用了0 一面内 的相邻网格作为插值依据，很有可能会带来远处非必要的信息，从而影响计算的 精度。为此，本文改进了插值方案，使得欧拉一拉格朗日方法的插值能够始终引 用距离相近的点作为插值依据，从而避免了引用垂向距离较远的。一面相邻网格。 在垂向采用。坐标，水平采用曲线坐标的三维e c o m 模式下，本文对欧拉一 拉格朗日方法中的插值方案进行了改进。在e c o m 模型中，每个网格点有3 个坐 标：水平l l 方向坐标i ( 自然数标号) 、水平方向坐标j ( 自然数标号) ，垂 向0 一坐标o ( 海表为0 ，海底为一1 ) 或者k ( 自然数标号) 。同样，在某个网格内的 任何一个点也可以用3 个相距离坐标来表示它在这个网格里的相对位置：水平e 1 方向相坐标d ，水平岛方向相坐标d 2 ，垂向相坐标d 3 。这3 个相坐标的取值都 介于一o 5 和0 5 之间。有了这3 个方向的坐标和相坐标，以及知道水平网格大 小和垂向。层分层比例，就可以通过3 个方向先后的插值来计算任何一个点的 值。改进前的方法就是先在各口一面里分两个方向先后插值，然后在垂向插值得 到最终结果的。在水深变化剧烈的地方，由于。一面是一个曲面，所以。一面里的 插值有时并不准确。故改进后的插值采用绝对坐标：水平t l 方向坐标i 和相坐 标d ，水平匕方向坐标j 和相坐标d 2 ，垂向绝对高度坐标z 。在改进后的新方案 里，先把网格点垂向插值到同一个绝对高度z 的水平面，然后再分两个水平方向 先后插值得到结果，这样就避免了由0 一面弯曲引起的误差。 改进前后所使用的欧拉一拉格朗日方法均采用预估修正法确定本文1 2 节所 述之p 点在n 时刻的位置，而插值则都采用防止越界出现的二阶插值法。 预估修正法的步骤如下： ( 1 ) 根据n + l 时刻位于网格点q 处的水质点速度订，以及单位时间步长出，初步 估算出该水质点在n 时刻的位置p ；( ，+ ：p 一弘r ) ( 2 ) 用插值的方法求出p 在n 时刻的速度，把它和q 在n + 】时刻的速度厅平均， 以估算出水质点从n 时刻到n + 1 时刻的平均速度i ，并以此为依据再次估算出 水质点在n 时刻的位置p ”，作为对p 的修正；( p ”= q i a t ) ( 3 ) 用同样方法，根据p “再次估算出p ”+ ，作为对尸“的修正；p “即作为所求 点p 的结果。 通过测试，可以发现修正2 次以后的位置已经非常稳定，所以没有必要继续 修正。 防止越界出现的二阶插值法是指：根据3 个已知点插值一个未知点( 速度或 物质浓度) ，插值采用二阶拉格朗日插值公式，但插值结果被限制在这3 个已知 点的最大值和最小值之间，即插值结果一旦大于最大值则取最大值，一旦小于最 小值则取最小值。这样就避免了越界的出现，从而消除了频散。 改进前后的插值方法的不同之处在于插值顺序和所选择的坐标依据。 改进前的插值法的具体过程如下： 先确定被插值的点p ( 这里记为p ，插值尸或p ”的速度时方法类似) 被夹 在哪2 个。一面内，把，竖直投影到这2 个。一面上，以便最后用这两个d 一面上的 投影点的值垂向一维插值p 的值；再把各。一面上的投影点沿一个水平网格方向 ( 不妨记为i 方向) 投影到最近的3 条另水平网格方向( 不妨记为 2 方向) 的网 格线上，以便用网格线上的投影点的值插值。一面上的投影点的值；而各网格线 上的投影点的值可以直接由该网格线上最近的3 个网格点的值插值得到。 图2 1 显示了上面描述的插值顺序：先通过每3 个网格点插值出1 个网格线 投影点得到6 个网格线投影点的值，进而通过每3 个网格线投影点插值出1 个。 一面投影点得到p 1 和p 2 的值，最后用p 1 和p 2 一维插值出p 的值。其中插值网 格线投影点和。一面投影点时都是利用防止越界出现的二阶插值法。 闻格 网格 a 面1 面投影点 一面2 图2 1 原来的插值过程 改迸前的方法在实旌过程中有时会遇到p 位于最上层的q - 一面之上或者最底 层的。一面之下的情况，这时，j p 就直接取最上层或最底层的0 一面投影点的值。 改进后的插值法的具体过程如下： 先把水平二维网格竖直投影到与点p 在同一绝对高度的水平面上，并在被投 影的水平面网格里找到与被插值的点p 最近的3 条e 1 方向的水平网格线和3 条 方向的水平网格线，这些网格线相交得到9 个网格点，其中每个网格点的值都可 以用垂向距离该点最近的2 个。一面上的网格点的值一维插值而得：然后，在水 平面网格里，把点p 沿e 。方向投影到最近的3 条 z 方向的网格线上，以便用网 格线上的投影点的值插值户的值；而其中各网格线上的投影点的值可以由该网格 线上最近的3 个网格点的值插值得到。 图2 2 显示了上面描述的插值顺序：先通过每2 个。一面上的网格点一维插 值出1 个水平面网格点得到9 个水平面网格点的值，进而通过每3 个水平面网格 点插值出1 个网格线投影点得到3 个网格线投影点的值，最后通过3 个冈格线投 影点插值出p 的值。其中插值网格线投影点和点尸时都是利用防止越界出现的二 阶插值法。 图2 , 2 新的插值过程 改进后的方法在实施过程中有时会遇到水平面网格点位于最上层的。一面之 上或者最底层的。一面之下的情况，这时，水平面网格点就直接取最上层或最底 层的。一面上的网格点的值。 简而言之，该进前的方法是先在。一面里插值，再在垂向插值；改进后的方 法是先在垂向插值，再在同一绝对高度的水平面里插值。 对比改进前后的两种插值方案：改进前的插值方案可能会由于。一面的弯曲 而引用垂向距离较远的o 面相邻网格点的值，并且在弯曲的线上插值也会造成 一定的不准确：而改进后的插值方案总是引用最近的点，同时所有插值都在直线 上进行，所以提高了准确性。以图2 2 为例，如果点尸在。一面1 和。一面2 之间， 那么在改进后的插值方案里最左边的3 个水平面网格点将引用。一面2 和0 一面3 上的网格点的值，而改进前的插值方案则会在相同水平坐标处引用。一面l 和。一 面2 上的网格点的值，有所不同；而如果点p 在。一面2 和。一面3 之间，那么在 改进后的插值方案里最右边的3 个水平面网格点将引用o 面1 和。一面2 上的网 格点的值，而改进前的插值方案则会在相同水平坐标处引用。一面2 和。一面3 上的网格点的值，其中。一面3 上相应位置的网格点低得在图中看不见，引用这 些点显然欠妥。 2 2 质点跟踪法 质点追踪的方法由来已久，比如早期c h o r i n 【1 2 】就作过这样的工作。而近期 江文胜、孙文心【1 3 , 1 4 1 则改进了悬浮物三维粒子追踪模式( s p m 模式) ，并用它模拟 了黄河排出的细颗粒物在渤海中的运动和大风情况下悬浮物的输运。不过总体来 说，质点追踪法目前在海洋研究中的应用并不是很多。而本文则要用类似的方法 改进对物质输运中平流项的计算。 本文跟据质点追踪的思想，设计了种在水平取曲线坐标、垂向取d 坐标的 三维e c o m 模式下的新的计算物质输运方程平流项的数值格式。由于这种格式采 用了在网格( 由水平网格和垂向。层分割而成的最小单元) 内设置一些质点并对 质点进行跟踪的方法，所以本文将之称为质点跟踪法或质点跟踪格式。以往的许 多数值格式之所以不能精确地模拟锋面的情形，很大程度上是由于信息的不完 备。比如大多数格式都只在网格点上定义了物质浓度值，如果要反映一小团盐度 为1 0 的不扩散盐水在淡水中的纯平移运动，网格点上就不得不用到中间过渡值， 即介于0 到1 0 之间的值，所以不可避免地产生耗散。所以，质点跟踪法的处理 方法与一般方法有所不同，它利用了更多的信息，通过质点把这些信息“安全准 时”地从一个网格送到另一个网格。质点跟踪法的计算过程如下： 一、初始化： 初始时刻在每个计算网格内均匀放置n 1 个质点，给每个质点赋予初始场的物 质浓度。 二、每次循环的操作： ( 一) 滩地与水体的交界处，在由于转湿而重新参与计算的网格点内重新均 匀放置m 个质点，并根据周围情况赋予质点新的物质浓度； ( 二) 把上一步中所有非平流因素造成的网格点物质浓度变化传递给网格内 的各质点： ( 三) 根据流场计算所有质点的新位置； ( 四) 把越过网格边界的质点转交给新的网格，其中如果某网格内出现质点 数减少到0 的情况，则在该网格内重新均匀放置m 个质点并赋予上游物质浓度信 息，如果某网格内出现质点多达2 n 个的情形，则把2 n 个质点两两合并成n 个质 点，新质点的位置和物质浓度值均为原来两点的平均； ( 五) 对网格内的所有质点的物质浓度值作平均得到网格点的物质浓度值， 以供其它过程( 如扩散) 和下一步计算。 其中m 和n 都是自然数常数。 原则上每个网格内的质点数越多，计算精度就越高，所以m 和1 1 应该取得尽 可能大。然而事实上只要质点能大致均匀分布到网格内的各个区域，即能具有相 当的代表性，所以未必需要非常多的质点，而且质点过多的话，计算量也会变得 非常大。所以在综合考虑计算精度和客观需要的情况下，本文把每个网格内的质 点数限制在8 个以内，而初始时只设置4 个，即m = n - - - 4 。 如果说欧拉一拉格朗日方法是通过逆向跟踪质点来计算平流项的话，那么质 点跟踪法则是通过正向跟踪来完成这一使命。质点跟踪法避免了欧拉一拉格朗目 方法的反复插值过程，从而极大她避免了耗散。但是质点跟踪法也不可避免地会 遇到个别网格中出现质点数减少到0 的情况，这时就需要用上游的信息补充，从 而造成局地物质浓度变化比实际情况略快的情况。不过在很多算例中，这种可能 性都仅仅大约为1 ，所以只造成了约1 的物质移动加快，对计算的影响很小。 此外，一h 游位置的计算采用了绝对位置，所以避免了因。层的波动带来的误差。 至于网格内质点达到2 n 个而合并的措麓虽然会造成耗散，但是由于网格已经是 计算的最小单位，所以这种耗散是极其微小的。质点跟踪法和欧拉一拉格朗日方 法一样无法产生新的极值，从而避免了频散。 质点跟踪法对网格的要求相对比较严格。如果网格在水平方向或者垂向的尺 度变化过大，那么就会增加网格内质点数减少到0 或增加到2 n 的可能，从而使 物质移动加快的现象发生得更加频繁，或者数值耗散更大。 第三章各种格式在理想试验下的比较 为了测试并比较各种平流项数值格式，本文设计了两个理想试验。第一个试 验是一维的，它测试各种格式在简单流动下的耗散和频散；第二个试验是二维的， 测试各种格式在水深变化剧烈处的垂向耗散和频散。 3 1 一维试验 在首先进行的一维试验中，共有1 0 0 个网格( 也可以看成是线段) ，它们从 左到右排成一行，分别标号为网格1 到网格1 0 0 。所有网格的长度都是1 0 0 m 。给 定这些网格的初始物质浓度值场：网格卜4 9 的初始物质浓度值为0 k g m 一，网格 5 0 一i 0 0 的初始物质浓度值为2 0 k g m 3 。给定按时间变化的均匀分布的流场： 流速u = c o s ( n 2 4 0 0 ) ，其中n 为时间步，单位时间步长a t = 2 5 s 。 如此，便形成了一个人工的以4 0 0 步即1 0 0 0 0 s 为周期的往复流。假定该物 质浓度变化只受平流影响，不受扩散等其它因素影响。在上述往复流下，用平流 格式对物质浓度变化进行计算。使用不同的平流格式将得到不周的结果。本试验 每1 0 0 步输出物质浓度值场一次，以观察各种格式的表现。理论上完美的格式应 该是没有频散或者耗散的，即始终能保持浓度0 k g m 3 与浓度2 0 k g m - 3 的明显分界， 当中应该没有中间数值。 这个试验共测试并比较了5 个平流格式，它们分别是迎风格式、 l a x w e n d r o f f 格式、t v d 格式、欧拉一拉格朗日方法和质点跟踪法。由于是一维 试验，所以改进前后的欧拉一拉格朗日方法是一样的，故不作区分。试验结果显 示，迎风格式耗散极其严重，3 个周期后，0 k g m 。3 到2 0 k g m - 。3 的锋面几乎完全耗尽： l a x w e n d r o f f 虽然耗散控制得较小，但是频散十分明显；t v d 格式和欧拉一拉格 朗日方法没有发生频散，且耗散较小，表现优于前两种格式，其中t v d 格式耗散 更小；质点跟踪法几乎完美，十分逼近真解。所以可以看出一维试验下质点跟踪 法保留了突变的锋面，是5 个格式中最好的格式。 以下几幅图为各方法在该试验中的表现： r 掌 ； 兽 g 8 g 0 2 0 1 0 0 2 0 o 2 0 1 0 0 2 0 1 0 0 气 j 一 j 一 一 j j 一 j一 j一 1 一 1 1 一 1 t i m es t e p 2 0 t i m es t e p 2 1 0 9 t i m es t e p = 2 0 0 t n es t e p 。3 0 0 t i m es t e p 2 4 0 0 t i m os l e p 2 s 0 0 t i m es t e p = 6 0 0 t i r i cs t e p = 7 0 0 t i m es t e p28 0 9 t i t i cs t e p 2 9 0 0 t i i 雠s t e p = 1 0 0 0 t i m es t e p = 1 1 0 0 t i m cs t e p 2 1 2 0 0 02 0 3 0 柏5 06 07 08 09 01 0 0 g r i dn u m b e r 图3 i 迎风格式在一维情形下的表现 _ f 一厂 _ 一 j _ _ 1 j 02 03 04 05 06 07 08 09 01 0 0 o r i dn u m b e r s t e p = 1 0 0 s t e p z 2 0 0 s t e p = 3 0 0 s t e p = 4 0 0 s t e p = 5 0 0 s t e p = 6 0 0 s t e p t 7 0 0 s t e p = 8 0 0 s t e p = 9 0 0 s t e p = 1 0 0 0 s t e p = 1 1 0 0 s t e p 21 2 0 0 图3 2l a x w e n d r o f f 格式在一维情形下的表现 0。00000。033393。o。030。o 加0曲0加0m0心0罅m0m0m0加m 一，+量一t2_glou j厂 j _ j j _ _ - _ _ _f_ _ -r_ j j j t l i n es k d = 0 t i s t e p 】0 0 t i m es m 2 0 8 t i m cs t e p = 3 0 0 t i n es t c d = 4 0 0 t i m es i c 口= 5 t i i ms t e p = 6 0 0 n m es t e p = 7 d o t i m es t 舭9 0 0 1 伯口s 讳= 9 0 0 t i m es t e p = 1 0 0 0 t i m cs t e p 21 1 0 0 t i m e s k 口= 1 2 。0 02 0 0加5 06 07 08 a9 01 0 0 0 n dn u m b e r 图3 3t v d 格式在一维情形下的表现 mz o3 ( 34 0s 06 07 0鼬钧i g r i dn u m b e r s m p = 0 s t e p = 1 0 0 s k d = 2 0 0 s k d = 3 0 0 s m p = 4 0 0 s t e p = 5 0 0 嘶= 6 0 0 s t e p = 7 0 0 s t e p = 8 0 0 s d = 9 d 0 s 唧= l o 。o s u e p = i i 0 0 s t e p = 1 2 0 0 图3 4 欧拉拉格朗目方法在一维情形下的表现 如m o 加m o加mom omome如momomo加旧。加omomo 一、5睾一0鲁鲁8e8 0珀0m0加0如0曲”o曲m o如0加0如o加0加加0加m0 一，量一菪；皇_850u 3 2 二维试验 02 03 05 06 07 08 09 01 0 0 g v 6n t l m b e 图3 5 质点跟踪法在一维情形下的表现 s t e p = 0 s t e p = 1 0 0 s t e p = 2 0 0 s t e p = 3 0 0 s t e p = 4 0 0 s t e p = 5 0 0 s t e p = 6 0 0 s t e p2 7 0 0 s 口= s 0 0 s m p = 9 0 0 s t e p = 1 0 0 0 s t e p = 1 1 0 0 s t e p = j 2 0 0 这是一个水平一维、垂向一维的二维试验。其中水平方向的网格数量和大小 与前面的试验一( 一维试验) 相同，垂向采用。坐标，均匀地分为2 5 层。水平 网格也像试验一那样从左到右标号为1 到1 0 0 ，其中1 到4 8 号网格设定水深为 2 0 m ，4 9 号为3 0 m ，5 0 号4 0 m ，5 1 号6 0 m ，5 2 号8 0 m ，5 3 到1 0 0 号则设定为1 2 0 m 。 这样，就人为地构造了一个水深剧烈变化的陡坡，陡坡两端水深相差6 倍。给定 这些网格的初始物质浓度值场：网格卜4 8 的初始物质浓度值为2 0 k g m - 3 ，网格 4 9 1 0 0 的初始物质浓度值为0 k g m 。给定按时间变化的均匀分布的流场：在每个 网格的水深2 0 m 以上部分，u = c o s ( n 2 兀4 0 0 ) ，其中i q 为时间步；在每个网格 的水深2 0 m 以下部分，流速u = 0 ：单位时间步长n c = 2 5 s 。 如此，便在水深2 0 m 以上区域形成了一个人工的以4 0 0 步即1 0 0 0 0 s 为周期 的往复流，水深2 0 m 以下部分则没有任何流动，其中水平网格1 - 4 8 由于设定水 深为2 0 m ，所以在往复流区域中。在假定物质浓度变化只受平流影响，不受扩散 等其它因素影响的基础上，分别用各种物质输运平流项数值格式对物质浓度变化 进行计算，并每1 0 0 步输出物质浓度值场一次，以观察各种格式的表现。在实际 情况条件下，初始场中在水平网格卜4 8 的2 0k g m 3 浓度值将无法影响水平网格 。o。m描蛐驯郅驺鲫闶。o口oo。 4 9 1 0 0 水深2 0 m 以下部分区域，2 0 m 以下部分区域将保持0k g m 3 的浓度不变。 由于在水深变化剧烈的地方，大多数格式都很容易出现垂向明显的数值耗散，所 以这里重点关注各种格式的垂向数值耗散情况。 这个试验测试并比较了迎风格式、l a x w e n d r o f f 格式、t v d 格式、改进前的 欧拉一拉格朗日方法、改进后的欧拉一拉格朗曰方法和质点跟踪法共6 种格式。这 6 种格式在该