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

この資料はオープンCAE勉強会＠岐阜で公開されているFS氏ご提供の「[SALOME-Mecaの使用法解説](http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?SALOME-Meca%A4%CE%BB%C8%CD%D1%CB%A1%B2%F2%C0%E2)：[6.0 接触解析の基本(1)](http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?plugin=attach&refer=SALOME-Meca%A4%CE%BB%C8%CD%D1%CB%A1%B2%F2%C0%E2&openfile=06-00.pdf)」を*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      | 荷重を負荷する面 |

![fig001.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/875ab843-8789-0a80-d4c6-533a39d48ce1.png)

## メッシュの作成

連結問題と同じ方法でメッシュを切る。
ツリーは、下記。
六面体の2次メッシュとし、ジオメトリからグループを作成する。

![fig002.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/5f15c16d-8bd3-55fe-492a-c233dd0c0668.png)

## 解析コードの作成

画面を*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）する。 下記参照。

```python:拘束条件
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に設定している。

```python:荷重条件
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*で追加する。この内容を以下で作成した。（ほとんどデフォルトのまま）

```python:接触条件
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のペアで入力する。

```python:増分関数
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倍）までを何分割して解析するのかを定義する。 下記参照。

```python:増分定義
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*の内容。

```python:非線形解析方法
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*で指定した名前を再利用する。

```python: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処理で再設定を行う。

```python:結果出力
IMPR_RESU(FORMAT='MED',              # 出力するバイナリ形式にMED形式を指定
          RESU=_F(NOM_CHAM=('DEPL', 'SIEQ_NOEU', 'SIGM_NOEU'), # DEPL（変位量）、SIEQ_NOEU（等価応力（節点））、SIGM_NOEU（応力（節点））
                  RESULTAT=resnonl), # 論理ユニット番号
          UNITE=80)
```

## 解析の開始

通常通り、解析をスタートさせる。

## 計算結果の確認

計算が終了したので、結果を確認する。以下が確認した結果になる。 うまく計算できている。

![fig003.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/585d2274-24fa-9bfd-7ea0-e2eed3d54192.png)

___

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

前記のモデルは、接触するものが四角柱のため、変形しても接触面積は変化しない。
このため、接触面のエッジに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      | 荷重を負荷する面 |

![fig004.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/c6e8dff0-6150-d25b-cfbd-03ea1ed4ee6e.png)


## メッシュの作成

メッシュは、六面体だと、エラーが発生し、メッシュが切れなかったので、四面体の2次メッシュとし、Automatic Length＝0.2とした。次図参照。

![fig005.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/8008942d-778a-541e-892e-35fcdd7ec792.png)

## 解析コードの作成

同じ方法で作成。

## 解析開始

通常通り、解析をスタートさせる。 エラーが発生し、解析が終了する。
メッセージファイルを確認すると、ニュートン法の最大繰り返し回数で収束しなかったために停止している。
ニュートン法の最大繰り返し回数を増やして再計算を実行する。

```python:非線形解析を修正
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$ 非線型になっている。）

![fig006.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/05c6e0be-530a-ab6d-195c-b3b454e6b97a.png "Z方向変位")
![fig007.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/a61c0ded-7add-711b-b2bd-cf57fbf32750.png "相当応力")

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

変形と共に接触面積が減少していくモデルを考える。下記のように変形と共に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      | 荷重を負荷する面 |

![fig008.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/92d3ebb0-72cf-57ce-ebf1-932b766435dd.png)

## メッシュの作成

四面体の2次メッシュとした。 細かさは、AutomaticLength＝0.1。
形状が複雑な分、メッシュが細かい。

![fig009.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/833ea426-e163-575d-21c3-0f200deffded.png)

## 解析コードの作成

コード自体は、ほとんど同じ。
press面の変位は、持ち上がりが確認できるように、大きな値（Y方向に-0.5mm）に設定。

## 計算

形状が複雑な分メッシュが多くなってしまっているので、計算時間は、長くかかる。
変位の結果を確認すると、Plateの両端が持ち上がっているのが確認できる。（下図参照。）
正しく接触判定をして、計算している。

![fig010.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/2e289386-3eda-0f97-743b-a825ebfa7932.png)![fig011.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/7de9b8b7-2459-6111-3e31-06b079ca8a16.png)

___

# ソースコード

```python: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()
```

```python: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()
```

```python: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版](https://github.com/JunTatsuno/LaTeX/raw/a5a30cddf1c197dbc06d4cca63aa3a818805e2d9/Contact%20Analysis%20Basics%20(Displacement%20Constraints)/How_to_use_Salome_Meca_6_0_Contact_Analysis_Basics__Displacement_Constraints.pdf)

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