Accueil | Cast3M

* fichier :  ortho-vari-coq4.dgibi
************************************************************************
************************************************************************
 
*** Options ...
 
    opti dime 3 mode trid ;
 
*** 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 UZ RX RY RZ ;
    motsdua = mots FX FY FZ MX MY MZ ;
 
*** 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 0 ;
    p2 = L 0 0 ;
    p3 = L L 0 ;
    p4 = 0 L 0 ;
 
*** Maillage ...
 
    su1 = manu qua4 p1 p2 p3 p4 ;
 
*** Champs de déplacement tests ...
 
    chpox chpoy chpoz = coor su1 ;
    chpo0 = manu chpo su1 1 scal 0 nature diffus ;
 
    tabdepl = table ;
    tabdepl . 1 = (epsilon * (nomc UX chpox)) et (nomc UY chpo0) et
                             (nomc UZ chpo0)  et (nomc RX chpo0) et 
                             (nomc RY chpo0)  et (nomc RZ chpo0) ;
    tabdepl . 2 = (epsilon * (nomc UX chpoy)) et (nomc UY chpo0) et
                             (nomc UZ chpo0)  et (nomc RX chpo0) et 
                             (nomc RY chpo0)  et (nomc RZ 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) 0. ;
    mo1 = mode su1 mecanique elastique orthotrope coq4 ;
    ma1 = mate mo1 DIRECTION dir1 PARALLELE YG1 val_e1
          YG2 val_e2  NU12 val_nu12 G12 val_g12 G23 val_g23
          G13 val_g13 EPAI 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.e-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 ;
* ATTENTION !!! Dans cette intégrale le champ des CARACTERISTIQUES
* n'est pas utilisé car les efforts sont déja intégrés dans l'épaisseur
* de la coque, il ne reste que l'intégrale sur la surface ...
    entot1 = intg dener mo1 ;
    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 ;