Accueil | Cast3M

* fichier : rotor6.dgibi
************************************************************************
************************************************************************
*
******************************************************************              
*                                                                *
*  Mots-clés : Vibrations, calcul modal, machines tournantes,    *
*              3D, flambage, precontrainte                       *
*                                                                *
*   rotor6.dgibi                                                 *
*                                                                *
*   Cas test basé sur rotor2.dgibi [Exemple Lalanne P.13]        *
*   Modélisation 3D dans le repère tournant                      *
*   Géo : épaisseur du disque + fine (hvol)                      *
*                                                                *
*   B. Prabel, janvier 2012 :                                    *
*   -teste les matrices en 3D massif créée par KENT KSIG         *
*    (CORI a ajouter a terme)                                    *
*   -teste l'algo de VIBR 'SIMUL' dans les cas ou :              *
*     + on cycle (grand nombre de modes)                         *
*     + K ou M n'est pas definie positive                        *
*                                                                *
******************************************************************              
 
 
******************************************************
*** Options ***
******************************************************
COMPLET = FAUX;
* COMPLET = VRAI;
 
opti dime 3 elem cub8 ;
*opti dime 3 elem cu20 ;
 
* tracés
OPTI TRAC PSC EPTR 4;
eye1 = 5. -10. 3. ;
 
*message dans VIBR
iimpi = 6;
* iimpi = 5;
* iimpi = 4;
iimpi = 0;
 
*test de DIAG
FLDIAG = FAUX;
 
 
******************************************************
*** Maillage ***
******************************************************
*
Ltot = 0.40 ;  
L1 = Ltot/3.;
L2 = 2.*Ltot/3.;
*                                                        
Re1 = 0.01 ;
Ri2 = 0.01 ;
Re2 = 0.15 ;
* hvol = 0.03;
* hvol = 0.005;
hvol = 0.002;
 
*
* Sections circulaires de l'arbre et du volant d'inertie
p0 = 0. 0. 0.;
vaxe = 1. 0. 0.;
p1 =  0. Re1 0. ;
p1e = 0. Re2 0. ;
p2 =  0. 0.  Re1 ;
p3 =  0. ((-1.)*Re1) 0. ;
p4 =  0. 0. ((-1.)*Re1) ;
* n1 = 2;
  n1 = 3;
* n1 = 4;
  n2 = 5;
*   naxe1 = 8; naxe2 = 2; naxe3 = 8; naxe4 = 8;
  naxe1 = 4; naxe2 = 2; naxe3 = 4; naxe4 = 4;
fins;
*
p1axe = p0;
p2axe = p0 plus ((L1 - (0.5*hvol))*vaxe);
p3axe = p0 plus ((L1 + (0.5*hvol))*vaxe);
p4axe = p0 plus ((L2)*vaxe);
p5axe = p0 plus ((Ltot)*vaxe);
pbal = p1e  plus ((L1 - (0.5*hvol))*vaxe);
 
*bp : modif du maillage de la section de l arbre
  p1p = p0 PLUS (0.2*(p1 MOIN p0)) ;
  p2p = p0 PLUS (0.2*(p2 MOIN p0)) ;
  p3p = p0 PLUS (0.2*(p3 MOIN p0)) ;
  p4p = p0 PLUS (0.2*(p4 MOIN p0)) ;
  cer1 = (cer3 (2*n1) p1 p2 p3) ;
  cer1p = (cer3 (2*n1) p1p p2p p3p) ;
  cer2 = (cer3 (2*n1) p3 p4 p1) ;
  cer2p = (cer3 (2*n1) p3p p4p p1p) ;
  secarbre = (REGL cer1 cer1p) et (COUT cer1p p0)
          et (REGL cer2 cer2p) et (COUT cer2p p0) ;
  elim secarbre (1.E-4*Re1) ;        
*
secvolan = (rota 360. (4*n1) (d n2 p1 p1e) p0 vaxe) et secarbre;
ELIM secvolan (1.E-4*Re1) ;
*
*opti elem cub8;
vol1 = secarbre volu naxe1 tran  ((L1 - (0.5*hvol))*vaxe);
vol2 = (secvolan plus ((L1 - (0.5*hvol))*vaxe))
        volu naxe2  tran (hvol*vaxe);
vol3 = (secarbre plus ((L1 + (0.5*hvol))*vaxe))
        volu naxe3  tran ((L2 - (L1 + (0.5*hvol)))*vaxe); 
vol4 = (secarbre plus (L2*vaxe))
        volu naxe4  tran ((Ltot - L2)*vaxe); 
voltot = vol1 et vol2 et vol3 et vol4;
elim 0.0001 (voltot et p1axe et p2axe et 
              p3axe et p4axe et p5axe et pbal);
*
* tracé
rp1 = @REPERE (prog 0.5 0.2 0.2) 'JAUN';
trac eye1 (voltot et rp1) cach;
trac eye1 (voltot) cach;
 
 
 
******************************************************************      
*** Données Matériaux et modèle ***
******************************************************************      
 
* arbre et disque en acier
Ey1 = 2.E+11 ;                
Nu1 = 0.3 ;                                                             
Rho1  = 7800. ;
Visc1 = 0.00015*Ey1;
 
* modele et materiau
modtot = MODE voltot 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE';
mattot = MATE modtot 'YOUN' Ey1 'NU' Nu1 'RHO' Rho1 'VISQ' Visc1;
 
*axe de rotation
* vaxe deja defini
 
 
***************************************************
*** Définition des matrices a l'aide de BRASERO ***
***************************************************
*matrices 
RIG1 = RIGI modtot mattot;
MAS1 = MASS modtot mattot;
AMO1 = AMOR modtot mattot;
KCE1 = KENT 'CENTRIFUGE' modtot mattot vaxe 'FLAM';
* on met 'FLAM' pour avoir le type MASSE sachant que ri1 ET ri2 (via fusrig)
* attribue le type RIGIDITE par defaut si ri1.mtymat different de ri2.mtymat
 
*CL
CL1 = BLOQ (p1axe ET p5axe) 'UX' 'UY' 'UZ';
*blocage de la rotation de corps rigide
pbal2 = voltot POIN PROC (pbal MOINS (0. (2.*Re2) 0.));
CL2 = rela 'UZ' pbal - 'UZ' pbal2;
 
*assemblage 
CLTOT  = CL1 et CL2;
RIGTOT = RIG1 et CLTOT;
 
 
***************************************************
*** Déformée centrifuge
***************************************************
 
* force centrifuge
XYZ = (CHAN COMP (COOR 1 voltot) 'UX') ET
      (CHAN COMP (COOR 2 voltot) 'UY') ET
      (CHAN COMP (COOR 3 voltot) 'UZ');
FCE1 = (-1.) * KCE1 * XYZ;
trac (vect FCE1 'FORC' 'ROUG') voltot;
 
*resolution
uce1 = RESO RIGTOT FCE1;
trac (defo uce1 voltot);
 
* precontrainte centrifuge
sig1 = sigm uce1 modtot mattot 'LINE';
ksi1 = ksig sig1 modtot mattot 'FLAM';
 
 
****************************************************
* Calcul des modes propres au repos
****************************************************
si(COMPLET);
  nmod1  = 100 ;
  fproc1 = 50.;
  nmax1 = 50;
* nmax1 = 50 permet d'avoir un petit trou a boucher...
sino;
  nmod1  = 16 ;
  fproc1 = 50.;
  nmax1 = 4;
fins;
si(FLDIAG); mess (diag RIGTOT) (diag MAS1); fins;
 
opti impi iimpi;
TMOD1 = VIBR 'SIMUL' fproc1 nmod1 RIGTOT MAS1  nmax1;
* on utilise txt1=4 a la fin pour limiter le nombre de modes calculé / cycle a 4 
*  et ainsi tester rapidement la stratégie de cycle de Lanczos utilisé par VIBR SIMUL 
opti impi 0;
 
opti echo 0;
mess '----------- Calcul des modes propres -----------';
 
*... post traitement "a la main"
mess 'mode   frequence (Hz)';
frq1p = prog;
imod1 = 0;
repe BMOD1 nmod1;  imod1 = imod1 + 1;
  frq1 = TMOD1 . 'MODES' . imod1 . 'FREQUENCE';
  umod1 = TMOD1 . 'MODES' . imod1 . 'DEFORMEE_MODALE';
  mess (chai imod1 '  ' frq1);
  toto1 = chai 'mode propre ' imod1 '  Freq=' frq1 'hz' ;
*   trac (defo umod1 voltot) 'TITR' toto1 'CACH' ;
  frq1p = frq1p et frq1;
fin  BMOD1;
 
* -> on remplace par
*... post traitement "via EXPLORER"
EXPLORER TMOD1;
opti echo 1;
 
 
****************************************************
* Calcul des modes propres en rotation
****************************************************
nmod2  = 4 ;
fproc2 = 50.;
* on se place apres la 1ere critique (>235Hz)
* de maniere a avoir des freq propres negatives
OMEG1 = 250.;
OMEG2 = (2.*pi*OMEG1)**2;
RIG2 = RIGTOT et (OMEG2 * (KCE1 et ksi1));
si(FLDIAG); mess (diag RIG2) (diag MAS1); fins;
 
opti impi iimpi;
TMOD2 = VIBR 'SIMUL' fproc2 nmod2 RIG2 MAS1 ;
* on ajoute VRAI pour avoir les modes dont la valeur propre est negative
opti impi 0;
 
opti echo 0;
mess '----------- Calcul des modes propres en rotation -----------';
mess 'mode   frequence (Hz)';
frq2p = prog;
imod2 = 0;
repe BMOD2 nmod2;  imod2 = imod2 + 1;
  frq2 = TMOD2 . 'MODES' . imod2 . 'FREQUENCE';
  umod2 = TMOD2 . 'MODES' . imod2 . 'DEFORMEE_MODALE';
  mess (chai imod2 '  ' frq2);
  toto2 = chai 'mode propre en rotation ( a ' OMEG1 'Hz)' imod2 
               '  Freq=' frq2 'hz' ;
  trac (defo umod2 voltot) 'TITR' toto2 'CACH' ;
  frq2p = frq2p et frq2;
fin  BMOD2;
opti echo 1;
 
****************************************************
* Calcul des modes critiques
****************************************************
nmod3  = 4 ;
fproc3 = 50.;
MAS3 =  -1.*(KCE1 et ksi1);
si(FLDIAG); mess (diag RIGTOT) (diag MAS3); fins;
* Attention !!! KCE1 et ksi1 matrices de type 'MASSE' (pour VIBR) grace a 'FLAM'
 
opti impi iimpi;
TMOD3 = VIBR 'SIMUL' fproc3 nmod3 RIGTOT MAS3 ;
* on ajoute VRAI pour avoir les modes dont la valeur propre est negative
opti impi 0;
 
opti echo 0;
mess '----------- Calcul des modes critiques -----------';
mess 'mode   frequence (Hz)';
frq3p = prog;
imod3 = 0;
repe BMOD3 nmod3;  imod3 = imod3 + 1;
  frq3 = TMOD3 . 'MODES' . imod3 . 'FREQUENCE';
  umod3 = TMOD3 . 'MODES' . imod3 . 'DEFORMEE_MODALE';
  mess (chai imod3 '  ' frq3);
  toto3 = chai 'mode critique ' imod3 '  Freq=' frq3 'hz' ;
  trac (defo umod3 voltot) 'TITR' toto3 'CACH' ;
  frq3p = frq3p et frq3;
fin  BMOD3;
opti echo 1;
 
 
* opti donn 5 trac x;
 
****************************************************
* Diagramme de Campbell dans le repère tournant
****************************************************
*
* a faire + tard pour verifier CORI ...
 
 
 
******************************************************************
*
*               TNR    
*
*****************************************************************
 
* valeur de ref calculées au 09/01/2012 en local (pclinux64)
*
xtol = 1.E-3 * 50.;
frq1 = 218.46804975;
frq2 = -81.9993;
frq3 =235.55230677;
frq3b =-292.46499848;
error1 = non (exis frq1p frq1 xtol);
error2 = non (exis frq2p frq2 xtol);
error3 = non (exis frq3p frq3 xtol);
error3b = non (exis frq3p frq3b xtol);
error = error1 ou error2 ou error3 ou error3b;
si(complet);
  fr1b = 34361.6036733;
  error1b = non (exis frq1p fr1b xtol);
  error = error ou error1b;
fins;
* Message
*       
TEMP 'IMPR' ;
SI (error) ;
  mess (chai error1  ' ' error2  ' ' error3 ' ' error3b);
  ERRE 5 ;
SINON ;
  MESS 'Aucun probleme detecte !' ;
FINSI ;
*
FIN ;