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 FricPressDropVars
    USE CMudPropertiesVariables
    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 (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 (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 (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 (ActiveRheologyModel == PowerLaw_RheologyModel .OR. 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 (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 (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 FricPressDropVars
    USE CMudPropertiesVariables
    USE Fluid_Flow_Startup_Vars
    USE KickVariables
    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 (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 (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