Question

计算程序

Answer 1

c name dfc1.for，无界承压含水层非稳定流抽水试验泰斯公式求T、S值

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

c 下一行定义计算泰斯井函数W（u）及自变量u中所用到的变量，需定义为双精度

double precision u，w，wu，y

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

c 用Q（m³/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

b＝0

c＝0

e＝0

f＝0

do 10 i＝1，n

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

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

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

if（u.gt.1）goto 20

c 利用（1-8-3）计算0＜u＜1时的W（u）：

wu＝（（（（0.00107857*u-0.00976004）*u+0.05519968）*u

*-0.24991055）*u+0.99999193）*u-0.57721566-dlog（u）

goto 30

c 以下4行利用（1-8-3）计算u≥1时的W（u）

20 y＝（（（u+8.5733287401）*u+18.059016973）*u+8.6347608925）*u

*+0.2677737343

w＝（（（u+9.5733223454）*u+25.6329561486）*u+21.0996530827）*u

*+3.9584969228

wu＝1/exp（u）/u*（y/w）

c 计算（1-8-1-3）式中的s₀（t）：

30 s10（i）＝q/（4*3.1415926*t0）*wu

c计算

；计算

：

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

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

c 以下5行计算（1-8-2）式中的各求和项

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

c 按（1-8-2）式计算ΔS、ΔT：

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

dt＝（e-b*ds）/a

t0＝t0+dt

s0＝s0+ds

c 检验误差是否满足要求：

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

goto 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