* fichier : four3.dgibi
************************************************************************
*
* Dynamique, calcul fréquentiel, reponse harmonique, modes de Fourier
*
* Etude de la flexion d'un cylindre sollicité par une excitation
* harmonique en mode de Fourier
* test des operateurs : IMPE, VIBR, RESO, PJBA
* creation : 2013-11-29, Benoit Prabel (#7870)
*
************************************************************************
*** OPTIONS ************************************************************
COMPLET = FAUX ;
GRAPH = VRAI; OPTI 'TRAC' 'PSC' 'POTR' 'HELVETICA_16' EPTR 5;
*GRAPH = VRAI;
*
OPTI DIME 2 ELEM QUA4 MODE FOUR NOHA ;
nhar = 1;
*** DONNEES GEO + DISCRETISATION ***************************************
*
R = 0.01 ; H = 0.5;
nr = 2 ; nz = 20;
moU = mots 'UR' 'UT' 'UZ';
moIU = mots 'IUR' 'IUT' 'IUZ';
*** MAILLAGE ***********************************************************
*
p1 = 0. 0.;
p2 = R 0.;
d12 = p1 droi nr p2;
d34 = d12 plus (0. H);
p3 = poin d34 initial; p4 = d34 poin final;
s1 = d12 regl nz d34;
*** MODELE MATERIAU ****************************************************
*
E1 = 2.E11; nu1=0.3; rho1=8000.;
visq1 = 1.E-1*((E1 * rho1)**0.5);
mod1 = mode s1 mecanique elastique isotrope;
mat1 = mate mod1 YOUN E1 NU nu1 RHO rho1 VISQ visq1;
*** MATRICES ***********************************************************
*
* matrices : partie symetrique (UR UZ UT seulement)
opti mode four nhar;
K1 = rigi (mod1) (mat1);
M1 = mass (mod1) (mat1);
C1 = amor mod1 mat1;
blo1 = BLOQ UR UT UZ d12;
Ktot1 = K1 et blo1;
*matrices impedances version 1 : on duplique la symetrique pour ddl IU
K2 = impe K1 'RAIDEUR';
M2 = IMPE M1 'MASSE' ;
C2 = impe C1 'AMORTISSEMENT';
* ! attention ! cela n'est pas equivalent a :
* IMPE C1 Cm1 'AMORTISSEMENT';
* ou Cm1 calculé depuis OPTI MODE FOUR (-1*nhar);
* car il ne s'agit de machines tournantes ici : on reste
* mono-directionnel : les ddl imaginaires correspondent au dephasage !
*
* blo2 = blo1 et (BLOQ IUR IUT IUZ d12);
* ci dessus avant evol #7774, ci dessous apres :
blo2 = IMPE blo1 'RAIDEUR';
Ktot2 = K2 et blo2;
*** CALCUL MODAL *******************************************************
*
nmod1 = 2;
fpro1 = 1.1;
* tmod1 = VIBR 'PROCH' Ktot1 (M1) (prog fpro1) (lect nmod1)
tmod1 = VIBR 'PROCH' Ktot2 (-1.*M2) (prog fpro1) (lect nmod1)
'TBAS' VRAI;
imod1 = 0;
repe bmod1 nmod1; imod1 = imod1 + 1;
frq1 = tmod1 . 'MODES' . imod1 . 'FREQUENCE' ;
u1 = tmod1 . 'MODES' . imod1 . 'DEFORMEE_MODALE' ;
tit1 = chai 'mode ' imod1 ' f=' frq1 'Hz';
mess tit1;
list (resu u1);
si GRAPH; TRAC (defo u1 s1) TITR tit1; finsi;
fin bmod1;
tstat = tabl 'LIAISONS_STATIQUES';
*** REPONSE **************************************************
* A UNE FORCE UNITAIRE
F2 = MANU 'CHPO' p4 2 'FR' 0.5 'FT' -0.5 NATURE DISCRET;
* Projection sur base modale
M3 = pjba M2 tmod1;
C3 = pjba C2 tmod1;
Ktot3 = pjba Ktot2 tmod1;
F3 = pjba F2 tmod1;
*** REPONSE STATIQUE **************************************************
* A UNE FORCE UNITAIRE sur base physique
u2 = reso Ktot2 F2;
urp4 = extr u2 p4 'UR';
utp4 = extr u2 p4 'UT';
mess urp4 utp4;
* A UN DEPLACEMENT IMPOSÉ
CLurp4 = BLOQ 'UR' p4;
CLutp4 = BLOQ 'UT' p4;
UIMP2 = (DEPI CLurp4 urp4) et (DEPI CLutp4 utp4);
uu2 = RESO (Ktot2 et CLurp4 et CLutp4) UIMP2;
ff2 = reac (CLurp4 et CLutp4) uu2;
* A UNE FORCE UNITAIRE sur base modale
u3 = reso Ktot3 F3;
u3 = reco u3 tmod1 tstat;
ur3p4 = extr u3 p4 'UR';
ut3p4 = extr u3 p4 'UT';
mess ur3p4 ut3p4;
* A UN DEPLACEMENT sur base modale
* UIMP3 = (DEPI CLurp4 ur3p4) et (DEPI CLutp4 ut3p4);
* UIMP3 = PJBA UIMP3 tmod1;
* CLur3p4 = PJBA CLurp4 tmod1;
* CLut3p4 = PJBA CLutp4 tmod1;
* !!! attention !!!
* A ce jour, il est impossible de projeter directement sur base modale :
* - un second membre de type FLX portant sur des relation non decrite
* dans cette base !
* - des relations et blocage non decrite dans cette base !
*==> pour l instant on fait "a la main" :
UIMP3 = VIDE 'CHPOINT'/'DIFFUS';
imod1 = 0;
repe bmod1 nmod1; imod1 = imod1 + 1;
umod1 = Tmod1 . 'MODES' . imod1 . 'DEFORMEE_MODALE';
ptrep = Tmod1 . 'MODES' . imod1 . 'POINT_REPERE';
UIMP3 = UIMP3 et (MANU 'CHPO' ptrep 1 'ALFA' (extr umod1 'UR' p4)
'NATURE' 'DIFFUS');
fin bmod1;
CL3 = RELA UIMP3;
list CL3;
UIMP3 = DEPI CL3 ur3p4;
list UIMP3;
uu3 = RESO (Ktot3 et CL3) UIMP3;
* verif du respect de la liaison :
uu3u = reco uu3 tmod1 tstat;
list (redu uu3u p4);
ERROR6u = maxi (prog ((extr uu3u 'UR' p4) - ur3p4)
((extr uu3u 'UT' p4) - ut3p4)) 'ABS';
ff3 = reac CL3 uu3;
*** REPONSE DYNAMIQUE **************************************************
* on trace les deformees pour les frequences suivantes
prtrac = prog 25. 50. 75. 100. PAS 100. 500.;
* fonction de transfert : on boucle sur les frequences
si (COMPLET);
promeg = prog 50. PAS 2. 70.
71. PAS 0.2 74.8
74.9 PAS 0.05 75.7
75.8 PAS 0.2 80.
82. PAS 2. 100.;
sino;
promeg = prog 70. PAS 0.2 80.;
fins;
* on cherche la resonance exacte (rem : on a de l amortissement)
promeg = ordo (promeg et frq1);
nomeg = dime promeg ;
* preparation des listreel resultats
ur2top = prog; ut2top = prog;
iur2top = prog; iut2top = prog;
ur3top = prog; ut3top = prog;
iur3top = prog; iut3top = prog;
iomeg= 0 ;
REPE BOMEG nomeg; iomeg = iomeg + 1 ;
omeg = extr promeg iomeg; w1 = 2. * pi * omeg;
chacha = chai iomeg ' : \W = ' omeg 'Hz';
mess chacha;
* formation de la rigidite dynamique
Kdyn2 = Ktot2 et (w1 * C2) et ((w1 ** 2) * M2);
* resolution
u2 = reso Kdyn2 F2;
* resultats
ur2top = ur2top et (extr u2 UR p4);
ut2top = ut2top et (extr u2 UT p4);
iur2top = iur2top et (extr u2 IUR p4);
iut2top = iut2top et (extr u2 IUT p4);
* trace
si (GRAPH);
si (exis prtrac omeg 0.0001);
trac ((defo u2 s1 BLEU) et (defo (u2 exco moIU moU) s1 OCEA))
'TITR' chacha;
fins;
fins;
* idem sur base modale
Kdyn3 = Ktot3 et (w1 * C3) et ((w1 ** 2) * M3);
u3 = reso Kdyn3 F3;
u3 = reco u3 tmod1 tstat;
ur3top = ur3top et (extr u3 'UR' p4);
ut3top = ut3top et (extr u3 'UT' p4);
iur3top = iur3top et (extr u3 'IUR' p4);
iut3top = iut3top et (extr u3 'IUT' p4);
FIN BOMEG;
*** POST TRAITEMENT ****************************************************
* evolutions brutes
ev2r = EVOL 'BLEU' 'MANU' 'OMEG' promeg 'UR' ur2top;
ev2t = EVOL 'BLEU' 'MANU' 'OMEG' promeg 'UT' ut2top;
ev2ir = EVOL 'OCEA' 'MANU' 'OMEG' promeg 'IUR' iur2top;
ev2it = EVOL 'OCEA' 'MANU' 'OMEG' promeg 'IUT' iut2top;
ev3r = EVOL 'ROUG' 'MANU' 'OMEG' promeg 'UR' ur3top;
ev3t = EVOL 'ROUG' 'MANU' 'OMEG' promeg 'UT' ut3top;
ev3ir = EVOL 'ROSE' 'MANU' 'OMEG' promeg 'IUR' iur3top;
ev3it = EVOL 'ROSE' 'MANU' 'OMEG' promeg 'IUT' iut3top;
* amplitude phase
ampx2 = ((ur2top**2) + (iur2top**2))**0.5;
ev2ax = EVOL 'BLEU' 'MANU' 'OMEG' promeg '|UX|' ampx2;
phix2 = atg iur2top ur2top;
ev2px = EVOL 'BLEU' 'MANU' 'OMEG' promeg 'phi' phix2;
ampx3 = ((ur3top**2) + (iur3top**2))**0.5;
ev3ax = EVOL 'BLEU' 'MANU' 'OMEG' promeg '|UX|' ampx3;
phix3 = atg iur3top ur3top;
ev3px = EVOL 'BLEU' 'MANU' 'OMEG' promeg 'phi' phix3;
* tracé des résultats
si GRAPH;
* table Tdess1
Tdess1 = tabl;
Tdess1 . 1 = mot ' MARQ PLUS REGU';
Tdess1 . 2 = mot ' MARQ S ROND REGU';
Tdess1 . 3 = mot 'TIRR MARQ PLUS REGU';
Tdess1 . 4 = mot 'TIRR MARQ S ROND REGU';
Tdess1 . 'TITRE' = tabl;
Tdess1 . 'TITRE' . 1 = mot 'UR';
Tdess1 . 'TITRE' . 2 = mot 'UT';
Tdess1 . 'TITRE' . 3 = mot 'IUR';
Tdess1 . 'TITRE' . 4 = mot 'IUT';
dess (ev2r et ev2t et ev2ir et ev2it)
TITX '\W (Hz)' POSX 'CENT'
TITY 'U (m)' POSY 'CENT'
TITR 'Fonction de transfert poutre Fourier base physique'
LEGE 'SO' Tdess1 ;
dess (ev3r et ev3t et ev3ir et ev3it)
TITX '\W (Hz)' POSX 'CENT'
TITY 'U (m)' POSY 'CENT'
TITR 'Fonction de transfert poutre Fourier base modale'
LEGE 'SO' Tdess1 ;
Tdess2 = tabl;
Tdess2 . 1 = mot ' MARQ PLUS REGU';
Tdess2 . 2 = mot 'TIRR MARQ S ROND REGU';
Tdess2 . 'TITRE' = tabl;
Tdess2 . 'TITRE' . 1 = mot 'base physique';
Tdess2 . 'TITRE' . 2 = mot 'base modale';
dess (ev2ax et ev3ax)
GRIL POINTILLE GRIS
TITX '\W (Hz)' POSX 'CENT'
TITY '|U| (m)' POSY 'CENT' LOGY
TITR 'Fonction de transfert poutre Fourier'
LEGE 'NO' Tdess2 ;
dess (ev2px et ev3px)
GRIL POINTILLE GRIS
TITX '\W (Hz)' POSX 'CENT' TITY '\j (°)' POSY 'CENT' YGRA 45.
TITR 'Fonction de transfert poutre Fourier'
LEGE 'SO' Tdess2 ;
finsi;
*** TESTS DE BON FONCTIONNEMENT ****************************************
* valeur de reference (obtenue depuis un calcul 3d avec des elements cub8
* et 24 elements sur la circonference)
wref3d = 75.354;
uref3d = 1.51858E-03;
* tolerance
XTOL1 = 1.E-2;
XTOL3 = 5.E-2;
XTOL4 = 0.5E-2;
XTOL5 = 1.E-9;
* frequence du 1er mode
ERROR1 = abs (1. - (frq1 / wref3d));
mess frq1 ERROR1;
ERROR1 = ERROR1 > XTOL1;
* frequence de resonance calculée depuis la phase = -90°
wres2 = ipol -90. phix2 promeg;
ERROR2 = abs (1. - (wres2 / wref3d));
mess wres2 ERROR2;
ERROR2 = ERROR2 > XTOL1;
* amplitude de la resonance
uamp2 = ipol frq1 promeg ampx2;
ERROR3 = abs (1. - (uamp2 / uref3d));
mess uamp2 ERROR3;
ERROR3 = ERROR3 > XTOL3;
* verif de la reponse projetee sur base modale
dur = (maxi(ur2top - ur3top) abs) / (maxi ur2top abs);
dut = (maxi(ut2top - ut3top) abs) / (maxi ut2top abs);
diur = (maxi(iur2top - iur3top) abs) / (maxi iur2top abs);
diut = (maxi(iut2top - iut3top) abs) / (maxi iut2top abs);
ERROR4 = maxi (prog dur dut diur diut);
mess dur dut diur diut ERROR4 XTOL4;
ERROR4 = ERROR4 > XTOL4;
* verif du deplacement impose sur base physique
ERROR5 = maxi (ff2 - F2) ABS;
ERROR5 = ERROR5 > XTOL5;
* verif du deplacement impose sur base physique
ERROR6 = maxi (ff3 - F3) ABS;
ERROR6 = (ERROR6 > XTOL5) ou (ERROR6u > XTOL5);
ERROR7 = (maxi (prog (ur3p4 - urp4) (ut3p4 - utp4) ) ABS)
/ (maxi (prog urp4 utp4) ABS);
mess ur3p4 urp4 ut3p4 utp4 ERROR7 XTOL3;
ERROR7 = ERROR7 > XTOL3;
TEMP 'IMPR' 'SGAC' 'IMPR';
si (ERROR1 ou ERROR2 ou ERROR3 ou ERROR4
ou ERROR5 ou ERROR6 ou ERROR7);
mess ERROR1 ERROR2 ERROR3 ERROR4 ERROR5 ERROR6 ERROR7;
ERRE 5;
sinon ;
MESS 'Aucun probleme detecte !' ;
finsi;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer