SimulationCore2/FluidFlow/Frictional_Press_Drop_Calc.i90

# 1 "/mnt/c/Projects/VSIM/SimulationCore2/FluidFlow/Frictional_Press_Drop_Calc.f90"
SUBROUTINE FricPressDrop(iloc)

    
    
!!      Record of revisions
!!          Date            Programmer          Discription of change
!!         ------          ------------        -----------------------
!!        1396/07/23        Sheikh              Original code
!!
    
    
    
    USE FricPressDropVarsModule
    USE SimulationVariables
    USE Fluid_Flow_Startup_Vars
    USE CError
    
    IMPLICIT NONE    

    INTEGER :: iloc
    REAL :: TauZero
    TauZero = 12.0
!ActiveRheologyModel = Herschel_Bulkley_RheologyModel

! 0 = Power Law   ,   1 = Bingham Plastic   ,   2 = Newtonian
!TotFricPressLoss = 0.0
    
        FlowEl(iloc)%alpha = 1                ! assume that all elements have annulus geometry
        FlowEl(iloc)%dPdLfric = 0.0
        FlowEl(iloc)%f = 0.0
        FlowEl(iloc)%FlowRate = ABS(FlowEl(iloc)%FlowRate)
        
        
    
    IF ((FlowEl(iloc)%FlowRate >= PressFlowrateTolerance) &
        .AND. (FlowEl(iloc)%MaterialType /= 1) &        ! not gas kick
        .AND. (ABS(FlowEl(iloc)%Length) >= PressLengthTolerance) &
        .AND. (FlowEl(iloc)%MaterialType /= 4)) THEN    ! not air
        
        IF (FlowEl(iloc)%Id==0) THEN
            FlowEl(iloc)%alpha = 0
        END IF
        
        FlowEl(iloc)%muPlastic = FlowEl(iloc)%Theta600 - FlowEl(iloc)%Theta300   ! cp
        FlowEl(iloc)%YieldP = 2.0 * FlowEl(iloc)%Theta300 - FlowEl(iloc)%Theta600   ! lbf/100ft**2
        FlowEl(iloc)%nIndex = 3.32 * log10(FlowEl(iloc)%Theta600 / FlowEl(iloc)%Theta300)
        FlowEl(iloc)%kIndex = 510.0 * FlowEl(iloc)%Theta300 / (511.0**FlowEl(iloc)%nIndex)    ! rabete fv2
        IF (data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel .AND. FlowEl(iloc)%alpha == 0) THEN
            FlowEl(iloc)%kIndex = 1.066 * FlowEl(iloc)%Theta300 / (511.0**FlowEl(iloc)%nIndex)
        ELSE IF (data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel .AND. FlowEl(iloc)%alpha == 1) THEN
            FlowEl(iloc)%nIndex = 3.32 * log10((FlowEl(iloc)%Theta600 - TauZero) / (FlowEl(iloc)%Theta300 - TauZero))
            FlowEl(iloc)%kIndex = 1.066 * (FlowEl(iloc)%Theta300 - TauZero) / (511.0**FlowEl(iloc)%nIndex)
        END IF
        
        
        
        
! Calculating velocity
        FlowEl(iloc)%vel = 0.408 * FlowEl(iloc)%FlowRate / (FlowEl(iloc)%Od**2 - FlowEl(iloc)%Id**2)   ! velocity in ft/s
!FlowEl(iloc)%vel = 24.51 * FlowEl(iloc)%FlowRate / (FlowEl(iloc)%Od**2 - FlowEl(iloc)%Id**2)   ! velocity in ft/min

!IF (FlowModel == Bingham_RheologyModel) THEN   ! Bingham Plastic
!    FlowEl(iloc)%Gf = (2. + FlowEl(iloc)%alpha) / 2.
!ELSE IF (FlowModel == PowerLow_RheologyModel) THEN
!    FlowEl(iloc)%Gf = ((3. - FlowEl(iloc)%alpha) * FlowEl(iloc)%nIndex + 1.) / FlowEl(iloc)%nIndex / (4. - FlowEl(iloc)%alpha) * (2. + FlowEl(iloc)%alpha) / 2.
!END IF
        
!FlowEl(iloc)%gammaW = 1.6 * FlowEl(iloc)%Gf * FlowEl(iloc)%vel / FlowEl(iloc)%Dhyd
        
!IF (FlowModel == Bingham_RheologyModel) THEN   ! Bingham Plastic
!    FlowEl(iloc)%tauW = 1.067 * ((4. - FlowEl(iloc)%alpha) / (3. - FlowEl(iloc)%alpha) * FlowEl(iloc)%YieldP + FlowEl(iloc)%muPlastic * FlowEl(iloc)%gammaW)
!    !FlowEl(iloc)%tauW = 1.067*(FlowEl(iloc)%YieldP+FlowEl(iloc)%muPlastic*FlowEl(iloc)%gammaW)
!ELSE IF (FlowModel == PowerLow_RheologyModel) THEN     ! Power law
!    FlowEl(iloc)%tauW = 1.067 * FlowEl(iloc)%kIndex * FlowEl(iloc)%gammaW**FlowEl(iloc)%nIndex
!END IF
        
! Calculating effective or apparent viscosity
        IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic
            FlowEl(iloc)%mueff = FlowEl(iloc)%muPlastic + 5. * FlowEl(iloc)%YieldP * FlowEl(iloc)%Dhyd / FlowEl(iloc)%vel
!write(*,*) 'pointer1' , FlowEl(iloc)%muPlastic ,  FlowEl(iloc)%YieldP , FlowEl(iloc)%Dhyd , FlowEl(iloc)%vel
        ELSE IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel .OR. data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel) THEN   ! Power Law
            FlowEl(iloc)%Gf = ((3. - FlowEl(iloc)%alpha) * FlowEl(iloc)%nIndex + 1.0) / FlowEl(iloc)%nIndex / (4.0 - FlowEl(iloc)%alpha) * (2.0 + FlowEl(iloc)%alpha) / 2.0
            FlowEl(iloc)%mueff = (FlowEl(iloc)%kIndex) / (1. + FlowEl(iloc)%alpha / 2.) * ((96. * FlowEl(iloc)%vel / FlowEl(iloc)%Dhyd)**(FlowEl(iloc)%nIndex - 1)) * FlowEl(iloc)%Gf**FlowEl(iloc)%nIndex
!write(*,*) 'pointer2' , FlowEl(iloc)%kIndex ,FlowEl(iloc)%alpha , FlowEl(iloc)%vel ,FlowEl(iloc)%Dhyd,FlowEl(iloc)%nIndex ,FlowEl(iloc)%Gf ,FlowEl(iloc)%nIndex
        END IF
        
        FlowEl(iloc)%gammaW = 96.0 * FlowEl(iloc)%Gf * FlowEl(iloc)%vel / FlowEl(iloc)%Dhyd
        
        FlowEl(iloc)%tauW = ((4.0 - FlowEl(iloc)%alpha) / (3.0 - FlowEl(iloc)%alpha))**FlowEl(iloc)%nIndex * TauZero + FlowEl(iloc)%kIndex * FlowEl(iloc)%gammaW**FlowEl(iloc)%nIndex
        
        
! Calculating Reynolds number
        IF (FlowEl(iloc)%Od == FlowEl(iloc)%Dhyd) THEN
            FlowEl(iloc)%GenRe = 928. * FlowEl(iloc)%density * FlowEl(iloc)%vel * FlowEl(iloc)%Dhyd / FlowEl(iloc)%mueff
        ELSE
            FlowEl(iloc)%GenRe = 757. * FlowEl(iloc)%density * FlowEl(iloc)%vel * FlowEl(iloc)%Dhyd / FlowEl(iloc)%mueff
        END IF
        
!FlowEl(iloc)%GenRe = 2997 * FlowEl(iloc)%density * FlowEl(iloc)%vel**2 / 19.36 / FlowEl(iloc)%tauW
        
! Calculating friction factor
        IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic
            IF (FlowEl(iloc)%GenRe <= 2000.0) THEN   ! laminar regime
                FlowEl(iloc)%f = 16.0 / FlowEl(iloc)%GenRe
            ELSE IF (FlowEl(iloc)%GenRe >= 4000.0) THEN   ! turbulent regime
                FlowEl(iloc)%a = 0.0791
                FlowEl(iloc)%b = 0.25
                FlowEl(iloc)%f = FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b
            ELSE    !! transition from laminar to turbulent regime
                FlowEl(iloc)%a = 0.0791
                FlowEl(iloc)%b = 0.25
                FlowEl(iloc)%f = (4000.0 - FlowEl(iloc)%GenRe) / 2000.0 * 16. / FlowEl(iloc)%GenRe &
                                + (FlowEl(iloc)%GenRe - 2000.0) / 2000.0 * FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b
            END IF
            
        ELSE IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel) THEN     ! Power law
            FlowEl(iloc)%ReCritLam = 3470. - 1370. * FlowEl(iloc)%nIndex
            FlowEl(iloc)%ReCritTurb = 4270. - 1370. * FlowEl(iloc)%nIndex
            
            IF (FlowEl(iloc)%GenRe <= FlowEl(iloc)%ReCritLam) THEN   ! laminar regime
                FlowEl(iloc)%f = 16.0 / FlowEl(iloc)%GenRe / (1 - 0.184 * FlowEl(iloc)%alpha)
            ELSE IF (FlowEl(iloc)%GenRe >= FlowEl(iloc)%ReCritTurb) THEN   ! turbulent regime
                FlowEl(iloc)%a = (log10(FlowEl(iloc)%nIndex) + 3.93) / 50.
                FlowEl(iloc)%b = (1.75 - log10(FlowEl(iloc)%nIndex)) / 7.
                FlowEl(iloc)%f = FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b
            ELSE
                FlowEl(iloc)%a = (log10(FlowEl(iloc)%nIndex) + 3.93) / 50.
                FlowEl(iloc)%b = (1.75 - log10(FlowEl(iloc)%nIndex)) / 7.
                FlowEl(iloc)%f = (FlowEl(iloc)%ReCritTurb - FlowEl(iloc)%GenRe) / 800.0 * 16. / FlowEl(iloc)%GenRe &
                                + (FlowEl(iloc)%GenRe - FlowEl(iloc)%ReCritLam) / 800.0 * FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b
            END IF
            
        END IF
        
!WRITE (*,*) 'fric press drop', iloc
!WRITE (*,*) 'Length', ABS(REAL(FlowEl(iloc)%Length))
!WRITE (*,*) 'FlowRate', FlowEl(iloc)%FlowRate
!WRITE (*,*) 'Theta600 , Theta300', FlowEl(iloc)%Theta600 , FlowEl(iloc)%Theta300
!WRITE (*,*) 'Dhyd', FlowEl(iloc)%Dhyd
!WRITE (*,*) 'GenRe', FlowEl(iloc)%GenRe
!WRITE (*,*) 'f', FlowEl(iloc)%f
        
        
        
        
    END IF

! Frictional pressure loss gradient calculation
! FlowEl(iloc)%dPdLfric = 1.076 * FlowEl(iloc)%f * FlowEl(iloc)%vel**2 * FlowEl(iloc)%density / 10**5 / FlowEl(iloc)%Dhyd
        FlowEl(iloc)%dPdLfric = FlowEl(iloc)%f * (FlowEl(iloc)%vel)**2 * FlowEl(iloc)%density / 25.81 / FlowEl(iloc)%Dhyd
        FlowEl(iloc)%FricPressLoss = FlowEl(iloc)%dPdLfric * ABS(REAL(FlowEl(iloc)%Length))
        
        IF (FlowEl(iloc)%FrictionDirection == -1) THEN
            FlowEl(iloc)%FlowRate = - FlowEl(iloc)%FlowRate
            FlowEl(iloc)%dPdLfric = - FlowEl(iloc)%dPdLfric
            FlowEl(iloc)%FricPressLoss = - FlowEl(iloc)%FricPressLoss
        END IF
        
        
            
!END DO


END SUBROUTINE FricPressDrop
    
    
SUBROUTINE PartialDerivativeFricToFlowRate(iloc)
    
    USE FricPressDropVarsModule
    USE SimulationVariables
    USE Fluid_Flow_Startup_Vars
    use KickVARIABLESModule
    USE CError
    
    
    IMPLICIT NONE    

    INTEGER :: iloc
    
    FlowEl(iloc)%FricToQPartialDiff = 0.0
!FlowEl(iloc)%FlowRate = ABS(FlowEl(iloc)%FlowRate)
    
    
    IF ((ABS(FlowEl(iloc)%FlowRate) >= PressFlowrateTolerance) &
        .AND. (FlowEl(iloc)%MaterialType /= 1) &        ! not gas kick
        .AND. (ABS(FlowEl(iloc)%Length) >= PressLengthTolerance) &
        .AND. (FlowEl(iloc)%MaterialType /= 4)) THEN    ! not air
        
        IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel) THEN     ! Power law
!IF (FlowEl(iloc)%Flowrate == 0.0) THEN
!    FlowEl(iloc)%Flowrate = 10.0
!    CALL FricPressDrop(iloc)
!END IF
            IF (FlowEl(iloc)%GenRe <= FlowEl(iloc)%ReCritLam) THEN     ! laminar flow
                FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate * FlowEl(iloc)%nIndex
                
            ELSE IF (FlowEl(iloc)%GenRe >= FlowEl(iloc)%ReCritTurb) THEN       ! turbulent flow
                FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &
                                                        * (2. - FlowEl(iloc)%b * (2. - FlowEl(iloc)%nIndex))
            
            ELSE        ! transition from laminar to turbulent
                FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &
                                                        * (2. + (2. - FlowEl(iloc)%nIndex) &
                                                            * ((FlowEl(iloc)%a * FlowEl(iloc)%GenRe**(1. - FlowEl(iloc)%b) - 16.) / 800. / FlowEl(iloc)%f - 1.))
            END IF
        
        ELSE IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic
            IF (FlowEl(iloc)%GenRe <= 2000.0 .OR. FlowEl(iloc)%f == 0.0) THEN     ! laminar flow    if f = 0.0, we have no flow in first time flowing
                FlowEl(iloc)%FricToQPartialDiff = (16. * FlowEl(iloc)%muPlastic * REAL(FlowEl(iloc)%Length) * 2.224 * (10.)**(-3)) &
                                                        / (25.81 * 928. * (1 - 0.184 * FlowEl(iloc)%alpha) * FlowEl(iloc)%Dhyd**2 * FlowEl(iloc)%Area)
        
            ELSE IF (FlowEl(iloc)%GenRe >= 4000.0) THEN       ! turbulent flow
                FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &
                                                        * (2. - FlowEl(iloc)%b * (2. - FlowEl(iloc)%muPlastic / FlowEl(iloc)%mueff))
                
            ELSE        ! transition from laminar to turbulent
                FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &
                                                        * (2. + (2. - FlowEl(iloc)%muPlastic / FlowEl(iloc)%mueff) &
                                                            * ((FlowEl(iloc)%a * FlowEl(iloc)%GenRe**(1. - FlowEl(iloc)%b) - 16.) / 2000. / FlowEl(iloc)%f - 1.))
            END IF
        END IF
        END IF
        
        
        IF (FlowEl(iloc)%FricToQPartialDiff < 0.0) THEN
!WRITE (*,*) ' iloc, Re, FricPressLoss, FricToQPartialDiff' , iloc, FlowEl(iloc)%GenRe, FlowEl(iloc)%FricPressLoss, FlowEl(iloc)%FricToQPartialDiff
!CALL ERRORSTOP('Error in Calculating FricToQPartialDiff')
        END IF
    
END SUBROUTINE PartialDerivativeFricToFlowRate