抽水蓄能電站水利樞紐工程邊坡穩(wěn)定分析

抽水蓄能電站水利樞紐工程邊坡穩(wěn)定分析1．1邊坡穩(wěn)定分析原理鋼筋混凝土面板堆石壩的穩(wěn)定分析根據(jù)工程經(jīng)驗(yàn)及設(shè)計(jì)理論,應(yīng)進(jìn)行施工期和蓄水發(fā)電期穩(wěn)定分析,對于本工程而言,由于壩體施工采用的是分層碾壓法,整個壩高在施工期都在不斷變化,而且壩體受的汽車荷載、施工機(jī)器荷載和人群荷載，具體考慮起來比較復(fù)雜，并且數(shù)量級相對較小。要準(zhǔn)確分析整個施工期壩體穩(wěn)定性比較困難。所以，在這里只計(jì)算竣工期的壩體穩(wěn)定，此時期的壩體穩(wěn)定分析應(yīng)分為上游坡、下游坡兩種情況，但對于兩種情況進(jìn)行比較可知，由于上游坡有鋼筋混凝土面板的作用，其整體性要好于干砌塊石及次堆石區(qū)組成的下游坡，再加上上游坡庫盆需用棄石碴回填至115m高程，由于這部分土的作用，上游坡的穩(wěn)定性也要好于下游坡。另外，對于鋼筋混凝土面板壩其防滲性能很好，其面板后的浸潤線很低，這對上、下游坡的影響都不大，不是控制性工況。而對蓄水期，在施工質(zhì)量得到保證的前提下，壩體的上游坡由于水的作用，其穩(wěn)定性會增大，對于下游坡的穩(wěn)定性變化會很小。根據(jù)上述分析：上游坡的穩(wěn)定性要強(qiáng)于下游坡，所以這里只對竣工期下游坡的穩(wěn)定性進(jìn)行分析。目前對于土石壩邊坡穩(wěn)定性分析，幾乎無例外的都是采用近似的方法——極限平衡法來進(jìn)行分析。常用的土石壩穩(wěn)定分析方法有：圓弧法和折線法。圓弧法：包括瑞典條分法和畢肖普法，多適用于粘土壩中，適合土壩。折線法：多適用于非粘土壩中，如薄心墻壩，斜墻壩壩體中。本工程為鋼筋混凝土面板堆石壩，故宜采用折線法。在進(jìn)行堆石體壩坡穩(wěn)定計(jì)算之前，要進(jìn)行大量的現(xiàn)場碾壓試驗(yàn)和室內(nèi)三軸試驗(yàn)，測得堆石體各區(qū)的容重和內(nèi)摩擦角。1.正常情況：圖4-3邊坡穩(wěn)定計(jì)算受力簡圖如圖所示說明折線滑動法計(jì)算壩坡安全系數(shù)的原理。假設(shè)ABC為可能的滑動面其中AB為下滑動面，BC為上滑動面，BD為錯動破裂面。由試驗(yàn)可知，當(dāng)滑動面BC上的摩擦角（）小于保護(hù)層土料的內(nèi)摩擦角（）時，BD面斜向上部，反之，則BD面可能斜向下部。設(shè)滑裂面BD與鉛直的夾角為0，滑動土體ABDA與BCDB之間的相互作用力為P，P與BD面法線的交角為（）。由滑動土體BCDB的平衡可知：當(dāng)時，當(dāng)時，從以上兩式聯(lián)立求解，可得式4-1同理，從土體ABDA的平衡得：式4-2利用式4-1與式4-2相等的條件，可得式4-3計(jì)算時，可先假設(shè)安全系數(shù)Kc,然后將滑動面上由試驗(yàn)得到的土體抗剪強(qiáng)度指標(biāo)c1,c2,，，除以Kc,即，，，，將求得的c1,c2,，，入式（4-3），若等式左右兩邊相等，則假設(shè)的Kc就是該滑動面的穩(wěn)定安全系數(shù)。如不相等，則需另行假設(shè)Kc，值，重新計(jì)算，直到方程兩邊相等為止。2.特殊情況下（即有水平地震力時）采用擬靜力法：從滑動土體BCED力的平衡可得：式4-4從滑動土體ABDA力的平衡可得：式4-5聯(lián)立兩式可得：式4-6式中：F1,F2分別為作用在第1，2兩塊上的水平地震力計(jì)算時，可先假設(shè)安全系數(shù)Kc,然后將滑動面上由試驗(yàn)得到的土體抗剪強(qiáng)度指標(biāo)c1,c2,，，除以Kc,即，，，，將求得的c1,c2,，，入式（4-6），若等式左右兩邊相等，則假設(shè)的Kc就是該滑動面的穩(wěn)定安全系數(shù)。如不相等，則需另行假設(shè)Kc，值，重新計(jì)算，直到方程兩邊相等為止。當(dāng)40<h,地震烈度為7度時，地震水平加速度分布系數(shù)如下：圖4-4地震水平加速度分布系數(shù)圖為了計(jì)算方便，并且有安全保障，可采用的地震水平加速度分布系數(shù)如下：圖4-5簡化的地震水平加速度分布系數(shù)圖分別計(jì)算出兩種工況下各種滑動面的安全系數(shù)，從中取出各工況下的最小的安全系數(shù)。最小的安全系數(shù)均要滿足各工況下的安全要求。５.2邊坡穩(wěn)定分析計(jì)算程序計(jì)算流程圖（一次循環(huán)）開始開始輸入壩體的數(shù)據(jù)輸入壩體的數(shù)據(jù)取取C點(diǎn)坐標(biāo)(x1,y1)取取B點(diǎn)坐標(biāo)(x2,y2)否計(jì)算否計(jì)算BC斜率K1AB斜率K2k1>k2k1>k2k2<1/m是是計(jì)算計(jì)算1,2兩塊重力G1,G2安全系數(shù)安全系數(shù)K=1否否K=K+0.01K=K+0.01滿足聯(lián)立方程滿足聯(lián)立方程k2<1/m是是輸出輸出K及各點(diǎn)坐標(biāo)結(jié)束結(jié)束應(yīng)用fortran語言編程如下：PROGRAMMAIN*********變量說明*********T-壩頂寬度，DELT---A的步長，DELB---B在Y方向的步長，DELA---B在X方向的步長*********A（50）存放A，B（50，500，2）存放B，K（50，500）存放K，*********ZXIA0K（50）存放每個A對應(yīng)下的最小K，BK（50，2）用來存放ZXIAOK對應(yīng)下B的座標(biāo)，*********（CX，CY）--壩的最低點(diǎn)的座標(biāo)H1，H2--壩頂點(diǎn)及壩下游最低點(diǎn)所對應(yīng)的高程****H3---壩底的控制高程***R1，R2---壩體材料的容重MOJIAO1，MOJIAO2--主堆石及次堆石的內(nèi)摩擦角***C----堆石的粘聚力GFK（10）每一分塊的重量***XFK（10），YFK（10）---分塊界點(diǎn)FKJIAO（10）--每分塊的傾角***FKMO（10）---每分塊的內(nèi)摩擦角LENGTH（10）--每分塊的底邊長度***AA（4），BB（4）-四條直線的的系數(shù)***Y1，Y2，Y3，Y4---通過B點(diǎn)的垂直線與四條直線的交點(diǎn)的縱座標(biāo)***N---分塊數(shù)JSK---用來計(jì)數(shù)***I，J---循環(huán)變量M---用來表示A的個數(shù)*** PO----下游坡度***KH，CZ，AI是用來表示地震力系數(shù)REALT,T1,DELT,DELB,DELA REALA(50),B(50,500,2),K(50,500) REALZXIAOK(50),BK(50,2) REALXFK(0:10),YFK(0:10) REALLENGTH(10),FKMO(10),FKJIAO(10),GFK(10) REALCX,CY,H1,H2,H3,PO REALR1,R2,MOJIAO1,MOJIAO2,C REALAA(4),BB(4) REALY1,Y2,Y3,Y4 REALKH,CZ,AI INTEGERI,J,M,JSK,N,KK***********************************************************如果KK等于0，則表示不考慮地震力，如果KK等于1，則表示考慮地震力。kk=1 IF(KK.EQ.0)THEN OPEN(0,FILE='W1.DAT')OPEN(1,FILE='W11.OUT') OPEN(2,FILE='W12.OUT') ENDIF IF(KK.EQ.1)THENOPEN(0,FILE='W2.DAT')OPEN(1,FILE='W21.OUT') OPEN(2,FILE='W22.OUT') ENDIF***************************************************** READ(0,*)T,DELT,DELA,DELB M=INT(T/DELT) T1=M*DELT IF(T.GT.T1)M=M+2 IF(T.EQ.T1)M=M+1 READ(0,*)(A(I),I=1,M) READ(0,*)H1,H2,H3 READ(0,*)R1,R2,MOJIAO1,MOJIAO2,C doi=1,2 READ(0,*)AA(I),BB(I) enddo READ(0,*)KH,CZ,AI************************************************************* PO=1.0/AA(1) CY=H1-H3 CX=PO*CY+T AA(3)=(H2-H3)/CX BB(3)=H1-H2 WRITE(*,*)(A(I),I=1,M),CX,CY MOJIAO1=MOJIAO1*3.1415926/180 MOJIAO2=MOJIAO2*3.1415926/180****************************************************** DO10I=1,M AA(4)=CY/(CX-A(I)) BB(4)=-A(I)*AA(4) JSK=1 B(I,JSK,1)=A(I)110B(I,JSK,2)=-DELB B(I,JSK,1)=B(I,JSK,1)+DELA IF(B(I,JSK,1).GE.CX)THEN JSK=JSK-1WRITE(*,*)I,JSK GOTO140 ENDIFB(I,JSK,2)=B(I,JSK,2)+DELB130Y1=AA(1)*B(I,JSK,1)+BB(1) Y2=AA(2)*B(I,JSK,1)+BB(2) Y3=AA(3)*B(I,JSK,1)+BB(3) Y4=AA(4)*B(I,JSK,1)+BB(4)100IF(B(I,JSK,2).LT.Y4)THEN B(I,JSK,2)=B(I,JSK,2)+DELB GOTO100 ELSE GOTO120 ENDIF120IF(B(I,JSK,2).GT.Y3)THEN GOTO110 ENDIF*************************************************************在FENKUAI子程序中，利用壩頂?shù)狞c(diǎn)A、折線與主次堆石分界線的交點(diǎn)（XJIAO，YJIAO）******折點(diǎn)B、壩腳點(diǎn)C將滑動體分塊，分別求出每分塊點(diǎn)的X、Y座標(biāo)、每分塊體底的長度、****每分塊體底的內(nèi)摩擦角及與水平線的夾角 CALLFENKUAI(B(I,JSK,1),B(I,JSK,2),A(I),CX,CY,N,AA,$BB,T,Y2,MOJIAO1,MOJIAO2,FKMO,FKJIAO,LENGTH,XFK,YFK)****在QIUZHOHGLI子程序中求出了每分塊體的重力 CALLQIUZHONGLI(GFK,XFK,YFK,R1,R2,B(I,JSK,1),B(I,JSK,2),$Y2,T,AA,BB,N,A(I))*****在QIUK子程序中求出了折點(diǎn)為B時的安全系數(shù) CALLQIUK(K(I,JSK),GFK,N,FKMO,FKJIAO,LENGTH,C,JSK,I,$B(I,JSK,1),B(I,JSK,2),KH,CZ,AI)************************************************************ JSK=JSK+1 B(I,JSK,1)=B(I,JSK-1,1) B(I,JSK,2)=B(I,JSK-1,2)+DELB GOTO130*****************************************************************此段用來求出在A點(diǎn)不動、B點(diǎn)循環(huán)所的的最小安全系數(shù)140 ZXIAOK(I)=K(I,1) BK(I,1)=B(I,1,1) BK(I,2)=B(I,1,2) DO20J=2,JSKIF(K(I,J).LT.ZXIAOK(I))THEN ZXIAOK(I)=K(I,J) BK(I,1)=B(I,J,1) BK(I,2)=B(I,J,2) ENDIF20CONTINUE WRITE(1,*)ZXIAOK(I),A(I),BK(I,1),BK(I,2)*******************************************************8810 CONTINUE CLOSE(1) CLOSE(2) CLOSE(3) END***********************************************************************分塊子程序*****變量說明（與主程序?qū)?yīng)的變量不說明）****（XJIAO，YJIAO）----折線與主次堆石分界線的交點(diǎn)座標(biāo)，DELX、DELY分塊點(diǎn)的步長****I-----循環(huán)變量 SUBROUTINEFENKUAI(X,Y,A,CX,CY,N,AA,BB,T,Y2,MOJIAO1,$MOJIAO2,FKMO,FKJIAO,LENGTH,XFK,YFK) REALX,Y,A,CX,CY,AA(3),BB(3),T,Y2,MOJIAO1,MOJIAO2REALFKMO(10),FKJIAO(10),LENGTH(10),XFK(0:10),YFK(0:10) REALXJIAO,YJIAO,DELX,DELY INTEGERN,I********************************************************************************************* IF(Y.GT.Y2.AND.A.LE.T)THEN CALLQIUZUOBIAO(XJIAO,YJIAO,FKMO,FKJIAO,N,A,X,Y$,CX,CY,MOJIAO1,MOJIAO2,Y2,T,AA,BB) XFK(0)=A YFK(0)=0 XFK(1)=X YFK(1)=Y XFK(2)=XJIAO YFK(2)=YJIAO DELX=(CX-XJIAO)/(N-2) DELY=(CY-YJIAO)/(N-2) DO210I=3,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY210 CONTINUE ENDIF********************************************************************* IF(Y.EQ.Y2.AND.A.LE.T)THEN CALLQIUZUOBIAO(XJIAO,YJIAO,FKMO,FKJIAO,N,A,X,Y$,CX,CY,MOJIAO1,MOJIAO2,Y2,T,AA,BB) XFK(0)=A YFK(0)=0 XFK(1)=X YFK(1)=Y DELX=(CX-XJIAO)/(N-1) DELY=(CY-YJIAO)/(N-1) DO220I=2,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY220CONTINUE ENDIF********************************************************************** IF(Y.LT.Y2.AND.A.LT.T)THEN CALLQIUZUOBIAO(XJIAO,YJIAO,FKMO,FKJIAO,N,A,X,Y$,CX,CY,MOJIAO1,MOJIAO2,Y2,T,AA,BB) XFK(0)=A YFK(0)=0 XFK(1)=XJIAO YFK(1)=YJIAO************************** IF(X.LE.2*XJIAO/3.0+CX/3.0)THEN XFK(2)=X YFK(2)=Y DELX=(CX-X)/(N-2) DELY=(CY-Y)/(N-2) DO230I=3,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY230CONTINUE ENDIF************************** IF(CX/3.0+2*XJIAO/3.0.LT.X.AND.X.LT.2*CX/3.0+XJIAO/3.0)THENDELX=(X-XJIAO)/3.0 DELY=(Y-YJIAO)/3.0 DO240I=2,N-3 XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY240 CONTINUEDELX=(CX-X)/3.0DELY=(CY-Y)/3.0DO250I=N-2,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY250 CONTINUE ENDIF*************************** IF(X.GE.2*CX/3.0+XJIAO/3.0)THEN DELX=(X-XJIAO)/4.0 DELY=(Y-YJIAO)/4.0 DO260I=2,N-1 XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY260 CONTINUE XFK(N)=CX YFK(N)=CY ENDIF***************************** ENDIF****************************************************************** IF(Y.LT.Y2.AND.A.EQ.T)THEN CALLQIUZUOBIAO(XJIAO,YJIAO,FKMO,FKJIAO,N,A,X,Y$,CX,CY,MOJIAO1,MOJIAO2,Y2,T,AA,BB) XFK(0)=A YFK(0)=0 IF(X.LE.2*XJIAO/3.0+CX/3.0)THEN XFK(1)=X YFK(1)=Y DELX=(CX-X)/(N-1) DELY=(CY-Y)/(N-1) DO270I=2,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY270CONTINUE ENDIF************************** IF(CX/3.0+2*XJIAO/3.0.LT.X.AND.X.LT.2*CX/3.0+XJIAO/3.0)THENDELX=(X-XJIAO)/3.0 DELY=(Y-YJIAO)/3.0 DO280I=1,N-3 XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY280 CONTINUEDELX=(CX-X)/3.0 DELY=(CY-Y)/3.0DO290I=N-2,N XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY290 CONTINUE ENDIF*************************** IF(X.GE.2*CX/3.0+XJIAO/3.0)THEN DELX=(X-XJIAO)/4.0 DELY=(Y-YJIAO)/4.0 DO200I=1,N-1 XFK(I)=XFK(I-1)+DELX YFK(I)=YFK(I-1)+DELY200 CONTINUE XFK(N)=CX YFK(N)=CY ENDIF******************************* ENDIF DO201I=1,N LENGTH(I)=SQRT((XFK(I)-XFK(I-1))**2+(YFK(I)-YFK(I-1))**2)201 CONTINUE END****************************************************************************************** 求交點(diǎn)，求每分塊的傾角，求每分塊的內(nèi)摩擦角*****變量說明與FENKUAI子程序所對應(yīng)的變量不說明*****JIAO1、JIAO2分別用來代表折線與水平線的交角*****DELX1、DELX2、DELY1、DELY2分別用來代表A點(diǎn)、B點(diǎn)、C點(diǎn)之間的X、Y座標(biāo)的差值*****DAIA、DAIX、DAICX分別表示用二分法求（XJIAO、YJIAO）時的中間變量 SUBROUTINEQIUZUOBIAO(XJIAO,YJIAO,FKMO,FKJIAO,N,A,X,Y$,CX,CY,MOJIAO1,MOJIAO2,Y2,T,AA,BB) REALXJIAO,YJIAO,FKMO(10),FKJIAO(10),A,X,YREALCX,CY,MOJIAO1,MOJIAO2,Y2,T,AA(3),BB(3) REALJIAO1,JIAO2,DELX1,DELX2,DELY1,DELY2,YY1,YY2 REALDAIA,DAIX,DAICX INTEGERN******************************************************************************* DELX1=X-A DELX2=CX-X DELY1=Y DELY2=CY-Y JIAO1=ATAN(DELY1/DELX1) JIAO2=ATAN(DELY2/DELX2) DAIA=A DAIX=X DAICX=CX****************************************************************************** IF(Y.GT.Y2.AND.A.LE.T)THEN N=6 FKMO(1)=MOJIAO1 FKMO(2)=MOJIAO1 FKJIAO(1)=JIAO1 FKJIAO(2)=JIAO2 DO310I=3,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2310 CONTINUE302 XJIAO=(DAICX+DAIX)/2.0 YY1=AA(2)*XJIAO+BB(2) YY2=TAN(JIAO2)*(XJIAO-CX)+CY IF((YY1-YY2).GE.0.001)THEN DAICX=XJIAO GOTO302 ENDIF IF((YY1-YY2).LE.-0.001)THEN DAIX=XJIAO GOTO302 ENDIF YJIAO=YY1 ENDIF*************************************************************** IF(Y.EQ.Y2.AND.A.LE.T)THEN N=5 XJIAO=X XJIAO=Y IF(A.EQ.T)FKMO(1)=(MOJIAO1+MOJIAO2)/2 IF(A.LT.T)FKMO(1)=MOJIAO1 FKJIAO(1)=JIAO1 DO320I=2,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2320 CONTINUE ENDIF*********************************************************** IF(Y.LT.Y2.AND.A.LT.T)THEN301 XJIAO=(DAIA+DAIX)/2.0 YY1=AA(2)*XJIAO+BB(2) YY2=TAN(JIAO1)*(XJIAO-X)+Y IF((YY1-YY2).GT.0.01)THEN DAIX=XJIAO GOTO301 ENDIF IF((YY1-YY2).LT.-0.01)THEN DAIA=XJIAO GOTO301 ENDIF YJIAO=YY1FKMO(1)=MOJIAO1FKJIAO(1)=JIAO1******************************* IF(X.LE.2*XJIAO/3.0+CX/3.0)THEN N=6 FKMO(2)=MOJIAO2 FKJIAO(2)=JIAO1DO330I=3,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2330 CONTINUE ENDIF**********************************IF(CX/3.0+2*XJIAO/3.0.LT.X.AND.X.LT.2*CX/3.0+XJIAO/3.0)THEN N=7 DO340I=2,N-3 FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO1340 CONTINUE DO350I=N-2,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2350CONTINUE ENDIF******************************* IF(X.GE.2*CX/3.0+XJIAO/3.0)THEN N=6 DO360I=2,N-1 FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO1360CONTINUE FKMO(N)=MOJIAO2 FKJIAO(N)=JIAO2 ENDIF******************************* ENDIF***************************************************************** IF(Y.LT.Y2.AND.A.EQ.T)THEN XJIAO=A YJIAO=0************************************ IF(X.LE.2*XJIAO/3.0+CX/3.0)THEN N=5 FKMO(1)=MOJIAO2 FKJIAO(1)=JIAO1DO370I=2,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2370 CONTINUE ENDIF***********************************IF(CX/3.0+2*XJIAO/3.0.LT.X.AND.X.LT.2*CX/3.0+XJIAO/3.0)THEN N=6 DO380I=1,N-3 FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO1380 CONTINUE DO390I=N-2,N FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO2390CONTINUE ENDIF********************************* IF(X.GE.2*CX/3.0+XJIAO/3.0)THEN N=5DO300I=1,N-1 FKMO(I)=MOJIAO2 FKJIAO(I)=JIAO1300CONTINUE FKMO(N)=MOJIAO2 FKJIAO(N)=JIAO2 ENDIF***************** ENDIF END*******************************************************************************求重力的子程序******S是求分塊重力時的分塊的面積，I是中間變量 SUBROUTINEQIUZHONGLI(G,XFK,YFK,R1,R2,X,Y,$Y2,T,AA,BB,N,A) REALXFK(0:10),YFK(0:10),G(10) REALR1,R2,X,Y,Y2,T,AA(3),BB(3),A INTEGERN,I REALS1,S11,S12,S21,S22,S23************************************************************************** IF(Y.GT.Y2.AND.A.LE.T)THEN IF(X.LT.T)THEN S1=(XFK(1)-XFK(0))*(YFK(1)-YFK(0))/2.0 S21=(YFK(1)+YFK(2))*(XFK(2)-XFK(1))/2.0-$(XFK(2)-T)*YFK(2)/2.0 S22=(XFK(2)-T)*YFK(2)/2.0-$(XFK(2)-T)*(AA(1)*XFK(2)+BB(1))/2.0 G(1)=R1*S1G(2)=R1*S21+R2*S22 ELSE S11=1/2.0*(XFK(1)-XFK(0))*(YFK(1)-YFK(0))-$1/2.0*(XFK(1)-T)*(AA(2)*XFK(1)+BB(2)) S12=1/2.0*(XFK(1)-T)*(AA(2)*XFK(1)+BB(2))-$1/2.0*(XFK(1)-T)*(AA(1)*XFK(1)+BB(1)) S21=1/2.0*(XFK(2)-XFK(1))*(YFK(1)+YFK(2))S22=1/2.0*(YFK(2)+AA(2)*XFK(1)+BB(2))*(XFK(2)-XFK(1)) S23=1/2.0*(AA(1)*XFK(2)+BB(1)+AA(1)*XFK(1)+BB(1))*$(XFK(2)-XFK(1)) G(1)=R1*S11+R2*S12 G(2)=R1*(S21-S22)+R2*(S22-S23)ENDIFDO410I=3,N G(I)=-1*R2*(XFK(I)-XFK(I-1))*(AA(1)*$(XFK(I)+XFK(I-1))/2+BB(1)-(YFK(I)+YFK(I-1))/2.0)410 CONTINUE ENDIF*******************888888*****************************8 IF(Y.EQ.Y2.AND.A.LE.T)THEN S11=1/2.0*(T-A)*YFK(1) S12=1/2.0*(XFK(1)-T)*(YFK(1)-AA(1)*XFK(1)-BB(1))G(1)=R1*S11+R2*S12 DO420I=2,N G(I)=-1*R2*(XFK(I)-XFK(I-1))*(AA(1)*$(XFK(I)+XFK(I-1))/2+BB(1)-(YFK(I)+YFK(I-1))/2.0)420 CONTINUE ENDIF******************************************************** IF(Y.LT.Y2.AND.A.LT.T)THEN S11=1/2.0*(T-XFK(0))*(YFK(1)-YFK(0)) S12=1/2.0*(XFK(1)-T)*(YFK(1)-YFK(0))-$1/2.0*(XFK(1)-T)*(AA(1)*XFK(1)+BB(1)-YFK(0)) G(1)=R1*S11+R2*S12 DO430I=2,N G(I)=-R2*(XFK(I)-XFK(I-1))*(AA(1)*$(XFK(I)+XFK(I-1))/2+BB(1)-(YFK(I)+YFK(I-1))/2.0)430 CONTINUE ENDIF******************************************************** IF(Y.LT.Y2.AND.A.EQ.T)THEN DO4