SimulationCore2/FluidFlow/Annulus_and_Openhole_Pressure_Distribution.f90

SUBROUTINE PressureAnnAndOHDistribution
    
    !!      Record of revisions
    !!          Date            Programmer          Discription of change
    !!         ------          ------------        -----------------------
    !!        1396/07/30        Sheikh              Original code
    !! 
    
    USE FricPressDropVarsModule
    USE MudSystemVARIABLES
use SimulationVariables !@@@
    use PressureDisplayVARIABLESModule
    USE GeoElements_FluidModule
    USE Fluid_Flow_Startup_Vars
    use KickVARIABLESModule
    USE CMudPropertiesVariables
    use SimulationVariables !@
    USE CReservoirVariables
    use MudSystemModule
    USE CHOKEVARIABLES
use SimulationVariables !@
    use CChokeManifoldVariables
    use SimulationVariables
    use SimulationVariables
    USE CError
    USE , INTRINSIC :: IEEE_ARITHMETIC

    
    IMPLICIT NONE
    
    INTEGER :: i , j , k , l
    INTEGER :: ifric
    REAL :: Fraction
    
    
    data%State%FricPressDrop%KBOP = 0.0
        

            
IF (KickVARIABLES%WellHeadOpen .OR. KickVARIABLES%NoGasPocket == 0) THEN    !! (mud circulation is normal wellhead may be open or closed) OR (kick is in the well and well head is open)
        
    
    !!!!!   Determining flow rate in each section
    i = data%State%FricPressDrop%AnnulusFirstEl
    j = data%State%FricPressDrop%OpenholeFirstEl - 1
        
    !!!!!!!!!!!!!!!!!!!!!!!!! flowrates due to external sources like pump and tripping
    FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = (ClingingFactor * FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%Area + FlowEl(data%State%FricPressDrop%StringFirstEl)%Area) * KickVARIABLES%DrillStringSpeed * ConvMintoSec * Convft3toUSgal     ! flowrate in annulus due to tripping
        
    FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate + REAL(data%State%MudSystem%MudVolume_InjectedToBH) * ConvMintoSec / dt       ! flowrate in annulus due to pump
                
    IF (data%State%MudSystem%ShoeFractured) THEN     ! reduction of flowrate due to formation fracture and lost circulation
            
        IF (ShoeFlowElNo > data%State%FricPressDrop%AnnulusLastEl) THEN      ! shoe is in openhole
            FlowEl(ShoeFlowElNo : data%State%FricPressDrop%NumbEl)%FlowRate = - data%State%MudSystem%Qlost
            FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate - data%State%MudSystem%Qlost
        ELSE        ! shoe is in annulus
            FlowEl(ShoeFlowElNo : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = FlowEl(ShoeFlowElNo : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate - data%State%MudSystem%Qlost
        END IF
    END IF
        
        
    !!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!! initial guess flowrates for opening BOP or choke line
    IF (KickVARIABLES%WellHeadWasOpen == .FALSE. .AND. KickVARIABLES%NoGasPocket > 0 .AND. KickIteration == 1) THEN
        IF (KickVARIABLES%ChokeKroneckerDelta == 1) THEN      ! flow on choke line
            IF (data%State%FricPressDrop%TotalOpenChokeArea < 0.01 * data%State%Choke%ChokeAreaFullyOpen) THEN
                WRITE (*,*) 'density , TotalOpenChokeArea' , DownHole%Density, data%State%FricPressDrop%TotalOpenChokeArea
                data%State%FricPressDrop%TotalOpenChokeArea = 0.01 * data%State%Choke%ChokeAreaFullyOpen
            END IF
            data%State%FricPressDrop%Kchoke = (KickVARIABLES%ChokeDensity / ((2.0 * 89158.0) * (0.26 * 0.61 * data%State%FricPressDrop%TotalOpenChokeArea)**2)) * 4.0    ! *4.d0: seyyed gofte
            GasPocketFlowInduced%Array(:) = MIN((0.6 / KickVARIABLES%NoGasPocket * SQRT(data%State%PressureDisplay%PressureGauges(2) / data%State%FricPressDrop%Kchoke)) , (0.05 * GasPocketNewVol%Array(:) * ConvFt3toUSgal / 60 / dt))
            WRITE (*,*) ' PressureGauges(2) , Kchoke' , data%State%PressureDisplay%PressureGauges(2) , data%State%FricPressDrop%Kchoke
            WRITE (*,*) 'Initial guess after opening choke =', GasPocketFlowInduced%Array(1)
            
            WRITE (*,*) ' valve 49 ', Manifold%Valve(49)%Status
            WRITE (*,*) ' valve 47 ', Manifold%Valve(47)%Status
            WRITE (*,*) ' valve 26 ', Manifold%Valve(26)%Status
            WRITE (*,*) ' valve 30 ', Manifold%Valve(30)%Status
            WRITE (*,*) ' valve 34 ', Manifold%Valve(34)%Status
            WRITE (*,*) ' valve 63 ', Manifold%Valve(63)%Status
            WRITE (*,*) ' valve 28 ', Manifold%Valve(28)%Status
            WRITE (*,*) ' valve 33 ', Manifold%Valve(33)%Status
            WRITE (*,*) ' valve 62 ', Manifold%Valve(62)%Status
            WRITE (*,*) ' valve 36 ', Manifold%Valve(36)%Status
            WRITE (*,*) ' valve 38 ', Manifold%Valve(38)%Status
            
        ELSE        ! flow through bell nipple
            k = data%State%FricPressDrop%NoHorizontalEl + data%State%FricPressDrop%NoStringEl + data%State%FricPressDrop%NoAnnulusEl
            data%State%FricPressDrop%KBOP = FlowEl(data%State%FricPressDrop%AnnulusLastEl)%Density / ((2.0 * 89158.0) * (0.26 * 0.61 * data%State%ShearRAM%MinimumOpenArea_InBOP)**2)
            GasPocketFlowInduced%Array(:) = MIN((0.1 / KickVARIABLES%NoGasPocket * SQRT(data%State%PressureDisplay%PressureGauges(6) / data%State%FricPressDrop%KBOP)) , (0.05 * GasPocketNewVol%Array(:) * ConvFt3toUSgal / 60 / dt))
            WRITE (*,*) 'PressureGauges(6), KBOP', data%State%PressureDisplay%PressureGauges(6), data%State%FricPressDrop%KBOP
            WRITE (*,*) 'Initial guess after opening BOP =', GasPocketFlowInduced%Array(1)
        END IF
    END IF
        !!!!!!!!!!!!!!!!!!!!!!!!!
        
        !!!!!!!!!!!!!!!!!!!!!!!!! flowrates due to expansion of gas pockets or kick influx
        !i = AnnulusFirstEl
        !j = OpenholeFirstEl - 1
        IF (KickVARIABLES%NoGasPocket > 0) THEN
            DO l = 1 , KickVARIABLES%NoGasPocket          !GasPocketFlowEl
                k = KickVARIABLES%GasPocketFlowEl(l , 1)
                !WRITE (*,*) 'KickVARIABLES%GasPocketFlowEl(l , 1)', l, k, j
                IF (k == 0)     CALL ERRORSTOP('GasPocketFlowEl(l , 1) == 0', l)
               
                IF (k >= data%State%FricPressDrop%OpenholeFirstEl) THEN       ! gas pocket is in open hole only
                    FlowEl(k : data%State%FricPressDrop%NumbEl)%FlowRate = FlowEl(k : data%State%FricPressDrop%NumbEl)%FlowRate + GasPocketFlowInduced%Array(l)     ! openhole elements above pocket
                    FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate + GasPocketFlowInduced%Array(l)     ! annulus and choke line elements
                ELSE IF (k < data%State%FricPressDrop%OpenholeFirstEl) THEN       ! gas pocket is in annulus ond/or choke line only
                    FlowEl(k : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate =  FlowEl(k : data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate + GasPocketFlowInduced%Array(l)     ! annulus or choke line elements above pocket
                END IF
            END DO
        END IF

        
        !!!!!!!!!!!!!!!!!!!!!!!!!
        !!!!!   END - Determining flow rate in each section
        
        !!!!!!!!!!!!!!!!!!!!!!!!! effect of surge and swab on frictional pressure drop direction
        DO l = data%State%FricPressDrop%AnnulusFirstEl , data%State%FricPressDrop%OpenholeFirstEl - 1
            IF (FlowEl(l)%FlowRate < 0.0) THEN
                FlowEl(l)%FrictionDirection = -1
                IF (FlowEl(l)%FlowRate > -1.0 * PressFlowrateTolerance .AND. ALLOCATED(GasPocketWeight%Array))      FlowEl(l)%FlowRate = - PressFlowrateTolerance
            ELSE
                FlowEl(l)%FrictionDirection = 1
                IF (FlowEl(l)%FlowRate < PressFlowrateTolerance .AND. ALLOCATED(GasPocketWeight%Array))             FlowEl(l)%FlowRate = PressFlowrateTolerance
            END IF
        END DO
        !!!!!!!!!!!!!!!!!!!!!!!!!

        !!!!!!!!!!!!!!!!!!!!!!!!! Calculating Back Pressure, in well to pit path back pressure = 0
                                ! in well to choke manifold path back pressure is equal to pressure before choke not casing pressure
        IF (KickVARIABLES%ChokeKroneckerDelta == 1) THEN

            IF (FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate < 0.0) THEN
                WRITE (*,*) ' Negative choke flowrate'
                FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate = MAX((REAL(data%State%MudSystem%MudVolume_InjectedToBH) * ConvMintoSec / dt) , 10.0)
            END IF
            data%State%MudSystem%deltaPchoke = (data%State%FricPressDrop%Kchoke * FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate * ABS(FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%FlowRate)) * 1.d0
            IF (data%State%MudSystem%deltaPchoke < 0.d0)      data%State%MudSystem%deltaPchoke = 0.d0
            data%State%FricPressDrop%BackPressure = REAL(data%State%MudSystem%deltaPchoke)
        ELSE
            data%State%FricPressDrop%BackPressure = 0.0
        END IF
        IF (IEEE_IS_NaN(data%State%FricPressDrop%BackPressure))      CALL ErrorStop('NaN in calculating back pressure' , FlowEl(j)%FlowRate)
        !write(*,*) 'BackPressure=' , BackPressure
        
        !!!!!!!!!!!!!!!!!!!!!!!!! when flow passes through choke manifold, solution process may be unstable
        IF (KickVARIABLES%ChokeKroneckerDelta == 1) THEN       ! thus we should stabilize solution
            IF (data%State%FricPressDrop%TotalOpenChokeArea > 0.5 * data%State%Choke%ChokeAreaFullyOpen) THEN
                KickVARIABLES%KickCorrectionUnderRelaxation = 0.6     
            ELSE IF (data%State%FricPressDrop%TotalOpenChokeArea > 0.1 * data%State%Choke%ChokeAreaFullyOpen) THEN
                KickVARIABLES%KickCorrectionUnderRelaxation = 0.5
            ELSE    ! TotalOpenChokeArea < 0.1 * ChokeAreaFullyOpen
                KickVARIABLES%KickCorrectionUnderRelaxation = 0.4
            END IF
        ELSE
            KickVARIABLES%KickCorrectionUnderRelaxation = 0.6
        END IF
        !!!!!!!!!!!!!!!!!!!!!!!!!
    
        !!!!!!!!!!!!!!!!!!!!!!!!!   calculating frictional pressure drop in annulus, chooke line and open hole elements
        DO ifric = data%State%FricPressDrop%AnnulusFirstEl , data%State%FricPressDrop%NumbEl
            CALL FricPressDrop(ifric)
            !WRITE (*,*) ' element No, FlowRate , Density, FricPressLoss', ifric, FlowEl(ifric)%FlowRate, FlowEl(ifric)%Density, FlowEl(ifric)%FricPressLoss
            IF (IEEE_IS_NaN(FlowEl(ifric)%FricPressLoss)) THEN
                WRITE (*,*) 'H, S, A, Ch, O', data%State%FricPressDrop%NoHorizontalEl , data%State%FricPressDrop%NoStringEl , data%State%FricPressDrop%NoAnnulusEl , data%State%FricPressDrop%NoWellToChokeEl , data%State%FricPressDrop%NoOpenHoleEl
                WRITE (*,*) 'Ann/Op start, end, density, Q, mu, Type' , FlowEl(ifric)%StartX, FlowEl(ifric)%EndX, FlowEl(ifric)%Density, FlowEl(ifric)%FlowRate, FlowEl(ifric)%mueff, FlowEl(ifric)%MaterialType
                CALL ErrorStop('NaN in calculating pressure drop' , ifric)
            END IF
        
        END DO

        
        !!!!!!!!!!!!!!!!!!!!!!!!! Pressure distribution in annulus
        j = data%State%FricPressDrop%OpenholeFirstEl - 1
        FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%EndPress = data%State%FricPressDrop%BackPressure
        FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%StartPress = FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%EndPress + FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%FricPressLoss + FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%StaticPressDiff
        
                        
        DO l = data%State%FricPressDrop%OpenholeFirstEl - 2 , data%State%FricPressDrop%AnnulusFirstEl , -1
            !WRITE (*,*) '123'
            FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
            FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%FricPressLoss + FlowEl(l)%StaticPressDiff
        END DO


        !!!!!!!!!!!!!!!!! Pressure distribution in Open Hole
        FlowEl(data%State%FricPressDrop%NumbEl)%EndPress = FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%StartPress
        FlowEl(data%State%FricPressDrop%NumbEl)%StartPress = FlowEl(data%State%FricPressDrop%NumbEl)%EndPress + FlowEl(data%State%FricPressDrop%NumbEl)%FricPressLoss + FlowEl(data%State%FricPressDrop%NumbEl)%StaticPressDiff
        
        DO l = data%State%FricPressDrop%NumbEl - 1 , data%State%FricPressDrop%OpenholeFirstEl , -1
                !WRITE(*,*) ' ope'
                FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
                FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%FricPressLoss + FlowEl(l)%StaticPressDiff
                !WRITE (*,*) ' Length, static, frictional open' , FlowEl(i)%Length, FlowEl(i)%StaticPressDiff, FlowEl(i)%FricPressLoss
                
                !END IF
        END DO

ELSE    ! wellhead is closed and kick is in the well
        !WRITE (*,*) ' well head is closed'
        k = KickVARIABLES%GasPocketFlowEl(KickVARIABLES%NoGasPocket , 1)
        !WRITE (*,*) 'k, Pocket Press', k, GasPocketOldPress%Array(KickVARIABLES%NoGasPocket) - StandardPress
        i = data%State%FricPressDrop%AnnulusFirstEl
        j = data%State%FricPressDrop%OpenholeFirstEl - 1
        FlowEl(k)%StartPress = GasPocketOldPress%Array(KickVARIABLES%NoGasPocket) - StandardPress
        FlowEl(k)%EndPress = GasPocketOldPress%Array(KickVARIABLES%NoGasPocket) - StandardPress
        IF (k > data%State%FricPressDrop%OpenholeFirstEl - 1) THEN     ! Top pocket StartX is in Open hole
            !WRITE (*,*) 'here 1'
            DO l = k - 1 , data%State%FricPressDrop%OpenholeFirstEl , -1           ! below elements in openhole
                !WRITE (*,*) 'here 1-1'
                FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
                FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%StaticPressDiff
            END DO
            
            DO l = k + 1 , data%State%FricPressDrop%NumbEl                   ! Above elements in openhole
                !WRITE (*,*) 'here 1-2'
                FlowEl(l)%StartPress = FlowEl(l - 1)%EndPress
                FlowEl(l)%EndPress = FlowEl(l)%StartPress - FlowEl(l)%StaticPressDiff
            END DO
            
                FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%StartPress = FlowEl(data%State%FricPressDrop%NumbEl)%EndPress
                FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%EndPress = FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%StartPress - FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%StaticPressDiff
            
            DO l = data%State%FricPressDrop%AnnulusFirstEl + 1 , data%State%FricPressDrop%OpenholeFirstEl - 1
                FlowEl(l)%StartPress = FlowEl(l - 1)%EndPress
                FlowEl(l)%EndPress = FlowEl(l)%StartPress - FlowEl(l)%StaticPressDiff
            END DO
            
        ELSE        ! Top pocket StartX is in annulus or choke line
            
            DO l = k - 1 , data%State%FricPressDrop%AnnulusFirstEl , -1           ! below elements in annnulus
                FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
                FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%StaticPressDiff
            END DO
            
            DO l = k + 1 , data%State%FricPressDrop%OpenholeFirstEl - 1                   ! Above elements in annulus
                FlowEl(l)%StartPress = FlowEl(l - 1)%EndPress
                FlowEl(l)%EndPress = FlowEl(l)%StartPress - FlowEl(l)%StaticPressDiff
            END DO
            
                FlowEl(data%State%FricPressDrop%NumbEl)%EndPress = FlowEl(data%State%FricPressDrop%AnnulusFirstEl)%StartPress
                FlowEl(data%State%FricPressDrop%NumbEl)%StartPress = FlowEl(data%State%FricPressDrop%NumbEl)%EndPress + FlowEl(data%State%FricPressDrop%NumbEl)%StaticPressDiff
            
            DO l = data%State%FricPressDrop%NumbEl - 1 , data%State%FricPressDrop%OpenholeFirstEl , -1
                FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
                FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%StaticPressDiff
            END DO
            
        END IF
        
END IF

    !!!!!!!!!!!!!!!!!!!!!! checking pressure for preventing NaN in pressures
    DO l = data%State%FricPressDrop%OpenholeFirstEl - 1 , data%State%FricPressDrop%AnnulusFirstEl , -1        ! annulus or choke elements
        !WRITE (*,*) 'start, end' , FlowEl(i)%StartX, FlowEl(i)%EndX
        IF (IEEE_IS_NaN(FlowEl(l)%EndPress)) THEN
            WRITE (*,*) 'H, S, A, Ch, O', data%State%FricPressDrop%NoHorizontalEl , data%State%FricPressDrop%NoStringEl , data%State%FricPressDrop%NoAnnulusEl , data%State%FricPressDrop%NoWellToChokeEl , data%State%FricPressDrop%NoOpenHoleEl
            WRITE (*,*) 'Ann/Ch start, end, density, Q, mu' , FlowEl(l)%StartX, FlowEl(l)%EndX, FlowEl(l)%Density, FlowEl(l)%FlowRate, FlowEl(l)%mueff, FlowEl(l)%MaterialType
            CALL ERRORSTOP('NaN in EndPress', l)
        END IF
    END DO

    DO l = data%State%FricPressDrop%NumbEl , data%State%FricPressDrop%OpenholeFirstEl - 1 , -1         ! op elements
        !WRITE (*,*) 'start, end' , FlowEl(i)%StartX, FlowEl(i)%EndX
        IF (IEEE_IS_NaN(FlowEl(l)%EndPress)) THEN
            WRITE (*,*) 'H, S, A, Ch, O', data%State%FricPressDrop%NoHorizontalEl , data%State%FricPressDrop%NoStringEl , data%State%FricPressDrop%NoAnnulusEl , data%State%FricPressDrop%NoWellToChokeEl , data%State%FricPressDrop%NoOpenHoleEl
            WRITE (*,*) 'Op start, end, density, Q, mu' , FlowEl(l)%StartX, FlowEl(l)%EndX, FlowEl(l)%Density, FlowEl(l)%FlowRate, FlowEl(l)%mueff, FlowEl(l)%MaterialType
            CALL ERRORSTOP('NaN in EndPress', l)
        END IF
    END DO
    !!!!!!!!!!!!!!!!!!!!!!
    
    !!!!!!!!!!!!!!!!!!!!!!
    KickVARIABLES%BottomHolePress = FlowEl(data%State%FricPressDrop%OpenholeFirstEl)%StartPress

    !!!!!!!!!!!!!!!!!!!!!!

END SUBROUTINE