裂隙岩体渗流场与应力场动态全耦合参数反演

1、第27卷第8期岩石力学与工程学报V ol.27 No.8 2008年8月Chinese Journal of Rock Mechanics and Engineering Aug.,2008 裂隙岩体渗流场与应力场动态全耦合参数反演王 媛1,2,刘 杰3(1. 河海大学岩土力学与堤坝工程教育部重点实验室,江苏南京 210098;2. 河海大学岩土工程研究所,江苏南京 210098;3. 广州市建设委员会,广东广州 510032摘要:考虑到裂隙岩体渗流和力学行为的发生与发展往往是一个动态过程,同时现场观测资料也是一个动态的时间序列,为能及时利用现场量测的新增信息使参数反演更为合理,基于求解非恒定

2、渗流场与弹性位移场动态全耦合正分析理论与方法,应用建立的混合遗传算法作为优化算法,同时利用水头、位移等多类型动态观测资料,建立了裂隙岩体渗流场与应力场动态全耦合的参数反演思路。为避免在耦合反问题中由于利用多类型量测资料所带来的量纲问题,采用了各时刻水头、位移的相对值来构造量纲一的目标函数。待反演参数同时考虑了力学参数与渗流参数两种类型,包括岩块的弹性模量、各组裂隙的切向与法向刚度系数、各组裂隙的初始等效渗透系数等。最后以一简单裂隙岸坡为算例,针对库水位快速上涨情况,以各时刻的动态全耦合正分析结果作为“假想”的实测数据,进行动态全耦合参数反演。反演结果表明,利用不断新增的实测资料可提高反演精度,

3、最终获得的参数反演解与理论解吻合很好。关键词:岩石力学;裂隙岩体;参数反问题;混合遗传算法;渗流与应力耦合分析中图分类号:TU 45;TD 32 文献标识码:A 文章编号:10006915(200808165207INVERSE ANALYSIS OF FULLY COUPLED DYNAMIC WATER FLOW AND STRESS IN FRACTURED ROCK MASSESWANG Yuan1,2,LIU Jie3(1.Key Laboratory for Geomechanics and Embankment Engineering of Ministry of Educati

4、on,Hohai University,Nanjing,Jiangsu 210098,China;2. Institute of Geotechnical Engineering,Hohai University,Nanjing,Jiangsu210098,China;3. Construction Committee of Guangzhou City,Guangzhou,Guangdong510032,ChinaAbstract:Due to the dynamic process of water flow and deformation in fractured rock masses

5、 as well as time-dependent characteristics of in-situ monitoring data,a dynamic inverse method for fully coupled problem of water flow and stress is presented,in which a hybrid genetic algorithm is used for optimization;and two different types of monitoring data about water head and displacement are

6、 taken into account. In order to avoid the dimension problem caused by different types of monitoring data,related values of water head and displacement at each time step are used in building objective function. In the coupling inverse analysis,both mechanical and seepage parameters are regarded as u

7、nknown variables,such as elastic modulus of rock block,shear stiffness and normal stiffness of each fracture set,and initial equivalent permeability coefficient of each fracture set. Finally,the presented inverse method is applied to a simple example of fractured rock bank slope in case that water l

8、evel of reservoir rises quickly,while the forward calculated results at each time step are regarded as the assumed收稿日期:20071210;修回日期:20080622基金项目:国家自然科学基金雅砻江水电联合研究基金项目(50579090,50639090;科技部“十一五”科技支撑计划项目(2006BAB04A06作者简介:王媛(1969,女,博士,1989年毕业于河海大学水工专业,现任教授、博士生导师,主要从事岩土工程渗流理论与测试、裂隙岩体应力和渗流耦合等方面的教学与研究工作。

9、E-mail:wangyuan第27卷第8期王媛,等. 裂隙岩体渗流场与应力场动态全耦合参数反演 1653 monitoring data. It is indicated that the accuracy of parameters identification can be improved if using the continuously increasing monitoring data in time;and the inverse results of parameters are in good agreement with theoretical solution.Key

10、 words:rock mechanics;fractured rock masses;parameter inverse problem;hybrid genetic algorithms; coupled analysis of water flow and stress1 引言岩体的渗流特性与裂隙的几何参数特别是隙宽密切相关,而裂隙的几何参数又会随岩体应力场变化而发生显著改变,因此裂隙岩体渗流场与应力场之间存在着强烈的耦合作用,已成为人们的共识。国内外学者18已进行了大量的两场耦合研究工作,然而如何确定渗流、应力及其耦合的参数一直是该领域研究的难点和关键问题,为此笔者9曾提出全耦合参数静

11、态反演的思路,基于水头、位移等多类型观测资料和全耦合正分析的方法,实现了渗流、应力的参数静态反演,如此不仅能够提高传统的参数反演方法的结果可信度,而且为确定这些不易直接测定的参数提供了一种有力的工具。但是在实际工程中,由于裂隙岩体渗流和力学行为的发生、发展往往是一个动态过程,它随着工程的施工、环境的变化和时间的持续在不断变化,而且有关岩体渗流和力学行为的观测资料往往也是若干时间序列,因此如果不考虑岩体的这种动态行为,只通过某一时刻的观测资料进行反演分析即静态反演,其反演结果在非恒定渗流场或者荷载与边界条件变化的情况下将会失真,而且也不便应用于岩体渗流和力学行为的动态预测。于是有必要开展基于动态

12、观测资料的两场耦合的动态反演方法研究,并建立裂隙岩体渗流和变形行为的动态预测模型。目前,有关渗流场和变形场单场动态反演的研究成果已有少量报道,而考虑两场全耦合作用的动态反演研究成果尚未见报道。岩土工程参数反分析方法1014很多,鉴于工程实践中裂隙岩体渗流场及其与应力场动态全耦合正问题比较复杂,往往具有明显的非线性特征,同时考虑到反演研究的工程实用性,本文仍采用新型混合遗传数值优化反演方法9开展裂隙岩体渗流场及其与应力场动态全耦合的参数反演研究。2 裂隙岩体渗流场与应力场动态全耦合方程笔者已建立了非恒定渗流场与应力场弹性动态全耦合分析的基本控制方程及增量形式的有限元方程组15。基本控制方程为(2

13、1=+iikjlljkijklfpuud,=lkji,(1,2,3 (1atuzpkH k kijij=+,w=ji,(1,2,3 (1b式中:iiu,(i = 1,2,3为位移分量;p为渗流水压力;ip,为渗流水压力对坐标x i的偏导数;ijkld为裂隙岩体弹性矩阵元素,可根据应变等效建立,与岩块的弹性模量rE、泊松比r以及各组裂隙的法向刚度系数lnK、切向刚度系数lsK、平均间距l d和法向余弦lin有关16;if为体积力;w为水的重度;z为位置坐标;t为时间;eH为改进初流量法处理自由面时以渗流水压力表示的区域识别函数16;ijk为裂隙岩体的渗透张量元素,可根据流量等效建立,与各组裂隙的

14、等效渗透张量与法向余弦lin有关,各组裂隙的等效渗透张量取决于隙宽l b、平均间距l d及流体的性质如密度、黏滞系数等16。考虑到隙宽受应力的影响,可将ijk表征为各组裂隙初始等效渗透系数与应力的函数。所建立的有限元支配方程为=FFpuKKKT(2 式中:pu,分别为未知节点位移和水压力增量列阵;KK,K,F及F的含义详见有关研究结果15,不再赘述。由于K与渗透张量有关,是应力张量的函数,同时K中包含自由面边界积分项,而对于某一级荷 1654 岩石力学与工程学报 2008年载下的每一个计算时步t 而言,自由面位置事先并不确定,因此K 也是随自由面位置的变动而改变的,因此以上建立的动态耦合有限元

15、方程组是一个强烈的非线性方程组,在每一个计算时步t 内都需要迭代求解。3 裂隙岩体渗流场与应力场动态全耦合反演方法建立王 媛和刘 杰9已建立了裂隙岩体渗流场与应力场静态全耦合参数反演的思路,与静态全耦合反演相比,动态全耦合参数反演的目标泛函是与时间相关的,实测点数目随时间也在不断增加,因此随着计算时间的增加,目标泛函对优化解的全局最优性要求越来越高,此外,动态耦合正分析过程需要的计算量比静态问题也大很多,因此动态耦合参数反演问题对优化算法的要求更高。于是本文仍采用对多变量、大空间优化问题适应性强、计算效率高的新型混合遗传算法9进行动态耦合参数反演的尝试性研究。基于最小二乘准则,考虑到水头和位移

16、量测信息是两个时间序列,本文构造如下动态全耦合参数反问题的目标泛函为=+=Tt O N i it it it itt H H H E 112obs obs 20(X =t O N i it it it u u u 12obs obs(3 式中:it H 为第i 个水头测点t 时刻的水头计算值,obs it H 为相应的实测值,t O N 为t 时刻水头的测点数,it 为相应的权函数,it u 为第i 个位移测点t 时刻的位移计算值,obsitu 为相应的位移实测值。 式(3所示的目标函数中同样采用了水头、位移的相对值,使目标函数成为量纲一的数值函数,从而可以避免在耦合反问题中利用多类型量测资料

17、所带来的量纲问题。为了更有效地利用水头实测点,采用类似静态的处理手段9,本文采用相应时刻的实测点水头对渗流参数的敏感性分析构造权函数it ,可表示为 2 1( |/|max |/|max 1M j x H x H t O N i j it jj it jit =,L (4式中:M 为待反演的渗流参数个数,j x 为待反演的渗流参数。至此,裂隙岩体渗流场与应力场动态耦合参数反问题即转化为式(3,(4构成的多变量约束非线性优化问题。动态全耦合反分析中待反演未知参数同时包括力学参数与渗流参数。力学参数的选取与静态反演类似,仍采用岩块的弹性模量及各组裂隙的切向、法向刚度系数。渗流参数虽然仍采用裂隙组的

18、初始等效渗透系数而不是宏观的等效渗透张量,但与静态反演不同的是,由于动态耦合分析中水头及位移的分布与各裂隙组等效渗透系数的绝对值有关,因此在动态全耦合反演中,需要以各裂隙组初始渗透系数的绝对值作为目标未知渗流参数,而在静态反演中则采用的是相对值。基于上述建立的目标函数表达式、待反演参数的类型及新型混合遗传的优化算法,考虑到裂隙岩体渗流场与应力场的动态变化特化特性,可建立如下的动态全耦合反演方法的思路:(1 针对具体的工程,首先依据相关方法进行参数敏感性分析9,确定力学目标未知参数。(2 根据试验成果或者工程经验,限定所有目标未知参数的变化范围,包括渗流和力学所有待反演的参数。(3 应用新型混合

19、遗传优化算法,结合动态正分析计算,根据已量测各时刻的水头及位移信息,反演获得待反演的渗流、力学参数。(4 以上一次反演得到的参数作为输入参数作动态耦合正分析计算,预测下一量测时刻的渗流及变形。(5 将下一量测时刻现场实测水头、位移实测值作为新增的量测信息,转入步骤(2重新做动态反演计算;对上一次的反演参数进行调整,并依此作为输入参数,预测下一量测时刻的渗流及变形行为;如此反复进行。鉴于裂隙岩体渗流场与应力场动态全耦合反分析问题的复杂性,在本文初步的研究中认为动态的变形过程是由外部荷载的变化引起的,暂不考虑岩体变形特性的时效性。4 渗流场与应力场动态全耦合参数反 演分析算例结合笔者曾研制的裂隙岩

20、体渗流场与应力场动态全耦合的三维有限元正分析程序DCou 15,依据上述所建立全耦合动态反演分析思路,编制了三维第27卷第8期王媛,等. 裂隙岩体渗流场与应力场动态全耦合参数反演 1655 动态耦合反分析程序DCouHGA。为验证本文所建立理论与研制程序的正确性,采用王媛和刘杰15的裂隙岩坡算例,进行渗流场与应力场的三维动态全耦合参数反演分析。裂隙岩坡中发育两组裂隙,采用等效连续介质模型模拟,岩坡几何尺寸及裂隙发育的几何参数同王媛和刘杰15的研究。根据静态力学敏感性分析结果,将水压力测点与位移测点按图1所示位置布置(节点563,565,579,581,603及605为位移测点,其余节点为水压力

21、测点,而且图1中测点号表示对应的有限元节点号。本文首先取用与王媛和刘杰15相同的材料计算参数,并针对其第1种工况库水位快速上涨(库水位以12 m/d速度上升、3 d时间可达到最高水位36 m情况,进行动态耦合有限元正分析,再利用计算出的各级荷载结束时的各测点位移及水压力值,作为假想“实测值”进行反演分析。表16分别给出了动态耦合有限元分析中各级荷载对应的测点位移及水压力假想“实测”值。其中,表2中存有负压,本文不以负压作为反演依据,因此这类测点属于“无效”测点。需要说明的是,与静态反演不同,动态反演中测点数目随加载过程不断增加,即本级荷载完成时,反演采用的实测值同时包括前面各级荷载的测点实测值

22、,本级荷载施加后获得的实测值仅是“新增”的实测值。 图1 裂隙岸坡断面( y = 3 m及测点布置示意图Fig.1 Cross-section of fractured rock slope( y = 3 m and location of measuring points表1 上游水位12 m时对应的动态耦合反演位移“实测”值Table 1 Assumed monitoring data of displacements when water level of upstream reaches 12 m mm 节点编号水平位移竖向位移21 6.37-1.57463 6.81-1.554657

23、.28-1.53461 6.37-0.13543 6.81-0.155457.30-0.17表2 上游水位12 m时对应的动态耦合反演水压力“实测”值Table 2 Assumed monitoring data of water pressures when water level of upstream reaches 12 m 104 Pa节点编号水压力节点编号水压力565*563 5.83581*579 5.74605*603 4.49注:断面位于y = 3.0 m处;“*”表示测点存在负压。表3 上游水位24 m时新增动态耦合反演位移“实测”值Table 3 Assumed moni

24、toring data of displacements when water level of upstream reaches 24 m mm节点编号水平位移竖向位移219.05-4.394639.43-4.3546510.49-4.324619.05-0.675439.43-0.6254510.79-0.67表4 上游水位24 m时新增动态耦合反演水压力“实测”值Table 4 Assumed monitoring data of water pressures when water level of upstream reaches 24 m 104 Pa节点编号水压力节点编号水压力5

25、6510.9656317.63581 6.4557913.10605 2.846038.61注:断面位于y = 3.0 m处。表5 上游水位36 m时新增的动态耦合反演位移“实测”值Table 5 Assumed monitoring data of displacements when water level of upstream reaches 36 m mm节点编号水平位移竖向位移2116.31-11.6346316.56-10.6346516.85-10.5346116.21-1.1954316.45-1.2954516.74-1.33 1656 岩石力学与工程学报 2008年表6

26、上游水位36 m 时新增的动态耦合反演水压力“实测”值Table 6 Assumed monitoring data of water pressures whenwater level of upstream reaches 36 m 104Pa节点编号 水压力节点编号水压力 565 22.52 563 29.77 581 17.13 579 23.54 6058.4760314.28注:断面位于y = 3.0 m 处。动态耦合参数反演的目标未知参数向量取为=02f 01f 2s 2n 1s 1n r k k K K K K E ,X (5式中:i iK K s n ,(1=i ,2分别为裂

27、隙的法向和切向刚度系数;0f i k (1=i ,2为初始等效渗透系数。动态全耦合参数反演分析中,以表7所示的参数取值范围作为HGA 的优化搜索范围,取种群个体总数=ps N 30,交叉概率=c p 0.9,变异概率=m p 0.1,加速循环的最优个体数目NA 取为ps N 的10%,需单纯形局部搜索的个体数目NS 取为ps N 的10%。各级加载结束时对应的参数反演结果及一些程序运行效率指标列于表810中。由表8可以看出,第一级荷载施加后(水位上升到12 m,利用获得的有效实测资料反演效果并不理想,目标函数虽然很小,但参数的反演解和理论解仍然相差很大,尤其是裂隙组的初始渗透系数误差较大,分析

28、出现这一现象的主要原因是由于渗流实测资料信息不足,致使反演解不惟一,另外第一级荷载作用下岸坡的位移和水压力很小,依据这样的实测资料反演,反演计算本身的误差也比较大,同样会造成目标函数很小,但反演解和理论解仍然相差较大。而由表9,10可以看到,随着荷载的施加,利用新增的实测资料参加反演之后,获得的参数反演解与理论解相当吻合,这也表明在测点有限的情况下,采用动态反演方法,利用不断增加的实测资料反演调整参数反演解是必要的,也是有效的。由表810提供的反演程序运行时间来看,动态反演的计算量巨大,反演耗时明显大于静态情况。但随着未来工程数值计算效率的提高,计算量的限制可望将逐渐消除。表7 裂隙岸坡动态耦

29、合目标未知参数限值Table 7 Limit values of unknown parameters for dynamic inverse analysis类别r E /MPa1n K /(MPa ·m -11s K /(MPa ·m -12n K /(MPa ·m -12s K /(MPa ·m -1f1k /(m ·s -1f2k /(m ·s -1X 下限值 50 1 0.05 1 0.05 1×10-5 1×10-5X 上限值100 5 0.50 5 0.50 1×10-31×10-

30、3表8 上游水位12 m 时对应的动耦合参数反演结果Table 8 Parameters by dynamic inverse analysis when upstream water level is 12 m类别r E /GPa 1n K /(GPa ·m -1 1s K /(GPa ·m -12n K /(GPa ·m -12s K /(GPa ·m -1f1k /(m ·s -1f2k /(m ·s -1目标函数值X 理论解 67.000 2.000 0 0.100 0 3.000 0 0.150 0 1.000 0×

31、;10-4 6.000 0×10-5 X 反演解67.0681.974 70.099 43.028 00.150 41.253 0×10-45. 635 0×10-54.3×10-5相对误差 1.0 13.3 6.2 9.3 2.4 253.0 60.8注:计算条件及效率:CPU :P4 2.1 GHz ;内存资源:DDR 256 M ;反演程序用时约 9.0 h 。表9 上游水位24 m 时对应的动耦合参数反演结果Table 9 Parameters by dynamic inverse analysis when upstream water lev

32、el is 24 m类别r E /GPa 1n K /(GPa ·m -1 1s K /(GPa ·m -12n K /(GPa ·m -12s K /(GPa ·m -1f1k /(m ·s -1f2k /(m ·s -1目标函数值X 理论解 67.000 0 2.000 0 0.100 0 3.000 0 0.150 0 1.000 0×10-4 6.000 0×10-5 X 反演解67.073 02.002 50.100 32.987 40.149 51.002 7×10-45.998 6×

33、;10-51.1×10-4相对误差 1.1 1.3 3.1 4.2 3.1 2.7 1.4注:计算条件及效率:CPU :P4 2.1 GHz ;内存资源:DDR 256 M ;反演程序用时约 21.0 h 。第 27 卷 第8期 王 媛, 裂隙岩体渗流场与应力场动态全耦合参数反演 等. 1657 表 10 Table 10 类别 上游水位 36 m 时对应的动耦合参数反演结果 Parameters by dynamic inverse analysis when upstream water level is 36 m - Er /GPa 67.000 0 67.065 0 1.0

34、1 K n /(GPam 1 1 K s /(GPam 1 - 2 K n /(GPam 1 - K s2 /(GPam 1 0.150 0 0.150 2 1.0 - 0 kf1 /(ms 1 - 0 kf2 /(ms 1 - 目标函数值 X 理论解 X 反演解 相对误差 2.000 0 2.004 2 2.1 0.100 0 0.099 8 2.0 3.000 0 2.998 7 0.4 1.000 0×10 1.003 8×10 3.8 -4 6.000 0×10 6.000 1×10 0.1 -5 6.7×10 -4 -5 -4 注:计

35、算条件及效率:CPU:P4 2.1 GHz;内存资源:DDR 256 M;反演程序用时约 33.0 h. 流固耦合数学模型J. 水利学报,2005,36(4:405412.(LIU 5 结 论 Xiaoli,LIANG Bing,WANG Sijing,et al. Fluid-solid coupled mathematical model for water-air two-phase infiltration and deformation (1 本文针对裂隙岩体渗流场与应力场两场动 态全耦合反演问题进行了尝试性研究,建立了动态 全耦合的反演分析方法,算例的应用结果表明,该 法是可行,有

36、效的. 3 of dual porous mediaJ. Journal of Hydraulic Engineering,2005, 36(4:405412.(in Chinese 孙 跃. 流形元与有限元变形和渗流的非同步耦合分析方法及其 应用J. 岩石力学与工程学报,2003,22(6:943950.(SUN Yue. Application research on non-simultaneous coupling calculation of deformation and seepage with manifold method and FEMJ. Chinese Journal

37、of Rock Mechanics and Engineering,2003,22(6:943 950.(in Chinese 4 张金才,王建学. 岩体应力与渗流的耦合及其工程应用J. 岩石力 学与工程学报, 2006, 25(10: 9811 989.(ZHANG Jincai, 1 WANG Jianxue. Coupled behavior of stress and permeability and its engineering applicationsJ. Chinese Journal of Rock Mechanics and Engineering,2006,25(10:1

38、 9811 989.(in Chinese 5 陈卫忠,杨建平,杨家岭,等. 裂隙岩体应力渗流耦合模型在压 力隧洞工程中的应用J. 岩石力学与工程学报,2006,25(12: 2 3842 391.(CHEN Weizhong,YANG Jianping,YANG Jialing, et al. Hydromechanical coupled model of jointed rock mass and its application to pressure tunnelJ. Chinese Journal of Rock Mechanics and Engineering,2006,25(1

39、2:2 3842 391.(in Chinese 6 殷黎明,杨春和,王贵宾,等. 地应力对裂隙岩体渗流特性影响 的研究J. 岩石力学与工程学报, 2005, 24(17: 0713 075.(YIN 3 Liming,YANG Chunhe,WANG Guibin,et al. Study on effect of (2 现场观测资料往往是一个动态的时间序 列,在参数反演过程中及时利用现场量测的新增信 息进行参数反演,可提高反演精度,动态反演在指 导工程施工反馈设计,进行工程安全预报等方面的 应用价值不容忽视. (3 动态反演的计算量明显大于静态反演,可 望随着工程数值计算效率的提高,计算量

40、的限制将 逐渐消除. (4 由于裂隙岩体渗流场与应力场动态全耦合分 析的复杂性,动态全耦合参数反问题对优化算法的 要求更高,本文虽采用多变量大空间优化问题适应 性强,计算效率高的新型混合遗传优化算法,但仍 存在计算工作量大的缺点,建议进一步开展并行算 法在反演分析中的应用,推动反演技术在实践中的 应用. 参考文献(References: 1 吉小明,杨春和,白世伟. 岩体结构与岩体水力耦合计算模型J. 岩土力学,2006,27(5:763768.(JI Xiaoming,YANG Chunhe, BAI Shiwei. Structure and hydromechanical coupled

41、 calculation model for rock massJ. Rock and Soil Mechanics,2006,27(5: 763768.(in Chinese 7 geostress on permeability of fractured rock massJ. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(17: 0713 075.(in 3 Chinese 陈卫忠,邵建富,DUVEAU G,等. 黏土岩饱和非饱和渗流应 力耦合模型及数值模拟研究J. 岩石力学与工程学报,2005, 24(17:3

42、0113 016.(CHEN Weizhong,SHAO Jianfu,DUVEAU 2 刘晓丽,梁 冰,王思敬,等. 水气二相渗流与双重介质变形的 1658 岩石力学与工程学报 2008 年 G,et al. Constitutive model of saturated-unsaturated clay and its numerical simulationJ. Chinese Journal of Rock Mechanics and Engineering,2005,24(17:3 0113 016.(in Chinese 8 何 翔,冯夏庭,张东晓. 岩体渗流应力耦合有限元计算的精

43、 神经网络在渗流反分析中的应用J. 岩土力学,2005,26(3: 404 408.(LIU Xianshan ,ZHOU Chuangbing,ZHANG Lijun. Application of neural network based on simulated annealing Gauss-Newton algorithm to seepage back analysisJ. Rock and Soil Mechanics,2005,26(3:404408.(in Chinese 13 刘先珊, 佘成学, 张立君. RBFNN 模型在渗透参数反演中的应用J. 岩土力学,2003,24

44、(6:1 0251 028.(LIU Xianshan,SHE Chengxue,ZHANG Lijun. Application of RBFNN model to inversion of percolation parameterJ. Rock and Soil Mechanics,2003,24(6: 1 0251 028.(in Chinese 14 王 媛, 刘 杰. 重力坝坝基渗透参数进化反演分析J. 岩土工程 细积分方法J. 岩石力学与工程学报, 2006, 25(10: 0032 008. 2 (HE Xiang,FENG Xiating,ZHANG Dongxiao. Precise integration algorithm for FEM simulation of coupled process of seepage field and stress field in rock massJ. Chinese Journal of Rock Mechanics and Engineering,2006,25(10:2 0032 008.(in Ch