版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、赵冲久1 ,秦崇仁2 ,杨华1 ,曹祖德1 ,蔡嘉熙1 ,冯玉林1 ,杨树森1(11 交通部天津水运工程科学研究所 ,天津 300456 ;21 天津大学建工学院 ,天津 300072)摘 要 :研究了波浪作用下悬移质含沙量的垂线分布规律 、底部高浓度含沙水体层的高度和含沙量分布以及该层中水体的运移速度 。从而得出了底部高浓度含沙水体层的输沙量计算公式 ,为粉沙质床面的 泥沙淤积计算提供了新的方法 。关键词 :粉沙 ;悬移质 ;底部高浓度含沙水体层中图分类号 : TV142 + 13文章编号 :1005 - 8443 (2003) 03 - 0101 - 08文献标识码 :A悬质分布规律的研究
2、1粉沙质海岸当风平浪静时水体含沙量较低 ,水体清澈 ,一般没有泥沙的剧烈运动 。当波浪增大时 ,由于水体紊动的增强 ,床面泥沙开始运动 ,并扩散到水体当中形成悬移质泥沙 。波浪较大时泥沙运动十分剧烈 ,水体含沙量急剧加大 。研究粉沙质床面的泥沙运动 ,须认真研究悬移质泥沙的运动规律 。 由于粉沙质海岸表现为大浪时泥沙运动为主 ,大浪时泥沙灾害较大 。结合现场观测每次大风 ( 浪) 过程中泥沙运动的变化过程 。悬移质泥沙从扬起到沉降可以分为 3 种型式 。第 1 种是悬扬型 ,在大浪形成时 ,水体运动剧烈 ,底部泥沙起动并扬动 ,泥沙进入悬浮状态 ,此时底部水体与上部水体含沙量呈下大上小 ,差别
3、较 大 ;第 2 种是稳定型 ,随着波浪持续时间的加长 ,泥沙不断地被悬扬 ,水体含沙量不断增大 ,并达到相对平衡 时期 ,此时上下层水体含沙量趋于均匀 ;第 3 种是沉降型 ,当波浪动力逐渐减弱时 ,悬沙开始沉降 ,此时上下水体的含沙量又呈现出上小下大的较大差别 。风平浪静后 ,最终恢复平静状态下清澈的水体和稳定的床面 。总体上讲 ,在一次大浪过程中 ,下层水体出现由不饱和挟沙到饱和挟沙的过程 ,而上层水体一直处于非饱和挟沙的过程 。111理论分析文献 11 概化了底部含沙量计算公式 :2uS = A ·s(1)ghu2 = u2 + u2w c由于底部掀沙的主体原因是波浪 ,忽略
4、 uc 项 、uw 项用近底水质运动的最大轨迹速度来表示 ,则式 (1) 变为 :u2mS = A ·s(2)ghD0F11, D 为 特 定 粒 径 , 取 0111 mm , 为 特 定 面 积 , 取引入泥沙粒径的修正 系 数, 并 取 F =0D K + D K01002 4 mm2 , D K 为泥沙粒径 。则底部水体含沙量为 :u21mS K = A ·s(3)gh F F收稿日期 :2003 - 06 - 27作者简介 :赵冲久 (1965 - ) ,男 ,辽宁盘锦人 ,博士 ,教授级高工 ,主要从事海岸河口泥沙研究 。uh = Tshkh221S K = A
5、 ·s于是 ,A 由实测资料所得 。(4)F FT2 ghsh2 kh下面讨论上部水体中的含沙量问题 。尽管 Nielson (1992) 15 已经指出 ,悬沙浓度的合理分布应该采用联合对流和扩散方程来描述 ,但大多数 关于悬移质泥沙含沙量分布的理论研究仍是基于泥沙的扩散方程进行的 。在泥沙粒径比较小的情况下 ,利 用扩散方程也确实能够得到较合理的结果 。取 x 轴和波向一致 、y 轴与波向垂直 、z 轴向上为正 ,波浪作用 下二维悬沙扩散的基本方程为 :d s 9 s 9 s 9 s 9 s 99 s 99 sd t = 9 t + u 9 x + 9 z = 9 z + 9x
6、(x 9 x ) + 9 z (z 9 z )(5)式中 : s 为空间某一点的含沙量 (或悬沙浓度) ,以单位体积水体内所含的泥沙重量计 ; 为泥沙的沉降速度 ;x 、z 为 x 、z 方向的泥沙扩散系数 。Nakatoo (1977) 16 的研究结果表明 ,水平方向的悬沙浓度梯度要比竖直方向的小的多 。因此 ,如果忽略 水平方向的泥沙扩散 ,且取其空间平均值并考虑时间平均的情况 ,二维扩散方程可化简为描述含沙量平均值 垂向分布的一维扩散方程 :z 9 s +¯s = 0(6)9 z式中 : s¯表示时均悬沙浓度 。求解的关键在于要已知泥沙扩散系数s 。在悬沙扩散系数的
7、取值方面 ,基本有两种做法 。其一是认为z 随空间坐标 Z 变化 ,如 Homma (1962) 17 、Bakker (1974) 18 基于混合长度假定导出悬沙分 布 ,Vonvisessomjai (1986) 19 对他提出的涡粘系数分布 (Vonvisessomjai ,198420 ) 进行修正得到悬沙扩散系数并导出 悬沙分布 ;另外一种做法是假定紊动与泥沙的交换沿垂向保持不变时 ,z 与 z 无关 ,由上式就可以得到 :¯s = s Kexp ( - z/ LS )(7)式中 : s k 为参考点时均含沙量 (从沙纹峰顶算起 z = 0) ; LS 为泥沙交换垂向尺度
8、, LS =z/ ,这时 ,只要 s k 或LS 和 s k 已知 ,就可得到含沙浓度的分布 。依 Nielson (1992) 15 研究结果 ,泥沙交换垂向尺度 LS 与沙纹高度密切相关 。其表达式为 : am amLS = 0 . 075 , < 18 amLS = 1 . 4,式中 :为波动圆频率 ;为沙纹高度 。 > 18(8)根据 Nielson (1981) 21 的研究结果 ,在实验室规则波条件下 ,沙纹波长 、高度 与沙纹陡度的经验关系式可分别表示为 :/ am = 2 . 2 - 0 . 3450134/ am = 01275 - 010220. 5(9)(10
9、) (11)= U2 / ( s21) gd 其中 :112水槽试验m50对于含沙量的测量 ,本试验采用虹吸采样称重法 。一般来说 ,虹吸方法只能测出时均含沙量分布 ,但根据Bosman (1987) 22 的研究 ,对于波动水流 ,当采样管口垂直于振荡流平面时 ,虹吸抽水的速度大于等于 A ·后 ( A为波动水质点振幅 ,为波动角频率) ,取样得到的泥沙含量与真实泥沙含量的比值就保持常数 ,约为 018 左右 。为了测量含沙量沿垂直方向的分布 ,共采用了 6 根吸管 ,吸管为内径 5 mm 的玻璃管 ,根据水深的不同 ,采用不同的距离固定在支架上 。测量时 ,把吸管连同支架一同放入
10、水中 ,吸管口方向与波浪振荡水流的平面 垂直 。最下面一根吸管的中心距床面 1115 cm (造波之前) 。吸管通过导管与取样瓶相联接 ,抽取的含沙浑 水由吸管口进入 ,经导管流入取样瓶 ,取样瓶为 50 ml 的广口瓶 ,瓶口用橡皮塞塞紧 。塞子联出两根管子 ,一为 6 个吸管提供动力 。缓冲瓶有防止浑水进入真空泵的作用 ,真空泵提供动力并控制采样时间 ,根据对实验条件的初步计算 ,波动水质点最大速度约为 0140 m/ s 。实验时采用了较大功率的真空泵 ,确保真空泵提供的 动力使 6 根吸管的抽水速度大于 0140 m/ s 。因此 ,可以认为实测含沙量与真实含沙量是比较接近的 。实验测
11、量阶段主要为波高的测量和浑水的抽取 ,浑水的抽取是在连续打波时间大于 1 h ,实验段中部悬 沙分布趋于稳定后进行的 。由于本实验是测量时间及空间平均含沙量 ,根据 Van Rijin 等 ( 1993) 23 的工作 , 移动距离有 5 个沙纹波长以上就可获得稳定的平均含沙量 。本实验小车移动距离为 50 cm ,约 10 个沙纹波 长左右 。实验时抽取浑水的时间为 12 min 左右 ,连续作用约 60 个波以上 。后期测量和整理工作主要是测量波浪作用下床面形成的沙纹高度和波长以及获得含沙量 。床面沙纹的高度用测针测量 ,测量精确到 015 mm ;沙纹波长用钢板尺测量 ,精确到 1 mm
12、 ,沙纹高度和波长均连续测量 6 个以上的数据 ,取其平均值而得到 。含沙浓度的获得需要如下步骤 :将抽取出的含沙浑水依次量出其体积 , 用已称重 、编号的滤纸分别进行过滤 ,然后烘干 、称重 ,可得浑水中泥沙的重量 ,根据下式即得最后所需的体 积含沙量 :体积比含沙量 = 泥沙重量/ 沙的比重浑水的体积本次实验共进行了两种粉沙的悬沙测量 ,即进行了维坊泥沙和去除部分粘性颗粒后的黄骅泥沙在波浪 作用下的悬沙实验 。表 1 中列出了所有实验组次的波浪条件和沙纹尺度测量结果 。在进行正式测量之前 , 先进行了两组波浪作用下的悬沙测量 ,以验证实验的可重复性 ,实验条件为表 1 中所列 A 组和 B
13、 组 。图 1 显示了这两组实验测得的悬沙分布结果 ,由图可以看出 ,对于这两组实验条件接近的情况 ,测得的悬沙分布也 是很接近的 。表 1 波浪作用下悬沙测量实验条件及测量结果实验组次水深(cm)周期( s)波高(cm)沙纹波长(cm)沙纹高度(cm)泥沙粒径(mm)skL s(cm)kg/ m3cm3 / cm3518 ×1024616 ×1024410 ×1024119 ×1024716 ×1024816 ×1024415 ×1024617 ×102411111101911311011311011311011
14、38128197118165186107148156167132182192123123123153154173175100125012501200130012501400160016501650170010730107301073010730107301073010360103601036010361155117501990151210221291119117821623137AB C D E F G H IJ2525252516162525161630321917253567816658919 ×10241217 ×1024图 2 分别列出了实验组次 C J 的悬沙测量
15、结果 。分析上述实验结果可以看出 ,当粉沙质床面上形成沙纹 时 ,波浪作用下沙纹上的悬沙分布可以分为两个部分 ,一 是靠近床面含沙量增加很快的高浓度部分 ,二是主水体 中的悬沙部分 。113 理论计算与试验结果的比较取维坊沙粒径 D50为 01073 mm ,静水沉速为318 mm/ s ;黄骅沙粒径为 D50 = 01036 mm ,静水沉速为 0182 mm/ s 。 根据试验结果看 ,作者对 A 、B 组试验适线 , 并做无量纲处理后得 :0111 (ln s ) 2 = z/ L S (开方取负值)(12)s k图 1 A 、B 组次实验测量结果及适线结果比较am表 2 理论计算公式而
16、 L S = ,沙纹高度用式 (4) (10) 计算 ,当 > 150 时 ,组次由式 12 得计算公式= D50 ×100 ,适线结果如图 1 。表 2 是 A - J 组理论计算公式 ,图 2 是 C - J 组计算结果与试 验结果比较图 。可以看出波浪动力强 、掺混强时公式与资料的吻 合好 。掺混越强含沙量的分布趋于均匀 ,反之 ,掺混较差时 ,上下层含沙量差别大 。114现场实测资料验证取 A = 010273 得 ,现场 D50 = 01036 mm , H = 2 m , h = 5 m ,0111 (ln¯s + 7145) 2 = z/ 300111
17、(ln¯s + 7182) 2 = z/ 19140111 (ln¯s + 8156) 2 = z/ 170111 (ln¯s + 7118) 2 = z/ 250111 (ln¯s + 7105) 2 = z/ 350111 (ln¯s + 7171) 2 = z/ 670111 (ln¯s + 7131) 2 = z/ 810111 (ln¯s + 6192) 2 = z/ 660111 (ln¯s + 6167) 2 = z/ 58A 、BC D E F G H IJT = 7 s , sk = 2 . 62
18、 kg/ m3 = 9 . 9 ×1024 cm3/ cm3 , L S = 5 . 7 m 。公式为0111 (ln ¯s + 6192) 2 = z/ 5 . 7 ,如图 3 示 。图 2 实测悬沙浓度分布及计算结果在 z = 0 处 , ¯s = 9 . 9 ×1024 cm3/ cm3 = 2 . 62 kg/ m3 ; z = 014 m 处 , ¯s = 4 . 4×1024 = 1 . 18 kg/ m3 ; z = 410 m 处 , ¯s = 1 ×1024 = 0 . 27 kg/ m3 。 其
19、含沙量分布的趋势与现场状况基本吻合 。2 底部高浓度含沙水体层的研究211底部高浓度含沙水体概念的提出研究表明 ,粉沙质海岸无风天泥沙运动微弱 ,水体清澈 ,而大风天水 体浑浊 ,泥沙运动剧烈 。研究海岸工程时 ,只有研究大风大浪时的泥沙运动才具有现实意义 。 在大风浪作用下粉沙质泥沙的运动形态仍为推移质和悬移质 。对于推移质泥沙有其特点 ,就是在大浪作用下床面泥沙易形成成层运动 。但图 3 现场水体含沙量计算结果此研究此部分泥沙的运动规律有其十分重要的意义 ,并因此提出底部高浓度水体的概念 。212底部高浓度含沙水体的高度由试验结果及曲线描述的情况 ,可以得出高浓度含沙水体的高度为 : am
20、hs = 011 ×L s = 011 (14)当 < 150 时 ,由 10 式计算 ;当 > 150 时 ,= 100 d50 。表 3 为 A - J 组试 验与现场状态下高浓度高浓度含沙水体层计算及实测厚度 (cm)表 3组次A 、BCDEFGHIJ现场1193119321533153617不明显811不明显616不明显518不明显计算实测3457无资料含沙水体层高度的计算值及部分实测值 。将 14 式代入 12 式得 , z = hs 处的含沙量为 :¯s s = 014 s k而由 12 式得出 ,当 z = 0 时 ,即底部最大含沙量为 :
21、5;s m = s k即最大含沙量与高浓度含沙水体上层含沙量之比为 :¯s m = 2 . 5 ¯s s式 15 和 17 说明用含沙量的数值也可以判定高浓度含沙水体的厚度 。(15)(16)(17)213高浓度含沙水体的含沙量与上部水体的含沙量的计算由式 12 求得 :z¯s = e Insks23(18)L上部水体平均含沙量为 :h1h - 0. 1 L s0. 1L¯s d zs u =s1 1332 L seln s k ( -3 1-z 2- z 2hL sL sz e- e=29L s3011L s 1 1011L s3 2 L s·
22、;s k ( 3 z 1 ez 2- z 2 -LL)+ e2ss= 9 ( h - 0 . 1 L s) Lsh2 L s·skh e - 3 h - e - 3h)(0175 - 3(19)= 9 ( h - 0 . 1 L s)L sL sL s下部高浓度水体层平均含沙量为 : 1 1 330 011L s31- z 2L- z 2L 1= 20 s k (5 s k011L s0 z 2 e+ e¯s d =¯s d zss=99L s011L s(20)表 4 为 A - J 组试验和现场状况下的上下层平均含沙量 。由表 4 可以看出 ,波浪动力越强 ,
23、掺混强度越大 ,则上下层含沙量趋于均匀 ,反之 ,上下层差别增大 。需 要特别指出的是 ,由于上层水体的含沙成分与底质不同 ,泥沙颗粒较细 ,沉降速度小 ,另外 ,由于水流作用下1h - 0 . 1L s横向掺混 ,因此 ,上层水体含沙量 s u =¯s u ,为大于 1 的系数 ,由实测资料与计算法比较确定 。表 4 A - J 组试验和现场条件下的上下层平均含沙量及比值 ( kg/ m3 )214近底水质点的运动速度影响近底水质点运动速度的 因素很多 , 包括波浪 、水流运动 、 地形地貌 、海岸工程以及异重流等 。因此 , 在工程中计算近底水 质点的运动速度 ,并加以概化是组次
24、A 、BCDEFGHIJ现场¯s u¯s d¯s d/ ¯s u012401863158011101555100010501285160013711123103015011272154012901662128014701992111017511451193019111872105013611454103 ¯sm/ ¯s u 6144 9100 10108 5145 4157 4110 3180 3147 3169 7125 一项十分复杂的工作 。目前的方法是采用数值模拟技术 ,并结合局部 、典型动力条件下的实测资料来概化 。为了研究问
25、题方便 ,暂考虑受波浪 、潮流动力影响下的理想状态 。在波浪 、潮流作用下 ,悬移质泥沙的运 动规律是“波浪掀沙 、潮流输沙”,在近底部泥沙的输移考虑潮流作用 。文献 10 在研究过程中引用了流速垂线分布的对数分布规律 ,表达为 : uc 1 z= k ln k s + B s()21u 3经过试算 k s = D50 ×104 ,此公式计算结果与现场实测资料对比如图 4 。图 5 是测站分布图和流速矢量图 。计算公式如表 5 。对比表明本公式可以较好地反映现场流速的变化规律 。图 4 某地近海潮流垂线分布 (虚线为实测数 、实线为计算结果)对式 21 积分得 h s 层 (高浓度含
26、沙水体层) 的平均流速 :1011L s ud zud ( t ) =0011 L s1u 3 ( t) ( 1 ln+ B s) d zz011L s=0011 L skk su 3 (t) ( zln - z) + B u ( t ) z011Lz=s 3kk s0 0 . 1 L s= u 3 ( t ) ( 1 ln 1)(22)+ B s -kk sk1011 L s图 5 近海区潮流流速矢量图表 5流速计算公式 (m/ s)站号潮况公式底层水深 (m)全部平均流速底层平均流速比值 uc z 涨潮0. 250. 380. 152. 5301039 = 2. 5ln 0. 36 + 7
27、. 31 #( - 215 m)ucz落潮0. 250. 260. 161. 63= 2. 5ln + 11. 6010190. 36 uc z 涨潮01023 = 2. 5ln 0. 36 + 14. 10. 570. 420. 291. 452 #( - 510 m)ucz落潮0. 570. 370. 251. 48= 2. 5ln + 12. 8010220. 36 uc z 涨潮01023 = 2. 5ln 0. 36 + 14. 00. 570. 440. 291. 513 #( - 810 m)ucz落潮0. 570. 370. 231. 61= 2. 5ln + 11. 7010
28、220. 36 uc z 涨潮01029 = 2. 5ln 0. 36 + 9. 00. 250. 380. 162. 384 #( - 215 m)ucz落潮0. 250. 310. 122. 58= 2. 5ln + 8. 0010260. 36 uc z 涨潮0. 570. 430. 281. 5401027 = 2. 5ln 0. 36 + 11. 95 #( - 510 m)ucz落潮0. 570. 350. 241. 46= 2. 5ln + 12. 8010210. 36 uc z 涨潮0. 570. 420. 241. 7501027 = 2. 5ln 0. 36 + 10.
29、26 #( - 810 m)ucz落潮= 2. 5ln + 11. 30. 570. 360. 221. 64 01022 0 . 36 采用表 3 中的现场资料 h s = 0 . 1 ×L s = 0 . 57 m ,求得图 5 所示各点的底部平均流速 ,见表 5 。由于 ud ( t ) 是时间 t 的函数 ,对于潮汐波可以假定 :ud ( t) = udmcost(23)对 t 在 - 区间积分得 :22ud = 2 udm(24)215底部高浓度含沙水体输沙量计算底部高浓度含沙水体的单宽输沙率 :(25)Q1 = sd·ud·hs其中 s d 由式 (2
30、0) 计算 , ud来源于涨急或落急流速断面 ,可以由实测得到 ,也可以由潮波推算 。参考文献12345678中国水利学会泥沙专业委员会泥沙手册 M . 北京 :中国环境科学出版社 ,1992 ,472101J TJ 213 - 98 ,海港水文规范 S 1赵冲久 . 波浪作用下粉沙质底沙运动特性的试验研究 J . 水道港口 ,1994 (1) :34391 赵冲久 ,等 . 粉沙质海岸泥沙运动特点的实验研究 J . 水道港口 ,2002 (4) :2592611 严恺 ,梁其荀 . 海港工程 M . 北京 :海洋出版社 ,1996 ,1712521曹祖德 ,等 . 波 、流共存时的水体挟沙力
31、 J . 水道港口 12001 ,22 (4) :1511551曹祖德 ,等 . 波 、流共存时的床面剪切力 J . 水道港口 12001 ,20 (2) :56601Liu Jia - ju. Determination of the silt concentration on shoal under the action of wind wave and tidal current , A Proc . of the 2nd COPEDEC C . Sept . 1987. 255 2591曹祖德 ,等 . 波 、流共同作用下的沙起动 R . 天津 :交通部天津水运工程科学研究所 ,200
32、21徐宏明 ,等 . 粉沙质海岸泥沙运移规律的研究 R . 天津 :交通部天津水运工程科学研究所 ,19971 南京水利科学研究院 . 粉沙质海岸建港中的泥沙问题 R . 南京 :南京水利科学研究院 ,19991 赵冲久 ,等 . 天然海况下推移质输沙规律的研究 R . 天津 :交通部天津水运工程科学研究所 ,19991 赵子丹 . 波浪作用下的泥沙起动 J . 海洋学报 ,1983 (1) :19231孙林云 ,等 . 粉沙质海岸建港中的泥沙问题研究 (1) 波浪作用下细沙及粉沙质泥沙运动问题研究 ( 初稿) R . 南京 : 南 科院 ,19961Nielson P. Coastal Bo
33、ttom Boundary Layers and Sediment Transport R . World Scientific . 1992.Nakato T ,Locher F A , Glover J R ,et al . Wave entrainment of sediment from ripple bedsA . Proc . ASCEC . ASCE ,1977. Homma M , Horkawa K. Suspended sediment due to wave actionA . Proc . 8 th Conf . Coastal EngC . ASCE ,1962.Ba
34、kker. W. T. Sand concentration in a oscillatory flowA . Proc . 14 th Int . Conf . Coastal EngC . 1974.Vongvisessomjai S. Profile of suspended sediment due to wave actionJ . J . Wtrwy ,Prot ,Coastal and Oc . Eng. ,1986 ,112 ( ) :8592. Vongvisessomjai S. Oscillatory boundary layer and eddy viscosityJ . J . Hydr. Eng. ,ASCE ,1984 ,110 (4) :7894.Nielson P. Dynamics and geometry of wave generated ripplesJ . Geophys. Res. ,1981 ,86 (c7) :
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 生产安全主题班会
- 内镜室医疗安全提醒
- 铝扣板候车室装修协议
- 前台礼仪规范塑造良好企业形象
- 展览活动吊车租赁合同样本
- 化学品安全异常处理办法
- 石家庄市文化交流中心租赁合同
- 制度评审在制造业的作用
- 市场合作经营协议
- 《人体早期发育》课件
- 金蝶精斗云进销存操作手册
- ug基础培训教材中文基础教程3
- 活动制度汇编15篇
- 《植物生理学》期末考试复习题库(500题)
- 畜禽粪污处理与资源化利用项目可行性研究报告
- GB/T 8170-2008数值修约规则与极限数值的表示和判定
- GB/T 33822-2017纳米磷酸铁锂
- GB/T 3217-1992永磁(硬磁)材料磁性试验方法
- GB/T 2467.2-1996硫铁矿和硫精矿中铅含量的测定第2部分:示波极谱法
- 第二章蚁群算法
- 关于水资源的开发利用
评论
0/150
提交评论