Télécharger hbm_duffing.dgibi

Retour à la liste

* fichier :  hbm_duffing.dgibi
******************************************************************      
*
*   Systeme a un degre de liberte
*
*   Oscillateur de Duffing calcule en k harmoniques :
*
*   m*x'' + c*x' + k*x + a*x^3 = f*cos(w*t)
*
*   Continuation + HBM + AFT
*
******************************************************************
 
*-----------------------------------------------------------------------
* CHARMECA.PROCEDUR 
* Calcul de la force non lineaire fnl(x,t) et de sa derivee dfnl(x,t)/dx
* dans le domaine temporel
*-----------------------------------------------------------------------
*
*  EN ENTREE : 
*
*  TAB1 
*  >> parametres génériques :
*  . 'DEP_NL'    : LISTCHPO DE DEPLACEMENTS TEMPORELLES
*                  ou (pas encore programmé)
*        . i . 'CHPOINT'  = CHPOINT unitaire des deplacements NL
*        . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
*  >> parametres propres a ce cas test :
*  . 'COEFF_FNL' : coefficient de non linearite
*  . 'POINT_FNL' : point de non linearite
*  . 'COMP_FNL'  : composante de non linearite
*
*
*  EN SORTIE :
*
*  TAB1
*  >> parametres génériques :
*   . 'FNL_T'  :   LISTCHPO DES FORCES NON LINEAIRES TEMPORELLES
*                  ou TABLE d'indice : (pas encore programmé)
*        . i . 'CHPOINT'  = CHPOINT unitaire des forces non lineaires
*        . i . 'LISTREEL' = evolution temporelle coefficient du chpoint
*   . 'KNL_T' :   DERIVEE DE FORCES NON LINEAIRES sous forme de TABLE :
*        . i . j = 'LISTREEL' evolution temporelle de dFi/dxj 
*
*-----------------------------------------------------------------------
DEBPROC CHARMECA TAB1*'TABLE';
 
********     parametres d'entree propres a ce cas test   **********
 
*coefficient de la force nl
SI (exis TAB1 . 'COEFF_FNL');
  AL = TAB1 . 'COEFF_FNL';
SINON;
  AL = 0.;
  MESS 'CHARMECA: PAS D EFFET NON LINEAIRE !';
FINSI;
 
SI (EXIS TAB1 . 'POINT_FNL');
  pNL = TAB1 . 'POINT_FNL';
SINON;
  MESS 'CHARMECA: IL MANQUE LA GEOMETRIE NL dans POINT_FNL'; ERRE 21;
FINSI;
 
SI (EXIS TAB1 . 'COMP_FNL');
  COMP_FNL = TAB1 . 'COMP_FNL';  
  n_CNL = DIME COMP_FNL;
SINON;
  MESS 'CHARMECA: IL MANQUE LA COMPOSANTE NL dans COMP_FNL'; ERRE 21;
FINSI;
 
********     parametres d'entree génériques   **********
 
*DEPNL : list de champ par points, deplacements temporels du point nl
*
SI (EXIS TAB1 'DEP_NL');  
  DEPNL = TAB1 . 'DEP_NL';
  IFCHPO = ega (type DEPNL) 'LISTCHPO';
  si (non IFCHPO);
    si (neg  (type DEPNL) 'TABLE');
      mess 'CHARMECA: DEP_NL doit etre de type LISTCHPO ou TABLE'; 
      ERRE 21;
    sinon; 
      si (neg (type DEPNL . 1) 'LISTREEL');
        mess 'CHARMECA: DEP_NL . i doit etre un LISTREEL'; 
        ERRE 21;
      finsi;
    finsi;
  finsi;
SINON; 
  MESS 'CHARMECA: IL MANQUE LA LISTE DES DEPLACEMENTS dans DEP_NL'; 
  ERRE 21;
FINSI;
 
* calcul de fnl ET dfnl/dx ? dans le doute on calcule tout...
SI (NEG (TYPE IFFNL) 'LOGIQUE'); IFFNL = VRAI; FINSI;
SI (NEG (TYPE IFKNL) 'LOGIQUE'); IFKNL = VRAI; FINSI; 
 
 
 
********  Calcul de FNL  ********
 
*---------------- Cas DEPNL LISTCHPO
SI IFCHPO;
 
* preparation des donnees de sortie
  FNLTEMP = VIDE 'LISTCHPO';
  FNLm = prog; 
  KNLp = prog;
  CNL = EXTR COMP_FNL 1;
 
* nbt  : nombre de pas de temps
  nbt = DIME DEPNL;
  ib = 0;
  REPE BFOR nbt; ib = ib + 1; 
 
    DEP_ib = extr DEPNL ib;
    DEP_i = DEP_ib ; 
    FNL_i = -1.*al* (DEP_i **3) ;
    FNLTEMP = FNLTEMP et FNL_i; 
    FNLm = FNLm et (MAXI FNL_i 'ABS');
 
*   pour le Calcul de DFNL :
    SI (IFKNL);
      xDEP_i = extr DEP_i pNL CNL;
      KNLp = KNLp et (3.*al* (xDEP_i **2));
    FINSI;
 
  FIN BFOR;
  FNLmax = maxi (ABS FNLm);  
 
*---------------- Cas DEPNL TABLE de LISTREEL
SINON; 
 
* DEPNL : table de listreel, indicée par les ddls NL
 
* nbt  : nombre de pas de temps
  nbt = DIME DEPNL . 1; 
  ib = 0;
  FNLTEMP = TABLE; FNLTEMP . 1 = PROG;
  KNLp = prog;
  REPE bfor nbt; ib = ib + 1; 
    DEP_i0 = EXTR DEPNL . 1 ib;
    DEP_i = DEP_i0 ; 
    FNLTEMP . 1 = FNLTEMP . 1 et (-1.*al* (DEP_i **3)); 
    SI (IFKNL);
      KNLp = KNLp et ( 3.*al* (DEP_i **2));
    FINSI;
  FIN bfor;
  FNLmax = maxi (ABS FNLTEMP . 1);
 
FINSI;
 
 
********  Calcul de DFNL  ********
 
SI (IFKNL);
*  uniquement sous forme de LISTREELs
   KNLT1A = TABL;
   KNLT1A . 1 = tabl;
   KNLT1A . 1 . 1 = KNLp;
FINSI;
 
TAB1 . 'FNL_T'  = FNLTEMP;
TAB1 . 'KNL_T' = KNLT1A;
TAB1 . 'FNLMAX' = FNLmax;
 
FINPROC;
*-----------------------------------------------------------------------
 
 
 
************************************************************************      
*
*  dynamique non lineaire 
*  reponse forcee d'un oscillateur de type Duffing
*  reponse par HBM 
*               
************************************************************************
 
*opti echo 0 ;
opti debug 1;
OPTI DIME 2 ELEM SEG2 MODE PLAN ;
 
GRAPH = FAUX ; SAVE = FAUX; TNR = VRAI;
* GRAPH = VRAI ; SAVE = VRAI; TNR = FAUX;
 
 
*-------- Definition des options du calcul :
 
* CALCUL NON LINEAIRE via charmeca?  (DEFAUT : VRAI)
Calc_nl = VRAI;
* si calcul par les listchpo dans AFT et DFT  (DEFAUT : FAUX -> listreel)
IFCHPO = FAUX;
* IFCHPO = vrai;
 
 
*-------- Definition des parametres du calcul :
 
*nombre d'harmoniques utilises
* nhbm = 1;
* nhbm = 3;
nhbm = 7;
* nhbm = 11;
 
*parametre de l'adimensionnement 
u1max = 1.;
 
*des iterations
itmoy = 6;
*ds = 0.1;
ds = 1.;
dsmax = 1.;
dsmin = 1.E-3 ;
 
* si AFT, OTFR tel que le nb de points pour TFR = 2**OTFR 
* OTFR = 7;
OTFR = ENTI 'SUPERIEUR' (log (2 * (2 * nhbm + 1)) / (log 2)); 
mess 'OTFR = ' OTFR;
 
 
*-------- Definition des donnees d'oscillateur ---------------
 
* M d2u/dt2 + C du/dt + K u + pfnl*u^3 = F(t)
*mass
pmass = 1.;
*amortissement
pamor = 0.1 ;
*rigidite
prigi = 1.;
*force d'exterieur
pfext = 0.5;
*coefficient de la force non lineaire
pfnl = 0.25;  
*pfnl = 1.;  
* pfnl = 5.;  
fprec = 1.E-3;
si (ega pfnl 0.25 fprec) ; mofnl = mot '025'; finsi;
si (ega pfnl 1.00 fprec) ; mofnl = mot '1'; finsi;
si (ega pfnl 5.00 fprec) ; mofnl = mot '5'; finsi;
 
* prefix :
prefix = chai 'hbm_duffing_' mofnl '_n' nhbm ;
si GRAPH ;
  ficps = chai prefix '.ps';
  opti TRAC PSC EPTR 6 POTR 'HELVETICA_16' 'FTRA' ficps;
finsi;
 
*plage d'etude de frequence en rad/s
* par defaut, CONTINU va supposer que la FREQUENCE[rad/s] vaut w(t)=t
t0 = 0.1;
si (pfnl <eg 1.00); t1 = 2.45; 
sinon;              t1 = 4.0; 
finsi;
LIS1T = prog t0 PAS 0.05 t1; 
 
 
*------------- geometrie  ---------------    
*
P1 = 0. 0.  ; 
P2 = 1. 0.  ;
ST = P1 D 1  P2 ;                                                                                                    
 
*----------- matrices caracteristiques ----------       
*                                                                               
LIM1 = MOTS UY ;
MASS1 = MASS 'UY' pmass P2 ;
RI1 = MANU 'RIGI' P2 LIM1 (prog prigi); 
AMOR1 = MANU 'RIGI' P2 LIM1 (prog pamor); 
 
 
*------- on définit une unique table pour tout (HBM et CONTINU) -------
TAB1 = TABLE ;
 
*--------- remplissage des matrices "temporelles" ---------
 
* ATTENTION AUX UNITES : ON VA TRAVAILLER EN rad/s ==> pas de 2pi
TAB1 . 'MASSE_CONSTANTE'        = MASS1 ;
TAB1 . 'AMORTISSEMENT_CONSTANT' = AMOR1 ;
TAB1 . 'RIGIDITE_CONSTANTE'     = RI1 ; 
TAB1 . 'N_HARMONIQUE'           = nhbm;
 
*------- remplissage des resultats attendus sur ddls "temporels" -------
mycoul12 = mots VIOL BLEU BLEU TURQ TURQ OCEA OCEA VERT VERT OLIV OLIV 
JAUN JAUN ORAN ORAN ROUG ROUG ROSE ROSE AZUR AZUR GRIS GRIS DEFA DEFA ;
 
TAB1 . 'RESULTATS' = tabl; 
TAB1 . 'RESULTATS' . 1 = tabl; 
TAB1 . 'RESULTATS' . 1 . 'POINT_MESURE' = P2; 
TAB1 . 'RESULTATS' . 1 . 'COMPOSANTE'  = mot 'UY'; 
TAB1 . 'RESULTATS' . 1 . 'TITRE'       = mot 'U'; 
TAB1 . 'RESULTATS' . 1 . 'COULEUR'     = mot 'BLEU'; 
TAB1 . 'RESULTATS' . 1 . 'COULEUR_HBM' = mycoul12; 
 
 
*--------- passage du probleme en frequentiel ---------
 
HBM TAB1;
 
MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP TEMPORELLES';
LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT';
MESS '>>>>>>>>>>>>>>> COMPOSANTES DEP FREQUENTIELLES';
LIST TAB1 .'COMPOSANTES' . 'DEPLACEMENT_HBM'; 
 
MESS 'RIGIDITE_HBM';            LIST TAB1 . 'RIGIDITE_HBM'       ; 
MESS 'BLOCAGES_HBM';            LIST TAB1 . 'BLOCAGES_HBM'       ; 
MESS 'AMORTISSEMENT_HBM';       LIST TAB1 . 'AMORTISSEMENT_HBM'  ; 
MESS 'MASSE_HBM';               LIST TAB1 . 'MASSE_HBM'          ; 
MESS 'MASSE_HBM_1';             LIST TAB1 . 'MASSE_HBM_1'        ; 
MESS 'MASSE_HBM_2';             LIST TAB1 . 'MASSE_HBM_2'        ; 
MESS 'AMORTISSEMENT_HBM_1';     LIST TAB1 . 'AMORTISSEMENT_HBM_1'; 
 
 
*----------- definition du chargement ----------------------------------
 
* ici, on definit directement le chargement en frequentiel
* F2 : composante en cos 1wt
FP11 = MANU 'CHPO' p2 1 'F2' pfext 'NATURE' 'DISCRET' ;
*                                                         
LIX1 = PROG 0. pas 0.1 (t1+10.) ;  
LIY1 = PROG (dime LIX1)*1.;   
EV1 = EVOL MANU 'w' LIX1 'F(w)' LIY1 ; 
CHA1 = CHAR MECA FP11 EV1 ; 
 
 
*----------- conditions initiales -----------------------  
 
* U0 = MANU CHPO p2 3 'U1' 0. 'U2' 0. 'V2' 0. ;
* en equilibre statique lineaire :
K0tot = TAB1 . 'RIGIDITE_HBM' et TAB1 . 'BLOCAGES_HBM';
U0 = RESO K0tot FP11;
*-> estimation lineaire l'equilibre n'est donc pas respecté !
 
 
*----------- OPTIONS DE CALCUL ----------------------
 
TAB1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
TAB1 . 'GEOMETRIE'  = P2;
TAB1 . 'CHARGEMENT' = CHA1;
TAB1 . 'PRECISION' = 1.E-5;
TAB1 . 'TEMPS_CALCULES' = LIS1T;
TAB1 . 'MAXIPAS' = 1000;
TAB1 . 'ACCELERATION' = 4;
TAB1 . 'DEPLACEMENTS' = TABLE;
TAB1 . 'DEPLACEMENTS' . 0 = U0 ;
TAB1 . 'NB_ITERATION' = itmoy;
TAB1 . 'MAXITERATION' = 30;
TAB1 . 'WTABLE' = tabl;
TAB1 . 'WTABLE' . 'DS' = ds;
TAB1 . 'WTABLE' . 'DSMAX' = dsmax;
TAB1 . 'WTABLE' . 'DSMIN' = dsmin;
TAB1 . 'MAXI_DEPLACEMENT' = u1max;
TAB1 . 'PAS_SAUVES' = 2;
 
* calcul des forces NL
TAB1 . 'PROCEDURE_CHARMECA'= Calc_nl ;
TAB1 . 'CALC_CHPO' = IFCHPO;
TAB1 . 'PROCEDURE_FREQUENCE_TEMPS' = mot 'AFT';
TAB1 . 'N_PT_TFR' = OTFR;
 
* parametres du cas test lié a AFT et CHARMECA
TAB1 . 'COEFF_FNL' = pfnl;
TAB1 . 'POINT_FNL' = P2;
TAB1 . 'COMP_FNL'  = mots 'UY' ;
 
* stabilité 
TAB1 . 'STABILITE' = MOTS 'DIAG' 'HILL';
 
 
*----------- resolution par la procedure CONTINU -----------------------  
 
*calcul avec continuation
CONTINU TAB1;
 
TEMP 'IMPR' 'MAXI';
 
*------------------------- post-traitement ----------------------------- 
 
 
********** tracé de toutes les harmoniques individuellement **********
 
wprog = TAB1  . 'TEMPS_PROG';
evtot = TAB1 . 'RESULTATS_HBM' . 'RESULTATS_EVOL';
 
TDESS1 = TABL;
TDESS1 . 'TITRE' = TAB1 . 'RESULTATS_HBM' . 'TITRE';
 
si GRAPH;
dess evtot 
     TITX '\w/\w_{0}' POSX CENT  
     TITY 'U_{k}' POSY CENT TDESS1 LEGE NO;
finsi;
 
*******************   uy2p  : amplitude en norme 2  ******************
 
ihbm = 0;
uy2p =  TAB1 . 'RESULTATS_HBM' . 1 . 'RESULTATS' ** 2;
repe BUYHBM nhbm; ihbm = ihbm + 1;
  ikcos = 2*ihbm;
  iksin = 2*ihbm + 1;
  uy2p = uy2p 
  + (0.5* (  (TAB1 . 'RESULTATS_HBM' . ikcos . 'RESULTATS' ** 2)
           + (TAB1 . 'RESULTATS_HBM' . iksin . 'RESULTATS' ** 2)) );
fin  BUYHBM;
 
*** evolution + tracé *** 
evuy2   = evol VIOL MANU wprog uy2p;
si GRAPH;
dess evuy2 
     TITX '\w/\w_{0}' POSX CENT  
     TITY '|U|_{2}' POSY CENT ;
finsi;
 
 
********** traduction des resultats frequentiels en temporels **********
 
HBM_POST TAB1;
 
************* tracé du max d'amplitude max_t(u(t)) *********************
 
evuyamp = TAB1 . 'RESULTATS' . 'RESULTATS_EVOL';
 
si GRAPH;
dess evuyamp 
    TITX '\w/\w_{0}'    POSX CENT  
    TITY 'max|U(t)|' POSY CENT ;
finsi;
 
************* tracé du max d'amplitude avec la stabilite ! *************
 
istab = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'STABILITE';
Tdess2 = TABL;
Tdess2 . 'LIGNE_VARIABLE' = TABL;
Tdess2 . 'LIGNE_VARIABLE' . 1 = istab;
si GRAPH;
dess evuyamp 
    TITX '\w/\w_{0}'    POSX CENT  
    TITY 'max|U(t)|' POSY CENT Tdess2;
finsi;
* zoom    
Tdess2 . 1 = mot 'MARQ S PLUS';
 
si GRAPH;
dess evuyamp 
    TITX '\w/\w_{0}'    POSX CENT  XBOR 1.55 1.70 YBOR 2.8 3.1
    TITY 'max|U(t)|' POSY CENT Tdess2;
dess evuyamp 
    TITX '\w/\w_{0}'    POSX CENT  XBOR 1.28 1.31 YBOR 0.8 1.5
    TITY 'max|U(t)|' POSY CENT Tdess2;
finsi;
 
* autres tracé de la stabilite : lambda_R en fonction de W
mycoul = mots 'VIOL' 'BLEU' 'AZUR' 'TURQ' 'VERT' 'OLIV' 
              'JAUN' 'ORAN' 'ROUG' 'ROSE';
ncoul = dime mycoul;
wprog2 = enle wprog (dime wprog);
evlrtot = VIDE 'EVOLUTIO' ;
evlitot = VIDE 'EVOLUTIO' ;
nl = dime TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL';
il = 0;
icoul = 0;
repe bl nl; il = il + 1; icoul = icoul + 1;
  si(icoul > ncoul); icoul = 1; finsi; coul1 = extr mycoul icoul;
  lr1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_REEL' . il;
  evlrtot = evlrtot  
  et (evol coul1 MANU '\w/\w_{0}' wprog2 '\l_{R}' lr1);
  li1 = TAB1 . 'RESULTATS_STABILITE' . 'FLOQ' . 'EXPOSANT_IMAG' . il;
  evlitot = evlitot  
  et (evol coul1 MANU '\w/\w_{0}' wprog2 '\l_{I}' li1);
fin  bl;
wprog3 = prog (mini wprog2) (maxi wprog2);
evzero = evol manu '\w/\w_{0}' wprog3 '\l_{R}' (0.*wprog3);
Tdess3 = tabl;
Tdess3 . 1 = mot 'TIRR';
 
si GRAPH;
dess (evzero et evlrtot) 
'TITY' '\l_{R} (s^{-1})' Tdess3;
dess (evlitot) 'TITY' '\l_{I} (Hz)' ;
finsi;
 
 
 
*-------------------------- sauvegarde -----------------------------
 
si SAVE ;
  ficxdr = chai prefix '.xdr'; 
  mess 'sauvegarde dans ' ficxdr;
  opti sauv ficxdr;
  sauv;
finsi;
 
 
 
*--------------------- test de non regression ---------------------
 
si TNR ;
 
* recup valeurs de la courbes de reponse
  wamp  = extr evuyamp 'ABSC';
  uyamp = extr evuyamp 'ORDO';
  namp = DIME uyamp;
 
* on teste la presence de 2 points limites
  dw1    = ( wamp enle 1) - ( wamp enle namp);
  dw1dw0 = (dw1 enle 1) * (dw1 enle (namp - 1));
  mptlim = (lect 0) et (enti 'PROCH' (masq dw1dw0 'EGINFE' 0.)) et 0;
  nptlim = somm mptlim;
  si (nptlim neg 2); 
    mess 'Il devrait y avoir 2 points limites !';
    list nptlim;
    ERRE 5; 
  finsi;
 
* on verifie que pour ces 2 points il y a chgt de stabilité
  iptlim = posi 1 DANS mptlim 'TOUS' ;
  sumstab = (extr istab (iptlim - 1)) + (extr istab (iptlim));
*  sumstab doit valoir 1 (ni 0 ni 2)
  si ( ((mini sumstab) neg 1) ou ((maxi sumstab) neg 1) ); 
    mess 'Les 2 points limites doivent correspondre a un changement '
         'de stabilite!';
    list sumstab;
    ERRE 5; 
  finsi;
 
* on termine par verifier l'amplitude en ces 2 points + l'amplitude max
  Utest = (extr uyamp iptlim) et (maxi uyamp) ;
* on compare a une precedente execution (lors de la creation du cas test)
  Uref = prog 3.0705       1.0831       3.0874;
  si ( (MAXI (Utest - Uref) 'ABS') > 0.01 ); 
    mess 'Erreur dans les amplitudes predites !';
    list Utest;
    ERRE 5; 
  finsi;
 
finsi;
 
 
fin ;
 
 
 
 
 
 
 
 
 
 

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