（计算机应用技术专业论文）多尺度非达西热流耦合并行算法研究.pdf

中国科学技术大学硕士学位论文 摘要 青藏铁路穿越6 3 2 公里的冻土区，如何保护冻土层不被破坏是青藏铁路安全 运营的关键。本文针对青藏铁路工程中碎石护坡，碎石铺垫两种情况，依据路基 的不同区域分别建立不同的数学模型，同时根据不同的路基截面进行多区域的网 格划分，采用有限差分对各区域中的方程进行离散化。由于其巨大的计算量，使 用并行计算方法计算不同碎石粒径、不同路基形式下的路基及地下的温度分布， 使用最小二乘法绘制多区域复杂边界下的等值线图。主要成果如下： 1 在块石及碎石区使用多孔介质渗流理论描述，考虑热交换的气体流动。并按其 雷诺数r e 的计算值，提出使用非达西渗流来校正运动方程，更加准确的描述了 其物理机制。有关冻土中多区域流动国内外未见有关非达西流的模型描述。 2 为使方程简洁并避免不同单位制带来麻烦，给出多区域无量纲化后的偏微分方 程，并给出了无量纲化后的偏微分方程二维问题的数值解法。因计算工作量打， 实现了并行化。国内外这种偏微分方程的数值解法还未见报道。 3 从网格划分、偏微分方程有限差分、并行算法设计及实现直至科学计算可视化 ( 如曲面拟合、温度等值线及速度矢量图等) 都全部采用v c + + 语言自己编写。 而国内外相关的计算中，一般采用如f l u e n t 等著名的计算力学商业软件计算， 而后处理软件采用t e c p l o t 或o r i g i n 。 本文的计算结果对青藏铁路设计及后期养管提供了一定科学依据；解决了多 孔介质中大孔道大雷诺数下空气的湍流及传热数值计算问题。数值计算结果证实 了青藏铁路提出的“气冷路基”理论，其计算结果也与部分实验相吻合。 关键词： 非达西热流耦合；偏微分方程组：c f d ；并行算法 i i 中国科学技术大学硕士学位论文 a b s t r a c t q i n g h a i t i b e tr a i l w a yg ot h r o u g ht h ep e r m a f r o s t 6 3 2k i l o m e t r e s t h er a i l w a y e m b a n k m e n th a st w os h a p e s ：g r a v e l s l o p ep r o t e c t i o na n dg r a v e lb e d d i n g p r e v e n t i n g p e r m a f r o s tf r o mt h a w i n g i st h ek e yf o r t h es a f e l yo p e r a t i n go f q i n g h a i t i b e tr a i l w a y t h i sp a p e rf o u n d e dd i f f i r e n tm a t h e m a t i cm o d e la n dd i v i d e dm e s h e sb a s e do nv a r i a n t s e c t i o na n de m b a n k m e n t i td i s p e r s e st h ee q u a t i o n sb yf i n i t ed i f f e r e n c e b e c a u s et h e c o m p u t a t i o ni sh u g e ，w en u m e r a t e dt h et e m p e r a t u r ef i e l do fv a r i a n tg r a i ns i z ea n d e m b a n k m e n tb yp a r a l l e l a l g o r i t h m ，a n d a l s o p l o t t e d i s o l i n ei n m u l t i p l e a r e a b y l e a s t s q u a r ep r o c e d u r e m a i nc o n t e n t ： 1t h i s p a p e r d e s c r i b e st h ea i rh e a t f l o wi n g r a v e l l y e m b a n k m e n ta s n o n d a r c y c o n v e c t i o ni np o r o u sm e d i a b a s e do nr e y n o l d sn u m b e rr e ，i tp r e s e n t sn o n d a r c y i n f i l t r a t i o nf l o wt oc o r r e c tm o t i o ne q u a t i o n s s oi td e s c r i b e st h ep h y s i c a lm e c h a n i s m m o r ew e l la n dt r u l y t h en o n d a r c yi n f i l t r a t i o nf l o wm o d e lf o rm u l t i p l ep e r m a f r o s t a r e a sh a v e n tb e e np r e s e n t e db yd o m e s t i ca n do v e r s e a s j o u r n a l s 2t h i sp a p e rw r i t e dn o n d i m e n s i o n a le q u a t i o n sf o rs i m p l i f y i n ge q u a t i o n sa n da v o i d i n g t r o u l e st a k e nb ym u l t i p l eu n i ts y s t e m s i ta l s op r e s e n t sn u m e r i c a ls o l u t i o n so ft h e s e n o n d i m e n s i o n a l e q u a t i o n s a n d p u tp a r a l l e la l g o r i t h m i n t oe f f e c tb e c a u s et h e c o m p u t a t i o ns c a l e i s h u g e t h e s en u m e r i c a l s o l u t i o n sh a v e n tb e e np r e s e n t e db y d o m e s t i ca n do v e r s e a sj o u r n a l s ， 3m e s hd i v i s i o n f i n i t ed i f f e r e n c eo fp d e ，i m p l e m e n t a t i o no fp a r a l l e la l g o r i t h ma n d v i s u a lc o m p u t a t i o n a lr e s u l t s ( s u c ha sf i t t i n gs u r f a c e ，i s o l i n eo ft e m p e r a t u r ef i e l da n d v e l o c i t yv e c t o rg r a p h ) a r e a l li m p l e m e n t e db yv c + + i ti su s e db yf l u e n tt oc o m p u t e a i r sh e a t f l o wa n dt e c p l o to ro r i g i nt od i s p l a yc o m p u t a t i o n a lr e s u l t si nd o m e s t i ca n d o v e r s e a s p a p c r s t h i sp a p e rs u p p l i e ds o m es c i e n t i f i cb a s e sf o rt h ed e s i g na n dm a i n t e n a n c eo f q i n g h a i - t i b e tr a i l w a y a n ds o l v e dt h e l a r g e r en o n - d a r f yh e a t - f l o w sn u m e r i c a l c o m p u t a t i o n i n b i g d r i l l p o r o r s m e d i u m t h en u m e r i e a lr e s u l t sa n d a n a l y s i s a c k n o w l e d g e dt h e “c o o l i n ge m b a n k m e m ”a n d m a t c h e ds o m e e x p e r i m e n t a ld a t a w e l l k e y w o r d s ：n o n d a r c yc o u p l e dh e a t f l o w ；p d e ；c f d ；p a r a l l e la l g o r i t h m 1 1 1 中国科学技术大学硕士学位论文 1 绪论 1 1 研究意义 在描述和量化某一自然现象的物理机制时，先要建立正确的物理和数学模型， 找到合适的边界和初始条件，然后求解该数学模型。遗憾的是，除了少数理想状态 外，自然界中大多数问题在用数学描述的时候都是非线性的，这些问题大都没有 解析解，只能求数值解。用人工的方法求数值解是不可想象的，这就需要人们利用 计算机的运算速度快来找到数值解。典型的是计算流体力学这门学科，通过差分 或其他方法把非线性偏微分方程转换为只包含加减乘除等运算的代数方程组，编 成程序后让计算机去按规定顺序重复的不厌其烦的完成这些代数运算，最后求得 数值解。即使这样，对于求解大型的非线性问题时，一台计算机的计算速度仍然显 得让人无法忍受，这就需要用多c p u 的并行计算机求解，并研究合适的并行算法。 现在，大规模科学计算已经成为人们认识和解释自然的强大工具。 正在兴建的国家特大型项目青藏铁路工程是“十五”期间西部大开发四大标 志性工程之一，具有重大的政治、经济、军事和社会意义。青藏铁路全长1 1 4 2 k m ， 穿越多年冻土5 3 2 k m 。该铁路普遍采用了碎石路基。本文以青藏铁路碎石路基中 的非达西自然对流为具体计算对象，研究了大雷诺数非线性渗流方程的并行算法 大大减少了计算时间，使定量了解碎石路基内的大雷诺数复杂热流耦合流动成为 可能，并给出了5 0 1 0 0 年间青藏铁路冻土区温度变化的计算结果，并将计算结 果可视化，为青藏铁路设计及后期养管提供了一定科学依据。 物理背景： b alc he l e c t 效应b a l c h 川于1 9 0 0 年提出，碎石堆可降低地温夏季时碎石堆 中上部气温高，下部气温低，无对流。由于空气导热率很小( k a y 0 0 2 w m k ) ， 上部向下部主要是热传导，。传输的热量很小。冬季上部低，下部高，产生对流， 使地面降温。有人在昆仑山垭口做实验，最大降温达7 1 。利用这个原理，可 采用碎石路基来保护青藏铁路的安全。 下面是自然对流的原理解释眨1 ： n s 方程的一维形式为 p 警+ p 一警一署+ 等一p g ， 石+ p7 ：瓦一2 一言 + 萨一p 。 当流动处于稳态( 定常) 时，左边第一项可略去，右边第二项是阻力项，用d 中国科学技术大学硕士学位论文 表不z 轴沿锚重万向同上为正。方程简化为 p y 婴：一罢一d p 。g k 【l p 仃一瓦) 】 ( 1 2 ) o zd z 。一 其中一詈是正值，p 是空气的热膨胀系数，岛，磊是空气密度和温度的某个参 考值，例如取地面t o 时刻的值。图中空气受四个力作力( 如图) 当 一窖+ p o g p ( r 一) d + p o g 有v o 且_ a v 0 梯度力一呈 0 它 译剪p ，g b 证一 叶、z i 阻力d 、r i 重力岛量 图1 1 烟囱中流动示意图 世界上在多年冻土区修筑铁路已有百年以上的历史。俄罗斯早在1 8 9 5 年就 开始修建第一天西伯利亚大铁路( 后贝加尔铁路) ，该线全长9 4 4 6 k m ，穿越多年 冻土2 2 0 0 k m 。2 0 世纪7 0 年代末，又开始修建第二条西伯利亚大铁路，该线全 长3 5 0 0 k r n ，穿越多年冻土2 5 0 0 k m 。美国于1 9 0 4 年修建阿拉斯加铁路，线路长 7 5 6 k m ，通过多年冻土3 7 8 k m 。加拿大在多年冻土区建有5 条铁路。最早的哈德 逊湾铁路建于1 9 1 0 年，全长8 2 0 k m ，穿越多年冻土6 1 l k m 。用碎石作为寒区铁 路的气冷路基，已经有很多年的历史。我国目前多年冻土区已兴建的道路工程有 大、小兴安岭地区的两条，穿越冻土8 0 0 k i n 左右。另外西北有两条铁路：其一 青海海西热水专用线：其二是穿越天山的南疆铁路。高原冻土区修筑的青藏公路 ( 跨越多年冻土5 5 0 k i n ) 和青康公路( 跨越多年冻土3 3 0 k i n ) 已有4 0 多年连续 通车运营的历史。 多年冻土由于反复的冻融作用，产生许多特殊的自然地质现象，如冻胀、融 沉、冻拔，冻裂，冰锥、冻融分选、冻融泥流等，对工程建筑有极大的影响。多 年冻土区常见的道路工程病害是融沉和冻胀，以青藏公路为例，8 5 的路基病害 是融沉造成的，1 5 为冻胀和翻浆所致。在寒区工程中，根据b a l c h 效应，碎石 中国科学技术大学硕士学位论文 路基被认为是气冷路基，具有保护路基下面冻土的作用。 中国科学院寒区与旱区工程研究所的程国栋院士领导着青藏铁路工程“高 原、环保、冻土”问题的科学论证组。本文物理背景方面的论述与分析引用了程 院士及其论证组成员已经发表的文献。本文并行计算方面的算法，引用了中国科 学技术大学陈国良院士编写的著作。中国科学技术大学孔祥言教授在课题报告 关于青藏铁路堤冷却的分析和计算方法( 修改中2 0 0 3 1 2 ) “1 中建立了数学物 理模型。本文的完成，得益于以上几位德高望重的师长所做的大量研究工作和成 果。 1 2 国内外研究趋势与状况 关于气冷路基的研究已经有几十年。中国科学院寒区与早区工程研究所的研 究人员做了大量的调研、论证、实验和计算工作。程国栋院士在中国铁道科学 杂志发表的用冷却路基的方法修建青藏铁路b 1 和在第四纪研究杂志发表 的青藏高原多年冻土区路基工程地质研究h 1 中详细说明了碎石路基的保护冻 土功能和其他一些相关的保护冻土措旌，如通风管、遮阳板、泥炭等。文中也详 细介绍了国内外研究的工作进展。 国外在这方面开展的研究比较早，俄罗斯、美国、国内的一些西北地区的研 究院所都曾做过物理实验。苏联铁路运输研究院沃罗丁冻土试验研究室和美国阿 拉斯加大学都曾做过实验和观测，原中科院冰川冻土所修筑过试验路堤。这些都 证实了碎石块石路堤是“气冷路堤”。应用计算机进行计算分析的报道近几年才 有。国内外一些学者在这方面做了大量的工作。阿拉斯加大学的g o e r i n gd j 发表数篇文章展示了用达西流来描述碎石内求流的计算结果，寒早所的赖远明、 俞文兵、张明义等也用f l u e n t 软件做了各种形状、各种颗粒粒径的路基计算分 析。 目前的计算工作概括起来基本有两种： 1 、用著名的计算流体商业软件( 如f l u e n t ，c f x 等) 计算 优点：可以方便快捷的计算，这种软件由众多力学和软件专家几十年研究和实践 经验锻造，功能强大，结果稳定，一般比较可靠。尤其是它们具有强大精美的可 视化功能，省去了自己编程序计算出结果后还要再找可视化软件( 如 t e c p l o t ，o r i g i n 等) 后处理的工序，所以大部分计算研究采用了商业软件。 缺点：对于商业计算流体软件来说，可靠性要求高于精确性和理论学术性要求。 中国科学技术大学硕士学位论文 为了使结果更加“合理”，一般会做增大粘性系数等一些处理，这样使流动过程 的一些细节难以表述。不过这个对于功能如此丰富强大的软件来说，是可以理解 的，毕竟计算流体力学是一门涉及数学非常深的复杂学科，商业软件也是建立在 理论研究的基础上的。对于用户的一些复杂要求，如特殊罕见的边界界面条件， 应用这些软件计算也是非常困难的。因为看面面俱到而又繁杂的帮助文档是一件 痛苦的事情，而且这些帮助文档也不一定能解决问题。因为这些软件的源代码是 不能修改也看不到的，所以只能按它要求的步骤和限制来工作。费用问题也是 一方面，这种商业软件价格昂贵，如果用它进行商业计算的话，价格更加昂贵。 2 、采用达西流模型，自己编制串行程序对于复杂边界的气体流动，由层流到湍 流的转捩雷诺数要比规则的管道流动中小的多。在青藏高原，常年有大风，碎石 堆里的气流边界十分复杂，属于多孔介质。根据后面的物理模型分析，碎石中会 有惯性一湍流效应，所以用达西流是不太精确的。而且对于求解复杂的非线性问 题( 如湍流) 串行的程序太费时，必须要编制并行算法，在并行机上计算。 1 3 论文研究的主要内容以及创新点 第一章是绪论，介绍了研究意义。第二章给出了偏微分方程的并行数值计算 的一些方法。第三章是核心部分，建立了数学物理模型，提出了大雷诺数渗流方 程组提出了正确的差分方法和数值算法，给出了计算结果和分析。第四章给出了 串行并行程序结构和加速比。第五章是计算结果可视化。第六章是结论与展望， 总结了已做工作和不足之处，还有以后需要做的工作。最后为参考文献等。 主要创新处是在对渗流雷诺数分析的基础上，以为认为非达蘧流的数学模型 更符合路堤碎石中空气渗流的自然对流物理机制( 此前国外均用达西流模型，计 算比较简单，但不尽合理) ，给出相应的非达西渗流微分方程、能量方程和复杂 边界条件：编制了并行计算程序，并进行了大量的计算分析，与现有的部分试验 数据进行了比较，结果比较符合。 4 中国科学技术大学硕士学位论文 2 并行背景一青藏铁路碎石路基数学物理模型 2 1 计算区域示意图 本文着重讨论两类路堤模型： 图2 一i 是碎石填铺的路堤。图2 2 是碎石护坡的路堤。其中区为多年冻土 区，区为亚粘土区，为碎石区，x 轴沿原始地面。空气的对流主要发生在 区，但计算温度分布要把各个区域联合起来，并且特别看重计算路堤范围原始地 面以下几米以内活动带的温度变化情况。所涉及的主要参数列于表2 1 ，l ， 鬲武n 。：一 袋聪鞭醐麓蔽繇瞄菌懑醪彀酣一1 。 醇土后 。 碎石护坡的路基示意图 图2 1 碎石铺垫的路基示意图 l a bc n 廑蕊 。 冻土 。e f o gh 图2 2 碎石护坡的路基示意图 中国科学技术大学硕士学位论文 表2 1 所涉及的材料的物理参数表 密度比热热导率孔隙度渗透率 c p 0 k g k )( w m 中 因子 p ( k g m )k ( 脚) b ( 1 m ) 空气0 6 4 l 1 0 0 40 0 2 抛石( 多孔 2 4 0 04 2 30 3 8 70 4 4 介质) l ，5 8 l o 曲1 0 5 可e 粘土( 多 1 9 0 01 2 4 01 3 5 1 o 4 3 0 1 0 ” 孔介质) 冻土( 固体) 1 9 3 09 5 61 8 2 4 2 2 雷诺数与瑞利数的分析 ， 定义渗流雷诺数r e 为 r e ：p y d( 2 1 ) 泓 其中v 是渗流速度( 或称比流量) ，是多孔介质的孔隙度，d 是固体颗粒的直径， 将有关数据代入可得：对于块石区 筠鬻粼硼 我们知道，达西( d a r c y ) 定律的适用范围为r e 1 0 ，流动为湍流，就需要考虑惯性一湍 流效应，特别是湍流效应。 渗流瑞利数r a 定义为 足口一p 鹰k h a t p 膪h 3 a t 善 ( 2 2 ) 耻a 。 u ah 其中卢为空气的热膨胀系数，q 2 夕红，) ，是总的有效热扩散系数a h 是路堤高 度，k 是多孔介质的渗透率，r 是上下温差。分析表明，对于等温和封闭边界， 达西流的临界瑞利数r a c = 4 厅2 z 3 9 4 7 8 ；非达西流稍大于4 0 。对于等温和定压边 界临界瑞利数约为1 2 2 0 。就是说当下边界温度高于上边界，且r e 大于这个 数值就会产生对流。 2 3 物理和数学模型的建立 2 3 1 数学控制方程 对于纯空气，动量方程要用纳维一司托克斯( n a v i e r - - s t o k e s ) 方程，其控制方 程为 v v = 0( 2 3 ) p ：兰+ ( v v ) v = 一v p + 胛2 v + , o o i 一( ，一撞 ( 2 4 ) 中国科学技术大学硕士学位论文 鲁+ ( v v f = 吒v 2 r 其中v 是流体流动的质点速度， 是空气的热扩散系数。 块石区( 区) 的数学模型 对于块石区，其控制方程为 v v = 0 即= 一簧v p s i v l v 协 1 一p ( r r ) o k 盯娶+ ( v v 弦：口，v ：r o t 其中b 是贝塔因子，盯和口，分别为比热比和总的热扩散系数 盯= 黯= 盟啬盟 ( 2 5 ) ( 2 6 ) ( 2 7 ) ( 2 8 ) 铲南= 铲 在动量方程( 2 4 ) 和( 2 7 ) 进行比较可以看到，在动量方程( 2 ，4 ) 中时间项 ，8 万v 被去掉，对流项( 硝，v ) v 被换成p 口i v i v 这种经验型的项。而在能量方程 ( 2 8 ) 和( 2 5 ) 的比较中，( 2 8 ) 中加入了比热比仃和总的熟扩散系数t 2 ，因为 碎石区有两种物质，碎石和空气。 亚粘土路堤和冻土区的数学模型 对于亚粘土路堤部分( 区) ，原则上属于达西流动。但在该特定条件下，其 瑞利数远小于临界瑞利数r a c 就是说在区多孔介质中的空气基本上处于静止 状态，是单纯的导热。区和区所适用方程形式相同，都是热传导方程 j 一0 ) 等+ bp ) f 譬= 女。v ：z ( i 钆2 ) ( 29 ) 其中l i 是区的冻融交变体积潜热，t i 是区的平衡温度。t f 是含盐冰的融化 温度，如果不考虑盐分的影响，就是0 c 。 2 3 2 无量纲化 对于开放边界 可在方程( 2 t 6 ) ( 2 8 ) 中消去变量v ， 量形式写成 盟+ 竺：o 融出 罢一等_ 一芦mm 。 化成只含p 和t 的两个方程。我们用分 ( 2 1 0 ) ( 2 1 1 ) 中国科学技术大学硕士学位论文 考= 一尝吃一膨吲眨一成g k n 一仃一t ) 】 盯詈+ 屹罢+ t 罢= q ( 警+ 窘) 引进所谓“非达西渗透率”k n 它是时间t 和空间变量的函数 k = kk 而2 忑磊 ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) 其中c f = b 4 k 是无量纲型阻系数 则由式( 2 1 1 ) 和( 2 1 2 ) 有 矿：一生望( 2 1 5 ) “a x _ ：一生l 罢+ p o g ( 1 一雕一t ) ) i ( 2 1 6 ) 将上二式分别对x 和z 求偏导数，再相加并利用式( 2 1 0 ) 可得 鲁+ 鲁+ 土f 警罢+ 监罢 鸭g pa ，t + 叙2瑟2 k ，氟良a z 出j 出 去, o g e l 一( r c ) 警+ g ( 丁一) 掣= 。( 2 1 7 ) 而( 2 1 3 ) 可写成 a 百o t 一百k u 磊劫面o t 一警 警+ 成g d 一卢仃一瓦) 】) 警= q ( 窘+ 窘 ( z s ) ( 2 1 7 ) 和( 2 1 8 ) 就是用于求解p 和t 的有量纲方程组。 在研制和开发软件时，为使运算通用化和规范化，通常是将方程( 和边界条 件) 无量纲化。引进无量纲量 x 2 万 z 2 百如= 暑， 世：监， “ k ：里屹，p 。：盟型，p q a o r t ：堕， 口 ：旦 口 等，= 警 亿1 9 ) 温度标尺t 应是有参考价值的恒定值。我们用地面平均温度一1 0 减去顶面 最低温度一1 3 在本向题中即取r = 1 2 k 于是式( 2 1 7 ) 和( 2 1 8 ) 化成无量纲方程组为 磐+ 磐+ 上f 堕监+ 坠丝1 _ r a 旦+ 麟磷k dl 橱d 硝da z do z dja z d 8 中国科学技术大学硕士学位论文 去( 肌脚) 等+ ( 1 一肚朋) 等= 。 鬻+ 骞+ k 。l ( o p d a o o p o 0 0 j + 1 岛( 科一r 甜) 嚣= 瓦0 0 而速度的无量纲形式由( 2 1 5 ) 和( 2 1 6 ) 化为 = 一磐 似d 吲= 晚+ 吃】 2 一玺一僻一r 日) ( 2 2 0 ) ( 2 2 1 ) = s i n a + c o s a ( 2 2 2 ) 其中a 是b d 与x 轴所夹锐角，而k d 、贝塔数b e 和瑞利数分别为 = 譬2 而习硼12 雨瓦眄1 丽( 2 2 3 ) b p ：a t , b k 啊 a t c f 4 k 啊 r d 。：p o g k h ，砌：p 。g k h f l a t ( 2 2 4 ) j a ，j c x r 对于区和区，引进无量纲量1 区【俨，) ，当作( p 勺) 1 ，所以o2 1 。 2 瓣l ，9 2 簪 ( 22 5 ) 则热体导方程的无量纲形式为 扛。子p 一巳) + - 尝= 骞+ 骞 ( z z e ) 边界条件 在a b 上：芦b 。：1 1 ：p f z ：o ) 一r 口o ( 1 一j 口a 丁口)( 2 2 7 ) - 1 ) = 筹+ 吕s i n ( 彘+ ) ( 2 2 8 ) 余类推。 方程( 2 2 0 ) 和( 2 2 1 ) 的数值解 方程( 2 2 0 ) 和( 2 2 1 ) 中，由于k d 是( x 。z 。t 。) 函数，所以方程是高度非线性 的，通常的算法是 作为零级近似，取k d 1 则方程退化为达西流情形，容易解出0 。口) 。 及其 导数。代入( 2 2 2 ) 和( 2 2 3 ) 求得v 。和及其导数。将它们代入( 2 2 0 ) 和( 2 2 1 ) 可求得1 级近似解0 。口) 。再迭代一次给出0 。臼) ：，就比较精确了。 对于封闭边界，由方程( 2 1 0 ) ，引进流函数比较方便，即 9 中国科学技术大学硕士学位论文 t = 詈，屹一业0 x ( 2 2 9 ) 将它们代入( 2 1 1 ) 和( 2 1 2 ) 后，用( 2 ，1 2 ) 对x 求导数减去( 2 1 i ) 对z 求导数，可 消去压力p ，而得只含流函数缈和温度函数t 的两个方程。 即 ( t + 半吲 v 2 v = - 等陪警扣a 彪i f , a 玎m 1 譬豢她。， 引进无量纲量 x d = 考，z d = 寺，t o = 嚣，妒。= 詈，k = 等，口= 等2 百2 百2 赢。2 i 2 百“2 百 则方程( 2 3 0 ) 化成无量纲形式为 ( t 他) v 2 舻也( 等掣+ 玺瓦o l v 。l 卜薏眨。t ， 其中 贝塔觏：等，州2 + ( 斟r 耻警眨。z ， 而能量方程的无量纲形式为 v 2 p 办篝 弦3 3 其中雅可比j ( o 。) 为 j ) ；瓦0 0 等瓦0 0 等( 2 3 4 ) 2 3 3 边界条件和界面条件( 以护坡示意图为例) ( 1 ) 开放边界条件 船处温度 t a b c = 2 0 + 1 5 s i n ( 彘屯+ 地面e f 处温度 t g v = - 1 0 + 1 2 s i n ( 淼” 侧面c f 处的温度可用和珞插值处理 其中t = 1 为一小时，8 7 6 0 为每年的小时数。 ( 2 3 5 ) ( 2 3 6 ) l o 中国科学技术大学硕士学位论文 冻土层内部g h 处的边界条件按q f 0 0 6 w m 2 和k l = 1 4 7 4 w m k 可得在g h 处 娶：争_ 0 0 4 0 7 k m 矾0 4 k 肌 ( 2 3 7 ) 盘丘 冻士中侧面f h 处因离受扰区较远，认为绝热，即热通量罢l 。= o ， a b 、 b c 、c e 、b d 、o d 、d e 按内部的界面条件处理。 在开放边界处的压力边界条件为 字= 一善e ( 2 3 8 ) 瓠k v 。 ” 罢= 一告圪一p o ( z ) g 1 一卢( r 一瓦) 】 ( 2 3 9 ) t y z a ( 2 ) 封闭边界 封闭边界的温度条件如开放边界一样，在具体计算时，把压力p 消掉，换 成流函数缈。b d 、d e 、c e 、b c 四边上v = 0 。 ( 3 ) 界面条件 b d 上的力学条件为压力连续和流量法向通量连续，热界面条件为假设所有 界面处接触热阻( 或膜系数) 为零，是良好的热接触，则有温度连续和热通量 连续。o d 和d e 上也是温度连续和热通量连续。 七i 罢“，罢 ( i ，j = l ，2 ，3 ，且i j )( 2 4 0 ) a 聆，a 珂， 2 4 控制方程的差分格式 在计算流体力学中，给出控制方程后，需要根据方程的形式和一些关键参 数的大小来选取合适的差分格式来保证格式的相容性、稳定性和收敛性。为方便 起见，在本文所做的网格划分一律取血= 母= a ，因为在本论文研究的实际物理 背景中，x 和y 方向的流动性质基本没有大的差异。二阶导数项一律采取中心差 分，如果气流流速不大的话一阶导数可为中心差分也可为迎风格式，如果流速较 大，最好采用迎风格式。网格划分如下： 中国科学技术大学硕士学位论文 图2 3 结构化的计算网格示意图 2 4 1 冻土区和亚粘土区控制方程的差分格式 冻土区和亚粘土区只有一个热扩散方程，熟扩散方程的数值解法相对简单， 本文采用最常用且稳定高效的a d i ( a l t e r n a t i n gd i r e c t i o ni m p l i c i tm e t h o d ) 方法求解，即把一个时间步长分作两个半步推进，第一个半步的推进是x 方向取 隐式，第二个半步的推进是y 方向采用隐式。离散化后如下 “= l l f ，0 + i 万 t i l t t + l ，2 2 虼”2 + 略? 2 + 嵋川一2 吧+ 虼一，) ( 2 4 1 ) 虼1 = y 孑“2 + ：轰( 睢盖。- 2 g 1 + 。n + 一l 。+ 妒= ：j 2 2 ：芦”+ n + 。l 2 ) ( 2 4 2 ) 该格式运算速度快，每个时间步解两个三对角矩阵，如果三对角矩阵的规模是 r ，那么算法的时间复杂度是0 ( ) 。对此差分格式用v o n n e u m a n n 方法进 行稳定性分析，可知本格式是无条件稳定的。 2 4 2 砰自区域径制万程 开放边界： 动量方程 磐+ 缉+ 一1f 坠亟+ 逊赴b 盟+ 瓠乞a z 乞k d 、a x da x d8 z da z o )o z d 去( 脒也口 ，i o 瓦k o + ( 一t o ) o p 。= 。 能量方程 鬻+ 鬻+ ( 薏嚣+ a p , , a o i j + 蹦砌吒脚) 薏= 薏 ( 2 4 3 ) ( 2 4 4 ) 动量方程中压力p 是未知量，由扩散项，对流项和“源项”组成，而没有时 间项。把相应的参数代入后发现，对流项的系数比扩散项的系数小1 到2 个数量 级左右，所以该方程近似于一个“变源”的泊松方程。分析“源项”后可知，源 的作用主要来自空气的热胀冷缩和自然对流。 中国科学技术大学硕士学位论文 能量方程中温度口是未知量，对流项的系数在1 0 4 左右，说明能量的传递主要 依靠气流来输运，而不是导热。对流项系数是变化的，说明气体的流速及流向随 时间发生变化。虽然导热项很小，但是不能把它忽略。忽略后的方程为欧拉能量 方程，这个方程并不适用高速渗流，因为在大雷诺数的气体流动中，流动是等熵 的，全场的温度基本上为一个恒定值，这不符合青藏铁路路基的实际情况。解决 这个问题的正是热扩散项，热扩散项系数里由孔隙度、总的比热比和热扩散系数 构成，说明了气体与多空介质( 石头) 、多孔介质中间的热传导。上述微分方程 非常恰当的描述了多空介质( 碎石) 中高速渗流的气体流动和能量传递物理机制。 困难在于这两个方程的数值解法，尤其是能量方程。如果是欧拉方程的话， 二维问题可以变成两个一维问题的叠加。热扩散方程的话可以采用a d i 方法。为 了稳定性和收敛性，我们采用了中心差分加全隐式格式。缺点是计算速度非常慢， 这一点也反应了物理实质：因为空气流动引起的热量传输速度远远高于热传导， 在本例中，气体对流引起的热量变化是以0 0 1 0 0 0 1 秒来计的，而碎石、冻土、 亚粘土的热传导是以小时来计的。在热量传输速度角度来讲，这是一个时间多尺 度问题。这个困难由并行和系数矩阵的应用而解决。并行的内容在其他章节介绍。 差分形式： 动量方程 p i ，- l 一2 p j + p 。j “p 一，j 一2 p f j + p i q j 5 x 。止 l f k h l k j ，lp w + l 一只，1 茂 ，一k “j 只t j 一屁川1 ( 2 4 5 ) k l 2 5 x2 a x 。 。 2 血2 a z j 也訾亡( r a 。- r a o , ，) 气+ ( t 一耻r ) 气= 。 能量方程 吲。一2 野1 + 吲，一2 嘭1 + 甜j 。 。夏2 1 。出2 1 上f 丛二垃! 型1 二型、纽二丛墅：二墅：1( 2 4 6 ) k ，i 2 a x2 a x2 a z2 a z j 蝇，( r a 。- r a o j f 1 ) 警= 华n + l n 其中每一时间步的k + 由上一步的值代入公式得到。 解上边两个方程，得到如下形式 动量方程组 4 ，。b ，= q 能量方程组坟。b ，= 6 j ， ( 刀为内网格点的数目) 4 。和最。都有如下的形式： 图中x 为各个点的系数，大部分非零。矩阵其他地方均为零，是个有着一定排列 中国科学技术大学硕士学位论文 规律的稀疏矩阵。 图2 4 动量和能量方程矩阵图 封闭边界： 动量方程 c 拙妒也( 玺掣+ 玺掣 一心嚣 t t ， 叫( 警) 2 + ( 斟r 亿t e ， 能量方程 v 2 0 _ 盟监+ 盟监：丝 瓠da 2 d融d 瓠do t o 差分腻”叫v 2 眈( 玺掣+ 玺掣) 一心薏 t 。， 动量方程 迭代在工程数值方法中经常用到，迭代的方法非常快，但是必须保证是对角 占优的情况下。为了简单的说明问题，我们先假设上面问题是一维问题。( 1 ) 上 面都是中心差分，可以看到 的系数是- 2 ( i + b e ) ，一。和y 。的系数为 l + ( 1 一型型) 和1 + b e ( 1 一掣型) ，肌和要型的值远远小于l ，这样的系数是可 以保证对角占优的。可以采用迭代的方法。 可以看到蚓噼( 玺 2 + ( 玺门m ，这个方程组是非线性的，不能直接 用代数方法求解。我们用迭代的方法求解，用右上标。似表示迭代i 次后的流函 数值。设当前时间步为第n 时间步，那么把上一时间步的”1 当做。o 并求出 1 4 中国科学技术大学硕士学位论文 ”，代入上式有 一见竺：f ( 4 v j o , v o ( ) ) d x d ( 2 5 0 ) 基于此式可以得到n 个以髫为未知数的线性代数方程组，解之得到。出。用 迭代法，并且左端保持不变： 一r 。0 0 ：，( ( 2 1 ，v d ( 1 ) ( y 。( 1 ) ) u x d 亿娑：，( p 。，垆( 。) ) u x d 一兄娑：，( ，矿- ( y 。) ) ( 2 5 1 ) g a d 一般的，l ；c ，。“会很快地向精确值”收敛。计算表明，迭代四次的情况下，误差 会一直保持在2 之内。 能量方程：能量方程，同样的对角占优分析。在f 足够小时，a d i 是个非常有效 的方法，显式也可以收敛。但是对于5 0 年的问题来说，秒或更小的时间单位的 显式或a b i 并不比以小时为时间单位的全隐式划算。能量方程可以显式，也可以 隐式。 2 5 计算结果和分析 2 5 1 几种情况的计算结果 l 开放边界条件下的粘土在上碎石在下 1 米厚的情况 1 碎石下0 高度线位置处一年内的温度变化是：秋天最高，在2 。c 3 之间。 春天最低，在一1 o 之间。 冬天基本都是o 左右。夏天在0 1 2 最高温度和最低温度不是出现在夏天和冬天，这是因为有一个滞后效应。就 如中午太阳直射对地表温度不是最高一样。 3 三十年的变化可看出此模型有一个降温效果，但不显著，在二十五年间，0 高度线处的温度大约降了0 2 o 4 c 碎石厚度2 米和1 米做了比较，可以看到，两米的降温效果比一米的好，结果显 著，同一位置的温度有o _ 3 0 5 的降温。这个明显是b a l c h 效应的结果。 中国科学技术大学硕士学位论文 图2 5 i 米 厚度 图2 6 1 6 中国科学技术大学硕士学位论文 图2 7 。，l f 。、。，。，。t i 抖 “u t y n v ，。，。，w “t 汰 1 米 厚度 夏天 速 度 矢 量 图 、w “l “t ，打v 打* v “伸_ _ - _ _ _ _ _ _ h _ n v ，f 怵苌k 、- - “i “l “ h 铀竹帕w _ 竹h _ _ _ _ 帕w “饵- k w r t i & i tl j | i i i i i i i a i | l i l a l l i “l 竹w h 竹_ - - 啊w ，劓n 叫_ h 图2 8 二二 中国科学技术大学硕士学位论文 图2 9 赫秋天 献藿 川。滁羹 。，。二，川。闷目h 图 铀- 一“y “镰 ，w ，_ 岬_ - _ _ 啊h 刊h 卅_ - ，r 坩v i 、叫a k w 1 “l “l “l 帕竹p w 忡 _ _ w _ _ _ _ h _ _ _ _ - _ - _ - _ h - v y j l q _ 图2 1 0 二二 中国科学技术大学硕士学位论文 图2 1 1 f 。l 。：，。，。h 。h ，。滢| ，筹薅 、。、。