Download ortho-vari-coq4.dgibi

Back to the list

* 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 ;

 

© Cast3M 2003 - All rights reserved.
Disclaimer