Télécharger con_calc.procedur

Retour à la liste

* CON_CALC  PROCEDUR  FANDEUR   22/03/01    21:15:03     11301          
*=======================================================
*
* CON_CALC
* ANALyse : calcul de la REPONSE_BALOURD NONlineaire par CONTINUATION
* calcule le Résidu R, la raideur K=dR/dU, les dérivées dR/dt et dK/dt
* creation : BP 28/01/2014
*
*=======================================================
*
************************************************************************
*
* CREATION : BP208322 28/01/2014
*
* OBJET    : Procedure de calcul de :
*            1. du Residu = F^ext - F^int
*            2. de la Raideur K = dResidu / du
*            3. de dt * la dérivée du Residu par rapport a t 
*            4. de dt * la dérivée de la Raideur par rapport a t 
*
* ENTREE   : 
* TAB1 = TABLE
*   . 'MODELE'
*   . 'CARACTERISTIQUES' 
*   . 'RIGIDITE_CONSTANTE'
*   . 'BLOCAGES_MECANIQUES'
*   . 'CHARGEMENT'
*   . 'TEMPS_CALCULES'        = plage indicative des temps a calculer et 
*                               de la finesse des pas
*   . 'GRANDS_DEPLACEMENTS'
*   . 'K_SIGMA'
*   . 'MAXITERATION'
*   . 'PROCEDURE_CHARMECA'    = vrai si chargement dependant de la confi
*                               -guration (pressions suiveuses par ex.)
*   . 'WTABLE'                = table avec infos de calcul interne
*       . 'CONTRAINTES'       = contraintes de Cauchy a t_n+1
* LMOT1 = LISTMOTS de mots-clés définissant les actions a effectuer
*         a choisir parmi { RESI(DU), RIGI, DRES(IDU), DRIG(I), ...}
*
* SORTIE   : 
* TAB1 . 'WTABLE'
*   . 'RESIDU'          = Residu = F^ext - F^int
*   . 'RAIDEUR'         = K = dResidu/du
*   . 'DERIVEE_RESIDU'  = dResidu/dt  * dt
*   . 'DERIVEE_RAIDEUR' = dK/dt * dt
*
* REM : 
* on suppose etre sur la config finale n+1
*
************************************************************************
 
DEBPROC CON_CALC  TAB1*'TABLE'    LMOT1*'LISTMOTS' 
                  time*'FLOTTANT' dtime*'FLOTTANT';   
 
 
************************************************************************
*                                                                      *
*         VERIFICATION ou simple recup DES DONNEES D'ENTREE            *
*                                                                      *
************************************************************************
* fldebug = (vale debug) >eg 1;
fldebug = faux;
 
SI (EXIS TAB1 'WTABLE');  
  WTAB = TAB1 . 'WTABLE';
SINON; 
  MESS 'IL MANQUE LA WTABLE'; ERRE 21;
FINSI;
 
FLCHAR1 = exis TAB1 'CHARGEMENT';
si (FLCHAR1); 
  CHAR1 = TAB1 . 'CHARGEMENT'; 
finsi;
 
FLMOD1  = exis TAB1  'MODELE'; 
si (FLMOD1);  
  MOD1  = TAB1 . 'MODELE';
  MAT1  = TAB1 . 'CARACTERISTIQUES';
  SIG1  = WTAB . 'CONTRAINTES' ;
finsi;
DEP1  = WTAB . 'DEPLACEMENTS';
 
* prendre les matrices constantes ou HBM ?
FLHBM = WTAB . 'HBM';
* Alternating Frequency / Time ?
IAFT =  WTAB . 'AFT';
* si AFT, Jacobienne calculee par TFR(Knl*gamma+Cnl*w*gamma°)
* ou par differentiation numerique ?
IDIFF = WTAB . 'DIFFERENTIATION';
* dans le cas dynamique frequentiel, on a besoin de la frequence w
* est-elle une inconnue ou fixee ?
IFREQ = WTAB . 'FREQ_INCONNUE';
si IFREQ;  w = EXTR DEP1 'FREQ'  WTAB . 'PT_FREQ';
           WTAB . 'FREQ' = w;
sinon;     w = IPOL TAB1 . 'FREQUENCE' time;
finsi;
 
* initialisations
TNL1  = mot 'NONDEFINI';
Fext1 = mot 'NONDEFINI';
TVARI1 = mot 'NONDEFINI';
 
SI FLHBM;
  morigi = mot 'RIGIDITE_HBM';
  momass = mot 'MASSE_HBM';
  moamor = mot 'AMORTISSEMENT_HBM';
  mocori = mot 'CORIOLIS_HBM';
  mokcen = mot 'CENTRIFUGE_HBM';
  mobloc = mot 'BLOCAGES_HBM';
SINON;
  morigi = mot 'RIGIDITE_CONSTANTE';
  momass = mot 'MASSE_CONSTANTE';
  moamor = mot 'AMORTISSEMENT_CONSTANT';
  mobloc = mot 'BLOCAGES_MECANIQUES';
FINSI;
 
KBLO1 = TAB1 . mobloc;
 
 
************************************************************************
*                                                                      *
*                     CALCUL DU RESIDU                                 *
*                                                                      *
************************************************************************
 
*** calcul des efforts Fext + Fnl ***
 
SI ((EXIS LMOT1 'RESI') ou (EXIS LMOT1 'DRES') ou IFREQ);
 
* Forces exterieure dues au chargement
  si (FLCHAR1); Fext1 = TIRE CHAR1 time;
  sinon;        Fext1 = VIDE 'CHPOINT'/'DISCRET';
  finsi;
*   mess '_Fext1='; list resum Fext1;
 
* Autres forces (types suiveuses par ex ou F^NL perso) 
  SI (TAB1 . 'PROCEDURE_CHARMECA');
    SI (EGA IAFT 0); TNL1 = CHARMECA TAB1 time;  FINSI;
    SI (EGA IAFT 1); TNL1 = DFT TAB1 time;       FINSI;
    SI (EGA IAFT 2); TNL1 = AFT TAB1 time LMOT1; FINSI;
    FNL1 = TNL1 . 'ADDI_SECOND';
*     mess 'concalc -> AFT 'time'=> FNL1='; list resum FNL1;
    Fext1 = Fext1 + FNL1;
  FINSI;
 
* Forces non-lineaires des structures variables BRASERO
  SI (EXIS TAB1 'STRUCT_VARIABLE');
*   attention a l'unite de w  !
    TVARI1 = BRASE403 TAB1 DEP1 w LMOT1;
    FVARI1 = TVARI1 . 'ADDI_SECOND';
    Fext1 = Fext1 + FVARI1;
  FINSI;
 
FINSI;
 
*** calcul des effort internes ***
 
SI ((EXIS LMOT1 'RESI') ou (EXIS LMOT1 'DRES'));
 
* Forces internes
* ! on suppose deja etre sur la config deformee !
  SI (FLMOD1);  Fint1 = BSIG SIG1 MOD1 MAT1;
  SINON;        Fint1 = VIDE 'CHPOINT'/'DISCRET';
  FINSI;
*   mess '_Fint1='; list resum Fint1;
 
* rappel : en frequentiel, K^dyn = K + t*C + t^2*M 
  FLdyn1 = faux;
  SI (EXIS TAB1 morigi);
    Kdyn1 = TAB1 . morigi;
*     mess '_on prend K^Harm';
    FLdyn1 = vrai;
  SINON;
    Kdyn1 = VIDE 'RIGIDITE'/'RIGIDITE';
  FINSI;
  SI (EXIS TAB1 moamor);
*   attention a l'unite de w  !
    Kdyn1 = Kdyn1 et (w * TAB1 . moamor);
*     mess '_on prend t * C^Harm';
    FLdyn1 = vrai;
  FINSI;
  SI (EXIS TAB1 momass);
*   attention a l'unite de w  !
    Kdyn1 = Kdyn1 et ((w**2) * TAB1 . momass);
*     mess '_on prend t^2 * M^Harm';
    FLdyn1 = vrai;
  FINSI;
  SI (FLdyn1); 
*     Fdyn1 = Kdyn1 * (DEP1 enle 'LX');
    Fdyn1 = Kdyn1 * DEP1;
*     mess '_Fdyn1='; list resum Fdyn1;
    Fint1 = Fint1 + Fdyn1;
  FINSI;
*dl + prise en compte des reactions
  Freac1 = REAC KBLO1 DEP1;
  WTAB . 'REACTION' = Freac1;
  Fint1 = Fint1 - Freac1;
*dl + :fin de l'ajout
 
* calcul et stockage effectif du Residu  
*bp-test   SI (EXIS LMOT1 'RESI');
    Res1 = Fext1 - Fint1;
*   on enleve les deplacements imposés deja vérifiés
    Res1 = Res1 - ((KBLO1*DEP1) EXCO 'FLX' 'NOID' 'FLX');
*   mess '_Res1='; list resum Res1;
    WTAB . 'RESIDU' = Res1;
*bp-test   FINSI;
 
* on va calculer ici une norme pour tester la convergence  
*     si (ega NF2 'NONDEFINI');
*    xden1 = maxi ((PSCA Res1 Res1 motCf motCf)**0.5);
*    xden2 = maxi ((PSCA Fext1 Fext1 motCf motCf)**0.5);
*    xden3 = maxi ((PSCA Fint1 Fint1 motCf motCf)**0.5);
*    mess 'xdeno = ' xden1 xden2 xden3;
* *   si (NF2 < 1.E-30); mess 'def de NF2 a revoir'; erre 21; finsi;
* *   WTAB . 'XDENO' = NF2;
*     finsi;
 
FINSI;
 
 
************************************************************************
*                                                                      *
*                     CALCUL DE LA RAIDEUR  (= - dResidu/dU)           *
*                                                                      *
************************************************************************
SI (EXIS LMOT1 'RIGI');
 
* Raideur elastique du modele
* ! on suppose deja etre sur la config deformee !
  SI (FLMOD1); Ktot1 = RIGI MOD1 MAT1;
  SINON;       Ktot1 = VIDE 'RIGIDITE'/'RIGIDITE';
  FINSI;
* relations et blocages
  Ktot1 = Ktot1 et KBLO1;
 
* Ksigma  
  SI (TAB1 . 'K_SIGMA');
    Ksig1 = KSIG SIG1 MOD1 MAT1;
    Ktot1 = Ktot1 et Ksig1;
  FINSI;
 
* Autres raideurs (dues au forces suiveuses ou autre K^NL perso) 
  SI (TAB1 . 'PROCEDURE_CHARMECA');
    si (neg (type TNL1) 'TABLE');
      SI (EGA IAFT 0); TNL1 = CHARMECA TAB1 time;  FINSI;
      SI (EGA IAFT 1); TNL1 = DFT TAB1 time;       FINSI;
      SI (EGA IAFT 2); TNL1 = AFT TAB1 time LMOT1; FINSI;
    finsi;
    KNL1 = TNL1 . 'ADDI_MATRICE';
    Ktot1 = Ktot1 et KNL1;
  SINON;
    KNL1 = VIDE 'RIGIDITE'/'RIGIDITE';
  FINSI;
 
* raideurs non-lineaires des structures variables BRASERO
  SI (EXIS TAB1 'STRUCT_VARIABLE');
    si (neg (type TVARI1) 'TABLE');
*   attention a l'unite de w  !
      TVARI1 = BRASE403 TAB1 DEP1 w LMOT1;
    finsi;
    KVARI1 = TVARI1 . 'ADDI_MATRICE';
    Ktot1 = Ktot1 et KVARI1;
    WTAB . 'RAIDEUR_NL' = KVARI1;
  FINSI;
 
* rappel : en frequentiel, K^dyn = K + w*C + w^2*M 
  SI (EXIS TAB1 morigi);
    Ktot1 = Ktot1 et TAB1 . morigi;
  FINSI;
  SI (EXIS TAB1 moamor);
*   attention a l'unite de w  !
    Ktot1 = Ktot1 et (w * TAB1 . moamor);
  FINSI;
  SI (EXIS TAB1 momass);
*   attention a l'unite de w  !
    Ktot1 = Ktot1 et ((w**2) * TAB1 . momass);
  FINSI;
 
* cas ou w est une inconnue indpte : il faut ajouter -dR/dw  a -dR/dU
  SI IFREQ;
    SI (TAB1 . 'PROCEDURE_CHARMECA');
*     on perturbe de 1.E-6  
      dw1c = 1.E-6;
      w1c = w + dw1c;
      WTAB . 'FREQ' = w1c;
*     Fext est considéré indpt de w, mais pas FNL
        SI (EGA IAFT 0); TNL1c = CHARMECA TAB1 time;  FINSI;
*       SI (EGA IAFT 1); TNL1c = DFT TAB1 time; FINSI;
        SI (EGA IAFT 2); TNL1c = AFT TAB1 time; FINSI;
        FNL1c = TNL1c . 'ADDI_SECOND';
*    -dR/dw = ...
      dRdw = (-1./dw1c) * (FNL1c - FNL1) ;
      dRdw = MANU 'RIGI' dRdw 'COLO' 'QUEL' WTAB . 'PT_FREQ' 'FREQ';
      Ktot1 = Ktot1 et dRdw;
*     on remet en l'etat initial    
      WTAB . 'FREQ' = w;
    FINSI;
*   sans oublier la partie constante     
    FLdyn1 = faux;
    SI (EXIS TAB1 moamor);
      dKdyn1 = TAB1 . moamor;
      FLdyn1 = vrai;
    SINON;
      dKdyn1 = VIDE 'RIGIDITE'/'RIGIDITE';
    FINSI;
    SI (EXIS TAB1 momass);
      dKdyn1 = dKdyn1 et ((2. * w) * TAB1 . momass);
      FLdyn1 = vrai;
    FINSI; 
    SI FLdyn1;
      dZUdw = dKdyn1 * DEP1;
      dZUdw = MANU 'RIGI' dZUdw 'COLO' 'QUEL' WTAB . 'PT_FREQ' 'FREQ';
      Ktot1 = Ktot1 et dZUdw;  
    FINSI; 
  FINSI;
 
* stockage effectif de la raideur
  WTAB . 'RAIDEUR'     = Ktot1;
 
FINSI;
 
 
************************************************************************
*                                                                      *
*                 CALCUL DE LA DERIVEE DU RESIDU avec t                *
*                                                                      *
************************************************************************
SI (EXIS LMOT1 'DRES');
 
* REM : quand on ne sait pas faire autrement, on calcule F(t+dt1b)-F(t)
*  dt1b = 0.5*dtime;
*  dt1b = 0.01*dtime;
  si(dtime > time); dt1b = 1.E-2 * dtime;
  sinon ;           dt1b = 1.E-6 * time;
  finsi;
  time1b = time + dt1b;
  si IFREQ;
*   cas ou w est une inconnue indpte: w ne varie pas explictement avec t
    w1b = w;
  sinon;
    w1b = IPOL TAB1 . 'FREQUENCE' time1b; 
    si (neg dt1b 0.); dw1b = (w1b - w) / dt1b; 
    sinon;            dw1b = 0.;
    finsi;
  finsi;
 
* Variation des Forces exterieures dues au chargement
  si (FLCHAR1);   Fext1b = TIRE CHAR1 time1b;
  sinon;          Fext1b = VIDE 'CHPOINT'/'DISCRET';
  finsi;
 
* Variation des Autres forces (types suiveuses par ex ou F^NL perso) 
* dependent-elles explicitement de time?
  SI (TAB1 . 'PROCEDURE_CHARMECA');
    SI (EGA IAFT 0); TNL1b = CHARMECA TAB1 time1b;   FINSI;
*     SI (EGA IAFT 1); TNL1b = DFT TAB1 time1b;        FINSI;
    LMOT1b = MOTS 'DRES';
    SI (EGA IAFT 2); TNL1b = AFT TAB1 time1b LMOT1b; FINSI;
    FNL1b = TNL1b . 'ADDI_SECOND';
*     mess 'concalc -> AFT 'time1b'=> FNL1b='; list resum FNL1b;
    Fext1b = Fext1b + FNL1b;
  FINSI;
 
* Variation des forces non-lineaires des structures variables BRASERO
* rem : calcul explicite pour l instant
*       on pourrait avoir directement les matrices ..? todo + tard
  SI (EXIS TAB1 'STRUCT_VARIABLE');
*  si w=0, avec structure variable amortissante, on n'a pas de 
*  "raideur" -> system irresoluble si second membre non nul  
   si (neg w 0.);
    TVARI1b = BRASE403 TAB1 DEP1 w1b (mots 'RESI');
    FVARI1b = TVARI1b . 'ADDI_SECOND';
    Fext1b = Fext1b + FVARI1b;
   finsi;
  FINSI;
 
* on calcule la variation des efforts en maintenant les u^imp a t
*   DFext1 = ((dtime / dt1b) * ((Fext1b - Fext1) ENLE 'FLX'))
*          + (EXCO Fext1 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET');  
* on calcule la variation des efforts et des u^imp 
  si (ega dtime 0.);
    DFext1 = 0.*(Fext1b - Fext1);
  sino;
    DFext1 = (dtime / dt1b) * (Fext1b - Fext1);
  finsi;
 
* Variation des Forces internes :
* on les suppose ne dependre que de u et pas de t -->  DFint1 = 0
  DFint1 = VIDE 'CHPOINT'/'DISCRET';
 
* rappel : en frequentiel, dK^dyn/dt = dw/dt * [C + 2*w*M]
* et donc dF = dtime * [dK^dyn/dt]*(u enle LX)
  SI (NON IFREQ);
    FLdyn1 = faux;
    SI (EXIS TAB1 moamor);
      dKdyn1 = TAB1 . moamor;
      FLdyn1 = vrai;
    SINON;
      dKdyn1 = VIDE 'RIGIDITE'/'RIGIDITE';
    FINSI;
    SI (EXIS TAB1 momass);
      dKdyn1 = dKdyn1 et ((2. * w) * TAB1 . momass);
      FLdyn1 = vrai;
    FINSI;
    SI ((neg dtime 0.) et (neg dw1b 0.) et FLdyn1); 
      DFdyn1 = dtime * dw1b * (dKdyn1 * DEP1);
      DFint1 = DFint1 + DFdyn1;
    FINSI;
  FINSI;
 
 
* calcul et stockage effectif de dResidu/dt*dt
  DRes1 = DFext1 - DFint1;
  WTAB . 'DERIVEE_RESIDU' = DRes1;
 
FINSI;
 
 
 
 
FINPROC ;
 
 
 
 
 

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