# 1 "/home/admin/SimulationCore2/FluidFlow/kick/Gas_Kick_Calculator.f90"
SUBROUTINE GasKickCalculator
use SimulationVariables
use KickVARIABLESModule
use SimulationVariables !@
Use CReservoirVariables
Use CFormationVariables
USE Fluid_Flow_Startup_Vars
use PressureDisplayVARIABLESModule
USE FricPressDropVarsModule
USE MudSystemVARIABLES
use SimulationVariables !@@@
USE CMudPropertiesVariables
USE CError
USE , INTRINSIC :: IEEE_ARITHMETIC
!! Note: a subject that may be confusing is that when we use gauge pressure and when're using absolute pressure?!
!! all stated pressure are gauge pressure, so I do like this.
!! only when we want to use a state equation (like ideal gas equation), we should use absolute equation and so I do this.
!! Thus gas pocket pressure are all absolute pressure.
IMPLICIT NONE
INTEGER :: i , j , k , l
KickVARIABLES%SolvingEquationError = .FALSE.
KickVARIABLES%KickVandPFunction(:) = 0.d0
KickVARIABLES%KickJacobian(: , :) = 0.d0
!====================================================
! Gas Properties (Methane Gas)
!====================================================
!GasPocketNewTemp%Array(1) = 600
KickVARIABLES%BottomHoleTemperature = 600
KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
KickVARIABLES%KickFluxAvgCompressibility = 0.98d0
KickVARIABLES%K_Aa = (5.8742362 * 10.**(-3) * KickVARIABLES%KickFluxAvgTemperature**1.2288) / (511.1728532 + KickVARIABLES%KickFluxAvgTemperature)
KickVARIABLES%K_Bb = 5.5565586 + (1000.01 / KickVARIABLES%KickFluxAvgTemperature)
KickVARIABLES%K_Cc = 2.47862 - 0.12294 * KickVARIABLES%K_Bb
KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
KickVARIABLES%GasKickViscosity = KickVARIABLES%K_Aa * EXP(KickVARIABLES%K_Bb * KickVARIABLES%GasKickSIDensity**KickVARIABLES%K_Cc)
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!! Calculating compressibility for bottom hole condition
KickVARIABLES%BottomHoleCompressibility = 0.98d0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
KickVARIABLES%GasKickBg = 0.00504 * KickVARIABLES%KickFluxAvgCompressibility * KickVARIABLES%KickFluxAvgTemperature / KickVARIABLES%KickFluxAvgPressure ![bbl/SCF]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!WRITE (*,*) 'Gas Kick Top'
!===> Kick Flow Rate Calculation
IF (KickVARIABLES%FormPressure > KickVARIABLES%BottomHolePress) THEN
KickVARIABLES%KickmdotACoef = 10.0**(-8) * 1.15741d0 * 7.080 * KickVARIABLES%FormPermeability * REAL(KickVARIABLES%KickFormLength) * data%State%GasType(KickVARIABLES%KickGasType)%StDensity / &
(KickVARIABLES%GasKickViscosity * KickVARIABLES%GasKickBg * LOG(10000.0))
!IF (WellHeadOpen .AND. NoGasPocket == 1) KickVARIABLES%KickmdotACoef = (1.d0 + 2.d0) * KickmdotACoef
ELSE
KickVARIABLES%KickmdotACoef = 0.0
END IF
i = data%State%FricPressDrop%StringLastEl
j = data%State%FricPressDrop%OpenholeFirstEl - 1
k = KickVARIABLES%GasPocketFlowEl(1 , 1)
KickVARIABLES%KickmdotBCoef = KickVARIABLES%FormPressure + StandardPress !! - Sum(static and friction pressure loss) of flow elements below gas pocket, see below
IF (KickVARIABLES%FormPressure > KickVARIABLES%BottomHolePress) THEN
!WRITE (*,*) 'k , i, j' , k , i, j
IF (k >= data%State%FricPressDrop%OpenholeFirstEl) THEN ! Bottom of active kick is in openhole
KickVARIABLES%KickmdotBCoef = KickVARIABLES%KickmdotBCoef - (SUM(FlowEl(data%State%FricPressDrop%OpenholeFirstEl : k)%StaticPressDiff)) !+ SUM(FlowEl(j + 1 : GasPocketFlowEl(1 , 1) - 1)%FricPressLoss
!WRITE (*,*) '1 SUM(FlowEl(j + 1 : k)%FricPressLoss', k, SUM(FlowEl(j + 1 : k)%FricPressLoss)
ELSE IF (k < data%State%FricPressDrop%OpenholeFirstEl) THEN ! bottom of 1st gas pocket (active kick) is in annulus ond/or choke line only
KickVARIABLES%KickmdotBCoef = KickVARIABLES%KickmdotBCoef - SUM(FlowEl(data%State%FricPressDrop%OpenholeFirstEl : data%State%FricPressDrop%NumbEl)%StaticPressDiff) &
- (SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : k)%StaticPressDiff) + SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : k)%FricPressLoss))
END IF
END IF
DO l = 1 , KickVARIABLES%NoGasPocket
KickVARIABLES%KickUnknownVector(2 * l - 1) = GasPocketNewVol%Array(l)
KickVARIABLES%KickUnknownVector(2 * l) = GasPocketNewPress%Array(l)
END DO
IF (KickVARIABLES%WellHeadOpen) THEN
!!!!!!!!!! Calculation of functions of pocket Pressure and gas Volumes
!IF (GasPocketElementNo(1) > 0) WRITE (*,*) ' GasPocketElementNo(1) ' , GasPocketElementNo(1)
!WRITE (*,*) ' Kick Unknown Vector' , KickUnknownVector!(1) , KickUnknownVector(2)
IF (KickIteration == 1) THEN ! updating initial guess based on previous time step data
DO l = 1 , KickVARIABLES%NoGasPocket
KickVARIABLES%KickUnknownVector(2 * l - 1) = KickVARIABLES%KickUnknownVector(2 * l - 1) + GasPocketDeltaVol%Array(l)
END DO
END IF
KickVARIABLES%KickVandPFunction(1) = KickVARIABLES%KickUnknownVector(1) - GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * & ! VandP(1) = V(1)
GasPocketNewTemp%Array(1) * (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * MAX(((KickVARIABLES%KickmdotBCoef - KickVARIABLES%KickUnknownVector(2)) * dt) , 0.0)) / KickVARIABLES%KickUnknownVector(2)
!WRITE (*,*) 'KickVandPFunction(1)',KickVandPFunction(1)
l = 2 * KickVARIABLES%NoGasPocket
CALL KickFunctionsCalculator(KickVARIABLES%KickVandPFunction(l) , KickVARIABLES%NoGasPocket , 2) ! VandP(last) = P(last)
!WRITE (*,*) 'KickVandPFunction(l)', l, KickVandPFunction(l)
DO l = 2 , KickVARIABLES%NoGasPocket ! VandP(Odd) = V(l, l > 1)
KickVARIABLES%KickVandPFunction(2 * l - 1) = KickVARIABLES%KickUnknownVector(2 * l - 1) - GasPocketCompressibility%Array(l) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * &
GasPocketNewTemp%Array(l) * GasPocketWeight%Array(l) / KickVARIABLES%KickUnknownVector(2 * l)
!WRITE(*,*) 'KickVandPFunction(V)', l, KickVandPFunction(2 * l - 1)
END DO
DO l = KickVARIABLES%NoGasPocket - 1 , 1 , -1
CALL KickFunctionsCalculator(KickVARIABLES%KickVandPFunction(2 * l) , l , 1) ! VandP(Even) = P(l, l < KickVARIABLES%NoGasPocket)
!WRITE(*,*) 'KickVandPFunction(P)', l , KickVandPFunction(2 * l)
END DO
!!!!!!!!!! END - Calculation of functions of pocket Pressure and gas Volumes
!!!!!!!!!! Calculation of Jacobian
DO k = 1 , 2 * KickVARIABLES%NoGasPocket ! Main Diagonal
KickVARIABLES%KickJacobian(k , k) = 1.d0
END DO
KickVARIABLES%KickJacobian(1,2) = (GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * GasPocketNewTemp%Array(1) &
* (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * KickVARIABLES%KickmdotBCoef * dt) / KickVARIABLES%KickUnknownVector(2)**2) ! Row 1 Finished
IF (KickVARIABLES%KickJacobian(1,2) == 0.d0) THEN
CALL Error('KickJacobian(1,2) = 0.0')
KickVARIABLES%KickJacobian(1,2) = KickVARIABLES%OldKickJacobian(1,2)
END IF
!WRITE(*,*) 'KickJacobian(1,2)', KickJacobian(1,2)
l = KickVARIABLES%NoGasPocket
CALL KickFunctionsCalculator(KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) , KickVARIABLES%NoGasPocket , 4) ! for last Row
IF (KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) == 0.d0) THEN
CALL Error ('KickJacobian(Last,Odd) = 0.0')
KickVARIABLES%KickJacobian(2 * l , 2 * l - 1) = KickVARIABLES%OldKickJacobian(2 * l , 2 * l - 1)
END IF
DO k = KickVARIABLES%NoGasPocket - 1 , 1 , -1
KickVARIABLES%KickJacobian(2 * l , 2 * k - 1) = KickVARIABLES%KickJacobian(2 * l , 2 * l - 1)
END DO ! Last Row Finished
!WRITE(*,*) 'KickJacobian(2,1)', KickJacobian(2,1)
DO k = 2 , KickVARIABLES%NoGasPocket
KickVARIABLES%KickJacobian(2 * k - 1 , 2 * k) = GasPocketCompressibility%Array(k) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * GasPocketNewTemp%Array(k) &
* GasPocketWeight%Array(k) / KickVARIABLES%KickUnknownVector(2 * k)**2
END DO ! Odd Rows (V equations) Finished
DO k = 1 , KickVARIABLES%NoGasPocket - 1
KickVARIABLES%KickJacobian(2 * k , 2 * k + 2) = -1.d0
END DO ! Even Rows (P equations) effect of upper pocket
DO k = 2 , 2 * (KickVARIABLES%NoGasPocket - 1) , 2
DO l = 1 , k - 1 , 2
CALL KickFunctionsCalculator(KickVARIABLES%KickJacobian(k , l) , k / 2 , 3)
IF (KickVARIABLES%KickJacobian(k , l) == 0.d0) THEN
WRITE (*,*) 'Jacobian Array = 0.0', k , l
CALL Error ('KickJacobian(k , l) = 0.0')
KickVARIABLES%KickJacobian(k , l) = KickVARIABLES%OldKickJacobian(k , l)
END IF
END DO
END DO
IF (ANY(IEEE_IS_NaN(KickVARIABLES%KickJacobian))) CALL ErrorStop ('NaN in calculating Kick Jacobian, Call your Service Provider')
!!!!!!!!!! Solving linear equation in order to finding correction vector for correcting pocket pressure and gas induced flowrates
KickVARIABLES%KickVandPFunction = -1.d0 * KickVARIABLES%KickVandPFunction
CALL SOLVE_LINEAR_EQUATIONS(KickVARIABLES%KickJacobian , KickVARIABLES%KickCorrectionVector , KickVARIABLES%KickVandPFunction , KickVARIABLES%SolvingEquationError, SIZE(KickVARIABLES%KickCorrectionVector))
IF (KickVARIABLES%SolvingEquationError) CALL ErrorStop( ' Error in solving kick equation ' )
KickVARIABLES%KickUnknownVector = KickVARIABLES%KickUnknownVector + KickVARIABLES%KickCorrectionUnderRelaxation * KickVARIABLES%KickCorrectionVector
DO l = 1 , KickVARIABLES%NoGasPocket
GasPocketNewVol%Array(l) = KickVARIABLES%KickUnknownVector(2 * l - 1)
IF (IEEE_IS_NaN(GasPocketNewVol%Array(l))) CALL ErrorStop('Volume of this pocket is Not a Number:', l)
IF (GasPocketNewVol%Array(l) <= 0.d0) CALL Error('Volume of this pocket is Negative or Zero:', l)
GasPocketNewPress%Array(l) = KickVARIABLES%KickUnknownVector(2 * l)
IF (IEEE_IS_NaN(GasPocketNewPress%Array(l))) CALL ErrorStop('Pressure of this Pocket is Not a Number:', l)
IF (GasPocketNewPress%Array(l) <= 0.d0) CALL ErrorStop('Pressure of this Pocket is Negative or Zero:', l)
END DO
ELSE ! well haed is closed, so build up process or migration occurs
!WRITE (*,*) 'GasPocketOldPress (before)' , GasPocketOldPress(1)
GasPocketNewPress%Array(1) = GasPocketOldPress%Array(1) * &
(REAL((GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * KickVARIABLES%KickmdotBCoef * dt) / (GasPocketWeight%Array(1) + KickVARIABLES%KickmdotACoef * GasPocketOldPress%Array(1) * dt)))
END IF
!DO l = 1 , KickVARIABLES%NoGasPocket
GasPocketDeltaVol%Array(:) = GasPocketNewVol%Array(:) - GasPocketOldVol%Array(:)
GasPocketFlowInduced%Array(:) = (GasPocketDeltaVol%Array(:)) / dt * 448.8 ! gpm
!END DO
KickVARIABLES%GasKickPumpFlowRate = 0.0
IF (NOT(KickVARIABLES%KickOffBottom) .AND. KickVARIABLES%WellHeadOpen) KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
END SUBROUTINE
SUBROUTINE KickFunctionsCalculator(ExitValue , GasPocketNo , CalcMode)
use SimulationVariables
use KickVARIABLESModule
USE FricPressDropVarsModule
USE Fluid_Flow_Startup_Vars
USE CError
USE , INTRINSIC :: IEEE_Arithmetic
IMPLICIT NONE
REAL(8) :: ExitValue
INTEGER :: GasPocketNo , CalcMode
INTEGER :: x
INTEGER :: y
INTEGER :: z , i , j
x = KickVARIABLES%GasPocketFlowEl(GasPocketNo , 1)
IF (GasPocketNo < KickVARIABLES%NoGasPocket) y = KickVARIABLES%GasPocketFlowEl(GasPocketNo + 1 , 1)
i = data%State%FricPressDrop%StringLastEl
j = data%State%FricPressDrop%OpenholeFirstEl - 1
! Case 1: gas pocket is completely in OP and STARTX of upper gas pocket is also
! Case 2: gas pocket is completely in OP and STARTX of upper gas pocket is above Bit
! Case 3: gas pocket is AROUNDBIT and so upper gas pocket is in ANN (or Choke line)
! Case 4: gas pocket is completely in ANN and upper gas pocket is also
! CalcMode 1: KickVandPFunction between 2 pocket
! CalcMode 2: KickVandPFunction for top gas pocket
! CalcMode 3: KickJacobian between 2 Pocket
! CalcMode 4: KickJacobian for top (last) gas pocket
! CalcMode 5: Calculating pressure of bottom pocket based on upper pocket
IF (CalcMode == 1) THEN ! calculating pressure difference between two pocket, include static pressure difference and frictional
! pressure difference, use in calculating 'KickVandPFunction'
ExitValue = KickVARIABLES%KickUnknownVector(2 * GasPocketNo) - KickVARIABLES%KickUnknownVector(2 * GasPocketNo + 2)
IF (x >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%OpenholeFirstEl) THEN ! Case 2 , Case 3
ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricPressLoss) &
- SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : y)%StaticPressDiff) - SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : y)%FricPressLoss)
ELSE ! Case 1 , Case 4
ExitValue = ExitValue - SUM(FlowEl(x : y)%StaticPressDiff) - SUM(FlowEl(x : y)%FricPressLoss)
END IF
ELSE IF (CalcMode == 2) THEN
ExitValue = KickVARIABLES%KickUnknownVector(2 * GasPocketNo) - StandardPress - data%State%FricPressDrop%Kchoke * FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%Flowrate**2
IF (x >= data%State%FricPressDrop%OpenholeFirstEl) THEN ! Gas Pocket is in Openhole
ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricPressLoss) &
- SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricPressLoss)
ELSE ! Gas Pocket is in Annulus
ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%OpenholeFirstEl - 1)%StaticPressDiff) - SUM(FlowEl(x : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricPressLoss)
END IF
ELSE IF (CalcMode == 3) THEN ! calculating derivative of pressure difference between two pocket, relative to change in flowrate
! use in calculating 'KickJacobian'
IF (x >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%OpenholeFirstEl) THEN ! Top kick STARTX is in Annulus
DO z = x , data%State%FricPressDrop%NumbEl ! open hole elements
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
DO z = data%State%FricPressDrop%AnnulusFirstEl , y ! Annulus elements
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
ExitValue = ExitValue - (SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricToQPartialDiff) + SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : y)%FricToQPartialDiff)) * 448.8 / dt
ELSE ! both pockets are one side of bit
DO z = x , y
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
ExitValue = ExitValue - SUM(FlowEl(x : y)%FricToQPartialDiff) * 448.8 / dt
END IF
ELSE IF (CalcMode == 4) THEN ! partial derivative of frictional pressure drop relative to flowrate for top gas pocket
ExitValue = - 2.d0 * data%State%FricPressDrop%Kchoke * FlowEl(data%State%FricPressDrop%OpenholeFirstEl - 1)%Flowrate * 448.8 / dt
IF (x >= data%State%FricPressDrop%OpenholeFirstEl) THEN ! kick STARTX is in openhole
DO z = x , data%State%FricPressDrop%NumbEl ! open hole elements
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
DO z = data%State%FricPressDrop%AnnulusFirstEl , data%State%FricPressDrop%OpenholeFirstEl - 1 ! Annulus elements
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
ExitValue = ExitValue - (SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%FricToQPartialDiff) + SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricToQPartialDiff)) * 448.8 / dt
ELSE
DO z = x , data%State%FricPressDrop%OpenholeFirstEl - 1 ! Annulus elements
CALL PartialDerivativeFricToFlowRate(z)
IF (IEEE_IS_NaN(FlowEl(z)%FricToQPartialDiff)) THEN
WRITE (*,*) ' FricToQPartialDiff , GenRe ' , x , FlowEl(z)%FricToQPartialDiff , FlowEl(z)%GenRe
WRITE (*,*) ' Op start, end, density, Q, mu' , FlowEl(z)%StartX, FlowEl(z)%EndX, FlowEl(z)%Density, FlowEl(z)%FlowRate, FlowEl(z)%mueff
CALL ErrorStop('NaN in calculating partial derivative')
END IF
END DO
ExitValue = ExitValue - SUM(FlowEl(x : data%State%FricPressDrop%OpenholeFirstEl - 1)%FricToQPartialDiff) * 448.8 / dt
END IF
ELSE IF (CalcMode == 5) THEN
IF (x >= data%State%FricPressDrop%OpenholeFirstEl .AND. y < data%State%FricPressDrop%OpenholeFirstEl) THEN ! Gas Pocket is in Openhole and upper pocket is in annulus
!WRITE (*,*) 'x , y 1' , x, y
ExitValue = GasPocketNewPress%Array(GasPocketNo + 1) + SUM(FlowEl(x : data%State%FricPressDrop%NumbEl)%StaticPressDiff) + SUM(FlowEl(data%State%FricPressDrop%AnnulusFirstEl : y)%StaticPressDiff)
ELSE ! Both gas pockets are in Annulus or openhole
!WRITE (*,*) 'x , y 2' , x, y
ExitValue = GasPocketNewPress%Array(GasPocketNo + 1) + SUM(FlowEl(x : y)%StaticPressDiff)
END IF
END IF
END SUBROUTINE
SUBROUTINE NewGasKick
use SimulationVariables
use KickVARIABLESModule
use SimulationVariables !@
Use CReservoirVariables
Use CFormationVariables
USE Fluid_Flow_Startup_Vars
use PressureDisplayVARIABLESModule
USE FricPressDropVarsModule
USE MudSystemVARIABLES
use SimulationVariables !@@@
USE CMudPropertiesVariables
USE CError
USE , INTRINSIC :: IEEE_ARITHMETIC
!! Note: a subject that may be confusing is that when we use gauge pressure and when using absolute pressure?!
!! all stated pressure are gauge pressure, so I do like this.
!! only when we want to use a state equation (like ideal gas equation), we should use absolute equation and so I do this.
!! Thus gas pocket pressure are all absolute pressure.
IMPLICIT NONE
INTEGER :: i , j , k , l
IF (NOT(ALLOCATED(GasPocketWeight%Array))) THEN ! 1st kick
WRITE (*,*) ' New Influx 1'
KickVARIABLES%NoGasPocket = 1
data%State%MudSystem%NewInfluxNumber = data%State%MudSystem%NewInfluxNumber + 1
data%State%MudSystem%NewInfluxElementCreated = 0
KickVARIABLES%KickOffBottom = .FALSE.
CALL GasPocketOldPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
CALL GasPocketNewPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
CALL GasPocketOldTemp%AddToFirst(600.0)
CALL GasPocketNewTemp%AddToFirst(600.0)
CALL GasPocketOldVol%AddToFirst(0.d0)
CALL GasPocketNewVol%AddToFirst(0.d0)
CALL GasPocketdeltaVol%AddToFirst(0.0)
CALL GasPocketFlowInduced%AddToFirst(0.0)
CALL GasPocketModifiedVol%AddToFirst(0.0)
CALL GasPocketWeight%AddToFirst(0.0)
CALL GasPocketDensity%AddToFirst(2.0)
CALL GasPocketCompressibility%AddToFirst(0.98)
ALLOCATE(KickVARIABLES%KickJacobian(2 , 2) , KickVARIABLES%OldKickJacobian(2 , 2) , KickVARIABLES%KickVandPFunction(2) , KickVARIABLES%KickUnknownVector(2) , KickVARIABLES%KickCorrectionVector(2))
KickVARIABLES%BottomHoleTemperature = 600
KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
KickVARIABLES%KickFluxAvgCompressibility = 0.98
KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
KickVARIABLES%GasKickDensity = KickVARIABLES%GasKickSIDensity * 8.3523
GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * KickVARIABLES%MinKickVol !1.0:seyyed gofte !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * &
GasPocketNewTemp%Array(1) * GasPocketWeight%Array(1) / GasPocketNewPress%Array(1)
GasPocketDeltaVol%Array(1) = 0.05 !GasPocketNewVol%Array(1)
GasPocketFlowInduced%Array(1) = (GasPocketDeltaVol%Array(1)) / dt * 448.8 ! gpm
KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
WRITE (*,*) ' FormPressure , BottomHolePress' , KickVARIABLES%FormPressure , KickVARIABLES%BottomHolePress, KickVARIABLES%GasKickDensity
WRITE (*,*) ' No Press(psia) Vol(gal) Weight(lbm) Flow Induced(gpm)'
DO i = 1 , KickVARIABLES%NoGasPocket
WRITE (*,102) i , GasPocketNewPress%Array(i), GasPocketNewVol%Array(i) * Convft3toUSgal, GasPocketWeight%Array(i), GasPocketFlowInduced%Array(i)
END DO
!ELSE IF (NoGasPocket < MaxGasPocket .AND. KickVARIABLES%KickOffBottom .AND. (GasPocketNewVol%Array(1) > KickVARIABLES%MinAllowableKickVol .OR. KickWasExitingThroughChoke)) THEN
ELSE IF (KickVARIABLES%NoGasPocket < KickVARIABLES%MaxGasPocket .AND. KickVARIABLES%KickOffBottom .AND. (GasPocketNewVol%Array(1) > KickVARIABLES%MinAllowableKickVol .OR. ANY(KickVARIABLES%GasPocketFlowEl(1 , :) == data%State%FricPressDrop%OpenholeFirstEl - 1))) THEN
WRITE (*,*) ' New Influx', KickVARIABLES%NoGasPocket + 1
102 FORMAT (I2, 4X, (F8.1), 3X, (F8.3), 2X, (F8.3), 8X, (F8.3))
KickVARIABLES%NoGasPocket = KickVARIABLES%NoGasPocket + 1
data%State%MudSystem%NewInfluxNumber = data%State%MudSystem%NewInfluxNumber + 1
data%State%MudSystem%NewInfluxElementCreated = 0
KickVARIABLES%KickOffBottom = .FALSE.
CALL GasPocketOldPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
CALL GasPocketNewPress%AddToFirst((KickVARIABLES%BottomHolePress + StandardPress) * 1.d0)
CALL GasPocketOldTemp%AddToFirst(600.0)
CALL GasPocketNewTemp%AddToFirst(600.0)
CALL GasPocketOldVol%AddToFirst(0.d0)
CALL GasPocketNewVol%AddToFirst(0.d0)
CALL GasPocketdeltaVol%AddToFirst(0.0)
CALL GasPocketFlowInduced%AddToFirst(0.0)
CALL GasPocketModifiedVol%AddToFirst(0.0)
CALL GasPocketWeight%AddToFirst(0.0)
CALL GasPocketDensity%AddToFirst(2.0)
CALL GasPocketCompressibility%AddToFirst(0.98)
DEALLOCATE(KickVARIABLES%KickJacobian , KickVARIABLES%OldKickJacobian , KickVARIABLES%KickVandPFunction , KickVARIABLES%KickUnknownVector , KickVARIABLES%KickCorrectionVector)
ALLOCATE(KickVARIABLES%KickJacobian(2 * KickVARIABLES%NoGasPocket , 2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%OldKickJacobian(2 * KickVARIABLES%NoGasPocket , 2 * KickVARIABLES%NoGasPocket))
ALLOCATE(KickVARIABLES%KickUnknownVector(2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%KickCorrectionVector(2 * KickVARIABLES%NoGasPocket) , KickVARIABLES%KickVandPFunction(2 * KickVARIABLES%NoGasPocket))
KickVARIABLES%BottomHoleTemperature = 600
KickVARIABLES%KickFluxAvgPressure = (KickVARIABLES%BottomHolePress + KickVARIABLES%FormPressure) / 2 + StandardPress
KickVARIABLES%KickFluxAvgTemperature = (KickVARIABLES%FormTemperature + KickVARIABLES%BottomHoleTemperature) / 2
KickVARIABLES%KickFluxAvgCompressibility = 0.98
KickVARIABLES%GasKickSIDensity = KickVARIABLES%KickFluxAvgPressure / (KickVARIABLES%KickFluxAvgCompressibility * &
KickVARIABLES%KickFluxAvgTemperature * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant) * Convpcftogpcm3
KickVARIABLES%GasKickDensity = KickVARIABLES%GasKickSIDensity * 8.3523
!GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * 0.05 !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
GasPocketWeight%Array(1) = KickVARIABLES%GasKickDensity * KickVARIABLES%MinKickVol !1.0:seyyed gofte !KickmdotACoef * (KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
GasPocketNewVol%Array(1) = GasPocketCompressibility%Array(1) * data%State%GasType(KickVARIABLES%KickGasType)%GasConstant * &
GasPocketNewTemp%Array(1) * GasPocketWeight%Array(1) / GasPocketNewPress%Array(1)
GasPocketDeltaVol%Array(1) = 0.05 !GasPocketNewVol%Array(1)
GasPocketFlowInduced%Array(1) = (GasPocketDeltaVol%Array(1)) / dt * 448.8 ! gpm
KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
WRITE (*,*) ' FormPressure , BottomHolePress' , KickVARIABLES%FormPressure , KickVARIABLES%BottomHolePress, KickVARIABLES%GasKickDensity
WRITE (*,*) ' No Press(psia) Vol(gal) Weight(lbm) Flow Induced(gpm)'
DO i = 1 , KickVARIABLES%NoGasPocket
WRITE (*,102) i , GasPocketNewPress%Array(i), GasPocketNewVol%Array(i) * Convft3toUSgal, GasPocketWeight%Array(i), GasPocketFlowInduced%Array(i)
END DO
ELSE ! no new kick, so mass of 1st kick should increase
GasPocketWeight%Array(1) = GasPocketweight%Array(1) + KickVARIABLES%KickmdotACoef * (KickVARIABLES%KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
KickVARIABLES%GasKickPumpFlowRate = GasPocketFlowInduced%Array(1)
IF (KickVARIABLES%NoGasPocket > 1 .OR. KickVARIABLES%SecondaryKickWeight > 0.0) THEN
KickVARIABLES%SecondaryKickWeight = KickVARIABLES%SecondaryKickWeight + KickVARIABLES%KickmdotACoef * (KickVARIABLES%KickmdotBCoef - GasPocketNewPress%Array(1)) * dt
KickVARIABLES%SecondaryKickVol = KickVARIABLES%SecondaryKickWeight / KickVARIABLES%GasReservoirDensity / 42.0 ! 42 USGal = 1bbl
END IF
END IF
END SUBROUTINE