（安全技术及工程专业论文）球面上偏微分方程的数值求解研究.pdf

s t u d yn u m e r i c a l l y i np a r t i a ld i f f e r e n t i a le q u a t i o no f s p h e r e x u y a b o s u p e r v i s o r ：p r o f c h a n gj i a n z h o n g a b s t r a c t i nt h i st h e s i s ，c u b e s p h e r em e t h o di sa d o p t e dt os o l v ep a r t i a ld i f f e r e n t i a le q u a t i o no f s p h e r e t h i sa v o i d st h ed i f f i c u l t yo f n u m e r i c a l l ys o l v i n gp a r t i a ld i f f e r e n t i a le q u a t i o n i ns p h e r e c o o r d i n a t i o na n dh a sg r e a tl a t e n tc a p a c i t yi nn u m e r i c a lc a l c u l a t i o n sa c c u r a c ya n dp o t e n c y e s p e c i a l l yi th a sg r e a ta d v a n t a g e o v e r s p e c t r a lm e t h o d i np a r a l l e ln u m e r i c a lc a l c u l a t i o n f i r s tt h es i xf o r m so ff o r m u l a t i o no fs h a l l o ww a t e r e q u a t i o n sa n das u i t eo fn u m e r i c a l t e s tm e t h o d sa r eg i v e n t h e nc o n f o r m a t i o nm e t h o d o f c u b e - s p h e r e c o o r d i n a t i o na n dt r a n s f o r m r u l e sb e t w e e ns p h e r ec o o r d i n a t i o na n d f i g h t - a n g l ec o o r d i n a t i o na n dc o r r e s p o n d i n gd i f f e r e n t i a l o p e r a t o r a r eg i v e ni nt h et h e s i s o nt h i sb a s i s ，t h ef o r m u l a t i o n so fs h a l l o ww a t e r e q u a t i o n s a r e d e d u c e df r o mt h ec u b e s p h e r ec o o r d i n a t i o n , u s i n gt h em o d i f i e dd i f f e r e n t i a lq u a d r a t u r e ( m d q ) s p e c i a ld i f f e r e n c es c h e m ea n df o u r t ho r d e rr u n g e - k u t t as c h e m et os o l v et h es h a l l o ww a t e r e q u a t i o n s t h e v e c t o r v e l o c i t ye x p r e s s i o n s o fs i x r e g i o n s i nt h el i n e a ra d v e c t i o no f r o s s b y h a u r w i t zw a v e c a s ea n dt h el i n e a ra d v e c t i o no fac o s i n eb e l lc a s ea r eg i v e nf r o mt h e c u b e s p h e r ec o o r d i n a t e a n de x p r e s s i o n so fs i x i n t e r v a l s v e l o c i t y v e c t o ra r e g i v e n i n c u b e - s p h e r ec o o r d i n a t i o n a ne x p l i c i tp s e u d o - v i s c o s i c t yt e r m ( - v 4w a su s e dt o f r e et h e n u m e r i c a ls o l u t i o n sf r o mt h en o i s ei n t r o d u c e da tt h ei n t e r n a lb o u n d a r i e sb yt h ei n t e r p o l a t i o n p r o c e d u r e t h em e t h o di s ( 1 ) f u l l ye x p l i c i t ( 2 ) f o u r t ho r d e ro fa c c u r a c yi nt i m e ( 3 ) h i g ho r d e r o f a c c u r a c y i ns p a c e b yc h a n g i n go n l yo n ep a r a m e t e r ( 4 ) e a s yt oo p e r a t eo nm a s s i v e l yp a r a l l e l c o m p u t e r se f f i c i e n t l y , a n di ti sr o b u s ti np r o c e e d i n gn u m e r i c a lp r o b l e m t l l i st h e s i si st h ef i r s tt i m et oa d o p tc u b e - s p h e r em e t h o di ns o l v i n gn u m e r i c a l l yp a r t i a l d i f f e r e n t i a l e q u a t i o ni no u rc o u n t r y t h ew o r k so ft h et h e s i sp u s ht h en u m e r i c a lw e a t h e r f o r e c a s t sb a s ew o r kf o r w a r d k e y w o r d s ：s h a l l o ww a t e r e q u a t i o n s ，c u b e s p h e r em e t h o d ，r u n g e - - - k u t t as c h e m e ， p s e u d o v i s c o s i c t y , m o d i f i e dd i f f e r e n t i a lq u a d r a t u r e ( m d q ) ，n u m e r i c a lw e a t h e rf o r e c a s t 本人声明 我声明，本论文及其研究工作是由本人在导师指导下独立完成的，在完成论文时 所利用的一切资料均已在参考文献中列出。 姓名：徐亚博 签字：徐埘 腑年弓月谚日 善：立方球体坐标下每个区间的角变量 ，7 ：立方球体坐标下每个区间的角变量 h ：流体深度 h ：自由面高出参考面( 海平面) 的高度 饥：被流体淹没的山高 v ：水平的速度矢量 f ：c o r i o l i s 参数 v ：梯度算予 v ：散度算子 x ，y ，万，c ，d ：辅助变量 “：掌方向的速度分量 v ：r l 方向的速度分量 ：辅助变量的流动元素 z ：空间辅助变量 c 。，一，：加权系数 m ：相邻的网格点 w ：空间的离散因子 n ：时间步长 e 。，p 2 ：扩散项系数 五：经度方向 0 ：纬度方向 “。：经度方向的速度分量 v 。：纬度方向的速度分量 。：平流的风速 2 ：固体物质环绕旋转轴和球坐标的极轴之间的夹角 g ：重力加速度g = 9 9 0 6 1 6 m s 一2 a r ：地球半径a r = 6 3 7 1 2 2 1 0 6 m q ：地球的旋转率q = 7 2 9 2 1 0 - 5 s - 1 中北大学学位论文 1 、绪论 1 1 研究的目的和意义 所谓气象是指大气的状态、变化与现象”】( 冷、热、雨、雪等等) 。随着社会经济的 快速发展，一方面气象对工农业生产、群众的生活和工作、商业销售、旅游开发等影响 越来越大，气象与人们的关系越来越密切：另一方面，经济建设和群众生活对自然界， 尤其是对社会生态环境的破坏，使自然灾害、尤其是气象灾害发生频繁，对人类造成的 灾难日渐增多。 天气预报就是对未来时期内天气变化的预先估计和预告。“天有不测风云”，这句话 充分说明了天气预报的难度。随着科学技术的发展，天气预报的准确率在不断提高，人 们根据天气预报，可以适时安排生产和生活，使气象为国民经济建设服务，减少气象灾 害的损失。天气预报是根据大气科学的基本理论和技术对某一地区未来的天气做出分析 和预测，这是大气科学为国民经济建设和人民生活服务的重要手段，准确及时的天气预 报对于经济建设、国防建设具有非常重要的作用。 大气的行为可以用一组以数学方程式表示的物理定理来表达，由这些物理定理可以 计算大气的量或场( 如温度、风向和风速及湿度等) 将如何改变，因此如果我们可以解出 这组代表物理定理的数学式，我们就可以由目前的天气状态推衍出对未来天气现象的描 述，这就是我们所谓的数值天气预报。数值天气预报以计算机为工具【2 一，通过解流体 力学，热力学，动力气象学组成的预报方程，来制作天气预报。在数值计算球面微分方 程时，由于球面固有的曲率问题，阻碍了以往在笛卡尔坐标系中发展的差分格式和其它 方法直接应用到球面微分方程的数值求解之中，这也是进一步发展数值求解球面微分方 程方法的原动力之所在。 1 2 国内外研究状况 由于球面固有的曲率问题，使得数值求解球面微分方程时，选择适当的方法( 坐标 系) 来保证数值求解的精度和效率就成了一个关键问题。虽然表述球面现象最自然的应 是球坐标系，但由于极点的问题，在数值求解球面微分方程上就带来了许多很难克服的 障碍。如坐标的单值对应关系在极点不能保证( 经度变成多个值) ，而且由于方程中某 中北大学学位论文 些项的值在极点变成无限值而出现了奇异点，球面速度矢量在其它点能很好的定义，但 在极点上却不行。这些特性都阻碍了以往在笛卡尔坐标系中发展的差分格式和其它方法 直接应用到球面偏微分方程的数值求解之中，这也是迸一步发展数值求饵球面偏微分方 程方法的原动力之所在。近年来，气象学和大气动力学在数值天气预报上的发展也更推 动了这一问题的研究n 5 j 。 前人为了求解球面上的偏微分方程，提出过一系列的方法。第一种方法是保角变换 的方法。该方法使两个相交的曲线之间的夹角保持不变。但是它并没有获得很普遍的应 用，主要是因为它的复杂性和费用i s 。在6 0 年代术，人们又尝试着使用几何分解的方法， 把整个球体表面分解到几个区问【7 柏i 。这种方法去除了几何的奇异点。但是，由于内部 边界处的数值的不稳定性和本身固有的复杂性，这种方法很快就被放弃了。在7 0 年代， 谱方法【10 1 和统一的经纬度嘲格技术【1 1 1 的进一步发展起来。目前，谱方法成为数值天气预 报应用的主要方法。但由于l e g e n d r e 变换【1 玛的计算所需时间过长和在计算中出现的非 物理现象的趋势f 1 3 j ，谱方法的计算量大大超过了有限差分和有限元方法，人们已对其在 数值预报中所占据的主导地位提出了疑问。另一方面，由于所使用的坐标系统的结构有 关的问题，使的该方法在精度和效率方面仍有限制。虽然可以采用对高纬处的不稳定的 短波进行过滤的方法来放松对时间的限制，但是同时这种过滤程序又严重的降低了在所 有纬度处的数值解的精度【1 4 1 。 如前所述上世纪六十年代有人提出用球面的几何分解和坐标变换的方法i 侧，求解球 面微分方程，但由于方法的过于复杂和数值计算的时效问题，未能达到成熟应用的水平。 在此基础上，上世纪末c r o n c h i 等人提出了“c u b e ds p h e r e ”方法求解球面微分方 程。立方球体技术作为一种新的有效的求解球面偏微分方程的方法，它是将球体表面分 成六个完全相同的区间，这些区间由球的外切立方体的中心投影丽获得。每个区间都可 以由一族等圆心角的大圆的弧所覆盖，这样就在每个区间定义了个准均的网格，且 在坐标系内没有奇异点。相应的计算网格是由角坐标瞄，叩，定义的六个完全相同的矩形 “f 一三7 z 网格所形成，每个区间跨越4q ，从而避免了在普通球坐标系中数值求解偏微分方 程时固有的许多困难。在六个区间的边界处，用网格组合技术的差值程序，将六个区间 2 中北大学学位论文 连接起来，可在较小的区间搭接量下获得稳定的数值解。该方法既避免了球坐标上数值 求解微分方程的难点，又在数值计算的精度和效能方面有很大的潜力，尤其是在并行数 值计算方面有谱方法不可比拟的优势。但在与数值求解格式的匹配与人工粘性的添加处 理方面和实际应用处理问题的经验方面，还需加以大量研究，在这些问题解决之后，可 以预计“c u b e ds p h e r e ”方法，在数值求解球面微分方程上将有重要的位置。 w i l l i a m s o n 给出了一套评价数值求解球面偏微分方程的算例”鲥8 】，即求解球面上 浅水波方程的一套算例，每一种算例都代表了在数值求解球面流动中所遇到的典型问 题，通过这几种算例可以检测出一种数值模型在求解球面流动中的能力，通过求解复杂 程度依次增加的几个算例，可以判断出该数值求解的精度及迅速检查所编程序的有效 性。计算流体力学的发展为数值求解偏微分方程提供了许多有效的数值求解格式和方 法，计算流体力学作为一个专门的学科是2 0 世纪6 0 年代形成的n ”2 1 ，它一出现就得到 了迅速的发展，解决了流体力学中的许多疑难问题，最引入注目的就是计算流体力学为 求解各种问题而发展的许多有效的求解方法和格式，不同的格式都是针对具体的物理问 题而设计的。本文在研究分析计算流体力学的一些格式之后，选择m d q 方法和四阶r r k 方法求解浅水波方程，它具有完全显示格式，时问上的四阶精度，仅仅通过改变一个参 数就可以获得空间上的高阶精度，意欲在并行的计算机上有效地执行，而且该数值方案 有较强的健壮性等特点。 1 3 本文的主要研究工作 正如前所述研究数值求解球面微分方程时，首先应选择合理的求解问题即测试方 案，其次应选择合适的表述球面的坐标系，以避免极点上的奇异和矢量方向不连续问题， 再有对具体的测试方案选择合理的数值求解格式。基于解决以上问题的考虑，本文首先 给出浅水波方程的几种形式，同时也给出4 个算例：其次选择用c u b e s p h e r e 方法来表 述球体几何上的偏微分方程：再有用m d q 方法和r r k 方法进行数值求解。用上述数 值求解方案对上面4 个算例中最有代表性的两个算例进行数值模拟。具体推导出浅水波 方程在立方球体坐标系中的形式和各个有关的参数在该坐标系下的形式，并分析模拟结 果，最后给出结论和建议。 总之本文的工作是以理论分析为指导，数值计算为基础，在前人已有成就的基础上， 中北大学学位论文 探索和研究数值求解球面偏微分方程的更为有效的新方法，为数值天气预报的发展打下 了更峰实的基础。 4 中北大学学位论文 2 1 引言 浅水波方程及其测试算例 在检测一个数值计算方法时，首先要确立该数值计算方法要解决的具体问题即所要 求解的具体方程。本文研究的是球面偏微分方程的数值求解方法，目的是为数值求解天 气预报等绕地球流动的数值模拟提供新的更有效的数值求解方法。在地球表面，包固地 球的大气主要集中在1 0 1 l k m 的对流层内，占地表约7 l 的海洋深度也不超过6 k m 上 下，与6 3 7 1 k m 的平均地球半径来比，大气层或者海洋带实在只是围绕地表的一个薄层， 其垂向运动速度比水平分速度小的多，一般情况下其垂直加速度也很小，从大尺度海气 运动而言，它们都是局部的( 微) 小尺度运动。因此绕地浅水波方程正恰好描述了如上 特征，因此该方程成了人们检测个求解球面微分方程的目标方程。本章将给出浅水波 方程的几种形式及其具体的四个算侧，为后面检测本文的数值模拟方案提供具体的求解 问题。 2 2 浅水波方程及其几种形式盼2 4 l 浅水波方程可以从对不可压缩流体的e u l e r 方程进行水深平均，利用垂向压强按 静水压强分布，自由面和海底起伏不随时间变化等条件，精确地推导出来。 h 为流体的深度，h 是自由面高出参考平面( 海平面) 的高度。如果 。代表被 流体淹没的山高，则h = h + k ，v 为水平的速度矢量，是c o f i o l i s 参数，g 是重 力加速度常数，t 是时间，霞分别为经度，纬度和向外方向的单位矢量。v 是梯 度算子，v 是散度因子，所，m 。是度量系数，距离增量讲。 通量形式 浅水波方程的动量和连续方程： 警+ v 坼矿) 一盛h v - 卵 ( 2 1 ) 中北大学学位论文 百a h * + v ) = 。 上面的通量形式的浅水波方程在球坐标系中可以写成如下形式 警m ) f + u 口t a h 口卜+ 面g h 瓦o h = 。 警m ) + f + u a t a h 臼) 如+ 譬等= 。 其中v ( ) ；去击( ) + 丢刍( ) v 一志c o s 0 良蜊0 0 d l 觑l c o r i o l i s 参数为f = 2 q s i n 0 平流形式 平流形式的水平动量和质量连续方程可写成如下形式： 警= 一露肌g v 厅 丝+ 矗v 矿：o d t 上面的方程在球坐标系中可写成如下形式： 孚c t + v v u - ( f + u 口t a h 口 v + 志嚣= 。 詈m v ( 卜u 口t a h 口+ 詈嚣= 。 百o h * m v 九面h 旧o u + 筹) = 。 其中鲁( ) z 言( ) + 缈v x ) 涡量散度形式 水平动量可用相应的涡量和龄序安亮诀 6 ( 2 2 ) ( 2 _ 3 ) ( 2 4 ) ( 2 5 ) ( 2 6 ) ( 2 7 ) ( 2 8 ) ( 2 9 ) ( 2 i o ) ( 2 1 i ) r 2 1 2 ) ( 2 1 3 ) 中北大学学位论文 涡量定义为f = k ( v v )( 2 1 4 ) 水平散度为j ；v 矿 ( 2 1 5 ) 应用下列矢量斌妒铲沙= v ( 半卜矿 ( 2 1 6 ) 可直接把涡量引进动量方程中去，得到 水平动量形式如下：詈= 一g + ，皿x v _ v ( g + t v o v ) ( 2 1 7 ) ol| 在球坐标中写成如下形式： 詈= g + 办一志刍阿三0 2 ) 坳 害砘m 丽1 f 1 g h + 圭o2 + v 2 ) 1 ( 2 1 9 ) 旋度算子被如下定义启( v 矿) = 夏岳| 言一百a u c o s 0 0 si ( 2 2 0 ) d cl 似 d i 应用旋度算子重v ( ) 和散度算子v ( ) ，动量方程变成如下形式： 坌a t = 一v g + i ) v ( 2 2 1 ) 警缘v 岍y ) v 卅h 半 ( 2 z z ) 在球坐标系中： 等一磊b 刍眙+ 办卜磊品眙+ 咖c o s 占】 ( z ) 警= 志刍眙+ ，m 一a c o l s o 刍眙+ ，k c 。s 口卜v 2 b + 圭0 2 + v ：) 1 c z 埘) 删= 志掣+ 志刍( 掣) 流函数。速度势形式 通过引入流函数和速度势z ，可以替代速度变量，速度和这两个标量变量有如下 姜絮 中北大学学位论文 v = 丘v + v z 速度的两个球面分量可以用如下式子表述： 。：一土丝+ j 一塑 。：l 塑+ 上亟 将旋度和散度算子带入以上两个速度分量中得到如下绝对涡量和散度。 绝对涡量即s f + ，= v 2 + 厂 散度艿= v 2 z 流函数和速度矢形式的动量与质量连续方程可被写成f 2 5 j ： 署+ v 慨) 一，白，；f ，) = o 警一v y ) 一，白，z ) = - v 2 ( k + 曲) 鲁+ v g ) 一j 缈) = 。 其y 水平流函数，z 速度势 正交曲线坐标系中的形式 g ，_ y ) 是正交曲线坐标系中的坐标，m ，为拉梅系数( 或称度规系数) ， 则相邻点的微分长度讲可表示为：( 讲) 2 = m ：出2 + m ：咖2 速度矢量v 在x 和y 方向的分量形式： u ：m 一! y ：棚生 动量方程： 警一 ，+ 去( 矿誓一唔肛砉挚 警+ ，+ 去( 矿誓一噜牡考挚 中北大学学位论文 连续方程： d h d t + 啬 鼽吼扣矿) ：0m 。朋yl 缸”7 移”7 j 对于球坐标， x= y = 口 r t l ，= a c o s o 笛卡尔形式 在笛卡尔坐标0 ，y ，= ) 下定义的速度矢量矿= ( x ，r ，z ) 7 笛卡尔坐标系下浅水波方程组中的动量方程为： 百c o v + c v + q 2 q + + j ) = o 其中c = q 7 s ( k ) 口= 9 r 2 3 9 ) ( 2 4 0 ) ( 2 4 1 ) f 2 4 2 ) ( 2 _ 4 3 ) ( 2 4 4 ) f 2 4 5 ) r 2 ，4 6 ) ( 2 4 7 ) f 2 4 8 ) ( 2 4 9 ) a 一砂 矿一 + a 一苏 旦 + a 一西 i i d 一曲 其 丝出 口 | ；塑出 d 口 1 1 | i 出一西砂出 一一 m 所 y = l l 晰 w p 科一出盯一瑟彪一瑟 橱咖订一砂凹一砂 翁一苦订缸彪一融 iii， z z z + + + y r r + + + 如舡舡 ，。l ，一矿 i i 口 9 q 7 占= 2 。o 。zli：一y：一ix茎y c z s 。， q 7 卢= p v ( h ( 2 5 1 ) 户= 乒 ；蔓x 三三羹：三i c 2 s 2 ， v 加( 罢，考，笔) 7 偿s ， 警p v “v 。儿。 ( 2 5 4 ) 其矩阵p 为任意的球体上的点g ，y ，z ) 正切的c a r t e s i a n 矢量 2 6 1 。 余弦钟每一次都绕着地球转一圈，规定了几种余弦钟的不同对流风场，包括环绕赤 1 0 中北大学学位论文 道，环绕南北极和在两种方向之间偏转一定角度的对流方向， 其对流速度由下式给出2 7 】 “= d o ( c o s l 7 c o s o t + s i n o c o s 2 s i n a l 以避免流动的对称性。 r 2 5 5 ) v = 一“os i n s i n 口( 2 5 6 ) u 和v 分别为经度五和纬度臼方向的速度，口为极地轴和余弦钟环绕轴之间的夹角。 “。2 天) 2 4 0 r e s ( 2 5 7 ) 平流余弦钟的初始形状由下式给出 以，p ) ：( h o j 2 x 1 + c o s ( , - ! ! 矿r 胄 (258)0fr 、。 ， 、。吖 其中h 。= 1 0 0 0 m ，为任一余弦钟范围内球面上的点以，口) 与余弦钟中心之间的大圆 距离。，= a a r c c o s s i n o cs i n 0 + c o s o cc o s o c o s ( 2 一以) 】，( 2 5 9 ) 其中q 。，扫。) = 0 , r 2 ，0 )( 2 6 0 ) 半径r = a 3 ( 2 6 1 ) 球带状稳定非线性地转流动 该方案对非线性的浅水波方程来说其解是一个稳定状态的。它由固体物质的旋转或 由有一定高度的绕地旋转的带状流动组成。c o r i o l i s 参数是经度和纬度的函数，所以描 述流动的球坐标的极点没有必要与地球的旋转轴一致。 其初始速度由下式给出 ”= u o ( c o s o c o s 口+ c o s 2 s i n o s i n t r )( 2 6 2 ) v = 一u o s i n2 s i n a f 2 6 3 ) 高度场由下式给出 咖= 硪一卜+ 卜c o s 2 c o s o s i n a + s i n o c o s t z ，2 偿e 4 ， 其中驴2 2 天) 4 0 m s ( 2 6 5 ) 中北大学学位论文 幽o = 2 9 4 x 1 0 4 m 2 s 2 ( 2 6 6 ) 其c o r i o l i s 参数为：f = 2 f 】( - c o s i c o s o s i n a + s i n o c o s a ) ( 2 6 7 ) 孤立山峰的带状流动 该方案由t a k a c s 利用来研究数值方案的整体的守恒性t 2 创，它由算例2 的带状流动 流过一个孤立的山峰。计算的风场和高度场都与第二种方案一致， 其中口= 0( 2 6 8 ) h o = 5 9 6 0 m “o = 2 0 m 5 ( 2 6 9 ) ( 2 7 0 ) 山高饥= h , o ( 1 r r ) ( 2 7 1 ) 其中 2 2 0 0 0m ( 2 7 2 ) r 2 州9 ( 2 7 3 ) ，2 = m i n k 2 ，0 一丸) 2 + p 一幺) zj ( 2 7 4 ) 中心取为以2 3 州2 ( 2 7 5 ) 眈2 州6 ( 2 7 6 ) r o s s b y h a u r w il ：z 波 r o s s b y h a u r w i t z 波是球体上非线性正压涡量方程的分析解洲。虽然不是浅水波方 程的分析解，自从p h i l l i p s 第一次作为测试方案以来刚】，它在气象学中经常被使用， 虽然不同的研究者选用的参数不同，但它实际上已经成为一个用来评估球体上的数值方 法的性能的标准测试方案。 初始风场是散度为零的，由以下流函数给出： i f f = - - a 2 w s i n p + 口2 芷c o s r 臼s i n 口c 。s r i ( 2 7 7 ) 这里w , k , r 是常数，h a u r w i t z 表明这一流动的形状。在以下角速和散度为零的条 中北大学学位论文 件下从西到东移动，其形状不会改变。 v ：r ( 3 + r ) w - 2 f 2 ( 1 + r ) ( 2 + r ) 速度分量如下： “，= a r w c o s 口+ a r k c o s 8 1p 伍s i n 20 一c o s 2o ) c o s r 五 v 。= 一a r k r c o s 8 10 s i n g s i n r 2 其w = k = 7 8 4 8 x 1 0 “j 一1 r = 4 高度为： ( 2 7 8 ) ( 2 7 9 ) ( 2 8 0 ) ( 2 8 1 ) ( 2 8 2 ) 砂= 劝。+ 胛4 ( 口) + 甜2 b ( o ) c o s r 2 + a r2 c ( o ) c 0 8 2 r a ( 2 8 3 ) 一p ) = 詈( 2 q + 国c o s 舢4 kc o $ 2 r 口k + 1 ) c o ：口+ 妇2 一r 一2 ) 一2 r 2 c o s - 2 口】( 2 8 4 ) b p ) 2 黼c 。s 。秽( r + 2 r + 2 ) 一( 胄+ ) 2 c 。s 2 目】 ( 。) c ( p ) = kc o s 2 r0 【( 胄+ 1 ) c o s 20 - ( 尺+ 2 ) 】 ( 2 8 6 ) 该算例的数值求解结果可以与给出的由谱方法计算而得到的结果进行比较。 24 太音，i 、结 本章通过分析绕地球流动的大气和海洋的情况，选定绕球旋转的浅水波方程为本文 求解的目标方程，即作为测试求解球面微分方程的数值方法的测试对象，并分析了浅水 波方程与无粘流动不可压缩欧拉方程的关系，给出了浅水波方程的六种形式，并介绍分 析了四个具体的浅水波方程的测试算例，为后面的工作提供了具体的求解对象。 中北大学学位论文 3 、c u b e - s p h e r e 技术 3 1 引言 由于球面固有的曲率问题，使得数值求解球面微分方程时，选择适当的方法( 坐标 系1 来保证数值求解的精度和效率就成了一个关键问题。虽然表述球面现象最自然的应 是球坐标系，但由于极点的问题，在数值求解球面微分方程上就带来了许多很难克服的 障碍。如坐标的单值对应关系在极点不能保证( 经度变成多个值) ，而且由于方程中某 些项的值在极点变成无限值而出现了奇异点，球面速度矢量在其它点能很好的定义，但 在极点上却不行。这些特性都阻碍了以往在笛卡尔坐标系中发展的差分格式和其它方法 直接应用到球面偏微分方程的数值求解之中，这也是进一步发展数值求解球面偏微分方 程方法的原动力之所在。上世纪六十年代有人提出用球面的几何分解和坐标变换的方 法，求解球面微分方程，但由于方法的过于复杂和数值计算的时效问题。未能达到成熟 应用的水平。在此基础上，上世纪末c r o n c | 珏等人提出了“c u b e s p h e r e ”方法求解球 面微分方程。立方球体技术作为一种新的有效的求解球面偏微分方程的方法。目前看来 能较高精度进行数值天气预报的数值模拟方案，它标注球体的坐标系和微分方程应满足 的六个标剐3 2 l ：1 ，图形必须免于任何的奇点。2 ，使求解方程保存有流体力学的一般形 式。3 ，定义在球体表面上的物理网格必须接近于规则的等距离网格。4 ，赤道应该是网 格系统的一条坐标线。5 ，所采用的数字方法与谱方法相比在计算时间和精确度上都要 有优势。6 ，所使用的数值算法应该是可升级的，并且应该允许在大量并行的计算机上 运算。本文在前人工作的基础上，经过深入的研究选择了c u b e s p h e r e 来作为标注球体 的坐标系，该方法既避免了球坐标上数值求解微分方程的难点，又在数值计算的精度和 效能方面有很大的潜力，尤其是在并行数值计算方面有谱方法不可比拟的优势。 3 2c u b e - s p h e r e 描述 我们采用的“c u b e s p h e r e ”技术来解决球体几何上的偏微分方程。该方法是在上 世纪末c ，r 0 一等人提出的。该方法是一种基于把球体分解成六个完全相同的区间的一种 新的网格技术，可以通过把球体的外切立方体投影到球体表面而获得。用这种方法定义 1 4 一一 生! ! 奎堂堂垡笙塞 一 ，- - _ 一 一 的六个部分避免了任何因坐标系引起的方程的奇异点，并且允许六个部分使用相同的数 值计算网格。正如图3 1 所示，两簇等距的大圆弧的交集定义了一个网格区间，这里大 圆的每个弧都可以联系着纵坐标或横坐标线。重复图3 1 中的办法六次，做出适当的选 择获得相同尺寸的区间，就可以得到精确的完整的球体表面( 如图3 2 ) 。从圈3 2 中很 显然看出，以这种方式定义的坐标系统是非垂直的，只有在赤道处圆弧才以9 0 。相交。 在图3 3 中的每个区间中，角变量f ，叩的范围为i 孚，手i 。该方法与以往的方法相比， l 叶叶j 最主要的创新就是：该方法可以给出坐标转换法则和相应的张量规则，能推导出新坐标 系下要求解方程的形式，免去数值计算时计算每个区间的张量元素，并且可以充分利用 六个区间的对称性。通过网格组合技术和合适的搭接区域插值方法控制数值计算得不稳 定性，特别是在六个区间的计算平面内数值网格是简单正交网格。通过与合理的数值离 散格式配合可使得我们获得在整个球体表面上的稳定和精确的数值解答。 愈饿 图3 1 两簇等距的大圆弧的交集定义的一个网格区间图3 2 立方球体的概图 i i v f 7 。 l 图3 3 立方球体的转换图 中北大学学位论文 为了使用立方球体坐标系的变换法则和微分算子，我们定义了下列的辅助变量。 x = t a n e )( 3 1 ) y = t a n b )( 3 2 ) j ：1 + x 2 + y 2 ( 3 ，3 ) c = 0 + z 2 y ”( 3 4 ) d = ( 1 + 】，2 y ”( 3 5 ) 六个区间中球坐标，直角坐标与c u b e s p h e r e 技术构成的坐标之间的坐标变换规则 如下： i ( 赤道) ： z = 么= t a n ( 3 ，6 ) y = 么= 秽。 ( 3 7 ) ，= g 2 + y 2 + z z y 2( 3 8 ) 鼢匕啪x y 8 v 凇2 4 , 。， i i ( 赤道) ： x = 一影y = 1 t a _ i l 声 ( 3 1 0 ) y = z y = 1 t a n s s i n ( 3 1 1 ) ，= g 2 + y 2 + z 2 y 2( 3 1 2 ) i i i ( 赤道) ： x 2 y x = t a i l 妒 ( 3 1 3 ) 】，2 一z x = 一1 u m o c o s ( 3 1 4 ) ，= e 2 + y 2 + 。2 严 ( 3 15 ) i v ( 赤：道) ： 1 6 中北大学学位论文 x = 一x y = 一1 t a n 矿 ( 3 1 6 ) y = 一z f y = 一l t a n o s i n 矿 ( 3 1 7 ) r = ( x 2 + ) ，2 + z 2y 2 ( 3 r e ) v ( 北极) ： x = y z = t a n o s i n 庐 ( 3 1 9 ) y = - x l z = - t a n o c o s 妒 ( 3 2 0 ) ，= b 2 + y 2 + z 2 ) 忙 ( 3 2 1 ) ( 挣由b 蒜2 ) ( 刳 扭z z ， v i ( 南极) ： x = y z = 一t a n o s i n o f 3 2 3 ) y = 一x z = - t a n o c o s ( 3 2 4 ) r = 0 2 + y 2 + z 2y ，2 ( 3 2 5 ) ：妨f 一似型黑p 1 ( 3 、 1 4 j 。矿矿l c y 一凹加水儿4 j 2 6 ) 这里弘，0 ，r 】为标准球坐标的经度，余纬和沿半径的坐标。x ，y ，z ) 为c a r t e s i a n 坐 标4 = 4 f + 彳。= a o e 。+ a f t 为球体上的普通矢量。在每个区间都有一套单元基本 矢量k ，；，j ，屯f = a ，= ，a ，= l 。 可以看到在所有的赤道区域内，都有相同的转换矩阵。所有的坐标转换都免于任何 的奇异点，因此在某种意义上说，六个区间之间就没有任何的差别，因此c u b e s p h e r e 方法可以充分利用六个区间的对称性。实际上，每个给定坐标系统的区域都可以由其他 区域在c a r t e s i a n 坐标轴上经过严格的旋转得到。在所有区间的张量都是相同的，为 1 7 一一 生i ! 奎堂堂垡鲨塞 f 1 一x y c d 0 1 g 2 i - 驯c d 10 j ( 3 2 7 ) 1 0 0 1 j 这一事实暗示出在每个区间内所有的微分算子都有相同的表达式。在每个区问内， 任何给定的p d e 都有相同的表达式，因此大大简化了编程和计算过程。 对于矢量的运算有，给定两个矢量a = k ，qj 和6 = 虹，岛，6 ，j 灿m 巾。坞b t + 口r 6 ，_ 茜眠乜如) ( 3 2 8 ) 砌2 嘉【c d ( 口，6 ，q ) + 灯- ，- a c b ) e ；+ 古【矾。印一) + c d 以弦。+ 筹k - a 。b 。， 。2 捆丁微分算于甭，给定标量厂皓，r l ，) ，矢量矿g ，r l ，) 则v 5 ；( d 面e f + 罟苗p + 三r f t 些c 笪a 掌+ c 面o f 一+ 等a ， ( s s 。) v 2 卜q 艿 c _ 2 叭oa o f 手卜苦( 豺器鲁愕甜善) 慨。z ， v 一丽 3 2 专 嘉) + 筹茜+ 砉杀) c 。一。， 弘阽；( 劳专r 比一罟专，_ + 警等k i l ( 万c d 万0 ，_ 一茜专，一警等k ( 3 。， + _ 胆- - l 旧c - d t h d 丝0 1 7 吉等) 一吉等+ 吉鲁一 在半径为口的球体上，沿f 和叩坐标方向的无限小的直线位移和相应的表面付穆南 峨= 业8 豢c o s 2 f 1 8 ( 3 3 4 ) 中北大学学位论文 ( 3 3 5 ) ( 3 3 6 ) 从上面的式子，很明显的看出，在每个区间的中心处，网格都是规则的，越往边界 处其变化就越大。最小的面元素( 在r = 0 ，善= 州4 或，7 = 州4 ，f = 0 时获得) 和中 心处( 玎= 孝；0 ) 的面元素的比值为2 2 * o 7 1 。在区域的四角处，该比值增加为0 7 7 。 在c u b e s p h e r e 方法中，网格组合，搭接区域及插值要求都并不十分苛刻，不同区 域的值连接并不要求非常光滑连接，即并不要求其导数连续。 从上可知，c u b e - s p h e r e 方法避免了球坐标系中以往求解球面微分方程的问题，如极 点、奇异、无界、矢量方向等难点，且计算网格为简单正交网格，可以方便的使用各种 计算流体力学的数值求解各种形式。 3 3c u b e - s p h e r e 技术中相邻区域的搭接及其节点的插值 当数值求解一个复杂区域的偏微分方程时，一个很可能的选择就是把原始区域分解 成有着简单结构的几个子域。这种方法的一个显著优点就是每个子域都可以被设计来遵 守特定的边界。但网格分解的方法有一个重要的问题就是如何把相邻的不同区域之间连 接起来，在差分格式中边界上的点也需要计算微分值【3 3 。5 1 。当处理多个坐标系统时，便 会出现如何得到与边界点相邻节点的值的问题，因为这些点并不被当地的坐标系所覆 盖。但由于设计了一个坐标系统覆盖了整个区域，那么这个点就成了相邻区域坐标系覆 盖的内部点。因此这个点的值可以由在其相邻区域的坐标系内插值而得。相邻坐标系被 连接在一起时，只有在边界处熏叠。然而，考虑到数值计算的稳定性和精确性，在复杂 网格之间有定的熏叠量是很有必要的。在实际中，重叠的数量只需考虑几个网格量即 可，因为在重叠点处的数值插值和计算，对计算时间和计算量上来讲是一个特别大的负 担。在插值时，插值计算的精度只要和采用的数值差分格式的精度相同就足以保证球面 数值模拟的精度和收敛。 在立方球体技术中，构成整个球体表面的六个区间的坐标线都是大圆的弧，因此， 1 9 击希 堕巧尘妒 勰 中北大学学位论文 一个给定区域的坐标线沿着一个方向延长与相邻区域的坐标线可连接在一起。如图3 4 所示。两个赤道区域的垂直坐标线都是由从北极到南极的大圆的弧组成，因此沿着赤道 方向是连续的。在这里我们仅考虑这样的情形，在每个区间内，在两个方向上都采用相 同的网格宽度，即n 。= n f = n ，n ，和f 是叩和f 方向的网格点数。我们采用这样的 方法，不是由于立方球体技术或可分解的网格方法而是为了使算法的效率达到最优的 需要。实际上，通过运用相同规则的网格，可以充分利用区域之间的对称性，这大大减 少了不同区间边界处网格相互连接时所需的计算量。 在图3 4 中，通过在两个区间内选择相同的网格宽度，重叠区间垂直网格线完全重 合。因此，为了在重叠区域内获得由善= k 4 + i a ) ( f - n 。，o + 。l ，n 。为重叠 区域的起始网格的节点编号，为网格间距，眠为重叠区域的网格数) 定义的坐标线 上的所有节点的值，只需沿垂直叩方向的一元差值即可获得。即= ( 叫4 + 弘) ， ( ，= n o ，n o + n s 1 ) 定义的坐标线上的节点值，也只需一元差值即可。 对于虬= 以= n 的网格，每一区间的水平坐标和垂直坐标的所有重叠点都是相邻 区域的内部网格点。因此，如果重叠只是局限于区域之间的边界线，对每个区间，只需 在4 n ，l 一2 ) 个节点上执行一元插值。 两个区间之间的重叠数有一个几何限制，因为一个区域的中心点是相邻的四个连续 坐标系统的奇异点。因此，为了保证精确性，n o + n 。 n 2 这严格条件必须保证。 虽然要求重叠的程度在l n o n 2 一n 。范围内，但是在n o = l 的情况下，同样可以获 得稳定的精确解，这种特性和在节点上只需一元插值，都大大减少了连接六个区间所需 的运算量【3 6 j 。 与其他的边界技术相比，c u b e s p h e