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