forma04aためしてみました

Hertzの接触応力の計算です。平面要素を使っています。

要素は次のようになっています。

commファイルです。

DEBUT(identifier='1:1',
      CODE=_F(NIV_PUB_WEB='INTERNET'),
      DEBUG=_F(SDVERI='OUI'))

E = 20000.

nu = 0.3

h = 2.

h_m = -1. * h

R = 50.

# -- Pression de hertz et demi-largeur de contact
p0 = -E / (pi * (1.0 - nu * nu)) * sqrt(2.0 * 2.0 * h / R)

a = sqrt(R * h)

MA = LIRE_MAILLAGE(identifier='12:1',
                   FORMAT='MED',
                   UNITE=20)

MA = MODI_MAILLAGE(identifier='13:1',
                   reuse=MA,
                   MAILLAGE=MA,
                   ORIE_PEAU=_F(GROUP_MA_PEAU=('CONT1', 'CONT2')))

# -- Groupe de noeuds oriente depuis le centre de la zone de contact esclave
MA = DEFI_GROUP(identifier='15:1',
                reuse=MA,
                CREA_GROUP_NO=_F(GROUP_MA='CONT1',
                                 GROUP_NO_ORIG='C1',
                                 NOM='CONT_ord',
                                 OPTION='NOEUD_ORDO'),
                MAILLAGE=MA)

MO = AFFE_MODELE(identifier='16:1',
                 AFFE=_F(GROUP_MA=('SPH_HAU', 'SPH_BAS'),
                         MODELISATION='AXIS',
                         PHENOMENE='MECANIQUE'),
                 MAILLAGE=MA)

MAT = DEFI_MATERIAU(identifier='17:1',
                    ELAS=_F(E=E,
                            NU=nu))

CHMAT = AFFE_MATERIAU(identifier='18:1',
                      AFFE=_F(GROUP_MA=('SPH_HAU', 'SPH_BAS'),
                              MATER=MAT),
                      MAILLAGE=MA)

RAMPE = DEFI_FONCTION(identifier='19:1',
                      NOM_PARA='INST',
                      VALE=(0, 0, 1, 1))

TPS = DEFI_LIST_REEL(identifier='20:1',
                     DEBUT=0,
                     INTERVALLE=_F(JUSQU_A=1,
                                   NOMBRE=20))

V_ABSC = DEFI_LIST_REEL(identifier='21:1',
                        DEBUT=0,
                        INTERVALLE=_F(JUSQU_A=10,
                                      NOMBRE=100))

CL = AFFE_CHAR_MECA(identifier='22:1',
                    DDL_IMPO=(_F(DX=0,
                                 GROUP_MA='A1A2'),
                              _F(DY=h_m,
                                 GROUP_MA='A2B2'),
                              _F(DY=h,
                                 GROUP_MA='A1B1')),
                    MODELE=MO)

# -- Definition du contact
CONT = DEFI_CONTACT(identifier='24:1',
                    FORMULATION='DISCRETE',
                    MODELE=MO,
                    ZONE=_F(ALGO_CONT='CONTRAINTE',
                            GROUP_MA_ESCL='CONT1',
                            GROUP_MA_MAIT='CONT2'))

STNL = STAT_NON_LINE(identifier='25:1',
                     CHAM_MATER=CHMAT,
                     COMPORTEMENT=_F(RELATION='ELAS'),
                     CONTACT=CONT,
                     EXCIT=_F(CHARGE=CL,
                              FONC_MULT=RAMPE),
                     INCREMENT=_F(LIST_INST=TPS),
                     MODELE=MO,
                     SOLVEUR=_F())

# -- Post-traitement
STNL = CALC_CHAMP(identifier='27:1',
                  reuse=STNL,
                  CONTRAINTE='SIEF_ELNO',
                  CRITERES=('SIEQ_ELGA', 'SIEQ_ELNO'),
                  RESULTAT=STNL)

STNL = CALC_CHAMP(identifier='28:1',
                  reuse=STNL,
                  CONTRAINTE='SIEF_NOEU',
                  CRITERES='SIEQ_NOEU',
                  RESULTAT=STNL)

# -- Pression de contact extraite sur le bord oriente
SIYY_abs = POST_RELEVE_T(identifier='30:1',
                         ACTION=_F(GROUP_NO='CONT_ord',
                                   INST=1.0,
                                   INTITULE='SIYY en fonction de ABSC_CURV',
                                   NOM_CHAM='SIEF_NOEU',
                                   NOM_CMP='SIYY',
                                   OPERATION='EXTRACTION',
                                   RESULTAT=STNL))

IMPR_TABLE(identifier='31:1',
           TABLE=SIYY_abs,
           UNITE=4)

# -- Trace de la pression et de la solution analytique
PRES = FORMULE(identifier='33:1',
               NOM_PARA='x',
               VALE='p0*sqrt(1.0-(x/a)*(x/a))',
               a=a,
               p0=p0)

PRES_sol = CALC_FONC_INTERP(identifier='34:1',
                            FONCTION=PRES,
                            LIST_PARA=V_ABSC,
                            NOM_PARA='x')

SIYY_f = RECU_FONCTION(identifier='35:1',
                       PARA_X='ABSC_CURV',
                       PARA_Y='SIYY',
                       TABLE=SIYY_abs)

# -- Impressions des resultats
IMPR_FONCTION(identifier='37:1',
              BORNE_X=(0.0, 20.0),
              BORNE_Y=(-3000.0, 0.0),
              COURBE=(_F(FONCTION=SIYY_f),
                      _F(FONCTION=PRES_sol)),
              FORMAT='XMGRACE',
              GRILLE_X=1.0,
              GRILLE_Y=500.0,
              LEGENDE_X='Abscisse curviligne (mm)',
              LEGENDE_Y='Pression (Mpa)',
              PILOTE='',
              TITRE='Pression de Hertz',
              UNITE=2)

IMPR_RESU(identifier='38:1',
          FORMAT='MED',
          RESU=_F(RESULTAT=STNL),
          UNITE=3)

# -- Validation des resultats
TEST_RESU(identifier='40:1',
          RESU=_F(CRITERE='RELATIF',
                  GROUP_NO='C1',
                  INST=1.0,
                  NOM_CHAM='SIEF_NOEU',
                  NOM_CMP='SIYY',
                  PRECISION=0.01,
                  REFERENCE='ANALYTIQUE',
                  RESULTAT=STNL,
                  VALE_CALC=-2773.5474109439,
                  VALE_REFE=-2798.3286697476105))

FIN(identifier='41:1',
    )