Accueil | Cast3M

* fichier :  rotor5.dgibi
************************************************************************
************************************************************************
*
******************************************************************
*                                                                *
*  Mots-cl�s : Vibrations, calcul modal, machines tournantes,    *
*              poutre, condensation statique, dynamique          *
*              reconstruction 3D, modes complexes                *
*                                                                *
*  Etude d'une machine tournante dans le repere inertiel (fixe)  *
*  Elements: POUT et volant d inertie                            *
*                                                                *
*  Etude d'un arbre lors d'une rupture d'aube                    *
*                                                                *
*    1- Calcul des inerties et des forces de balourd a l aide    *
*       d un modele 3D massif                                    *
*    2- Etude de la reponse a une rupture accidentelle d'aube    *
*                                                                *
*  Auteur: Didier Combescure Juillet 2003                        *
*  Modif : Benoit Prabel, 06/2014 ajout test DSPR et TFR         *
*                                                                *          
******************************************************************
*COMPLET = VRAI;
 COMPLET = FAUX;
*GRAPH = VRAI; opti trac PSC EPTR 6 POTR HELVETICA_16;
GRAPH = FAUX;
FLAVS = faux;
FLVTK = faux;
*
OPTI 'DIME' 3 'ELEM' SEG2;                                                                
 
******************************************************************
*
* Mod�lisation (maillage, modeles, matrices ...)
*
******************************************************************
*
* Caract�ristique de l'arbre
*
L1 = 10.0;                                                
Re1 = 0.2 ; 
Re1S = 0.22 ;                                                     
E1 = 2.e+11 ;                                                           
Nu1 = 0.2 ;                                                             
Rho1  = 7800. ;
RhoFl = 1000.;
*amortissement faible?
Visc1 = E1/(10000.); 
*amortissement eleve?
*Visc1 = E1/(100.); 
Sarbre = pi*(Re1**2);
Inxarb = pi*(Re1**4)/2.;
Inyarb = pi*(Re1**4)/4.;
*                                                                     
P0 = 0. 0. 0.;
P1 = (0.75*L1)  0.  0. ;
P2 = L1  0.  0. ;
*
Pb0 = 0. 0. 0.;
Pb2 = L1 0. 0.;
GeoRot = (D 75 P0 P1) et (D 25 P1 P2); 
*
* G�om�trie de la turbine
*
OPTI 'DIME' 3 ELEM  CU20;
*
Re2 = 1.0;
Ep2 = 0.05;
vh2 = Ep2 0. 0.;
*
PM1=  P1 plus (0. Re1 0.);
PM2 = P1 plus (0. 0. Re1); 
PM3=  P1 plus (0. ((-1.)*Re1) 0.);
PM4 = P1 plus (0. 0. ((-1.)*Re1));
PMrot =  P1 plus (0. Re1 Re1);
*
PM1e=  P1 plus (0. Re2 0.);
PM2e = P1 plus (0. 0. Re2); 
PM3e=  P1 plus (0. ((-1.)*Re2) 0.);
PM4e = P1 plus (0. 0. ((-1.)*Re2));
*
* Lig2 = (CER3 4 PM1 PM2 PM3) ET (CER3 4 PM3 PM4 PM1);
* Sur2 = SURF LIG2 'PLANE';
* Vol2 = (VOLU Sur2 TRAN (0.5*VH2)) et (VOLU Sur2 TRAN ((-0.5)*VH2));
* elim 0.00001 Vol2;
*
*
PM11e = d 2 PM1 PM1e;
Sur3a = PM11e ROTA 2 90. P0 P2;
Vol3a =  (VOLU 1 Sur3a TRAN (0.5*VH2))
      et (VOLU 1  Sur3a TRAN ((-0.5)*VH2));
DEPL Vol3a TOUR 25. P1 PMrot;
elim 0.001 Vol3a;
Vol3b = Vol3a TOUR 90. P0 P2;
Vol3c = Vol3b TOUR 90. P0 P2;
Vol3d = Vol3c TOUR 90. P0 P2;
VOL3 = Vol3a et Vol3b et Vol3c et Vol3d;
VOL3Bal = Vol3a et Vol3b et Vol3c;
*
* Modele de l'arbre
*        
Mod1 = MODELI GEOROT 'MECANIQUE' POUT;                                          
Mat1 = MATE Mod1 'YOUN' E1 'NU  ' Nu1 'RHO ' Rho1
       'SECT' Sarbre 'INRY' Inyarb 'INRZ' Inyarb
*        'TORS' Inxarb 'OMEG' (1.0) 'VISQ' (Visc1);    
       'TORS' Inxarb  'VISQ' (Visc1);    
*
MESH3D = POUT2MAS Mod1 Mat1 'MASSIF';
SI (GRAPH);
  OEIL3d = 100. 100. 100.;
  ELIM 0.0001 MESH3D;
  Titre 'Maillage poutre';
  trac OEIL3d MESH3D CACH;
  trac OEIL3d (MESH3D et GeoRot);
  Titre 'Maillage poutre + volant';
  trac OEIL3d (MESH3D et VOL3) CACH;
  oeilx = 100. 0. 0.;
  Titre 'Volant';
  trac oeilx VOL3 cach;
  trac oeilx VOL3Bal cach;
FINSI;
*
* Modele de la turbine
*
Mod2 = MODELI Vol3 'MECANIQUE' 'ELASTIQUE'  ;                                        
Mat2 = MATE Mod2 'YOUN' E1 'NU  ' Nu1 'RHO ' Rho1;
*
Mod2B = MODELI Vol3Bal 'MECANIQUE' 'ELASTIQUE'  ;                                        
Mat2B = MATE Mod2B 'YOUN' E1 'NU  ' Nu1 'RHO ' Rho1;
*
*
Blo2 = (BLOQ 'UZ' P0) ET
       (BLOQ 'UY' P0) ET
       (BLOQ 'UX' P0);   
 
Blo3 = (BLOQ 'UZ' P2) ET
       (BLOQ 'UY' P2) ET
       (BLOQ 'UX' P2);                             
Blo1 = BLOQ 'RX' Georot;
*
* Matrice de masse
*                                                                      
MMass1 = MASS Mod1 Mat1 ;                                   
MRigi1 = RIGI Mod1 Mat1 ;   
*
* Condensation de la matrice de masse du volant d'inertie sur 1 point
*
MMass2 = MASS Mod2 Mat2 ;
MRigi2 = RELA CORI DEPL ROTA (P1 et Vol3);
*
SupK2 = SUPER 'RIGIDITE' 'NOMU' MRigi2 P1;
SupM2 = SUPER 'MASSE' SupK2 MMass2;
*
InxVol = EXTR (PSCAL (SupM2 *(MANU CHPO P1 1 RX 1.)) 
       (MANU CHPO P1 1 RX 1.) (MOTS MX) (MOTS RX)) 'SCAL' P1;
InyVol = EXTR (PSCAL (SupM2 *(MANU CHPO P1 1 RY 1.)) 
       (MANU CHPO P1 1 RY 1.) (MOTS MY) (MOTS RY)) 'SCAL' P1;
*  
*
MMassT = MMass1 et SupM2;
*
Geovolan = MANU POI1 P1;
CG2 = MANU RIGI 'TYPE' 'AMORTISSEMENT' Geovolan
      ('MOTS' 'RX  ' 'RY  ' 'RZ  ') 'QUEL'
      (prog 0. 0. 0.
            0. 0. (InyVol)
            0. ((-1.)*InyVol) 0.);
*
* Matrice de couplage gyroscopique de l'arbre
*
 CG1 = GYROSCOPIQUE Mod1 Mat1;
 CGTOT = CG1 et CG2;
* 
* Matrice d'amortissement
*
 CDAM1 CDAH1 = 'AMOR' Mod1 Mat1 'COROTATIF';   
*
*
* Calcul des forces de balourd (pour une vitesse de rotation de 1 tour/s)
*
*   D�finition du champs d'acceleration centrifuge
*
XYZ =  (0.*(CHAN COMP (COOR 1 Vol3Bal) 'UX')) ET
       (CHAN COMP (COOR 2 Vol3Bal) 'UY') ET
       ((CHAN COMP (COOR 3 Vol3Bal) 'UZ'));
*
XYZ0 = MANU CHPO Vol3Bal 3 'UX' (0.) 'UY' (0.)
                         'UZ' (0.);
*
* Calcul avec la matrice de raideur centrifuge
*
OmegRot = 1.0;
Vrota = OmegRot 0. 0.;
KCE2 ='KENT' 'CENTRIFUGE' Mod2b Mat2b Vrota;
FCENT2 = KCE2*(XYZ - XYZ0);
*
* Calcul avec la matrice de masse
*
FCENT = (OmegRot**2)*MMass2*(XYZ - XYZ0);
*
SI (GRAPH);
  TRAC oeilx Vol3Bal  (vecteur FCENT FX FY FZ 0.01);
FINSI;
*
* Condensation des forces centrifuges pour le calcul du balourd
*
Fbalourd = SUPER 'CHARGE' SupK2 FCENT;
FbalR = Fbalourd;
mess 'FbalR';
list FbalR;
Vecrot = MANU CHPO (extr FbalR mail) 3 
              RX 1. RY 0. RZ 0.;
FbalI = (FbalR '*'VecRot ('MOTS' FY FZ FX) 
                ('MOTS' RZ RX RY) ('MOTS' FX FY FZ)) 
     -  (FbalR '*'VecRot ('MOTS' FZ FX FY) 
                ('MOTS' RY RZ RX) ('MOTS' FX FY FZ));
FbalI = CHAN ATTRIBUT FbalI NATURE DISCRET;
 
******************************************************************
*
*   Calcul des modes propres
*
******************************************************************    
TBasR1 = VIBR 'PROCHE' (prog 10.) (lect 2)
           (MRigi1 ET Blo1 ET Blo2 et Blo3) 
             MMassT  ;
*  
TMod1 = TBasR1 . 'MODES' ;
NbMod1 = (DIME TMod1) - 2 ;
TInd1 = INDEX (TMod1 . 1) ;
*
SI (GRAPH);
 REPETER bbb1 NbMod1 ;
  i = &bbb1 ;
  Frq1 = TMod1 . i . 'FREQUENCE' ;
  Def1 = TMod1 . i . 'DEFORMEE_MODALE' ;
  TITRE Frq1 ;
  trac oeil3d (defo 0.5 georot def1);
 FIN bbb1 ;
FINSI;
Frq1 = TMod1 . 1 . 'FREQUENCE' ;
mess Frq1;
*
*
  MAmor1 = CDAM1;
*
*-Projection des matrices assembl�es sur la base r�elle
*
  MMasP1 = PJBA TBasR1 MMassT ;
*   MRigP1 = PJBA TBasR1 (MRigi1 ET Blo1 ET Blo2 et Blo3)  ;
  MRigP1 = PJBA TBasR1 MRigi1  ;
  MAmoP1 = PJBA TBasR1 MAmor1 ;
  MAmHP1 = PJBA TBasR1 CDAH1;
  Cg1P = PJBA TBasR1 CgTot ;
*
******************************************************************
*
* Calcul des modes complexes avec VIBC 
*
******************************************************************
*
* Vitesse de rotation (10 Hz)
* 
* Omeg1Hz = 0.;
Omeg1Hz = 10.;
Omeg1 =Omeg1Hz*2.*pi;
*
   BasCi1 = VIBC MMasP1 (MRigP1 et (Omeg1*MAmHP1)) 
                        (MAmoP1 et (Omeg1*Cg1P)) TBasR1 ;
*
  ModBasC = BasCi1.'MODES';
  NbMod1 = (DIME ModBasC) - 2;
*
FRQR = TABLE;
FRQI = TABLE;
FRQ = TABLE;
AMOR= TABLE;
*
mess 'Frequences et amortissement a ''Hz';
mess '__________________________________________________';
mess 'mode  |     w_R     |     w_I     |     xi';
* 
  REPETER bbb2 NbMod1 ;
   i2 = &bbb2 ;
   FrqR1 = BasCi1.'MODES'.i2.'FREQUENCE_REELLE' ;
   FrqI1 = BasCi1.'MODES'.i2.'FREQUENCE_IMAGINAIRE' ;
   DefR1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_REELLE';
   DefI1 = BasCi1.'MODES'.i2.'DEFORMEE_MODALE_IMAGINAIRE';
   FRQR.i2 = FrqR1;
   FRQI.i2 = FrqI1;
   FRQ.i2 = ((FrqR1 ** 2) + (FrqI1 ** 2)) ** 0.5 ;
   AMOR.i2   = FrqI1 /(FRQ.i2 );
   mess i2 (FrqR.i2) (FrqI.i2) (AMOR.i2);
  FIN bbb2;
mess '__________________________________________________';
*
 
 
******************************************************************
*
* Calcul du diagramme de Campbell 
*
******************************************************************
 
TAB1 = TABLE;
TAB1.'BASE_MODALE' = TBasR1;
TAB1.'MASS_PROJ' = MMasP1;
TAB1.'RIGI_PROJ' = MRigP1;
TAB1.'GYRO_PROJ' = 2.*pi*Cg1P;
TAB1.'AMOR_PROJ' = MAmoP1;
TAB1.'KROT_PROJ' = (2.*pi*MAmHP1); 
TAB1.'AFFICHAGE' = VRAI;
 
PRFREQ = prog 0. 0.1 0.5 1. pas 1. 50.;
 
CAMPBELL TAB1 PRFREQ;
*
EVFRb1 = (TAB1. 1). 'FREQUENCE_REELLE' coul TURQ;
EVFRb2 = (TAB1. 2). 'FREQUENCE_REELLE' coul ROSE;
EVFRb3 = (TAB1. 3). 'FREQUENCE_REELLE' coul BLEU;
EVFRb4 = (TAB1. 4). 'FREQUENCE_REELLE' coul ROUG;
 
EVAMb1 = (TAB1. 1). 'AMORTISSEMENT' coul TURQ;
EVAMb2 = (TAB1. 2). 'AMORTISSEMENT' coul ROSE;
EVAMb3 = (TAB1. 3). 'AMORTISSEMENT' coul BLEU;;
EVAMb4 = (TAB1. 4). 'AMORTISSEMENT' coul ROUG;
*
EvFreOme = evol manu (prog 0. 50.) (prog 0. (50.));
EvFreZer = evol manu (prog 0. 50.) (prog 0. 0.);
*
SI (GRAPH);
 DESS (EVFRb1 et EVFRb2 et EVFRb3 et EVFRb4 et EvFreOme)
 'TITX' '\W (Tr/s)' POSX CENT
 'TITY' 'f (Hz)' POSY CENT YBOR 0. 20.
 'TITR' 'Campbell diagram (red:Forward - blue:Backward)';
 DESS (EVAMb1 et EVAMb2 et EVAMb3 et EVAMb4 et EvFreZer)
 'TITX' '\W (Tr/s)'  POSX CENT
 'TITY' '\z'  POSY CENT
 'TITR' 'Campbell diagram (red:Forward - blue:Backward)';
finsi;
 
wr50F = ipol 50. EVFRb3;
wr50B = ipol 50. EVFRb4;
mess wr50F wr50B ;
 
******************************************************************
*
* Calcul de la reponse stationnaire au balourd
*
******************************************************************
 
SI COMPLET;
 promeg = (2.*pi) * (prog 1. pas 0.1 50.);
SINON;
 promeg = (2.*pi) * (prog 1. pas 1. 6 pas 0.2 10 pas 1. 50.);
FINSI;
*
TAB2 = TABLE;
TAB2.'BASE_MODALE' = TBasR1 ;
TAB2.'MASS_PROJ' = MMasP1;
TAB2.'RIGI_PROJ' = MRigP1;
TAB2.'GYRO_PROJ' = Cg1P;
TAB2.'AMOR_PROJ' = (MAmoP1);
TAB2.'KROT_PROJ' = (MAmHP1); 
TAB2.'AFFICHAGE' = VRAI;
TAB2.'FBALOURD' = Fbalourd;
TAB2.'SAUVDEFO' = FAUX;
TAB2.'REPONSE' = TABLE;
TAB2.'REPONSE' . 1  = TABLE;
TAB2.'REPONSE' . 1 . 'POINT' = P1;
 
BALOURD TAB2 promeg;
*
EvUYbal = EVOL VIOL MANU (promeg/(2.*pi))
        (extr (((TAB2.'REPONSE'). 1) . 'UY' ) ordo);
EvUZbal =  EVOL BLEU MANU (promeg/(2.*pi))
       (extr  (((TAB2.'REPONSE'). 1) . 'UZ' ) ordo);
EvGabal =  EVOL MANU (promeg/(2.*pi))
        (extr (((TAB2.'REPONSE'). 1) . 'GAXE') ordo);
EvPabal =  EVOL MANU (promeg/(2.*pi))
        (extr (((TAB2.'REPONSE'). 1) . 'PAXE') ordo);
EvAngEl =  EVOL MANU (promeg/(2.*pi))
        (extr (((TAB2.'REPONSE'). 1) . 'ANGLEGAXE') ordo);
*
SI (EGA GRAPH 'Y');
 DESS (EvUYbal et EvUZbal)
 GRIL 'POIN' GRIS
 TITY '|u|' LOGY POSY CENT
 TITX '\W (tr/s)' POSX CENT
 XBORD 0. 50.;
FINSI;
 
uybal5 = ipol 5. EvUYbal;
uzbal5 = ipol 5. EvUZbal;
mess uybal5 uzbal5;
 
 
*******************************************************
*
*               Calculs transitoires
*
*******************************************************
 
* vitesse de rotation en Hz
  fr1 = 5.0;
  Tperiod = 1./fr1;
* nombre de periodes
  SI COMPLET;  Nperiod = 20;
  SINON;       Nperiod = 7;
  FINSI;
 
*******************************************************
*
*                     DYNAMIC
*
*******************************************************
dt = (1./(50.*fr1));
tdebut = (20.*dt);
* tduree = Nperiod * 30. * dt;
tduree = Nperiod * Tperiod - tdebut;
 
prtime2 = prog 0. pas dt (tduree + (100.*dt));
prtime = (prog 0.) et
         (prtime2 + (prog (dime prtime2)*tdebut));
prcren = (prog 0.) et (prog (dime prtime2)*1.);
prcos = (prog 0. ) et (cos (fr1*prtime2*360.));
prsin = (prog 0. ) et (sin (fr1*prtime2*360.));
prtimef = prog 0. pas dt  (tdebut+tduree);
 
* cahrgement de balourd
EvReel = evol manu prtime 
   ((prcos*prcren)*((2.*pi*fr1)**2));
EvImag = evol manu prtime
   ( (((-1.)*prsin*prcren)*((2.*pi*fr1)**2)));
*
ChaBal1 = (CHAR FbalR EvReel);
Chabal2 = (CHAR FbalI EvImag);
Chabal = Chabal1 et Chabal2;
*
LIS1 = prtimef;
NTT = ((dime prtimef) - 5);
TAB3in = TABLE ;
TAB3in.'CHAR' = ChaBal ;
TAB3in.'DEPL' = MANU CHPO GEOROT 3 UX 0. UY 0. UZ 0. ;
TAB3in.'VITE' = MANU CHPO GEOROT 3 UX 0. UY 0. UZ 0. ;
TAB3in.'RIGI' =  (MRigi1 ET Blo1 ET Blo2 et Blo3) 
            et (0.*CDAH1*2.*pi*fr1);
TAB3in.'MASS' =  MMassT;
TAB3in.'AMOR' =  MAmor1 et (CgTot*2.*pi*fr1);
TAB3in.'FREQ' =  (1./dt);
TAB3in.'INST' = LIS1 ;
TAB3 = DYNAMIC TAB3in ;
*
Time = PROG  ;
PRUY1 = PROG  ;
PRUZ1 = PROG  ;
I = 0 ;
REPE  BB1 NTT  ;
  I = I + 1 ;   
  T_I = TAB3 . I . 'TEMP' ;
  Time = Time ET T_I;                                           
  DEP1 = (TAB3 . I . 'DEPL') exco (mots 'UX' 'UY' 'UZ') ;
  PRUY1 = PRUY1 ET( PROG (EXTR DEP1 UY  P1) );                            
  PRUZ1 = PRUZ1 ET( PROG (EXTR DEP1 UZ  P1) ); 
* sortie paraview
  si (FLAVS);
    ficsort = chai 'AVS/rotor5_DYNAMIC_' I '.inp';
    OPTI SORT ficsort;
    mess ficsort;
    SORT 'AVS' GEOROT  DEP1 'TEMP' T_I;                        
  fins;
* sortie paraview
  si (FLVTK);
    mess 'I=' I '  T_I=' T_I;
    ficsort = chai 'AVS/rotor5_DYNAMIC'  ;
    OPTI SORT ficsort;
    mess ficsort;
    SORT 'VTK' GEOROT  DEP1 'TEMP' T_I;                        
  fins;
FIN BB1;
* opti donn 5;
*
EvUY1 = evol VIOL manu Time PRUY1;
EvUZ1 = evol BLEU manu Time PRUZ1;
EvUYUZ1 = EVOL BLEU MANU PRUY1 PRUZ1;
*
SI (GRAPH);
 tit1 = mot 'transitoire DYNAMIC';
 tdess1 = tabl;
 tdess1 . 1 = mot 'MARQ PLUS REGU'; 
 tdess1 . 2 = mot 'TIRR MARQ CARR REGU'; 
 tdess1 . 3 = mot 'MARQ CROI REGU'; 
 tdess1 . 4 = mot 'TIRR MARQ ROND REGU'; 
 tdess1 . 'TITRE' = tabl;
 tdess1 . 'TITRE' . 1 = mot 'UY (implicit)';
 tdess1 . 'TITRE' . 2 = mot 'UZ (implicit)';
 tdess1 . 'TITRE' . 3 = mot 'UY (explicit)';
 tdess1 . 'TITRE' . 4 = mot 'UZ (explicit)'; 
 DESS (EvUY1 et EvUZ1)
 TITY 'U'        POSY CENT
 TITX '\W (tr/s)' POSX CENT
 'TITR' tit1 tdess1 LEGE NO;
 DESS EvUYUZ1
  TITY 'UZ'        POSY CENT
  TITX 'UY'        POSX CENT
  'CARRE' 'TITR' tit1 ;
FINSI;
 
 
*** calcul de la DSP ***
 
*... via l'operateur DSPR
nduree = (tdebut+tduree) / dt ;
nduree = enti ((log nduree) / (log 2));
DSP_UY1 = DSPR nduree EvUY1 'VIOL';
DSP_UZ1 = DSPR nduree EvUZ1 'BLEU';
 
SI (GRAPH);
  dess (DSP_UY1 et DSP_UZ1) 
  GRIL POIN GRIS
  XBOR 0 50. TITX 'f(Hz)'  POSX CENT
  LOGY TITY 'DSP(U)' POSY CENT 
  TITR 'Implicit + physical base' tdess1 LEGE NE; 
FINSI;
 
*... via la Transformee de Fourier Rapide de l'evolution
Tduree = dt * (2**nduree); mess Tduree;
TF_UY1 = TFR nduree EvUY1 'REIM' 'VIOL';
prfreq = TF_UY1 EXTR 'ABSC' ;
prUY1R = TF_UY1 EXTR 'ORDO' 1;
prUY1I = TF_UY1 EXTR 'ORDO' 2;
S_UY1 = (2./Tduree) * ((prUY1R**2) + (prUY1I**2));
S_UY1 = EVOL 'VIOL' 'MANU' 'Hz' prfreq '|UY|^{2}' S_UY1;
 
*comparaison
DSP_som = (INTG ((DSP_UY1 extr 'AVAN' 10.) extr 'APRE' 2.5)) ;
S_som   = (INTG ((S_UY1 extr 'AVAN' 10.) extr 'APRE' 2.5))  ;
mess DSP_som S_som;
 
 
**********************************************************
*
*           DYNE
*
**********************************************************
TRIG = TABLE 'RAIDEUR_ET_MASSE';
TRIG.'RAIDEUR' = MRigP1;
TRIG.'MASSE'  = MMasP1;
*
TAMOR = TABLE 'AMORTISSEMENT';
TAMOR . 'AMORTISSEMENT' = MAmoP1;
*
*
CGP = (PJBA TBasR1 (CgTot*2.*pi*fr1));
*
FLPLEIN = vrai ;
si (FLPLEIN);
 
* Couplage gyroscopique => mis au 1er membre
*
TAMOR . 'AMORTISSEMENT' = MAmoP1 et CGP;
TAMOR . 'NATURE'        = mot 'PLEINE';
 
 
sinon;
 
* Couplage gyroscopique => mis au second membre
*
TLIAI = TABLE 'LIAISON'  ;
TTLA = TABLE 'LIAISON_A' ;
TLIAI.'LIAISON_A' = TTLA ;
*
Nmode = (DIME (TBASR1.'MODES')) - 2;
I=0;
Iliai = 0;
REPETER BLOC1 Nmode;
 I=I+1;
 J=0;
REPETER BLOC2 Nmode;
 J=J+1;
 Iliai = Iliai + 1;
*
 PtI = ((TBASR1.'MODES'). I). 'POINT_REPERE'; 
 PtJ = ((TBASR1.'MODES'). J). 'POINT_REPERE'; 
*
 Vunit = MANU CHPO PtI 1 'ALFA' 1.;
 CHPOM1 = CGP*Vunit;
 CoefCG1 = EXTR CHPOM1 'FALF' PtJ;
 list CoefCG1;
 TLIAI.'LIAISON_A'.Iliai=TABLE 'LIAISON_ELEMENTAIRE';
 TLIAI.'LIAISON_A'.Iliai.'TYPE_LIAISON'= MOT 'COUPLAGE_VITESSE';
 TLIAI.'LIAISON_A'.Iliai.'SUPPORT'=PtJ;
 TLIAI.'LIAISON_A'.Iliai.'ORIGINE'=PtI;
 TLIAI.'LIAISON_A'.Iliai.'COEFFICIENT'= (-1.)*CoefCG1;
*
FIN BLOC2;
FIN BLOC1;
*
finsi;
 
TCHAR = TABLE 'CHARGEMENT' ; 
CHAP1 = PJBA  TBasR1 ChaBal;     
TCHAR.'BASE_A' = CHAP1;
*
TSORT = TABLE 'SORTIE'   ;                           
TSORV = TABLE 'VARIABLE' ;                           
TSORT.'VARIABLE' = TSORV ;                           
TSORV.'VITESSE'     = VRAI;                     
TSORV.'DEPLACEMENT' = VRAI ;                     
TSORV.'ACCELERATION'= VRAI ;
* bp:
TSORV . 'TYPE_SORTIE' =  mot 'LISTREEL';
*
DTEX = DT/2;
si (FLPLEIN);
TRESU = DYNE DE_VOGELAERE TBasR1 TAMOR TCHAR
             (NTT*2) DTEX TSORT  ;
sinon;
TRESU = DYNE DE_VOGELAERE TBasR1 TAMOR TCHAR TLIAI
             (NTT*2) DTEX TSORT  ;
finsi;
 
******************************************************************
*         TRACE DES DEPLACEMENTS EN FONCTION DU TEMPS
******************************************************************
 
Ev2UY1 = EVOL ROSE RECO TRESU TBasR1 'DEPL' P1  UY  ;
Ev2UZ1 = EVOL ROUG RECO TRESU TBasR1 'DEPL' P1  UZ  ;
Ev2UYUZ1 = EVOL ROUG MANU (extr Ev2UY1 ordo) (extr Ev2UZ1 ordo);
*
SI (GRAPH);
 tit2 = mot 'transitoire DYNAMIC-DYNE';
 tdess2 = tabl;
 tdess2 . 1 = mot 'MARQ PLUS REGU'; 
 tdess2 . 2 = mot 'TIRR MARQ CROI REGU'; 
 tdess2 . 'TITRE' = tabl;
 tdess2 . 'TITRE' . 1 = mot 'DYNAMIC (implicit)';
 tdess2 . 'TITRE' . 2 = mot 'DYNE (explicit)';
 DESS (EvUY1 et EvUZ1 et Ev2UY1 et Ev2UZ1)
 TITY 'U'        POSY CENT
 TITX '\W (tr/s)' POSX CENT
 'TITR' tit1 tdess1 LEGE NO;
 DESS (EvUYUZ1 et Ev2UYUZ1)
  TITY 'UZ'        POSY CENT
  TITX 'UY'        POSX CENT
  'CARRE' 'TITR' tit1 tdess2 LEGE NO;
FINSI;
*
UyMa= maxi (extr EvUY1 ordo) abs;
Uy2Ma= maxi (extr Ev2UY1 ordo) abs;
mess UyMa Uy2Ma;
 
*** calcul de la DSP ***
 
NTT2 = enti ((log (NTT*2)) / (log 2));
DSP_UY2 = DSPR NTT2 Ev2UY1 'ROUG';
DSP_UZ2 = DSPR NTT2 Ev2UZ1 'ROSE';
DSP_som2 = (INTG ((DSP_UY2 extr 'AVAN' 10.) extr 'APRE' 2.5)) ;
mess DSP_som S_som DSP_som2;
 
tdess3 = tabl;
tdess3 . 'TITRE' = tabl;
tdess3 . 'TITRE' . 1 = mot 'UY (explicit)';
tdess3 . 'TITRE' . 2 = mot 'UZ (explicit)';
 
SI (GRAPH);
  dess (DSP_UY2 et DSP_UZ2) 
  GRIL POIN GRIS
  XBOR 0 50. TITX 'f(Hz)'  POSX CENT
  LOGY TITY 'DSP(U)' POSY CENT
  TITR 'Explicite base modale' tdess3 LEGE NE; 
FINSI;
 
 
*** On va vérifier que EVOL RECO avec un maillage fonctionne aussi ***
 
evgeoUY = EVOL VERT 'RECO' TRESU TBasR1 'DEPL' GEOROT  'UY'  ;
GEO1 = CHAN GEOROT 'POI1';
iP1 = POSI P1 'DANS' GEO1;
evgeoUY1 = EXTR evgeoUY 'COUR' iP1;
*test
TestEVOLRECO = MAXI 'ABS' ((EXTR Ev2UY1 'ORDO') - (EXTR evgeoUY1 'ORDO'));
SI (TestEVOLRECO > 1.E-12); ERRE 5 ; FINSI; 
 
* on recommence avec le point suivant
iP1next = iP1 + 1;
P1next = POIN GEO1 iP1next;
Ev2UY1next = EVOL ROSE RECO TRESU TBasR1 'DEPL' P1next  UY  ;
evgeoUY1next = EXTR evgeoUY 'COUR' iP1next;
*test
TestEVOLRECO = MAXI 'ABS' ((EXTR Ev2UY1next 'ORDO') - (EXTR evgeoUY1next 'ORDO'));
SI (TestEVOLRECO > 1.E-12); ERRE 5 ; FINSI; 
 
 
 
******************************************************************
*                    Test de bon fonctionnement                  *
******************************************************************
 
***  valeurs theoriques  ***
 
*mode au repos (1ere frequence)
 FR1Th = 7.496;
*campbell (frequence a 50Hz)
 wr50Fth = 7.3838;  wr50Bth = 7.6115;
*balourd stationnaire a 5Hz
 uy5th = 1.34618E-02  ; uz5th = 1.37074E-02;
*transitoire : comparaison DYNAMIC, DYNE
 UyTh = 1.902E-02; Uy2Th = 1.820E-02;
*verif DSP entre 2.5 et 10Hz 
 DSPth = 1.E-04;
 
***  ecarts relatif et test d'erreur ***
TEST1 = (ABS (Frq1 - FR1Th)) >eg (0.01*FR1Th);
TEST2 = ( (ABS (wr50F - wr50Fth)) >eg (0.01*wr50Fth) )
     ou ( (ABS (wr50B - wr50Bth)) >eg (0.01*wr50Bth) );
TEST3 = ( (ABS (uybal5 - uy5th)) >eg (0.01*uy5th) )
     ou ( (ABS (uzbal5 - uz5th)) >eg (0.01*uz5th) );
TEST4 = ( (ABS (UyMa  - UyTh )) >eg (0.01*UyTh ) )
     ou ( (ABS (Uy2Ma - Uy2Th)) >eg (0.01*Uy2Th) );
si (COMPLET);     
  TEST5 = ( (ABS (DSP_som - S_som)) >eg (0.0001*DSPth) );
sino;
  TEST5 = ( (ABS (DSP_som - S_som)) >eg (0.0001*DSPth) )
      ou ( (ABS (DSP_som - DSPth)) >eg (0.10*DSPth) )
      ou ( (ABS (DSP_som2 - DSPth)) >eg (0.10*DSPth) );
finsi;
 
* opti donn 5 trac X;
SI (TEST1 ou TEST2 ou TEST3 ou TEST4 ou TEST5);
   ERRE 5;
SINON;
   ERRE 0; 
FINSI;
*
 
 
*opti donn 5 trac X;
 
FIN ;