DOI: 10.11817/j.issn.1672-7207.2016.08.024

基于颗粒离散元的抗滑桩土拱效应分析

韩同春^{1, 2}，邱子义^{1, 2}，豆红强³

(1. 浙江大学 滨海与城市岩土工程研究中心，浙江 杭州，310058；

2. 浙江大学 软弱土与环境土工教育部重点实验室，浙江 杭州，310058；

3. 福州大学 环境与资源学院，福建 福州，350116)

摘要：运用 YADE离散元软件对三维空间下的土拱效应进行研究。探讨在不同土层埋深、桩间距及土颗粒不均匀分布等影响因素下的土拱效应发展及变化。研究结果表明：浅层土中的土拱较之深层土中的土拱强度较差，破坏也有所提前；随着桩间距的增大，土拱效应的影响将会减弱，并且减弱的速度逐渐变慢；土颗粒粒径分布的不均匀性增强会使得土拱的形成提前，但是对加载过程中，桩上水平力分担比峰值的影响不大。

关键词：土拱效应；YADE离散元；三维模型；土层深度；桩间距变化；土体不均匀性

中图分类号：TU473             文献标志码：A         文章编号：1672-7207(2016)08-2715-08

Soil arching effect between anti-slide piles based on YADE discrete element method

HAN Tongchun^{1, 2}, QIU Ziyi^{1, 2}, DOU Hongqiang³

(1. Research Center of Coastal and Urban Geotechnical Engineering, Zhejiang University, Hangzhou 310058, China;

2. Key Laboratory of Soft Soils and Geoenvironmental Engineering, Ministry of Education,

Zhejiang University, Hangzhou 310058, China;

3. College of Environment and Resources, Fuzhou University, Fuzhou 350116, China)

Abstract: YADE discrete element method was used to analyze the soil arching effect between anti-slide piles. The depth of soil layers, space between anti-slide piles and variation of soil particles diameters were given into consideration. The results show that at the different depths of soil layers, the developing of soil arching is different. When the depth of soil layer is shallow, the strength of soil arching is lower and the soil arching is destroyed early. When the depth of soil layer is deep, the strength of soil arching is higher and the soil arching is destroyed later. The arch rise decreases with the increase of the depth of soil layers. Space increasing between anti-slide piles weakens soil arching effect and the speed of weakness of soil arching effect is slower and slower. The increase of range of soil particles diameters distribution lets soil arching take shape in advance, but has little influence on the maximum horizontal load sharing ratio of piles in the loading process.

Key words: soil arching effect; YADE discrete element method; three-dimensional numerical model; depth of soil layers; variation of space between anti-slide piles; variation of soil particles diameters

土拱效应是岩土工程中一种常见的现象，TERZAGHI^[1]在1943年就通过“活动门”试验证实了土拱效应的存在。近些年来，众多学者分别从数值模拟、理论推导和室内试验等角度对土拱效应进行了深入研究^[2-7]。李忠诚等^[8]利用三维有限元结合ITOS等^[9]提出的土压力理论对成拱效应中桩间距的影响进行了系统的分析；杨明等^[10]使用离心机对桩间的土拱效应进行了研究；向先超等^[11]利用二维的离散元也进行了部分的研究。上述研究主要是建立在连续介质的基础之上，事实上土体是一种非均质各向异性的散粒体，利用有限元连续介质模型的研究必然存在一定的误差。并且，除此之外，由土拱效应的触发机理可知，土拱效应的产生过程即为桩后土体及桩间土和桩后土的水平位移发展不协调及运动土体的“楔紧”过程，土体的压密发生在三维空间中，仅仅在二维平面内考虑土拱问题与实际不符。随着离散元在岩土工程中的广泛应用^[12-14]，该方法对于大变形和非连续介质模拟的适用性良好的优点越来越突出。为此，本文作者采用离散元软件对三维空间中抗滑桩土拱效应进行研究，探讨不同深度土层的成拱形式及桩间距、土颗粒粒径分布等因素对土拱效应的影响。

1  土拱效应形成机理简述

土拱的形成主要是土体具有颗粒间的黏结力和摩阻力^[2]。就抗滑桩而言，土拱效应产生的机理为桩后土体在土压力或者滑坡的作用下产生位移，而抗滑桩的存在致使桩后土位移的发展受阻。与此同时，土体为典型的非均质各向异性散粒体材料，表现为桩后土体位移的发展在土压力作用下并不均匀，进而使得土体间发生“楔紧”作用。在土体自身存在抗剪强度及桩侧的摩阻力作用下，桩后土体形成以相邻两桩为拱脚的土拱，且桩后的土压力p随着土拱传递到拱脚及其周围土体上(见图1，其中D₁为桩径，D₂为桩间距)。所以，土拱效应的实质是应力在土体和抗滑桩之间的转移过程，即从运动的土体转移到抗滑桩上，使桩土共同分担荷载。

图1  桩间土成拱分析示意图

Fig. 1  Sketch map of soil arching effect

2  YADE离散法基本理论简介

离散元法的主要思想是将整个研究的物体离散为一系列刚性小球，小球几何特性(直径、空间位置、排列等)由人为控制产生，并且可以赋予刚性小球一定的物理特性(重度、摩擦角等)或是化学特性。小球运动受经典力学理论控制，整个研究物体的宏观物理性质由小球的运动和相互作用描述。

就YADE离散元软件而言，它是一个内核基于C++语言编写，外部的控制代码为Python语言编写的开源代码非商业三维离散元软件，方便用户基于此软件进行二次开发和更加深入的研究，亦可以根据自己的需要调用软件中相应的模块来实现对工程问题的模拟。而且，在离散元法的计算效率上，接触面的发现算法是影响计算效率的重要因素。YADE离散元程序在接触面的发现算法中采用“边界盒法”，具有较高的效率^[15]。

参考该离散元软件的教程^[16],两刚性小球的接触刚度分为法向刚度K_N和切向刚度K_T。切向刚度K_T由法相刚度K_N乘上K_T/K_N得到。其中，K_T/K_N决定宏观上泊松比的取值，对于法向刚度K_N，用2个长度等于球体半径的弹簧链表示，结合图2得

               (1)

式中：E₁和E₂分别为两接触小球的弹性模量；和在本次模拟采用的接触特性中取=2r₁，=2r₂。

图2  两小球法向接触刚度

Fig. 2  Normal stiffness of contact between two spheres

对于应变的计算，YADE中将应变分为法向应变和切向应变。法向应变的计算如图3所示。其中，和分别为两小球初始球心坐标，定义：

              (2)

           (3)

图3  两小球接触初始状态图

Fig. 3  Initial contact state of two spheres

              (4)

法向位移u_N及法向应变ε_N的计算式如下：

           (5)

         (6)

式中：和分别为两小球发生位移后的球心坐标。

对于压缩大变形的情况，由式(5)和(6)可知：其应变将趋向于负无穷，但是小球间的力还保持为一个有限的值，会发生球体间彼此穿透的现象。所以，对法向应变的计算式进行如下调整：

      (7)

剪切应变的计算有全量算法和增量算法2种，增量算法的剪切位移为

       (8)

式中：为前一计算步的剪切位移；为投影引起的剪切位移增量；为旋转引起的剪切位移增量；为两小球相对运动引起的剪切位移增量。

对于应力的计算，在已知的法向和切向位移条件下，正应力F_N和切应力的计算式如下：

                (9)

                (10)

切应力仅仅是计算最终切应力的中间过程，最终切应力F_T的计算需要经过下式的计算获得：

      (11)

式中：φ为土体的内摩擦角。

3  数值模型

参照杨明等^[10-11]建立的二维离散元模型,建立三维离散元模型。根据对称性，取排桩中相邻两桩中心线土体为研究对象，所研究的土体限于桩后及两桩中间的土体。由于所研究的土体为某一深度处的土层，在上部土层的约束作用下，土体竖向的位移相对于水平向位移较小，并且，整个土层的竖向位移为零并不会导致土颗粒间的“楔紧”作用仅发生在水平方向，土拱效应的产生依然是一个三维的过程。所以，作近似处理，假定土体的位移仅仅发生在水平方向上，在矩形盒的顶部和底部施加竖直向的刚性约束。

包围土体的矩形盒表面光滑，前端开口。为保证土体因水平方向上产生的不均匀位移引发的“楔紧”作用不是人为施加，而是抗滑桩和土颗粒间及土颗粒和土颗粒间相互作用产生的结果，使用竖向加载墙对各土层施加以相同的初始水平位移。加载的速度为0.05 m/s，所施加的速度方向为沿矩形盒长度方向。矩形盒的长度为5 m，高度为3 m，宽度依据不同的桩间距调整。

抗滑桩为直径0.9 m的圆截面桩，桩土摩擦因数为0.5，忽略桩的弹性变形，视之为刚体且桩端完全固定，无水平方向的位移。同时忽略桩后土体的土压力作用以形成临空面。

采用刚性小球来模拟离散土体，若是按照实际的尺寸来模拟土颗粒，则所建立的模型中小球将多达数十万个，以目前计算机的计算速度无法满足要求。而且，本文所研究的是土拱效应在不同条件下的变化趋势，数值仅需比较相对大小，并不需与实际的工程做比较。所以，模拟中对土颗粒的尺寸进行适当的放大，土体的孔隙率为15%，小球直径为0.02 m。根据室内三轴试验和YADE离散元模拟实验之间的参数对比分析，确定土体的其他细观物理力学参数见表1。

表1  土体细观参数

Table 1  Meso-parameters of soil

依据以上边界条件及土体的各项物理指标，利用YADE开源代码离散元软件，在程序中定义符合上述物理指标的材料直接赋予所产生的刚性小球，并且，控制小球产生的数量来满足初始孔隙率的要求。建立数值模型如图4所示。

图4  离散元三维数值模型

Fig. 4  Numerical model of DEM^3D

4  土拱效应影响因素的数值分析结果

4.1  不同深度土层的成拱形式

取抗滑桩的净间距为2.7 m，程序中虚拟时间每过1 s保存1次加载过程中各离散元小球的物理状态。保存的数据利用paraview后处理软件进行处理，即可得到离散元小球接触面上正应力力传递图。通过观察离散元小球接触面上正应力力传递图的变化，来间接得知抗滑桩间土拱的形成、发展和破坏过程。所进行的模拟步为25 000步左右，所用时间约为1.5 h。

对于整个模型所模拟的土体，在加载墙不同的推进位移情况下，如图4所建立坐标系，取不同h处的正应力传递图的水平剖面，可得结果如图5所示的土颗粒小球接触面间正应力传递图(d为加载墙位移)。从图5可以看出：处于不同埋置深度h的土层，所成的土拱形状也不相同。随着土层深度的增加，土拱的矢高f逐渐降低。对于同一h处的土层水平剖面，当h= 0.5 m时，随着加载墙上加载位移的发展，土拱在早期就可形成，并且在加载位移d=5~15 cm时，土拱的形状变得更为清晰，在d=40 cm时，两相邻抗滑桩间的应力分布图已出现不连续的现象，土拱发生破坏；当h=1.5 m时，土拱在d=5 cm时，土拱已然形成，在d= 15 cm时，土拱形状已经很完整，在d=15~40 cm的过程中，土拱已经开始产生破坏，但是破坏的程度不如h=0.5 m时明显；当h=3.0 m时，在整个加载过程中，土拱形状始终完整。

因此，可以得知：随着土体的埋置深度的增加，土拱的矢高高度降低，土拱的形状变得扁平，其破坏的过程亦会发生滞后。此现象与二维平面中单一土层成拱和破坏的过程有所不同。

考虑到三维空间与二维平面在几何方面的不同点在于三维空间中需要考虑深度方面的影响，并且，在本次模拟中，沿深度h方向施加的外力为重力，所以，土拱形状沿深度的变化及破坏的滞后性均与重力有关。

若是忽略重力的影响，其余条件同前，则可以得出如图6所示的整个土体正应力分布俯视图。

结合图6和图5可以得知：重力对于不同h处土层的成拱形式有着很大的影响。在无重力条件下，各层土体的成拱和土拱的破坏均按照同一种形式完成，并且土拱的破坏也有所提前，在d=30 cm时就发生了土拱的破坏；然而在重力的影响下，各层土体的成拱和土拱的破坏过程各不相同，下层土的破坏较之于上层土表现出一种滞后性。

进一步探究重力对于不同h处土层成拱形式的影响原因，图7所示为在d=15 cm时，常重力场和无重力场条件下的加载墙上水平推力分布竖直剖面图。由图7可知：加载墙上的推力分布沿深度增加。所以，重力的存在使得不同h处产生土拱时，所受到的推力不同，近似可以认为沿土层深度成梯形分布。研究表明^[17]：随着桩后土压力的增加，土拱矢高f呈陡降形对数函数的变化规律。所以，在土层推力因重力随深度增加时，土拱的矢高f会减小。并且，根据Mohr-Coulomb破坏准则，土体间的有效正应力增加可以提高土体的抗剪强度。重力的存在使得土体间有效正应力增加，提高了抗剪强度。所以，深层土的成拱发生时间较早，土拱也较为稳定。

综合上述分析可知：不同深度的土层在重力场的影响下，所受到的加载墙水平推力不同，从而土拱的形状表现出不同的形式。而且，重力的存在有利于土拱的形成。

4.2  桩间距对土拱效应的影响

保持其他的控制条件不变，改变桩间的间距，将桩径记为D₁，桩间净间距记为D₂，桩上水平力记为F，加载墙上水平推力记为p，定义桩上水平力分担比为F/p.分别对D₂/D₁为1，2，3，5和7的情况进行模拟计算。图8所示为在不同D₂/D₁下，各加载过程中，F/p的最大值。

图5  不同深度土层土颗粒间正应力传递图

Fig. 5  Pathways of normal stress of soil layers with different depths

图6  无重力条件下整个土层土颗粒间正应力传递图

Fig. 6  Pathway of normal stress for whole soil without gravity acceleration

图7  加载墙水平推力分布

Fig. 7  Distributions of horizontal force on load wall

图8  不同D₂/D₁下，桩体水平荷载最大分担比变化

Fig. 8  Variation of maximum load sharing ratio of piles with different D₂/D₁

从图8可以看出：随着D₂/D₁的增大，桩对于加载墙上水平推力分担比减小；并且，D₂/D₁较小时，分担比的衰减速度较快，而D₂/D₁较大时，分担比的衰减速度较慢。

在加载墙位移d=10 cm和土层深度h=1.5 m的条件下，进一步考察不同D₂/D₁条件下，抗滑桩间离散元小球接触面上正应力传递图的应力传递变化，结果如图9所示。

从图9可以看出：当D₂/D₁为1，2和3时，抗滑桩间的土拱明显，两桩间的应力线可以连为一体，土拱的整体性表现明显；但是当D₂/D₁为5时，两桩周围的应力线出现连接稀疏，土拱的整体性也不是十分明显；当D₂/D₁为8时，两桩周围的应力线已经基本没有连接，土拱也没有出现完整的形状。

结合图8与图9可知：抗滑桩间土拱是否能完整形成与桩间距有很大关系。当桩间距较小时，桩后土体会明显受到桩的约束作用，土体与桩之间的相互作用将促使土拱效应的发生，土体所受的力将经过土拱传递给抗滑桩，表现为抗滑桩上水平推力的集中。成拱效应越显著，水平推力的集中度越高。若是桩间距较大时，桩后土体受到桩的约束将大幅减弱。在加载墙水平推力的作用下，桩后土体将从两抗滑桩中间滑出，土体上所受到的力有相当大的一部分不会传递到桩上，仅仅是靠近桩的部分土体与桩之间存在力的传递。另外，在桩间距较小时，桩上水平力的集中程度主要由成拱效应的显著程度来控制，所以，当桩间距增加时，成拱效应的显著程度将下降，桩上水平力的集中程度也将明显下降；在桩间距较大时，桩上水平力的集中程度实际上由靠近桩的部分土体与桩之间相互作用控制，受桩间距的变化的影响较小，所以表现出桩上水平力分担百分比的衰减速度减慢的现象。

4.3  土颗粒粒径分布的影响

土体由于形成原因的不同，土颗粒的粒径也会有所不同，即使在同一地理位置的土体，其土颗粒粒径也会有所差异。YADE离散元程序可以对通过调整生成刚性小球的直径，来模拟土体中土颗粒粒径的差异。取桩间距为2.7 m，其余条件均与初始相同，保持生成的刚性小球平均直径为0.02 m，调整生成小球时直径的方差分别为0，0.3，0.6和0.8，以此来达到模拟土体中土颗粒粒径不同的目的。所生成的小球颗粒级配曲线如图10所示。

图9  不同D₂/D₁条件下，应力传递图

Fig. 9  Pathways of normal stress for different D₂/D₁

图10  不同粒径分布下级配曲线

Fig. 10  Grading curves of different particle size distributions

图11所示为在不同方差下，整个加载过程中，桩上水平力占加载墙上水平推力的变化。由图11可以看出：随着土体不均匀性的加强，桩上水平力占加载墙上推力比值的峰值有所提前，在产生小球方差为0时，在加载墙位移d=15 cm的情况下，桩上水平力与加载墙上水平力比值F/p达到峰值，在产生小球方差为0.3，当d=10 cm时，F/p即可十分接近峰值。在产生小球方差为0.6和0.8，当d=5 cm时，F/p即为峰值。这表明土体的不均匀，使得土拱在加载的早期形成，并且不均匀性越强，土拱形成的时间越早。这是因为随着土体中土颗粒的不均匀性增强，在土拱形成的过程中，土颗粒间彼此的相对位置调整加剧，完成“楔紧”过程所需要的时间缩短，使得土拱的形成时间提前。

图11  不同粒径分布下桩体水平荷载分担比

Fig. 11  Horizontal load sharing ratios of piles with different particle size distributions

若是从峰值方面来看，土体不均匀性较大时，桩上水平力占加载墙上水平推力比值的峰值略有上升，土体不均匀性并没有对桩上水平力占加载墙上水平推力比值的峰值产生较大的影响。

5  结论

1) 三维情况下土拱的形成和矢高与土层的深度有关，浅层土土拱强度较差，破坏早于埋深较深的土层形成的土拱；土拱的矢高随着埋置深度的增加而减小，原因是重力场的存在使得各土层的的水平推力不同，从而对土拱矢高产生影响。所以，在实际的抗滑桩设计中，应当明白桩后土体水平推力沿深度的分布形式，结合土拱矢高及强度随桩后土体水平推力变化的规律综合考虑来确定合理的桩间距。

2) 抗滑桩桩间距对土拱的形成有很大的影响，桩间距增大，土拱效应将会减弱，并且，成拱效应随桩间距增加时的衰减程度并不相同，在桩间距较小时衰减较为明显。

3) 土体分布的不均匀性对于土拱效应的发展有较明显的影响，土体的不均匀性越大，土拱效应的形成会越早发生。但是，土体的不均匀性对于桩上水平力占加载墙上水平推力比值的峰值影响不大，仅仅是略有上升。

参考文献：

[1] TERZAGHI K. Theoretical soil mechanics[M]. New York: John Wiley and Sons, 1943: 66-76.

[2] LIANG R, ZENG S. Numerical study of soil arching mechanism in drilled shafts for slope stabilization[J]. Soil and Foundation, 2002, 42(2): 83-92.

[3] BOSSCHCR J, CRAY H. Soil arching in sandy slopes[J]. Journal of Geotechnical Engineering, 1986, 112(6): 626-645.

[4] 张建华, 谢强, 张照秀. 抗滑桩结构的土拱效应及其数值模拟[J]. 岩石力学与工程学报, 2004, 23(4): 699-703.

ZHANG Jianhua, XIE Qiang, ZHANG Zhaoxiu. Arching effect of anti-slide pile structure and its numerical simulation[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(4): 699-703.

[5] 周德培, 肖世国, 夏雄. 边坡工程中抗滑桩合理桩间距的探讨[J]. 岩土工程学报, 2004, 26(1): 132-135.

ZHOU Depei, XIAO Shiguo, XIA Xiong. Discussion on rational spacing between adjacent anti-slide piles in some cutting slope projects[J]. Chinese Journal of Geotechnical Engineering, 2004, 26(1): 132-135.

[6] 刘钦, 李地元, 刘志祥, 等. 水平推力作用下抗滑桩间土拱效应影响因素的数值分析[J]. 中南大学学报(自然科学版), 2011, 42(7): 2071-2077.

LIU Qin, LI Diyuan, LIU Zhixiang, Numerical analysis of influence factors on soil arching effect between anti-sliding piles under horizontal pushing loads[J]. Journal of Central South University (Science and Technology), 2011, 42(7): 2071-2077.

[7] 杨雪强, 吉小明, 张新涛. 抗滑桩桩间土拱效应及其土拱模式分析[J]. 中国公路学报, 2014, 27(1): 30-37.

YANG Xueqiang, JI Xiaoming, ZHANG Xintao, Analysis of arching effect between anti-slide piles and different arch body modes[J]. Chinese Journal of Highway and Transport, 2014, 27(1): 30-37.

[8] 李忠诚, 梁志荣. 侧移土体成拱效应及被动桩计算模式分析[J]. 岩土工程学报, 2011, 33(Suppl 1): 106-111.

LI Zhongcheng, LIANG Zhirong. Soil arching and calculation model for soil pressure of passive piles[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(Suppl 1): 106-111.

[9] ITO T, MATSUI T. Methods to estimate lateral force acting on stabilizing piles[J]. Soil and Foundations, 1975, 14(4): 43-59.

[10] 杨明, 姚令侃, 王广军. 桩间土拱效应离心模型试验及数值模拟研究[J]. 岩土力学, 2008, 29(3): 817-822.

YANG Ming, YAO Lingkai, WANG Guangjun. Study of centrifuge model tests and numerical simulation on soil arching in space of piles[J]. Rock and Soil Mechanics, 2008, 29(3): 817-822.

[11] 向先超, 张华, 蒋国盛, 等. 基于颗粒流的抗滑桩土拱效应研究[J]. 岩土工程学报, 2011, 33(3): 386-391.

XIANG Xianchao, ZHANG Hua, JIANG Guosheng, et al. Soil arching effect of anti-slide piles based on particle flow method[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(3): 386-391.

[12] 周健, 亓宾, 曾庆有. 被动侧向受荷桩模型试验及颗粒流数值模拟研究[J]. 岩土工程学报, 2007, 29(10): 1449-1454.

ZHOU Jian, QI Bin, ZENG Qingyou, Model tests and PFC^2D numerical analysis on laterally loaded passive piles[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(10): 1449-1454.

[13] 韩高孝, 宫全美, 周顺华. 摩擦型岩土材料土拱效应微观机制颗粒流模拟分析[J]. 岩土力学, 2013, 34(6): 1791-1798.

HAN Gaoxiao, GONG Quanmei, ZHOU Shunhua. Analysis of microcosmic mechanism of soil arching in frictional geotechnical material by particle flow simulating[J]. Rock and Soil Mechanics, 2013, 34(6): 1791-1798.

[14] 许朝阳, 周健, 完绍金. 桩承式路堤承载特性的颗粒流模拟[J]. 岩土力学, 2013, 34(1): 501-507.

XU Chaoyang, ZHOU Jian, WAN Shaojin. Simulation of bearing characteristics of pile-supported embankments by particle fiow code[J]. Rock and Soil Mechanics, 2013, 34(1): 501-507.

[15] 刘红, 刘军. 离散元法中三种球形颗粒接触发现算法的比较[J]. 黑龙江科技学院学报, 2007, 16(6): 360-363.

LIU Hong, LIU Jun. Arithmetic compare of three contact detection methods for spherical particles in DEM[J]. Journal of Heilongjiang Institute of Science & Technology, 2007, 16(6): 360-363.

[16] V, CATALANO E, CHAREYRE B, et al. Yade Documentation[EB/OL].https://yade-dem.org/doc/, 2013-10-17.

[17] 吕庆, 孙红月, 尚岳全. 抗滑桩桩后土拱形状及影响因素[J]. 哈尔滨工业大学学报, 2010, 42(4): 629-633.

L Qing, SUN Hongyue, SHANG Yuequan. Shape of soil arch behind anti-slide piles and its major influence factor[J]. Journal of Harbin Institute of Technology, 2010, 42(4): 629-633.

(编辑  杨幼平)

收稿日期：2015-08-17；修回日期：2015-10-11

基金项目(Foundation item)：国家自然科学基金资助项目(51178423)；浙江省自然科学基金资助项目(LY15E080010)(Project (51178423) supported by the National Natural Science Foundation of China; Project(LY15E080010) supported by the Natural Science Foundation of Zhejiang Province)

通信作者：韩同春，博士，副教授，从事边坡稳定、地基处理和土工数值分析等研究；E-mail：htc@zju.edu.cn