河流模拟课程设计—水库一维泥沙淤积计算

1、水库一维泥沙 淤积计算课程设计目录一、目的与要求1二、根本原理11、根本方程12、方程离散13、公式补充2三、计算步骤2四、计算框图3五、计算结果41、历年输沙量特征值 42、各年淤积总量53、各年水位库容关系 64、水面线的变化 75、深泓变化86、坝前断面变化 9六、结果分析121、剖面形态分析 122、库容损失合理性分析 12七、 计算程序13、目的与要求具备一定并编制出完整通过课程设计，初步掌握一维数学模型建立数学模型的根本过程和计算方法, 的解决实际问题的能力。以水流、泥沙方程为根底，构建恒定流条件下的河道一维水沙数学模型, 的计算程序，并以某个水库为实例，进行水库泥沙淤积计算。水流

2、条件：恒定非均匀流。泥沙条件：包括悬移质，推移质的均匀沙模型，推移质计算模式为饱和输沙，悬移质计 算模式为不饱和输沙，水流泥沙方程采用非耦合解。二、根本原理1、根本方程水流连续方程:AQx0QQ2txAQQ2txAtSAcxy°Gb水流运动方程泥沙连续方程河床变形方程xgA0 thgAgA i。 ifx2 2Q n c g4-0AR3(S S*)(S S)推移质平衡输沙方程G=G水流挟沙力公式采用张瑞瑾公式，推移质输沙率公式采用Mayer-_Peter公式,MAYER-PETER公式中的能坡J按均匀流曼宁公式近似计算每个断面不同。2、方程离散方程在恒定流情况下有q 0，离散为：Q=c

3、onstx方程变形为£_ QgA x AQ2n2a2r"3Q2x 2gA2 2Q n24A R3上式离散为1 Lj22gQj2A2Zj 1 Zj2xn(AjQj4R/3j 1方程S4去掉时间项得到Sj i S j i SjS*jexpS*jSji-1 exp xBQx对2号断面以下，上式可以离散为:由方程4-5 可得QS' y00B 0BG bBGb QS QSxx'Bjyoj(10)B j 1 y0j 1t一S S q该方程的解析解为：对于进口断面，推移质不考虑，悬移质采用单点离散方程5可离散为：y。(S SJ t103、公式补充S3mk gRK取0.12

4、4 , m取 1.05，干密度0 取 1.3恢复饱和系数0.25均匀沙粒径为d=0.041mm 悬移质，d=2 mm 推移质二、计算步骤1、输入河床地形糙率等数据求得断面面积与水位的关系A Z,进而求得断面平均流速Q，水力学半径A2、读入一个时段的水沙数据特别注意，不要一次性将数据全部读入读入第一时段Q,S值3、计算水面线，同时得到各断面的水力要素求得各个断面的河宽、断面面积、水深、平均流速等值计算前要注意在坝前输入水位，各断面均应对流量赋值4、计算悬移质水流挟沙力gRK取 0.124，m取 1.05。5、计算推移质输沙率采用 mayer-peter公式3/2n 3/2hJ 0.047 s d

5、 ngb0.125 1/2gs6、计算各断面含沙量公式Sj1 S*j1 Sj Sj exp 2S*j Sj1expq x7、计算各断面冲淤厚度进口断面y。一S1 SJ t08、修改水各断面水下河床高程H H。 丫9、重新进入2进行下一循环10、计算10年河床变形，计算时段为一天，单位为秒 s11、淤积总量年输出一次，其余每两年输出一次计算结果四、计算框图读入地形资料和糙率输出初始库容、深泓等YQ=0计算S和Gb计算水流子程序计算S(j)和dy(j)1r4i=i+11F读入一个判断库尾冲刷与否YJvnpnxt重新计算S(j)修改河床地形五、计算结果1、历年输沙量特征值1。其中推移质输沙量是按照计

6、算各年悬疑质输沙量、推移质输沙量、年均流量见表悬疑质输沙量的1.5%计算的。表1龙开口坝址历年流量输沙量特征值代表年年平均悬移质年平均推移质年平均流量输沙量万t输沙量万t万t1984294144.1213301985412061.8017501986322148.3213601987482872.4217601988254338.1514401989642696.3919001990438565.7818001991597289.5818701992198229.7312901993645696.841900取大值645696.841900最小值198229.731290多年平均值428764

7、.3116402、各年淤积总量表2历年淤积总量的变化年悬移质输沙 量万吨年均流量万吨总淤积量亿 立方米000011821.494845560.1425291.20510994180.387316993.5419658070.879420027.3424697751.085521668.3327656171.204623853.9830875691.368730155.0537460621.701834879.5943121571.987936754.1846403382.1171039289.6650197622.243S 2 5 1 S 甲2 100 2 4 6 8 10图1历年淤积总量的变

8、化3、各年水位库容关系表3各年水位库容关系水位m0年2年4年6年8年10年2732.5742.2071.5721.3130.7650.5692712.3852.0181.3841.1320.6160.4382692.2021.8351.2020.980.4930.3332672.0241.6611.0440.8590.3980.262651.8521.510.9280.7730.3460.232631.6881.3720.830.6980.3050.2042611.5321.2450.7450.6280.2690.1792591.3841.1260.6680.5630.2370.1562571

9、.2471.0140.5940.5010.2070.1332551.1170.9070.5230.4420.1790.1122530.9950.8080.4570.3870.1540.0922510.880.7170.3990.3340.1310.0742490.7710.6320.3440.2870.110.0592470.670.5530.2930.2440.0890.0452450.5760.4780.2440.2040.0710.0342430.4920.4080.1990.1680.0530.0252410.4170.3470.1590.1370.0380.01727U2US2602

10、0245图2各年水位库容关系曲线4、水面线的变化表4水面线变化年距坝里程水位年距坝里程水位年距坝里程水位(m)(m)(m)(m)(m)(m)034920267.6434920267.88834920274.43033810267.16433810267.7683381027401431350267.53831350273.43029280267429280267.37829280272.93027000267427000267.21827000272.4024050267424050267.07824050271.5021846267421846267.0282184

11、6270.53019896267419896267819896269.58018270267418270267818270268.92015790267415790267815790268.1013170267413170267813170267.41010170267410170267810170267.05078202674782026787820267049502674495026784950267021002674210026782100267002674026780267年距坝里程水位年距坝里程水位年距坝里程水位(m)(m)(m)(m)(m)(m)234920267.14634920

12、272.011034920273.99233810267.12633810271.751033810273.78231350267.06631350271.011031350273.49229280267.02629280270.371029280273.18227000267627000269.631027000272.79224050267624050268.641024050272.22221846267621846267.941021846271.71219896267619896267.391019896271.07218270267618270267.111018270270.48

13、2157902676157902671015790269.72131702676131702671013170268.942101702676101702671010170267.882782026767820267107820267.19249502676495026710495026722100267621002671021002672026760267100267图3水面线变化曲线5、深泓变化表5逐年深泓变化年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)年距坝里程(m)深泓(m)034920264.1434920264.77834920271.75033810261.7433810

14、264.12833810270.46031350253.3431350262.29831350268.95029280249.4429280261.92829280268.6027000246.6427000263.23827000269.21024050241.3424050261.99824050267.71021846239.2421846261.5821846267.87019896237.9419896261.498198962681418270252.83818270266.84015790230.1415790250.82815790264.620131

15、70228.2413170240.94813170264.74010170224.4410170239.22810170262.5907820222.147820230.9887820253.7604950219.844950228.2984950237.2802100215.9242100219.9782100233.9900213.1740221.0480224.22年距坝里程深泓年距坝里程深泓年距坝里程深泓(m)(m)(m)(m)(m)(m)234920262.79634920269.961034920271.76233810262.41633810269.081033810270.93

16、231350259.92631350267.251031350269.37229280259.53629280266.71029280269.62227000259.48627000267.211027000270.09224050250.62624050265.391024050268.98221846247.11621846265.191021846269.24219896242941019896269.7218270242.84618270264.711018270268.61215790233.11615790258.711015790267.03213170

17、232.41613170249.521013170266.88210170226.49610170242.18101017026527820224.1167820234.07107820264.7924950220.6564950228.12104950245.6322100216.8662100222.71102100236.5420213.7860219.7100229.24280 r260220T0年潇泓T-丄年琛叛 亠4年擁涨召年济亂8毎儼泓迪年探泓1Q年水厨线200 11010QD0200005000040000图4深泓变化曲线6、坝前断面变化该水库为三角形淤积，最靠近坝前断面淤积并

18、不会特别严重，为使断面变化情况表现的更加直观，这里采用IP-17断面作为典型断面并分析其淤积情况。表6坝前断面变化年起点距高程年起点距高程年起点距高程(m)(m)(m)(m)(m)(m)00300.220300.240300.2010.5272.29210.5272.29410.5272.29015.8266.93215.8267415.8267.01018.2260.24218.2261.18418.2264.29018.9259.62218.9260.56418.9263.67019.4259.88219.4260.82419.4263.93024.8259.91224.8260.8542

19、4.8263.96026.4259.82226.4260.76426.4263.87044.2246.15244.2247.09444.2250.2061.2234.47261.2235.41461.2238.52078.6218.17278.6219.11478.6222.22091.5215.92291.5216.86491.5219.97097.3216.1297.3217.04497.3220.150106.8216.982106.8217.924106.8221.030121.6217.472121.6218.414121.6221.520132.4218.052132.4218.9

20、94132.4222.10145.6218.452145.6219.394145.6222.50173.8218.962173.8219.94173.8223.010191.7219.432191.7220.374191.7223.480203.6230.782203.6231.724203.6234.830216.5242.882216.5243.824216.5246.930228.7255.042228.7255.984228.7259.090241.1267.442241.1267.444241.1267.440253.7279.812253.7279.814253.7279.8102

21、65.3291.22265.3291.24265.3291.20279.1304.522279.1304.524279.1304.520293.1300.32293.1300.34293.1300.3年起点距高程年起点距高程年起点距高程(m)(m)(m)(m)(m)(m)60300.280300.2100300.2610.5272.29810.5272.291010.5272.29615.8267.01815.8267.011015.8267.01618.2267818.22671018.2267618.9266.41818.92671018.9267619.4266.67819.426710

22、19.4267624.8266.7824.82671024.8267626.4266.61826.42671026.4267644.2252.94844.2264.221044.2267661.2241.26861.2252.541061.2255.09678.6224.96878.6236.241078.6238.79691.5222.71891.5233.991091.5236.54697.3222.89897.3234.171097.3236.726106.8223.778106.8235.0510106.8237.66121.6224.268121.6235.5410121.6238.

23、096132.4224.848132.4236.1210132.4238.676145.6225.248145.6236.5210145.6239.076173.8225.758173.8237.0310173.8239.586191.7226.228191.7237.510191.7240.056203.6237.578203.6248.8510203.6251.46216.5249.678216.5260.9510216.5263.56228.7261.838228.726710228.72676241.1267.448241.1267.4410241.1267.446253.7279.8

24、18253.7279.8110253.7279.816265.3291.28265.3291.210265.3291.26279.1304.528279.1304.5210279.1304.526293.1300.38293.1300.310293.1300.3图5坝前断面变化曲线六、结果分析1、剖面形态分析A、纵断面淤积情况分析：由深泓变化曲线图可以看出，随着时间的增加，河道深泓点高程成增高趋势。该水库运行10年后，在距坝约10000米处，深泓点高程增加 40多米，计算结果说明，该水库泥沙以 淤积为主，致使水库库容逐年损失，且情况较为严重。B坝前断面变化曲线分析：由坝前断面变化曲线图可以看出

25、， 随着时间的增加，坝前断面各点的总体高程也成增高 趋势。该水库运行10后，该断面河床抬高 20多米，说明在该处河道逐年以淤积为主，且淤 积情况较为严重。2、库容损失合理性分析由水位库容曲线可以看出，水库在运行10年后，经计算可知，在273m水位条件下，库容的损失率为77.89%，在241m的水位条件下，库容的损失率高达95.92%，计算结果说明，经过十年的运行，水库库容根本成全部淤满状态，所以在10年后将不能再继续使用。结果说明，该水库的调节能力较差。上述计算结果与现实不太相符,其原因是在计算该水库库容的泥沙损失时,计算冲往下游的泥沙量不合理。在给定条件下，该水库并没有考虑水位的变化及坝前冲

26、沙情况。该水库的计算条件为正常蓄水位保持不变,在高水位运行条件下， 必然会导致水库的严重淤积， 而实际情况是，在洪水来临时，坝前水位要降至防洪限制水位，以增大防洪库容，该条件下,水流的速度及水流挟沙力较大，会造成库区泥沙处于冲刷状态，以增加水库库容。七、计算程序program mainparameter nn=16,mm=60, ndisp=0, npxt=16,ND=3653,NY=10,T=15)paramete (rs=26.5,r=10.0,d=0.000041,failev=0.5,faidy=0.5,dt=86400)dime nsion x(mm,2 ,nn ),rough( n

27、n), npoint(nn ),B( nn ),A (nn ),Gb( nn ),dy (nn ),Nday(NY) dime nsion dx( nn ),Xw( nn ),dxa( nn ),alow( nn ),alowO( nn ),Q( nn ),Sx( nn ),S( nn ),Zlevel( nn)data(Nday(m),m=1,10)/366,731,1096,1461,1827,2192,2557,2922,3288,3653 /I*翻开输出文件*open(10,file='河床地形.txt',status='old')open(11,fil

28、e='流量沙量.txt',status='old')open(12,file='逐年累计淤积量.txt',status='unknown')open(13,file='逐年深泓变化.txt',status='unknown')open(14,file='逐年水面线变化.txt',status='unknown')open(15,file='逐年坝前断面淤积变化.txt',status='unknown')open(16,file=

29、9;逐年水位库容关系变化.txt',status='unknown')读入地形数*call Qbed(m m,n pxt,x,dxa,alow0 ,npoint)算沉降速*W=FW(T,r,rs,d,0)write(*,*)W算断面间住巨 *do i=1, npxt-1dx(i)=dxa(i)-dxa(i+1)en ddo!*计算初始水面线、深泓、库容水位关系*dWQST=0.0dWQS=0.0dWGb=0.0AVQ=0.0do l=1,NPXT rough(l)=0.035 en ddoQ(1)=295.do I=2,NPXTQ(I)=Q(I-1)en ddoZlev

30、el( npxt)=267.0CALL level(x,rough,nn,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) !初始水面线输出到文件 *CALLFILE(12,13,14,15,16,0,x ,n poi nt,Zlevel,alow,dxa,dx,mm, nn ,dWQST,AVQ,dWQS,dWGb, ndisp)!*计算各个时段水库的淤积情况*DO 40,I=1,NDWRITE(*,*)'读入天数:',I!*读入流量、沙量数据*READ(11,*)Q(1),S(1)IF(Q(1)=0)THENGOTO 41ENDIF

31、DO J=2,NPXTQ(J)=Q(1)ENDDOI*CALL level(x,rough, nn ,zlevel,dx,Q, npoi nt,b,a,xw, nn ,mm,failev,0)!*计算各断面的水流挟沙力和推移质输沙率*CALL SGB( nn ,a,b,Q,Sx,Gb,rough,w, ndisp)!*计算各断面的含沙量和冲淤厚度*CALL DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, npoi nt,faidy, ndisp)I *验质量守丿恒 *CALL CHECK(dy,B,S,Gb,dt, npxt,Zlevel,x

32、, npoi nt,Q,dx, ndisp)41CONTINUE!* 计算悬移质、推移质年均流量*WQS 仁Q(1)*SGb(1)=WQS1*0.015WQSN=Q( nn 广S( nn)dWQST=dWQST+WQS1*dt/10000000 !悬移质累计输沙量dWQS=dWQS+(WQS1-WQSN)*dt/1300.0/100000000!悬移质累计冲淤量dWGb=dWGb+(Gb(1)-Gb( nn )*dt/1300.0/100000000!推移质累计淤积量AVQ=AVQ+Q(1)*dt/10000!累计流量!*输出相关文件*DO m=1,NYIF(l=Nday(m)thenCALL

33、FILE(12,13,14,15,16,m,x ,n poi nt,Zlevel,alow,dxa,dx,mm, nn ,dWQST,AVQ,dWQS,dWGb ,n disp)ENDIFENDDO40 CONTINUEclose(10)close(11)close(12)close(13)close(14)close(15)close(16) endI *读入地形数据*subrouti ne Qbed(m m,nn, x,dxa,alowO ,npoint)dime nsion x(mm,2, nn ),dxa( nn ),alowO( nn), npoi nt( nn)do i=1, nn

34、read(10,*)read(10,*) dxa(i), npoi nt(i)read(10,*)do j=1, npoin t(i)read(10,*) x(j,1,i),x(j,2,i)en ddo计算深泓*alow0(i)=x(1,2,i)do m=2, npoin t(i)if(alow0(i)>x(m,2,i)thenalow0(i)=x(m,2,i)en difen ddoen ddoend!*泥沙沉降速度的计算，采用张瑞瑾公式*FUNCTION FW(T,gama,gamas,d, ndisp)IF(NDIS 圧Q.-1)WRITE(*,*)'INT0 FWcall

35、 VISCOS(T,CMU)A仁CMU*13.95/dA2=1.09*9.8*d*(gamas-gama)/gamaFW=(A1*2.0+A2)*0.5-A1RETURNENDSUBROUTINE VISCOS(T,CMU)X=1.775E-06A=1+0.0337*T+0.000221*T*TCMU=X/ARETURNEND水面线及各断面水力要素计算函数*subroutine level(x,rough,npxt,zlevel,dx,q,npoint,b,a,xw,nn,mm,failev,NDISP) dime nsion x(mm,2, nn ),rough( npxt),dx( npx

36、t),zlevel( npxt) dime nsion q(n pxt), npoint(n pxt),b( npxt),a( npxt),xw (n pxt) IF(NDISP=-1)WRITE(*,*)'INTO LEVEL' nc=1000dz=0.5dz1= 0.1317dz2=2.079call area( npoint(n pxt),x(1,1, npxt),x(1,2, npxt),b (n pxt),a (n pxt),xw (n pxt),zlevel( npxt), ndisp) do ip=n pxt-1,1,-1NR=0zmin=zlevel(ip+1)

37、+dz30call area( npoin t(ip),x(1,1,ip),x(1,2,ip),b(ip),a(ip),xw (n pxt),zm in,n disp)if(A(i p)<=0)the nzmin=zmin+0.173goto 30en diffmin=flevel(zlevel(ip+1),zm in, dx(ip),q(ip+1),q(ip),rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),failev, ndisp)if(fmin> O)the nfr=(q(ip)/a(ip)*2.0*b(ip)/(9.8*a(ip)if(fr&l

38、t;1)the nzmax=zm in+dz计算断面要素*20callarea( npoin t(ip),X(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(ip),zmax ,n disp)fmax=flevel(zlevel(ip+1),zmax,dx(ip),q(ip+1),q(ip),rough(ip),b(ip+1),b(ip),a(ip+1),a(ip),failev ,ndisp)if(fmax>0)the nzmax=zmax+dzgoto 20en difelsezmin=zmin+dz2NR=NR+1IF(NR>NC)THENzmin=zmin

39、-dz1write(*,*)'the loop is death,pause'WRITE(*,*)' nr,ip,zmi n,fmin', nr,ip,zmi n, fminread(*,*)en difgoto 30en difelsezmin=zmin-dz1goto 30en difcallbisec(zm in, zmax,fm in, zlevel(ip), npoin t(ip),X(1,1,ip),x(1,2,ip),b(ip),a(ip),xw(ip),q(ip), ndisp,zlevel(ip+1),b(ip+1),a(ip+1),q(ip+

40、1),dx(ip),rough(ip),failev)en ddoreturnendFUNCTION flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev, ndisp)IF(NDISPEQ.-1)WRITE(*,*)'INTO FLEVEL'HLOWER=ALOWER/BLOWERH=A/BAA=ZLELO-ZLEVELB仁failev*(QLOWER/BLOWER)*2.0/HLOWER*(10/ 3.0)B2=(1-failev)*(Q/B)*2.0/H*(10/3.0)BB=DX*ROUGH*2.0

41、*(B1+B2)C仁(QLOWER/ALOWER)*2.0C2=(Q/A)*2.0CC=(C1-C2)/(2*9.8)FLEVEL=AA+BB+CCRETURNEND计算断面要素*SUBROUTINE AREA(NPOINT,X,Z,B,A,XW,ZLEVEL,NDISP)DIMENSION X(NPOINT),Z(NPOINT)IF(NDISP<=-1)WRITE(*,*)'INTO AREA'IF(NPOINT<=2)THENWRITE(*,*)'IN AREA THE INPUT DATA ARE FALSE N=',Npoi ntSTOPEN

42、DIFIF(Z(1)<ZLEVEL.OR.Z(NPOINT)<ZLEVEL)THENWRITE(*,*)'IN AREA THE WATER LEVEL IS TOO HIGH OR THE HEIGHT OF THE SECTIONAT EDGE IS TOO LOWSTOPENDIFA=0.0B=0.0XW=0.0DO l=1,NPOINT-1ZMIN=AMIN1(Z(I),Z(I+1)IF(ZMIN<ZLEVEL)THENZMAX=AMAX1(Z(I),Z(I+1)IF(ZMAX<ZLEVEL)THENDB=X(I+1)-X (I)DH=ZLEVEL-0.

43、5*(Z(l)+Z(l+1)DX=(Z(l)-Z(l+1)*2.0+DB*DB)*0.5ELSEDB=(ZLEVEL-ZMIN)/(ZMAX-ZMIN)*(X(I+1)-X(I)DH=0.5*(ZLEVEL-ZMIN)DX=(2*DH)*2.0+DB*DB)*0.5ENDIFIF(DB<0.0)THENWRITE(*,*)'THE DISTANCE FOR NODE ',I,'AND',I+1,'ARE FALSE'WRITE(*,*)'IN AREA',I,I+1,X(I),X(I+1)STOPENDIFDA=DB*DHB

44、=B+DBA=A+DAXW=XW+DXENDIFENDDORETURNEND二分法计算*subrouti nebisec(zmi n,zmax,fmi n,zlevel, npoi nt,x ,z, b,a,xw,q, ndisp,zlelo,blower,alower,qlower,dx,rough,failev) dime nsion x(npoin t),z( npoint)IF(NDISFEQ.-1)WRITE(*,*)'INTO BISEC'ERR=0.001 10 ddz=zmax-zm inzlevel=0.5*(zm in+zmax)call area( npo

45、in t,x,z,b,a,xw,zlevel, ndisp)f=flevel(zlelo,zlevel,dx,qlower,q,rough,blower,b,alower,a,failev, ndisp) if(ddz>err)the nff=f*fminif(ff>0)the nzmin=zlevelelsezmax=zlevelen dif ddz=zmax-zm in goto 10en difreturnend!悬移质挟沙力和推移质输沙率计算，悬移质采用张瑞谨挟沙力公式唯移质采用mayer-peter公式SUBROUTINE SGB( nn ,a,b,Q,Sx,Gb,rou

46、gh,w, ndisp)DIMENSION A( nn ),B( nn ),Q( nn ),Sx( nn ),Gb( nn ),rough( nn)IF(n disp=-1)WRITE(*,*)'INTO SGB'H=A(I)/B(I)U=Q(I)/A(I)VA仁 U*3.0VA2=9.8*H*WVA3=(VA1VA2)*1.05Sx(l)=VA3*0.124!悬移质挟沙力!推移质输沙率，断面输沙率Gb(l)=FGb(rough(l),10.0,H,Q(l),A(l),26.5,0.002,1.0,2.65,0)*B(i)ENDDOENDFUNCTION FGb(rn,r,H,

47、Q,A,rs,d,rou,rous, ndisp)IF(ndisp=-1)WRITE(*,*)'INTO FGb'A1=(rn*2.0)*(Q*2.0)A2=(A*2.0)*(H*(4 /3.0)RJ=A1A2rnn=(d*(1 /6.0)/26.0AA=(rn n/rn)*1.5)*rs*H*RJBB=0.047*(rs-r)*dCC=(rous-rou)/rous)*9.8*(rou*0.5)*0.125IF(AA-BB)<O)THENFGB=0ELSEFGB=(AA-BB)*1.5)/CCENDIFRETURNEND!*断面冲淤厚度计算*SUBROUTINE DEPO( nn ,mm,w,x,S,Sx,Gb,dt,B,Q,dx,dy,alow,alowO, npoi nt,faidy, ndisp)DIMENSIONx(mm,2 ,nn), npoint(nn ),B( nn ),Gb( n