Télécharger la plaquette Cast3M


fr en

Télécharger hbmhill.eso

Retour à la liste

Numérotation des lignes : 

hbmhill

  1. C HBMHILL   SOURCE    OF166741  26/05/11    21:15:11     12538          
  2.  
  3. *=======================================================================
  4. *                          Calcul de stabilite
  5. *=======================================================================
  6. * On utilise la methode de Hill (domaine frequentiel) pour une solution
  7. * periodique (Q1,OMEG) des equations d'equilibre dynamique R(X,w) = 0.
  8. *
  9. *                      (Rx + µ*D1 + µ**2*D2)*Phi = 0
  10. * µ : valeurs propres
  11. * Phi: vecteurs propres
  12. * ==> 2*n*(2*H+1) paires (µ,Phi) pour H harmoniques et n ddl, dont 2*n
  13. *     correspondent aux exposants de Floquet.
  14. * Dans cette sous-routine, on utilise une procedure simplifiee pour
  15. * la caracterisation des bifurcations.
  16. * -----------
  17. *   Entrees
  18. * -----------
  19. *         OMEG       : frequence du cycle
  20. *     NT,NHBM,NDDL   : taille du systeme, # d'harmoniques, # de DDL
  21. *        XM,XASM     : matrices de masse, amortissement
  22. *          RX        : matrice Jacobienne
  23. *       dv0,dv1      : variations dans le parametre de continuation
  24. *                      aux pas precedent et actuel, respectivement
  25. *        CPOSRE      : indicateur de bifurcation au pas precedent
  26. *        ZSTAB       : stabilite du pas precedent (0 si initialisation)
  27. * -----------
  28. *   Sorties
  29. * -----------
  30. *        CPOSRE      : indicateurs de bifurcation au pas actuel
  31. *         FLAG       : caracterise le type de bifurcation
  32. *       ZSTABN=ZSTAB : stabilite de la solution (Q1,OMEG)
  33. *=======================================================================
  34.  
  35.       SUBROUTINE HBMHILL(NT,NHBM,NDDL,OMEG,Q1,RX,XM,XASM,dv0,dv1,ZINIT,
  36.      &                  LHILL,ZSTABN,CPOSRE,FLAG,ZSTAB)
  37.  
  38.       IMPLICIT INTEGER(I-N)
  39.       IMPLICIT REAL*8(A-H,O-Z)
  40.  
  41. -INC CCREEL
  42.  
  43.       INTEGER NT,NHBM,NDDL,CPOSRE
  44.       REAL*8 OMEG,dv0,dv1
  45.       REAL*8 Q1(*), RX(NT,*),XM(NDDL,*),XASM(NDDL,*),LHILL(2,*)
  46.       LOGICAL ZSTABN,ZSTAB,ZINIT
  47.       CHARACTER*8 FLAG
  48.  
  49.       SEGMENT mwrkml
  50.         INTEGER INDEIG(2*NT)
  51.         REAL*8 Rw(NT),MATJ(NT+1,NT+1),t0(NT+1)
  52.         REAL*8 EXPIM(2*NT),EXPRE(2*NT)
  53.         REAL*8 VR(2*NT,2*NT),VL(2*NT,2*NT),WORK(8*NT),BB
  54.         REAL*8 DEL1(NT,NT),DEL2(NT),Mi,Ci,JJ(2*NT,2*NT)
  55.         REAL*8 mumx(2*NDDL),ZIINDM
  56. cbp       REAL*8 ZR(2*NDDL),ZI(2*NDDL)
  57.       ENDSEGMENT
  58.  
  59.       INTEGER INFO,INDM,I,J,CPOSREN
  60.       REAL*8 ZERO,ONE,XTOL,DEUXPI
  61.       PARAMETER (ZERO=0.D0,ONE=1.D0,XTOL=1.D-8,
  62.      &           DEUXPI=2.D0*XPI)
  63.  
  64.       SEGINI,mwrkml
  65. *-----------------------------------------------------------------------
  66. *  1. Construction de la matrice de Hill
  67. *     DEL2
  68.       DO I=1,2*NHBM+1
  69.          i_z = NDDL*(I-1)
  70.          DO J=1,NDDL
  71.             DEL2(i_z+J) = ONE/XM(J,1)
  72.          ENDDO
  73.       ENDDO
  74. *     DEL1
  75.        DO J = 1,NT
  76.          DO I = 1,NT
  77.             DEL1(I,J)=ZERO
  78.           ENDDO
  79.        ENDDO
  80.        DO I=1,NDDL
  81.           DEL1(I,I) = XASM(I,1)
  82.        ENDDO
  83.        DO J = 2,2*NHBM,2
  84.          DO I=1,NDDL
  85.              Ci = XASM(I,1)
  86.              Mi = XM(I,1)
  87.              BB = OMEG*J*Mi
  88.              DEL1(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = Ci
  89.              DEL1(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = BB
  90.              DEL1(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -BB
  91.              DEL1(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = Ci
  92.          ENDDO
  93.        ENDDO
  94.  
  95. *-----------------------------------------------------------------------
  96. *  2. Assemblage de JJ
  97. *
  98. *     JJ = [DEL2\DEL1   DEL2\Rx]
  99. *          [     -I         0  ]
  100. *
  101.        DO  I=1,NT
  102.          DO  J=1,NT
  103. *  JJ_11 = DEL2\DEL1
  104.             JJ(I,J) = -DEL1(I,J)*DEL2(J)
  105. *  JJ_12 = DEL2\Rx
  106.            JJ(I,NT+J) = -RX(I,J)*DEL2(J)
  107. *  JJ_21 = [I]
  108.            IF (I.EQ.J) THEN
  109.              JJ(NT+I,J) = ONE
  110.            ELSE
  111.              JJ(NT+I,J) = ZERO
  112.            ENDIF
  113. *  JJ_22 = [0]
  114.            JJ(NT+I,NT+J) = ZERO
  115.          ENDDO
  116.        ENDDO
  117.  
  118. *-----------------------------------------------------------------------
  119. *  3. Calcul des valeurs/vecteurs propres de JJ
  120.       CALL DGEEV('N','V',2*NT,JJ,2*NT,EXPRE,EXPIM,VL,1,VR,2*NT,
  121.      &                                             WORK,8*NT,INFO)
  122.  
  123. *-----------------------------------------------------------------------
  124. *  4. Tri des valeurs propres
  125.       CALL HBMORDO(2*NT,NDDL,OMEG,EXPIM,INDEIG)
  126.  
  127. *-----------------------------------------------------------------------
  128. *  5. Evaluation de stabilite
  129. *     On recupere les exposants de Floquet
  130.       ZSTABN=.TRUE.
  131.       DO  I=1,2*NDDL
  132. cbp : TRIVALS a permute les elements de EXPIM, mais pas EXPRE
  133.          LHILL(1,I) = EXPRE(INDEIG(I))
  134.          LHILL(2,I) = EXPIM(I)
  135.          IF (LHILL(1,I).GT.ZERO) THEN
  136.             ZSTABN=.FALSE.
  137.          ENDIF
  138. *        calcul des Multiplicateurs de Floquet --> deplace par bp
  139.       ENDDO
  140.  
  141. *-----------------------------------------------------------------------
  142. *  6. Detection des bifurcations
  143. * On utilise comme indicateur le nombre d'exposants de Floquet dont la
  144. * partie reelle est positive.
  145.  
  146. *      COMPTAGE -> CPOSREN=?
  147.        CPOSREN = 0
  148.        DO I = 1,2*NDDL
  149. c          IF (LHILL(1,I).LT.ZERO)THEN
  150.          IF (LHILL(1,I).GE.ZERO)THEN
  151.            CPOSREN = CPOSREN+1
  152.          ENDIF
  153.        ENDDO
  154.  
  155.        FLAG = 'S'
  156.        IF (.NOT.ZINIT) THEN
  157. *        WRITE(*,*) 'FTEST=',ABS(CPOSREN-CPOSRE)
  158.  
  159. *         Bifurcations statiques:
  160. *         un seul exposant traverse l'axe imaginaire
  161.           IF (ABS(CPOSREN-CPOSRE).EQ.1) THEN
  162.  
  163.             IF (dv0*dv1.LT.ZERO) THEN
  164. *             Point Limite
  165.               FLAG = 'L'
  166.             ELSE
  167. *             Branchement
  168.               FLAG = 'B'
  169.             ENDIF
  170.  
  171. *         Bifurcations dynamiques:
  172. *         une paire d'exposants traversent l'axe imaginaire
  173.           ELSEIF (ABS(CPOSREN-CPOSRE).EQ.2) THEN
  174.  
  175. c *           Multiplicateurs de Floquet
  176. c             DO I=1,2*NDDL
  177. c cbp           ZR(I) = EXP(DEUXPI*LHILL(1,I)/OMEG)*COS(DEUXPI*LHILL(2,I)/OMEG)
  178. c cbp           ZI(I) = EXP(DEUXPI*LHILL(1,I)/OMEG)*SIN(DEUXPI*LHILL(2,I)/OMEG)
  179. c               mumx(I) = EXP(2.D0*XPI*LHILL(1,I)/OMEG)
  180. c             ENDDO
  181. c             INDM = MAXLOC(mumx,2*NDDL)
  182. cbp: ci-dessus, remplace par (car EXP est une fonction croissante) :
  183.             INDM=1
  184.             XINDM=LHILL(1,INDM)
  185.             DO I=1,2*NDDL
  186.               IF(LHILL(1,I).GT.XINDM) THEN
  187.                 INDM=I
  188.                 XINDM=LHILL(1,I)
  189.               ENDIF
  190.             ENDDO
  191. cbp: calcul de INDM douteux si ce n'est pas la 1ere instabilite -> a verifier + tard
  192.             ZIINDM = EXP(DEUXPI*LHILL(1,INDM)/OMEG)
  193.      &             * SIN(DEUXPI*LHILL(2,INDM)/OMEG)
  194. cbp            IF (ABS(ZI(INDM)).LT.1.E-8) THEN
  195.             IF (ABS(ZIINDM).LT.XTOL) THEN
  196. *             Doublement de periode
  197.               FLAG = 'P'
  198.             ELSE
  199. *             Neimark-Sacker (Hopf secondaire)
  200.               FLAG = 'N'
  201.             END IF
  202.  
  203.           ELSE
  204. c           cas non traite
  205.           ENDIF
  206.  
  207.        ENDIF
  208.  
  209. *-----------------------------------------------------------------------
  210. * 7.  enregistrement pour le pas suivant
  211.       CPOSRE = CPOSREN
  212.       ZSTAB = ZSTABN
  213.  
  214.       END
  215.  
  216.  
  217.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales