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