Accueil | Cast3M

* 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 1 ;
 
*** 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
 
    MOP = 'MODE' SU1 'CHARGEMENT' 'PRESSION' elemact ;
    MAP = 'PRES' MOP 'PRES' 10000. ;
    fopres = 'BSIG' MOP MAP ;
 
    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 . TEMPS_CALCULES = prog 0.0 pas 1.e-1 1.0 ;
       tabnl . TEMPS_SAUVES = lr1 ;
       tabnl . GRANDS_DEPLACEMENTS = vrai ;
       tabn1 . 'HYPOTHESE_DEFORMATIONS'='LINEAIRE' ;
 
       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 'LINE' ;
    sinon ;
       si1 = sigm mo1 ma1 de1 'QUAD' ;
    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 'LINE' ;
    sinon ;
       ep1 = epsi mo1 de1 ma1 'QUAD' ;
    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 ;