Computation of the steady state






    Before doing any calculation of flood, it is necessary to have the steady-state results, which corresponds to the river in everyday conditions. During these calculations, we realized that the first boundary conditions were not  fitting our problem.

            *    First of all, the flow of water came from the whole plain border instead of the minor bed of the river. So the inlet flow rate was fixed to 10 m3/s for only the minor bed of the river.
        That implies that the first meters of the computation field do not represent anything. Indeed, the boundary condition for the major bed will be dry for the whole simulation, which is wrong during the flood. The solution would be to extend the computation field beyond the inlet boundary. In that way, it would be possible to model the fact that the inlet section of the river changes with the flow rate.

            *    Moreover, the height of water given for the outlet conditions had to be calculated from the Strickler law:
 

                                                              with:       K:    Stricler coefficient
                                                                                                                           R:    equivalent radius
                                                                                                                                                      I:    constant = 1e-3
                R is defined for the minor bed by the quotient of the area over the wet perimeter: 
                K represents the rubbing and  is equal to 30 for the river and 20 for the plain. But K will be taken constant for both the river and the plain and equal to 30.

NB: This equation is wrong but we realized it later ....

This calculation of H gave us different values because some terms in the equation were so little that the resolution of this equation needs a good precision of the results. The most realistic value was found with EXCEL Software and was 2.51 m. For the first calculations, the given height of water was 2m. In this condition, the outlet flow rate was much higher than 10 m3/s. That shows that 2 m was too low. That is why we decided to run some calculations with a fixed height of water at the outlet, equal to 2.5m.
The steady state for the free surface profile is displayed by the following picture.

    The two blue lines correspond to results at 700 and 1300 s. Since they are quite similar, we considered that the steady state was reached.

    Since the boundary condition "Given Height Of Water" cannot be applied for the flood computation, we decided to change the outlet boundary condition. The flow rate was given by the Strickler law. The inlet flow rate was still equal to 10 m3/s.


 
 

We can notice that the height of water is about higher than expected. That is due to a mistake in the computation of the Strickler law. Indeed, the flow rate we calculated with the Strickler law was actually Q/S. So we over-estimated H.

                *    The initial condition for all the previous computations was "Constant altitude of water".  We have chosen this one because the CONDIN subroutine was wrong and did not consider the initial condition "Constant height of water".  Indeed, here is the warning we received when we put "Constant height":

    First we thought that it would have not had a huge consequence on the results of our computation. But after all the problems of convergence we have met, we tried to change this condition in order to facilitate the convergence of the software. We changed the subroutine CONDIN but we did not have the time to compute some realistic results.
 
 

How to solve these problems.

  • The main error comes from the Q(H) formula . the right formula sould be multiply by S hence :
  • Q=S*K*I^0.5*S*R^(2/3)
  • The problems of initialization can be resolved understanding the CONLIM fortran file. It is there that H(x,y) can be specified. We specify a 1.2 m  constant height in the river bed.

  • Here is the source file :

    Doing like that we reach the steady state very quicly and the solution semms more realistic. e.g. the in the medium plan y-ward we have:


     
     
     
     

    C-----------------------------------------------------------------------
    C
          STOP
          END
    C                       *****************
                            SUBROUTINE CONDIN
    C                       *****************
    C
         *(U,V,H,T,AK,EP,VISC,AUX3,X,Y,ZF,NPOIN,TRAC,TRAC0,ITURB,PROPNU,
         * TEMPS,AKEP,CDTINI,COTINI,HAUTIN,NBOR,NPTFR,PRIVE)
    C
    C***********************************************************************
    C PROGICIEL : TELEMAC 2D        31/10/91  J-M HERVOUET TEL: 30 87 80 18
    C
    C***********************************************************************
    C
    C     FONCTION  : INITIALISATION DES TABLEAUX DES GRANDEURS PHYSIQUES
    C
    C-----------------------------------------------------------------------
    C                             ARGUMENTS
    C .________________.____.______________________________________________
    C |      NOM       |MODE|                   ROLE
    C |________________|____|______________________________________________
    C |    U , V       |<-- | COORDONNEES DES VECTEURS VITESSE.
    C |    H           |<-- | HAUTEUR D'EAU VARIABLE
    C |    T           |<-- | TRACEUR.
    C |    AK          |<-- | ENERGIE TURBULENTE.
    C |    EP          |<-- | DISSIPATION.
    C |    VISC        |<-- | VISCOSITE TURBULENTE
    C |    AUX3        |<-- | VARIABLE SUPPLEMENTAIRE PREVUE.
    C |    X,Y         | -->| COORDONNEES DES POINTS DU MAILLAGE .
    C |    ZF          | -->| COTE DE LA SURFACE LIBRE .
    C |    LIUBOR      | -->| INDICATEURS DE CONDITIONS AUX LIMITES SUR U.
    C |    LIVBOR      | -->| INDICATEURS DE CONDITIONS AUX LIMITES SUR V.
    C |    LIHBOR      | -->| INDICATEURS DE CONDITIONS AUX LIMITES SUR H.
    C |    IBOR,JBOR   | -->| INDICES DES POINTS FRONTIERES.
    C |    INDIC       | -->| INDICATEURS DE SITUATION DES POINTS
    C |    UBOR        | -->| VALEURS DE U IMPOSEES AUX BORDS.
    C |    VBOR        | -->| VALEURS DE V IMPOSEES AUX BORDS.
    C |    HBOR        | -->| VALEURS DE H IMPOSEES AUX BORDS.
    C |________________|____|______________________________________________
    C MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
    C***********************************************************************
    C
          IMPLICIT NONE
          INTEGER LNG,LU
          COMMON/INFO/LNG,LU
    C
     INTEGER IPOIN,NPOIN,ITURB,NPTFR,NBOR(NPTFR)
          DOUBLE PRECISION U(NPOIN)  , V(NPOIN) , H(NPOIN)   ,T(NPOIN)
          DOUBLE PRECISION AK(NPOIN) , EP(NPOIN), VISC(NPOIN),AUX3(NPOIN)
          DOUBLE PRECISION X(NPOIN)  , Y(NPOIN) , ZF(NPOIN)  ,PRIVE(NPOIN,*)
    C
          DOUBLE PRECISION PROPNU,TEMPS,Z,C,COTINI,HAUTIN,COTE,TRAC0,EIKON
    C
     
    C
          LOGICAL TRAC,AKEP
    C
          CHARACTER *(*) CDTINI
    C
          INTRINSIC EXP
    C
    C-----------------------------------------------------------------------
    C
    C   INITIALISATION DU TEMPS
    C
          TEMPS = 0.D0
    C
    C-----------------------------------------------------------------------
    C
    C   INITIALISATION DES VITESSES : VITESSES NULLES
    C
          CALL OV( 'X=C     ' , U , Y , Z , 0.D0 , NPOIN )
          CALL OV( 'X=C     ' , V , Y , Z , 0.D0 , NPOIN )
    C
    C-----------------------------------------------------------------------
    C
    C   INITIALISATION DE H , LA HAUTEUR D'EAU
    C
          IF(CDTINI(1:10).EQ.'COTE NULLE') THEN
            COTE = 0.D0
            CALL OV( 'X=C     ' , H , Y  , Z , COTE , NPOIN )
            CALL OV( 'X=X-Y   ' , H , ZF , Z , C    , NPOIN )
          ELSEIF(CDTINI(1:14).EQ.'COTE CONSTANTE') THEN
            COTE = COTINI
            CALL OV( 'X=C     ' , H , Y  , Z , COTE , NPOIN )
            CALL OV( 'X=X-Y   ' , H , ZF , Z , C    , NPOIN )
          ELSEIF(CDTINI(1:13).EQ.'HAUTEUR NULLE') THEN
            CALL OV( 'X=C     ' , H , Y  , Z , 0.D0 , NPOIN )
          ELSEIF(CDTINI(1:13).EQ.'PARTICULIERES') THEN
          DO 10 IPOIN=1,NPOIN
            IF((X(IPOIN).GE.-5).and.(X(IPOIN).LE.5)) THEN
     
            H(IPOIN) = 1.2D0
     
            ELSE
            H(IPOIN) = 0.D0
     
            ENDIF
     
    10    CONTINUE
          ELSE
            PRINT *,'CONDIN : CONDITION INITIALE NON PREVUE : ',CDTINI
            CALL PLANTE(0)
            STOP
          ENDIF
    C
    C-----------------------------------------------------------------------
    C
    C   INITIALISATION DU TRACEUR
    C
          IF(TRAC) CALL OV( 'X=C     ' , T , Y , Z , TRAC0 , NPOIN )
    C
    C-----------------------------------------------------------------------
    C
    C INITIALISATION DE LA VISCOSITE
    C
          CALL OV( 'X=C     ' , VISC , Y , Z , PROPNU , NPOIN )
    C
    C-----------------------------------------------------------------------
    C
          RETURN
          END