# 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