friction_lindner.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\friction_lindner.f
00002 !
00060                      SUBROUTINE FRICTION_LINDNER
00061 !                    ***************************
00062 !
00063      &(VA,HA,CF,VK,G,DP,SP,CP)
00064 !
00065 !***********************************************************************
00066 ! TELEMAC2D   V6P1                                   21/08/2010
00067 !***********************************************************************
00068 !
00069 !
00070 !
00071 !
00072 !
00073 !
00074 !
00075 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00076 !| CF             |-->| FRICTION COEFFICIENT FOR BOTTOM ROUGHNESS
00077 !| CP             |<--| FRICTION COEFFICIENT FOR NON-SUBMERGED VEGETATION
00078 !| DP             |-->| DIAMETER OF ROUGHNESS ELEMENT
00079 !| G              |-->| GRAVITY ACCELERATION
00080 !| HA             |-->| FLOW DEPTH
00081 !| SP             |-->| SPACING OF ROUGHNESS ELEMENT
00082 !| VA             |-->| VELOCITY
00083 !| VK             |-->| KINEMTIC VISCOSITY
00084 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00085 !
00086       IMPLICIT NONE
00087 !
00088 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00089 !
00090       DOUBLE PRECISION, INTENT(IN)  :: VA,HA,CF,VK,G,DP,SP
00091       DOUBLE PRECISION, INTENT(OUT) :: CP
00092 !
00093 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00094 !
00095       INTEGER,          PARAMETER :: KMAXITER = 200
00096       DOUBLE PRECISION, PARAMETER :: KPRECISION = 1.0D-3
00097       INTEGER                     :: CWRAV
00098       INTEGER                     :: CWRMAX
00099       INTEGER                     :: CWRCOUNT
00100       INTEGER                     :: ANLAV
00101       INTEGER                     :: ANLMAX
00102       INTEGER                     :: ANLCOUNT
00103       INTEGER                     :: ITERR
00104 !
00105       INTEGER :: I, J
00106       INTEGER :: ICWR               ! ITERATION COUNTER: CWR
00107       INTEGER :: IANL               ! ITERATION COUNTER: ANL
00108       INTEGER :: REALROOTS
00109       LOGICAL :: LCWR
00110 !
00111       DOUBLE PRECISION :: CW, CWR, RCWR, DCWR, ANL, ANB, RANL
00112       DOUBLE PRECISION :: CWR1, CWR2, DCWR1, DCWR2
00113       DOUBLE PRECISION :: LAMBDA, FR
00114       DOUBLE PRECISION :: X(3), VRATIO, HRATIO
00115       DOUBLE PRECISION :: ALFA, ACOF, BCOF, CCOF, DCOF
00116 !
00117       DOUBLE PRECISION :: TMP1, TMP2
00118 !
00119 !=======================================================================!
00120 !=======================================================================!
00121 !                               PROGRAMME                               !
00122 !=======================================================================!
00123 !=======================================================================!
00124 !
00125       LCWR     = .TRUE.
00126       CWRAV    = 0
00127       CWRMAX   = 0
00128       CWRCOUNT = 0
00129       ANLAV    = 0
00130       ANLMAX   = 0
00131       ANLCOUNT = 0
00132       ITERR    = 0
00133       IF ((DP < 1.0E-3)     .OR.(SP < 1.0E-2).OR.
00134      &    (ABS(VA) < 1.0E-3).OR.(HA < 1.0E-3)     ) THEN
00135         CP = 0.D0
00136       ELSE
00137         ! INITIALIZATION
00138         ! --------------
00139         CWR   = 1.0      ! DRAG COEFFICIENT
00140         ANL   = SP/2.0   ! WAKE LENGTH OF A CYLINDER
00141         CWR1  = 1.0
00142         CWR2  = 1.0
00143         DCWR1 = 0.0
00144         DCWR2 = 0.0
00145         ! START OF ITERATION FOR CWR
00146         ! --------------------------
00147         DO ICWR = 1, KMAXITER
00148           ! SUPERPOSED FRICTION COEFFICIENT
00149           ! -------------------------------
00150           LAMBDA = 8.0D0*CF  +  4.0D0*CWR*HA*DP/SP/SP
00151           ! DRAG COEFFICIENT CW FOR ONE CYLINDER
00152           ! ------------------------------------
00153           CALL DRAGCOEFF(VA, DP, VK, CW)
00154           ! WAKE LENGTH OF A CYLINDER (ITERATIVE COMPUTATION)
00155           ! -------------------------------------------------
00156           DO J=1, KMAXITER
00157             TMP1 = 1.0D0  +  ANL*LAMBDA/4.0D0/HA
00158             TMP2 = 30.0D0/ABS(TMP1)**(1.5)
00159             RANL = CW*DP*ABS(TMP2)**(1.429)
00160             ! TEST FOR CONVERGENCE
00161             ! --------------------
00162             IF (ABS((RANL-ANL)/RANL) < KPRECISION) THEN
00163               ANL  = RANL
00164               IANL = -1*J
00165               EXIT
00166             ENDIF
00167             ANL = 0.5 * (RANL + ANL)
00168           ENDDO
00169           ! STATISTICS OF CWR ITERATION
00170           ! ---------------------------
00171           IF ( IANL > 0 ) THEN
00172             ANL = SP/2.0D0
00173           ELSE
00174             IANL = ABS(IANL)
00175             ANLCOUNT = ANLCOUNT + 1
00176             ANLAV = IANL + ANLAV
00177             IF (IANL > ANLMAX) ANLMAX = IANL
00178           ENDIF
00179           ! WAKE WIDTH
00180           ! ----------
00181           ANB = 0.24 * ABS(ANL)**(0.59) * ABS(CW*DP)**(0.41)
00182           ! RATIO OF VELOCITY IN FRONT OF AND BEHIND CYLINDER
00183           ! -------------------------------------------------
00184           VRATIO = 1.151 * ABS(ANL/SP)**(-0.483)
00185      &           +   0.5 * ABS(ANB/SP)**(1.1)
00186           ! RATIO OF FLOW DEPTH
00187           ! -------------------
00188           FR = VA / SQRT( G * HA ) ! FROUDE NUMBER
00189           ALFA = DP / SP
00190           ACOF =  FR * FR * (1.0D0 - ALFA * CWR/2.0D0)
00191           BCOF = -FR * FR - (1.0D0 - ALFA) / 2.0D0
00192           CCOF =  0.0D0
00193           DCOF = (1.0D0 - ALFA) / 2.0D0
00194           HRATIO = 1.0D0
00195           IF (ABS(ACOF) < 1.0E-10) THEN
00196             HRATIO = SQRT( -DCOF / BCOF)
00197           ELSE
00198             CALL CUBEEQUATION(ACOF, BCOF, CCOF, DCOF, REALROOTS, X)
00199             DO I = 1, REALROOTS
00200               IF (X(I) > 0.0  .AND.  X(I) < 1.0)  THEN
00201                 HRATIO = X(I)
00202                 EXIT
00203               ENDIF
00204             ENDDO
00205           ENDIF
00206           ! REVISE DRAG COEFFICIENT CWR
00207           ! ---------------------------
00208           RCWR = 1.3124D0*CW*VRATIO + 2.0D0*(1.0D0-HRATIO)/FR/FR
00209           ! TEST FOR CONVERGENCE
00210           ! --------------------
00211           IF ( ABS((RCWR-CWR)/RCWR) < KPRECISION ) THEN
00212             LCWR = .FALSE.
00213 !           ICWR = -1/ICWR
00214             EXIT
00215           ENDIF
00216           ! USE PEGASUS ALGORITHM FOR CWR ITERATION
00217           ! ---------------------------------------
00218           DCWR = RCWR - CWR
00219           IF ((ICWR >= 3) .AND. (DCWR1*DCWR2 < 0.0D0)) THEN
00220             IF (DCWR2*DCWR < 0.0D0) THEN
00221               DCWR1 = DCWR2/(DCWR2+DCWR)*DCWR1
00222             ELSE
00223               CWR1  = CWR2
00224               DCWR1 = DCWR2
00225             ENDIF
00226             CWR2  = CWR
00227             DCWR2 = DCWR
00228             CWR   = CWR2 - DCWR2*(CWR2-CWR1)/(DCWR2-DCWR1)
00229           ELSE
00230             CWR1 = CWR2
00231             DCWR1 = DCWR2
00232             CWR2 = CWR
00233             DCWR2 = DCWR
00234             IF ((ICWR >= 2) .AND. (DCWR1*DCWR2 < 0.0 )) THEN
00235               CWR = CWR2 - DCWR2*(CWR2-CWR1)/(DCWR2-DCWR1)
00236             ELSE
00237               CWR = RCWR
00238             ENDIF
00239           ENDIF
00240         ENDDO !ICWR = 1, KMAXITER
00241 !
00242         IF (LCWR) THEN
00243           ITERR = ITERR + 1
00244           CP = -1.D0
00245         ELSE
00246           ! STATISTICS OF CWR ITERATION
00247           ! ---------------------------
00248           ICWR = -1/ICWR ! AS THE PROGRAM RISMO2D FROM THE BAW
00249           ICWR = -ICWR
00250           CWRCOUNT = CWRCOUNT + 1
00251           CWRAV = ICWR + CWRAV
00252           IF (ICWR > CWRMAX) CWRMAX = ICWR
00253           CP = LAMBDA/8.0D0 - CF
00254         ENDIF
00255       ENDIF
00256 !
00257 !=======================================================================!
00258 !=======================================================================!
00259 !
00260       RETURN
00261       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0