计算程序

c name taisi-TS.for，无界承压含水层非稳定流抽水试验泰斯公式用收敛级数求解，计算T、S值

c计算公式：s＝Q/（4πT）W（u）u＝Srr/（4Tt）

c 求W（u）的子程序改为级数，下一行泰斯井函数W（u）及自变量u，需定义为双精度

double precision u，wu

dimension t（500），s1（500），s10（500），ss（500），st（500）

c 用Q（m3/d）的定流量抽水，在距主孔r（m）处观测孔中观测水位降深值，s0：待求参数S赋初值，t0：待求参数T赋初值

open（1，file＝’qr431.txt’）

read（1，*）q，r，s0，t0

close（1）

s00＝s0

t00＝t0

open（1，file＝’st431.txt’）

cn：水位降深观测次数，读进观测孔中观测到的水位降深值n对：

read（1，*）n

read（1，*）（t（i），s1（i），i＝1，n）

close（1）

50 a＝0.0

b＝0.0

c＝0.0

e＝0.0

f＝0.0

do 10 i＝1，n

c 下一行计算t（i）时刻的u值：

u＝s0*r*r/t0/t（i）/4.0

write（*，*）’u＝’，u

call taisi（u，wu）

s10（i）＝q/（4.0*3.1415926*t0）*wu

ss（i）＝-q*exp（-u）/s0/t0/4/3.1415926

st（i）＝-s10（i）/t0+q*exp（-u）/t0/t0/4.0/3.1415926

a＝a+st（i）**2

b＝b+st（i）*ss（i）

c＝c+ss（i）**2

e＝e+st（i）*（s1（i）-s10（i））

f＝f+ss（i）*（s1（i）-s10（i））

10 continue

ds＝（a*f-b*e）/（a*c-b**2）

dt＝（e-b*ds）/a

t0＝t0+dt

s0＝s0+ds

if（abs（dt/t0）.lt.0.001）goto 40goto 50

40 write（*，*）’t＝’，t0，’（m2/d）s＝’，s0

open（1，file＝’ts431.txt’）

write（1，*）’Q（m3/d）r（m）T0（m2/d）S0’，’T（m2/d）S’

write（1，60）q，r，t00，s00，t0，s0

close（1）

60 format（2f10.2，f10.2，f10.5，f10.2，f10.5）

open（1，file＝’st0431.txt’）

do 80 i＝1，n，3

write（1，90）t（i），s1（i），t（i+1），s1（i+1），t（i+2），s1（i+2）

80 continue

close（1）

90 format（3（f10.4，f7.3，9x））

stop

end

subroutine taisi（u，w）

double precision u，w

dimension ni（100）

c 泰斯井函数值W为指数积分，用收敛级数前10项之和表示：

w＝0

ni（1）＝1

do 30m＝2，10

30 ni（m）＝ni（m-1）*m

do 60 n＝1，10

60 w＝w+（-1）**（n+1）*u**n/（n*ni（n））

w＝w-0.577216-dlog（u）

10 continue

if（w.lt.0）w＝0

return

end

一、动态载荷的添加

通用程序控制中的外载荷为恒定载荷，对于材料在动态载荷作用下的破坏过程，所受的外载荷一般都是随着时间变化的载荷，例如三角波载荷等。通过对通用程序控制进行改进添加冲击三角波载荷，以适应对冲击问题的模拟。三角波载荷的示意图如图13-1所示。

图13-1 三角波载荷示意图

其中：Fmax为动态三角波冲击载荷的峰值；t1和t2分别为动态载荷的峰值时间和截止时间。

二、模型计算流程

岩石是一种典型的d脆性材料，在外力作用下，一般不会发生明显的塑性变形，因此可以不考虑材料的塑性变形。由于岩石的抗压强度远大于其抗拉强度，可以认为岩石在受压时，材料内部不会产生新的损伤，即不考虑体积压缩时的损伤累积。尽管在采用静态d性模型计算时，可以不考虑材料的损伤演化，但是由于在其他时步的计算过程中，产生的损伤是不会恢复的，所以在进行d性计算时，材料的d性模量将不同于最初的d性模量，变为损伤后的d性模量。根据以上的思路，可按照以下的计算流程来进行计算：

1.总体计算流程（针对单个块体）

（1）若块体单元损伤变量D为1，则不再对该单元进行计算而直接进入下一个单元的计算，否则进入（2）。

（2）计算体积应变εv，如果εv＜0即单元处于体积压缩状态，则进入（3），否则进入（4）。

（3）按照d性损伤模型，由应变计算应力，并且保持D值不变（不考虑体积压缩时的损伤积累），进入（5）。

（4）当岩石处于受拉状态时，材料的应力应变关系并不都是按照应力波衰减模型来计算的，而且还与应变变化率有一定的关系，仅当应变率大于10时，才能按照动态的应力波衰减模型，由应变计算应力值，进入（5）；而当应变率小于10时，按照岩石受压状态下的d性损伤模型来计算（高尔新等，1999），然后进入（5）。

（5）如果D＞1或D=1，则令D=1，并置当前所计算块体的压力和偏应力为零，进入下一步的计算。

在每一步的计算中需要对所有单元进行一一计算，并保存当前步的应力及应变计算结果以供下一步计算使用。

2.具体的计算模型和公式

第一步计算：

（1）由于对D值的计算需要用到D的变化率，因此本步计算中保持D值不变。

（2）应变变化率

岩石断裂与损伤

（3）体积应变

岩石断裂与损伤

（4）动态应力波衰减模型

σij=3K（1-D）εδij+2G（1-D）eij

即

岩石断裂与损伤

（5）d性损伤模型

岩石断裂与损伤

在以上各式中：E、G、K、μ分别为材料的d性模量、剪切模量、体积模量和泊松比。

第二步及其以后的计算：

（1）D值的计算：

岩石断裂与损伤

（2）应变变化率的计算：

岩石断裂与损伤

（3）体积应变的计算：

岩石断裂与损伤

应力波衰减模型与d性损伤模型同第一步相同。

以上的讨论和建立的公式都是针对于平面应力问题而言的，如果要求解的问题属于平面应变问题，需要把以上方程式中的E换为，μ换为。

---