1.1 冲积河渠床面阻力试验研究(1)
王士强
摘要:本节分析总结了在可调坡循环长槽内进行的冲积床面阻力试验资料,同时分析了国内外其他大量水槽试验及天然冲积河渠的床面形态及阻力资料;结合对冲积床面阻力影响因素的物理分析,总结提出了预报冲积河渠床面阻力的规律及公式,同时提出了区分床面形态低能态和高能态的判别公式,初步搞清和提出了过渡区的阻力关系。根据本节提出的公式,可以计算预报不同水沙条件下长江三峡库区淤积河段、长江中下游、黄河下游及各种冲积河渠的床面阻力。
1.1.1 前言
冲积河渠阻力是河渠水力要素及挟沙力计算的基础,是河流防洪及工程设计的重要依据。冲积河流总阻力包括床面阻力、岸壁阻力及河段局部阻力。因大多数天然冲积河段河宽远大于水深,故床面阻力一般占据了总阻力的绝大部分。床面阻力随着水流、床沙条件及相应的床面形态而变化,规律十分复杂。
冲积河渠床面阻力是泥沙运动机理研究的重要方面。自1914年Gilbert进行冲积床面水槽试验以来,国外已积累了大量冲积床面形态及阻力的水槽试验的资料,国内则测验和分析了大量天然河渠的水流阻力资料。国内外对于明槽冲积床面的形态及阻力已进行过大量的研究。在众多的研究中,值得提出已为国内外广泛引用的Engelund的研究成果,他提出了低能态和高能态两种不同的冲积阻力θ*与θ'∗的关系。其中
D为床沙代表粒径,Rb、R'b分别为整个床面的和沙粒的水力半径,S为水流能坡,γ和γs分别为水和泥沙的容重。新近很多学者的研究表明,Engelund的阻力关系尚存在三个问题:①θ'∗不仅随θ'∗而变化,还受其他因子的影响;②过渡区的阻力关系还是空白;③用弗劳德数Fr大于1或小于1区分床面形态高、低能态与实际不符。
Hayashi根据实际资料分析,提出θ∗=f(θ'∗,S)的阻力关系,即θ∗尚与能坡有关。W. R. White等根据实际资料分析确定的阻力关系中,在与θ∗、θ'∗相类似的两个参数,即总剪切可动性和有效剪切可动性之间,无论在低能态或高能态,都另外受到无量纲粒径D∗的影响。
本节针对国内外已有研究中存在的问题安排了实验,着重研究不同水深、比降及粒径在低、高能态和过渡区对床面阻力变化的影响规律。
1.1.2 长槽水流阻力试验
1. 试验设备及方法
本试验在清华大学水利系新建的可调坡循环水槽内进行。水槽总长64m,有效段长60m,测验段长50m,槽宽1.2m,高0.62m,边壁全部为玻璃,水槽内铺床沙厚15cm,为天然泥沙,γs=2.65。床沙中径D50为0.076mm,D65为0.086mm,几何标准差σ=0.5(D84/D50-D50/D16)=1.51,级配曲线如图1所示。水槽纵比降可自动调节。水槽尾部设调节水位的闸门,在尾部加设了一个15cm高的坎。试验中流量Q和水深H为自变量,随着床面形态的自动调整,能坡S和挟沙浓度为从变量。试验运行时间以床面形态发展达到平衡、水流达到均匀流为准,即水流阻力和挟沙能力沿程、历时不变,达到平衡状态。θ'∗越小,达到平衡的调整时间越长,此时间一般为4~8h,长则超过20h。
图1 试验床沙级配曲线
每组试验在平衡状态下测量诸水力要素和床面形态要素(含沙浓度、分布及级配也进行了测验,另文总结)。流量由经过量水堰率定的电磁流量计测定。水槽沿程水面高程由测压管、测针及跟踪式水位计三种方法测定。长50m测验段内沿程两壁共布置了36个测压孔。床面高程由水尺或测淤高程仪测定。因本试验水槽比国外绝大多数水槽长得多,故提高了至关重要的水面比降的测量精度。另外,在长水槽内,每组试验是否达到平衡状态的判断比短水槽要准确。
2. 试验成果及分析
(1)试验成果
本试验共进行了64组,主要水力要素和床面形态要素测验数据及计算成果列于表1。试验流量Q范围为24.2~410L/s,单宽流量为20.2~342L/s·m,水深H范围为6.2~36.5cm,能坡S为0.01×10-3~3.05×10-3,无量纲床面剪切力θ∗为0.012~3.07,弗劳德数Fr为0.073~1.05。
试验中床面形态随着水流强度增大,由静平整发展至沙纹,该临界条件θ∗=θ'∗=0.04。在这种粉细砂质河床,水深只能达到几十厘米的水槽内,床面形态发展过程中不出现沙垄(相同粒径泥沙在大水深的天然河流内出现沙垄)。在水深变化6倍的情况下,最大平均沙纹高度Δ和波长λ并无明显变化,Δ为2.5cm,λ为18~20cm。当θ'∗增大到一定值以后,沙纹开始急剧趋于消失,在θ∗-θ'∗图上出现明显转折,θ∗随θ'∗增长反而减小。这就是沙波从低能态向高能态发展的过渡区。水深越小,这种转折在较小θ'∗时开始,随着水深增大,转折发生时的θ'∗增大。沙纹高度和波长在此过渡区急剧减小,但陡度减小不太多,背水坡角α减小甚为明显。沙纹完全消失达动平衡以后,随着θ'∗的继续增长,床面开始出现逆行沙坡,床面形态极不稳定,相伴的水面波动越来越烈。
表1 60m长槽阻力试验数据(D50=0.076mm)
续表
∗ 1:低能态平整床面;2:沙纹;4:过渡区,急剧减小的沙纹;5:高能态动平整床面;7:逆行沙波
(2)试验成果初步分析
表1中所列θ'∗除D50以外,取决于能坡S及沙粒水力半径R'b,后者由Engelund及很多国外研究者使用式(1)计算确定:
经过对本实验及其他大量各种粒径床沙的平整床面阻力的分析发现,对于D50<0.1~0.11mm的粉细砂情况,取沙粒当量粗糙度K's=0.5D65,对其他情况取K's=D65比较符合实际。粉细砂床沙的沙粒阻力比较小的原因,尚待进一步研究清楚。
表1所列θ∗中包括的床面水力半径Rb由H. A. Einstein的处理方法计算确定。取玻璃边壁糙率nw=0.008。达西-韦斯巴赫阻力系数f=8(V∗/V)2,V∗=,V为平均流速。
将根据试验成果计算的θ∗及θ'∗点绘于图2,且区分不同的Rb/D50情况。由图可见,在同一种床沙粒径时,除过渡区外,θ∗主要取决于θ'∗,与水深或Rb/D50无明显的单独关系,因小水深相应大比降,故与S也无明显的单独关系。在过渡区,不同的Rb/D50相应地存在不同的θ∗-θ'∗关系,θ∗随θ'∗增加而减小;在相同θ'∗时,θ∗随Rb/D50增大而增大。图2内也点绘了黄河下游及Chitale的渠道[1]资料,其床沙粒径与本试验床沙相似,它们的θ∗与θ'∗关系与本试验点据完全一致。在过渡区,天然河流的θ∗随θ'∗增大连续变化而无突变,这是因为天然河流资料的Rb/D50均大于10000,这是水槽试验中随着Rb/D50的增大θ∗转折推后发生,即低能态转为高能态的临界θ'∗不断增大的必然结果。
图2 θ∗与θ'∗、Rb/D50关系(D∗=1.9)
阎颐在16m长可调坡循环水槽内以D50=0.22mm的天然泥沙作床沙进行了床面阻力试验,将其成果,还有Guy等[1]、Barton、林秉南[1]的水槽阻力试验资料及长江中下游的实测资料点绘于图3。这些资料的床沙中径相近。图3中θ∗随θ'∗的变化规律与图2完全一致,在过渡区,对于给定的Rb/D50情况,θ∗也随θ'∗增加而减小;在相同的θ'∗时,Rb/D50越大,θ∗越大。对比图2与图3可以看出,在相同θ'∗时,不同粒径D时θ∗明显不同。图5将各种不同粒径的水槽及野外资料统一点绘,这点尤为显著。
将上述点据以阻力系数f与θ'∗关系点绘于图4。由图可见,f随着θ'∗的增加,先增加后下降,床面处于动平整时f呈极小值,随着θ'∗的继续增加,f又再度增加。在同一粒径及相同Rb/D50时,f-θ'∗呈一连续变化曲线,这与θ∗-θ'∗关系中高、低能态时为二条不连续的曲线是不同的。在同一粒径时,水深越大,阻力系数越小,但当Rb/D50大于10000以后,随着水深或Rb/D50增大,f减小甚微,可以视为不变。随着水深增大,本试验和文献[7]中的沙波相对高度Δ/H及陡度Δ/λ均随之减小,天然河流中沙垄背水坡坡角均比小水深的水槽中的沙垄背水坡角小得多,波峰平顶段又长得多,这些因素都使f值在相同D和θ'∗时随水深增大而呈指数递减。图4只列出了床沙中径为0.076mm及0.2mm左右两种情况,所有其他粒径f随θ'∗及Rb/D50的变化规律都是类似的,只不过资料不够完整。
图3 θ∗与θ'∗、Rb/D50关系(D∗=5)
1.1.3 冲积床面阻力影响因素分析
作者认为,无量纲床面剪切力θ∗及阻力系数f分别包括以下四个主要部分。
式中,θ'∗为无量纲沙粒剪切力,也是决定床面剪切力其他部分的主要因子;θ″∗为沙波阻力,由床面形态产生;θb为运动推移质颗粒对水流形成的阻力;θw为水波破碎及水流紊乱直接形成的阻力;f'、f″、fb及fw则分别为相应的阻力系数。θb和θw并非床面直接造成的对水流的阻力,但因推移质靠近床面、与床面沙粒不断交换且与床沙组成密切相关,而逆行沙波阶段的水波破碎及水流紊乱与沙浪相伴而产生,另外,影响这些阻力的主要因子相同,故宜于将它们归入床面阻力一并分析计算。
f'=8(V'∗/V)2,由式(1)可知,其值取决于相对沙粒粗糙度K's/Rb,在沙纹形成以前的低能态平整床面阶段,f=f'或θ∗=θ'∗。
推移质运动颗粒对水流造成阻力的主要原因是颗粒运动速度ub小于当地水流速度uy,这部分阻力换算成单位床面上的剪切力Fb为
图4 两种床沙f与θ'∗、Rb/D50关系
对于跃移质,根据作者推导,Bagnold的分析及实测资料,(uy-ub)~ω,则水流运动强度较大时有
式中,ω为颗粒沉速;CD为阻力系数常数;N为单位床面上运动的推移质颗粒总数。根据推移质输沙率qb的函数关系,qb~θ'p∗D3/2,ub~θ'∗D1/2,N~qb/(ubD3)~θ'm∗D-2,式中m>0,因θb~Fb/D,故
由上式可见,θb随θ'∗增大而增大。因随粒径D的增大,沉速ω由与D2成正比逐渐改变到与D1/2成正比,故对于D<0.076~0.1mm情况,θb~D3。当雷诺数Re=ωD/ν>103以后,即D约大于3mm以后,θb就与粒径无关。在动平整床面及逆行沙波初期,θ″∗及θw基本为零,θ∗-θ'∗=θb,图5中θ∗-θ'∗随θ'∗及D∗变化的规律与上述分析完全一致。
在逆行沙波阶段,因沙浪形状圆滑呈流线型,故θ″*甚小,此阶段阻力除θ'∗以外,主要包括θb及θw。在此阶段中,水流非常紊乱,因水波波速大于沙浪波速,故逆行沙波一般很不稳定,床面冲淤变化频繁,沙浪时现时隐,水面波涛汹涌,水波时大时小,经常破碎,水面时有巨大漩滚。这些能量损失更多地取决于绝对水流强度ρv'2∗,粒径影响较小,即θw~θ'a∗Db∗(a>0, 0<b<1),故相同θ'∗时D大者θw较大,同时θb较大。此阶段θ∗-θ'∗=θb+θw,故在高能态时相同θ'∗下θ∗随D增大而增大,式中无量纲粒径D∗定义见式(4)。
在低能态沙纹或沙垄阶段,沙波阻力是床面阻力的主要部分。以往常以水流突然扩散模式分析沙波阻力影响因子,从而得出f″~(Δ/λ)(Δ/H)。根据观察分析,作者认为,从物理本质上看,f″主要取决于单位床面上沙波背水坡后面产生的漩涡总体积与水流体积之比,实际情况更接近于渐变扩散,故除了Δ/λ及Δ/H以外,背水坡角α及Δ1/Δ也对f″产生重大影响,其中Δ1为沙波顶高程与平均河床高程之高差。
作者分析了水力学管道中渐变扩散的水头损失系数与扩散角的关系,发现在大部分情况下损失系数与tan α接近成正比,这一分析与Klaasen等在水槽内固定沙波情况下测得的α对阻力系数的影响结果完全一致。过渡区的阻力系数所以迅速减小,除了Δ/λ、Δ/H减小之外,α的减小也是一个显然的原因。分析长江汉口及黄河花园口河段沙波测验资料发现,这些沙垄背水坡角α一般不超过5°,远比水槽试验中相同的θ'∗时的α值为小,这说明在相同D及θ'∗时,水深很大的天然河流比小水深水槽中的f值所以要小的一个重要原因就是背水坡角的调整减小。
进一步分析床面阻力系数和沙波要素的关系看出,Δ1/Δ也是不可忽略的一个重要的调整因素。Δ1/Δ越小,意味着沙波波峰平顶段长度和面积越大,漩涡总体积所占比例越小。分析沙波纵剖面及平面形态照片发现,当沙纹、沙垄从三角形发展为梯形时,阻力系数一般总是减小。长江、黄河等天然河流中的沙垄一般呈梯形,故比水槽中三角形沙垄的阻力系数为小。
综上分析,当水力和床沙条件改变时,沙波自动调整其大小、形状,主要通过Δ/H、Δ/λ、tan α及Δ1/Δ这四个要素的变化,改变床面形态阻力系数。作者以θ∗或f或Ks/D65为纵坐标,以θ'∗为横坐标点绘了大量水槽和天然河渠冲积床面的阻力资料,对比了以S或Rb/D或D作参变数的不同情况。由于水槽试验中细粒径往往伴随着小比降,故若以S为参数[5]确实也得出了相同θ'∗时S越大θ∗越小的大体趋势。但合并分析大量天然冲积河渠资料时就会发现,在θ∗~θ'∗图上,虽然S和Rb/D50相差很大,但只要D相近,则除过渡区外,点群都相当集中。相反,粒径相差很大的资料,即使S相近,但点群仍以D的大小有规律地分散成不同的带群。
1.1.4 冲积河渠床面阻力变化规律及预报
图5以θ∗及θ'∗为纵、横坐标,以D50为参数点绘了D50从0.018mm至28.65mm各种粒径的水槽试验及天然河渠阻力资料近1000组,Rb/D50的范围为10~120000。图中点据区分了高能态和低能态。分析了现有的所有过渡区的资料,过渡区的θ∗-θ'∗关系与第二节的分析完全一致。因图上点据太多,故略去了水槽试验中过渡区的点据以免混淆低、高能态的点据。根据图5大量点据适线确定的水槽及天然河渠冲积床面阻力关系如下。(2)
图5 冲积河渠θ∗与θ'∗关系
在低能态区(对天然河渠,过渡区也包括在内),
式中,x=lg(θ'∗/0.04),
D∗=[g(γs/γ-1)/ν2]1/3 D50;ν为水的运动黏滞系数;当D∗>80以后,θ∗与θ'∗的关系变化甚微,可按D∗=80计算。
在高能态区,
式中,A=1.4/lg(θ'1/0.04);G=1+4.874exp(-0.79D∗);θ'1=0.68+0.32exp(-0.1D∗)。
由图5可见,在低能态,在相同θ'∗时,θ∗随D∗增大而减小,在高能态则反之,在式(5a)中θ∗与θ'∗的差值主要是θb造成的,因为此时θ'∗不太大,θw比较小。式中θ'1与D∗呈指数函数关系,这就反映了实际资料中θ∗随D∗增大而增大的趋势迅速减弱直至与D∗基本无关这一事实,这与式(3b)中D对θb影响的物理分析是完全一致的。
对于天然冲积河流、渠道、水库,一般Rb/D50均大于10000~15000,此时可用式(4)和式(5)分别计算求出二对θ∗、θ'∗值,其中一对较大的θ∗及相应较小的θ'∗值即为需求预报值。这对数值若是据式(4)得出的,即处于低能态,若是据式(5)得出的,即为高能态,不必首先判断属哪种能态。
对于Rb/D50<10000的情况,在θ∗-θ'∗关系图上会出现转折过渡区,如图2、图3所示。为建立计算预报过渡区的θ∗-θ'∗关系,首先需要找到低能态转为高能态即沙纹、沙垄完全消失转为动平整床面的判别式。沙纹、沙垄的消长主要取决于水槽床底流速及粒径,前者又取决于平均流速V和相对水深Rb/D50,故沙粒弗劳德数Frd=V/及Rb/D50是床面形态发展、判别高低能态的主要参数。图6绘出了区分低能态和高能态两类床面形态的大量水槽试验及天然河渠实测资料,后者包括长江中下游、黄河下游、Missouri河及Rio Grande Conveyance Channel等天然河渠床面形态资料。根据图6,当Rb/D50≤10000时,得到如下区分高低能态的判别式,当Rb/D50>10000以后,临界Frd不再随Rb/D50变化。
图6 床面形态高低能态判别
根据式(6)、式(5)及式(1),即可求得给定Rb/D50情况下相应于临界动平整床面的θ'c及相应θc,此点即为过渡区θ∗-θ'∗转折曲线与高能态公式(5)曲线的交点。Rb/D50越小,则此点沿式(5)曲线越向左下方移动(图6上)。根据已有的水槽试验资料,在求得临界θ'c及θc以后,可初步按式(7)确定过渡区的阻力。
式中,
图7点绘了Rb/D50>10000的天然河渠的阻力系数f与θ'∗的关系。在相同θ'∗时,在低能态,f随D∗增大而减小,对高能态则反之。对于D50>0.15mm的沙质河床,f随D减小甚微。根据图7实测点据适线确定式(8)。图中线6为式(8)计算值而无实测点据,但该种粒径的Rb/D50较小的水槽试验点据说明线6也是合理的。
式中,x=lg(θ'∗/0.04);
图7 天然河渠f与θ'∗关系
式(8)包括了天然河渠高、低能态及过渡区的阻力关系。由图7可见,在f-θ'∗图上过渡区并不明显,只是在低能态的后半段f随θ'∗增长减少较快。
式(8)与式(4)、式(5)是基本等价的,但对于D50>0.7mm或Rb/D50<10000的情况,式(8)未曾包括,可按式(4)~式(7)计算。对于天然沙质河渠,利用式(8)预报阻力系数极其简便,适用范围已相当大。
众所周知,黄河下游河床糙率比一般河流的糙率小。其实按阻力系数规律,若床沙较细的黄河与床沙较粗的长江床面形态均处于低能态具有相同θ'∗时,黄河下游的f值反而大于长江中下游。但由于以下几点原因使黄河下游糙率普遍较小。
(1)黄河下游沙细、比降陡,故流量不大时θ'∗就比较大,使大多数情况床面处于动平整左右。黄河下游沙垄一般都处于消失过渡区。沙垄波长λ很大,背水坡角α很小,这样就使处于低能态的f值比床沙较它为粗的河流大得不多。对于较大的洪水,当θ'∗>1.3以后,由图7可见,在相同θ'∗时沙细的黄河下游的f值均比床沙较粗的河流小得多。
(2)黄河下游沙细、冲泻质多,沙粒阻力系数f'较小,沙细使fb也较小,这就使动平整左右的f值特别小,因为此时的f值主要只包括f'和fb。
(3)在水文测验和工程设计中普遍采用糙率n值来表示阻力系数,n与f值的关系为n=R1/6×(f/8g)1/2,糙率是有量纲的,除与f1/2成正比外,尚与R1/6成正比。黄河下游与长江中下游相同θ'∗时,前者的水深一般只为后者的1/5,这就使两者f相同时,前者的n值只是后者的约3/4。因此高能态床面形态情况下黄河下游的n值比长江的显得小得多,而低能态情况下,f值稍高的黄河下游,其n值反显得比长江的还小。即使在θ'∗=0.2左右两者f值差别最大的阶段,黄河下游的糙率也不比长江为大。上述分析尚未包括测验高含沙浓度及冲淤影响。
1.1.5 结语
(1)低能态时冲积床面阻力主要随床面形态而变化。沙波形态阻力主要取决于沙波背水坡后的漩涡的相对体积。除沙波相对高度及陡度以外,沙波背水坡坡角及沙波形状对形态阻力也有重要影响。高能态时阻力主要因推移质运动滞后于水流及水流漩涡紊乱、水波破碎而造成。
(2)影响床面阻力的参数主要有θ'∗、D∗及Rb/D。达西-韦斯巴赫阻力系数f(或f″)在相同D∗及Rb/D情况下,随θ'∗增大而先增后减再增大,呈一连续变化曲线。对于每一种给定粒径的床沙来说,在相同θ'∗时,Rb/D越大,f值越小,但当Rb/D>10000以后,f值变化甚微,故对于一般天然河渠来说,f主要只取决于θ'∗及D∗。可据式(8)预报天然沙质河床f值。对于天然河渠或Rb/D相同时,在θ'∗相同时,低能态时f随D∗增大而减小,在高能态时反之,但D50>0.15mm以后的沙质河渠,不同粒径的f-θ'∗曲线相差甚小。
(3)在θ∗-θ'∗阻力关系图上,对于Rb/D50>10000的天然河渠来说,θ∗为θ'∗及D∗的函数,θ∗随θ'∗增长而增长。当θ'∗固定时,θ∗随着粒径增大,在低能态时减小,在高能态时增大,但这种增大趋势呈指数迅速递减而趋于一定值。可据式(4)、式(5)计算预报Rb/D50>10000的所有天然冲积河渠的床面剪切力θ∗值。
(4)在低能态向高能态转变时,沙纹或沙垄消失减弱较快,在θ∗-θ'∗关系图上呈现一个明显的转折过渡区。在此区固定Rb/D时,就有一个确定的θ∗-θ'∗关系,Rb/D越大,转折越晚发生。对于给定的Rb/D,θ∗随θ'∗增长而减小。当Rb/D<10000时,可按高、低能态的临界值判别式(6)求得给定Rb/D50及D50时的临界流速Vc,再据式(5)、式(1)确定θ'c及相应的θc值,按式(7)预报过渡区的θ∗及θ'∗值。其他区域的阻力关系仍按式(4)、式(5)预报。
参考文献
[1] Brownlie William R. Compilation of alluvial channel data[R]. California Institute of Technology, 1981.
[2] IRTCES. Sediment transport data in laboratorial flumes[R]. circular No. 2, Beijing, China. 1987.
[3] Engelund F. Hydraulic resistance of alluvial stremns[J]. Proc. ASCE, JHD, Vol. 92, No. HY2, 1966.
[4] Engelund F. Hydraulic resistance of alluvial sireams (closure to discussion) [J]. Proc. ASCE, JHD, Vol. 93, No. HY4, 1967.
[5] Hayashi Taizo. Alluvial bed forms and roughness [R]. NSF Sediment Research Workshop, San Francisco, USA. 1986.
[6] White W R, et al. Frictional characteristics of alluvial streams in the lower and upper regimes [J]. Proc. Instn. Civ. Engrs, part2, Vol. 83. Dec. 1987.
[7] 阎颐(指导教师:王士强).冲积床面形态及阻力水槽试验研究[D].清华大学硕士学位论文,1989.
[8] Wang Shiqiang, Zhang Ren. A New equation of bedload transport [C]. Proc. 22 Congress IAHR, 1987.
[9] Klaassen G J, et al, DHL—Research on bedforms, resistance to flow and sediment transport[C]. Proc. 3rd ISRS, USA. 1986.
Experimental study on hydraulic resistance of alluvial streams
Wang Shiqiang
Abstract: The investigation of hydraulic resistance of alluvial streams was conducted in a 60-meter recirculating flume at the Sediment Research Laboratory of Tsinghua University in Beijing. A great number of laboratory test results and field data from other investigators have also been analyzed in this paper.
The result of study shows that the dimensionless resistance parameter θ∗ depends not only on θ'∗ but also on the diameter of the bed material in lower and upper regimes.
In transition zone from ripple or dune to flat bed, θ∗ decreases with the increase of θ'∗ for a given D and Rb/D and it increases with the increase of Rb/D for given D and θ'∗.
The friction factor f depends on θ'∗, D and R/D. f changes with θ'∗ continuously for given D and Rb/D and decreases with the increase of Rb/D but is maintained constant when Rb/D is larger than 10000. Friction factor f decreases in lower regime and increases in upper regime with the increase of D for given θ'∗.
θ∗ and θ'∗ can be computed by the formulae (1), (4)~(7) suggested in this paper. The value of f on sand bed for natural streams can be predicted by formula (8).