土石坝渗流

土石坝渗流分析，采用非饱和土渗流参数，迭代计算浸润线，根据前次计算结果，

不断修改单元的渗透系数和浸润逸出点位置，直到满足精度要求。

本算例的土石坝体型比较简单.采用非饱和渗流计算.即渗透系数为空隙压力的函数.

首先建立一个数据文件PPPP.TXT，存储渗透系数函数关系，如下。

第一列为空隙压力值（水头M)，第二列为渗透系数指数，渗透系数等于10^A(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+00

-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,STATIC ! a static analysis

ET,1,PLANE55 !2D温度场单元

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,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,LINE,,6

NSLL,S,1

*GET,N_NUM2,NODE,,COUNT !获得下游所有结点总数，保存在N_NUM2中

*DIM,N_NO2,ARRAY,N_NUM2 !定义数组N_NO2，大小为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,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 ! 临时变量，用于存储浸润线出口坐标

*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_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 !计算单元中心点平均温度

PRESS_T=TEMP_Y-CENTRY(I)

*IF,PRESS_T,GT,0,THEN

PRESS_T=0

MPCHG,1,I

*ELSEIF,PRESS_T,LT,-10,THEN

PRESS_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(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

*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