* fichier : ktest-calp.dgibi
***********************************************************
* *
* ktest-calp.dgibi *
* *
* Test de fonctionnement de l'opérateur CALP sur un *
* simple exemple de plaque carrée en flexion pure. *
* Plusieurs types de coques sont testées en même temps. *
* *
* Auteur : Michel Bulik *
* *
* Date : Août 1995 *
* *
***********************************************************
*** Options ...
opti dime 3 mode trid elem seg2 echo 0 ;
* ajout de option epsilon lineaire pour la precision des test!
OPTION epsilon lineaire;
*** Paramètres ...
* longueur et largeur de la plaque
L = 1. ;
* nombre de divisions le long de ses côtés
nbdiv = 10 ;
* module de YOUNG
valyoun = 2.e+11 ;
* coeff. de Poisson nul pour qu'il n'y ait pas
* de courbure anticlastique
valnu = 0. ;
* épaisseur de la plaque
ep = 1.e-2 ;
* calcul linéaire ou non
calclin = vrai ;
* affichage graphique ou non
graph = faux ;
* types d'éléments testés
typdel = mots COQ4 DKT DST COQ3 ;
*** Points ...
dens (L / nbdiv) ;
p1 = 0 0 0 ;
p2 = L 0 0 ;
vectz = 0 0 L ;
vectxz = L 0 L ;
*** Lignes ...
li1 = p1 d nbdiv p2 ;
*** Surface ...
opti elem qua4 ;
su1 = li1 tran nbdiv vectz ;
*** Début de la boucle sur les types d'éléments ...
nbfois = dime typdel ;
repeter surelem nbfois ;
elemact = extr typdel &surelem ;
si(ega &surelem 2) ;
mess 'SU1 passe au type TRI3' ;
su1 = chan tri3 su1 ;
finsi ;
mess '--------------------------------------' ;
mess 'On analyse actuellement l élément ' elemact ;
*** Modèle ...
mo1 = mode su1 mecanique elastique elemact ;
*** Matériau ...
ma1 = mate mo1 YOUN valyoun NU valnu ;
ca1 = cara mo1 EPAI ep ;
ma1 = ma1 et ca1 ;
*** Rigidité + Blocages ...
si(calclin) ;
ri1 = rigi mo1 ma1 ;
finsi ;
cl1 = (bloq uy uz li1) et
(bloq rx li1) et
(bloq ux p1 ) ;
*** Chargement ...
* un champ de forces obtenus par l'intégration avec les
* fonctions de forme - on va le modifier pour obtenir
* le chargement réparti uniformément sur le bord supérieur
* de la plaque
fopres = pres coqu mo1 10000. NORM ;
mafo1 = su1 poin droi vectz vectxz ;
fo1 = exco fy (redu fopres mafo1) mx ;
* moment appliqué par unité de longueur
momtot = (abs (maxi (resu fo1))) / L ;
mess 'Moment appliqué par (par unité de longueur) = '
momtot ;
*** L'affichage éventuel du chargement ...
vecchar = faux ;
si(exis (extr fo1 comp) FX) ;
titr 'Chargement : forces nodales' ;
vecfo1 = vect fo1 FX FY FZ (1 / (maxi (abs fo1))) ;
vecchar = vrai ;
finsi ;
si (exis (extr fo1 comp) MX) ;
titr 'Chargement : couples nodaux' ;
vecfo1 = vect fo1 MX MY MZ (1 / (maxi (abs fo1))) ;
vecchar = vrai ;
finsi ;
si(graph et vecchar) ;
trac vecfo1 su1 ;
finsi ;
*** Résolution ...
si(calclin) ;
de1 = reso (ri1 et cl1) fo1 ;
sinon ;
tabnl = table ;
tabnl . MODELE = mo1 ;
tabnl . CARACTERISTIQUES = ma1 ;
tabnl . BLOCAGES_MECANIQUES = cl1 ;
lr1 = prog 0. 1. ;
ev1 = evol manu 't' lr1 'f(t)' lr1 ;
tabnl . CHARGEMENT = char MECA fo1 ev1 ;
tabnl . TEMPS_CALCULES = prog 0.0 pas 1.e-1 1.0 ;
tabnl . TEMPS_SAUVES = lr1 ;
tabnl . GRANDS_DEPLACEMENTS = vrai ;
tabnl . MTOL = 1.e+30 ;
pasapas tabnl ;
de1 = tabnl . DEPLACEMENTS . 1 ;
finsi ;
*** Affichage éventuel de la déformée et des DDL ...
defo0 = defo su1 de1 0. ;
titr 'La deformee' ;
defo1 = defo su1 de1 1. ;
si(graph) ;
trac (defo0 et defo1) ;
* lescomp = extr (redu de1 su1) comp ;
lescomp = mots UY RX ;
nbcomp = dime lescomp ;
repe surcomp nbcomp ;
lacomp = extr lescomp &surcomp ;
titr 'Element' elemact ' : Degre de liberte'
lacomp ;
trac (exco lacomp de1) su1 ;
fin surcomp ;
finsi ;
*** Contraintes ...
si(calclin) ;
si1 = sigm mo1 ma1 de1 ;
sinon ;
si1 = sigm mo1 ma1 de1 II ;
finsi ;
si1 = rten si1 mo1 (1 0 0) (0 0 1) ;
*** L'affichage éventuel des contraintes ...
si(graph) ;
lescomp = extr si1 comp ;
nbcomp = dime lescomp ;
repe surcom2 nbcomp ;
lacomp = extr lescomp &surcom2 ;
titr 'Element' elemact ' : Contrainte generalisee'
lacomp ;
trac (exco lacomp si1) mo1 ;
fin surcom2 ;
finsi ;
*** Vérif de la valeur de la contrainte ...
valmax = maxi (exco m22 si1) ;
valmin = mini (exco m22 si1) ;
si( (valmax - valmin) > 1.e-6 ) ;
mess 'Attention ! Résultat non homogène (M22)' ;
* valm22 = (valmin + valmax) / 2 ;
valm22 = (intg (exco m22 si1) mo1) / (L * L) ;
sinon ;
valm22 = valmax ;
finsi ;
ecartrel = 100 * ((valm22 - momtot) / momtot) ;
mess 'Ecart relatif (M22) = ' ecartrel '%' ;
si(ecartrel > 1.e-5) ; erreur 5 ; finsi ;
*** Application de CALP ...
si1sup = calp si1 ca1 mo1 SUPE ;
si1inf = calp si1 ca1 mo1 INFE ;
valtheor = (6 * momtot) / (ep * ep) ;
valsup = (intg (exco smyy si1sup) mo1) / (L * L) ;
valinf = (intg (exco smyy si1inf) mo1) / (L * L) ;
ecsup = 100 * ((valsup - valtheor) / valtheor) ;
mess 'Ecart relatif (peau supérieure) = ' ecsup '%' ;
si(ecsup > 1.e-5) ; erreur 5 ; finsi ;
ecinf = 100 * ((valinf + valtheor) / valtheor) ;
mess 'Ecart relatif (peau inférieure) = ' ecinf '%' ;
si(ecinf > 1.e-5) ; erreur 5 ; finsi ;
*** Résultats analytiques pour le déplacement
*** et la rotation ...
* inertie par unité de longueur ...
I = (ep * ep * ep) / 12 ;
depanal = (momtot * L * L) / (2 * valyoun * I) ;
depcal = (maxi (exco uy de1)) ;
ecdep = 100 * ((depcal - depanal) / depanal) ;
mess 'Écart relatif (déplacement maxi) = ' ecdep '%' ;
si(ecdep > 0.5) ; erreur 5 ; finsi ;
rotanal = (momtot * L) / (valyoun * I) ;
rotcal = abs (mini (exco rx de1)) ;
ecrot = 100 * ((rotcal - rotanal) / rotanal) ;
mess 'Écart relatif (rotation maxi) = ' ecrot '%' ;
si(ecrot > 2.0) ; erreur 5 ; finsi ;
*** Déformations généralisées ...
si (calclin) ;
ep1 = epsi mo1 de1 ma1 ;
sinon ;
ep1 = epsi mo1 de1 ma1 II ;
finsi ;
ep1 = rten ep1 mo1 (1 0 0) (0 0 1) ;
*** et sur les peaux ...
ep1sup = calp ep1 ca1 mo1 SUPE ;
ep1inf = calp ep1 ca1 mo1 INFE ;
*** Comparaison avec les valeurs analytiques ...
* RTTT
valmax = maxi (exco rttt ep1) ;
valmin = mini (exco rttt ep1) ;
si( (valmax - valmin) > 1.e-6 ) ;
mess 'Attention ! Résultat non homogène (RTTT)' ;
valrttt = (intg (exco rttt ep1) mo1) / (L * L) ;
sinon ;
valrttt = valmax ;
finsi ;
valanal = momtot / (valyoun * I) ;
ecartrel = 100 * ((valrttt - valanal) / valanal) ;
mess 'Ecart relatif (RTTT) = ' ecartrel '%' ;
si(ecartrel > 1.e-5) ; erreur 5 ; finsi ;
* EPYY (peau supérieure)
valanal = (ep * valrttt) / 2 ;
valmax = maxi (exco epyy ep1sup) ;
valmin = mini (exco epyy ep1sup) ;
si( (valmax - valmin) > 1.e-6 ) ;
mess 'Attention ! Résultat non homogène (EPYY sup)' ;
valepyys = (intg (exco epyy ep1sup) mo1) / (L * L) ;
sinon ;
valepyys = valmax ;
finsi ;
ecartyys = 100 * ((valepyys - valanal) / valanal) ;
mess 'Ecart relatif (EPYY sup) = ' ecartyys '%' ;
si(ecartyys > 1.e-5) ; erreur 5 ; finsi ;
valmax = maxi (exco epyy ep1inf) ;
valmin = mini (exco epyy ep1inf) ;
si( (valmax - valmin) > 1.e-6 ) ;
mess 'Attention ! Résultat non homogène (EPYY inf)' ;
valepyyi = (intg (exco epyy ep1inf) mo1) / (L * L) ;
sinon ;
valepyyi = valmax ;
finsi ;
ecartyyi = 100 * ((valepyyi + valanal) / valanal) ;
mess 'Ecart relatif (EPYY inf) = ' ecartyyi '%' ;
si(ecartyyi > 1.e-5) ; erreur 5 ; finsi ;
*** Fin de la boucle sur les types d'éléments ...
fin surelem ;
*** Bye ...
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer