Procédures | Cast3M

* FLOQUET   PROCEDUR  BP208322  16/09/01    21:15:03     9010           
*
*=======================================================
*
* FLOQUET
* ANALyse de la stabilite : calcul des exposants de Floquet
*                           (= valeurs propres du determinant de Hill)
* appelé par CONTINU
* CREATION : BP 30/07/2015 (recopie des travaux de these de Lihan Xie)
*
*=======================================================
*
DEBP FLOQUET TAB1;
 
*   MESS '~~~~~  CALCUL DES EXPOSANTS DE FLOQUET  ~~~~~';
 
 
  nhbm = TAB1 . 'N_HARMONIQUE';
  WTAB = TAB1 . 'WTABLE';
  SI (EXIS WTAB 'TIME0'); TIME0 = WTAB . 'TIME0'; FINSI;
* sinon TIME0 est passsé par transitivité depuis continu !!!
  wHz0 = IPOL TAB1 . 'FREQUENCE' TIME0;
 
 
* RECUPERATION DES MATRICES DEJA CONSTRUITES ***************************
 
* KTOT  = (K + KNL)  + (wHz0 * ATOT) ET (wHz0**2 * MTOT);
  KTOT  = WTAB . 'RAIDEUR' ;
  ATOT  = TAB1 . 'AMORTISSEMENT_HBM' ;
  MTOT  = TAB1 . 'MASSE_HBM';
  ATOT1 = TAB1 . 'AMORTISSEMENT_HBM_1' ;
  MTOT1 = TAB1 . 'MASSE_HBM_1';
  MTOT2 = TAB1 . 'MASSE_HBM_2';
 
* * matrices specifiques machines tournantes  
*   OMEG0 = (2.*PI) * wHz0;
*   GTOT  = OMEG0 * (TAB1 . 'CORIOLIS_HBM');
*   KCTOT = OMEG0**2 * (TAB1 . 'CENTRIFUGE_HBM');
*   GTOT1 = OMEG0 * (TAB1 . 'CORIOLIS_HBM_1');
 
 
* CONSTRUCTION DE :                        *****************************
 
* - DELTA0
    DELTA0 = KTOT ;
 
* - DELTA1
    DELTA1 = ATOT1 ET (wHz0 * MTOT1);
 
* - DELTA2
    DELTA2 = MTOT2;
 
    DELTA0 = CHAN DELTA0 'TYPE' 'RIGIDITE';
    DELTA1 = CHAN DELTA1 'TYPE' 'AMORTISS';
    DELTA2 = CHAN DELTA2 'TYPE' 'MASSE';
 
 
 
* CALCUL DES VALEURS PROPRES *******************************************
 
*bp, 2016-06-13 : on pourrait ne sortir que les 2 NDDL premieres v.p.
  BAS_J = VIBC DELTA2 DELTA0 DELTA1 ;
  MBAS_J = BAS_J . 'MODES';
  nmodC = (DIME MBAS_J) - 2;
 
 
* ON RECUPERE LES VALEURS PROPRES DANS 2 LISTREELS *********************
*           + LES DEFORMEE MODALES (normees) DANS 2 LISTCHPO
 
* l = lR + i lI = 2*pi *  i (wR + i wI)
  j = 0;
  lrprog = prog;
  liprog = prog;
  sq_R = suit;
  sq_I = suit;
  tcomp = TAB1 . 'COMPOSANTES'; 
  mophys = tcomp . 'DEPLACEMENT';
  mocomp = tcomp . 'DEPLACEMENT_HBM';
  REPE BMODC nmodC; j = j + 1;
    lrprog = lrprog et (-2.*pi * MBAS_J . j . 'FREQUENCE_IMAGINAIRE');  
    liprog = liprog et ( 2.*pi * MBAS_J . j . 'FREQUENCE_REELLE');
    q_R  = MBAS_J . j . 'DEFORMEE_MODALE_REELLE';
    q_I  = MBAS_J . j . 'DEFORMEE_MODALE_IMAGINAIRE';
    xq2 = (XTX q_R) + (XTX q_I);
    q_R = q_R * (1./(xq2**0.5));
    q_I = q_I * (1./(xq2**0.5));
    sq_R = sq_R et q_R;
    sq_I = sq_I et q_I;
  FIN  BMODC ;
 
* Calcul du cout = distance ente 2 modes complexes + distance entre lambda
  FLORDO = FAUX;
 
* * Calcul de la distance ente 2 modes complexes
*   SI (EXIS WTAB 'PHI_R');
*     Coutphi = prog;
*     REPE BPHI nmodC;
*       p_R = EXTR WTAB . 'PHI_R' &BPHI;
*       p_I = EXTR WTAB . 'PHI_I' &BPHI;
*       REPE BPHJ nmodC;
*         q_R = EXTR sq_R &BPHJ;
*         q_I = EXTR sq_I &BPHJ;
*         xpq_R = (XTY p_R q_R mocomp mocomp) 
*               + (XTY p_I q_I mocomp mocomp);
*         xpq_I = (XTY p_R q_I mocomp mocomp) 
*               - (XTY p_I q_R mocomp mocomp);
* *         mess 'i,j='&BPHI &BPHJ '--> x='xpq_R xpq_I;
*         xpq2 = ((xpq_R**2) + (xpq_I**2))**0.5;
*         Coutphi = Coutphi et xpq2;
*       FIN  BPHJ;
*     FIN  BPHI;
*     Coutphi = 1. - Coutphi;
*     FLORDO=VRAI;
* * *   on ordonne pour minimiser le cout
* *     lpermu = lect 1 pas 1 nmodC;
* *     Xcout lrprog liprog lpermu 
* *     = ORDO 'COUT' Coutphi lrprog liprog lpermu;
*   SINON;
*     Coutphi=0.;
*   FINSI;
 
* on cree la distance Dij = |l_i - l_j|^2 / lref^2
  SI (EXIS WTAB 'LAMBDA_R');
     lrp = WTAB . 'LAMBDA_R';
     lip = WTAB . 'LAMBDA_I';
*    si on peut (pas de point limite),
*    on fait une extrapolation prediction lineaire
     si (exis  WTAB 'LAMBDA_R-1');
       DTM1 = WTAB . 'TIME0-1' - WTAB . 'TIME0-2' ;
       DTM0 = WTAB . 'TIME0'   - WTAB . 'TIME0-1' ;
       si (((ABS DTM1) > 1.E-6) et ((DTM1*DTM0) > 0.));
        DLR1  = lrp -  WTAB . 'LAMBDA_R-1';
        DLI1  = lip -  WTAB . 'LAMBDA_I-1';
        lrp = lrp + (DLR1 * (DTM0/DTM1));
        lip = lip + (DLI1 * (DTM0/DTM1));
       finsi;
     finsi;
*      lref2 = (somm ((lrprog**2) + (liprog**2))) / nmodC;
     Dprog = prog;
*    --> boucle sur les lignes
     repe b_i nmodC;
       xr2 = (lrprog - (extr lrp &b_i))**2;
       xi2 = (liprog - (extr lip &b_i))**2;
       lref2 = (lrprog**2) + (liprog**2) 
             + ((extr lrp &b_i)**2) + ((extr lip &b_i)**2);
       Dprog = Dprog et ((xr2 + xi2) / lref2);
     fin  b_i;
     FLORDO=VRAI;
* *    on ordonne pour minimiser le cout
*      Xcout lrprog liprog = ORDO 'COUT' Dprog lrprog liprog;
  SINON;
    Dprog = 0.;
  FINSI;
 
*   Coutot = Coutphi + Dprog;
  Coutot = Dprog;
 
  SI FLORDO;
*   on ordonne pour minimiser le cout
    lpermu = lect 1 pas 1 nmodC;
    Xcout lrprog liprog lpermu 
    = ORDO 'COUT' Coutot lrprog liprog lpermu;
  SINON;  
*   on ordonne selon reel croissant et |imag| croissant (pour centrer)
    lrprog2 liprog2 = ORDO 'DECROISSANT' lrprog  liprog;
    toto liprog lrprog = ORDO 'CROISSANT' (ABS liprog2) liprog2 lrprog2;
* * imedian pour methode de LAZARUS & THOMAS
*   lrprog2 liprog2 imedian2 = ORDO 'DECROISSANT' lrprog  liprog imedian;
*   toto liprog  lrprog  = ORDO 'CROISSANT' imedian2 liprog2 lrprog2;
  FINSI;
* *   on n'oublie pas d'ordonner les listchpo
*     sq_R0 = COPI sq_R;
*     sq_I0 = COPI sq_I;    
*     REPE BPHI nmodC;
*       REMP sq_R (EXTR lpermu &BPHI) (EXTR sq_R0 &BPHI);
*       REMP sq_I (EXTR lpermu &BPHI) (EXTR sq_I0 &BPHI);
*     FIN BPHI;
 
* on stocke pour prochain pas
  SI (EXIS WTAB 'LAMBDA_R');
    WTAB . 'LAMBDA_R-1' = WTAB . 'LAMBDA_R';
    WTAB . 'LAMBDA_I-1' = WTAB . 'LAMBDA_I';  
    WTAB . 'TIME0-2'    = WTAB . 'TIME0-1';  
  FINSI;
  WTAB . 'LAMBDA_R' = lrprog;
  WTAB . 'LAMBDA_I' = liprog;
  WTAB . 'TIME0-1'  = TIME0;  
*   WTAB . 'PHI_R' = sq_R;
*   WTAB . 'PHI_I' = sq_I;
 
 
 
* on ne conserve ainsi que les (2*nddls) valeurs centrées autour de 0
  si (faux);
  NDDL = (nmodC / 2) / (2*nhbm + 1);
  lprem  = lect 1 pas 1 (2 * NDDL);
  liprog = EXTR liprog lprem;
  lrprog = EXTR lrprog lprem;
  finsi;
 
 
FINP lrprog  liprog;