（流体力学专业论文）热带太平洋赤道附近数值模拟与海温同化研究.pdf

摘要 本文研究了利用m i t g c m 模式和最伉插馕法相结合，同化热带太平洋券道附 近的海温数据的方法。主蒙工作有以下几个方面： 翁先，介绍了m i t g e m 模式的基本方程、离散方法帮代码结稳。目隧简单奔 缨了阈化方法的意义、必要性和阅化的发展，详细分绍了最忧捶僮方法的基本原 理和优点。 其次，采取了剿用奇羿值分解法解决方程中的瘸态现象，和蔫逐点遂次插佼 法，掇赢嗣德效果，改进桷关性模型等方法来宠善最忧捶值法。著利用将热繁太 平洋赤道附i 匠的固定观测站点的观测资料和m i t g c m 模式计算值，进行同化试 验，骏证了方法静霄效性。 褥次，农本论文的基谶上，对建立的同化平台可以改谶蛇地方傲了说明。 关键词：最优插值m i t g c m 模式 热带太平洋 a b s t r a c t t h i st e x tc o n s t r u c t sa i la s s i m i l a t i o ns y s t e mw i t hm i t g c mm o d e la n do p t i m a l i n t e r p o l a t i o nm e t h o dt ov e r i f yw i t ho c e a nt e m p e r a t u r ed a t aa r o u n dt h ee q u a t o ri n t r o p i co c e a n t h ec o n t e x to f t h et e x ti n c l u d e st h r e ea s p e c t s ： f i r s t ，i ti n t r o d u c e st h eb a s i cd y n a m i c se q u a t i o n s 、m e t h o d so fd i f f e r e n c ea n d c o d es t r u c t u r e i ta l s os i m p l yi n t r o d u c e st h em e a n i n ga n dn e c e s s i t yo f a s s i m i l a t i o na n dt h ee v o l u t i o no f a s s i m i l a t i o nm e t h o d s 。t h eb a s i ct h e o r ya n di t s s u p e r i o r i t yo f o p t i m a li n t e r p o l a t i o nm e t h o da r ee x p l a i n e de m p h a t i c a l l ) ： s e c o n d ，s o m em e t h o d sa r ee x p l o i t e di nt h eo it oi m p r o v et h ef i n a lr e s u l t ， s u c ha ss v d 、s t e p w i s ei n t e r p o l a t i o n 、t h eu p s w i n go fr e l a t i v 盼m o d e l 。a n a s s i m i l a t i o ne x p e r i m e n t a t i o nu s i n go b s e r v a t i o nd a t aa r o u n dt h ef i x e dp o r ti nt h e e q u a t o ro f t h ep a c i f i co c e a n a n dt h er e s u l to f m o d e li sd o n e i tf a l l so u tt h a tt h e a s s i m i l a t i o ns y s t e mc a nr e f l e c tt h ed i s t r i b u t i o no ft h eo c e a nt e m p e r a t u r eb e t t e r a n dc a nb eu s e dt of o r e c a s te f f e c t i v e l y t h i r d ，t h et e x ti n d i c a t e ss o m ea s p e c t st h a tc a l lb ei m p r o v e do nc o n t i n u a l l y k e yw o r d s ：o p t i m a li n t e r p o l a t i o nm e t h o dm i t g c m t h et r o p i cp a c i f i co c e a n 4 第一章绪论 1 1 论文的研究意义 热带海洋区域一直是海洋科学家重点关注的西太平样暖池和台风的发源地， 也是流经我国东海的著名的黑潮的起源海域。西太平洋暖池不仅是全球海气相互 作用中最典型的海区之一，而且与其相邻的其他大洋所生成的台风等也经常侵袭 我国东南沿海乃至内陆地区并造成重大的经济损失。因此，开展对热带西太平洋 暖池变化的研究以及毗邻洋区的海洋观测具有重要的科学意义，对于预测我国大 陆的气候波动也具有重要的现实意义。 1 9 9 7 年爆发了百年来最强的一次热带中，东太平洋表层海水温度大范围增 温事件一厄尔尼诺，给世界许多国家造成了严重的自然灾害【2 】，我们国家也受 到严重的影响吲。厄尔尼诺事件最典型的基本特征是赤道东太平洋大范围海 水的异常增温1 6 1 ”，而其先兆现象如赤道偏西风的爆发和次表层暖水的东移等却 发生在西太平洋。研究人员通过对热带太平洋海洋变量异常的分析，以期望获得 厄尔尼诺事件发生的原因与必然联系。同时，通过建立海洋数值模型，来努力进 行数值模拟和准确预报，来减少人类的损失。 随着科学技术的发展，人类观测海洋的手段也越来越多。卫星观测、建立固 定观测站点、施放浮标、利用船舶收集等为人类提供了越来越多的观测资料。但 是相对于占地球表面7 0 的海洋来说，通过现有的观测资料难以达到正确认识海 洋的目的。现有的海洋数值模式也无法准确的描述复杂多变的海洋运动。同化方 法 8 j 是通过一定的数学模型和优化标准，将不同空间，不同时间，采用不同手段 获得的观测资料有机的结合在一起，建立相互协调的分析或预报优化系统，确定 那些不能直接观测的量，以及没有观测到的地方的相关信息，同时模式本身也可 以得到优化。因此利用数据同化技术可以充分的利用仍然不很充分的海洋观测资 料，加深对海洋动力和热力过程的认识和理解。同化技术在海洋和气象研究上有 着广泛应用的发展前景。通过同化技术将数值模式和观测数据结合起来，能够获 得更准确的数值计算结果和更精确的数值预报结果。 1 2 本文的工作 本文结合m i t g c m 模式帮最饶舔德方法建立了一个有效酶海溢数据目耗平 台。利用赤道附近的固定站点的观测斑料，对其进行试验。试验绺果表明，同化 结果有效的兼顾了数值模拟值和观测德，纠正了模式模拟值出现黝误差情况，更 好鲍反映了赤遂辩近的海滋场的势京特铤，可戳有效躲进行海澄数据数值颟掇。 1 2 1 主要内容 本文的主要内容分为以下几个方面； 1 、第一章是主要简单介绍了本文研究的意义。阐述了研究热带太平洋赤道 瓣避熬意义，势霭革说明矮篾强化技零黝必要瞧。 2 、第二章简单的介绍了m i t g c m 数值计算模黧。介绍了模烈的基本方程、 离散方法以及模式的结构特点。 3 、第三搴阐述了最毯矮篷夔基本琢灌。逶逮对足静数据露纯方法夔毙较， 说明了最优捕值法的优点，并详细的介绍了最优插德原理。 4 、第四鞲试验结果的分析与比较。通过对几个闷内的赤道附孟匠温度的模式 诗算结栗、阏纯缝象襄褒溅结栗豹遴孬 l 较，分援了该系统兹套效毽窝实露瞧。 5 、第五章展望。概括了本文，给蹦了本文的绪论。 1 2 2 本文的主要特色秘鲢裁点 1 、根据海洋流动的特点和观测站点的分布特点，使用椭圆形数据相关性。 2 、对予计算过程中出现的病态矩黪，采用奇异镶分解法( s v d ) ，以获得精确 解。 3 、对预报误差协方藏矩阵和观测误差协方差矩阵，每半个月爨新次，以 获想更为可纛豹精瘦。 4 、锋对诗算数据量避大导致计算内存不够帮计箨不稳定的情况，采用逐点 逐次插值法。 6 第二章m i t g c m 模式介绍 本文采用的是美国麻省理工学院的m i t g c m 模式m 1 ( m i tg e n e r a lc i r c u l a t i o n m o d e lo 该模式具有一下几个方面的主要特征： ( 1 ) 该模式具有一个统一的核心动力学方程，既可以用来模拟研究大气物理 运动现象，也可以用来模拟研究物理海洋的运动过程。参见图( 2 1 ) 。 ( 2 ) 可以用来模拟研究非静态动力过程，对于大小尺度的动力过程均可以模 拟。 ( 3 ) 使用有限元体积法来处理非规则边界和非规则几何体，生成正交曲线网 格。图( 2 2 ) 。 ( 4 ) 利用该模式可以进行灵敏度分析和最优化研究。 ( 5 ) 模式可以有效的在不同的平台运行等。 2 1 删t g c m 模式的核心动力学基本方程 由于m i t g c m 模式对大气物理运动的计算和海洋物理运动的计算采用的是 同一个核心动力方程( 参见图1 - 1 ) ，所以采用同构方法来表达其核心方程( 在模式 中，海洋变量和气象变量的对应关系如图1 - 3 所示1 。 m i t 模式的动力学基本方程： 2 1 1 水平动力方程 盟+ ( 2 n d t 舌) + v 妒= r v h ， ， 2 1 2 垂向动力方程 华+ j c ( 2 n i ) + 掣+ b ：0 戚 o r ( 2 1 ) ( 2 2 ) 2 1 3 连续方程 v 碗+ 要= o 泖 2 1 4 状态方程 b = b ( o ，s ，r ) 2 1 5 位温方程 岛一。 百一啪 2 1 6 盐度方程 鲁呱 ( 2 3 ) ( 2 4 ) ( 2 5 ) ( 2 6 ) 其中r 是垂向坐标，旦d t = 鲁+ 市v 是全导公式，v = v h + 喏是梯度算子， v h 表示水平方向的梯度，昙表示垂直方向的梯度( i c 是垂向矢量) ，f 是时间， 哥= ( “，v ，) = ( 瓦，) 是速度矢量，妒是压力重力势，西是地球旋转速度，b 是 压强梯度力，口是位温，s 是盐度，l 是压力项和i 的耗散项，3 。是压力项和 位温目的耗散项，同样3 。表示压力项和盐度s 的耗散项。 2 2 海洋模式边界条件 为了封闭海洋动力基本方程，必然用到初始条件和边界条件。这些边界条件 包括： 在，= r 。= 一h ( x ，y ) 即海洋底面 o d l h a t m o s p h e r i c p h y s i c s 图2 1 s t r y n1 1 i n n “* n 田 图2 2 t i - , n , o 盱p tt ( 】3 9 z - pi s o m o r p h i s m 趴e 排l i c “) o 幽i t a t e s ) z pa t i i 删p l 僻雌f pc o o r d i n a t e s ) d 占+ r ”：+ 里：p = p h 昏 ( 证+ f 。+ 已雷 g p + a 7 p = op # q+ 屯小= o 玉“! a = - = 0 w m l 川d 。：o 却= q 0 c l f b = ( 邮一ss ”qd ：q = s a ，i lf 里r nl m x = p e n h p 。 a 。r i ! 、! = o r 胡胡1 肛一p i w 可j 在，= 如。= 玎即在海平面 图2 3 一鲁+ 【“象+ v 努j 西【苏却j 其线性形式为： 。：塑 西 固体侧边界：法向速度为零 “h = 0 ( 2 7 ) ( 2 8 ) ( 2 9 ) f 2 1 0 ) 亓是固体边界的法向量。 在海洋中我们假定：( 1 ) 高度：，= z ，( 2 ) 垂直方向的速度为：，= 面0 2 = w ， ( 3 ) 压强：矽= 旦，( 4 ) 浮力：b ( o ，s ，) = g - - - ( p ( o ，s ，r ) 一n ) 。其中肛是海水参 p cp c 照密度，g 为重力加速度，如。，。= r + r l ，r o ( x ，y ) 是海平面的平均高度( 通常 假定为r o = 0 即r m o v 。= r 1 ) ，叩是超出或者低于海平面平均高度的值。 2 3 静态动力项、准静态动力学项、非静力动力项 首先将压力西分解成海表压力，静态压力和非静态压力 妒【x ，y ，) = 允( 工，y ) + h 【x ，y ，) + 九( x ，y ，)( 2 1 1 ) 将方程( 2 1 ) 写成以下形式： 盟3 t + v h 睡+ v + 巳h v h 谚o 2 嚷 ( 2 1 2 ) 同时有以下公式： 孥：一6( 2 1 3 ) 鲁+ 等= q ( 2 1 4 ) 其中毛。是非静力参数，( 龟，g ) 分别代表动力方程中不同方向的平流项，量纲 项，耗散项和科氏力项之和。在球坐标系下，它们具有以下的形式： 瓯一加。一 等一半 _ 一2 n 伊瑚咖s 卅l ( 2 1 5 ) g ，：一v v v - ( f f 一竺堕) 一 一2 q “s i n 妒) + r ， ( 2 1 6 ) q ；_ v v + 兰旦) + ( 2 n “c 。s 四+ r p( 2 1 7 ) 其中，r 表示距离地球地心的距离，妒代表纬度，一哥v 。，一v v ，一v v ，分别 代表各个方向的平流项，一 堡u v t a n q ，一 v i u 2t a n 伊 ， 兰三_ 上) ，分别代表各 个方向的量纲项，一 - 2 q vs i n 妒+ 2 q ，c o s 讲，一( - 2 f l u s i n o ) ， 2 q c o s 纠分别代表 各个方向的科氏力项，f ，r ，r ，分别代表各个方向的耗散项。 2 3 1 游态动力项和准静态动力项组成部分 在静力学方程中，公式( 2 1 5 - - 2 1 7 ) 变成以下形式： g 。：一哥v 。一 一u v t a n o 一 2 f 2 v s i m p + f 。 娃 鼠。一v v v - 一u 2 t a n g o 一 一2 f 2 u s i n 尹) + r ， a & 一0 日表示玩球的平均半径。 在凇静力学方程中，公式 瓯洳。一 等一u v t ，a n 妒 一 ( 2 1 5 2 1 7 ) 变成以下形式： ( 2 1 8 ) ( 2 1 9 ) f 2 。2 0 ) - 2 f 2 v s i n 妒+ 2 f c o s f a f 。( 2 2 1 ) g ，：十v ，一芦一u 2 t a n o 一 2 妨村s i n 四十r ， rr g ： 生芝) + 2 q “c _ 0 s 妒 公式( 2 。1 4 ) 变成以下形式： 婺l ：2 q “c o s 伊 o r 2 。3 2 滞静力学形式 q ：一节v 。一 堡一u v t a n e 一 - 2 n v s i m p + 2 n c o s 妒) + r 。 ， 瓯：母v ，一芒u 2 t a n e 一 - 2 f 2 嚣s 遗癣+ f ， r， g ：十v ，+ 生旦 + 2 n u c o s g , + 0 r ( 2 ，2 2 ) r 2 2 3 ) ( 2 。2 4 ) ( 2 2 5 ) f 2 + 2 6 ) f 2 2 7 ) 1 2 2 4 压力求解 2 4 1 静态压力求解 静态压力的求解是通过对方程( 2 1 3 ) 在垂直方向上积分得到的 r 挚= 蛎 ? = r a 咖 最终有： 蛎( w ，r ) = r 6 咖 2 4 2 海表压力求解 r 2 2 8 ) r 2 2 9 ) 拇回压力足遇趣料连续刀栏( 2 3 ) 垂且积分得到嗣： e ( v 瓦垣e ) a r = o ( 2 3 0 ) 所以有： 署w v 呷+ e dv h 瓦扣o ( 2 3 1 ) 这里有叩= 如。，。一r o ，也可以运用莱布尼兹法则变换得到如下计算公式： 署钾一”死拈舳眦e ( 2 3 2 ) 2 4 3 非静态压力求解 将方程( 2 1 4 ) 两边求散度，对方程( 2 1 3 ) 两边求微昙，然后代入连续方程( 2 3 ) 中 d ， 我们就得到如下方程： v ；如= v 岛一( v ：癣+ v 2 ) = v 户 ( 2 3 3 ) 2 5 耗散项求解 2 5 1 动量方程中的耗散项求解 b = 碱2 v + 4 窘+ 锅4 v r 2 3 4 ) 其中，4 ，4 分别表示水平涡粘性系数和垂向涡粘性系数。4 是双谐摩擦水 平系数。 2 5 2 轨迹方程中的耗散项求解 a t ，s = v ；! ，( r ，s ) + 墨v ：( r ，s ) 其中丝是扩散张量。 j 2 6 海洋模式的动力学基本方程 2 6 1b o u s s i n e s 假定和不可压缩假定下的海洋动力学基本方程 鲁+ 五哪去v 扩于 o e 。d 肪w 。+ 百g p + 去老= l v ：吒+ 8 w = 0 o z p = p ( 曰，s ) 鲁喝 导电 d 。 ( 2 3 5 ) f 2 3 6 ) ( 2 3 7 ) r 2 3 8 ) r 2 3 9 ) r 2 4 0 ) ( 2 4 1 ) 2 6 2 “半压缩”b o u s s i n e s q 假定下的海洋动力学基本方程 不可压缩假定是指在流动和密度方面都是不可压缩的，识是很多情况f ，我 们在计簿的过程中，需要考虑到密度变化的影响，所以对密度p 做以下处理： p 2 p o + p p g 2 p e p p ( o ，s ，p o ( z ) ) - p o ( 2 4 1 ) r 2 4 2 ) ( 2 4 3 ) 假p o 是不变化的，密度的大小变化主要体现在p 上，风为参照密殿。这 耱馕瑟一f 海洋动力学鏊本方程裁交瓷驻下形式： 鲁+ 矗时1 见v , p = 于 等+ 菪罢+ 去鲁= r ，瓦+ 菪磊+ 瓦嵩r w v ：魄十娑= o o 暑 p = p ( o , s ，p o ( z ) ) 一风 鲁喝 鲁喝 2 7 嫩标变换 f 2 4 4 ) ( 2 4 5 ) ( 2 4 6 ) ( 2 ，4 7 ) f 2 4 8 ) r 2 4 9 ) 童予焱海洋系统豹硬究过程中，缀多情凌下瑟在球坐标下( 参觅蚕1 - 4 ) 求鼹 变量，丽“ v ，w ，梯魔算子，散度辣子等在球坐标系下具有其特殊的表达式： d z 舻n 0 8 鲈瓦 ，：，堡 d t f 2 + 5 0 ) f 2 5 1 ) 肼 d t v ；1 8 1 面0 ，刳 v 一上r c o s e p 翥+ 杀泓。嘲 专警l 蕊却、，。毋 2 8 差分方法及算法 2 8 1 时问步上的差分 f 2 5 2 ) f 2 5 3 ) 控5 0 ) 穰式审稷努豹动力方程包括溷令媛缀方程( 群葶委v ，湿度0 ，蕊度s ) 懿三兮 诊断方程( 1 4 ，密度压强梯度力p t b ，压强画。) 。模式通过一定的计算方法来 保证在满足诊断方程约束的情况下，向前推进预报方程中的变霪，以获得最终的 结果。 2 8 1 1 压强项的差分 2 8 1 1 1 刚盖假定压强格式 水平动力方程和连续方程可以写成以下的形式： 香g o ，簪= 瓯 2 5 1 ) a ，v + g o 。叩= g v ( 2 5 2 ) 蚤：“+ 艿，v + o ：w = 0 ( 2 。5 3 ) 将式子( 2 5 3 ) 在深度方向积分，同时运用刚盖假定( 在刚盖处，w = 0 ) 和边界条 件，变成以下形式： 8 ；壬菇+ 矛，饼= 0 2 5 4 ) 其中，脚= e “砌是“在深度方向上的积分，脚一l 毗足v 在深度方向上 的积分。 则方程( 2 5 1 ，2 ，5 2 ，2 5 4 ) 可以在时间和空间上离散成以下形式： 1 6 图2 4球坐标示意图 五为经度，妒为纬度，为距离地球中心的位置 w 羔，u 1 7 “”“+ a t g o 。叮”t 。= “”+ f q ”川心 v 肿1 + a t g c 3 。呷”制= v ”+ “z “卅坨 a 。，描”“+ a 。f 茆”“= 0 f 2 5 5 ) f 2 5 6 ) ( 2 5 7 ) 将公式2 5 5 和公式2 5 6 代入公式2 5 7 就可以得到一个求解卵“椭圆方程，可以 将公式2 5 5 2 5 7 做以下安排： “= 矿+ q ”“2 v + = v ”+ r d ”“2 佗5 8 ) r 2 5 9 ) o ，a t g h o 。r l ”“+ a y a t g h o y 叩”“= a ，上磁+ a ，日 ( 2 6 0 ) u ”1 = 村a r g o 。矿+f 2 6 1 ) v ”+ 1 = v + - a t g o 。1 7 ”+ 1( 2 6 2 ) 公式2 5 8 2 6 2 的顺序执行求解的过程就是模式中运用的刚盖假定下的压力方 法的原理过程，图2 1 描述了压力方法在时间步上的求解过程。 t i n t e n a t ( n + 1 ) a t f - ；u ( 1 1 v + ) 州n u f v 8 丫2 q 2v ( ：) j 。什。 图2 - 5 2 8 1 1 2 线性自由表面的隐式压强格式 刚盖假定过虑了外部重力波的影响，并因此改变了r o s s b y 波的频散关系。 椭圆方程的离散形式具有一些特征值为零，这就潜在地使得方程本身难以处理或 者求解效率降低。将刚盖假定用自由表面地线性化方程来代替，则公式2 4 可以 写成： 0 , r + 0 x 日舀+ 0 y 爿蚕= p e + r ( 2 6 3 ) 其中p e + r 是淡水源项。在时间上对公式2 6 3 进行离散得到： 刁“+ a t o ，。，m ”卅+ & t o 。舟移”= 叩”+ r ( p e + r )( 2 6 4 ) 将公式2 5 5 ，2 5 6 ，2 6 4 按照一定地次序排列分解成以下： ”= u ”+ a t q ”2 ( 2 6 5 ) v = v ”+ a t g _ ：”+ 1 7 2 ( 2 6 6 ) r = ( 玎”+ 尸一e + 五) 一& t o ，m + + f a y m ( 2 6 7 ) 蛐t g h b x r l “+ o y & t g h 0 7 1 - 等一 ( 2 6 8 ) “= & t g o ，r l ( 2 6 9 ) v “= v + & t g o 。矿“( 2 7 0 ) 为选择开关( 当= 1 表示刚盖假定的压力方法，当= o 表示隐式线性化自 由表面的压力方法) 。对方程2 6 5 2 7 0 顺序执行求解的过程就是模式中隐式线性 化自由表面的压力方法的工作原理。 2 8 1 1 3 显式a d a m s b a s h f o r t h 格式 对于压力强格式中的q “”，q ”“，我们采用的是准二阶 a d a m s b a s h f o r t h 格式来处理，这种时间步的处理方法适合动量方程和轨迹方程 中的任何显式项。假设f 代表预报方程中的任一变量甜，v ，0 ，s ，则有： q ”1 ”= ( 3 2 + 毛日) 畔一( 1 2 + ) e 。( 2 7 1 ) 这是一种向前外差的方法。当e a 。= 0 时公式( 2 7 1 ) 具有二阶精度，但是在计算的 过程中容易产生不稳定因素，所以通过对颤。赋值来保证该方法的稳定性。 2 8 1 1 4 隐式后差格式 对于垂向扩散项和粘性项，可以对它们运用隐式后差格式。方程如下： f “1 一a t a ，k a ，f ”1 = f ”+ r q ”2 ( 2 7 2 ) 对于式中的g “可采用公式( 2 7 1 ) 1 拘1 方法来处理。对于公式( 2 7 2 ) i h 以改写成如 下形式： f = 矿+ f d 7 2 ( 2 7 3 ) f ”1 = 1 ( f + )( 2 7 4 ) 其中，1 为算子l = 【1 + a t o r t co r 的逆算子。 2 8 1 1 5 时间同步格式 当所有的状态变量在时间上同步的话，那么，a d a m s b a s h f o r t h 显式外推 差分法就可以应用到压力方法中去。 睇，s = g o ，s ( ，伊，s ”) ( 2 - 7 5 ) g 器“2 = ( 3 2 + ) g ：，s 一( 1 2 + 矗口a n ，- i ( 2 7 6 ) ( 口，s + ) = ( 口”，s ”) + f g 嚣“2 ( 2 7 7 ) ( 矿“，s ”1 ) = 鹾。( 矿，s + ) ( 2 7 8 ) 铂= j 6 ( 矿，s ”) d r ( 2 7 9 ) 留= 岛( i “，硝b ) ( 2 8 0 ) 印删= ( 3 2 + ) 它一( 1 2 + ) 留一1 ( 2 8 1 ) 哥= 矿+ 厢：n + l 2 )( 2 8 2 ) 矿= 三- 1 ( 矿)( 2 8 3 ) 叩+ = 占店( 矿+ j p e + r ) 一a t v 月舌 ( 2 8 4 ) v g h v r f l + l _ 丛a f t = 丢 ( 2 8 5 ) 矿“= 矿一a t g v 矿“ ( 2 8 6 ) 图2 - 6 展示了所有状态变量在时间上同步的计算过程。虚箭头表示运用的是 显式a d a m s b a s h f o r t h 外推差分，实箭头表示的是从n 层值对”+ 1 层进行估值。 2 8 1 1 6 斜压交错网格格式 霞秀分鼷豹原溺，内部夔力渡虿戆在霹瘸步上受翔一定豹袋铡， 这种情况下，将热力举变量和流变嫩在时间步，t 交错开来更为有效，这样在 本质上允许内部重力波在时间步上变化以满足= 阶精度和稳定性。相对于时 滴丽多丽言，交舞步豹关键变健在予对静态蓬力豹求簿，将公式( 2 。7 努交为 如下： 颤= f 6 ( ，沙 7 ) _ 举茁，s l 威s ；p 纛如丰，1 lu ，3 g ? i j 弋fg n m f 一_ g ! 均 u , v n - l 1 龄“l 矿疆i v ? 等； 图2 - - 6 则交绩时间步的计算顺序如下： 硪“= g e 。s ( 扩，8 - u 2 , s “) 础= ( 3 2 + ) 喇”一( 1 2 + g a b ) 睇坨 ( 痧，s ) = ( 移”，s 4 ) + 矗f 磺 f 2 ，8 8 ) f 2 8 9 ) f 2 ，9 0 ) 2 1 0 “+ 1 2 , s + j n ) = 嗣，( 九s ) 掰“= 舷扩“2 ，s ”“2 涉 g ：= g ；( i ”) 窿“”= ( 3 2 + s a 8 ) 霹一( 1 2 + c a # ) 霹一 事= 矿+ 嫡2 - v 蜊n + l7 2 ) = 掣( 矿) 邛+ = ( 坪”+ p 一嚣十露) 一a t v 日季 v g 日v 零“一寺gn + l = 一立a t 2 口”“= i 一a t g v t ”+ 1 图2 7 展示了_ 交错时间步的计算过程。 1 1 1 1i 2 ) at( r ) - - h ) a t n + ) a t t i m e 秘一1 a t n a t n + 1 ) a t g 3 i 嗽一g ：爹卜g 罐 0 , s n - i a ，“也一 参9 s ；0t,s“+ i i 啦h ； g 黔一- g 嚣 v 1 夕” 一 车uv t 1 1 + i 图2 7 ( 2 9 1 ) ( 2 ，9 2 ) f 2 9 3 ) ( 2 。9 4 ) ( 2 9 5 ) ( 2 9 6 ) f 2 9 7 ) ( 2 9 8 ) f 2 9 9 ) 一 、 一、 、 。+ 一 一 一一一一一 ，叫 n u v g u 2 8 1 2 非静态动力颈的差分 菲静态确力学方程渗及弱垂离动象方程，对予j 静态动力顼霭矮求解一个三 维的椭圆方秘。对于静态渤力项，我们仍然采用对其进行垂向积分并求解二维椭 圆方程的处避方法。 动量方稳在对阔步上瓣离散魏下： l f u n + l + g d 。+ a 。= 1 + 瑚( 2 1 0 0 ) 一& 1v n + l + 9 8 y 艿y 2 岔1i 口 , n + 7 耸( 2 1 0 i ) 古w + 。，簖1 = 石1 + 础删( 2 1 0 2 ) 以上方程必须满足连续方穰，连续方稷离数翔下： t 甜”“+ 0 。v ”“+ 0 ，矿轧= o( 2 1 0 3 ) 幼量方程中变量的估计值计算如下： 帮= 扩+ a t 0 7 ” f 2 。1 0 4 ) v + = v ”+ f g ? “7 2( 2 1 0 5 ) = 矿出域删2 在计算u n + 1 ，v ”1 的时候，我们引入一个中间步； 甜“m “”一a t 8 ：蝣1# ”= u 。一a t g g ，l “ 稼。1 0 6 ) ( 2 1 0 7 ) v “= v ”一a t 3 ，簖1 v ”= v 4 - - a t g o y r i ”1 ( 2 1 0 8 ) 将以上方程代入公式( 2 6 3 ) 中，并与公式( 2 1 0 3 ) 联合求得近似方程： 色h 8 匙矿+ 笳1 ) + 8 y h c 3 ，( 册”l + 嚣卜等寺 ( 2 l 。9 ) 流动必须满足连续方程和遵续方程的熏向积分方程，用以上两个约束条件来形成 求解蕺1 蕊三维糖瑟方程； a 。谚嚣1 十a w 谚嚣1 + a ，谚等1 = a ，“”+ a ，v ”+ a ，w ( 2 1 1 0 ) 整个求解过程按热下方程蹶序执行： “= “”+ r q ”“2 v + = v ”+ f g ( 川7 2 w + = w ”+ f g ? 2 卵+ = e 正( r ”+ p e + 月) 一a t 0 ，m + a 。月矿 e l g h 8 川n + l + a y g h 0 7 1 一- a 扩t - - - - 广= 一百1 7 甜”= “+ a r g o ，矿“ v ”= v - a r g o 。r l ”+ 1 a 。簖1 + a w 簖1 + a ，簖1 = o x 甜”+ a y v ”+ a ，w + ”1 = “”一a t o ，簖1 v ”1 = v ”- a t 0 ，簖1 a ，w ”“= 一a j ”“一o y v ”+ 1 2 8 1 3 自由表面的差分 自由表面的差分方法包括：非线性格式，隐式c r a n k - n i c h o l s o 格式，显式分 裂格式。下面介绍以下c r a n k - n i c h o l s o 格式。 在= 0 即不包含非静态动力项的情况下，自由表面的椭圆方程如下： s 瑁“1 一v a t 2 ( r o r m d ) v 以，7 ”1 = r ( 2 1 2 2 ) 刁= 产a t v 皂；咖+ ( 尸咽“( 2 1 2 3 ) ：；a t v 以 ( 2 1 2 4 ) 在= 1 即包含非静态动力项的情况下 颤矿1 - v a t 2 ( r o r 胁d ) v 趣矿1 = r ( 2 1 2 5 ) 叩= 啄r - a t v 。- 免；咖+ f ( j p 咽“( 2 1 2 6 ) m 哟 四 m 四 岣 聊 坳 聊 卿 哪 l 1 l 1 l l 1 1 1 l e 口 伫 q 亿 亿 亿 亿 亿 伍 亿 ( v ：+ a 。编“= v 。- 州，+ ) ( 2 1 2 7 ) i 一 v = v 一矗h a t v 九“1 ( 2 1 2 8 ) c r a n k - n i c k e l s o n 正压格式 完全隐式格式是无条件稳定的，但是却过虑了快重力波，导致了位势能的损失。 改进的方式是通过将隐式部分( f l ，) 和显式部分( 1 一，1 一，) 有机的结合起来 ，y 分别代表海表压力梯度和正压流发散。 当= y = 1 时，是全隐式格式，当= ，= 1 2 时，是能量守恒，无条件稳定的 c r a n k - n i c k e l s o n 格式，当( ，) = ( 1 ，0 ) 或者( f l ，) = ( o ，1 ) 时，分别代表只有在小的 时间步上才能稳定的前差格式和后差格式。 ；+ = ；”+ r 瓦”“2 + 一1 ) a t v h 幼”+ a t v h 磊w 7 2(2129)(fl1 ) a t v ha t v h1 2 9 )v = v + r q +一 ” 饵”比 ( 2 r + = 啄矿+ 姒p 呐一a t v 。皂【y “( 1 训；”妙 ( 2 1 3 0 ) 罟+ v 瑚矿1 + ( 1 一f 1 ) 棚+ e h v 砧“= 吾 ( 2 1 3 1 ) 古+ v 一良【叮”1 + ( 1 一矿】+ 一矿12 古 ( 2 华氓t 昕”1 + ( 1 训强= ( p 咽( 2 1 3 2 ) 如果是静态动力方程的情况下( ：0 ) ，我们按照以下方法计算r n + l 和；”1 呷”“一v f l r z x t 2 b a & 一哆，d ) v r ”1 = r + ( 2 13 3 ) ：；+ f l a t v 良( 2 1 3 4 ) 2 8 2 空问上的差分 动力方程空间上的离散是用有限元体积法来实现的。也就是说在流场内部 用的是格点法( 也称谓二阶中心有限差分法) ，同时允许在边界上使用不规则网 格。我们对水平方向和垂向方向分别进行不同的处理。 公式中符号的意义以及算子表达公式： 缸，缈，分别代表x ，】，r 方向上网格间距。 4 ( d = “，v ，w ) 平面在垂直于0 方向上的投影面积。 圪，圪，分别表示包含“，v ，w ，0 点的网格体的体积。 ，j ，k 分别表示对应于x ，y ，r 方向的网格点的下标标识符。 4 中= 巾，+ j ，2 巾，2 。 巾。= ( 巾，2 + 中i - 1 1 2 ) 2 。 中= 二一点巾。 1 血。 审o = 瓯巾，占。卿梯度算子表达式。 可- y - = 去 巧坝+ 6 ：x f y 散度算子表达式。 审2 中= 一v 审。拉普拉斯算子表达式。 有限体积法用来对方程在空间上进行离散，“有限体积”实际上有两方面的 意义：一方面表示该方法可以对边界进行不规则处理，另外表示可以对像振动等 这样的不平滑现象的进行处理的非线性插值法。 有限元体积法是基于流场中的流在任何时候都是不分叉的基础上被用来推 动动力方程在空间上向前运行，也就是说流的分量在空间上是交错的( 参见图 2 4 1 图2 8 图2 8 展示了流的三个方向”m w 上交错的情况，从图中可以看出，水流在u 方向上流经的面处在v ，w 方向上流经的面之间，它们之间是相互交错的关系。 图2 9 展示了网格在水平方向上的交错关系。 图2 1 0 展示了网格在垂直方向上的交错关系。 2 8 2 1 连续方程和水平压力梯度项 即+ 瓦1q i 3 0 1 1 。，t 石c h 。4 中- = 瓯一击4 中： 即+ 毒t 詈j 。叩+ 轰t 吒咆一毒t 中： ( b w + 古。嚷中_ ) = q + 酽一击瞑巾： 巧妙g q 九”+ t 哦嘶嚏v + 8 k a 。w = a , s k ( p e ) “ 其中，蝇，蚬，峨，哦参照图2 - - 9 ，a r c ，q 参照图2 1 0 ( 2 1 3 5 ) r 2 1 3 6 ) ( 2 1 3 7 ) 陀1 3 8 ) 0q 0 图2 9 jf 了 一o 一 l j 廿一 一g 一 v 2 8 _ 2 2 静态动力学项 图2 1 0 j 仁叶， 一o 一 1 = 小 一o 一 e 一 垂向动量方程中有静态动力学项和准静态动力学项， 在z 坐标系下，静力动力学项被离散成： c p + g l + 毛驰h = ”( 2 1 3 9 ) 2 8 2 3 流量形式的动力方程 g u = g 苫“+ g ：4 + g ：。9 + g i 。”+ g ：。c 七g 0 一”“ g ，= g + g 警r + g ：叫枷+ g 了d m + g ：。r “+ g 警m 曲慢 g w = g 0 + g ? + g ：m “+ g ：4 “+ g ：”“+ g 0 一”“ 2 8 2 4 动力方程中的平流项 平流算子在空间上是二次精度的： a 。m 如警= 6 0 t 丽七8 帚l 破1 + 6 雾t 蛩 a 5 rr h ；g = 6 霜i 可t + 6 妒j - _ + 6 帝i 可4 4 g = 4 盯访。+ 玩旷订+ 皖矿。订 其中的u ，v ，w 定义如下： u = a y , 6 r ：h 。u v = 缸g 心r h w = 一w 2 8 2 5 科氏力项 4 奶九6 ：o = f a 声：h y m f a c f h 矽 a , a r ：h ，g 7 r = 一，a c 心f h 罗 a c r c g 警= nf j a c 心f h f 其中科氏力参数，和，定义为 厂= 2 q s i n 缈 ，。= 2 q c o s 伊 r 2 1 4 0 ) ( 2 1 4 1 ) ( 2 1 4 2 ) r 2 1 4 3 ) f 2 1 4 4 ) r 2 1 4 5 ) f 2 1 4 6 r 2 1 4 7 ) f 2 1 4 8 ) f 2 1 4 9 ) f 2 1 5 0 ) ( 2 1 5 1 ) f 2 1 5 2 ) r 2 1 5 3 ) 2 8 2 6 动量方程中的向量不定式 由于球坐标系f 的方程出现量纲项，使得方程表达不具有普遍性。利用向量 不定式方程可以方便的在不同的正交曲线坐标系( 如球坐标，边界流和等角球立 方体系统) 下转换。 非静力学向量不定式方程表达为： a ，矿+ ( 2 q + f ) a i 一6 尹+ v b ；v f ( 2 1 5 4 ) 该方程是正交曲线坐标系下的通用公式，b 为b e r n o u l l i 方程，f = v a i 为涡 矢量。 对于向量不定式方程的离散方法与流量形式的方程是一样。唯一有区别的是 g 。g 的表达式如下： q = 印+ 货”+ g ”+ 四，。+ 霹7 + 哦一“9 + q 一如9 ( 2 1 5 5 ) g v = g 拿+ g ：9 + g + 声+ g ：一+ g ：”。s ”+ g l 。s s 4 q 1 5 6 ) g w = g 0 + g 譬+ g 铲+ g 鲁b + g ：矗蛳p + g i 。啦p q 1 5 7 ) ( 1 ) 相对涡度 相对涡度的垂向组成部分被显性计算，显性离散。 相关涡度定义为： 炙2 石f2 毒( 坝v 一缸“) ( 2 1 5 8 ) 其中，4 是垂向方向的涡度单元的面积，f 是该单元的流通量。 ( 2 ) 动能 动能用k e 表示，定义如下： k e ：昙( 7 + + ) ( 2 1 5 9 ) ( 3 ) 科氏力项 线性科氏力项的位涡度拟能守恒形式定义如下： 钟= 忐吾雨7 亿， g 净毒吾雨7 ( 2 蚓) 科氏力参数被定义在涡度中心点。 非线性科氏力项的位涡度拟能守恒形式定义如下： 弘击吾丽。 抛， g 昝一毒吾雨 亿s ， 科氏力项也可以用绝对涡度表达式+ 来求得： 四+ 啡”= 瓦1t f + ：3 - - m 。v f 。 ( 2 1 6 4 ) 础科一击孚菰 亿， 在求解高阶平流尺度的时候，绝对涡度和相对涡度是不相同的。绝对涡度的 单调平流形式和相对涡度的单调平流形式是不同的。 ( 4 ) 剪切应力项 印”2 赤万( 耻口w ) ( 2 1 6 6 ) 嘶”2 赤万( 驴t w ) ( 2 1 6 7 ) 掣2 击删瑚) 2 专删。+ k e ) f 6 ) 7 j ( 平耗散项 肚4 1 - - ! - - ( 4 a y 7 “+ 缸嘻v ) q h - d i s s i p 一- - 壶4 ( 如。一如4 d ) 一面( 吩f 一 。f + ) ( 2 1 6 8 ) r 2 1 6 9 ) r 2 1 7 0 ) f 2 1 7 1 ) 3 2 口一p = 面击4 瞳( 4 f 一4 f +