この資料はオープンCAE勉強会＠岐阜で公開されているFS氏ご提供の「SALOME-Mecaの使用法解説：6.0 接触解析の基本(1)」をSalome-Meca 2019（Code_Aster 14.4）を用いて実行しました。
はじめに
接触解析において、Salome-Meca 2010.1以降からCode_Asterのコマンド体系が変更され、従来の解析コードそのままでは、エラーが発生し、接触解析の計算ができなくなっている。このため、新しいコマンド体系で接触解析をしてみる。
2部品同士を変位拘束（1部品を固定、1部品を変位させる）したときの接触解析を行なってみる。この場合は、1部品とも変位拘束されているので、剛体移動は発生せず、普通に解ける。この問題は、接触面の状態が安定（変形と共に接触位置や面積が変化しない）している場合は、Solid同士を連結させて解析することと同じになる。しかし、実際は、接触面ですべりが生じたり、変形と共に接触位置が変わっていく。
このような問題を解くためには、負荷（荷重や変位）を少しずつ掛けていき、その都度、解を求めて最終的な解を求める方法をとる。$\rightarrow$ 非線形解析
ここでの問題は、接触面の摩擦がなく、すべりが発生する問題を考えている。
接触判定（接触しているかどうかの判定）は、Toleranceを確認した結果、デフォルトで「5e-3」に設定されていた。つまり5μm以下で接触していると判定する。
モデルの読み込み
モデルは、連結問題で使用したモデルをそのまま使う。「multi-bar-1.stp」をmm単位で読み込む。
モデルを読み込んだ後は、モデルサイズを、「Inspection」＞「Dimensions」＞「Bounding Box」で確認しておく。解析は、Barの上面（press面）を-0.2mmZ方向に変位させる接触問題として解析してみる。
解析は、~/CAE/contact-bar/というフォルダーを作り、この中で解析する。
Entityの作成
連結問題と同様に解析で使用するVolumeやFaceをグループ化しておく。
ツリーの構造は下記。
また、前項でモデルサイズを確認しており、モデルはmm単位で作成されているので、変位の境界条件は、mm単位で入力することになる。
| グループ名 | 備考 |
|---|---|
| Base | Solid1(Base) |
| fix | 固定面 |
| contBase | Baseの接触面 |
| Bar | Solid2(Bar) |
| contBar | Barの接触面 |
| press | 荷重を負荷する面 |
メッシュの作成
連結問題と同じ方法でメッシュを切る。
ツリーは、下記。
六面体の2次メッシュとし、ジオメトリからグループを作成する。
解析コードの作成
画面をAsterStudyモジュールに変えて、アシスタント（ウィザード）を使って、通常通りCode_Asterの解析コードを作成する。このとき、固定面はfix面（X=Y=Z=0）、荷重面はpress面で0.1MPaとしておく。
材料定数は、ベリ銅の値をそのまま使用。
- ヤング率：130,300MPa
- ポアソン比：0.343
Code_Asterの結果ファイルは、フォルダー~/CAE/contact-bar/内に「multi-bar.rmed」として保存されるようにする。
解析コードの編集
AsterStudyモジュールを使って、作成された解析コードを接触問題が解けるように編集する。
従来までのCode_Asterは、接触のコマンドが境界条件を設定するAFFE_CHAR_MECAコマンドの下にContactコマンドがあったが、Salome-Meca 2010からは、最上位にDEFI_CONTACTコマンドが準備されるようになった。
境界条件の編集
境界条件は、
- 通常の境界条件
- 負荷を少しづつ変化させる条件
の2種類の条件に分けて設定する。以下に各々の境界条件設定法方について示す。
通常の境界条件
通常の境界条件は、fix面を固定する条件となる。
この境界条件は、アシスタント（ウィザード）が作成した境界条件を編集して、作成する。
press面をZ方向に-0.2mm変位させるが、XY方向の拘束がないので、XY方向も拘束する必要がある。ここでXY方向を拘束（X=Y=0）する。 下記参照。
mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
DY=0.0,
DZ=0.0,
GROUP_MA=('fix', )), # 固定する面(fix)を固定
_F(DX=0.0,
DY=0.0,
GROUP_MA=('press', ))), # 負荷をかける面(press)のXY方向を固定
MODELE=model)
少しづつ負荷させる境界条件作成
press面をZ方向に-0.2mm変位させるが、この変位が接触面に直接影響を与えるので、この変位を少しづつ変化させていくようにする必要がある。
このため、この境界条件を独立させて定義する。
また、アシスタント（ウィザード）で入力した圧力の条件も削除する。
現在設定されているAFFE_CHAR_MECAの後に、以下を追加する。
DZは、モデルの大きさに合わせて、設定する。
今回のモデルは、mm単位で作成されていたので、変位DZは、-0.2に設定している。
mecach = AFFE_CHAR_MECA(DDL_IMPO=_F(DZ=-0.2, # Z方向に-0.2mm変位させる
GROUP_MA=('press', )), # 負荷をかける面(press)を
MODELE=model)
接触の定義
ここは、Salome-Meca 2010 で新しく設定されたコマンドになる。従来は、境界条件（AFFE_CHAR_MECA）内で設定していた。この接触の定義をDEFI_CONTACTで追加する。この内容を以下で作成した。（ほとんどデフォルトのまま）
contact = DEFI_CONTACT(FORMULATION='CONTINUE', # 拡張ラグランジュ法（摩擦接触に対しては安定しており、推奨設定）
MODELE=model,
ZONE=_F(CONTACT_INIT='OUI', # 接触想定領域、初期条件
GROUP_MA_ESCL=('contBase', ), # スレイブ面
GROUP_MA_MAIT=('contBar', ))) # マスター面
接触のためのコード追加
引き続き、次の行に、接触問題を解くためのファンクションを追加する。
press面の変位を0から0.2mmまで徐々に変位させていく方法を取るため、0〜0.2mmまでの中間の値をどのように設定するか（線形or非線形で回帰）を設定する。普通に線形で回帰させる（ramp制御）方法とする。
このためのファンクションを下記のように定義する。
値は、倍率を表しており、「1」は、-0.2mmを示している。
座標の入力は、X,Yの形式でXYのペアで入力する。
func = DEFI_FONCTION(NOM_PARA='INST', # 名称は任意で可。NOM_PARAは「INST(Time)」を選択。変数は、VALEで入力
VALE=(0.0, 0.0, 1.0, 1.0)) # 原点(0,0)から(1,1)までを線形で回帰する
次に1.0（1.0倍）までを何分割して解析するのかを定義する。 下記参照。
listr = DEFI_LIST_REEL(DEBUT=0.0, # 初期値を設定
INTERVALLE=_F(JUSQU_A=1.0, # 0～1秒までを
PAS=0.2)) # 0.2秒毎に5分割する。
非線形解析方法の設定
ここで今までに設定した条件、ファンクションを使って、非線形（接触）問題を解く方法を設定する。
アシスタント（ウィザード）で設定したMECA_STATIQUEは削除する。
Salome-Meca 2010では、solver（STAT_NON_LINE）内にcontactコマンドが追加されているので、以下の様に追記した。（必要最小限の変更にした。）
以下のコードがSTAT_NON_LINEの内容。
resnonl = STAT_NON_LINE(CHAM_MATER=materfl, # 材料を指定
CONTACT=contact, # 接触を読み込む
EXCIT=(_F(CHARGE=mecabc), # 通常の境界条件
_F(CHARGE=mecach, # 少しづつ負荷させる条件
FONC_MULT=func)), # 中間の変位を線形で求める
INCREMENT=_F(LIST_INST=listr), # 0.2秒づつ増える
MODELE=model, # モデルを指定
NEWTON=_F(REAC_ITER=1)) # 接線剛性を更新する「ニュートン・ラプソン法」を指定
Post処理の追加
Post処理を追加する。
CALC_CHAMPの名前は、STAT_NON_LINEで指定した名前を再利用する。
resnonl = CALC_CHAMP(reuse=resnonl,
CONTRAINTE=('SIGM_ELNO', 'SIGM_NOEU'), # CONTRAINTE（応力）
CRITERES=('SIEQ_ELNO', 'SIEQ_NOEU'), # CRITERES（基準）
MODELE=model,
RESULTAT=resnonl)
結果出力の修正
MECA_STATIQUEを削除したとき、これにリンクされている結果出力側（IMPR_RESU）がエラーになるので、この再設定が必要になる。追加したPost処理で再設定を行う。
IMPR_RESU(FORMAT='MED', # 出力するバイナリ形式にMED形式を指定
RESU=_F(NOM_CHAM=('DEPL', 'SIEQ_NOEU', 'SIGM_NOEU'), # DEPL（変位量）、SIEQ_NOEU（等価応力（節点））、SIGM_NOEU（応力（節点））
RESULTAT=resnonl), # 論理ユニット番号
UNITE=80)
解析の開始
通常通り、解析をスタートさせる。
計算結果の確認
計算が終了したので、結果を確認する。以下が確認した結果になる。 うまく計算できている。
接触面積が増加するモデルの場合
前記のモデルは、接触するものが四角柱のため、変形しても接触面積は変化しない。
このため、接触面のエッジにR面取りを施し、変形と共に接触面積が増加するモデルを作って解析してみる。
モデルの読み込み
モデルは、「multi-bar-1-R.stp」を読み込む。接触面のエッジをR面取りしたモデル。
解析は、~/CAE/contacr-R/のフォルダーを作り、この中で解析する。
Entityの作成
グループ化は、まったく同じように実施。
ただし、Barの接触面（contBar）は、接触する平面と変形と共にR面にも接触することになるので、平面＋R面をcontBarとしている。
| グループ名 | 備考 |
|---|---|
| Base | Solid1(Base) |
| fix | 固定面 |
| contBase | Baseの接触面 |
| Bar | Solid2(Bar) |
| contBar | Barの接触面 |
| press | 荷重を負荷する面 |
メッシュの作成
メッシュは、六面体だと、エラーが発生し、メッシュが切れなかったので、四面体の2次メッシュとし、Automatic Length＝0.2とした。次図参照。
解析コードの作成
同じ方法で作成。
解析開始
通常通り、解析をスタートさせる。 エラーが発生し、解析が終了する。
メッセージファイルを確認すると、ニュートン法の最大繰り返し回数で収束しなかったために停止している。
ニュートン法の最大繰り返し回数を増やして再計算を実行する。
resnonl = STAT_NON_LINE(CHAM_MATER=materfl,
CONTACT=contact,
CONVERGENCE=_F(ITER_GLOB_MAXI=30), # ニュートン法の最大繰り返し回数を追加（デフォルト値は10）
EXCIT=(_F(CHARGE=mecabc),
_F(CHARGE=mecach,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr),
MODELE=model,
NEWTON=_F(REAC_ITER=1))
結果の確認
最大応力は、3,400MPaでRの根元部で発生。
変形の様子を確認すると、R面に沿ってBase側が変形していることがわかる。（変形と共に接触面積が増える。 $\rightarrow$ 非線型になっている。）
接触面積が減少するモデルの場合
変形と共に接触面積が減少していくモデルを考える。下記のように変形と共にPlateの両端が持ち上がり、接触面積が減少していく場合を考える。
下記モデルで解析し、Plateの両端が持ち上がっている（接触面積が減少する）かどうかを確認する。
モデルの読み込み
モデルは、「test-contact-1.stp」を読み込む1。
~/CAE/contact-plate/というフォルダーを作り、この中で解析する。
Entityの作成
基本的には、前記したモデルと同じ。違いは、contBaseが2ヶ所ある（2平面をグループ化してcontBaseとした）ことと、press面はPlate中央の円柱の上面としていること。
| グループ名 | 備考 |
|---|---|
| Base | Solid1(Base) |
| fix | 固定面 |
| contBase | Baseの接触面 |
| Plate | Solid2(Plate) |
| contPL | Plateの接触面 |
| press | 荷重を負荷する面 |
メッシュの作成
四面体の2次メッシュとした。 細かさは、AutomaticLength＝0.1。
形状が複雑な分、メッシュが細かい。
解析コードの作成
コード自体は、ほとんど同じ。
press面の変位は、持ち上がりが確認できるように、大きな値（Y方向に-0.5mm）に設定。
計算
形状が複雑な分メッシュが多くなってしまっているので、計算時間は、長くかかる。
変位の結果を確認すると、Plateの両端が持ち上がっているのが確認できる。（下図参照。）
正しく接触判定をして、計算している。
ソースコード
DEBUT(LANG='EN')
mesh = LIRE_MAILLAGE(FORMAT='MED',
UNITE=20)
model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
PHENOMENE='MECANIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater = DEFI_MATERIAU(ELAS=_F(E=130300.0,
NU=0.343))
materfl = AFFE_MATERIAU(AFFE=_F(MATER=(mater, ),
TOUT='OUI'),
MODELE=model)
func = DEFI_FONCTION(NOM_PARA='INST',
VALE=(0.0, 0.0, 1.0, 1.0))
listr = DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=1.0,
PAS=0.2))
mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
DY=0.0,
DZ=0.0,
GROUP_MA=('fix', )),
_F(DX=0.0,
DY=0.0,
GROUP_MA=('press', ))),
MODELE=model)
mecach = AFFE_CHAR_MECA(DDL_IMPO=_F(DZ=-0.2,
GROUP_MA=('press', )),
MODELE=model)
contact = DEFI_CONTACT(FORMULATION='CONTINUE',
MODELE=model,
ZONE=_F(CONTACT_INIT='OUI',
GROUP_MA_ESCL=('contBase', ),
GROUP_MA_MAIT=('contBar', )))
resnonl = STAT_NON_LINE(CHAM_MATER=materfl,
CONTACT=contact,
EXCIT=(_F(CHARGE=mecabc,
FONC_MULT=func),
_F(CHARGE=mecach,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr),
MODELE=model,
NEWTON=_F(REAC_ITER=1))
resnonl = CALC_CHAMP(reuse=resnonl,
CONTRAINTE=('SIGM_ELNO', 'SIGM_NOEU'),
CRITERES=('SIEQ_ELNO', 'SIEQ_NOEU'),
MODELE=model,
RESULTAT=resnonl)
IMPR_RESU(FORMAT='MED',
RESU=_F(NOM_CHAM=('DEPL', 'SIEQ_NOEU', 'SIGM_NOEU'),
RESULTAT=resnonl),
UNITE=80)
FIN()
DEBUT(LANG='EN')
mesh = LIRE_MAILLAGE(FORMAT='MED',
UNITE=20)
model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
PHENOMENE='MECANIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater = DEFI_MATERIAU(ELAS=_F(E=130300.0,
NU=0.343))
materfl = AFFE_MATERIAU(AFFE=_F(MATER=(mater, ),
TOUT='OUI'),
MODELE=model)
func = DEFI_FONCTION(NOM_PARA='INST',
VALE=(0.0, 0.0, 1.0, 1.0))
listr = DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=1.0,
PAS=0.2))
mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
DY=0.0,
DZ=0.0,
GROUP_MA=('fix', )),
_F(DX=0.0,
DY=0.0,
GROUP_MA=('press', ))),
MODELE=model)
mecach = AFFE_CHAR_MECA(DDL_IMPO=_F(DZ=-0.2,
GROUP_MA=('press', )),
MODELE=model)
contact = DEFI_CONTACT(FORMULATION='CONTINUE',
MODELE=model,
ZONE=_F(CONTACT_INIT='OUI',
GROUP_MA_ESCL=('contBase', ),
GROUP_MA_MAIT=('contBar', )))
resnonl = STAT_NON_LINE(CHAM_MATER=materfl,
CONTACT=contact,
CONVERGENCE=_F(ITER_GLOB_MAXI=30),
EXCIT=(_F(CHARGE=mecabc),
_F(CHARGE=mecach,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr),
MODELE=model,
NEWTON=_F(REAC_ITER=1))
resnonl = CALC_CHAMP(reuse=resnonl,
CONTRAINTE=('SIGM_ELNO', 'SIGM_NOEU'),
CRITERES=('SIEQ_ELNO', 'SIEQ_NOEU'),
MODELE=model,
RESULTAT=resnonl)
IMPR_RESU(FORMAT='MED',
RESU=_F(NOM_CHAM=('DEPL', 'SIEQ_NOEU', 'SIGM_NOEU'),
RESULTAT=resnonl),
UNITE=80)
FIN()
DEBUT(LANG='EN')
mesh = LIRE_MAILLAGE(FORMAT='MED',
UNITE=20)
model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
PHENOMENE='MECANIQUE',
TOUT='OUI'),
MAILLAGE=mesh)
mater = DEFI_MATERIAU(ELAS=_F(E=130300.0,
NU=0.343))
materfl = AFFE_MATERIAU(AFFE=_F(MATER=(mater, ),
TOUT='OUI'),
MODELE=model)
func = DEFI_FONCTION(NOM_PARA='INST',
VALE=(0.0, 0.0, 1.0, 1.0))
listr = DEFI_LIST_REEL(DEBUT=0.0,
INTERVALLE=_F(JUSQU_A=1.0,
PAS=0.2))
mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
DY=0.0,
DZ=0.0,
GROUP_MA=('fix', )),
_F(DX=0.0,
DY=0.0,
GROUP_MA=('press', ))),
MODELE=model)
mecach = AFFE_CHAR_MECA(DDL_IMPO=_F(DZ=-0.2,
GROUP_MA=('press', )),
MODELE=model)
contact = DEFI_CONTACT(FORMULATION='CONTINUE',
MODELE=model,
ZONE=_F(CONTACT_INIT='OUI',
GROUP_MA_ESCL=('contBase', ),
GROUP_MA_MAIT=('contPL', )))
resnonl = STAT_NON_LINE(CHAM_MATER=materfl,
CONTACT=contact,
CONVERGENCE=_F(ITER_GLOB_MAXI=30),
EXCIT=(_F(CHARGE=mecabc),
_F(CHARGE=mecach,
FONC_MULT=func)),
INCREMENT=_F(LIST_INST=listr),
MODELE=model,
NEWTON=_F(REAC_ITER=1))
resnonl = CALC_CHAMP(reuse=resnonl,
CONTRAINTE=('SIGM_ELNO', 'SIGM_NOEU'),
CRITERES=('SIEQ_ELNO', 'SIEQ_NOEU'),
MODELE=model,
RESULTAT=resnonl)
IMPR_RESU(FORMAT='MED',
RESU=_F(NOM_CHAM=('DEPL', 'SIEQ_NOEU', 'SIGM_NOEU'),
RESULTAT=resnonl),
UNITE=80)
FIN()
-
モデルを読み込むと、モデルの上下方向がY軸方向となっているので、で、OX軸に対して90度回転させている。 ↩