# 1 "/mnt/c/Projects/VSIM/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 (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 = 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) * 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 (WellHeadWasOpen == .FALSE. .AND. NoGasPocket > 0 .AND. KickIteration == 1) THEN IF (ChokeKroneckerDelta == 1) THEN ! flow on choke line IF (data%State%FricPressDrop%TotalOpenChokeArea < 0.01 * data%State%Choke%ChokeAreaFullyOpen) THEN WRITE (*,*) 'density , TotalOpenChokeArea' , "data%Equipments%DownHole%Density", data%State%FricPressDrop%TotalOpenChokeArea data%State%FricPressDrop%TotalOpenChokeArea = 0.01 * data%State%Choke%ChokeAreaFullyOpen END IF data%State%FricPressDrop%Kchoke = (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 / 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 / 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 (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 >= 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 (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 (ChokeKroneckerDelta == 1) THEN ! thus we should stabilize solution IF (data%State%FricPressDrop%TotalOpenChokeArea > 0.5 * data%State%Choke%ChokeAreaFullyOpen) THEN KickCorrectionUnderRelaxation = 0.6 ELSE IF (data%State%FricPressDrop%TotalOpenChokeArea > 0.1 * data%State%Choke%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 = 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 = GasPocketFlowEl(NoGasPocket , 1) !WRITE (*,*) 'k, Pocket Press', k, GasPocketOldPress%Array(NoGasPocket) - StandardPress i = data%State%FricPressDrop%AnnulusFirstEl j = data%State%FricPressDrop%OpenholeFirstEl - 1 FlowEl(k)%StartPress = GasPocketOldPress%Array(NoGasPocket) - StandardPress FlowEl(k)%EndPress = GasPocketOldPress%Array(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 !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! BottomHolePress = FlowEl(data%State%FricPressDrop%OpenholeFirstEl)%StartPress !!!!!!!!!!!!!!!!!!!!!! END SUBROUTINE