mahmood19227
/
SimulationCore2
1
0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
Simulation Core
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
69 Commits
2 Branches
134 MiB
 
 
 
 
 
 
Tree: a57ac8edd6
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'a57ac8edd6'
${ noResults }
SimulationCore2/FluidFlow/Frictional_Press_Drop_Calc.i90

232 lines
12 KiB
Raw Blame History 

  1. # 1 "/mnt/c/Projects/VSIM/SimulationCore2/FluidFlow/Frictional_Press_Drop_Calc.f90"

  2. SUBROUTINE FricPressDrop(iloc)

  3. 

  4.     

  5.     

  6. !!      Record of revisions

  7. !!          Date            Programmer          Discription of change

  8. !!         ------          ------------        -----------------------

  9. !!        1396/07/23        Sheikh              Original code

  10. !!

  11.     

  12.     

  13.     

  14.     USE FricPressDropVarsModule

  15.     USE SimulationVariables

  16.     USE Fluid_Flow_Startup_Vars

  17.     USE CError

  18.     

  19.     IMPLICIT NONE    

  20. 

  21.     INTEGER :: iloc

  22.     REAL :: TauZero

  23.     TauZero = 12.0

  24. !ActiveRheologyModel = Herschel_Bulkley_RheologyModel

  25. 

  26. ! 0 = Power Law   ,   1 = Bingham Plastic   ,   2 = Newtonian

  27. !TotFricPressLoss = 0.0

  28.     

  29.         FlowEl(iloc)%alpha = 1                ! assume that all elements have annulus geometry

  30.         FlowEl(iloc)%dPdLfric = 0.0

  31.         FlowEl(iloc)%f = 0.0

  32.         FlowEl(iloc)%FlowRate = ABS(FlowEl(iloc)%FlowRate)

  33.         

  34.         

  35.     

  36.     IF ((FlowEl(iloc)%FlowRate >= PressFlowrateTolerance) &

  37.         .AND. (FlowEl(iloc)%MaterialType /= 1) &        ! not gas kick

  38.         .AND. (ABS(FlowEl(iloc)%Length) >= PressLengthTolerance) &

  39.         .AND. (FlowEl(iloc)%MaterialType /= 4)) THEN    ! not air

  40.         

  41.         IF (FlowEl(iloc)%Id==0) THEN

  42.             FlowEl(iloc)%alpha = 0

  43.         END IF

  44.         

  45.         FlowEl(iloc)%muPlastic = FlowEl(iloc)%Theta600 - FlowEl(iloc)%Theta300   ! cp

  46.         FlowEl(iloc)%YieldP = 2.0 * FlowEl(iloc)%Theta300 - FlowEl(iloc)%Theta600   ! lbf/100ft**2

  47.         FlowEl(iloc)%nIndex = 3.32 * log10(FlowEl(iloc)%Theta600 / FlowEl(iloc)%Theta300)

  48.         FlowEl(iloc)%kIndex = 510.0 * FlowEl(iloc)%Theta300 / (511.0**FlowEl(iloc)%nIndex)    ! rabete fv2

  49.         IF (data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel .AND. FlowEl(iloc)%alpha == 0) THEN

  50.             FlowEl(iloc)%kIndex = 1.066 * FlowEl(iloc)%Theta300 / (511.0**FlowEl(iloc)%nIndex)

  51.         ELSE IF (data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel .AND. FlowEl(iloc)%alpha == 1) THEN

  52.             FlowEl(iloc)%nIndex = 3.32 * log10((FlowEl(iloc)%Theta600 - TauZero) / (FlowEl(iloc)%Theta300 - TauZero))

  53.             FlowEl(iloc)%kIndex = 1.066 * (FlowEl(iloc)%Theta300 - TauZero) / (511.0**FlowEl(iloc)%nIndex)

  54.         END IF

  55.         

  56.         

  57.         

  58.         

  59. ! Calculating velocity

  60.         FlowEl(iloc)%vel = 0.408 * FlowEl(iloc)%FlowRate / (FlowEl(iloc)%Od**2 - FlowEl(iloc)%Id**2)   ! velocity in ft/s

  61. !FlowEl(iloc)%vel = 24.51 * FlowEl(iloc)%FlowRate / (FlowEl(iloc)%Od**2 - FlowEl(iloc)%Id**2)   ! velocity in ft/min

  62. 

  63. !IF (FlowModel == Bingham_RheologyModel) THEN   ! Bingham Plastic

  64. !    FlowEl(iloc)%Gf = (2. + FlowEl(iloc)%alpha) / 2.

  65. !ELSE IF (FlowModel == PowerLow_RheologyModel) THEN

  66. !    FlowEl(iloc)%Gf = ((3. - FlowEl(iloc)%alpha) * FlowEl(iloc)%nIndex + 1.) / FlowEl(iloc)%nIndex / (4. - FlowEl(iloc)%alpha) * (2. + FlowEl(iloc)%alpha) / 2.

  67. !END IF

  68.         

  69. !FlowEl(iloc)%gammaW = 1.6 * FlowEl(iloc)%Gf * FlowEl(iloc)%vel / FlowEl(iloc)%Dhyd

  70.         

  71. !IF (FlowModel == Bingham_RheologyModel) THEN   ! Bingham Plastic

  72. !    FlowEl(iloc)%tauW = 1.067 * ((4. - FlowEl(iloc)%alpha) / (3. - FlowEl(iloc)%alpha) * FlowEl(iloc)%YieldP + FlowEl(iloc)%muPlastic * FlowEl(iloc)%gammaW)

  73. !    !FlowEl(iloc)%tauW = 1.067*(FlowEl(iloc)%YieldP+FlowEl(iloc)%muPlastic*FlowEl(iloc)%gammaW)

  74. !ELSE IF (FlowModel == PowerLow_RheologyModel) THEN     ! Power law

  75. !    FlowEl(iloc)%tauW = 1.067 * FlowEl(iloc)%kIndex * FlowEl(iloc)%gammaW**FlowEl(iloc)%nIndex

  76. !END IF

  77.         

  78. ! Calculating effective or apparent viscosity

  79.         IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic

  80.             FlowEl(iloc)%mueff = FlowEl(iloc)%muPlastic + 5. * FlowEl(iloc)%YieldP * FlowEl(iloc)%Dhyd / FlowEl(iloc)%vel

  81. !write(*,*) 'pointer1' , FlowEl(iloc)%muPlastic ,  FlowEl(iloc)%YieldP , FlowEl(iloc)%Dhyd , FlowEl(iloc)%vel

  82.         ELSE IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel .OR. data%Configuration%Mud%ActiveRheologyModel == Herschel_Bulkley_RheologyModel) THEN   ! Power Law

  83.             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

  84.             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

  85. !write(*,*) 'pointer2' , FlowEl(iloc)%kIndex ,FlowEl(iloc)%alpha , FlowEl(iloc)%vel ,FlowEl(iloc)%Dhyd,FlowEl(iloc)%nIndex ,FlowEl(iloc)%Gf ,FlowEl(iloc)%nIndex

  86.         END IF

  87.         

  88.         FlowEl(iloc)%gammaW = 96.0 * FlowEl(iloc)%Gf * FlowEl(iloc)%vel / FlowEl(iloc)%Dhyd

  89.         

  90.         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

  91.         

  92.         

  93. ! Calculating Reynolds number

  94.         IF (FlowEl(iloc)%Od == FlowEl(iloc)%Dhyd) THEN

  95.             FlowEl(iloc)%GenRe = 928. * FlowEl(iloc)%density * FlowEl(iloc)%vel * FlowEl(iloc)%Dhyd / FlowEl(iloc)%mueff

  96.         ELSE

  97.             FlowEl(iloc)%GenRe = 757. * FlowEl(iloc)%density * FlowEl(iloc)%vel * FlowEl(iloc)%Dhyd / FlowEl(iloc)%mueff

  98.         END IF

  99.         

  100. !FlowEl(iloc)%GenRe = 2997 * FlowEl(iloc)%density * FlowEl(iloc)%vel**2 / 19.36 / FlowEl(iloc)%tauW

  101.         

  102. ! Calculating friction factor

  103.         IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic

  104.             IF (FlowEl(iloc)%GenRe <= 2000.0) THEN   ! laminar regime

  105.                 FlowEl(iloc)%f = 16.0 / FlowEl(iloc)%GenRe

  106.             ELSE IF (FlowEl(iloc)%GenRe >= 4000.0) THEN   ! turbulent regime

  107.                 FlowEl(iloc)%a = 0.0791

  108.                 FlowEl(iloc)%b = 0.25

  109.                 FlowEl(iloc)%f = FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b

  110.             ELSE    !! transition from laminar to turbulent regime

  111.                 FlowEl(iloc)%a = 0.0791

  112.                 FlowEl(iloc)%b = 0.25

  113.                 FlowEl(iloc)%f = (4000.0 - FlowEl(iloc)%GenRe) / 2000.0 * 16. / FlowEl(iloc)%GenRe &

  114.                                 + (FlowEl(iloc)%GenRe - 2000.0) / 2000.0 * FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b

  115.             END IF

  116.             

  117.         ELSE IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel) THEN     ! Power law

  118.             FlowEl(iloc)%ReCritLam = 3470. - 1370. * FlowEl(iloc)%nIndex

  119.             FlowEl(iloc)%ReCritTurb = 4270. - 1370. * FlowEl(iloc)%nIndex

  120.             

  121.             IF (FlowEl(iloc)%GenRe <= FlowEl(iloc)%ReCritLam) THEN   ! laminar regime

  122.                 FlowEl(iloc)%f = 16.0 / FlowEl(iloc)%GenRe / (1 - 0.184 * FlowEl(iloc)%alpha)

  123.             ELSE IF (FlowEl(iloc)%GenRe >= FlowEl(iloc)%ReCritTurb) THEN   ! turbulent regime

  124.                 FlowEl(iloc)%a = (log10(FlowEl(iloc)%nIndex) + 3.93) / 50.

  125.                 FlowEl(iloc)%b = (1.75 - log10(FlowEl(iloc)%nIndex)) / 7.

  126.                 FlowEl(iloc)%f = FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b

  127.             ELSE

  128.                 FlowEl(iloc)%a = (log10(FlowEl(iloc)%nIndex) + 3.93) / 50.

  129.                 FlowEl(iloc)%b = (1.75 - log10(FlowEl(iloc)%nIndex)) / 7.

  130.                 FlowEl(iloc)%f = (FlowEl(iloc)%ReCritTurb - FlowEl(iloc)%GenRe) / 800.0 * 16. / FlowEl(iloc)%GenRe &

  131.                                 + (FlowEl(iloc)%GenRe - FlowEl(iloc)%ReCritLam) / 800.0 * FlowEl(iloc)%a / FlowEl(iloc)%GenRe**FlowEl(iloc)%b

  132.             END IF

  133.             

  134.         END IF

  135.         

  136. !WRITE (*,*) 'fric press drop', iloc

  137. !WRITE (*,*) 'Length', ABS(REAL(FlowEl(iloc)%Length))

  138. !WRITE (*,*) 'FlowRate', FlowEl(iloc)%FlowRate

  139. !WRITE (*,*) 'Theta600 , Theta300', FlowEl(iloc)%Theta600 , FlowEl(iloc)%Theta300

  140. !WRITE (*,*) 'Dhyd', FlowEl(iloc)%Dhyd

  141. !WRITE (*,*) 'GenRe', FlowEl(iloc)%GenRe

  142. !WRITE (*,*) 'f', FlowEl(iloc)%f

  143.         

  144.         

  145.         

  146.         

  147.     END IF

  148. 

  149. ! Frictional pressure loss gradient calculation

  150. ! FlowEl(iloc)%dPdLfric = 1.076 * FlowEl(iloc)%f * FlowEl(iloc)%vel**2 * FlowEl(iloc)%density / 10**5 / FlowEl(iloc)%Dhyd

  151.         FlowEl(iloc)%dPdLfric = FlowEl(iloc)%f * (FlowEl(iloc)%vel)**2 * FlowEl(iloc)%density / 25.81 / FlowEl(iloc)%Dhyd

  152.         FlowEl(iloc)%FricPressLoss = FlowEl(iloc)%dPdLfric * ABS(REAL(FlowEl(iloc)%Length))

  153.         

  154.         IF (FlowEl(iloc)%FrictionDirection == -1) THEN

  155.             FlowEl(iloc)%FlowRate = - FlowEl(iloc)%FlowRate

  156.             FlowEl(iloc)%dPdLfric = - FlowEl(iloc)%dPdLfric

  157.             FlowEl(iloc)%FricPressLoss = - FlowEl(iloc)%FricPressLoss

  158.         END IF

  159.         

  160.         

  161.             

  162. !END DO

  163. 

  164. 

  165. END SUBROUTINE FricPressDrop

  166.     

  167.     

  168. SUBROUTINE PartialDerivativeFricToFlowRate(iloc)

  169.     

  170.     USE FricPressDropVarsModule

  171.     USE SimulationVariables

  172.     USE Fluid_Flow_Startup_Vars

  173.     use KickVARIABLESModule

  174.     USE CError

  175.     

  176.     

  177.     IMPLICIT NONE    

  178. 

  179.     INTEGER :: iloc

  180.     

  181.     FlowEl(iloc)%FricToQPartialDiff = 0.0

  182. !FlowEl(iloc)%FlowRate = ABS(FlowEl(iloc)%FlowRate)

  183.     

  184.     

  185.     IF ((ABS(FlowEl(iloc)%FlowRate) >= PressFlowrateTolerance) &

  186.         .AND. (FlowEl(iloc)%MaterialType /= 1) &        ! not gas kick

  187.         .AND. (ABS(FlowEl(iloc)%Length) >= PressLengthTolerance) &

  188.         .AND. (FlowEl(iloc)%MaterialType /= 4)) THEN    ! not air

  189.         

  190.         IF (data%Configuration%Mud%ActiveRheologyModel == PowerLaw_RheologyModel) THEN     ! Power law

  191. !IF (FlowEl(iloc)%Flowrate == 0.0) THEN

  192. !    FlowEl(iloc)%Flowrate = 10.0

  193. !    CALL FricPressDrop(iloc)

  194. !END IF

  195.             IF (FlowEl(iloc)%GenRe <= FlowEl(iloc)%ReCritLam) THEN     ! laminar flow

  196.                 FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate * FlowEl(iloc)%nIndex

  197.                 

  198.             ELSE IF (FlowEl(iloc)%GenRe >= FlowEl(iloc)%ReCritTurb) THEN       ! turbulent flow

  199.                 FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &

  200.                                                         * (2. - FlowEl(iloc)%b * (2. - FlowEl(iloc)%nIndex))

  201.             

  202.             ELSE        ! transition from laminar to turbulent

  203.                 FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &

  204.                                                         * (2. + (2. - FlowEl(iloc)%nIndex) &

  205.                                                             * ((FlowEl(iloc)%a * FlowEl(iloc)%GenRe**(1. - FlowEl(iloc)%b) - 16.) / 800. / FlowEl(iloc)%f - 1.))

  206.             END IF

  207.         

  208.         ELSE IF (data%Configuration%Mud%ActiveRheologyModel == Bingham_RheologyModel) THEN   ! Bingham Plastic

  209.             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

  210.                 FlowEl(iloc)%FricToQPartialDiff = (16. * FlowEl(iloc)%muPlastic * REAL(FlowEl(iloc)%Length) * 2.224 * (10.)**(-3)) &

  211.                                                         / (25.81 * 928. * (1 - 0.184 * FlowEl(iloc)%alpha) * FlowEl(iloc)%Dhyd**2 * FlowEl(iloc)%Area)

  212.         

  213.             ELSE IF (FlowEl(iloc)%GenRe >= 4000.0) THEN       ! turbulent flow

  214.                 FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &

  215.                                                         * (2. - FlowEl(iloc)%b * (2. - FlowEl(iloc)%muPlastic / FlowEl(iloc)%mueff))

  216.                 

  217.             ELSE        ! transition from laminar to turbulent

  218.                 FlowEl(iloc)%FricToQPartialDiff = FlowEl(iloc)%FricPressLoss / FlowEl(iloc)%FlowRate &

  219.                                                         * (2. + (2. - FlowEl(iloc)%muPlastic / FlowEl(iloc)%mueff) &

  220.                                                             * ((FlowEl(iloc)%a * FlowEl(iloc)%GenRe**(1. - FlowEl(iloc)%b) - 16.) / 2000. / FlowEl(iloc)%f - 1.))

  221.             END IF

  222.         END IF

  223.         END IF

  224.         

  225.         

  226.         IF (FlowEl(iloc)%FricToQPartialDiff < 0.0) THEN

  227. !WRITE (*,*) ' iloc, Re, FricPressLoss, FricToQPartialDiff' , iloc, FlowEl(iloc)%GenRe, FlowEl(iloc)%FricPressLoss, FlowEl(iloc)%FricToQPartialDiff

  228. !CALL ERRORSTOP('Error in Calculating FricToQPartialDiff')

  229.         END IF

  230.     

  231. END SUBROUTINE PartialDerivativeFricToFlowRate