* fichier : ortho-vari-2D.dgibi
*** Options ...
opti dime 2 mode plan cont ;
opti epsilon lineaire;
*** Paramètres ...
L = 1.0 ;
epsilon = 0.001 ;
val_ep = 0.01 ;
val_e1 = 1.e+11 ;
val_e2 = 1.e+10 ;
val_e3 = 5.e+10 ;
val_nu12 = 0.3 ;
val_nu23 = 0.2 ;
val_nu13 = 0.25 ;
val_g12 = 1.e+9 ;
val_g23 = 5.e+9 ;
val_g13 = 3.e+9 ;
motspri = mots UX UY ;
motsdua = mots FX FY ;
*** Coefficients de la matrice de Hooke en repère d'orthotropie ...
denomc = 1.0 - (val_nu12*val_nu12*val_e2/val_e1) ;
Q11 = val_e1 / denomc ;
Q12 = val_nu12 * val_e2 / denomc ;
Q22 = val_e2 / denomc ;
Q66 = val_g12 ;
*** Points ...
p1 = 0 0 ;
p2 = L 0 ;
p3 = L L ;
p4 = 0 L ;
*** Maillage ...
su1 = manu qua4 p1 p2 p3 p4 ;
*** Champs de déplacement tests ...
chpox chpoy = coor su1 ;
chpo0 = manu chpo su1 1 scal 0 nature diffus ;
tabdepl = table ;
tabdepl . 1 = (epsilon * (nomc UX chpox)) et (nomc UY chpo0) ;
tabdepl . 2 = (epsilon * (nomc UX chpoy)) et (nomc UY chpo0) ;
*** Boucle sur différents déplacements imposés ...
opti echo 0 ;
repeter surdepl (dime tabdepl) ;
mess '===================================' ;
mess 'Champ de déplacement N°' &surdepl ;
depl1 = tabdepl . &surdepl ;
*** Boucle sur différents angles du repère d'orthotropie ...
angle = 0. ;
repeter surangl 12 ;
angle = angle + 7. ;
mess '---------------------------------' ;
mess ' Angle =' angle ;
*** Modèle + matériau ...
dir1 = (cos angle) (sin angle) ;
mo1 = mode su1 mecanique elastique orthotrope ;
ma1 = mate mo1 DIRECTION dir1 PARALLELE YG1 val_e1
YG2 val_e2 YG3 val_e3 NU12 val_nu12 NU23 val_nu23
NU13 val_nu13 G12 val_g12 DIM3 val_ep ;
mh1 = hook mo1 ma1 ;
ri1 = rigi mo1 ma1 ;
*** "Calculs" ...
epsi1 = epsi mo1 ma1 depl1 ;
sigma1s = sigm mo1 ma1 depl1 ;
sigma1e = elas mo1 epsi1 ma1 ;
sigma1h = mh1 * epsi1 mo1 ;
test12 = abs (sigma1s - sigma1e) ;
test23 = abs (sigma1e - sigma1h) ;
si((maxi test12) > 1.5e-8) ;
mess 'Erreur !!! Différence entre opérateurs ELAS et SIGM' ;
erre 5 ;
finsi ;
si((maxi test23) > 1.e-8) ;
mess 'Erreur !!! Différence entre opérateurs ELAS et HOOK' ;
erre 5 ;
finsi ;
dener = (ener mo1 epsi1 sigma1s) / 2 ;
entot1 = intg dener mo1 ma1 ;
mess 'Énergie de déformation : /// sigma * epsilon =' entot1 ;
fr_kxu = ri1 * depl1 ;
entot2 = (xty depl1 fr_kxu motspri motsdua) / 2 ;
mess 'Énergie de déformation : (K * u) . u =' entot2 ;
fr_bsig = bsig mo1 sigma1s ma1 ;
entot3 = (xty depl1 fr_bsig motspri motsdua) / 2 ;
mess 'Énergie de déformation : (/// B^T * sigma) . u =' entot3 ;
test_fr = maxi (abs (fr_kxu - fr_bsig)) ;
si(test_fr > 1.e-5) ;
mess 'ERREUR !!! Forces de réaction incorrectes !!!' ;
erre 5 ;
finsi ;
cosa1 = cos angle ;
sina1 = sin angle ;
cosa2 = cosa1 * cosa1 ;
sina2 = sina1 * sina1 ;
cosa4 = cosa2 * cosa2 ;
sina4 = sina2 * sina2 ;
si(ega &surdepl 1) ;
Q11bar = (Q11*cosa4) + (2*cosa2*sina2*(Q12+(2*Q66))) +
(Q22*sina4) ;
en_ana = Q11bar * epsilon*epsilon * L*L * val_ep / 2 ;
finsi ;
si(ega &surdepl 2) ;
Q66bar = (Q66*(cosa4+sina4)) + (sina2*cosa2*(Q11+Q22-(2*Q12)
-(2*Q66))) ;
en_ana = Q66bar * epsilon*epsilon * L*L * val_ep / 2 ;
finsi ;
mess 'Énergie donnée par la formule analytique =' en_ana ;
test12 = abs (entot1 - entot2) ;
test23 = abs (entot2 - entot3) ;
test34 = abs (entot3 - en_ana) ;
* list test12 ;
* list test23 ;
* list test34 ;
testgl = (test12 > 1.e-10) ou (test23 > 1.e-10) ou
(test34 > 1.e-10) ;
si(testgl) ;
mess 'ERREUR !!! Écart intolérable !!!' ;
erre 5 ;
finsi ;
*** Fin de la boucle sur les angles ...
fin surangl ;
*** Fin de la boucle sur les déplacements ...
fin surdepl ;
*** Bye ...
saut 1 lign ;
opti echo 1 ;
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer