* 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 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 ;