この資料はオープン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コマンドが準備されるようになった。

境界条件の編集

境界条件は、
1. 通常の境界条件
2. 負荷を少しづつ変化させる条件
の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で指定した名前を再利用する。

Post処理

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」を読み込む¹。
~/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の両端が持ち上がっているのが確認できる。（下図参照。）
正しく接触判定をして、計算している。

ソースコード

multi-bar.commの場合

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()

multi-bar-R.commの場合

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()

contact-plate.comm（R面取りしたモデル）場合

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()

PDF版

モデルを読み込むと、モデルの上下方向がY軸方向となっているので、で、OX軸に対して90度回転させている。 ↩

Salome-Meca 2019を用いた接触解析の基本（1）（SSNV）

はじめに

モデルの読み込み

Entityの作成

メッシュの作成

解析コードの作成

解析コードの編集

境界条件の編集

通常の境界条件

少しづつ負荷させる境界条件作成

接触の定義

接触のためのコード追加

非線形解析方法の設定

Post処理の追加

結果出力の修正

解析の開始

計算結果の確認

接触面積が増加するモデルの場合

モデルの読み込み

Entityの作成

メッシュの作成

解析コードの作成

解析開始

結果の確認

接触面積が減少するモデルの場合

モデルの読み込み

Entityの作成

メッシュの作成

解析コードの作成

計算

ソースコード