基于水动力学模型冻土冻胀数值模拟的改进.pdf

第3 5 卷第6 期 2 0 0 6 年1 1 月 中国矿业大学学报 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 文章编号1 0 0 0 - 1 9 6 4 2 0 0 6 0 6 - 0 7 6 2 - 0 5 V 0 1 ．3 5N o ．6 N O V ．2 0 0 6 基于水动力学模型冻土冻胀数值模拟的改进 商翔宇，周国庆，周金生 中国矿业大学建筑工程学院，江苏徐州2 2 1 1 1 6 摘要基于水动力学冻胀模型，对变化的边界条件下土体的冻胀量进行了大量计算．结果表明 由于水动力模型的控制方程中耦合源项在计算迭代步之内总是被不合理地估算，易于产生数值 振荡以致整个计算难以收敛；利用土壤冻结特性曲线及计算出的负温和未冻水含量之间的特定 关系，可以明确判断出迭代步计算所得耦合源项变化是偏大或偏小．据此，本文提出在迭代步内 直接对含冰量变化进行合理修正的改进算法，并编制了基于有限体积法的自调节时间步长计算 程序．该算法能够加速收敛、消除数值振荡，且与实测结果相当吻合，适用于复杂多变边界条件下 土体的冻胀计算． 关键词冻土；冻胀；数值模拟；算法改进 中图分类号T U4 4 5 文献标识码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 o z e n S o i l ’SF r o s tH e a v ew i t hH y d r a u l i c sM o d e l S H A N GX i a n g y u ，Z H O UG u o q i n g ，Z H O UJ i n - s h e n g S c h o o lo fA r c h i t e c t u r e8 LC 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 g T e c h n o l o g y ， X u z h o u ，J i a n g s u2 2 1 1 1 6 ，C h i n a A b s t r a c t B a s e do nt h eh y d r a u l i c sm o d e l ，al a r g en u m b e ro fc o m p u t a t i o n so ff r o z e ns o i l ’Sf r o s t h e a v ew i t hc h a n g e f u lb o u n d a r yc o n d i t i o n sw e r ec a r r i e do u t ．T h er e s u l t ss h o wt h a tb e c a u s et h e c o u p l i n gs o u r c et e r mi nt h em o d e lg o v e r n i n ge q u a t i o ni se s t i m a t e di m p r o p e r l y ，n u m b e rs u r g e s u n a v o i d a b l ya n dt h ec o m p u t a t i o ni sd i f f i c u l tt oc o n v e r g e ．T h em a g n i t u d eo fs o u r c et e r mc a nb e ju d g e dd e f i n i t e l yt ob eo v e r e s t i m a t e do ru n d e r e s t i m a t e du s i n gt h es p e c i f i cr e l a t i o nb e t w e e nm i n u st e m p e r a t u r e ，u n f r o z e nw a t e rc o n t e n ta n ds o i l f r e e z i n gc h a r a c t e r i s t i cc u r v e ．O nt h eb a s i so f t h i s ，a ni m p r o v e da l g o r i t h m ，w h i c ha d ju s t e dt h ei c ec o n t e n td i r e c t l yi ne a c hi t e r a t i v es t e pw a s p r o p o s e d ，a n dt h ec o r r e s p o n d i n gp r o g r a mw i t hs e l f - a d ju s tt i m es t e pa n df i n i t ev o l u m em e t h o d i sp r o g r a m m e d ．T h i sa l g o r i t h mc a nn o to n l ya c c e l e r a t et h ec o n v e r g e n ts p e e da n dr e m o v et h e n u m e r i c a lo s c i l l a t i o n ，b u ta l s oa p p l yt om a n yk i n d so fc h a n g e f u lb o u n d a r yc o n d i t i o n s ． k e yw o r d s f r o z e ns o i l ；f r o s th e a v e ；n u m e r i c a ls i m u l a t i o n ；a l g o r i t h mi m p r o v e m e n t 不论是寒区建设工程还是运用人工冻结施工 方案的地下工程都不可避免地遇到冻胀危害的威 胁．前者如我国东北和西北地区的各种土木工程建 设活动，例如青藏公路，每年都会因冻胀和融沉造 成道路翻浆致使交通受阻．后者如采用人工地层冻 结技术的隧道、城市地铁、超深基坑等工程‘卜2 1 ，如 瑞士苏黎世M i l e h b u c h 公路隧道第一冻结区测得 的冻胀量高达1 0 5m m [ 3 | ，日本大田干线隧道底面 收稿日期2 0 0 6 一0 1 1 2 基金项目国家重点基础研究发展规划 9 7 3 项目 2 0 0 2 C t M l 2 7 0 4 ，国家自然科学基金重点项目 5 0 5 3 4 0 4 0 ；国家自然科学基金项目 4 0 4 7 1 0 2 1 ；高等学校博士学科点专项科研基金项目 2 0 0 4 0 2 9 0 5 0 2 作者简介商翔宇 1 9 7 7 一 ，男，河南省洛阳市人，博士研究生，从事冻土工程和岩土力学方面的研究． E - m a l l x y s h a n g c u m t ．e d u ．c n T e l 0 5 1 6 - 8 3 9 9 5 0 7 8 万方数据 第6 期 商翔宇等基于水动力学模型冻土冻胀数值模拟的改进 实测最大冻胀量达到1 2 0m m ．而在冻结凿井工程 中冻胀量甚至达到10 0 01 T I l l 3 ．以上．这些过量的冻 胀会使地表产生不均匀变形，造成工程上部或邻近 建筑物、构筑物基础及煤气、电力、通讯及供排水等 管线的破坏H ] ，甚至影响到工程的成败．因此在进 行这些工程建设之前必须尽可能准确地预测可能 冻胀量，是保证工程成功实施的先决条件之一． 国内外学者为预测冻胀量，不断深化对冻胀的 研究水平，提出了定量描述和解释冻胀现象主要特 征的各种变量之间的关系式即许多冻胀预测模型， 并逐渐得到验证[ 5 ] ．根据所依据的冻胀理论的不 同，可以将冻胀预测模型分为4 大类第一冻胀理 论模型、第二冻胀理论模型、热力学模型和水动力 模型．这里限于篇幅不再对各种模型做详尽描述， 只对本文研究对象水动力冻胀模型作简单介绍． 1 水动力冻胀模型简介 水动力模型是R ．L ．H a r l a n 基于土体即使在 冰点以下数十度依然有未冻水存在的认识于1 9 7 3 年提出的，直到上世纪末关于水动力学模型的研究 依然在进行[ 6 。8 ] ，这足见其生命力之强大．该模型认 为土壤冻结过程主要发生3 种物理过程热量的传 递、水分的迁移和水分的液一固相变．这3 个过程 同土壤中的水分场和温度场相互作用、相互制约、 相互影响．模型不考虑冰透镜体的形成与冻结缘的 发展，亦不考虑上覆荷载的影响，假设冻胀形成的 临界值．该模型能很好地模拟土壤中的水热状况． 其控制方程如下 c ，等一麦 Aa ⋯te L 鲁， 1 孥一麦 脚。卺 警一瓦P i 鲁，c 2 ， 乱≤9 T T 0 ， 3 式中C 。为土壤体积热容量，J ／ m 3 K ；A 为热导 率，W ／ m K ；L 为融化潜热，J ／m 3 ；T 为温度，K ； D 为土壤水扩散率，1 3 3 ．2 ／s ；K 为导水率，m ／s ；9 。为 水含量；馈为冰含量；2 为空间坐标，m ；t 为时间，s ； p T 是与T 相对应的未冻水含量． 水动力模型旨在解决含相变的水热耦合问题， 其数学表示是非线性的抛物型耦合方程组．求解该 方程组只能采取数值方法，求解非常复杂，关于该 问题的求解特点在现有文献中已有较多论述[ 9 ] ．认 为相变界面上吸收和释放的热量由于相变潜热的 量级比相邻单元传递的热量大的多，该热量如同一 个内部脉冲源，造成界面上能量的突变和能量守恒 条件的强非线性，极易造成数值振荡，是该问题求 解困难的主要原因．换言之，含相变的水热耦合问 题是由于水热方程中的含冰量即源项 也即方程组 的耦合项 的存在使得计算不易达到收敛口0 | ．每一 计算时间步内最困难的是源项具体数值的确定，因 此不少改进的算法都是围绕如何确定含冰量变化 进行的．表1 简要总结评价现有文献中关于土壤冻 标志为孔隙冰含量超过一个不超出土体孔隙度的 结条件下水热耦合问题的算法． 表1现有水动力模型算法评价 T a b l e1C o m m e n to nt h er e p o r t e dh y d r a u l i c sm o d e l s 由表1 可见，现有的算法改进几乎都是朝着降 低方程耦合性的方向进行的，而方程耦合性的降低 则是以对另外一些要求的增高为代价的，比如对土 壤冻结曲线的连续性从一阶提高到了两阶．因为目 前研究冻胀控制及人工冻结技术研究领域，经常遇 到边界温度变化幅度和速度俱快的情况，这种情况 下，震荡更加剧烈，而且计算精度要求亦不能降低． 因此，能在不降低方程本身耦合性的同时，寻求到 万方数据 中国矿业大学学报 第3 5 卷 提高计算收敛速度和降低算法震荡的方法无疑是 有价值的． 2 算法改进 2 ．1 方程离散 首先采用控制容积法对控制微分方程进行离 散，对离散方程则按有限差分方程来讨论其稳定性 和精度．这里略去离散公式的具体推导，只作简要 说明．空间离散采用“先单元后结点”的方法，该法 能保证结点一定位于所代表单元的几何中心；时间 域离散采用稳定精度比较好的全隐格式．得到的是 全隐格式的空间二阶精度格式．下列各式中下标 J ，挖分别为空间和时间坐标标记． 整理得到热方程的离散方程 A ，T V - } B ，T 广1 G T j 计 1 1 一E j ， 4 式中A i 一土堕；c i 一土堕； B ，一一A i C f 一 c ， 广1 ； E j 一一 c ， 广1 T ；一L [ 纯 广1 一 吼 ；] ； A 鹞一鲣j - - 1 嘻业；A 鹞鳢≥业． 水分方程的离散方程 J i 9 声} G P 广1 H j P 搿一I j ， 5 E D i 式中J J 一_ ‘- ；G 一一- ，J 日』一1 ； r 一，．， z rE D 葛F A t 马 』g j l - - g j ；E 一万1 等1 一丽2 oJ ； 。 。 z ， 一一2 f 一一 LJ L 一一孵 E [ K 鑫一K 蔓] 譬[ 仍 尹1 一 9 i ] ． 本文研究的水热耦合方程属于初边值问题，其 定解条件包括初始条件和边界条件2 类．本文的初 始温度条件可以是均匀分布的，也可以是非均匀分 布的；初始含水量只考虑均匀分布的情况．边界条 件包括热边界和水分边界条件，热边界条件仅考虑 第1 类边界条件，即边界温度已知的情况；水分边 界条件考虑第1 和第2 两类边界，即边界上含水量 已知和通过边界的水分流量已知的情况． 2 ．2 改进思路 作者大量的计算实践表明每一迭代步之内含 冰量的变化很难确定，其主要原因是迭代过程中温 度在冻结点上下跳动的单元含冰量的变化是质的 变化，即从有到无或从无到有，而不像其它点 一般 是远离冻结锋面的点 含冰量的变化仅仅是从多到 少或相反，是量的变化．尽管采用超松驰迭代这样 具有平均光滑意义的加速收敛方法可以使得方程 收敛，但数值振荡依然无法完全避免．因此如果能 为含冰量的修改像最小梯度法那样指定一个明确 的修改方向，收敛可能就会加快，而且振荡会大大 降低． 事实上，产生计算振荡的直接原因是含冰量修 正不当引起的．这种不当与计算出来的未冻水含量 和温度值之间有明确的关系，而这种关系提供了含 冰量修改方向的线索．可以设想，如果本计算时间 步得出的含冰量偏大 与真实合理的情况相比较 ， 则由水分方程计算出的未冻水含量则必然偏小，由 热方程计算出的温度值则会偏高；反之，则计算出 的未冻水含量和温度则分别会偏大和偏低．基于这 种关系，可以得出如下结论若按照冻结特性曲线 计算出的未冻水含量比由水分方程计算出的大，则 意味着原来计算的含冰量偏大，反之则偏小．迭代 过程中以由冻结特性曲线计算出的未冻水含量与 由水分方程计算出的未冻水含量之差作为控制误 差，进行含冰量修正；该误差值大于零就减小原来 的冰含量，反之就增大之．在迭代计算过程中根据 上述关系不断修正含冰量，直至控制误差达到允许 的范围之内．可以看到上述这种直接修改含冰量的 做法，由于基于含冰量和所属控制误差之间的明确 关系，为迭代过程中的含冰量的修改指出了确定的 方向，有利于迭代计算的收敛和稳定．另外，如果在 第一个迭代步内不考虑本时间步可能产生的新含 冰量，则根据以上关系计算出的含冰量肯定是偏大 的，那么该计算就估计出了本时间步可能产生的最 大量的含冰量和可能的相变发生范围，余下的迭代 步就可以根据上述关系不断地修正含冰量，直至控 制误差达到允许的范围之内．本文将基于这一思路 对现有算法进行修正． 2 ．3 算法实现 本文的计算实践表明空间步长加密到一定程 度只能增多计算量，而无助于计算精度的大幅提 高，为了兼顾计算效率和计算精度，本文选取等间 距空间步长0 ．5c m 次最后步长为0 ．1 8e m ．时间 步长的选择从两方面考虑进行自调节时间步长，一 是边界条件的考虑，另一是计算过程中达到收敛的 迭代次数．如果边界条件在某个时间段内变化剧 烈，那么时间步长就不能过大，这是首先要保证的； 如果收敛较慢的话，则需要适当减小时间步长．如 果边界条件变化缓慢，且收敛情况较好，则时间步 长可考虑大幅的增大．因此基于以上考虑的时间步 长选择方法就可以兼顾计算速度和准确性2 个方 万方数据 第6 期商翔宇等基于水动力学模型冻土冻胀数值模拟的改进 面在变化缓慢的时段计算时间步长不会过短以至 于影响计算速度，在需要时间步长小以保证计算精 度的时段时间步长不会过长以至于影响计算精度． 算法用M A T L A B 语言进行编程实现 见图1 ． 输入基本参数初始条件；水热参数；时空参数 N 计算时间步乒O ；计算总时间t t O J ；i ， 1 f 户“ △t 1 1 k f T ．1 t t t t - A t 减小时间步长 设△妒， 0 ，进行初算求出T l 及q J ．1 ，并由冻结特性曲线求出△仍 Y { | 燮耋堡攀篓 尚∑ 1 T ≤2 0 I 根据P 的大小对△卿进行调节 赋时间步末值为下一时间步的初值 ≮堡塑璺卿 翻 图1计算流程 F i g ．1 F l o w c h a r to ft h ea l g o r i t h m N 2 ．4初步验证 图2 比较了用松弛迭代法直接进行耦合计算 和用本文方法计算出的冻结深度随时间的变化曲 线．采用的算例是文献[ 1 1 ] 提供的，这里给出初始 含水量为0 ．2 ，模拟3 ．0m 垂直土柱冻结的情况． 1 2 g 1 0 萤 暴4 U1 0 02 0 03 0 04 0 05 0 06 0 07 0 08 0 0 t f m i n 图2冻结锋面随时间变化曲线 F i g ．2C h a n g i n gc u r v e so ff r e e z i n gf r o n tw i t ht i m e 从图2 可以看出本文提出的算法把数值振荡 基本消除了．而且，由于收敛性能好，该算法使得模 型对于高含水量和低温度边界的适应性大大提高 了，表现出较好的计算性能． 3算法验证 为了验证本文提出算法的可靠性，这里采用了 文献[ 1 2 ] 提供的室内试验原始数据进行了计算，并 与该文献提供的试验结果和计算结果进行了比较． 比较结果表明本文算法是可靠的． 该计算仅经过了2 5 1 个时间步，历时不足1 m i n ．计算结果给出模拟试验结束时，含水量的分 布情况和试验过程中3 4 ，5 2 8 和28 3 0m i n 时刻的 温度分布情况，见图3 ．结果表明计算结果与实测 结果接近且趋势一致．表2 比较了计算结束时本文 计算结果与实测结果及文献E 1 2 ] 计算结果关于最 大冻结深度和土柱内最大与最小含水量值． 深度，m 深度／m a 计算与实测土样内含水量分布对比 b 计算与实测温度场对比 图3水分场与温度场的计算结果与实测数据对比 F i g ．3C o m p a r i s o no fc o m p u t e da n de x p e r i m e n t a ld a t ao fm e i s t u r ec o n t e n ta n dt e m p e r a t u r ef i e l d s 裹2试验结束时计算结果的对比 T a b l e2 C o m p a r i s o no fn u m e r i c a la n de x p e r i m e n t a lr e s u l t sa tt h ee n do ft h et e s t 数据源最大冻结深鼽m 土譬警量麓揪土翟筒量麓裂焉 作者还做了封闭和开放系统2 种情况下的数 值模拟，模拟结果较好地反映了2 类问题的冻胀特 征m ] ，限于篇幅这里不再给出计算结果．作此说 明，以表明该算法对不同边界条件的良好适应性． 万方数据 7 6 6中国矿业大学学报 第3 5 卷 4 结论 本研究表明水动力冻胀计算模型求解难以收 敛的主要原因在于方程中的源项，即每个时间步的 体积含冰量难以给出合理的估计；根据计算过程中 未冻水含量、负温及体积含冰量大小之间的动态联 系，可以给出合理估计或修正体积含冰量大小的思 路．基于此提出了在初定的相变范围之内直接不断 调整体积含冰量大小的改进算法，计算模拟与实测 结果比较吻合．计算实践表明该方法不但可以加速 收敛过程，消除数值振荡，而且对于多种不同的边 界条件均具有良好的适应性和敏感性．该方法不仅 可以适宜于复杂边界条件下的冻土冻胀数值模拟， 而且也会有助于改进其它类型包含相变问题的数 值模拟． 致谢感谢清华大学水利系尚松浩副教授和广东省 水利水电科学研究院王海丽老师的指导；感谢中国 矿业大学青年科研基金项目 2 0 0 5 A 0 0 5 对本文研 究的资助． 参考文献 [ 1 3 金永军，杨维好．冻土墙用于基坑支护的力学设计方 法探讨[ J 1 ．中国矿业大学学报，2 0 0 3 ，3 2 2 1 2 3 1 2 7 ． J I NY o n g ju n ，Y A N GW e i h a o ．D i s c u s s i o no fd e s i g n m e t h o do ff r o z e nw 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 ． [ 2 ] 王衍森，蒋武军，庞法扬，等．某深基坑冻土墙应力 与变形的数值计算[ J ] ．中国矿业大学学报，2 0 0 2 ，3 1 3 2 8 5 2 8 8 ． W A N G , Y a n - s e n ，J I A N GW u j u n ，P A N GF a - y a n g ， e ta 1 ．N u m e r i c a la n a l y s i so fs t r e s sa n dd e f o r m a t i o no f f r o z e nw a l la r o u n dd e e pf o u n d a t i o np i t [ J ] ．J o u r n a lo f C 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 2 ，31 3 2 8 5 - 2 8 8 ． [ 3 ] M E T T I E RK ．瑞士苏黎世M i l c h b n c k 公路隧道的冻 结施工[ C ] ／／K I N O S I T AS ，F U K U D AM ．P r o c e e d i n go f4 t hI n t e r n a t i o n a lS y m p o s i u monG r o u n dF r e e z i n g ．R o t t e r d a m ，N e t h e r l a n d s AAB A L K E M A ， 1 9 8 5 8 3 8 8 ． [ 4 ] 陈瑞杰，程国栋，李述训，等．人工地层冻结应用研 究进展和展望[ J ] ．岩土工程学报，2 0 0 0 ，2 2 1 4 0 4 4 ． C H E NR u i j i e ，C H E N GG u o d o n g ，L iS h u - x u n ，e t a 1 ．D e v e l o p m e n ta 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 lg r o u n df r e e z i n g [ J ] ．C h i n e s eJ o u r n a l o fG e o t e c h n i c a lE n g i n e e r i n g ．2 0 0 0 ，2 2 I 4 0 4 4 ． [ 5 ] 徐学祖．冻土物理学[ M ] ．北京科学出版社，2 0 0 1 ． [ 6 ] 王海丽．冻土水热运动的数值模拟[ J ] ．内蒙古农牧学 院学报，1 9 9 8 ，9 1 9 9 1 0 5 ． W A N GH a i _ l i ．N u m e r i c a ls i m u l a t i o no fm o i s t u r ea n d h e a tt r a n s f e rd u r i n gs o i lf r e e z i n g [ J ] ．J o u r n a lo fI n n e r M o n g o l i ai n s t i t u t eo fa g r i c u l t u r e ＆a n i m a lh u s b a n d r y ，1 9 9 8 ，9 1 9 9 1 0 5 ． [ 7 ] S U K N A MK ．C o u p l e dh e a ta n dm o i s t u r ef l o wa n a l y s i si nu n s a t u r a t e ds o i l [ D ] ．D e p a r t m e n to fC i v i lE n g i n e e r i n g ，U n i v e r s i t yo fT o l e d oo fC a n a d a ，2 0 0 2 ． [ 8 3 s A L L YAS 。S U A NRB ．M o i s t u r em i g r a t i o nd u r i n g f r e e z ea n dt h a wo fu n s a t u r a t e ds o i l s m o d e l i n ga n l a r g es c a l ee x p e r i m e n t s [ J ] ．C o l dR e g i o n sS c i e n c ea n d T e c h n o l o g y ，1 9 9 7 ，2 5 3 3 4 5 ． [ 9 3 陶文铨．数值传热学[ M ] ．西安西安交通大学出版 社，2 0 0 1 ． [ 1 0 ]尚松浩，雷志栋，杨诗秀．冻结土壤冻结条件下水热 耦合运移的数值模拟的改进[ J ] ．清华大学学报， 1 9 9 7 ，3 7 8 6 2 6 4 ． S H A N GS o n g h a o ，L E IZ h i d o n g ，Y A N GS h i x i u ． 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 fc o u p l e dm o i s t u r ea n dh e a tt r a n s f e rd u r i n gs o i lf r e e z i n g [ J ] ．J o u r n a lo fT s i n g h u aU n i v e r s i t y ，1 9 9 7 ，3 7 8 6 2 6 4 ． [ 1 1 ] 杨诗秀，雷志栋，朱强．土壤冻结条件下水热耦合运 移的数值模拟[ 刀．清华大学学报，1 9 8 8 ，2 3 增1 11 2 1 2 0 ． Y A N GS h i - x i u ，L E IZ H I - d o n g ，Z H UQ i a n g ．N u m e r i c a ls i m u l a t i o nf o rc o u p l e dm o i s t u r ea n dh e a t t r a n s f e rd u r i n gs o i lf r e e z i n g [ J ] ．J o u r n a lo fT s i n g h u a U n i v e r s i t y ，1 9 8 8 ，2 3 s u p p ．1 1 1 2 1 2 0 ． [ 1 2 ] 胡和平，杨诗秀，雷志栋．土壤冻结时水热迁移规律 的数值模拟[ J ] ．水利学报，1 9 9 2 7 1 - 8 ． H uH e - p i n g ，Y a n gS h i x i u ，L E IZ h i d o n g ．An u m e r i c a ls i m u l a t i o nf o rh e a ta n dm o i s t u r et r a n s f e r d u r i n gs o i lf r e e z i n g [ J ] ．J o u r n a lo fh y d r a u l i ce n g i n e e r i n g ，1 9 9 2 7 1 - 8 ． [ 1 3 ] 商翔宇．冻土冻胀与冻结模式关系的试验与数值模 拟研究[ D ] ．徐州中国矿业大学建筑工程学院， 2 0 0 4 ． 责任编辑 王继红 万方数据