* fichier : vpparf1.dgibi
*
complet = faux;
*
*************************************************
* pour calcul complet mettre complet à : vrai;
*************************************************
**
*
* FICHIER GIBIANE POUR TESTER L'IMPLANTATION
* DU MODELE VISCO PLASTIQUE PARFAIT
*
* On impose une force sur l'extremite d'un barreau cylindrique
* La force evolue dans le temps
* La température varie également et modifie les carcatéristiques
* mécaniques
*
* Le resultat théorique est obtenu par integration
*
* MAILLAGE
*
opti dime 2 mode axis elem qua4 ;
* ajout de option epsilon lineaire pour la precision des test!
OPTION epsilon lineaire;
graph = 'N';
*
p1 = 0. 0. ;
p2 = 5.d-3 0. ;
p3 = 5.d-3 100.d-3 ;
p4 = 0. 100.d-3 ;
*
n1 = 1 ;
n2 = 2 ;
*
li1 = d n1 p1 p2 ;
li2 = d n2 p2 p3 ;
li3 = d n1 p3 p4 ;
li4 = d n2 p4 p1 ;
*
su1 = dall li1 li2 li3 li4 'PLAN' ;
si ( ega graph 'O');
trace su1 ;
finsi;
*
* MODELE MECANIQUE
*
mod1 = MODE su1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'VISCOPLASTIQUE'
'PARFAIT' ;
*
* CARACTÉRISTIQUES DU MATÉRIAU
*
ltemps = prog 0. 200. 400. 600. 800. 1000. 1200. 1400. 1450.
1600;
* young
lyoun = prog 200.d9 200.d9 180.d9 160.d9 140.d9 80.d9 20.d9 1.d9 1.d6
1.d6 ;
evyoun = evol manu 'T' ltemps 'YOUN' lyoun ;
si ( ega graph 'O');
dess evyoun titre 'Module de Young';
finsi;
*
* nu
*
lnu = prog 0.3 0.32 0.34 0.36 0.38 0.40 0.44 0.49 0.495
0.495 ;
evnu = evol manu 'T' ltemps 'NU' lnu ;
si ( ega graph 'O');
dess evnu titre 'Coefficient de Poisson' ;
finsi;
*
* alpha
*
lalpha = prog 0.d-6 0d-6 0.d-6 0.d-6 0.d-6 0.d-6 0.d-6 0.d-6 0.d-6
0.d-6;
evalpha = evol manu 'T' ltemps 'ALPH' lalpha ;
si ( ega graph 'O');
dess evalpha titre 'Coefficient de dillatation thermique' ;
finsi;
*
* rho
*
lrho = prog 7.9d3 7.9d3 7.8d3 7.8d3 7.75d3 7.7d3 7.7d3 7.65d3 7.d3
7.d3;
evrho = evol manu 'T' ltemps 'RHO' lrho ;
si ( ega graph 'O');
dess evrho titre 'Densite' ;
finsi;
*
* N
*
lden = prog 2. 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2
1.;
evden = evol manu 'T' ltemps 'N' lden ;
si ( ega graph 'O');
dess evden titre 'Exposant N' ;
finsi;
*
* K
*
lk = prog 1.d7 1.d7 1.d7 1.d7 1.d7 1.d7 1.d7 1.d7 1.d7
1.d7;
evk = evol manu 'T' ltemps 'K' lk ;
si ( ega graph 'O');
dess evk titre 'Coefficient K' ;
finsi;
*
* SIGY
*
ly = prog 150d6 150.d6 120.d6 80.d6 60.d6 20.d6 10.d6 1.d6 0.d6
0.d6;
evy = evol manu 'T' ltemps 'SIGY' ly ;
si ( ega graph 'O');
dess evy titre 'Limite elastique' ;
finsi;
*opti donn 5;
*
*
mat1 = MATE mod1 'YOUN' evyoun 'NU' evnu 'ALPH' evalpha 'RHO' evrho
'N' evden 'SIGY' evy 'K' evk;
*
*
*
chpt20 = manu chpo su1 1 'T' 20. ;
chelt20 = changer CHAM chpt20 mod1 'RIGIDITE' ;
mat2 = vari 'NUAG' mod1 mat1 chelt20 ;
*
* CALCUL DE LA SOLUTION THEORIQUE
*
* liste des temps
linstant = prog 0. 50. 100. ;
* evolution de la contrainte
lsigma = prog 140.d6 30.d6 140.d6 ;
evsig = evol manu 'Temps' linstant 'Contrainte' lsigma ;
si ( ega graph 'O');
dess evsig titr 'Evolution de la contrainte axiale' ;
finsi;
* evolution de la température
lTdet = prog 20. 1000. 20. ;
evdet = evol manu 'Temps' linstant 'temperature' lTdet ;
si ( ega graph 'O');
dess evdet titr 'Evolution de la temperature' ;
finsi;
*opti donn 5;
*
* deplacement initial (sigma/E)
*
ezz = 0.00075 ;
lztheo = prog (ezz*0.1);
*
* pas de temps
*
dt = 0.1 ;
t1 = 0. ;
*
si complet ;nfois = 1000; sinon; nfois=450; finsi;
repeter bou1 nfois ;
t2 = dt * (&bou1) ;
t15 = t1 +t2 /2.d0 ;
TP15 = ipol t15 linstant lTdet ;
sigy15 = ipol TP15 ltemps ly ;
k15 = ipol TP15 ltemps lk ;
n15 = ipol TP15 ltemps lden ;
sigma15 = ipol t15 linstant lsigma ;
sige1 = (ipol t1 linstant lsigma)/(ipol t1 ltemps lyoun);
sige2 = (ipol t2 linstant lsigma)/(ipol t2 ltemps lyoun);
* deformation viscoplastique
syk = (sigma15-sigy15 + (abs (sigma15-sigy15)))/2.d0/k15 ;
dezp = dt*(0.6666*(syk**n15)) ;
* deformation elastique
dezz = dezp + sige2 - sige1 ;
* mess 'Pas ' &bou1 'ezz=' ezz 'dezzp='dezp ;
ezz = ezz +dezz ;
lztheo = lztheo et (prog (ezz *0.1));
t1 = t2 ;
fin bou1 ;
mess 'EZZ= theorique' ezz 'DZ theorique' (EZZ*0.1);
*opti donn 5 ;
*
* CALCUL AVEC PASAPAS
*
tab1 = table ;
*
* CONDITIONS AUX LIMITES
*
cl1 = bloq li1 'UZ' ;
cl2 = bloq li4 'UR' ;
tab1.'BLOCAGES_MECANIQUES' = cl1 et cl2 ;
*
tab1.'CARACTERISTIQUES' = mat1 ;
*
* FORCES APPILQUEES
*
f1 = press mass mod1 li3 (-1.) ;
titr 'Maillage et forces appliquees (mode axisymmetrique)';
vf1 = vect f1 'FR' 'FZ' 5.d3;
si ( EGA graph 'O');
trace vf1 su1 ;
finsi;
evf = evol manu 'T' linstant 'Intensite' lsigma ;
charf = charge evf f1 'MECA' ;
*
* CHARGEMENT THERMIQUE
*
evt = evol manu 'TPS' linstant 'T' lTdet ;
chpt = manu chpo su1 1 'T' 1. ;
chart = charg evt chpt 'T' ;
tab1.'CHARGEMENT'= chart et charf ;
*
si complet;
lcal = prog 0. pas 10. 40. 44. pas 0.25 47.5
pas 0.1 52.5 pas 0.25 55.5 pas .1
56. pas 4. 60. pas 10. 100. ;
sinon;
lcal = prog 0. pas 10. 40. 44. pas 0.25 45.;
finsi;
tab1.'TEMPS_CALCULES' = lcal ;
tab1.'TEMPS0' = 0. ;
tab1.'MODELE' = mod1 ;
tab1.'TALPHA_REFERENCE' = chpt20 ;
tab1.'MAXITERATION' = 100 ;
*
TMASAU=table;
tab1 . 'MES_SAUVEGARDES'=TMASAU;
TMASAU .'DEFTO'=VRAI;
TMASAU .'DEFIN'=VRAI;
pasapas tab1 ;
*
* POST TRAITEMENT
*
tdep1 = tab1.'DEPLACEMENTS' ;
depfin = tdep1 . ((dime tdep1) - 1 ) ;
dz = extr depfin p3 'UZ' ;
mess 'DZ calcule' dz ;
*
erre = abs ( ezz*0.1 -dz )/ ezz *0.1 ;
sauter 1 ligne ;
mess 'Erreur=' erre ;
si ( erre >eg 0.02 ) ;
erreur 5 ;
finsi ;
*
*opti sauv 'res' ;
*sauver tab1 ;
*fin ;
*------------------------------------------------------------
*opti rest 'res' ;
*rest ;
*
* COURBE TEMPS-DEPLACEMENT
*
lz = prog ; lxx = prog ;
tdep1 = tab1.'DEPLACEMENTS' ;
*
repeter bou3 (dime tdep1 ) ;
depfin = tdep1 . (&bou3 - 1 ) ;
dz = extr depfin p3 'UZ' ;
lz = lz et (prog dz );
fin bou3 ;
*
evz = evol manu 'TEMPS' lcal 'Deplacment vertical' lz ;
si complet ; nft=100.; sinon; nft=45.; finsi;
lttheo = prog 0. pas 0.1 nft ;
evztheo = evol manu 'TEMPS' lttheo 'Deplacment vertical' lztheo ;
dTAB1=TABLE;
dTAB1.1='MARQ PLUS NOLI';
dTAB1.'TITRE'=TABLE;
dTAB1.'TITRE'. 3 = MOT 'Essai de Traction' ;
si ( ega graph 'O');
dess (evz et evztheo) logo dtab1 ;
finsi;
*
*
*___________________________________________________________________
* SECONDE PARTIE
*
*
* cas test pour le modele visco plastique parfait
*
* on impose un champ de deplacement en cisaillement pur MODE PLAN DEFO
* et on regarde evoluer les contraintes
* On compare aux resultats analytiques
*
* la deformation inelastique doit varier comme 0.5*(1.-exp(-t))
* le vmises doit valoir Y = 3 ** 0.5 * 50 = 86.603
*
opti dime 2 mode plan defo elem seg3;
p1 = 0. 0. ;
p2 = 1. 0. ;
p3 = 1. 1. ;
p4 = 0. 1. ;
elem1 = 'MANU' 'QUA4' P1 P2 P3 P4 ;
mode1 = 'MODE' elem1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'VISCOPLASTIQUE'
'PARFAIT' ;
mate1 = 'MATE' mode1 'YOUN' 200. 'NU' 0. 'SIGY' (3. ** 0.5 * 50.)
'N' 1. 'K' 200. ;
bloc1 = 'BLOQ' elem1 'DEPL' ;
depimp1 = ('NOMC' ('COOR' 1 elem1 ) 'UY' ) 'ET'
('MANU' 'CHPO' elem1 1 'UX' 0. 'NATURE' 'DIFFUSE' ) ;
fldep1 = 'DEPI' bloc1 depimp1 ;
defor_0 = 'EPSI' mode1 depimp1 ;
sig_0 = 'ELAS' mode1 defor_0 mate1 ;
varin_0 = 'ZERO' mode1 'VARINTER' ;
defor_0 = 'ZERO' mode1 'DEFORMAT' ;
si complet;
ltemps = 'PROG' 0. 'PAS' 0.05 5. 'PAS' 0.5 10. ;
sinon;
ltemps = prog 0. pas 0.05 0.5;
finsi;
lchar = 'PROG' ('DIME' ltemps) * 1. ;
tab1 = 'TABLE' ;
tab1 . 'DEPLACEMENTS' = 'TABLE' ;
tab1 . 'CONTRAINTES' = 'TABLE' ;
tab1 . 'VARIABLES_INTERNES'= 'TABLE' ;
tab1 . 'DEFORMATIONS_INELATSIQUES' = 'TABLE';
tab1 . 'DEPLACEMENTS' . 0 = depimp1 ;
tab1 . 'CONTRAINTES' . 0 = sig_0 ;
tab1 . 'VARIABLES_INTERNES' . 0 = varin_0 ;
tab1 . 'DEFORMATIONS_INELATSIQUES' . 0 = defor_0 ;
ev1 = 'EVOL' 'MANU' 'TEMPS' ltemps 'INTENSITE' lchar ;
char1 = 'CHAR' 'DIMP' ev1 fldep1 ;
tab1 . 'BLOCAGES_MECANIQUES' = bloc1 ;
tab1 . 'CHARGEMENT' = char1 ;
tab1 . 'MODELE' = mode1 ;
tab1 . 'CARACTERISTIQUES' = mate1 ;
tab1 . 'TEMPS_CALCULES' = ltemps ;
TMASAU=table;
tab1 . 'MES_SAUVEGARDES'=TMASAU;
TMASAU .'DEFTO'=VRAI;
TMASAU .'DEFIN'=VRAI;
PASAPAS tab1 ;
*** solution numérique
levp = 'PROG' ;
'REPETER' bou1 ('DIME' ltemps) ;
evp = 'EXTR' ( tab1 . 'DEFORMATIONS_INELASTIQUES' . (&bou1 - 1))
'GIXY' 1 1 1 ;
levp = levp 'ET' ('PROG' evp ) ;
'FIN' bou1 ;
*** solution theorique
levpth = ('PROG' ( 'DIME' ltemps) * 1. ) - ('EXP' (ltemps * -1.))
* 0.5 ;
ev_th = 'EVOL' 'MANU' 'TEMPS' ltemps 'GIXY' levpth ;
ev_num = 'EVOL' 'MANU' 'TEMPS' ltemps 'GIXY' levp ;
tabgr = table ;
tabgr . 1 = 'MARQ CROI TIRR ' ;
tabgr . 2 = 'MARQ TRIA ' ;
tabgr.'TITRE' = table ;
tabgr .'TITRE'. 1 = MOT 'analytique' ;
tabgr .'TITRE'. 2 = MOT 'numerique' ;
si ( ega graph 'O');
dess (ev_th et ev_num ) tabgr LEGE LOGO 'YBOR' 0. 1. ;
finsi;
* contrainte equivalent à l'issue
n1 = ('DIME' ltemps) - 1 ;
vmis_fin = 'EXTR' ('VMIS' mode1 (tab1 . 'CONTRAINTES' . n1 ))
'SCAL' 1 1 1 ;
si complet ;
valthe = 3 **0.5 * 50;
sinon;
valthe = 139.78;
finsi;
err1 = abs (vmis_fin / valthe - 1. );
err2 = 'MAXI' 'ABS' (levpth - levp ) ;
mess 'ERREUR = ' err1 err2 ;
'SI' (( err1 >EG 0.01 ) 'OU' ( err2 >EG 0.01)) ;
ERREUR 5 ;
'FINSI' ;
*
*
*___________________________________________________________________
* TROISIEME PARTIE
*
*
*
* cas test pour le modele visco plastique parfait
*
* on impose un champ de deplacement en cisaillement pur MODE AXIS
* et on regarde evoluer les contraintes
* On compare aux resultats analytiques
*
* la deformation inelastique doit varier comme 0.5*(1.-exp(-t))
* le vmises doit valoir Y = 3 ** 0.5 * 50 = 86.603
*
opti dime 2 mode axis elem seg3;
p1 = 0. 0. ;
p2 = 1. 0. ;
p3 = 1. 1. ;
p4 = 0. 1. ;
elem1 = 'MANU' 'QUA4' P1 P2 P3 P4 ;
mode1 = 'MODE' elem1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' 'VISCOPLASTIQUE'
'PARFAIT' ;
mate1 = 'MATE' mode1 'YOUN' 200. 'NU' 0. 'SIGY' (3. ** 0.5 * 50.)
'N' 1. 'K' 200. ;
bloc1 = 'BLOQ' elem1 'DEPL' ;
depimp1 = ('NOMC' ('COOR' 1 elem1 ) 'UZ' ) 'ET'
('MANU' 'CHPO' elem1 1 'UR' 0. 'NATURE' 'DIFFUSE' ) ;
fldep1 = 'DEPI' bloc1 depimp1 ;
defor_0 = 'EPSI' mode1 depimp1 ;
sig_0 = 'ELAS' mode1 defor_0 mate1 ;
varin_0 = 'ZERO' mode1 'VARINTER' ;
defor_0 = 'ZERO' mode1 'DEFORMAT' ;
si complet;
ltemps = 'PROG' 0. 'PAS' 0.05 5. 'PAS' 0.5 10. ;
sinon;
ltemps = prog 0. pas 0.05 0.5;
finsi;
lchar = 'PROG' ('DIME' ltemps) * 1. ;
tab1 = 'TABLE' ;
tab1 . 'DEPLACEMENTS' = 'TABLE' ;
tab1 . 'CONTRAINTES' = 'TABLE' ;
tab1 . 'VARIABLES_INTERNES'= 'TABLE' ;
tab1 . 'DEFORMATIONS_INELATSIQUES' = 'TABLE';
tab1 . 'DEPLACEMENTS' . 0 = depimp1 ;
tab1 . 'CONTRAINTES' . 0 = sig_0 ;
tab1 . 'VARIABLES_INTERNES' . 0 = varin_0 ;
tab1 . 'DEFORMATIONS_INELATSIQUES' . 0 = defor_0 ;
ev1 = 'EVOL' 'MANU' 'TEMPS' ltemps 'INTENSITE' lchar ;
char1 = 'CHAR' 'DIMP' ev1 fldep1 ;
tab1 . 'BLOCAGES_MECANIQUES' = bloc1 ;
tab1 . 'CHARGEMENT' = char1 ;
tab1 . 'MODELE' = mode1 ;
tab1 . 'CARACTERISTIQUES' = mate1 ;
tab1 . 'TEMPS_CALCULES' = ltemps ;
TMASAU=table;
tab1 . 'MES_SAUVEGARDES'=TMASAU;
TMASAU .'DEFTO'=VRAI;
TMASAU .'DEFIN'=VRAI;
PASAPAS tab1 ;
*** solution numérique
levp = 'PROG' ;
'REPETER' bou1 ('DIME' ltemps) ;
evp = 'EXTR' ( tab1 . 'DEFORMATIONS_INELASTIQUES' . (&bou1 - 1))
'GIRZ' 1 1 1 ;
levp = levp 'ET' ('PROG' evp ) ;
'FIN' bou1 ;
*** solution theorique
levpth = ('PROG' ( 'DIME' ltemps) * 1. ) - ('EXP' (ltemps * -1.))
* 0.5 ;
ev_th = 'EVOL' 'MANU' 'TEMPS' ltemps 'GIRZ' levpth ;
ev_num = 'EVOL' 'MANU' 'TEMPS' ltemps 'GIRZ' levp ;
tabgr = table ;
tabgr . 1 = 'MARQ CROI TIRR ' ;
tabgr . 2 = 'MARQ TRIA ' ;
tabgr.'TITRE' = table ;
tabgr .'TITRE'. 1 = MOT 'analytique' ;
tabgr .'TITRE'. 2 = MOT 'numerique' ;
si ( ega graph 'O');
dess (ev_th et ev_num ) tabgr LEGE LOGO 'YBOR' 0. 1. ;
finsi;
* contrainte equivalent à l'issue
n1 = ('DIME' ltemps) - 1 ;
vmis_fin = 'EXTR' ('VMIS' mode1 (tab1 . 'CONTRAINTES' . n1 ))
'SCAL' 1 1 1 ;
mess vmis_fin;
si complet ;
valthe = 3 **0.5 * 50;
sinon;
valthe = 139.78;
finsi;
err1 = abs (vmis_fin / valthe - 1. );
err2 = 'MAXI' 'ABS' (levpth - levp ) ;
mess 'ERREUR = ' err1 err2 ;
'SI' (( err1 >EG 0.01 ) 'OU' ( err2 >EG 0.01)) ;
ERREUR 5 ;
'FINSI' ;
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer