板（plate、ソリッド）とブロック（stamp、ソリッド）の接触解析です。
概要はこんな感じです。突っ込まれる側のplateを少し傾けて、片あたり状態にしています。
stampの材料特性は、バイリニアとしています。
plate側は弾性材料です。
Geomは下図のようになります。
Meshは下図のとおり。plateはmesh1、stampはmesh2です。
commです。まずはplateのmeshを読み込みます。
DEBUT(LANG='FR')
pmesh = LIRE_MAILLAGE(identifier='0:1',
UNITE=20)
次にstampのmeshを読み込みます。
mesh = LIRE_MAILLAGE(identifier='1:1',
UNITE=2)
stampは宙に浮いていて、解析開始時に不安定なので、上の四隅にばねを追加します。
smesh = CREA_MAILLAGE(identifier='2:1',
CREA_POI1=_F(GROUP_NO=('RBT', 'LBT', 'RFT', 'LFT'),
NOM_GROUP_MA=('spElmt', )),
MAILLAGE=mesh)
plateとstampのmeshを統合します。
nmesh = ASSE_MAILLAGE(identifier='3:1',
MAILLAGE_1=pmesh,
MAILLAGE_2=smesh,
OPERATION='SUPERPOSE')
モデルを定義します。
model = AFFE_MODELE(identifier='4:1',
AFFE=(_F(MODELISATION='3D',
PHENOMENE='MECANIQUE',
TOUT='OUI'),
_F(GROUP_MA=('spElmt', ),
MODELISATION='DIS_T',
PHENOMENE='MECANIQUE')),
MAILLAGE=nmesh)
安定用のばねの定義を行います。
softSp = AFFE_CARA_ELEM(identifier='5:1',
DISCRET=_F(CARA='K_T_D_N',
GROUP_MA=('spElmt', ),
VALE=(50.0, 50.0, 50.0)),
MODELE=model)
材料の定義です。
steel = DEFI_MATERIAU(identifier='6:1',
ELAS=_F(E=206800.0,
NU=0.29))
ssteel = DEFI_MATERIAU(identifier='7:1',
ECRO_LINE=_F(D_SIGM_EPSI=3283.0,
SY=155.8),
ELAS=_F(E=206800.0,
NU=0.29))
plate用がsteel、stamp用がssteelです。これをモデルに割り振ります。
fieldmat = AFFE_MATERIAU(identifier='8:1',
AFFE=(_F(GROUP_MA=('plate', ),
MATER=(steel, )),
_F(GROUP_MA=('stamp', ),
MATER=(ssteel, ))),
MODELE=model)
非線形解析用の関数を定義します。funcが刻み。listr0がfuncに対する荷重増加の関数、listrは出力するステップの定義です。
func = DEFI_FONCTION(identifier='9:1',
INTERPOL=('LIN', ),
NOM_PARA='INST',
VALE=(0.0, 0.0, 1.0, 1.0))
listr0 = DEFI_LIST_REEL(identifier='10:1',
DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=1.0,
NOMBRE=1600))
listr = DEFI_LIST_REEL(identifier='11:1',
DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=1.0,
NOMBRE=40))
今回初トライの最大せん断応力を出力するために関数です。trescaを1/2倍するだけですけど。
MAXTAW = FORMULE(identifier='12:1',
NOM_PARA=('TRESCA', ),
VALE='TRESCA/2.0')
境界条件です。plateを四辺支持状態にしています。あと、stampがZ方向に下がる（分力でX方向やY方向に動かない）ために拘束しています。
fix = AFFE_CHAR_MECA(identifier='13:1',
DDL_IMPO=(_F(DX=0.0,
DY=0.0,
DZ=0.0,
GROUP_NO=('bedge', )),
_F(DX=0.0,
DY=0.0,
GROUP_NO=('CF', )),
_F(DX=0.0,
GROUP_NO=('RF', )),
_F(DX=0.0,
GROUP_NO=('LF', ))),
MODELE=model)
外力として、stamp上面に80MPa加えます。
load = AFFE_CHAR_MECA(identifier='14:1',
FORCE_FACE=_F(FZ=-80.0,
GROUP_MA=('FTOP', )),
MODELE=model)
接触の定義です。
contact = DEFI_CONTACT(identifier='15:1',
ALGO_RESO_CONT='NEWTON',
ALGO_RESO_GEOM='NEWTON',
CONT_STAT_ELAS=0,
FORMULATION='CONTINUE',
MODELE=model,
RESI_CONT=-1.0,
RESI_GEOM=0.0001,
ZONE=_F(ALGO_CONT='STANDARD',
GROUP_MA_ESCL=('sbottom', ),
GROUP_MA_MAIT=('ptop', )))
非線形解析の定義です。完走させるため、収束判定を２段階にしています。
１段目。
resnonl = STAT_NON_LINE(identifier='16:1',
ARCHIVAGE=_F(LIST_INST=listr),
CARA_ELEM=softSp,
CHAM_MATER=fieldmat,
COMPORTEMENT=(_F(GROUP_MA=('stamp', ),
RELATION='VMIS_ISOT_LINE'),
_F(GROUP_MA=('plate', ),
RELATION='ELAS'),
_F(GROUP_MA=('spElmt', ),
RELATION='ELAS')),
CONTACT=contact,
CONVERGENCE=_F(ITER_GLOB_ELAS=50,
ITER_GLOB_MAXI=200),
EXCIT=(_F(CHARGE=fix),
_F(CHARGE=load,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr0,
NUME_INST_FIN=800,
NUME_INST_INIT=0),
METHODE='NEWTON',
MODELE=model,
NEWTON=_F(MATRICE='TANGENTE',
REAC_INCR=1,
REAC_ITER=1))
収束条件を変えて（甘くして）、２段目です。
resnonl = STAT_NON_LINE(identifier='17:1',
reuse=resnonl,
ARCHIVAGE=_F(LIST_INST=listr),
CARA_ELEM=softSp,
CHAM_MATER=fieldmat,
COMPORTEMENT=(_F(GROUP_MA=('stamp', ),
RELATION='VMIS_ISOT_LINE'),
_F(GROUP_MA=('plate', ),
RELATION='ELAS'),
_F(GROUP_MA=('spElmt', ),
RELATION='ELAS')),
CONTACT=contact,
CONVERGENCE=_F(ITER_GLOB_ELAS=50,
ITER_GLOB_MAXI=200,
RESI_GLOB_RELA=2e-06),
ETAT_INIT=_F(EVOL_NOLI=resnonl,
INST_ETAT_INIT=0.5),
EXCIT=(_F(CHARGE=fix),
_F(CHARGE=load,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr0,
NUME_INST_FIN=1600,
NUME_INST_INIT=800),
METHODE='NEWTON',
MODELE=model,
NEWTON=_F(MATRICE='TANGENTE',
REAC_INCR=1,
REAC_ITER=1),
RESULTAT=resnonl)
Mises応力を計算したりする後処理の定義です。
resnonl = CALC_CHAMP(identifier='18:1',
reuse=resnonl,
CONTRAINTE=('SIGM_NOEU', ),
CRITERES=('SIEQ_NOEU', ),
RESULTAT=resnonl)
最大せん断応力の計算です。
resnonl = CALC_CHAMP(identifier='19:1',
reuse=resnonl,
CHAM_UTIL=_F(FORMULE=(MAXTAW, ),
NOM_CHAM='SIEQ_NOEU',
NUME_CHAM_RESU=1),
RESULTAT=resnonl)
結果の出力です。
IMPR_RESU(identifier='20:1',
RESU=_F(RESULTAT=resnonl,
TOUT='OUI'),
UNITE=3)
FIN()
Mises応力のコンターです。stampの半分くらいに155MPa（赤色）の応力が生じています。
原点をとおるXZ面のMisesのコンターはつぎのとおり。
ほぼ原点（x=-0.05、ｙとｚは０）をとおるＹＺ面です。
今回のファイルはこちら。
コメント