libfrodtube.bas

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

Attribute VB_Name = "LibFRodTube"
Option Explicit
Rem 在光杆与泵示功图仿真结果的基础上，迭加定向井杆管之间的摩擦力
Sub PRL_Adding_Frt(WellType As String, NcalWell As Integer, XWell() As Single, SetWell() As Single, FaiWell() As Single, _
                   nrod As Integer, drod() As Single, lrod() As Single, jcal As Integer, _
                   llpump() As Single, pr() As Single, vpr() As Single, apr() As Single, prl() As Single)
    Rem 120为杆柱上的离散点数
    Dim j As Integer
    Dim NcalRod(5) As Integer, NcalXpoint As Integer
    Dim mrodpoint(110) As Single, wrodpoint(110) As Single, alrod(110) As Single
    Dim Xpoint(110) As Single, qpoint(110) As Single
    Dim setpoint(110) As Single, faipoint(110) As Single
    Dim vsetpoint(110) As Single, vfaipoint(110) As Single
    
    ReDim PRLpoint(110, jcal) As Single
    ReDim frt(110, jcal) As Single
    ReDim frtsum(jcal) As Single
    
    Dim JprMax As Integer, JprMin As Integer, JDeltpr As Integer, JSpr_Deltpr As Integer
    Dim DeltDeltpr As Single, DeltSpr_Deltpr As Single
    Dim prMax As Single, prMin As Single, Deltpr As Single
    
    Call Set_Fai_Q_Cal(WellType, NcalWell, XWell, SetWell, FaiWell, _
                   nrod, drod, lrod, NcalRod, NcalXpoint, Xpoint, qpoint, setpoint, faipoint, vsetpoint, vfaipoint)
  
    
    Call MWpoint_cal(WellType, NcalWell, XWell, SetWell, FaiWell, nrod, drod, lrod, mrodpoint, wrodpoint, alrod)
    
    Call PRLpoint_Cal(nrod, drod, lrod, jcal, llpump, vpr, apr, NcalXpoint, _
                       mrodpoint, wrodpoint, alrod, PRLpoint)
    Call Frt_Cal(WellType, NcalWell, XWell, SetWell, FaiWell, _
                NcalXpoint, jcal, Xpoint, qpoint, setpoint, faipoint, vsetpoint, vfaipoint, PRLpoint, frt, frtsum)
    
    
    JprMin = 1
    JprMax = 1
    prMin = Abs(pr(1))
    prMax = Abs(pr(1))
    
    For j = 1 To jcal
         If Abs(pr(j)) < prMin Then
            prMin = pr(j)
            JprMin = j
        End If
        
        If Abs(pr(j)) > prMax Then
            prMax = pr(j)
            JprMax = j
        End If
    Next j
    
    Deltpr = 0.05 * prMax
    JDeltpr = 1
    JSpr_Deltpr = 1
    DeltDeltpr = Abs(pr(1) - Deltpr) '上冲程
    DeltSpr_Deltpr = Abs(pr(1) - (prMax - Deltpr)) '下冲程
    
    For j = 1 To jcal
        If vpr(j) > 0 Then
            If Abs(pr(j) - Deltpr) <= DeltDeltpr Then
                DeltDeltpr = Abs(pr(j) - Deltpr)
                JDeltpr = j
            End If
        End If
        
        If vpr(j) < 0 Then
            If Abs(pr(j) - (prMax - Deltpr)) <= DeltSpr_Deltpr Then
                DeltSpr_Deltpr = Abs(pr(j) - (prMax - Deltpr))
                JSpr_Deltpr = j
            End If
        End If
    Next j
    
    For j = 0 To jcal
        If vpr(j) > 0 Then
            If Abs(pr(j)) >= Deltpr Then
                prl(j) = prl(j) + frtsum(j)
            Else
                prl(j) = prl(j) - frtsum(JprMin) + (frtsum(JprMin) + frtsum(JDeltpr)) * Abs(pr(j)) / Deltpr
            End If
        Else
            If Abs(pr(j)) >= prMax - Deltpr Then
                prl(j) = prl(j) - frtsum(JSpr_Deltpr) + (frtsum(JSpr_Deltpr) + frtsum(JprMax)) * (pr(j) - (prMax - Deltpr)) / Deltpr
            Else
                prl(j) = prl(j) - frtsum(j)
            End If
        End If
    Next j
End Sub

Rem 定向井杆柱各离散点下部杆柱的质量、重量在轴线上的投影
Sub MWpoint_cal(WellType As String, NcalWell As Integer, XWell() As Single, SetWell() As Single, FaiWell() As Single, _
                nrod As Integer, drod() As Single, lrod() As Single, mrodpoint() As Single, wrodpoint() As Single, _
                alrod() As Single)
          
    Dim qrod(5) As Single, arod(5) As Single
    Dim Wrod As Single, wrod2 As Single, wrod3 As Single
    Dim i As Integer, j As Integer
    Dim NcalXpoint As Integer '整个抽油杆柱上所取的计算点数。由set_fai_q_cal（）返回，101个点
    Dim Xpoint(110) As Single, qpoint(110) As Single, NcalRod(5) As Integer
    Dim setpoint(110) As Single, faipoint(110) As Single
    Dim vsetpoint(110) As Single, vfaipoint(110) As Single
    Dim pi As Single
    
    pi = 3.14159265
    
    Call Set_Fai_Q_Cal(WellType, NcalWell, XWell, SetWell, FaiWell, _
            nrod, drod, lrod, NcalRod, NcalXpoint, Xpoint, qpoint, setpoint, faipoint, vsetpoint, vfaipoint)
    
    For i = 1 To nrod
        qrod(i) = pi / 4 * drod(i) ^ 2 * 7850 * 9.8 * 1.05
        arod(i) = pi / 4 * drod(i) ^ 2
    Next i
    
    If nrod = 1 Then
        For i = 0 To NcalXpoint
            mrodpoint(i) = (lrod(1) - Xpoint(i)) * qrod(1) / 9.81
            alrod(i) = (lrod(1) - Xpoint(i)) * arod(1)
            wrodpoint(i) = 0
            For j = i + 1 To NcalXpoint
                wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                      (Xpoint(j) - Xpoint(j - 1))
            Next j
        Next i
    End If
        
    If nrod = 2 Then
        wrod2 = 0
        For j = NcalRod(1) + 1 To NcalXpoint
           setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
           qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
           wrod2 = wrod2 + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                      (Xpoint(j) - Xpoint(j - 1))
        Next
        
        For i = 0 To NcalRod(1)
            mrodpoint(i) = (lrod(1) - Xpoint(i)) * qrod(1) / 9.81 + lrod(2) * qrod(2) / 9.81
            alrod(i) = (lrod(1) - Xpoint(i)) * arod(1) + lrod(2) * arod(2)
            wrodpoint(i) = wrod2
            For j = i + 1 To NcalRod(1)
                wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + _
                        qpoint(j) * Cos(setpoint(j - 1))) / 2 * (Xpoint(j) - Xpoint(j - 1))
            Next j
        Next i
         
        For i = NcalRod(1) + 1 To NcalXpoint
           mrodpoint(i) = (lrod(1) + lrod(2) - Xpoint(i)) * qrod(2) / 9.81
           alrod(i) = (lrod(1) + lrod(2) - Xpoint(i)) * arod(2)
           wrodpoint(i) = 0
           setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
           qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
           For j = i + 1 To NcalXpoint
               wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                      (Xpoint(j) - Xpoint(j - 1))
           Next j
        Next
    End If
        
    If nrod = 3 Then
        wrod2 = 0
        For j = NcalRod(1) + 1 To NcalRod(1) + NcalRod(2)
            setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
            qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
            wrod2 = wrod2 + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                      (Xpoint(j) - Xpoint(j - 1))
        Next
    
        wrod3 = 0
        For j = NcalRod(1) + NcalRod(2) + 1 To NcalXpoint
            setpoint(NcalRod(1) + NcalRod(2)) = setpoint(NcalRod(1) + NcalRod(2) + 1)
            qpoint(NcalRod(1) + NcalRod(2)) = qpoint(NcalRod(1) + NcalRod(2) + 1)
            wrod3 = wrod3 + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                      (Xpoint(j) - Xpoint(j - 1))
        Next
        
        For i = 0 To NcalRod(1)
            mrodpoint(i) = (lrod(1) - Xpoint(i)) * qrod(1) / 9.81 + lrod(2) * qrod(2) / 9.81 + lrod(3) * qrod(3) / 9.81
            alrod(i) = (lrod(1) - Xpoint(i)) * arod(1) + lrod(2) * arod(2) + lrod(3) * arod(3)
            wrodpoint(i) = wrod2 + wrod3
            For j = i + 1 To NcalRod(1)
                wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + _
                        qpoint(j) * Cos(setpoint(j - 1))) / 2 * (Xpoint(j) - Xpoint(j - 1))
            Next j
        Next i
        
        For i = NcalRod(1) + 1 To NcalRod(1) + NcalRod(2)
            mrodpoint(i) = (lrod(1) + lrod(2) - Xpoint(i)) * qrod(2) / 9.81 + lrod(3) * qrod(3) / 9.81
            alrod(i) = (lrod(1) + lrod(2) - Xpoint(i)) * arod(2) + lrod(3) * arod(3)
            wrodpoint(i) = wrod3
            setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
            qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
            For j = i + 1 To NcalRod(1) + NcalRod(2)
                wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + _
                        qpoint(j) * Cos(setpoint(j - 1))) / 2 * (Xpoint(j) - Xpoint(j - 1))
            Next j
        Next i
         
        For i = NcalRod(1) + NcalRod(2) + 1 To NcalXpoint
           mrodpoint(i) = (lrod(1) + lrod(2) + lrod(3) - Xpoint(i)) * qrod(3) / 9.81
           alrod(i) = (lrod(1) + lrod(2) + lrod(3) - Xpoint(i)) * arod(3)
           wrodpoint(i) = 0
           setpoint(NcalRod(1) + NcalRod(2)) = setpoint(NcalRod(1) + NcalRod(2) + 1)
           qpoint(NcalRod(1) + NcalRod(2)) = qpoint(NcalRod(1) + NcalRod(2) + 1)
           For j = i + 1 To NcalXpoint
               wrodpoint(i) = wrodpoint(i) + (qpoint(j - 1) * Cos(setpoint(j - 1)) + _
                  qpoint(j) * Cos(setpoint(j - 1))) / 2 * (Xpoint(j) - Xpoint(j - 1))
           Next j
        Next
    End If
End Sub

Rem 已知泵负荷仿真结果，计算抽油杆柱离散点的轴向负荷（静载+惯性负荷+杆液摩擦负荷，用于计算杆管摩擦力）
Sub PRLpoint_Cal(nrod As Integer, drod() As Single, lrod() As Single, _
                 jcal As Integer, llpump() As Single, vpr() As Single, apr() As Single, _
                 NcalXpoint As Integer, mrodpoint() As Single, wrodpoint() As Single, alrod() As Single, _
                 PRLpoint() As Single)
  
     Dim wr As Single, arod As Single, srod As Single, dprlj As Single
     Dim i As Single, j As Integer
     
     For i = 0 To NcalXpoint
        For j = 0 To jcal
            PRLpoint(i, j) = llpump(j) + mrodpoint(i) * apr(j) + wrodpoint(i)
        Next j
     Next i
    
     For i = 0 To NcalXpoint
        For j = 0 To jcal
            srod = 7850
            dprlj = 0.5 * srod * alrod(i) * vpr(j)
            PRLpoint(i, j) = PRLpoint(i, j) + dprlj
        Next j
     Next i
End Sub
Rem 计算杆柱上任意一点单位长度摩擦力和杆柱附加的总摩擦力
Sub Frt_Cal(WellType As String, NcalWell As Integer, XWell() As Single, SetWell() As Single, FaiWell() As Single, _
            NcalXpoint As Integer, jcal As Integer, Xpoint() As Single, qpoint() As Single, setpoint() As Single, _
            faipoint() As Single, vsetpoint() As Single, vfaipoint() As Single, _
            PRLpoint() As Single, frt() As Single, frtsum() As Single)
    
     Dim i As Single, j As Integer
     Dim Nrt(110, 15000) As Single, mu_ft As Single
     
     mu_ft = 0.08
     For i = 0 To NcalXpoint
        For j = 0 To jcal
            Nrt(i, j) = (PRLpoint(i, j) * vsetpoint(i) - qpoint(i) * Sin(setpoint(i))) ^ 2 + _
                         (PRLpoint(i, j) * Sin(setpoint(i)) * vfaipoint(i)) ^ 2
            Nrt(i, j) = Sqr(Nrt(i, j))
            frt(i, j) = Nrt(i, j) * mu_ft
        Next j
     Next i
     For j = 0 To jcal
         frtsum(j) = 0
         For i = 1 To NcalXpoint
             frtsum(j) = frtsum(j) + (frt(i - 1, j) + frt(i, j)) / 2 * (Xpoint(i) - Xpoint(i - 1))
         Next i
     Next j
End Sub

Rem 直井与定向井杆柱空气中重量在井眼轴线上的投影
Sub Wrod_Cal(WellType As String, NcalWell As Integer, XWell() As Single, SetWell() As Single, FaiWell() As Single, _
             nrod As Integer, drod() As Single, lrod() As Single, Wrod As Single)
         
     Dim i As Integer, j As Integer, NcalXpoint As Integer
     Dim qrod(5) As Single
     Dim NcalRod(5) As Integer
     Dim Xpoint(110) As Single, qpoint(110) As Single
     Dim setpoint(110) As Single, faipoint(110) As Single
     Dim vsetpoint(110) As Single, vfaipoint(110) As Single
     Dim pi As Single
    
     pi = 3.14159265
   
     If WellType = "定向井" Then
         Call Set_Fai_Q_Cal(WellType, NcalWell, XWell, SetWell, FaiWell, _
                    nrod, drod, lrod, NcalRod, NcalXpoint, Xpoint, qpoint, setpoint, faipoint, vsetpoint, vfaipoint)
         
         Wrod = 0
         If nrod = 1 Then
             For j = 1 To NcalXpoint
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next j
         End If
         
         If nrod = 2 Then
             For j = 1 To NcalRod(1)
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next j
              
             For j = NcalRod(1) + 1 To NcalXpoint
                setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
                qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next
         End If
         
         If nrod = 3 Then
             For j = 1 To NcalRod(1)
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next j
              
             For j = NcalRod(1) + 1 To NcalRod(1) + NcalRod(2)
                setpoint(NcalRod(1)) = setpoint(NcalRod(1) + 1)
                qpoint(NcalRod(1)) = qpoint(NcalRod(1) + 1)
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next
             
             For j = NcalRod(1) + NcalRod(2) + 1 To NcalXpoint
                setpoint(NcalRod(1) + NcalRod(2)) = setpoint(NcalRod(1) + NcalRod(2) + 1)
                qpoint(NcalRod(1) + NcalRod(2)) = qpoint(NcalRod(1) + NcalRod(2) + 1)
                Wrod = Wrod + (qpoint(j - 1) * Cos(setpoint(j - 1)) + qpoint(j) * Cos(setpoint(j - 1))) / 2 * _
                           (Xpoint(j) - Xpoint(j - 1))
             Next
         End If
     Else
         Rem 直井
         For i = 1 To nrod
              qrod(i) = pi / 4 * drod(i) ^ 2 * 7850 * 9.8 * 1.05
         Next i
         Wrod = 0
         For i = 1 To nrod
            Wrod = Wrod + qrod(i) * lrod(i)
         Next i
     End If
End Sub