第37卷第5期2023年10月水土保持学报J o u r n a l o f S o i l a n d W a t e rC o n s e r v a t i o nV o l .37N o .5O c t .,2023收稿日期:2023-03-07资助项目:国家自然科学基金项目(U 20A 20110 01);中国地质调查局地质调查项目(D D 20221746) 第一作者:杨昶(1997 ),男,硕士,主要从事泥石流数值模拟技术研究㊂E -m a i l :y a n g c h a n g21@m a i l s .u c a s .a c .c n 通信作者:铁永波(1979),男,博士,教授级高级工程师,博士生导师,中国地质调查局杰出地质人才, 重点地区特大地质灾害链调查评价 工程首席专家,主要从事地质灾害评价与防治研究㊂E -m a i l :t y b 2009@q q.c o m 我国泥石流数值模拟研究历程及发展趋势杨昶1,2,铁永波2,3,4,张宪政2,3,4,向炳霖2,3,4,5(1.中国地质科学院,北京100037;2.中国地质调查局成都地质调查中心(西南地质科技创新中心),成都610081;3.自然资源部地质灾害风险防控工程技术创新中心,成都611734;4.自然资源部成都地质灾害野外科学观测研究站,成都610000;5.重庆交通大学河海学院,重庆400074)摘要:针对泥石流流变特性和动力学过程复杂的特点,在系统收集国内泥石流数值模拟研究相关文献的基础上,按初步探索㊁逐渐完善和成熟多元3个阶段对我国泥石流数值模拟研究的历程及现状进行了总结㊂根据物质组成和动力学特征来看,泥石流的动力学模型可划分为连续介质㊁离散介质和混合介质模型,通过分析不同模型的特点及应用场景,对比总结各类模型和数值计算方法的适用性及不足,在此基础上针对当前泥石流模拟存在的难点,对泥石流数值模拟的发展进行了展望,旨在为泥石流数值模拟研究及应用软件开发提供参考㊂关键词:泥石流;动力学模型;数值模拟;研究综述中图分类号:P 694 文献标识码:A 文章编号:1009-2242(2023)05-0001-11D O I :10.13870/j.c n k i .s t b c x b .2023.05.001R e s e a r c hH i s t o r y a n dD e v e l o pm e n t T r e n do fN u m e r i c a l S i m u l a t i o no fD e b r i s F l o w i nC h i n aY A N GC h a n g 1,2,T I EY o n g b o 2,3,4,Z HA N G X i a n z h e n g 2,3,4,X I A N GB i n gl i n 2,3,4,5(1.C h i n e s eA c a d e m y o f G e o l o g i c a lS c i e n c e s ,B e i j i n g 100037;2.C h e n g d uC e n t e r o f Ch i n a G e o l o g i c a lS u r v e y (S o u t h w e s tG e o l o g i c a lS c i e n c e a n dT e c h n o l o g y I n n o v a t i o nC e n t e r ),C h e n gd u 610081;3.Te c h n o l o g y I n n o v a t i o nC e n t e rf o rR i s kP r e v e n t i o na n d M i t ig a t i o no f G e oh a z a r d ,Mi n i s t r y o f Na t u r a l R e s o u r c e s ,C h e n g d u 611734;4.Ob s e r v a t i o na n dR e s e a rc hS t a t i o no f C h e n gd uGe o l o g i c a lH a z a r d s ,M i n i s t r y o fN a t u r a lR e s o u r c e s ,C h e n g d u 610000;5.H e h a iC o l l e g e ,C h o n g q i n g J i a o t o n g U n i v e r s i t y ,C h o n g q i n g 400074)A b s t r a c t :A i m i n g a t t h e c o m p l e x c h a r a c t e r i s t i c s o f r h e o l o g i c a l p r o p e r t i e s a n d d yn a m i c p r o c e s s o f d e b r i s f l o w ,o n t h eb a s i s o f s y s t e m a t i c a l l y c o l l e c t i n g re l e v a n t l i t e r a t u r e s o nn u m e r i c a l s i m u l a t i o nof d e b r i s f l o wi nC h i n a ,t h e h i s t o r y a n d c u r r e n t s i t u a t i o no f n u m e r i c a l s i m u l a t i o no f d e b r i s f l o w i nC h i n a a r e s u mm a r i z e d a c c o r d i ng t o th r e es t a g e s :p r e li m i n a r y e x p l o r a t i o n ,g r a d u a li m p r o v e m e n t a n d m a t u r i t y .A c c o r d i n g t o t h e m a t e r i a l c o m p o s i t i o na n dd y n a m i cc h a r a c t e r i s t i c s ,t h ed y n a m i c m o d e lo fd e b r i s f l o wc a nb ed i v i d e d i n t oc o n t i n u o u s m e d i u m ,d i s c r e t e m e d i u m a n d m i x e d m e d i u m m o d e l .B y a n a l y z i n g t h ec h a r a c t e r i s t i c sa n d a p pl i c a t i o n s c e n a r i o s o f d i f f e r e n tm o d e l s ,t h e a p p l i c a b i l i t y a n d s h o r t c o m i n gs o f v a r i o u sm o d e l s a n dn u m e r i c a l c a l c u l a t i o n m e t h o d s a r e c o m pa r e d a n d s u mm a r i z e d .O n t h i sb a s i s ,i n v i e wo f s o m e d i f f ic u l t i e s i n t h ede b r i sf l o ws i m u l a t i o n ,t h e f u t u r e d e v e l o p m e n t o f d e b r i s f l o wn u m e r i c a l s i m u l a t i o n i s p r o s p e c t e d ,a i m i ng t o p r o v i d i n g r e f e r e n c e f o r d e b r i s f l o wn u m e r i c a l s i m u l a t i o n r e s e a r c ha n d a p p l i c a t i o n s o f t w a r e d e v e l o p m e n t .K e y w o r d s :d e b r i s f l o w ;d y n a m i cm o d e l ;n u m e r i c a l s i m u l a t i o n ;r e v i e w 泥石流是一种常见于山区由重力作用驱动的典型自然灾害,具有突发性强㊁流动速度快和破坏能力显著的特点㊂随着我国经济水平的不断发展,人们开始投身于山区的开发利用,与之而来的泥石流灾害问题也进入人们的视野[1]㊂早期对泥石流的研究过程通常为野外观测ң数据采集ң经验总结[2]㊂但这一思路通常只能对某些问题进行定性的分析,无法定量说明泥石流运动过程中的特有现象,如阵性来流㊁高Copyright ©博看网. All Rights Reserved.速低阻等[3]㊂为深入研究泥石流的运动规律㊁侵蚀堆积过程及冲击破坏能力,我国学者着眼于宏观和细观尺度,采用相应的数值技术对泥石流开展从单尺度到多尺度㊁从连续介质模型到耦合模型的研究历程,推动泥石流的研究体系逐步由定性向定量转变㊂泥石流的数值模拟最初采用单流体模型和多流体模型㊂单流体模型是指将泥石流视为整体,建立单一组分的动力学模型㊂多流体模型考虑泥石流固液两相之间的相互作用关系,通过两相动量间的动量交换建立起多组分的动力学模型㊂常用的单流体模型有B i n g h a m 模型[4]㊁V o e l l m y模型[5]和膨胀体模型[6],多流体模型有库伦混合流模型[7]和广义两相流模型[8]㊂单流体模型虽然计算参数少㊁计算量小且使用方便,但其假设具有局限性,即认为泥石流体的固㊁液两相之间速度基本一致,忽略两相间的相对运动㊂这在黏性泥石流中大致成立,但对稀性泥石流和水石流则明显与实际不符㊂因此,为研究固㊁液两相之间的相互作用关系,有学者借鉴水和细颗粒泥沙的两相流模型,提出浆体和粗颗粒泥沙/石块的结构两相流模型,并进一步对高浓度和低浓度的泥石流运动机理进行分析[9-12]㊂该模型考虑固液两相之间的动量交换,采用两组分的流体动力学方程进行描述,相比于单流体模型更加适用于黏性泥石流,本质上结构两相流模型是多流体模型的一种㊂随着研究的深入,学者们发现无论是单流体模型还是多流体模型,对固相的处理都存在不足,尤其是对宽级配泥石流体的粗颗粒组分难以进行定量分析㊂尽管多流体模型将固㊁液两相分别采用2组方程进行描述,但连续介质模型从本质上来说并不适用于对粗颗粒物质进行描述[13]㊂因此,适用于宽级配固相颗粒的基于离散介质模型的数值计算方法引起广泛关注㊂离散介质模型大体可分为从微观角度出发建立的离散元模型和基于介观建立的动力学模型,常用的数值计算方法有离散单元方法(D E M)[14-17]㊁格子玻尔兹曼方法(L B M)[18-19]和元胞自动机方法(C A M)[20]㊂离散元模型适用于颗粒流和碎屑流的动力学机理研究,动力学模型适用于固液碰撞和多孔渗流等问题的研究㊂总体来讲,2种模型均适用于细 微观尺度下的动力学研究,但针对的研究问题各有侧重㊂对宏观尺度的流动问题,因计算效率和收敛性等原因,通常采用连续介质模型进行处理㊂为了兼顾连续介质模型和离散介质模型的优点,混合模型由此产生㊂这一模型不仅解决连续介质模型在描述粗颗粒及相间作用力上的困难,同时也避免离散介质模型对宏观流动问题计算效率的不足㊂相比于前2类模型,混合模型更为符合泥石流的物理本质,目前已有的混合模型计算方法有S P H-D E M[21]㊁L B M-D E M[22]和C F D-D E M[23]㊂尽管混合模型在物理上更符合泥石流的特性,但在实际使用过程中存在较多难点,比如,各相间相互作用力难以精确量化㊁宽级配泥石流体的多尺度问题及耦合算法的界面信息传递等问题,其理论和计算方法都有待进一步完善㊂本文从我国泥石流数值模拟的研究历程出发,聚焦不同类型的泥石流动力学问题,对比分析不同数值模拟方法的优缺点,总结目前取得的成果,并对下一步的研究提出建议,旨在为快速㊁全面了解泥石流动力学模型和数值计算方法提供参考㊂1研究历程1.1初步探索阶段国内对泥石流数值模拟的研究始于20世纪90年代㊂这一时期,由于人类经济和社会的快速发展,所面临的自然灾害问题也日益严重㊂了解泥石流的运动规律㊁分析其冲出堆积过程㊁预测其危险程度是当时亟待解决的科学问题㊂为研究上述问题,唐川等[24-26]基于守恒定理推导出一维非定常流的控制方程,进而将其扩展至二维㊂通过引入隐式剖开算子,对方程组进行求解,并结合云南怒江的芭蕉河泥石流进行实例验证㊂数值和实测结果表明,这一方法能有效模拟泥石流的冲出堆积过程(图1)㊂基于上述的数值结果,以泥深和流速为评价指标,提出泥石流堆积扇的危险度评价体系,开创性地采用数值模拟手段进行泥石流堆积区危险性评价㊂图1芭蕉河泥石流泛滥的模拟范围与实际范围比较[24]泥石流由高浓度的粗颗粒和水或浆体组成,其内部存在由粗颗粒运动导致的碰撞和摩擦力,粗颗粒和浆体之间的湍流力,孔隙流体或颗粒和浆体混合物的黏性力,以及两相之间相对运动产生的作用力[27-29]㊂随着流体中平均粗颗粒体积浓度(C)的变化,其主要作用力不同,结合雷诺数㊁B a g n o l d数和相对深度,可以对流体进行相应划分[27](图2),其流变特性及相应的力学模型见表1㊂2水土保持学报第37卷Copyright©博看网. All Rights Reserved.注:附图修改自参考文献[27]㊂图2 固液两相流体运动形式判据表1 不同流体的剪切力模型流体本构模型模型名称模型方程式方程编号水流牛顿流变模型τ=μγ[30-31](1)挟沙水流低含沙水流采用公式(1)[32]高含沙水流采用公式(2),但需要对μ进行修正[33]B i n g h a m 模型τ=τ0+μγ[34-35](2)V o e l l m y 模型τ=σf +ρg v 2ξ[36-37](3)泥石流膨胀模型τ=τ0+μγ+ξγ2[38-39](4)H e r s c h e l -B u l k l e y 模型τ=τ0+μγn [40-41](5)双线性流变模型τ=τ0+μγ+τ0γ0γ+γ0[40-41](6)库伦混合流模型τ=c +σt a n φ+μγn [42](7) 注:τ为剪切力;τ0为屈服应力;μ为黏性系数;γ为剪应变;γ0为临界剪应变;σ为正应力;f 为摩擦系数;φ为内摩擦角;ρ为密度;g 为重力加速度;v 为流向速度;ξ为湍流系数;n <1㊂ 泥石流动力学方程大多采用不可压缩流体的N S 方程,或将N S 方程在深度上平均简化后得到S a i n t -V e n a n t 方程㊂上述方程组的求解过程可以借鉴不可压缩流体C F D 的计算思路,通过有限差分/有限体积法/有限元方法对N S 方程进行离散,而后采用S I M P L E /P I S O /P I M P L E 算法求解计算区域的流场和压力场㊂基于上述思路,可以分析出不同边界条件下,泥石流区别于水流和挟沙水流的物理特性及流场差异,有助于理解泥石流的运动特点[43]㊂单流体模型便于使用且应用广泛,但其假定泥石流固液两相之间强耦合,不具有相对流动速度,忽略泥石流的固㊁液两相流特征㊂一种更为合理的方法是将泥石流视为一种特殊的结构两相流体,固相由粗颗粒部分组成,液相由细颗粒和水组成[9-10]㊂在此基础上针对泥石流运动堆积过程的特点,将其简化为平面二维流动㊂泥石流的停止是重力势能转化为阻力做功的结果,并忽略密度在流动过程中的沿程变化和动量扩散影响,假定速度在垂向上变化很小,可忽略,两相之间的相互作用通过阻力项进行描述,进而体现相间动量交换过程,最后可以得到结构两相流的动力学方程组,并用于重现1991年沙湾大沟泥石流的运动堆积过程(图3)[44-46]㊂上述的结构两相流方程不仅可以模拟泥石流的运动堆积过程,通过特定的简化和采取合适的阻力模型,结合P I C 算法可模拟稀性和黏性泥石流的阵性运动和堆积过程[47-48],也可用于分析不同流量补给条件下的泥石流堆积特点[49]㊂初期泥石流的数值模拟从动力学模型到数值求解方法,基本上沿袭水力学和计算流体力学2个学科的思路,对泥石流物理模型上的认知更多处于单流体和多流体模型阶段,即总体上认为泥石流的运动以液相运动为主,符合连续介质假设㊂尽管单流体和多流体模型无法描述以固相颗粒运动为主的泥石流体,但其简化方便㊁容易理解,从实际使用效果来看,足以满足工程实践的需求,直至今日仍广泛应用㊂以上工作在我国泥石流数值模拟发展的初期起到奠基作用,对泥石流的动力学方程及其计算方法进行了有益探索㊂3第5期 杨昶等:我国泥石流数值模拟研究历程及发展趋势Copyright ©博看网. All Rights Reserved.图31991年沙湾大沟泥石流运动堆积过程的数值重现[46] 1.2逐渐完善阶段从2000年开始,三峡库区及震区的灾害防治理论与技术研究,极大地推进防灾减灾技术的发展[50]㊂基于长期的实践经验发现,泥石流的运动过程难以用简单的流动ң堆积进行概述㊂由于不规则地形的影响,泥石流的运动过程常伴随着汇聚与分流㊁滞留与侵蚀等复杂现象[51]㊂其运动过程包含启动㊁沿程运动(侵蚀放大㊁颗粒分选)和汇流3个阶段(图4),因此需针对各阶段的特点进行相应研究㊂注:(a)修改自参考文献[52];(b)修改自参考文献[53];(c)修改自参考文献[54];(d)修改自参考文献[55]㊂图4泥石流启动到入汇主河的动力学过程示意泥石流的启动通常与地震和降雨密切相关,其触发过程通常可以用驱动力(地震㊁降雨㊁人类活动等)和抗力(土体的抗剪强度)的相互作用来描述[56]㊂泥石流的启动机理通常可按土力和水力两类泥石流进行研究㊂土力泥石流多为降雨入渗下,土体饱和失稳液化形成的一类泥石流,其启动因子为土体水分饱和度和孔隙水压力(图4c);水力泥石流是受水流作用控制的一类泥石流,由地表径流对沟道物源进行冲刷形成,启动因子为流体的冲击力和拖拽力[57]㊂土力泥石流的数值研究方法通常采用离散元方法,通过研究不同含水率和降雨强度条件下的泥石流启动过程发现,土体失稳时的含水率往往在饱和含水率附近[54,58-59]㊂水力4水土保持学报第37卷Copyright©博看网. All Rights Reserved.泥石流的启动则可参考沟道侵蚀过程㊂泥石流沿沟道运动过程中,常伴随着颗粒分选和沿程侵蚀过程㊂泥石流在沟道中的运动过程可以用明渠流进行类比㊂已有研究[11-12]表明,管道和明渠中的固液两相流在近壁区存在显著的滞后现象,为探究其深层原因,有学者从动理学角度出发,采用F o k k e r-P l a n c k扩散模型描述湍流与颗粒间的相互作用,进一步引入修正的B G K模型对颗粒的碰撞效应进行近似处理,采用C h a p m a n-E n s k o g展开将复杂的动力学方程转化为线性方程形式,并建立描述颗粒相运动的动力学模型㊂采用上述模型,对高/低浓度的固液两相流进行研究发现,颗粒脉动能与作用在颗粒上的外力是影响颗粒运动的主要影响因素[11-12,60-61]㊂上述研究表明,泥石流在运动过程中存在颗粒分选效应,粗颗粒逐渐向泥石流浅表层及前缘聚集,细颗粒则迁移至泥石的下层及底部(图5b和图5c),前者是粗颗粒的向上分选过程,后者为细颗粒的向下筛分过程,最终形成泥石流的堆石垄结构(图5a)[16,62-63]㊂此外,Z h o u等[64]将泥石流视作二元颗粒流结合C F D-DE M方法提出颗粒分选模型,研究完全浸没在不同环境流体中的颗粒流分选情况表明,黏性流体的存在有效地降低颗粒分选效应,但并非直接进行影响,而是通过降低剪切速率和压力的形式表现㊂注:(a)修改自参考文献[62];(b)和(c)修改自参考文献[63];箭头方向指标重力方向㊂图5运动泥石流形成堆石垄过程中的颗粒分选机理泥石流的沿程侵蚀过程(图4d)常被描述为 滚雪球 效应[65],如何考虑泥石流在运动过程中由于侵蚀作用导致的体积放大效应是一大难点㊂解决上述问题,通常有3种思路:(1)基于试验结果和自然物理现象(例如,一次泥石流事件发生前后的体积变化)总结得到的经验公式;(2)基于平衡状态假设推导得出的侵蚀模型;(3)基于守恒方程推导得出的侵蚀模型㊂例如,较多研究[66-73]基于上述思路开发出的侵蚀模型,更多国内外学者提出的不同侵蚀模型见表2㊂在复杂地形条件下,局部的地形变化可能对侵蚀作用产生显著的影响㊂已有研究[74-76]表明,地形起伏强烈影响着流体的运动,主要表现为剪切作用和动量传递过程在地形阻碍或重定向运动的地方存在突变,其本质是复杂地形引起的垂直加速度和曲率效应的影响[77]㊂总体趋势为地形整体起伏变化较小,流体速度越高,侵蚀越强;地形整体起伏变化较大,流体能量耗散率越高;速度越小,更易堆积在沟谷中[78-79];另一大难点是数值算法的稳定性,解决这一问题通常可采用激波捕捉格式的有限体积格式对复杂地形的流动过程进行模拟[80-82]㊂泥石流的入汇是十分复杂的动力学过程,涉及固相颗粒的运移及水流对泥石流的多场耦合作用,在计算过程中需确定泥石流和水流的相界面㊂采用常规的数值计算方法难以解决这一问题㊂基于上述考量,可从流体力学理论出发结合非牛顿流体的本构关系得出泥石流与水流交汇区的统一方程㊂进一步考虑交汇区流体存在从牛顿体向非牛顿体转变的过程,结合前人研究成果,采用经验关系式求出浆体屈服应力,最后通过标志网格法(MA C)求解交汇区的控制方程,得到泥石流堵江的动态过程[93]㊂泥石流的堆积过程通常伴随着流体的体积质量减少和能量损耗[46]㊂在控制方程组中通常作为源项进行描述,源项为正表示侵蚀作用,源项为负表示堆积作用㊂表示堆积作用的源项通常以阻力为主的形式出现在方程右侧[24,26,46,94-95];也可以采用由现场资料得出的经验公式,如假定泥石流体积与沟谷横截面积和平面堆积面积存在关系,结合非定常明渠流理论得到的L A H A R Z模型[96]㊂对于颗粒流的堆积过程,可采用离散单元法对泥石流的堆积过程进行研究,得到有源和无源条件下的固体颗粒运移过程[14]㊂这一阶段主要将初期简化的流动ң堆积过程进行更为深入的考虑,将泥石流的流动过程细化为启动ң沿程运动(侵蚀放大㊁颗粒分选)ң汇流 固体物质运移㊁堆积㊂对固相颗粒动力学机理的研究中引入离散介质模型(D E M㊁L B M),对泥石流的物理建模更为充分且能够捕捉到颗粒在整个过程中的运移机制,为深入研究泥石流的动力学机理开拓了新的思路㊂5第5期杨昶等:我国泥石流数值模拟研究历程及发展趋势Copyright©博看网. All Rights Reserved.表2 侵蚀率公式的总结及说明[83]来源侵蚀率公式备注文献[42]E =α1c e -c c *-c e h 1u 1d ,i f (c <c e )α2c e -c c *h 1u 1d,e l s e ìîíïïïï(8)基于流体的泥沙浓度(c )和平衡浓度(c e )的对比得到文献[84]E =c *u 1t a n (θ-θe )(9)基于平衡状态下的坡度(θe )和实际环境坡度(θ)的对比得到文献[85]E =αu 21+v 21,i f (τ1z x b o t ȡτt h r e s h o l d )0,e l s e{(10)考虑侵蚀与流速线性相关的经验公式文献[66,86-89]E =φθ-θc h 1u 1d -0.2,i f (θȡθc )0,e l s e {(11)考虑溃坝洪水对河床侵蚀响应的经验公式文献[71]∂z b∂t=E M h 1u 21+v 21E M =l n (V f /V 0)/S(12)基于侵蚀前后体积变化求得平均侵蚀率并采用类似公式(11)思想得到的经验公式文献[90]E =σ2t o p -σ1t o pρw 1(z b ),i f (σ2t o p >σ1t o p )0,e l s e{(13)基于守恒方程推导得出的描述雪崩对基底侵蚀的公式文献[83]-∂z b ∂t =τ1b o t -τ2t o p ρu 1-u 1σ2t o p -σ1b o t τ1b o t -τ2t o p (14)不考虑河床地形影响的侵蚀率通用形式文献[91]-∂z b ∂t =τ1b o t -τ2t o pρu 1τ1b o t =ρs C f u 1u 1τ2t o p =t a n (φ2)(σz -p )u 1/u 1(15)不考虑河床地形及剪胀影响的侵蚀率通用形式文献[92]-∂z b ∂t =τ1b o t -τ2t o pρu 21+v 21τ1b o t 可用B i n g h a m ,H e r s c h e l -B u l k e y ,或V o e l l m y 模型τ2t o p =c +h 1ρg (1-λ)c o s θt a n φ2(16)公式(16)的变式文献[73]-∂z b ∂t =τ1b o t -τ2t o pρu 21+v 21τ1b o t =m a x [(ρg h t a n φ+(u 2+v 2)/C 2z ),(ρ(1-s )g h t a n δ)]τ2t o p =c +h 1ρg (1-λ)c o s θt a n φ2(17)公式(16)的变式 注:α1和α2为经验系数;c *为静止河床的泥沙浓度;C f 为无量纲摩擦系数;d 为粒径;φ为经验侵蚀系数;V a 和V f 为侵蚀发生前后的体积;h 为层厚;u ㊁v ㊁w 分别为速度沿x ㊁y ㊁z 的分量;下标1,2分别表示移动的上层流体和静止的下层基底㊂1.3 成熟多元阶段泥石流的动力致灾过程具有空间尺度大㊁时间跨度长的特点,因而传统的计算流体力学(C F D )和离散元(D E M )技术难以解决这一问题㊂我国学者在此基础上,提出广义深度积分模型[97]㊂考虑流动在流向和垂向上的巨大差异,对垂向进行深度积分,进而将三维流动问题简化为伪三维流动问题,因而大大缩短计算时间,提高计算效率㊂由于泥石流固㊁液两相流的特点,使得单一的连续介质模型或离散介质模型都不能很好地体现其物理本质,故发展结合两者优势的混合模型㊂构建混合模型普遍采用的方法是对一定粒径范围内的细颗粒和流体描述为泥石流液相浆体部分,对超出临界粒径的固体颗粒描述为泥石流的固相部分,固相部分采用离散介质模型,液相部分采用连续介质模型来刻画[98-100]㊂混合模型适用于固 液两相作用形式复杂,固相颗粒间作用力无法忽略的泥石流体,通常的技术手段是对液相部分用N S 方程进行模化,固相部分采用D E M 进行描述,通过一定的耦合技术可以实现泥石流对拦挡物的实时动力过程监测,如冲击破坏[98,101]和坝体截留[99]等㊂然而,在C F D D E M 耦合中,流体与颗粒之间的相互作用基于一个粗糙的网格,该网格无法精确地描述单个颗粒周围的流动或微尺度的流动,因此,有学者采用L B M D E M 的6水土保持学报 第37卷Copyright ©博看网. All Rights Reserved.耦合形式,以更加精细地模拟流㊁固两相之间的动力作用过程[22,102]㊂求解混合模型的难点在于如何耦合流体与固体部分的计算,若考虑固体变形对流体的影响则称为双向耦合,反之则为单向耦合[103](图6)㊂流固耦合技术的难点在于界面信息的传递,已有的流固耦合技术通过协调网格法(c o n f o r m i n g -m e s h m e t h o d s )或浸入边界/域(i mm e r s e db o u n d a r /d o m a i n m e t h -o d)的方法来实现流㊁固之间的信息传递㊂然而,协调网格的流固耦合技术在流体 固体界面的数据处理困难及迭代收敛性问题上存在问题,目前更多采用浸入式方法[104]㊂ 注:附图修改自参考文献[103]㊂图6 流固耦合的2种计算流程与前2个阶段研究成果相比,该阶段整体上体现我国泥石流的数值研究已整体迈入世界先进行列㊂对泥石流的动力致灾机理及物理性质理解更为深入,采用的建模方式更具有针对性和科学性,数值计算方法也更为稳定㊁高效㊂通过3类模型(连续介质模型㊁离散介质模型㊁耦合模型)可以对各种性质的泥石流进行物理描述,进而采用合理的数值计算方法研究泥石流在不同尺度下的实际问题㊂2 讨论经过近30年的发展,我国泥石流的数值模拟研究取得系统性的研究成果,体现在从宏观到微观㊁从灾后反演到临灾预报㊂随着计算机技术的发展,后处理技术由原先的平面二维可视化演变为如今的时空四维可视化,进而结合G I S 平台为我国大型工程项目建设㊁山区城镇规划以及防灾减灾工作实施,提供强有力的科技支撑㊂针对不同的泥石流动力学问题,模型的选取和数值处理方式不尽相同㊂泥石流运动模型的适用性与泥石流剪切率有关,如I v e r s o n 的混合流理论可归结为准静态模型,其应力模型与剪切率无关;T a k a h a s h i 的水石流膨胀体模型可以认为是高剪切率状态下的研究;而中剪切率状态下泥石流应力 应变过程与剪切速率有关㊂数值方法上,连续介质模型适用于宏观流动过程的研究,离散介质模型适用于研究微观颗粒的动力作用和多孔渗流问题㊂通过耦合离散介质与连续介质并且与室内试验相结合,可研究多尺度下的泥石流动力学机理㊂通过数值模拟技术可以更好地帮助研究人员探究各尺度下的泥石流动力学机理㊂根据不同的研究问题,选择合适的模型与方法尤为重要㊂通过对已有文献进行汇总,给出不同泥石流数值模拟方法的优缺点及其适用性(表3)㊂表3 不同数值模拟方法的优缺点及其适用性汇总数值模型计算方法优点缺点适用性有限差分法(F DM )简单且易于构造高精度格式复杂网格处理困难液相作用为主的泥石流有限体积法(F VM )守恒性好,适于处理复杂网格精度不易提高连续介质模型有限元法(F EM )基于变分原理,守恒性好复杂方程不易处理粒子类方法(S PH ㊁MPM 等)无网格依赖性,易于处理复杂外形精度较低㊁计算量大E-L 方法(P I C ㊁MA C 等)易于捕捉界面流动计算量大㊁稳定性较差离散介质模型格子玻尔兹曼法(L BM )并行性好,边界条件处理方便复杂形状问题求解困难,方程推导过程复杂离散单元法(D EM )直接描述颗粒的相互作用,尺度跨度大(n m~m )计算量大㊁参数标定困难固相颗粒作用为主的泥石流混合介质模型C F D-D EM ㊁L BM-D EM ㊁S PH-D EM 等能够处理复杂的流㊁固相互作用问题界面信息交互困难固液两相共同作用的泥石流7第5期杨昶等:我国泥石流数值模拟研究历程及发展趋势Copyright ©博看网. All Rights Reserved.。
