高心墙堆石坝心墙水力劈裂的颗粒流模拟

杨 艳 ，周 伟，常晓林，花俊杰

（1. 武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2. 西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100）

1 引 言

在众多坝型中，土石坝因其具有广泛的适应性、可就地取材、施工方法简单且抗震性能良好等优点而成为高坝的首选坝型，在坝工建设中得到了极其广泛的应用。特别是随着我国西南水电开发进程的加快，在建和拟建的200～300 m 级高土石坝的数量在逐渐增加。但在高土石坝的建设中，心墙堆石坝作为土石坝中常用的坝型之一，还存在着一些急需解决的重要问题，其中心墙的水力劈裂问题是土石坝设计和建设中倍受关注的焦点之一。

目前，已有很多学者针对心墙的水力劈裂问题进行了研究，并取得了不少成果。黄文熙[1]指出，心墙中任何一点处的孔隙水压力如果使该点处的最小主应力的有效值降低至心墙土料的抗拉强度，心墙就会沿着这个最小主应力面产生水力劈裂。殷宗泽等[2]分析了常用的有效应力法和总应力法计算水力劈裂时结果产生较大差异的原因，认为应该用心墙外水压力是否超过心墙上游面处土中的中主应力来判别水力劈裂发生的可能。沈珠江等[3]模拟了Teton 坝水力劈裂的离心模型试验，并得出结论：深截水槽心墙未必会发生水力劈裂破坏，常用的总应力法过大地估计了水力劈裂发生的可能性。张丙印等[4]研制了一种新型的水力劈裂试验装置，采用糯扎渡高心墙堆石坝心墙混合土料进行了水力劈裂试验，结果证实了土石坝心墙中可能存在的渗水弱面以及水库在快速蓄水的过程中产生的弱面水压楔劈效应是导致水力劈裂发生的重要条件。李全明 等[5]采用弥散裂缝理论描述水力劈裂裂缝的发展过程，建立了用于描述水力劈裂发生和扩展过程的数学模型及有限元计算模式。朱俊高等[6]研究认为，水库蓄水初期是水力劈裂的危险期；完全均质的心墙内不会发生水力劈裂；裂缝或局部的缺陷及迅速蓄水的初期是土石坝心墙发生水力劈裂的两个重要条件，水力劈裂发生的根本原因是局部高水力梯度的存在。陈五一等[7]指出，目前水力劈裂分析常用的有效应力法存在不足，并提出改进方法，定义了一种判定水力劈裂的安全系数。曹雪山等[8]提出了研究心墙水力劈裂问题的非饱和土固结简化计算的有效应力分析方法，研究发现，提高心墙的渗透系数和心墙填筑土的初始饱和度及在初次蓄水时放慢蓄水速度等均可防止心墙水力劈裂的发生。毕庆 涛[9]提出了用非饱和心墙料做固结不排水试验得出总应力法计算模型参数的改进方法以及一种新的水力劈裂判断标准，用心墙内紧靠上游表面单元的组合应力(32σ -1σ )与心墙前库水压力进行比较来判定水力劈裂的发生。

综上所述，可以看出现有的成果大部分是通过试验和有限元数值模拟等手段，从宏观的角度研究水力劈裂的发生、发展过程及其力学机制，至今尚未形成一致的观点，仍存在许多争议。正是由于从宏观的角度难以完全揭示水力劈裂的破坏机制，因此，本文在考虑流-固耦合的前提下，采用PFC2D软件的颗粒流方法，从细观的角度初步探究了心墙水力劈裂发生的机制，并将该成果与现有的室内试验成果[4]进行比较分析，证实了该方法的可行性，对今后进一步研究心墙的水力劈裂具有重要的参考价值。

2 颗粒流流-固耦合的基本原理

颗粒流流-固耦合的基本原理[10]基于2 个假设：①假想固体中流体的渗流路径是由颗粒间接触处的平行板通道组成，该通道称为“管道”（如图1 中黑色线段所示，图中灰色圆形表示颗粒，白色线段表示颗粒间的接触），其管径的大小与接触处颗粒间的法向距离成正比；②本文采用的PFC2D颗粒流程序没有模拟流体的存在，仅仅通过流体压力来体现流体的作用，于是假设计算模型中存在可以存储压力的单元，该单元在PFC 中用域（如图1 中由白色线条围成的一个闭合多边形区域，黑色圆点表示域的中心）来定义，相邻域之间通过管道连接，于是可以根据域间的压力差实现流体流动的模拟，而流体的流动将会引起域内压力的改变。

图1 管道和域的示意图 Fig.1 Sketch of pipes and domains

PFC2D颗粒流程序在模拟固体和流体的耦合作用时，是通过下面的流体计算公式和压力计算公式实现的。

2.1 流体计算公式

在PFC 中用管道概化流体的渗流路径，该管道相当于一个平行板通道，其长度为L，管道孔径为a，在垂直于平面的方向上取单位厚度；另外，域将流体计算区域离散后，计算域内的“流体”是通过相邻域之间的压力差产生流动。于是，由平行板均匀流的立方定律可以得到管道内的流量公式为

式中：k 为渗透系数；a 为管道的孔径；P2-P1为两个相邻域的压力差；L 为管道的长度。

对于具有一定渗透性的无裂缝材料而言，其PFC 模型中的渗流通道用管道概化后，a 的大小将影响模型的渗透性。假设当颗粒间的接触力为0 时对应的管道孔径为a0，并称a0为残余孔径。当给模型施加实际的应力边界条件后，接触处法向可能的受力情况分别是受压、受拉和不受力。

当法向接触力为压力时，a 将随着接触处法向压力的增大而逐渐减小，此时a 与a0间的经验关系如下：

式中：F0为a 减小到a0/2 时的法向压力；F 为荷载作用下的法向接触压力。

当法向接触力为拉力或者0 时：

式中：g 为两颗粒表面间的法向距离；m 为两颗粒表面间法向距离的缩放因子，其作用是调节a 的大小。

当接触力为0 时，两接触颗粒表面间的法向距离g 为0，从式（3）中可以得出此时颗粒间接触处的a =a0，这与前面的假设是一致的。

2.2 压力计算公式

在一个时间步 tΔ 内，域中流体压力的变化可以通过下式计算：

式中：Kf为流体的体积模数；Vd为域的表观体积；∑q 为每个域从周围管道中获得的总流量，以流入为正；dVΔ 是考虑力的作用引起域的体积改变。

2.3 流-固耦合方式

在PFC2D中，流-固耦合的实现可以采用以下3种方式：（1）通过接触的张开与闭合或者接触力的改变实现管道孔径的变化；（2）作用在颗粒上的力改变了域的体积，从而引起域内压力的变化；（3）域内的压力差导致颗粒上作用着渗透体积力。

以上3 种形式中，前两点已经通过式（2）～（4）实现。针对第3 点，为了简化施加在颗粒上的力，假设域内压力施加在周围颗粒上的作用力沿其周围接触颗粒间的连线均匀分布（如图2 所示，图中黑色圆点表示域的中心，白色线段表示相邻接触颗粒间的连线，黑色箭头表示均匀分布的水压力P）。

图2 域内压力分布 Fig.2 Distribution of pressure in domain

由此可知，作用在颗粒上的力矢量可以表示为

式中：ni为相邻接触颗粒间连线的法向单位矢量；s为颗粒圆心到接触点的距离。

2.4 求解方法

PFC2D进行流-固耦合计算时采用显示求解方法，对每个管道和域分别进行流量和压力计算，并且在整个求解过程中两者的计算是交替循环进行的。

然而，为了保证渗流计算的稳定性，渗流计算时步不应大于临界时步，下面推导临界时步的计算公式。假设某个域存在扰动压力 pPΔ ，由于扰动流入域内的流量可以根据式（1）计算得出：

式中：R 为该域周围颗粒的平均半径；N 为连通该域的管道个数。

由式（4）可以计算出该流量在域内产生的压力变化 rPΔ 为

要保证渗流计算稳定，那么扰动流量引起的压力变化 rPΔ 必须小于扰动压力 pPΔ ，当两者相等时，由式（6）、（7）可求出临界时间步长为

在计算中，整体的时间步长必须取所有局部时间步长中的最小值，同时需乘以一个小于1 的安全系数。

3 颗粒流水力劈裂模拟

由于土石或堆石坝的心墙近似满足平面应变的特性，在模拟时可以按照平面应变问题考虑，因此，本文采用PFC2D进行颗粒流水力劈裂的二维模拟。

3.1 颗粒流模型的建立

本文参考张丙印等[4]的室内试验数据资料，将本次颗粒流模型的尺寸取为20 cm×20 cm，颗粒半径的取值范围为2.0～3.0 mm，半径取值的分布形式采用均匀分布，孔隙率为0.15。根据以上参数，本次模拟一共生成颗粒1 665 个（见图3，图中灰色圆形表示颗粒，四周的黑色线段表示墙）。 另外，考虑到心墙土为黏性土，本次模拟在颗粒间设置了接触黏结。但不足之处是，目前PFC 中的细观和宏观参数之间还没有建立起对应的关系，然而黏结强度（包括法向和切向）的取值将直接影响数值试验中水力劈裂的起始压力，因此，在模拟过程中，为了使得数值计算结果符合室内试验的规律，经多次调试，最终选取的法向和切向黏结强度n_bond 和s_bond 均为1.05 kN。当接触力超过黏结强度时，相应的接触黏结发生破坏，从而可以模拟裂缝的形成。

图3 颗粒流模型 Fig.3 Particle flow model

为了保证计算模型的边界条件符合心墙坝上游面实际的应力条件，本文采用PFC2D中的伺服控制机制控制边界墙的速度，从而使得计算模型达到给定的应力边界条件，如图3 所示。

3.2 心墙水力劈裂发生的必备条件

目前很多学者对于产生水力劈裂的条件有着一致的看法，认为产生水力劈裂必须具备以下条件[11]：①心墙中存在初始裂缝或缺陷；②心墙材料具有较低的渗透性；③能够产生楔劈效应；④水库快速蓄水。本文为了能够较真实地模拟水力劈裂的发生和发展过程，尽量满足上述的4 个必备条件。

本文生成模型时，在上游侧设置了长为0.08 m，厚为0.005 m 的水平软弱带，并通过在模型的上游面和水平软弱带区域内设置相应的库水压力来实现水库的快速蓄水过程，如图4 中深灰色区域所示（为了更加清晰地显示上游水压力和软弱带的布置，图中只标出了颗粒圆心的位置并未显示颗粒，浅灰色线条表示颗粒间的接触，圆圈1 表示测量圆1，本文中设置该测量圆是用于测量模型内部的应力值）。如果在初设的上游库水压力作用下未发生水力劈裂，则逐渐增加库水压力，直到水力劈裂发生。本文在模拟过程中设置水平初始缺陷的主要原因是，与竖向初始缺陷相比土石或堆石坝在逐层碾压的施工过程中较易产生水平缺陷。

图4 数值模型 Fig.4 Numerical model

3.3 计算参数的选取

采用颗粒流模拟土的实际宏观力学行为，如心墙的水力劈裂，在模拟中能否实现细观参数与宏观力学参数的统一将直接影响结果的正确与否。然而，这是一个十分复杂的过程，目前还没有统一的方法。本文为了减少在选取参数时的盲目性，首先对流-固耦合中的流体计算参数逐一进行敏感性分析，结果见图5。

图5 中的5 组曲线图均显示的是计算模型中平均应力(1σ +2σ )/2 随时步的变化过程，此处的1σ 和2σ 是指通过模型内部的测量圆1 得到的主应力。由于在水力劈裂研究中十分关注裂缝的继续扩展，因此，本文只跟踪记录了水平软弱带之后模型内部平均应力的变化过程，将测量圆1 布置在如图4 所示的位置。而本文之所以将心墙中的平均应力作为衡量指标，主要是考虑到目前很多学者认为，当心墙前的外水压力值大于心墙内的某一主应力（1σ 、2σ或3σ ）时，将导致水力劈裂的发生。由此可见，心墙内的应力值是判断水力劈裂发生的重要数据之一。另外，图5 中只显示了5 000 步之后平均应力的变化过程，这是因为在此之前是模型生成过程以及在给定的应力边界条件下达到平衡所经历的时步，而本文研究的重点是快速蓄水之后模型内部的应力变化及裂缝产生和发展的过程，因此，图5 仅从施加水荷载之后开始记录平均应力的变化过程。

图5 平均应力-时步关系曲线 Fig.5 Curves of average stress-step

从图5(a)、(b)可以看出，残余孔径a0和渗透系数k 取不同值对试样平均应力有较大影响。当a0值在0.010 0～0.000 5 mm 之间变化时，试样平均应力的最大差值达到0.16 MPa，是a0=0.010 0 mm 时的1.3 倍；当k 值在1.0×10-6～1.0×10-9mm/s 之间变化时，试样平均应力的最大差值达到0.13 MPa，是k = 1.0×10-6mm/s 时平均应力值的0.79 倍。这主要是由于上述2 个参数直接反映了心墙的渗透性，而心墙发生水力劈裂的实质是在高水力梯度的作用下，裂缝形成并发展的过程。因此，在模拟过程中应选取适当的渗透参数，确保心墙材料的低渗透性，从而模拟高水力梯度的作用。而a =a0/2 时的法向压力F0、流体体积模数kf以及颗粒表面间法向距离的缩放因子m 分别取不同值时，试样平均应力的最大差值仅为0.04 MPa，由此可见，上述3 个参数在合理地取值范围内变化时，对平均应力的结果影响较小。因此，参数选取时应重点调整a0和k 的取值。

此外，图5 中5 组曲线的变化规律基本一致，大约在5 500 步之前，压应力随时间步的增加而逐渐减小，这是因为在水压力作用下，初始的水平软弱带内形成了水平裂缝，在裂缝尖端附近区域的接触黏结受拉（见图6，图中黑色线段表示颗粒间的接触黏结承受拉力，线段越粗表示其对应的拉力值越大，灰色线段表示颗粒间的接触），这部分拉力抵消了一部分竖向初始压力，这种效应被称为水压楔劈效应[3]。而随后压应力又随时间步的增长而逐渐增加，其原因是当接触黏结破坏导致颗粒间产生裂缝后，这部分拉力被释放，作用在颗粒上的水压力使得颗粒向裂缝两侧移动，增加了裂缝两侧颗粒间的接触压力。但是当a0和k 取值较大时，由于模型的渗透性较大未能形成高水力梯度的作用，因此，没有发生水力劈裂，之前产生的拉应力也就无法释放。

图6 裂缝尖端区接触拉力分布 Fig.6 Distribution of contact tension around crack tip

本文在敏感性分析的基础之上，经过反复调整参数取值，最终选取的流体计算参数值见表1。

表1 流体计算参数 Table 1 Computational parameters of fluid

3.4 数值模拟结果及分析

本文在模拟心墙水力劈裂过程时，首先取一初始水压力作用在模型的上游表面和水平弱面内，如果在给定的水压力值下运行2 000 步后未能在模型内部产生新裂缝（本文以接触黏结发生破坏作为相应接触处产生裂缝的判定标准），则在计算过程中逐渐增加水压力值，直至在其内部产生连通上下游的贯穿性裂缝（如图7 所示），并在数值模拟中以此作为判别水力劈裂发生的依据，此时对应的库水压力值即为劈裂水压力Pf。

图7 PFC2D 模拟结果 Fig.7 Simulation result of PFC2D

本文针对表2 中1σ 、3σ 的3 种组合方案分别进行了颗粒流模拟，并将PFC2D的模拟结果与室内试验成果[3]做了比较（见表2 和图8）。从图可以看出，由数值模拟得到的结果，其规律与室内试验的基本一致，即劈裂水压力Pf随着竖向应力的增大而增大，且两者基本成线性关系。但由于在PFC2D的数值模拟过程中还存在着许多不确定的因素，如计算结果受颗粒的大小以及颗粒的组装形式等因素的影响较大。另外，在数值模拟过程中还做了一些简化和假定，从而导致计算结果均较室内成果偏大。

当竖向应力为0.2 MPa 时，本文采用颗粒流数值模拟得到的劈裂水压力为0.21 MPa，比室内试验的结果大0.03 MPa，偏大部分占室内结果的16.7%；当竖向应力采用0.3 MPa 时，本文得到的劈裂水压力为0.29 MPa，比室内试验的结果大0.02 MPa，偏大部分占室内结果的6.9%；在第3 种方案中竖向应力取为0.4 MPa，此时本文计算得到的劈裂水压力为0.38 MPa，与室内试验的结果相比大0.03 MPa，偏大部分占室内结果的8.6%。经过上述的比较可以看出，采用PFC2D进行水力劈裂的模拟，其结果虽与室内试验的成果有一定出入，但偏大的范围仅在6.9%～16.7%之间，这说明采用颗粒流方法研究心墙的水力劈裂是一种可行的方法。

表2 PFC2D 数值计算结果 Table 2 Numerical results of PFC2D

图8 数值结果与室内试验成果 Fig.8 Numerical and laboratory results

产生水力劈裂后的模型见图7，从该图中可以清晰地看到模型在高水力梯度作用下，由初始水平软弱带向下游发展形成的水力劈裂面，见图中近似水平的白色连通区域（其内部用黑色的线段表示颗粒间产生的裂缝）。此外，图中还显示了模型内部速度场的分布，从局部放大图中可以看出，裂缝周围颗粒的速度方向（黑色箭头所示）基本垂直于裂缝的延伸方向，从而初步证明了心墙水力劈裂破坏应属于颗粒间的张拉破坏，其破坏的主要力学原因是由于心墙中的张拉应力超过了土体的抗拉强度。

然而，目前关于水力劈裂发生的力学原因，主要有2 种观点，即拉裂破坏和剪切破坏[11]。本文在模拟水力劈裂的发展过程中，除了根据裂缝周围颗粒的速度方向定性地证明了水力劈裂破坏应属于颗粒间的张拉破坏外，还通过PFC2D程序分别统计并记录了法向和切向接触黏结破坏数目的变化过程，图9 显示的是初始竖向应力为0.2 MPa 时的结果，另外2 种方案的结果规律与之类似。图中初始的21条法向黏结破坏裂缝是由于在水平软弱带瞬时施加的水压力产生的。当运行到6 622 步时，模型内部形成了贯穿性裂缝，总的法向黏结破坏裂缝数目达到57 条。从该图可以看出，随着水力劈裂过程的发展，颗粒之间的接触黏结逐渐破坏，其中法向黏结破坏的数目随着时步在逐渐增加，而切向黏结破坏的个数却始终保持0，再一次证明水力劈裂现象的发生主要是由心墙内部法向拉裂破坏导致。

此外，从图9 中还可以看出，形成的裂缝数目与时步之间并不呈线性关系。其原因是，只有在裂缝尖端处积聚的能量足够使裂缝进一步扩展时，才能产生新的裂缝，因此，图中线段斜率较小的时步段对应的是裂缝尖端处能量积聚以及内部调整的时段，一旦能量足够大，就会形成大量新裂缝。

图9 法向裂缝数目-时步关系曲线 Fig.9 Curves of number of normal crack-step

在本次数值模拟中，从裂缝的发展方向（见图7）可以看出，水力劈裂面近似垂直于1σ 的作用方向，这表明在满足高水力梯度作用下，产生的水力楔劈效应降低了裂缝尖端区附近的1σ ，从图8 和表2 中可以看出，对于3 组不同的（1σ ，3σ ）组合，当1σ 小于或接近相同高程处心墙上游面的外水压力时，就会发生水力劈裂。

4 结 论

（1）本文运用PFC2D模拟了心墙水力劈裂的发生和发展过程，并将模型中产生连通上下游的贯穿性裂缝作为判别水力劈裂发生的依据。文中将数值模拟的结果和室内试验的成果做了对比，两者的规律基本一致，即随着竖向应力的增加所需的劈裂水压力值大致呈线性增加趋势。

（2）通过数值模拟的结果可以看出，水力劈裂发生的主要原因是心墙在高水力梯度作用下形成的水楔效应导致了心墙内裂缝尖端区的张拉破坏，并且破坏面是近似垂直于1σ 的作用方向，这表明水楔效应降低了心墙内原有的1σ ，当该值小于或接近心墙的外水压力时将会发生水力劈裂。

（3）由于模型在发生水力劈裂破坏的过程中法向黏结破坏的数目随着时步在不断增加，切向黏结破坏却始终未发生，进一步证明了水力劈裂的破坏形式属于法向张拉破坏。

（4）本文的研究成果还存在不足之处，尚需进一步改进。如：颗粒的大小及组装形式对计算结果有较大影响；计算参数的选取还没有较严格的理论依据；流-固耦合中还没有考虑到应力对渗透特性的影响；研究成果未能深入揭示水力劈裂发生的机制等。但本文的研究成果为今后进一步研究心墙的水力劈裂提供了新思路。

[1] 黄文熙. 对土石坝科研工作的几点看法[J]. 水利水电技术, 1982, (4): 23－27.

[2] 殷宗泽, 朱俊高, 袁俊平, 等. 心墙堆石坝的水力劈裂分析[J]. 水利学报, 2006, 37(11): 1348－1352. YIN Zong-ze, ZHU Jun-gao, YUAN Jun-ping, et al. Hydraulic fracture analysis of rock-fill dam with core wall[J]. Journal of Hydraulic Engineering, 2006, 37(11): 1348－1352.

[3] 沈珠江, 易进栋, 左元明. 土坝水力劈裂的离心模型试验及其分析[J]. 水利学报, 1994, (9): 67－77. SHEN Zhu-jiang, YI Jin-dong, ZUO Yuan-ming. Centrifuge model test of hydraulic fracture of earth dam and its analysis[J]. Journal of Hydraulic Engineering, 1994, (9): 67－77.

[4] 张丙印, 李娜, 李全明, 等. 土石坝水力劈裂发生机理及模型试验研究[J]. 岩土工程学报, 2005, 27(11): 1277－1281. ZHANG Bing-yin, LI Na, LI Quan-ming, et al. Mechanism analysis and model test of hydraulic fracturing in embankment dams[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(11): 1277－1281.

[5] 李全明, 张丙印, 于玉贞, 等. 土石坝水力劈裂发生过程的有限元数值模拟[J]. 岩土工程学报, 2007, 29(2): 212－217. LI Quan-ming, ZHANG Bing-yin, YU Yu-zhen, et al. Numerical simulation of the process of hydraulic fracturing in earth and rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(2): 212－217.

[6] 朱俊高, 王俊杰, 张辉. 土石心墙水力劈裂机制研究[J]. 岩土力学, 2007, 28(3): 487－492. ZHU Jun-gao, WANG Jun-jie, ZHANG Hui. Study on mechanism of hydraulic fracturing in core of earth-rockfill dam[J]. Rock and Soil Mechanics, 2007, 28(3): 487－492.

[7] 陈五一, 赵颜辉. 土石心墙水力劈裂计算方法研究[J]. 岩石力学与工程学报, 2008, 27(7): 1380－1385. CHEN Wu-yi, ZHAO Yan-hui. Study of calculation method of hydraulic fracturing for core of earth-rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(7): 1380－1385.

[8] 曹雪山, 殷宗泽. 土石坝心墙水力劈裂的非饱和土固结方法研究[J]. 岩土工程学报, 2009, 31(12): 1851－1856. CAO Xue-shan, YIN Zong-ze. Consolidation method of unsaturated soils for hydraulic fracturing of core walls of rock-fill dams[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(12): 1851－1856.

[9] 毕庆涛. 心墙水力劈裂的总应力法分析[J]. 水力发电, 2010, 36(3): 51－54. BI Qing-tao. Hydraulic fracturing analysis of core wall by total stress method[J]. Water Power, 2010, 36(3): 51－54.

[10] NG A K L, SMALL J C. A case study of hydraulic fracturing using finite element methods[J]. Canada Geotechnical Journal, 1999, 36(5): 861－875.

[11] Itasca Consulting Group Inc. PFC2D User’s Manual(version 3.1)[M]. Minneapolis: Itasca Consulting Group Inc., 2004.

[12] Itasca Consulting Group Inc. PFC3D verification problems and example applications[M]. Minneapolis: Itasca Consulting Group Inc., 2005.

[13] 王俊杰, 朱俊高, 张辉. 关于土石坝心墙水力劈裂研究的一些思考[J]. 岩石力学与工程学报, 2005, 24(增刊2): 5664－5668. WANG Jun-jie, ZHU Jun-gao, ZHANG Hui. Some ideas of study on hydraulic fracturing of core of earth-rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(Supp.2): 5664－5668.