libwellcurve.bas

本系统是给大庆油田做的一个示例程序

        faixpoint_cal = F(1)
        vfaixpoint_cal = F1(1)
    Else
         setxpoint_cal = SetWell(1)
         faixpoint_cal = FaiWell(1)
         vsetxpoint_cal = (SetWell(2) - SetWell(1)) / (XWell(2) - XWell(1))
         vfaixpoint_cal = (FaiWell(2) - FaiWell(1)) / (XWell(2) - XWell(1))
    End If
      
    setxpoint_cal = setxpoint_cal * pi / 180
    faixpoint_cal = faixpoint_cal * pi / 180
    vsetxpoint_cal = vsetxpoint_cal * pi / 180
    vfaixpoint_cal = vfaixpoint_cal * pi / 180
    
    'vfaixpoint_cal = 0.003
    
    K_Curve = Sqr(vsetxpoint_cal ^ 2 + vfaixpoint_cal ^ 2 * (Sin(setxpoint_cal)) ^ 2)
    
    If K_Curve >= 0.02 Then K_Curve = 0.02
    
    If K_Curve <= 0.00001 Then
        R_Curve = 10000
    Else
        R_Curve = 1 / K_Curve
    End If
End Sub

Rem 井眼轴线上任意一点的全角变化率
Sub VAnglePerm_Cal(WellType As String, NcalWell As Integer, XWell() As Single, _
                SetWell() As Single, FaiWell() As Single, _
                xpoint_cal As Single, VAnglePerm As Single)
  
    Dim i As Integer, j As Integer, deltx As Single
    Dim n As Integer, M As Integer
    Dim setxpoint_cal As Single, faixpoint_cal As Single
    Dim vsetxpoint_cal As Single, vfaixpoint_cal As Single
    Dim x(300) As Single, Y(300) As Single, F(300) As Single
    Dim F1(300) As Single, F2(300) As Single, t(300) As Single
    Dim pi As Single
    
    pi = 3.14159265
    
    If xpoint_cal > XWell(1) And xpoint_cal < XWell(NcalWell) Then
        n = NcalWell
        M = 1
        
        For j = 1 To n
            x(j) = XWell(j)
            Y(j) = SetWell(j)
        Next j
        
        t(1) = xpoint_cal
        Call ChZh(n, x, Y, M, t, F, F1, F2)
        
        setxpoint_cal = F(1)
        vsetxpoint_cal = F1(1)
        
        For j = 1 To n
            x(j) = XWell(j)
            Y(j) = FaiWell(j)
        Next j
        
        t(1) = xpoint_cal
        
        Call ChZh(n, x, Y, M, t, F, F1, F2)
        
        faixpoint_cal = F(1)
        vfaixpoint_cal = F1(1)
    Else
         setxpoint_cal = SetWell(1)
         faixpoint_cal = FaiWell(1)
         vsetxpoint_cal = (SetWell(2) - SetWell(1)) / (XWell(2) - XWell(1))
         vfaixpoint_cal = (FaiWell(2) - FaiWell(1)) / (XWell(2) - XWell(1))
    End If
      
    setxpoint_cal = setxpoint_cal * pi / 180
    faixpoint_cal = faixpoint_cal * pi / 180
    vsetxpoint_cal = vsetxpoint_cal * pi / 180
    vfaixpoint_cal = vfaixpoint_cal * pi / 180
    'vfaixpoint_cal = 0.003
    
    VAnglePerm = Sqr(vsetxpoint_cal ^ 2 + (Sin(setxpoint_cal) * vfaixpoint_cal) ^ 2)
    If VAnglePerm >= 0.02 Then VAnglePerm = 0.02
End Sub
Rem 井眼轴线上任意一点曲率和曲率半径
Sub RCurve_Cal(WellType As String, NcalWell As Integer, XWell() As Single, SetWell() As Single, FaiWell() As Single, _
               xpoint_cal As Single, K_Curve As Single, R_Curve As Single)
  
    Dim i As Integer, j As Integer, deltx As Single
    Dim n As Integer, M As Integer
    
    Dim x(300) As Single, Y(300) As Single, F(300) As Single
    Dim F1(300) As Single, F2(300) As Single, t(300) As Single
    Dim setxpoint_cal As Single, faixpoint_cal As Single, vsetxpoint_cal As Single
    Dim vfaixpoint_cal As Single
    Dim pi As Single
    
    pi = 3.14159265
    
    If xpoint_cal > XWell(1) And xpoint_cal < XWell(NcalWell) Then
        n = NcalWell
        M = 1
        
        For j = 1 To n
             x(j) = XWell(j)
             Y(j) = SetWell(j)
        Next j
        
        t(1) = xpoint_cal
        Call ChZh(n, x, Y, M, t, F, F1, F2)
        
        setxpoint_cal = F(1)
        vsetxpoint_cal = F1(1)
        
        For j = 1 To n
            x(j) = XWell(j)
            Y(j) = FaiWell(j)
        Next j
         
        t(1) = xpoint_cal
        
        Call ChZh(n, x, Y, M, t, F, F1, F2)
        
        faixpoint_cal = F(1)
        vfaixpoint_cal = F1(1)
    Else
        setxpoint_cal = SetWell(1)
        faixpoint_cal = FaiWell(1)
        vsetxpoint_cal = (SetWell(2) - SetWell(1)) / (XWell(2) - XWell(1))
        vfaixpoint_cal = (FaiWell(2) - FaiWell(1)) / (XWell(2) - XWell(1))
    End If
      
      setxpoint_cal = setxpoint_cal * pi / 180
      faixpoint_cal = faixpoint_cal * pi / 180
      vsetxpoint_cal = vsetxpoint_cal * pi / 180
      vfaixpoint_cal = vfaixpoint_cal * pi / 180
       
      K_Curve = Sqr(vsetxpoint_cal ^ 2 + vfaixpoint_cal ^ 2 * (Sin(setxpoint_cal)) ^ 2)
      
      If K_Curve <= 0.00001 Then
          R_Curve = 10000
      Else
          R_Curve = 1 / K_Curve
      End If
End Sub


Rem 三次样条插值程序
Sub Curve_3splin(x() As Single, Y() As Single, n As Integer, dy1 As Single, dyn As Single, _
                 XX() As Single, M As Integer, Dy() As Single, DDy() As Single, s() As Single, _
                 Ds() As Single, DDs() As Single, t As Single)
              
    Dim h(300) As Single
    Dim i As Integer, j As Integer
 
    Dim h0 As Single, H1 As Single, beta As Single, ALPHA As Single
    Dim pi As Single
    
    pi = 3.14159265
        
    Dy(1) = 0#
    h(1) = dy1
    h0 = x(2) - x(1)
    
    For j = 2 To n - 1
        H1 = x(j + 1) - x(j)
        ALPHA = h0 / (h0 + H1)
        beta = (1# - ALPHA) * (Y(j) - Y(j - 1)) / h0
        beta = 3# * (beta + ALPHA * (Y(j + 1) - Y(j)) / H1)
        Dy(j) = -ALPHA / (2# + (1# - ALPHA) * Dy(j - 1))
        h(j) = (beta - (1# - ALPHA) * h(j - 1))
        h(j) = h(j) / (2# + (1# - ALPHA) * Dy(j - 1))
        h0 = H1
    Next j
    
    Dy(n) = dyn
    For j = n - 1 To 1 Step -1
        Dy(j) = Dy(j) * Dy(j + 1) + h(j)
    Next j
    
    For j = 1 To n - 1
        h(j) = x(j + 1) - x(j)
    Next j
    
    For j = 1 To n - 1
        H1 = h(j) * h(j)
        DDy(j) = 6# * (Y(j + 1) - Y(j)) / H1 - 2# * (2# * Dy(j) + Dy(j + 1)) / h(j)
    Next j
    
    H1 = h(n - 1) * h(n - 1)
    DDy(n) = 6# * (Y(n - 1) - Y(n)) / H1 + 2# * (2# * Dy(n) + Dy(n - 1)) / h(n - 1)
    
    t = 0#
    For i = 1 To n - 1
        H1 = 0.5 * h(i) * (Y(i) + Y(i + 1))
        H1 = H1 - h(i) * h(i) * h(i) * (DDy(i) + DDy(i + 1)) / 24#
        t = t + H1
    Next i
    For j = 1 To M
        If XX(j) >= x(n) Then
           i = n - 1
        Else
           i = 1
           Do
               If XX(j) > x(i + 1) Then
                   i = i + 1
               Else
                   Exit Do
               End If
           Loop
        End If
        
        H1 = (x(i + 1) - XX(j)) / h(i)
        s(j) = (3# * H1 * H1 - 2# * H1 * H1 * H1) * Y(i)
        s(j) = s(j) + h(i) * (H1 * H1 - H1 * H1 * H1) * Dy(i)
        Ds(j) = 6# * (H1 * H1 - H1) * Y(i) / h(i)
        Ds(j) = Ds(j) + (3# * H1 * H1 - 2# * H1) * Dy(i)
        DDs(j) = (6# - 12# * H1) * Y(i) / (h(i) * h(i))
        DDs(j) = DDs(j) + (2# - 6# * H1) * Dy(i) / h(i)
        H1 = (XX(j) - x(i)) / h(i)
        s(j) = s(j) + (3# * H1 * H1 - 2# * H1 * H1 * H1) * Y(i + 1)
        s(j) = s(j) - h(i) * (H1 * H1 - H1 * H1 * H1) * Dy(i + 1)
        Ds(j) = Ds(j) - 6# * (H1 * H1 - H1) * Y(i + 1) / h(i)
        Ds(j) = Ds(j) + (3# * H1 * H1 - 2# * H1) * Dy(i + 1)
        DDs(j) = DDs(j) + (6# - 12# * H1) * Y(i + 1) / (h(i) * h(i))
        DDs(j) = DDs(j) - (2# - 6# * H1) * Dy(i + 1) / h(i)
    Next j
End Sub