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    PV090527  26/09/16    21:15:06     12645          
  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,mwrkml)
  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. C  segini passe a l'etage superieur
  66. *-----------------------------------------------------------------------
  67. *  1. Construction de la matrice de Hill
  68. *     DEL2
  69.       DO I=1,2*NHBM+1
  70.          i_z = NDDL*(I-1)
  71.          DO J=1,NDDL
  72.             DEL2(i_z+J) = ONE/XM(J,1)
  73.          ENDDO
  74.       ENDDO
  75. *     DEL1
  76.        DO J = 1,NT
  77.          DO I = 1,NT
  78.             DEL1(I,J)=ZERO
  79.           ENDDO
  80.        ENDDO
  81.        DO I=1,NDDL
  82.           DEL1(I,I) = XASM(I,1)
  83.        ENDDO
  84.        DO J = 2,2*NHBM,2
  85.          DO I=1,NDDL
  86.              Ci = XASM(I,1)
  87.              Mi = XM(I,1)
  88.              BB = OMEG*J*Mi
  89.              DEL1(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = Ci
  90.              DEL1(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = BB
  91.              DEL1(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -BB
  92.              DEL1(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = Ci
  93.          ENDDO
  94.        ENDDO
  95.  
  96. *-----------------------------------------------------------------------
  97. *  2. Assemblage de JJ
  98. *
  99. *     JJ = [DEL2\DEL1   DEL2\Rx]
  100. *          [     -I         0  ]
  101. *
  102.        DO  I=1,NT
  103.          DO  J=1,NT
  104. *  JJ_11 = DEL2\DEL1
  105.             JJ(I,J) = -DEL1(I,J)*DEL2(J)
  106. *  JJ_12 = DEL2\Rx
  107.            JJ(I,NT+J) = -RX(I,J)*DEL2(J)
  108. *  JJ_21 = [I]
  109.            IF (I.EQ.J) THEN
  110.              JJ(NT+I,J) = ONE
  111.            ELSE
  112.              JJ(NT+I,J) = ZERO
  113.            ENDIF
  114. *  JJ_22 = [0]
  115.            JJ(NT+I,NT+J) = ZERO
  116.          ENDDO
  117.        ENDDO
  118.  
  119. *-----------------------------------------------------------------------
  120. *  3. Calcul des valeurs/vecteurs propres de JJ
  121.       CALL DGEEV('N','V',2*NT,JJ,2*NT,EXPRE,EXPIM,VL,1,VR,2*NT,
  122.      &                                             WORK,8*NT,INFO)
  123.  
  124. *-----------------------------------------------------------------------
  125. *  4. Tri des valeurs propres
  126.       CALL HBMORDO(2*NT,NDDL,OMEG,EXPIM,INDEIG)
  127.  
  128. *-----------------------------------------------------------------------
  129. *  5. Evaluation de stabilite
  130. *     On recupere les exposants de Floquet
  131.       ZSTABN=.TRUE.
  132.       DO  I=1,2*NDDL
  133. cbp : TRIVALS a permute les elements de EXPIM, mais pas EXPRE
  134.          LHILL(1,I) = EXPRE(INDEIG(I))
  135.          LHILL(2,I) = EXPIM(I)
  136.          IF (LHILL(1,I).GT.ZERO) THEN
  137.             ZSTABN=.FALSE.
  138.          ENDIF
  139. *        calcul des Multiplicateurs de Floquet --> deplace par bp
  140.       ENDDO
  141.  
  142. *-----------------------------------------------------------------------
  143. *  6. Detection des bifurcations
  144. * On utilise comme indicateur le nombre d'exposants de Floquet dont la
  145. * partie reelle est positive.
  146.  
  147. *      COMPTAGE -> CPOSREN=?
  148.        CPOSREN = 0
  149.        DO I = 1,2*NDDL
  150. c          IF (LHILL(1,I).LT.ZERO)THEN
  151.          IF (LHILL(1,I).GE.ZERO)THEN
  152.            CPOSREN = CPOSREN+1
  153.          ENDIF
  154.        ENDDO
  155.  
  156.        FLAG = 'S'
  157.        IF (.NOT.ZINIT) THEN
  158. *        WRITE(*,*) 'FTEST=',ABS(CPOSREN-CPOSRE)
  159.  
  160. *         Bifurcations statiques:
  161. *         un seul exposant traverse l'axe imaginaire
  162.           IF (ABS(CPOSREN-CPOSRE).EQ.1) THEN
  163.  
  164.             IF (dv0*dv1.LT.ZERO) THEN
  165. *             Point Limite
  166.               FLAG = 'L'
  167.             ELSE
  168. *             Branchement
  169.               FLAG = 'B'
  170.             ENDIF
  171.  
  172. *         Bifurcations dynamiques:
  173. *         une paire d'exposants traversent l'axe imaginaire
  174.           ELSEIF (ABS(CPOSREN-CPOSRE).EQ.2) THEN
  175.  
  176. c *           Multiplicateurs de Floquet
  177. c             DO I=1,2*NDDL
  178. c cbp           ZR(I) = EXP(DEUXPI*LHILL(1,I)/OMEG)*COS(DEUXPI*LHILL(2,I)/OMEG)
  179. c cbp           ZI(I) = EXP(DEUXPI*LHILL(1,I)/OMEG)*SIN(DEUXPI*LHILL(2,I)/OMEG)
  180. c               mumx(I) = EXP(2.D0*XPI*LHILL(1,I)/OMEG)
  181. c             ENDDO
  182. c             INDM = MAXLOC(mumx,2*NDDL)
  183. cbp: ci-dessus, remplace par (car EXP est une fonction croissante) :
  184.             INDM=1
  185.             XINDM=LHILL(1,INDM)
  186.             DO I=1,2*NDDL
  187.               IF(LHILL(1,I).GT.XINDM) THEN
  188.                 INDM=I
  189.                 XINDM=LHILL(1,I)
  190.               ENDIF
  191.             ENDDO
  192. cbp: calcul de INDM douteux si ce n'est pas la 1ere instabilite -> a verifier + tard
  193.             ZIINDM = EXP(DEUXPI*LHILL(1,INDM)/OMEG)
  194.      &             * SIN(DEUXPI*LHILL(2,INDM)/OMEG)
  195. cbp            IF (ABS(ZI(INDM)).LT.1.E-8) THEN
  196.             IF (ABS(ZIINDM).LT.XTOL) THEN
  197. *             Doublement de periode
  198.               FLAG = 'P'
  199.             ELSE
  200. *             Neimark-Sacker (Hopf secondaire)
  201.               FLAG = 'N'
  202.             END IF
  203.  
  204.           ELSE
  205. c           cas non traite
  206.           ENDIF
  207.  
  208.        ENDIF
  209.  
  210. *-----------------------------------------------------------------------
  211. * 7.  enregistrement pour le pas suivant
  212.       CPOSRE = CPOSREN
  213.       ZSTAB = ZSTABN
  214.  
  215.       END
  216.  
  217.  
  218.  
  219.  

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