s080020478.frm

用matlab语言编写的试井分析的程序

VERSION 5.00
Object = "{F9043C88-F6F2-101A-A3C9-08002B2F49FB}#1.2#0"; "comdlg32.ocx"
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3060
   ClientLeft      =   120
   ClientTop       =   420
   ClientWidth     =   4560
   LinkTopic       =   "Form1"
   ScaleHeight     =   3060
   ScaleWidth      =   4560
   StartUpPosition =   3  '窗口缺省
   Begin MSComDlg.CommonDialog CommonDialog1 
      Left            =   1800
      Top             =   1920
      _ExtentX        =   847
      _ExtentY        =   847
      _Version        =   393216
   End
   Begin VB.CommandButton Command1 
      Caption         =   "计算"
      Height          =   735
      Left            =   960
      TabIndex        =   0
      Top             =   720
      Width           =   1335
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Const M = 8               '（M=4~18,取偶数）在试井分析中通常取M=8比较合适
Private Function BessI0(X As Double) As Double   '计算I0的子函数
Dim ax As Double, ax1 As Double, y, Bessel As Double
If Abs(X) < 3.75 Then                     '当X小于3.75时的计算公式
   y = X * X / (3.75 * 3.75)
   Bessel = 1# + y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.0360768 + y * 0.0045813)))))
Else                                     '当X大于3.75时的计算公式
   ax = Abs(X)
   ax1 = Exp(ax)
   y = 3.75 / ax
   Bessel = (ax1 / Sqr(ax)) * (0.39894228 + y * (0.01328592 + y * (0.00225319 + y * (-0.00157565 + y * (0.00916281 + y * (-0.02057706 + y * (0.02635537 + y * (-0.01647633 + y * 0.00392377))))))))
End If
BessI0 = Bessel                          '返回函数值
End Function
Private Function BessK0(X As Double) As Double   '计算K0的子函数
Dim temp As Double, y, Bessel As Double
If X <= 2 Then                           '当参数小于2时计算公式
    y = X * X / 4
    Bessel = (-Log(X / 2#) * BessI0(X)) + (-0.57721566 + y * (0.4227842 + y * (0.23069756 + y * (0.0348859 + y * (0.00262698) + y * (0.0001075 + y * 0.0000074)))))
Else                                    '当参数大于2时计算公式
   y = 2# / X
   temp = Exp(-X)
   Bessel = (temp / Sqr(X)) * (1.25331414 + y * (-0.07832358 + y * (0.02189568 + y * (-0.01062446 + y * (0.00587872 + y * (-0.0025154 + y * 0.00053208))))))
End If
BessK0 = Bessel
End Function
 '计算I1的子函数
Private Function BessI1(X As Double) As Double
Dim ax As Double, ax1 As Double, y As Double
If Abs(X) < 3.75 Then
   y = X * X / (3.75 * 3.75)
   ans = y * (0.15084934 + y * (0.02658733 + y * (0.00301532 + y * 0.00032411)))
   ans = y * (0.87890594 + y * (0.51498869 + ans))
   ans = X * (0.5 + ans)
   'ans = X * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934 + y * (0.02658733 + y * (0.00301532 + y * 0.00032411))))))
Else
    ax = Abs(X)
    ax1 = Exp(ax)
    y = 3.75 / ax
    Bessel = 0.02282967 + y * (-0.02895312 + y * (0.01787654 - y * 0.00420059))
    Bessel = 0.39894228 + y * (-0.03988024 + y * (-0.00362018 + y * (0.00163801 + y * (-0.01031555 + y * BessI1))))
    Bessel = (ax1 / Sqr(ax)) * Bessel
End If
BessI1 = Bessel
End Function
'计算K1的子函数
Private Function BessK1(X As Double) As Double
  Dim temp As Double, y, Bessel As Double
  If X <= 2 Then
    y = X * X / 4#
    Bessel = (Log(X / 2#) * BessI1(X)) + (1# / X) * (1# + y * (0.15443144 + y * (-0.672778579 + y * (-0.18156897 + y * (-0.01919402 + y * (-0.00110404 + y * (-0.00004686)))))))
  Else
    y = 2# / X
    temp = Exp(-X)
    Bessel = (temp / Sqr(X)) * (1.25331414 + y * (0.23498619 + y * (-0.0365562 + y * (0.01504268 + y * (-0.00780353 + y * (0.00325614 + y * (-0.00068245)))))))
  End If
  BessK1 = Bessel
End Function
'阶乘函数（多次递归调用影响效率）
Private Function fact(n As Integer) As Double
  If n = 0 Then
     fact = 1
  Else
     fact = n * fact(n - 1)
  End If
End Function
 '三点法求导数（试井分析意义下的压力导数）
Private Sub derivative(tD() As Double, ft() As Double, D_ft() As Double, point_num As Integer)
'ReDim tD(point_num) As Double, ft(point_num) As Double, D_ft(point_num) As Double
  Dim i As Integer
  D_ft(0) = tD(0) * (-3 * ft(0) + 4 * ft(1) - ft(2)) / (4 * tD(1) - 3 * tD(0) - tD(2))
  D_ft(1) = tD(1) * (ft(2) - ft(0)) / (tD(2) - tD(0))
  For i = 2 To point_num - 1
    D_ft(i) = tD(i) * (ft(i - 2) - 4 * ft(i - 1) + 3 * ft(i)) / (tD(i - 2) - 4 * tD(i - 1) + 3 * tD(i))
  Next i

End Sub


Private Sub Command1_Click()
Const point_num = 100   '设置tD~ft的总点数
Const CDE_2S_num = 20    '设置参数的总数
Dim i As Integer
Dim tD(point_num) As Double, ft(point_num) As Double, D_ft(point_num) As Double, CDE_2S(CDE_2S_num) As Double
    'tD()无因次时间函数      ft()无因次拟压力函数     D_ft()无因次拟压力导数函数 CDE_2S（）CDE_2S()的取值函数

For i = 0 To point_num - 1
  tD(i) = 0.01 * 10 ^ (9 * i / 99)           '无因次时间的取值
Next

For i = 0 To 20
  CDE_2S(i) = 0.01 * 10 ^ i                       'CDE(2s)的取值
Next
'输出数组
Dim fileName As String
On Error GoTo errhandle
CommonDialog1.InitDir = App.Path '
CommonDialog1.Filter = "*.txt|*.txt"
CommonDialog1.ShowSave '显示"打开"对话框
fileName = CommonDialog1.fileName '保存文件名
If Len(CommonDialog1.fileName) > 0 Then
    'File = FreeFile() '获得可用文件号
    Open fileName For Output As #1 '打开文件
End If

MousePointer = 11
For i = 0 To CDE_2S_num - 1
  Print #1, "CDE_2S:", CDE_2S(i)
  Call stehfest(tD, CDE_2S(i), ft, point_num)       '调用拉氏数值逆变换"Stehfest算法"
  Call derivative(tD, ft, D_ft, point_num)                '拟压力导数
  Print #1, "无因次时间tD", "拟压力", "拟压力导数"
  For j = 0 To point_num - 1
    Print #1, tD(j), ft(j), D_ft(j)
  Next
Next
Close #1
errhandle:
MousePointer = 0

End Sub


Private Sub stehfest(tD() As Double, CDE_2S As Double, ft() As Double, point_num As Integer)
'ReDim tD(point_num) As Double, ft(point_num) As Double
Dim i As Integer, j As Integer, k As Integer, min As Integer, flag As Integer
Dim u As Double, s As Double, temp As Double
Dim v(M + 1) As Double

If (M / 2 + 1) / 2 = 0 Then
   flag = 1
Else
   flag = -1
End If

 '求反解系数v(1:n)
For j = 1 To M
  v(j) = 0
   '判断min(j,N/2)的值
  If j <= M / 2 Then
    min = j
  Else
    min = M / 2
  End If
  For k = Int((j + 1) / 2) To min
   '计算系数V(j)
     v(j) = v(j) + k ^ (M / 2) * fact(2 * k) / (fact(M / 2 - k) * fact(k) * fact(k - 1) * fact(j - k) * fact(2 * k - j))
  Next k
  v(j) = flag * v(j)
  flag = -flag
  'MsgBox ("v=" & v(j))
Next j

 '反解laplace空间压力到实空间
 For i = 0 To point_num - 1 '0
   ft(i) = 0
   u = Log(2) / tD(i)
   For j = 1 To M
     s = j * u
     temp = Sqr(s / CDE_2S)
     ft(i) = ft(i) + (v(j) / s) * BessK0(temp) / (s * BessK0(temp) + temp * BessK1(temp))
   Next
   ft(i) = ft(i) * u
   'MsgBox ("ft=" & ft(i))
 Next
 End Sub