冻结土壤温度场数值模拟的改进.pdf

第3 4 卷第2 期中国矿业大学学报v 0 1 ．3 4N o ．2 2 0 0 5 年3 月J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g T e c h n o l o g yM a r ．2 0 0 5 文章编号1 0 0 0 1 9 6 4 2 0 0 5 0 2 0 1 7 9 ．0 5 冻结土壤温度场数值模拟的改进 商翔宇1 ，周国庆h 2 ，别小勇1 ’2 1 ．中国矿业大学建筑工程学院，江苏徐州2 2 1 0 0 8 ； 2 ．中国科学院寒区旱区环境与工程研究所冻土工程国家重点实验室，甘肃兰州 7 3 0 0 0 0 摘要针对土壤冻结过程的数值分析，基于显热容法用有限差分法构造出热传导方程变空间步 长的半隐格式和全隐格式．结合实例分析了F E M 和F D M 方法造成的‘相变遗漏’，提出一种自 中图分类号T U2文献标识码A 模型求解迭代大为简化。有效避免了‘相变 解冻土温度场，是对传统方法的有效改进， N u m e r i c a lS i m u l a t i o nI m p r o v e m e n to fF r e e z i n g S o i l ’ST e m p e r a t u r eF i e l d S H A N GX i a n g y u l 。Z H O UG u o q i n 9 1 ”，B I EX i a o y o n 9 1 ’2 1 ．S c h o o lo fA r c h i t e c t u r ea n dC i v i lE n g i n e e r i n g ，C h i n aU n i v e r s i t yo fM i n i n g8 LT e c h n o l o g y ， X u z h o u ，J i a n g s u2 2 1 0 0 8 ，C h i n a ； 2 ．S t a t eK e yL a b o r a t o r yo fF r o z e nS o i lE n g i n e e r i n g ，C A R W W E I ，C A S 。L a n z h o u ，G a n s u7 3 0 0 0 0 ，C h i n a A b s t r a c t B a s e do nt h en u m e r i c a lc a l c u l u so fs o i lf r e e z i n ga n dt h ea p p a r e n th e a tc a p a c i t ym e t h o d ， c o m p l e t e l ya n dp a r t l y ‘i m p l i c i td i f f e r e n c es c h e m e so ft h eh e a tt r a n s f e re q u a t i o nw e r ed e d u c e db yt h e f i n i t ed i f f e r e n c em e t h o dw i t hv a r i a b l es p a c es t e p ．T h e nt h ep h e n o m e n o no f ‘p h a s ec h a n g el o s t ’， r e s u l t e df r o mt h ef i x e dt i m es t e p ，w a sa n a l y z e dt h r o u g hat y p i c a ln u m e r i c a le x a m p l e ．A n dt h e na m e t h o dw i t hs e l f a d j u s t e dt i m es t e pi nw h i c hc o m p l e t e l ya n dp a r t l yi m p l i c i td i f f e r e n c es c h e m e sw e r e u s e da l t e r n a t e l yw a sp u tf o r w a r d ．T h i sm e t h o dc a na v o i dt h e ‘p h a s ec h a n g el o s t ’s u c c e s s f u l l y ，a n d a sar e s u l t ，t h ei t e r a t i o np r o c e d u r ec a nb es i m p l i f i e dg r e a t l y ．A tt h es a m et i m e ，i tc a nr e d u c et h e c a l c u l a t i o nt i m ee f f i c i e n t l y ．T h em e t h o dp r o p o s e di sa ni m p o r t a n ti m p r o v e m e n tt ot h et r a d i t i o n a l n u m e r i c a ls i m u l a t i o n ，w h i c hc a nh e l pt oa n a l y z et h ep r o c e s so fs o i lf r e e z i n g ，f r o s th e a v ea n dt h a w s u b s i d e n c ea c c u r a t e l y ． K e yw o r d s s o i lf r e e z i n g ；n u m e r i c a lc a l c u l u s ；p h a s ec h a n g el o s t 我国东北、西北地区的寒区土木工程方兴未 艾，人工地层冻结技术在土木工程建设中也越来越 发挥重要的作用‘1 。4 ‘．要分析天然冻土和人工冻土 的冻结工程，必须准确确定土体的温度变化规律． 常见的冻胀模型如K o n r a d 的分凝冰模型‘5 3 以及 收稿日期 基金项目 作者简介 D a i y uW a n g 的冻胀模型嘲在计算冻胀量之前都要 求准确求出冻土的温度场．然而，在用有限差分法 求解温度场过程中，由于时间步长取法不当很容易 造成‘相变遗漏，．本文提出一种自调节时间步长、 半隐和全隐格式交替使用的方法，该方法能够有效 2 0 0 4 0 2 1 6 国家9 7 3 项目 2 0 0 2 C B 4 1 2 7 0 4 ；中国科学院知识创新工程重大项目 K Z C X l 一S W 一0 4 ；中国科学院冻土工程国家重点实验室 开放基金项目 S K L F S E 2 0 0 3 0 4 商翔宇 1 9 7 7 一 ，男，河南省洛阳市人，博士研究生，从事岩土工程方面的研究． 遗也 蒸一 黼撇荆龋 一一～一 隐差臻扮 全误皎借 和析蚴数 一一一一 翱堋涨潞～一一{ 蓦 万方数据 1 8 0中国矿业大学学报第3 4 卷 防止传统计算过程中出现的相变遗漏，并且计算速 度较之固定时间步长方法有了很大的提高． 1 数学模型 常见的一维土壤冻结过程热传导方程为 c y d d T fd d 烈{ 。3 a 7 z “ L 肛d 面O i ， 1 式中C y ，A 为土壤体积热容量 J ／ m 3 K 和热导 率 J ／ m S K ；L ，T 为融化潜热 J ／m 3 和温度 ‘C ；胁，岛为冰密度 k g ／m 3 和冰含量 m 3 ／m 3 ；z ， t 为空间坐标 m 和时间 s ． 求解上述方程的一种有效方法是显热容法，即 在发生相变的一个小的温度范围之内构造比热函 数，以代替以上方程的右端第二项，只以温度为待 求函数‘7 8 1 c d T d z r d td誓dj ， 2 。y z ＼“zJ ’ Ⅵ7 其中 f C f 丁 T f △丁 ， A 。 f T 丁f A T ， 式中c ；为等效热容；A ’为等效导热系数；c ，， 分 别为已冻土的体积热容和导热系数；C 。，A 。分别为 未冻土的体积热容和导热系数；乃一舒≤丁≤丁， A T 为假设小的相变温度范围． 2 差分格式与相变遗漏 方程 2 是与时间和空间都有关系的发展方 程，采用有限差分法能对其有效求解．在常用的有 限差分格式中，显格式精度最差，隐格式次之，而 C r a n k N i c o l s o n 格式精度较高踟；但是从解的稳定 性角度分析，则隐格式的稳定性最好，因此本文采 用隐格式进行求解．根据方程系数的取法不同，隐 格式又有全隐格式和半隐格式之分．全隐格式公式 如下 C 1 ； ； 1 丁； 1 7 1 ； 一j 生一 z 汁1 2 r 1 [ ∽卅n 忱l 鼍著卅Ⅵ地笨署] ， 3 ， 式中J ，船分别为空间和时间坐标标记。 所谓半隐格式是将式 3 中系数c 0 和A 的时 间步n l 换成竹，通过半隐格式只需求解线性方程 组即可．事实上，半隐格式和全隐格式对本方程的 求解精度影响不大可忽略不计[ 4 ] ．作者也分别编制 了两种不同格式的求解程序进行了验证． 由式 2 ， 3 可见，相变的影响在模拟中是通 过在假设相变范围的等效热容来实现的，因此要全 面考虑相变的影响，必须保证在每一计算时间步内 所有发生相变的空间单元的等效热容能够包含相 变潜热项．而土体空间单元的等效热容是否应该包 含相变潜热项，则是通过该单元的温度值来判断 的．由式 2 知，只有其温度落在T { 一△T ≤T ≤T f 凹范围之内，该单元的相变影响才能被考虑．这 样，就不免会出现‘相变遗漏，的情况．比如空间某 单元在某一时间步之内温度由T 1 n △丁变为 T 2 兀 △丁而乃一△丁≤T 2 ≤T f △丁 时段初 未冻而时段末处于过渡期 则相变作用就会遗漏； 采用全隐格式，若T f - A T T 1 T f 凹而T 2 兀 △丁 且丁2 丁f 一△丁则不论采用半隐或是全隐格式相 变作用都会被遗漏掉． 可见如果仅仅采用一种差分格式，那么无论是 半隐格式还是全隐格式，对于本文所述问题都会不 可避免地存在‘相变遗漏’这一事实，因此是不正确 和不可取的．由于冻结过程伴随着大量相变潜热的 释放，所以原来采用固定格式算法得出温度将会偏 低，而冻结锋面推进速度则会偏快． 3 修正方法 为了避免数值分析中相变作用的遗漏，必须有 效地控制时间步长．本文根据半隐格式和全隐格式 的不同对方程的求解精度影响不大这一特点，提出 万方数据 第2 期 商翔宇等冻结土壤温度场数值模拟的改进 1 8 1 了自调节时间步长、半隐和全隐格式交替使用的方 法 M o d i f i e dF D M 以下简称M F D M ，其中半隐格 式直接求解线性方程组即可，全隐格式则要采用迭 代法求解．由于本文提出的方法采用的是稳定的计 算格式，而且采用了自调节时间步长的技术，所以 该方法具有稳定和计算速度快的特点． 采用该方法进行数值分析的步骤如下 1 对第一时间步的计算进行特殊处理第一 时间步采用较小的步长，由于第一时间步内温度变 化比较剧烈，而且土体从初始的正温开始变化，所 以先采用半隐格式进行计算，然后令所有温度落在 丁 乃 △T 的空间点的温度为n ，再对第一时间 步进行全隐格式计算； 2 将上一时间步末的空间点按其温度值分为 两类一类是温度落在丁≤T f △丁的空间点，另一 类是温度大于丁r 凹的点；同时采用与上一时间 步相同步长的时间步，用半隐格式对本时间步进行 试算； 3 对第2 步中划分出的温度大于丁r △丁的 点的试算结果进行判断如果这些点中的温度最小 值 在自上而下冻结的情况下，此点应为这些点中 最上端的点 小于T r A T ，就减小时间步长，再用 半隐格式对此时间步进行试算；若此最小值大于 丁f △丁，则增大时间步长，再用半隐格式对该时间 步进行试算，直至该最小值落在E T t A T ，T t △7 1 ] 之中； 4 用第3 步中确定下的最终时间步长对2 中划分出的温度大于7 ’f △T 的点进行全隐格式计 算； 5 重复2 ～4 步，直至计算时间结束． 可以设想，采用了上述改进办法，相变遗漏避 免了，对于冻结问题也就是说释放出来的相变潜热 被全面地考虑了．采用固定格式算法得出的温度值 和冻结锋面都会得到相应的纠正，温度会相应升 高，锋面推进会相应减缓． 4 计算实例及分析 采用文献[ 1 0 ] 提供的基于F E M 的经典算例． 空间总长为1m ，计算总时间为2 ．8 8 1 0 5S ，采用 第一类边界条件，T c 一一2 0 ‘C ，h 一1 0 ℃；初始温 度T 。一1 0 C ．参数如下 T I 一0 ‘C ，A T - - - - 1 ‘C ， L 一3 3 8 0 0 0 0 0 0 J ／m 3 ， , i f 一2 ．2 2 J ／ m s ‘C ， 。一0 ．5 5 6 J ／ m S ℃ ， C f 一1 7 6 2 0 0 0 J ／ m 3 ‘C ， C 。 4 2 2 6 0 0 0 J ／ m 3 ℃ ． 1 采用固定时间步长半隐差分格式 F D M 求 解．将空间等步划分为4 0 0 个单元，时间等步划分 为1 9 2 个时间步长．计算结果与文献[ 6 ] 提供的结 果几乎相吻合．通过比较可以看出多处存在相变遗 漏，其中计算时间初始阶段出现较多，随着计算时 间的增加相变遗漏现象逐渐变少，见表1 和表2 ． 表1 中在1 ．5 1 0 3 ～3 1 0 3s 时间步之内，0 ．0 2 ～ 0 ．0 2 75m 处4 个空间点显然发生了相变，且它们 在时间步初的温度值都大于，也即超出了考虑相变 的温度范围，由于采用半隐格式，根据时间步初的 温度值判断的结论是这4 个空间点都没有发生相 变 表中黑体数字 ，因此它们的等效热容将不会包 含相变潜热项，也就是说，本时间步之内这4 个点 的相变都被遗漏掉了．同样在3 1 0 3 ～4 ．5 1 0 3S 时间步之内，0 ．0 31 T I 和0 ．0 3 25m 两点的相变被 遗漏掉了；在4 ．5 1 0 3 ～6 1 0 3S 时间步之内， 0 ．0 3 5m 和0 ．0 3 75m 两点的相变被遗漏掉了；在 6 1 0 3 ～7 ．5 1 0 3S 时间步之内，点0 ．0 4I n 处的相 变被遗漏掉了．在表2 中因为相变不像表1 所示的 时间初那样剧烈，因此相变遗漏相应减少，只有黑 体数字表示的两点发生相变遗漏． 2 用M F D M 求解将空间变步划分为1 1 3 个 单元，初始时间步长设为0 ．2 5s ，采用此种方法只 芸焉囊至一芸百蝴墨嘉篡篡嚣羔懈墨乒攥m凳霎誊喜再祗|||一一一一一一一一一一邺埘猢一一雌一一一一一一吣蝴㈨ 一～一一一一一～一一一一一一一一～一一一一一一一～一一『一一～一一一一}量|一一一一一一『一一一 一 一 5 5 5 5 一l l ，- R _ 一 一 5 5 1 一 n 一 2 5 7 2 5 7 一摹} o ．、一 n 一 2 5 7 k一矶一∞睨∞∞∞∞∞∞叫一∥表 №一矶一匏盟毖毖 洲一湘一毗一一一吣一一一一一标 ～一湘一毗一一一 万方数据 1 8 2 中国矿业大学学报第3 4 卷 计算了6 8 个时间步计算就结束了，计算结果见表 3 ．由于对每一时间步内黑体数据以上的点采用半 隐格式，而黑体数据及其以下的点则采用全隐格 式，所以‘相变遗漏’被有效避免了．比如在表3 中 的3 ．1 7 X1 0 3 ～3 ．6 8 1 0 3S 时间步之内，0 ．0 2 5 ～ 0 ．0 2 8m 处3 个空间点采用了半隐格式，点0 ．0 3 m 处则采用了全隐格式，它们相应的温度判据都 落在了] - T f A T ，T f △丁] ，因此实际发生的相变都 被考虑到了．在3 ．1 7 1 0 3 ～3 ．6 8 1 0 3S 时间步之 内，点0 ．0 2 5m 处已不存在相变，点0 ．0 2 6m 和 0 ．0 2 8m 处均存在相变，对这两个点必须用半隐格 表3 采用了本文方法计算后的相变情况 T a b l e3‘P h a s ec h a n g e ’p r o f i l ea tt i m es t e p sw i t hM F D M ，t／103 S z ／m～ 3 ．1 7 3 ．1 73 ．6 84 ．1 95 ．2 25 ．2 5 6 ．2 7 0 ．0 2 5 0 ．0 2 6 0 ．0 2 8 0 ．0 3 0 0 ．0 3 1 0 ．0 3 3 0 ．0 3 5 0 ．0 3 7 0 ．0 3 9 式才能考虑相变的影响，而点0 ．0 3m 和0 ．0 3 1m 处因为时间步初末的温度值都在[ 7 1 r 一△T ，丁r A T 3 ，采用半隐或全隐格式都能考虑到相变影响， 点0 ．0 3 3m 处也存在相变但只有时间步末的温度 值在相变范围内，因此只能采用全隐格式．根据设 计的算法在0 ．0 3 3m 点以上各点采用半隐格式，在 0 ．0 3 3m 点 包括该点 以下各点采用全隐格式，是 可以满足全部考虑相变影响要求的． 图1 横轴表示距土体低温端的空间距离，纵 轴表示温度值 比较了用F E M ，F D M 和本文方法 计算所得结果中的两个时刻的温度场．可以看出， F E M 和F D M 计算结果几乎相同，而用本文方法 计算的土体中同一点的温度值则比F E M 和F D M 的均要高．这符合本文前述的分析和设想传统 F E M 和传统的固定格式的F D M 均存在相变遗 漏，致使没能全面考虑释放的相变潜热，因此计算 出的温度偏低．从图1 还可看出传统算法和修正算 法之间最大的误差出现在冻结锋面附近，这也是符 合实际情况的．因为实际冻结过程中相变发生最剧 烈的地方就是冻结锋面附近，固定格式算法遗漏相 变最多的地方自然也在此处． M F D M F D M F E M b t 1 4 ．4 t 0 4S 图1各种计算方法的温度场比较 F i g ．1C o m p a r i s o no ft h et e m p e r a t u r ef i e l dr e s u l t sb yt h eF E M ，F D Ma n dM F D M 由于本文提出的方法考虑了原来固定时间步 长法没有考虑的‘相变遗漏’的影响，所以比传统方 法更加全面地考虑了水相变成冰释放出来的潜热， 这必然延缓冻结锋面的推进速度．表4 比较了采用 半隐格式固定时间步长方法和采用自调节时间步 长方法求出的冻结锋面的变化情况．从图2 横轴 表示冻结时间，纵轴表示冻结锋面距土体低温端的 表4 本文方法和传统方法冻结锋面位置比较 T a b l e4C o m p a r i s o no ft h ef r e e z i n gf r o n t ’p o s i t i o n o b t a i n e db ya u t h o r ’sm e t h o da n dt h et r a d i t i o n a lm e t h o d 空间距离 可以看出，随着计算时间的增加，两者冻 结锋面位置差距越大，即用传统固定格式算法没能 考虑的相变潜热量积累越多．这也进一步验证了本 文前述的分析和设想． 图2 本文方法与传统方法冻结锋面位置比较 F i g ．2C o m p a r i s o no ft h ef r e e z i n gf r o n t ’p o s i t i o no b t a i n e d b ya u t h o r ’Sm e t h o da n dt h et r a d i t i o n a lm e t h o d 6 7 5 7 O 6 5 l 4 9 O l 9 9 1 3 7 3 ● ● ● ● ● 4 4 3 l O O 0 O l 一 一 一 一 一 一 1 5 7 6 4 l 9 9 5 6 6 6 4 1 5 9 6 3 3 2 1 O 0 O O l 2 一 一 一 一 1 2 7 6 6 2 6 9 1 5 3 1 3 l 5 1 7 4 ● ● ● ● ● ● ● ● ● 3 2 1 O 0 0 l 1 2 一 一 一 一 7 2 9 7 8 6 2 7 1 4 8 6 O 5 9 6 2 9 5 ● ● ● ● ● ● ● 1 O O O O l Z 2 3 一 一 1 6 7 2 5 5 3 1 7 2 4 7 O 5 9 6 3 9 6 2 ● ● ● ● ● ● ● ● ● O O O O 1 2 2 3 4 53 5 9 4 5 4 1 4 5 0 5 9 7 4 1 8 4 O ● ● ● ● ● ● ■ ● ● O O O l 2 3 3 4 5 万方数据 第2 期商翔宇等冻结土壤温度场数值模拟的改进1 8 3 5 结论 本文采用的自调节时间步长、半隐全隐格式交 替求解的方法能全面考虑冻结过程中的相变作用， 使温度场分析结果更加符合实际．实际上，这种方 法修正了传统算法的计算方法，不仅对土壤冻结问 题适用，而且也能够为其它发生相变的传热问题所 用． 参考文献 [ 1 3 陈瑞杰，程国栋，李述训，等．人工地层冻结应用研究 进展和展望[ J ] ．岩土工程学报，2 0 0 0 ，2 2 1 4 0 4 4 ． C h e nRJ ，C h e n gGG ，L iSX ，e ta 1 ．D e v e l o p m e n t a n dp r o s p e c to fr e s e a r c ho na p p l i c a t i o no fa r t i f i c i a l g r o u n df r e e z i n gE J ] ．C h i n e s eJ o u r n a lo fG e o t e c h n i c a l E n g i n e e r i n g ，2 0 0 0 ，2 2 1 4 0 4 4 ． [ 2 3 程国栋．冻土力学与工程的国际研究新进展2 0 0 0 年国际地层冻结和土冻结作用会议综述[ J ] ．地球科 学进展，2 0 0 1 ，1 6 3 2 9 3 2 9 9 ． C h e n gG D ．I n t e r n a t i o n a lA c h i e v e m e n t so fs t u d yo n f r o z e ns o i lm e c h a n i c sa n de n g i n e e r i n g - - S u m m a r yo f t h ei n t e r n a t i o n a ls y m p o s i u mo ng r o u n df r e e z i n ga n d f r o s ta c t i o ni ns o i l s [ J ] ．A d v a n c ei nE a r t hS c i e n c e s ， 2 0 0 1 ，1 6 3 2 9 3 2 9 9 ． [ 3 3 金永军，杨维好．冻土墙用于基坑支护的力学设计方 法探讨[ J ] ．中国矿业大学学报，2 0 0 3 ，3 2 2 1 2 3 1 2 7 ． J i nYJ ，Y a n gW ．D i s c u s s i o no fd e s i g nm e t h o do f f r o z e n w a l la ss u p p o r t i n go ff o u n d a t i o np i t s [ J ] ． J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g T e c h n o l o g y ， 2 0 0 3 ，3 2 2 1 2 3 1 2 7 ． [ 4 3周国庆．间歇冻结抑制人工冻土冻胀机理分析[ J ] ．中 国矿业大学学报，1 9 9 9 ，2 8 5 4 1 3 4 1 6 ． Z h o uGQ ．A n a l y s i so fm e c h a n i s mo f r e s t r a i n i n go f s o i l f r e e z i n gs w e l l i n gb yu s i n gi n t e r m i s s i o nm e t h o d [ J ] ．J o u r n a l o fC h i n a U n i v e r s i t yo fM i n i n g ＆ T e c h n o l o g y 。1 9 9 9 ，2 8 5 4 1 3 - 4 1 6 ． [ 53 K o n r a dJM ，M o r g e n s t e r nNRA ．M e c h a n i s t i c t h e o r yo fi c el e n sf o r m a t i o ni nf i n e g r a i n e ds o i l s [ J 3 ． C a n ．G e o t e c h ，1 9 8 7 ， 1 7 4 7 3 4 8 6 ． [ 6 3R a z a q p u rAG ，W a n gDW ．F r o s t i n d u c e dd e f o r m a t i o n sa n ds t r e s s e si np i p e l i n e s [ J ] ．I n tP r e sV e s8 L P i p i n g 。1 9 9 6 6 9 1 0 5 1 1 8 ． [ 7 3陶文铨．数值传热学[ M ] ．西安西安交通大学出版 社，2 0 0 1 ． [ 8 3 郭宽良，孔祥谦。陈善年．计算传热学[ M ] ．合肥中国 科学技术大学出版社，1 9 8 8 ． [ 9 ]陆金甫．偏微分方程数值解法[ M ] ．北京清华大学出 版社，1 9 8 7 ． [ 1 0 3 孔祥谦．有限单元法在传热学中的应用[ M ] ．北京 科学出版社。1 9 9 8 ． 责任编辑陈其泰 万方数据