1、土石坝渗流分析，采用非饱和土渗流参数，迭代计算浸润线，根据前次计算结果，不断修改单元的渗透系数和浸润逸出点位置，直到满足精度要求。本算例的土石坝体型比较简单.采用非饱和渗流计算.即渗透系数为空隙压力的函数.首先建立一个数据文件 PPPP.TXT，存储渗透系数函数关系，如下。第一列为空隙压力值（水头 M)，第二列为渗透系数指数，渗透系数等于 10A(M/D)。 -10.00 -4.0E+00 -9.00 -3.6E+00 -8.00 -3.2E+00 -7.00 -2.8E+00 -6.00 -2.4E+00 -5.00 -2.0E+00 -4.00 -1.6E+00 -3.00 -1.2E+0
2、0 -2.00 -8.0E-01 -1.00 -4.0E-01 0.00 0.0E+00 土坝顶宽 4M，上下游坡比均为 1:2，总高 12M，底宽 52M。上游水深 8M，下游无水。 FINI /TITLE, EARTHDAM SEEPAGE /FILNAME,SEEPAGE5 /PLOPTS,DATE,0 *DIM,TPRE,TABLE,11,1,1,PRESS,KKPE ! 定义水压与渗透系数的关系数组 *TREAD,TPRE,PPPP,TXT ! 读入数组 *DIM,NCON,ARRAY,4 ! 定义数组，用于存贮单元四个节点号 /PREP7 SMRT,OFF ANTYPE,STATI
3、C ! THERMAL ANALYSIS ET,1,PLANE55 MP,KXX,1,1 ! 饱和状态下的渗透系数 MP,KXX,2,1E-4 ! 完全干燥下的渗透系数，假设空隙水压力小于-10M 时 K,1,24,12 K,2,24,0 K,3,0,0 K,4,28,12 K,5,28,0 K,6,52,0 L,1,3 L,3,2 L,1,2 L,4,5 L,5,6 L,4,6 LESIZE,ALL,24 A,1,3,2 A,1,2,5,4 A,4,5,6 MSHK,2 ! MAPPED AREA MESH IF POSSIBLE MSHA,0,2D ! USING QUADS AMESH,
4、ALL ! MESH AREAS NUMMRG,NODE ! MERGE NODES AT BOTTOM OF CAISSON *GET,N_MAX,NODE,NUM,MAX ! 获得最大节点号 *GET,E_MAX,ELEM,NUM,MAX ! 获得最大单元号 *DIM,N_TEMP,ARRAY,N_MAX ! 定义节点温度变量-总水头 *DIM,N_PRE,ARRAY,N_MAX ! 定义节点压力水头变量 !定义上游面总水头值 LSEL,S,LINE,1 NSLL,S,1 NSEL,R,LOC,Y,0,8 D,ALL,TEMP,8 !定义上游面总水头值 !定义下游面总水头值 LSEL,S,
5、LINE,6 NSLL,S,1 *GET,N_NUM2,NODE,COUNT *DIM,N_NO2,ARRAY,N_NUM2 II=0 *DO,I,1,N_MAX *IF,NSEL(I),EQ,1,THEN ! 判断节点是否选中 II=II+1 N_NO2(II)=I ! 存储渗流可能逸出点节点编号 *ENDIF *ENDDO NSEL,R,LOC,Y,0,8 ! 第一次计算，假设浸润线逸出点在 8M 高位置,与上游同高 *GET,N_NUM,NODE,COUNT ! 获得渗流出口节点总数 *DIM,N_NO,ARRAY,N_NUM ! 定义变量,存储渗流出口节点编号 II=0 *DO,I,1
6、,N_MAX *IF,NSEL(I),EQ,1,THEN ! 判断节点是否选中 II=II+1 N_NO(II)=I ! 存储渗流出口节点编号 *ENDIF *ENDDO *DO,I,1,N_NUM D,N_NO(I),TEMP,NY(N_NO(I) ! 定义下游面总水头值 *ENDDO ALLSEL,ALL FINISH /SOLU SOLVE FINISH !第一次计算完毕 !- !迭代计算 CONUTT=20 ! 最大循环次数 DD_HEAT=0.001 ! 前后两次计算，总水头最大允许计算差 CHUK_ST=3 ! 出口边界条件重新设定的起始点 CHUK_MAXY2=10E5 ! 临时
7、变量，用于存储浸润线出口坐标 *DO,COM_NUM,1,CONUTT DD_H=0 /POST1 SET,1 *DO,I,1,N_MAX *IF,COM_NUM,GT,CHUK_ST+1,THEN DD1=N_TEMP(I) *IF,ABS(DD1-TEMP(I),GT,DD_H,THEN DD_H=ABS(DD1-TEMP(I) *ENDIF *ENDIF N_TEMP(I)=TEMP(I) ! 计算节点温度（总水头） N_PRE(I)=N_TEMP(I)-NY(I) ! 计算节点压力，总水头-Y 坐标 *ENDDO *IF,COM_NUM,GT,CHUK_ST+1,THEN *IF,DD
8、_H,LE,DD_HEAT,THEN *EXIT *ENDIF *ENDIF /PREP7 ! 重新给每个单元设定材料 MATNUM=2 *DO,I,1,E_MAX *DO,KK,1,4 *GET,NCON(KK),ELEM,I,NODE,KK ! 获取单元四个节点编号 *ENDDO TEMP_Y=(N_TEMP(NCON(1)+N_TEMP(NCON(2)+N_TEMP(NCON(3)+N_TEMP(NCON(4)/4 !计算单元中心点平均温度 RESS_T=TEMP_Y-CENTRY(I) *IF,PRESS_T,GT,0,THEN RESS_T=0 MPCHG,1,I *ELSEIF,P
9、RESS_T,LT,-10,THEN RESS_T=-10 MPCHG,2,I *ELSE MP,KXX,MATNUM+1,10*TPRE(PRESS_T) MPCHG,MATNUM+1,I MATNUM=MATNUM+1 *ENDIF *ENDDO ! 重新设定出口边界条件 *IF,CONUTT,GT,CHUK_ST,THEN !前 CHUK_ST 次采用原边界条件 LSEL,S,LINE,6 NSLL,S,1 DDELE,ALL,TEMP ! 删除原边界条件 II=0 CHUK_MAXY=0 *DO,JJ,1,N_NUM2 *IF,N_TEMP(N_NO2(JJ),GE,NY(N_NO2(
10、JJ),THEN D,N_NO2(JJ),TEMP,NY(N_NO2(JJ) ! 总水头=Y 坐标 *IF,NY(N_NO2(JJ),GT,CHUK_MAXY,THEN CHUK_MAXY=NY(N_NO2(JJ) *ENDIF *ENDIF *ENDDO *IF,CHUK_MAXY2,NE,CHUK_MAXY,THEN ! 判断前后两次计算的浸润线出口位置是否相同 NSEL,R,LOC,Y,CHUK_MAXY ! 选择最高节点 *IF,CHUK_MAXY,GT,0,THEN DDELE,ALL,TEMP ! 删除出口最高节点边界条件 *ENDIF CHUK_MAXY2=CHUK_MAXY *
11、ENDIF *ENDIF ALLSEL,ALL FINI /SOLU SOLVE FINISH *ENDDO SAVE !迭代计算完毕，进入后处理 FINISH /POST1 /CLABEL,1 /EDGE,0 /CONTOUR,8,0,1,8 PLNSOL,TEMP ! 显示总水头云图 PLVECT,TF, , , ,VECT,ELEM,ON,0 PLVECT,TF, , , ,VECT,NODE,ON,0 LSEL,S,LINE,6 NSLL,S,1 PRRSOL,HEAT ! PRINT FLOWRATE THROUGH SOIL FSUM,HEAT ! 计算渗流量 *GET,Q_DAY,FSUM,0,ITEM,HEAT ALLSEL,ALL SAVE *DO,I,1,N_MAX N_TEMP(I)=TEMP(I) ! 计算节点总水头（温度） N_PRE(I)=N_TEMP(I)-NY(I) ! 计算节点压力，总水头-Y 坐标 DNSOL,I,TEMP,N_PRE(I) ! 将压力水头值复制到节点 *ENDDO PLNSOL,TEMP ! 显示压力水头云图 FINI
