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