SUBROUTINE PressureAnnAndOHDistribution
!! Record of revisions
!! Date Programmer Discription of change
!! ------ ------------ -----------------------
!! 1396/07/30 Sheikh Original code
!!
USE FricPressDropVars
USE MudSystemVARIABLES
USE PressureDisplayVARIABLES
USE GeoElements_FluidModule
USE Fluid_Flow_Startup_Vars
USE KickVariables
USE CMudPropertiesVariables
USE TD_WellGeometry
USE CReservoirVariables
USE MudSystem
USE CHOKEVARIABLES
USE CChokeManifoldVariables
USE VARIABLES
USE CError
USE , INTRINSIC :: IEEE_ARITHMETIC
IMPLICIT NONE
INTEGER :: i , j , k , l
INTEGER :: ifric
REAL :: Fraction
KBOP = 0.0
IF (WellHeadOpen .OR. 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 = AnnulusFirstEl
j = OpenholeFirstEl - 1
!!!!!!!!!!!!!!!!!!!!!!!!! flowrates due to external sources like pump and tripping
!WRITE (*,*) 'StringFlowRate', StringFlowRate
FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate = (ClingingFactor * FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%Area + FlowEl(StringFirstEl)%Area) * DrillStringSpeed * ConvMintoSec * Convft3toUSgal ! flowrate in annulus due to tripping
FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate = FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate + REAL(MudVolume_InjectedToBH) * ConvMintoSec / dt ! flowrate in annulus due to pump
!WRITE (*,*) 'Drillstring speed (ft/s)' , FlowEl(j)%FlowRate
!IF (NoWellToChokeEl > 0) THEN ! flowrate in choke line
! FlowEl(NoHorizontalEl + NoStringEl + NoAnnulusEl + 1 : NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWellToChokeEl)%FlowRate = AnnulusFlowRate + (DeltaVolumePipe * ConvMinToSec / dt)
!END IF
IF (ShoeFractured) THEN ! reduction of flowrate due to formation fracture and lost circulation
!WRITE (*,*) ' SHoe fractured', PressureGauges(5), FlowEl(ShoeFlowElNo)%FlowRate
IF (ShoeFlowElNo > AnnulusLastEl) THEN ! shoe is in openhole
FlowEl(ShoeFlowElNo : NumbEl)%FlowRate = - QLost
FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate = FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate - QLost
ELSE ! shoe is in annulus
FlowEl(ShoeFlowElNo : OpenholeFirstEl - 1)%FlowRate = FlowEl(ShoeFlowElNo : OpenholeFirstEl - 1)%FlowRate - QLost
END IF
END IF
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!! initial guess flowrates for opening BOP or choke line
IF (WellHeadWasOpen == .FALSE. .AND. NoGasPocket > 0 .AND. KickIteration == 1) THEN
IF (ChokeKroneckerDelta == 1) THEN ! flow on choke line
IF (TotalOpenChokeArea < 0.01 * ChokeAreaFullyOpen) THEN
WRITE (*,*) 'density , TotalOpenChokeArea' , DownHole%Density, TotalOpenChokeArea
TotalOpenChokeArea = 0.01 * ChokeAreaFullyOpen
END IF
Kchoke = (ChokeDensity / ((2.0 * 89158.0) * (0.26 * 0.61 * TotalOpenChokeArea)**2)) * 4.0 ! *4.d0: seyyed gofte
GasPocketFlowInduced%Array(:) = MIN((0.6 / NoGasPocket * SQRT(PressureGauges(2) / Kchoke)) , (0.05 * GasPocketNewVol%Array(:) * ConvFt3toUSgal / 60 / dt))
WRITE (*,*) ' PressureGauges(2) , Kchoke' , PressureGauges(2) , 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 = NoHorizontalEl + NoStringEl + NoAnnulusEl
KBOP = FlowEl(AnnulusLastEl)%Density / ((2.0 * 89158.0) * (0.26 * 0.61 * MinimumOpenArea_InBOP)**2)
GasPocketFlowInduced%Array(:) = MIN((0.1 / NoGasPocket * SQRT(PressureGauges(6) / KBOP)) , (0.05 * GasPocketNewVol%Array(:) * ConvFt3toUSgal / 60 / dt))
WRITE (*,*) 'PressureGauges(6), KBOP', PressureGauges(6), 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 (NoGasPocket > 0) THEN
DO l = 1 , NoGasPocket !GasPocketFlowEl
k = GasPocketFlowEl(l , 1)
!WRITE (*,*) 'GasPocketFlowEl(l , 1)', l, k, j
IF (k == 0) CALL ERRORSTOP('GasPocketFlowEl(l , 1) == 0', l)
IF (k >= OpenholeFirstEl) THEN ! gas pocket is in open hole only
FlowEl(k : NumbEl)%FlowRate = FlowEl(k : NumbEl)%FlowRate + GasPocketFlowInduced%Array(l) ! openhole elements above pocket
FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate = FlowEl(AnnulusFirstEl : OpenholeFirstEl - 1)%FlowRate + GasPocketFlowInduced%Array(l) ! annulus and choke line elements
ELSE IF (k < OpenholeFirstEl) THEN ! gas pocket is in annulus ond/or choke line only
FlowEl(k : OpenholeFirstEl - 1)%FlowRate = FlowEl(k : OpenholeFirstEl - 1)%FlowRate + GasPocketFlowInduced%Array(l) ! annulus or choke line elements above pocket
END IF
END DO
END IF
!IF (ChokeKroneckerDelta == 1 .AND. ABS(FlowEl(i + NoAnnulusEl)%FlowRate / 600.0 - Ann_Saved_MudDischarged_Volume_Final) > 0.05) THEN
! WRITE (*,*) 'Difference between flowrates', FlowEl(i + NoAnnulusEl + 1)%FlowRate / 600.0, Ann_Saved_MudDischarged_Volume_Final
!END IF
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! END - Determining flow rate in each section
!!!!!!!!!!!!!!!!!!!!!!!!! effect of surge and swab on frictional pressure drop direction
DO l = AnnulusFirstEl , 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 (ChokeKroneckerDelta == 1) THEN
IF (FlowEl(OpenholeFirstEl - 1)%FlowRate < 0.0) THEN
WRITE (*,*) ' Negative choke flowrate'
FlowEl(OpenholeFirstEl - 1)%FlowRate = MAX((REAL(MudVolume_InjectedToBH) * ConvMintoSec / dt) , 10.0)
END IF
!Kchoke = ChokeDensity / ((2. * 89158.0) * (0.26 * 0.61 * TotalOpenChokeArea)**2)
deltaPchoke = (Kchoke * FlowEl(OpenholeFirstEl - 1)%FlowRate * ABS(FlowEl(OpenholeFirstEl - 1)%FlowRate)) * 1.d0
!WRITE (*,*) '**deltaPchoke , Kchoke, choke flowrate' , deltaPchoke , Kchoke, FlowEl(i)%FlowRate
!WRITE (*,*) '**TotalOpenChokeArea , Total Open Choke Area Percent' , TotalOpenChokeArea , TotalOpenChokeArea / 4.0 * ChokeAreaFullyOpen
IF (deltaPchoke < 0.d0) deltaPchoke = 0.d0
BackPressure = REAL(deltaPchoke)
!WRITE (*,*) ' Choke inlet FlowRate, Density, pressure' , FlowEl(j)%FlowRate, FlowEl(j)%Density, FlowEl(j)%StartPress
!WRITE (*,*) ' Choke outlet Density' , FlowEl(i)%Density
!WRITE (*,*) ' deltaPchoke , choke flowrate' , deltaPchoke , FlowEl(i)%FlowRate
!WRITE (*,*) 'Total Open Choke Area Percent' , TotalOpenChokeArea / 4.0 * ChokeAreaFullyOpen
ELSE
BackPressure = 0.0
END IF
IF (IEEE_IS_NaN(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 (ChokeKroneckerDelta == 1) THEN ! thus we should stabilize solution
IF (TotalOpenChokeArea > 0.5 * ChokeAreaFullyOpen) THEN
KickCorrectionUnderRelaxation = 0.6
ELSE IF (TotalOpenChokeArea > 0.1 * ChokeAreaFullyOpen) THEN
KickCorrectionUnderRelaxation = 0.5
ELSE ! TotalOpenChokeArea < 0.1 * ChokeAreaFullyOpen
KickCorrectionUnderRelaxation = 0.4
END IF
ELSE
KickCorrectionUnderRelaxation = 0.6
END IF
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!! calculating frictional pressure drop in annulus, chooke line and open hole elements
DO ifric = AnnulusFirstEl , 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', NoHorizontalEl , NoStringEl , NoAnnulusEl , NoWellToChokeEl , 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
!IF (ChokeKroneckerDelta == 1) THEN
!WRITE (*,*) ' velocity and flowrate', FlowEl(i)%vel, FlowEl(i)%flowrate
!WRITE (*,*) ' Theta600, Theta300', FlowEl(i)%Theta600 , FlowEl(i)%Theta300
!WRITE (*,*) ' kIndex , nIndex', FlowEl(i)%kIndex, FlowEl(i)%nIndex
!WRITE (*,*) ' last el. mueff, gen. Rey.', i, FlowEl(i)%mueff, FlowEl(i)%GenRe
!END IF
!!!!!!!!!!!!!!!!!!!!!!!!! Pressure distribution in annulus
j = OpenholeFirstEl - 1
FlowEl(OpenholeFirstEl - 1)%EndPress = BackPressure
FlowEl(OpenholeFirstEl - 1)%StartPress = FlowEl(OpenholeFirstEl - 1)%EndPress + FlowEl(OpenholeFirstEl - 1)%FricPressLoss + FlowEl(OpenholeFirstEl - 1)%StaticPressDiff
!write(*,*) 'FlowEl(j)%StartPress=' ,j, FlowEl(j)%StartPress
!write(*,*) 'FlowEl(j)%Length=' ,j, FlowEl(j)%Length, FlowEl(j)%EndX
!write(*,*) 'FlowEl(i)%dPdLFric=' ,i, FlowEl(i)%dPdLFric
DO l = OpenholeFirstEl - 2 , AnnulusFirstEl , -1
!WRITE (*,*) '123'
FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%FricPressLoss + FlowEl(l)%StaticPressDiff
!WRITE(*,*) "ANNULUS: bottom , top Pressure", l , FlowEl(l)%StartPress , FlowEl(l)%EndPress , FlowEl(l)%fricPressLoss
!WRITE(*,*) "ANNULUS: Start , End X", FlowEl(l)%StartX , FlowEl(l)%EndX
!write(*,*) 'FlowEl(i)%StartPress=' ,i, FlowEl(i)%StartPress
!WRITE (*,*) ' FlowEl(i)%GenRe, FlowEl(i)%ReCritLam ' , FlowEl(i)%GenRe , FlowEl(i)%ReCritLam
END DO
!!!!!!!!!!!!!!!!! Pressure distribution in Open Hole
FlowEl(NumbEl)%EndPress = FlowEl(AnnulusFirstEl)%StartPress
FlowEl(NumbEl)%StartPress = FlowEl(NumbEl)%EndPress + FlowEl(NumbEl)%FricPressLoss + FlowEl(NumbEl)%StaticPressDiff
!WRITE (*,*) 'op top and op down' , FlowEl(NumbEl)%EndPress, FlowEl(j + 1)%StartPress
!write(*,*) 'FlowEl(NumbEl)%dPdLFric=' , FlowEl(NumbEl)%dPdLFric
!write(*,*) 'FlowEl(NumbEl)%dPdLGrav=' , FlowEl(NumbEl)%dPdLGrav
DO l = NumbEl - 1 , OpenholeFirstEl , -1
!WRITE(*,*) ' ope'
FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
!IF (FlowEl(i)%FlowRate < 0.0d0) THEN
! FlowEl(i)%StartPress = FlowEl(i)%EndPress - FlowEl(i)%FricPressLoss + FlowEl(i)%StaticPressDiff
!ELSE
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 = GasPocketFlowEl(NoGasPocket , 1)
!WRITE (*,*) 'k, Pocket Press', k, GasPocketOldPress%Array(NoGasPocket) - StandardPress
i = AnnulusFirstEl
j = OpenholeFirstEl - 1
FlowEl(k)%StartPress = GasPocketOldPress%Array(NoGasPocket) - StandardPress
FlowEl(k)%EndPress = GasPocketOldPress%Array(NoGasPocket) - StandardPress
IF (k > OpenholeFirstEl - 1) THEN ! Top pocket StartX is in Open hole
!WRITE (*,*) 'here 1'
DO l = k - 1 , 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 , 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(AnnulusFirstEl)%StartPress = FlowEl(NumbEl)%EndPress
FlowEl(AnnulusFirstEl)%EndPress = FlowEl(AnnulusFirstEl)%StartPress - FlowEl(AnnulusFirstEl)%StaticPressDiff
DO l = AnnulusFirstEl + 1 , 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 , 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 , 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(NumbEl)%EndPress = FlowEl(AnnulusFirstEl)%StartPress
FlowEl(NumbEl)%StartPress = FlowEl(NumbEl)%EndPress + FlowEl(NumbEl)%StaticPressDiff
DO l = NumbEl - 1 , OpenholeFirstEl , -1
FlowEl(l)%EndPress = FlowEl(l + 1)%StartPress
FlowEl(l)%StartPress = FlowEl(l)%EndPress + FlowEl(l)%StaticPressDiff
END DO
END IF
!
! !WRITE (*,*) ' first annulus bottom pressure ' , FlowEl(NoHorizontalEl + NoStringEl + 1)%StartPress
! !WRITE (*,*) ' last OpenHole bottom pressure' , FlowEl(NumbEl)%StartPress
! !WRITE (*,*) ' Gas Pocket pressure' , GasPocket%NewPress
END IF
!!!!!!!!!!!!!!!!!!!!!! checking pressure for preventing NaN in pressures
DO l = OpenholeFirstEl - 1 , 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', NoHorizontalEl , NoStringEl , NoAnnulusEl , NoWellToChokeEl , 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 = NumbEl , 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', NoHorizontalEl , NoStringEl , NoAnnulusEl , NoWellToChokeEl , 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
!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!
BottomHolePress = FlowEl(OpenholeFirstEl)%StartPress
!DO i = 1 , NoGasPocket
! WRITE (*,*) ' Pocket, Pressure, Vol, Flow Induced, FlowElPress', i, REAL(GasPocketNewPress%Array(i)), REAL(GasPocketNewVol%Array(i)), GasPocketFlowInduced%Array(i), FlowEl(GasPocketFlowEl(i , 1))%StartPress
!END DO
!WRITE (*,*) ' BottomHolePress =' , BottomHolePress
!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!IF (ChokeKroneckerDelta == 1) THEN
! WRITE (*,*) ' ChokeLine flowrate' , FlowEl(NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWellToChokeEl)%FlowRate , stringflowrate
! !i = NoHorizontalEl + NoStringEl + NoAnnulusEl
! !j = NoHorizontalEl + NoStringEl + NoAnnulusEl + NoWellToChokeEl
! !WRITE (*,*) ' Well Outlet and Chokeline Outlet Pressure' , FlowEl(i)%EndPress, FlowEl(j)%EndPress
!END IF
!IF (GasPocket%ElementNo == 0) THEN
! KickUnknownVector(2) = BottomHolePress
!!ELSE
!! KickUnknownVector(2) = FlowEl(GasPocket%ElementNo)%StartPress
!END IF
!IF (WellHeadOpen)
! GasPocket%NewPress = KickUnknownVector(2)
!END IF
!WRITE (*,*) 'Ann End'
END SUBROUTINE