W a t e r R e s o u r c e s a n dh y d r o p o w e r E n g i n e e r i n gV o l .36N o .2中国水利水电科学研究院的面板堆石坝数值分析研究的回顾与展望徐泽平(中国水利水电科学研究院,北京 100044)【摘 要】 主要介绍中国水利水电科学研究院在面板堆石坝数值计算分析方面的一些经验,以及有关峡谷地区高面板堆石坝、利用软岩修筑面板堆石坝和深覆盖层上面板堆石坝等方面的研究成果,同时还对面板堆石坝数值计算分析中需进一步研究的问题进行了讨论.【关键词】 面板堆石坝;数值分析;应力;变形;中国水利水电科学研究院中图分类号:T V 641.43 文献标识码:A 文章编号:1000-0860(2005)02-0039-05R e s e a r c h w o r k s o f I WH Ro nn u m e r i c a l a n a l y s i s o f C F R D :r e v i e wa n dp r o s p e c tX UZ e -p i n g(C h i n a I n s t i t u t e o f W a t e r R e s o u r c e s a n dH y d r o p o w e r R e s e a r c h ,B e i j i n g 100044,C h i n a )A b s t r a c t :S o m e e x p e r i e n c e s o f I WH R o nt h en u m e r i c a l a n a l y s i s o f C F R D a r ei n t r o d u c e dh e r e i n .I na d d i t i o n ,s e v e r a l r e -s e a r c hp r o j e c t s s u c h a s :h i g hC F R Db u i l t i nn a r r o wv a l l e y ,C F R Dc o n s t r u c t e dw i t hs o f t r o c k f i l l a n dC F R Dc o n s t r u c t e do n d e e p a l l u v i a l f o u n d a t i o na r e b r i e f l y i n t r o d u c e d .F u r t h e r m o r e ,s o m e p r o b l e m s t h a t n e e d t o b e f u r t h e r s t u d i e d a r e d i s c u s s e da s w e l l .K e yw o r d s :C F R D ;n u m e r i c a l a n a l y s i s ;s t r e s s ;d e f o r m a t i o n ;I WH R收稿日期:2004-07-11作者简介:徐泽平(1963—),男,教授级高级工程师,博士.1 概 述混凝土面板堆石坝是以堆石体为支承结构、并在其上游表面设置混凝土面板作为防渗结构的一种堆石坝形式.采用薄层碾压技术的现代面板堆石坝始于20世纪70年代,而中国于1985年也开始采用现代技术修建混凝土面板堆石坝.20多年来面板堆石坝的设计和施工技术取得了迅速的发展.在面板堆石坝的设计中,数值分析方法是一种有效的分析工具.通过数值计算分析,可以深入了解面板堆石坝在各种工况下坝体和面板的应力变形特性,预测坝体可能的变形与破坏形式,从而为坝体的结构设计与施工提供指导.伴随着中国面板堆石坝技术的发展,中国水利水电科学研究院(以下简称“中国水科院”)在面板堆石坝的数值计算方面进行了大量的研究工作.从中国第一座面板堆石坝———西北口面板堆石坝到目前在建的世界最高的面板堆石坝———水布垭面板堆石坝,中国水科院完成了大量实际工程的计算分析工作.表1所列为中国水科院参与计算分析的部分面板堆石坝工程.近些年来,在国家科技攻关项目和国电公司重点项目的资助下,中国水科院先后参加、完成了“200m 级高混凝土面板堆石坝关键技术研究”、“峡谷地区高混凝土面板堆石坝关键技术应用研究”、“利用软岩筑面板堆石坝关键技术应用研究”、“深覆盖层上面板堆石坝关键技术研究”、“超硬岩筑面板堆石坝关键技术应用研究”等项目的研究工作.同时,仍在继续进行“面板堆石坝长期变形特性研究”、“面板堆石坝湿化变形特性研究”、“面板堆石坝非连续变形分析方法研究”等项目的研究工作.39表1 中国水利水电科学研究院参与计算分析的部分面板堆石坝工程编号坝 名地点坝高/m筑坝材料1西北口湖北95灰岩2水布垭湖北233灰岩3天生桥一级贵州178灰岩4洪家渡贵州179.5灰岩5三板溪贵州185.56紫坪铺四川159砂岩、灰岩7吉林台新疆152凝灰岩8公伯峡青海130砂砾石9街面福建126重结晶坚硬泥岩、石英砂岩10芹山福建122凝灰岩11积石峡青海100黑云角闪片岩12龙首二级(西流水)甘肃146.5辉绿岩13九甸峡甘肃136.5灰岩14舍网浙江35.515溪口抽水蓄能电站上库浙江38.116溪口抽水蓄能电站下库浙江43.117张河湾抽水蓄能电站上库河北18西龙池抽水蓄能电站上库山西5219泰安抽水蓄能电站上库山东99.820十三陵抽水蓄能电站上库北京75安山岩21察汗乌苏新疆107.6砂砾石22云荞云南82.2灰岩23宜兴抽水蓄能电站上库江苏75石英砂岩夹泥岩24鱼跳重庆110砂岩25盘石头河南106砂岩、页岩26大坳江西90砂岩27柴石滩云南103白云岩2 数值计算分析方法2.1 堆石体的本构模型在面板堆石坝的设计施工中,堆石体的变形控制是影响坝体和面板工作性状的关键问题.就目前而言,有限元计算分析是进行面板堆石坝应力变形分析的主要手段,而材料的的本构关系又是有限元计算分析的核心问题之一,因此,堆石材料的本构模型对于面板堆石坝数值计算分析的结果有着至关重要的作用.在多年的计算分析实践中,中国水科院开发了基于多种本构模型的分析方法,其中主要包括:非线性弹性K-G模型、试验曲线插值模型和D u n c a n-C h a n g双曲线模型,目前正在进行基于剑桥弹塑性模型分析方法的开发.对于试验曲线插值模型,其模量的计算基于常规三轴试验的实测数据,通过插值计算的方法,由应力不变量、轴应变和八面体剪应变算出.模型的主要计算公式为 G=(Δτo c t/Δε1)/(Δγo c t/Δε1)(1)ν=(3/22)(Δγo c t/Δε1)-1(2)K=G(2/3)(1+ν)/(1-2ν)(3)σo c t=(σ1+2σ3)/3(4)J3=σ1(σ3)2(5)τo c t=(2/3)(σ1-σ3)(6)γo c t=(2/3)(3ε1-εν)(7)对于D u n c a n-C h a n g双曲线模型,在实际的计算中曾尝试了使用其E-ν模式和E-B模式.经计算分析比较,E-ν模式由于模量估算失真,计算结果与实测相差较大,因此,在堆石坝的分析中一般宜采用E-B模式.D u n c a n模型E-B模式的主要计算公式如下.切线弹性模量E t=K·P a·(σ3P a)n[1-R f(1-s i nφ)(σ1-σ3)2c c o sφ+2σ3s i nφ]2(8)切线体积模量B=K b·P a·(σ3P a)m(9)卸荷时,采用卸荷弹性模量E u r=K u r·P a·(σ3P a)n(10)式中,σ1和σ3为最大和最小主应力;P a为大气压力; c和φ为强度指标;R f为破坏比;K为弹性模量数;n为弹性模量指数;K b为体积模量数;m为体积模量指数; K u r为卸荷弹性模量数.在采用非线性弹性模型进行面板堆石坝的计算分析时,程序的处理上要注意各弹性常数间的相关关系.对于D u n c a n E-B模型,尤其要注意泊桑比ν的控制,过大或过小的ν值都将对计算结果的精度产生不利的影响.从大量计算分析的结果上看,由于堆石材料应力应变关系具有明显的非线性,因此,其本构模型必须准确反映这种非线性关系,线弹性模型对于堆石变形的计算是不适用的.另外,堆石料的剪缩特性对面板的应力有一定的影响,对于堆石剪缩特性的合理考虑,宜采用弹塑性模型.采用多屈服面、非关联流动法则的弹塑性模型将是未来面板堆石坝计算分析的发展方向,但这一类模型目前仍面临着试验方法特殊、计算参数类40水利水电技术 第36卷 2005年第2期水利水电技术 第36卷 2005年第2期比性差、以及计算复杂等问题.就目前而言,采用邓肯的E -B 模型并结合一些适当的修正,其计算分析的结果基本上是可信的,就工程实用的角度而言,其计算结果的精度也可以满足工程的要求.2.2 接触面及接缝系统的模拟在面板堆石坝的结构中,涉及到刚性的混凝土面板与散粒体堆石的相互作用面,以及面板纵缝、面板周边缝等接缝系统,因此,在面板堆石坝的数值计算分析中,需对坝体结构中各类不同材料的接触面和面板的分缝进行有效地处理.由于接触面两侧材料性质相差悬殊,在外力作用下,通常都会表现出与连续体不同的剪切滑移、脱开分离等特殊的变形特征,因此,在计算分析中需要采用特殊的单元来加以模拟,以准确、真实地反映坝体各部位相互作用的特性.在以往的计算分析中,接触面的模拟常采用无厚度的G o o d m a n 单元来模拟,其单元的受力变形与接触面两边节点的相对位移直接相关,对于G o o d m a n 单元,由于其无厚度的特性,在实际计算分析中,很难保证单元两侧的节点不发生相互嵌入的现象,而且法向劲度取值过大或过小对于计算的精度也会产生不利的影响.事实上,从实际工程的观测和试验室的试验中均可发现,在两种材料性质相差悬殊的介质之间,一般都会在材料性质相对较弱的一面形成一个薄层的剪切带,因此,采用薄层接触面单元来模拟不同材料之间的接触情况可能会更接近实际.对于薄层接触面单元,接触面上的变形可以分为基本变形和破坏变形两部分.在正常受力情况下,单元产生基本变形(ε′),其材料的本构关系与垫层料一致,薄层单元在计算过程中按普通实体单元参与计算.当剪应力达到抗剪强度产生了沿接触面的滑动破坏或接触面受拉产生了拉裂破坏时,单元产生破坏变形(ε″),破坏变形采用刚塑性假定,即假定接触面单元破坏前,接触面上无相对位移,当接触面发生张裂和剪切变形时,则相对位移将不断发展.对于接触面上的破坏变形(ε″),可以用下式表示Δε″s Δε″n Δγ″s n=0001E ″001G ″ΔσsΔσnΔτs n=C ″(Δσ)(11)式中,E ″和G ″分别是反映拉裂破坏变形和滑动破坏变形的模量参数.平行于接触面方向上的正应变由于受到混凝土的约束不会发生破坏,因此可取Δεs″=0,对应的,C ″矩阵中的相应元素取为0.接触面的总变形为基本变形和破坏变形的叠加(Δε)=(Δε′)+(Δε″)=C ′(Δσ)+C ″(Δσ)=C (Δσ)(12)对于面板之间的接缝,采用在分缝处设置双节点的分离缝单元模拟,缝单元采用无厚度形式.当接缝呈受拉趋势时,缝节点分离,两缝节点各自随垂直缝两边的面板单元位移;当接缝受压时,缝节点重合.对于趾板与面板之间的周边缝,则采用软单元的方式进行模拟.即在面板单元与趾板单元之间设置一个薄层软单元,当单元受压时,此软单元取混凝土材料的力学特性;当接缝受拉或受剪时,软单元则取为较低模量的柔性材料特性.需要指出的是,尽管薄层接触面单元在模拟接触面的剪应力传递和避免接触面两边介质的互相嵌入上具有一定的优势,但是,它对于接触面相互脱开的模拟仍存在缺陷.在未来的计算分析中,我们将进一步研究能够反映接触面非连续变形特性的界面单元的应用.2.3 施工及蓄水过程的仿真从大量的计算分析中可以发现,面板堆石坝的施工顺序及蓄水过程对于坝体的变形以及面板的应力有着明显的影响,因此,在计算分析中应真实地模拟坝体的实际填筑过程和面板的分期浇筑,尤其是坝体临时挡水断面及其后续填筑层的施工过程模拟.计算中每一填筑施工层作为一个增量荷载步进行迭代计算,水库蓄水的过程也应分成多个荷载步进行计算.计算分析的结果表明,在每一填筑施工层计算结束后,采用填筑层表面位移归零的处理方式,将会得到与观测结果较为一致的成果.在蓄水期,坝体堆石的上游侧由于水荷载的作用,主应力方向发生偏转,σ3数值明显增加,(σ1-σ3)数值减小,产生偏应力卸载现象.在蓄水期的计算中,应充分考虑这一卸荷过程,对相应的单元采用卸荷模量进行计算,以准确计算蓄水期面板的应力.3 工程应用研究在长期的面板堆石坝数值计算分析的实践过程中,中国水科院逐步形成并完善了一套完整的面板堆石坝应力变形分析技术,并进行了大量实际工程的计算分析和方案论证工作.其中包括面板堆石坝施工过程的仿真分析、岸坡地形对面板坝应力变形特性的影响分析、次堆石区坝料特性对坝体和面板应力变形特性的影响分析、面板堆石坝与坝基覆盖层的相互作用41分析、面板的温度应力与干缩应力分析、特殊结构形式(高趾墙、挤压式挡墙)的面板堆石坝应力变形分析等.以下所列为中国水科院近些年完成的部分面板堆石坝数值分析项目成果简介.3.1 利用软岩筑面板堆石坝对于坝高小于120m的面板堆石坝,采用软岩堆石(饱和抗压强度小于30M P a)作为筑坝材料,充分利用开挖料并提高料场的利用率,是面板坝工程中节约投资、缩短工期的重要途径.这项研究工作主要是通过对江西大坳面板堆石坝、重庆鱼跳面板堆石坝和河南盘石头面板堆石坝的数值计算研究,就软岩堆石分区对坝体应力变形特性的影响进行分析,并在此基础上提出软岩堆石分区布置和优化的一般原则.研究结果表明,对于软岩堆石料的利用,应保证软岩料区的下边界线在大坝运行时处于干燥区,以便坝体排水畅通,并避免软岩遇水产生湿化变形等;上边界线应保证其上游有不小于2m的新鲜硬岩填筑层;下游边界线应保证坝体下游边坡的稳定,且在其外侧留有不小于2m新鲜硬岩填筑区,以防止软岩料的继续风化;软岩堆石区的上游边界线应通过应力变形计算分析,在保证坝体施工期、运行期的沉降量以及面板的应力在合理范围内的前提下,尽量往坝体上游侧靠近,以期能够最大限度的利用软岩材料.3.2 峡谷地形下的面板堆石坝应力、变形特性研究对于峡谷地形条件下修建的高面板堆石坝,岸坡的约束和坝肩边坡的地形对于坝体和面板的应力和变形特性有着直接的影响,因此,在计算分析中,必须在综合考虑坝址地形条件的情况下,对坝体进行整体三维有限元计算分析.这项研究主要是通过数值计算的方法,结合179.5m高的洪家渡面板堆石坝工程,研究峡谷地形条件下高混凝土面板堆石坝的应力、变形特性.在计算分析中,通过考虑河谷地形和施工顺序对坝体及面板应力变形的影响,深入研究了洪家渡面板堆石坝在狭窄、不对称河谷地形条件下的应力、变形分布规律,以及面板周边缝的变形特点,并为坝体材料的合理分区及断面的优化设计提供依据.3.3 深覆层上的面板堆石坝对于修建于深覆盖层上的面板堆石坝,采用混凝土防渗墙构筑一道垂直防渗体系是工程中经常采用的有效方法.采用这样的防渗体系,防渗墙与上部坝体防渗体(面板)的连接将是整个防渗体系的关键部位.在这项专题的研究中,主要是通过对坝体和坝基整体结构在不同工况下应力变形的数值计算分析,深入研究覆盖层上面板堆石坝的应力变形特性、深覆盖层对上部坝体变形的影响规律以及坝体—面板—趾板—连接板—防渗墙—覆盖层之间的静力相互作用性态,对工程的结构形式和设计方案作出论证.4 今后的研究工作4.1 面板堆石坝的长期变形研究近些年来,随着面板堆石坝坝高的不断增长,以及所采用的筑坝材料日益广泛,坝体的长期变形问题也日趋突出.对于面板堆石坝,其长期变形主要包括水库运行期蓄水荷载变化引起的变形以及堆石体的流变变形等.目前,关于坝体堆石流变变形的机理研究尚不成熟.一般认为,在堆石坝中,堆石颗粒之间接触紧密,颗粒接触点处的局部应力集中将可能导致颗粒的局部破坏,从而引起周围堆石体的应力重新分布,同时也可能会引起颗粒的重新排列.所有这些调整的过程,均需要一定的时间完成,因而就会造成变形随时间不断发展的过程.另外,随着时间的推移,堆石体内理、化条件的变化,也可能会导致堆石变形的发展.就工程实践而言,目前最为关心的还是工程竣工后,坝体的后期变形量及其对面板应力和位移形态的影响.对于蓄水荷载变化引起的变形,可以通过常规的方法进行计算分析.对于坝体堆石流变变形量的确定,则取决于计算分析中采用什么样的流变模型.目前,对于堆石的流变分析,主要有两种方法,一种是采用理论模型,如M a x w e l l模型、V o g i t模型和M e r c h a n t模型等;另一种是采用经验函数模型,如指数衰减函数、幂函数、双曲函数和对数函数等.采用不同的流变模型可以得出不同的流变变形,即使是对于同一种流变模型,不同的参数取值也会导致不同的流变计算结果.因此,在堆石的流变变形研究中,材料特性试验不可或缺.目前,中国水科院的大型高压蠕变试验仪即将投入使用.今后的数值计算分析研究中,将充分结合试验研究的相关成果,研究开发堆石流变分析模型及相关参数的测定方法,对面板堆石坝的长期变形特性进行深入地研究.4.2 堆石材料的湿化变形研究堆石材料由干燥状态遇水变为湿态时所产生的变42水利水电技术 第36卷 2005年第2期形称为堆石的湿化变形.堆石在湿化过程中所产生的附加变形将会使坝体位移和应力重新分布,在大多数情况下,它将会导致坝体的不均匀变形,从而对面板的应力产生影响.由于面板堆石坝的防渗层位于坝体上游面,在面板未发生大的裂缝和止水系统完好的情况下,水库蓄水的升降不会造成堆石的湿化,因此,常规的数值计算分析中一般不考虑堆石的湿化变形.但是,当因面板裂缝或止水破坏而导致库水渗漏时,在一定的条件下,堆石的湿化变形将有可能导致面板的进一步损坏.另外,对于高面板堆石坝,由于雨水入渗和下游尾水位变化所引起的堆石湿化变形也有可能对坝体的整体变形和面板的应力造成不利的影响.在今后的研究中,中国水科院将通过应力控制式的三轴试验仪,对各种堆石料的湿化变形特性进行研究,并在此基础上建立堆石坝湿化变形的计算模型和分析方法,以深入研究堆石湿化变形对面板堆石坝(尤其是高面板堆石坝)应力变形特性的影响.5 结 语本文结合中国水科院从事面板堆石坝数值计算分析的实践,针对面板堆石坝应力变形分析中的主要问题进行了简要的分析.随着面板堆石坝技术的发展,未来的工程设计与施工中所面临的问题也日趋复杂.尤其是随着坝高的增长,对以往的一些工程经验和设计准则需要予以重新考虑与分析.在今后的研究工作中,对于复杂工程条件下的面板堆石坝工程,数值分析方法将在迎接新的挑战的基础上不断发展,与此同时,数值分析方法与物理模型分析方法的结合,也将会成为新的发展方向(限于篇幅,本文中未包括动力数值计算方面的内容).参考文献:[1] 蒋国澄,傅志安,凤家骥.混凝土面板坝工程[M].武汉:湖北科学技术出版社.[2] 钱家欢,殷宗泽.土工原理与计算(第二版)[M].北京:水利电力出版社,1994.[3] I n t e r n a t i o n a l C o m m i t t e e o nL a r g eD a m.R o c k f i l l D a m s w i t hC o n c r e t eF a c i n g,S t a t eo f t h e A r t[R].I C O L DB u l l e t i n70,1989.[4] J u s t oJ L,e t a l.T h e U p s t r e a mF a c i n g o f M a r t i n G o n z a l o R o c k f i l l D a m[A].T r a n s a c t i o n s16t h C o n g r e s s o n L a r g e D a m s[C].V o l.2,Q.61,R47,1988.[5] P i n t o AV,e t a l.D e s i g n C r i t e r i a o f U p s t r e a m C o n c r e t e F a c i n gR o c k-f i l l D a m s[A].T r a n s a c t i o n s16t h C o ng r e s so nL a r g eD a m s[C].V o l.2,Q.61,R47,1988.[6] S h e r a r dJ L,C o o k eJ B.T h ec o n c r e t ef a c er o c k f i l l d a m:I a s s e s s-m e n t[J].J o u r n a l o f G e o t e c h n i c a l E n g i n e e r i n g,A S C E,N e wY o r k,1985.[7] L e p s TM.R e v i e wo f s h e a r i n gs t r e n g t ho f r o c k f i l l[J].J o u r n a l o f s o i lm e c h a n i c s a n d f o u n d a t i o nd i v i s i o n,A S C E,1983,96(S M4).[8] D u n c a nJ M,C h a n gCY.N o n l i n e a r A n a l y s i s o f S t r e s s a n dS t r a i ni nS o i l s[A].P r o c.A S C E,J S M F D[C].V o l.96,S M5,1970.[9] D u n c a nJ M,B y r n e P,Wo n gKS,e t a l.S t r e n g t h,S t r e s s-S t r a i na n dB u l kM o d u l u s P a r a m e t e r s f o r F i n i t e E l e m e n t A n a l y s i s o f S t r e s s e s a n dM o v e m e n t si nS o i lM a s s e s,G e o t e c h n i c a lE n g i n e e r i n gR e p o r tN o.U C B/G T/80-01[R].U n i v e r s i t y o f C a l i f o r n i a,1980.[10] M a r t i nH L.A t h r e e-d i m e n s i o n a l d e f o r m a t i o na n a l y s i so f s t o r v a s sd a m[J].I n te r.J o u r n a l N u m e r i c a l A n a l y t i c a l M e t h o d si nG e o m e-c h a n i c s,V o l.2,p p3-17,1978.[11] D u n c a nJ M,Wo n g KS,O z a w a Y.F E A D A M:AC o m p u t e r P r o g r a mf o r F i n i t e E l e m e n t A n a l y s i s o f D a m s[R].U n i v e r s i t yo f C a l i f o r n i a,G e o t e c h n i c a l E n g i n e e r i n g R e p o r t N o.U C B/G T/80-02,1980.[12] M o r i R T,P i n t oN L.A n a l y s i s o f D e f o r m a t i o n si nC o n c r e t eF a c eR o c k f i l l D a m s t oI m p r o v eF a c e M o v e m e n t P r e d i c t i o n[A],T r a n s-a c t i o n s16t hC o n g r e s so nL a r g eD a m s[C].V o l.2,Q.61,R2,1988.[13] I m a i z u m i H,S a r d i n h aAE.A S t u d yo f D e f o r m a t i o n s i nC o n c r e t eF a c e d R o c k f i l l D a m s,C o n c r e t e F a c eR o c k f i l l D a m—D e s i g n,C o n-s t r u c t i o n a n d P e r f o r m a n c e[M].E d i t e db yJ.B.C o o k ea n dJ.L.S h e r a r d,1985,528-540.[14] W o o dD M.S o i l b e h a v i o r a n dc r i t i c a l s t a t es o i l m e c h a n i c s[M].C a m b r i d g e U n i v e r s i t y P r e s s.1990.[15] G o o d m a nRE,T a y l o r RL,B y e k k e TL.Am o d e l f o r t h e m e c h a n i c so f j o i n t e d r o c k[J].J o u r n a l o f s o i l m e c h a n i c sa n df o u n d a t i o nd i v i-s i o n,A S C E,V o l.94,S M31968.[16] D e s a i CS,e t a l.T h i n l a y e r e l e m e n t f o r i n t e r f a c e a n d j o i n t s[J].I n-t e r.J o u r n a l f o r N u m.＆M e t h.I nG e o m e c h,1984,8(1):19-43.[17] N o g u e r a G,P i n i l l aL,S a nM a r t i nL.2000,C F R DC o n s t r u c t e do nD e e pA l l u v i u m,J.B A R R Y C O O KE R V O L U M E—C o n c r e t ef a c eR o c k f i l l D a m s[M].[18] I m a i z u m i H,e t a l.As t u d yo f d e f o r m a t i o ni nc o n c r e t e f a c e dr o c k f i l ld a m s,d i s c u s s i o n b yW.H.S c h u m a n n,e ta l[J].J o u r n a lo fG e o t e c h n i c a l E n g i n e e r i n g.A S C E,V o l.113,N o.10,1264-1268,1987.[19] ADM P e n m a n.T h ew a t e r p r o o f e l e m e n t f o r e m b a n k m e n t d a m s[J].Wa t e r P o w e r a n d D a mC o n s t r u c t i o n,V o l.37,N o.7,1985.(责任编辑 欧阳越)植树造林 美化环境43水利水电技术 第36卷 2005年第2期。
