Salome-Meca 2020を用いた線形熱応力解析（SSLV）

この資料は[オープンCAE勉強会＠岐阜](http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?OpenCAE%20Users%20Wiki)で公開されている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)：5.0 線形熱応力、5.1 線形熱応力（2）」をSalome-Meca 2020を用いてモデル作成、メッシュ生成、解析条件設定、解析実行、結果処理を行いました。

# 1. はじめに

Salome-Mecaを使って、静的に熱応力を計算します。
温度は、均一に分布しているものとして、線形の弾性解析をします。
まずは、単純な細長い棒のモデルで、温度変化を与えて応力が計算どおりになっているかどうかを確認します。
その後、Bi-Metalのモデルを作り、温度変化による変位と応力を確認します。

# 2. 単純モデルの場合

100×20×10mmのモデルで、長手方向を拘束、20℃から120℃の温度変化を与えて、熱膨張分の応力がかかるようにする。
20℃のときが応力0の状態で、この状態から温度を120℃まで上昇させ熱応力をかけます。
~/CAE/thermal-bar/というフォルダーを作りこの中で解析します。

## 2-1. モデルの読み込み
「Bar-100x20x10.stp」のモデルを読み込みます。
このモデルは、100×20×10mmの四角柱となっています。

## 2-2. Entityの作成

モデルを拘束させる部位を定義（グループ化）します。
単純化する為に、片側端面固定、反対側端面の軸方向（X方向）のみ固定させます。
固定面を「fix0」、反対側端面を「fix100」でグループ化します。

![2021-05-29-08-38-04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/675abba2-2f2d-0cb5-d2b3-85606d9c5f7a.png)

## 2-3. メッシュの作成

Mesh画面に変え、「Mesh」「Create Mesh」を選択。
Create Mesh画面で「Assign a set of hypothees」ボタンをクリックして、「3D:Automatic Hexahedralization」（四角形のメッシュ）を選択。
線の分割数を指定する「Hypothesis Construction」画面が表示されますので、キャンセルします。
1Dタグを選択して、「Hypothesis」の歯車アイコンをクリックします。
メッシュの粗さは「Automatic Length」を選択し、「Fineness」は「0.2」を設定。

「Add. Hypothesis」の歯車アイコンをクリックして「Quadratic Mesh」（2次メッシュ）を選択。
（四角形のメッシュで2次メッシュにした。精度がよくなる為。）

「Apply and Close」を選択し、Create Mesh画面を閉じます。
Object Browseツリー上の「Mesh_1」を右クリックして、「Compute」でメッシュを作成します。
次に「Mesh_1」を右クリックして、「Create Groups from Geometry」でメッシュにグループを設定します。

![2021-05-29-08-55-58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/108c4c9f-fc83-f0db-5a76-07a4b462f305.png)

## 2-4. Aster Codeの作成

AsterStudy画面に変えて、「Isotropic liner elasticity」ウイザードでAster Codeを作成します。
fix0は全面固定の為、XYZ各方向0を入力、「＋」アイコンをクリックして、項目を追加してfix100に設定、fix100もXYZ各方向0を入力します。
圧力は、適当な面（fix0 or fix100）に適当な値を入力します。（後で圧力の境界条件は削除するので適当で可。）
結果ファイルは、~/CAE/thermal-bar-true/フォルダー内に「bar-true.rmed」として保存します。

## 2-5. Aster Codeの編集

編集内容は、

* 材料の定義（DEFI_MATERIAU）
  * ヤング率
  * ポアソン比
  * 線膨張係数
* 設定温度の定義（CREA_CHAMP）
  * 温度を120℃に設定
* 材料の設定（AFFE_MATERIAU）
  * 定義した材料を設定
  * 参照温度（20℃）を設定
  * 設定温度（120℃）を材料に関連つけ
* 境界条件（AFFE_CHAR_MECA）
  * fix0の設定：ウィザードで設定済み
  * fix100の設定：ウィザードで設定済み
  * 圧力の設定：ウィザードで設定したものを使用しない為、削除

となります。
AFFE_MATERIAUコマンド内ですべて定義し、温度の設定は、材料Fieldで定義します。

### 材料の定義

ここでは、方法を検証するためなので、線膨張係数の大きい下記の材料を選びました。

* 材料名：Aluminum
* ヤング率：70,600 MPa
* ポアソン比：0.345
* 線膨張係数：23.0e-6 1/℃

この値を入力します。
ツリーの構造は下記。

```python
Aluminum = DEFI_MATERIAU(ELAS=_F(ALPHA=2.3e-05,
                                 E=70600.0,
                                 NU=0.345))
```

### 設定温度の定義

新たに、ドロップダウンリストの「Post Processing」から「CALC_CHAMP」を作り、設定温度を定義します。
ツリーの構造は下記。

```python
tempS = CREA_CHAMP(AFFE=_F(NOM_CMP=('TEMP', ),
                           TOUT='OUI',
                           VALE=(120.0, )), # 120℃で定義
                   MODELE=model,
                   OPERATION='AFFE',
                   TYPE_CHAM='NOEU_TEMP_R') # 節点温度に設定
```

### 材料と温度の設定

定義した材料と設定温度をモデルに設定します。
ツリーの構造は下記。

```python
materfl = AFFE_MATERIAU(AFFE=_F(MATER=(Aluminum, ), # 材料の設定
                                TOUT='OUI'),
                        AFFE_VARC=_F(CHAM_GD=tempS, # このオペランドを追加し、温度Fieldを作ります、設定温度を関連つけ
                                     NOM_VARC='TEMP',
                                     TOUT='OUI',
                                     VALE_REF=20.0), # 参照温度を設定
                        MODELE=model)
```

### 境界条件の設定

まず、圧力の設定は、使わないのでコンセプト名mecachのAFFE_CHAR_MECAを削除します。
モデルの両端を固定する。ツリーの構造は、下記。

```python
mecabc = AFFE_CHAR_MECA(DDL_IMPO=_F(DX=0.0,
                                    DY=0.0,
                                    DZ=0.0,
                                    GROUP_MA=('fix0', 'fix100')),
                        MODELE=model)
```

コンセプト名mecachのAFFE_CHAR_MECAの削除後、コンセプト名resultのMECA_STATIQUEでエラーが発生するので、EXCITのCHARGE=mechachを削除します。
通常通り[^1]後処理の「CALC_CHAMP」の追加と、結果出力の「IMPR_RESU」を修正します。

## 2-6. 計算開始

通常通りに計算させます。

変位
![2021-05-29-10-24-13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/6958c4cd-89f7-78dc-5731-3a15f37c043b.png)

相当応力
![2021-05-29-10-24-26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/aab9793f-d3e2-8e19-c2c7-fd4c9402b71f.png)

最大主応力
![2021-05-29-10-24-41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/88e3f580-8242-546c-11ae-55764b4f8a73.png)

### 端面の拘束無しの場合

今回検証する為に、Code_Asterを編集して、fix100のX方向の拘束を無くして、変位のみを確認します。
この結果、端面コーナ部の変位は、

* X方向：0.2351
* Y方向：0.0115
* Z方向：0.0230

となります。

![2021-05-29-21-45-30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/4dac90fb-7f72-934e-3885-63216ad0b92f.png)

X方向の伸びは、100mm×100℃×23.0e-6＝0.230mm
となるはずであり、Salomeの答えは、0.2351mmなのでほぼ合致しています。
fix0を全面固定した為に、この部分で歪が発生しており、X方向の伸びが若干狂っていると思われます。
Y方向の伸びは、10mm×100℃×23.0e-6＝0.0230mmとなります。
Salomeの計算は、Y方向にまったく拘束していない為、コーナ部の伸びは片側の伸びを表しています。
したがって、SalomeのY方向の伸びは、0.0115mm×2＝0.0230mmとなり、計算どおりの答えになります。
同様にして、Z方向の伸びは、20mm×100℃×23.0e-6＝0.046mmとなるはずです。
Salomeの計算は、0.0230mm×2＝0.0460であり、計算通りの答え。
この為、温度の計算は、正しく行われています。

# 3. Bi-Metalのモデルの場合

線膨張係数の異なる2種類の金属を張り合わせたモデルを考え、このモデルを温度変化させたとき、変位や応力がどうなるか確認する問題を解きます。

## 3-1. モデルの読み込み

「circle2.stp」を読み込みます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/f2c086f5-d73e-8d2b-37ec-8f390801bd03.png)

このモデルは、円板が2枚あるモデルを中心で1/4にしたモデル。

* R40mm：t3mm
* R50mm：t3mm

モデルが対称なので、1/4のモデルで考えます。
解析は、~/CAE/thermal-circle/というフォルダーを作りこの中で解析します。

## 3-2. モデルの再構築

読み込んだモデルは、2ヶのSolidで構成されている為、いったん、Fuseコマンドで一体のモデルにした後、Partitionコマンドで、分割してメッシュを斬ります。
こうすることで、2ヶのSolidの境界面でで節点が共有できるモデルを作ることができます。
（部品を連結させても良いが、誤差が大きくなってしまうので、複合材として解析します。）
方法については、複合モデルの解析方法を参照。
ツリーの構造は、下記。

| ツリーの構造  | グループ名 | 備考                                                    |
|---------------|---------|---------------------------------------------------------|
| Circle2.stp_1 |         |                                                         |
|               | Lower   | 下側のSolid                                             |
|               | Upper   | 上側のSolid                                             |
|               | Sepa    | 分割面（Lowerの上面）を定義                             |
| Fuse_1        |         | LowerとUpperをFuseで一体化                              |
|               | *Lower  |                                                         |
|               | *Upper  |                                                         |
| Partition_1   |         | 一体化したFuse_1をSepa面で分割（Partition）             |
|               | *Fuse_1 |                                                         |
|               | *Sepa   |                                                         |
|               | Lower   | 下側のSolid（解析に使用するため再定義）                 |
|               | Upper   | 上側のSolid（解析に使用するため再定義）                 |
|               | fixY    | XZ平面でカットしたLower、Upperの面（Y方向に拘束する面） |
|               | fixX    | YZ平面でカットしたLower、Upperの面（X方向に拘束する面） |
|               | Cline   | Z軸に沿った円板（Lower、Upper）の中心線（XY方向に拘束） |
|               | Cpoint  | Lower下面の円板中心点（XYZ方向を拘束）                  |                                             

解析を1/4のモデルで実施する為、カットした面は、対称面として垂直方向に拘束する必要があります。
また、この他に移動・回転しないように、CLine、CPointを定義して、変形に影響を与えない方向で拘束します。
これら解析に使用するVolume、Face、Pointは、最終的なモデル（Partition_1）の下に上記のように再定義します。

![2021-05-29-11-48-39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/b10fee5d-4198-bee2-b5ee-af875698153b.png)

## 3-3. Meshの作成

画面をMeshに変えます。
メッシュは、「Automatic Tetrahedralization」でlengthは0.2、2次メッシュでメッシュを斬りました。
ツリー構造は、下記。

![2021-05-29-11-52-37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/946a7497-f493-c41b-c4dc-6a35d2de983a.png)

次に「Mesh_1」を右クリックして、「Create Groups from Geometry」でメッシュにグループを設定します。

![2021-05-29-12-21-40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/9fc46cf7-6b28-09ae-f67a-278ab6793b00.png)

## 3-4. Aster Codeの作成

AsterStudy画面に切り変えます。
ウイザードを使ってAster Codeを作ります。
材料や境界条件は、あとで編集する為、適当に入力しておいても問題ありません。
結果ファイルの保存は、~/CAE/thermal-circle/circle2.rmedとしました。

## 3-5. Aster Codeの編集

編集する項目は、

* 材料の定義
* 設定温度の定義（温度を120℃に設定）
* 材料の設定（参照温度を20℃に設定）
* 境界条件

となります。
最終的に、AsterStudyのツリーの構造は、下記となります。

```python
steel = DEFI_MATERIAU(ELAS=_F(ALPHA=1.18e-05, # steelの材料定義
                              E=212000.0,
                              NU=0.293))

Aluminum = DEFI_MATERIAU(ELAS=_F(ALPHA=2.3e-05, # Aluminumの材料定義
                                 E=70600.0,
                                 NU=0.345))

mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
                                     GROUP_MA=('fixX', )), # X方向を拘束
                                  _F(DY=0.0,
                                     GROUP_MA=('fixY', )), # Y方向を拘束
                                  _F(DX=0.0,
                                     DY=0.0,
                                     DZ=0.0,
                                     GROUP_NO=('CPoint', )), # 固定
                                  _F(DX=0.0,
                                     DY=0.0,
                                     GROUP_MA=('CLine', ))), # XY方向を拘束
                        MODELE=model)

tempS = CREA_CHAMP(AFFE=_F(NOM_CMP=('TEMP', ),
                           TOUT='OUI',
                           VALE=(120.0, )), # 温度を120℃に設定
                   MODELE=model,
                   OPERATION='AFFE',
                   TYPE_CHAM='NOEU_TEMP_R') # 節点温度

materfl = AFFE_MATERIAU(AFFE=(_F(GROUP_MA=('Lower', ), # Lowerにsteelをセット
                                 MATER=(steel, )),
                              _F(GROUP_MA=('Upper', ), # Upperにaluminumをセット
                                 MATER=(Aluminum, ))),
                        AFFE_VARC=_F(CHAM_GD=tempS,
                                     NOM_VARC='TEMP',
                                     TOUT='OUI',
                                     VALE_REF=20.0), # モデル全体の温度を20℃にセット（参照温度）
                        MODELE=model)
```

## 3-6. 計算開始

通常通りに計算開始します。

## 3-7. 結果の確認

変位を確認すると、温度上昇によって、下側へ反り返った状態になっています。aluminumのほうが伸びるので、下側に0.30mm反り返っています。
最大の相当応力は、120.6MPaと確認できます。

変位
![2021-05-29-13-06-03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/7930e43a-d4bc-84c9-4dfa-75c3b40b0e1e.png)

相当応力
![2021-05-29-13-06-26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/c113c6a7-23ce-c8a1-11a8-748790cdeb41.png)

最大主応力
![2021-05-29-13-07-07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/201403/04fd2be9-93eb-046a-7209-43ba671120ff.png)

# 4. Aster Codeの内容

以下に今回解析したCodeの内容を示します。

```python:bar-true.commの内容
DEBUT(LANG='FR')

mesh = LIRE_MAILLAGE(FORMAT='MED',
                     UNITE=20)

model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
                            PHENOMENE='MECANIQUE',
                            TOUT='OUI'),
                    MAILLAGE=mesh)

Aluminum = DEFI_MATERIAU(ELAS=_F(ALPHA=2.3e-05,
                                 E=70600.0,
                                 NU=0.345))

mecabc = AFFE_CHAR_MECA(DDL_IMPO=_F(DX=0.0,
                                    DY=0.0,
                                    DZ=0.0,
                                    GROUP_MA=('fix0', 'fix100')),
                        MODELE=model)

tempS = CREA_CHAMP(AFFE=_F(NOM_CMP=('TEMP', ),
                           TOUT='OUI',
                           VALE=(120.0, )),
                   MODELE=model,
                   OPERATION='AFFE',
                   TYPE_CHAM='NOEU_TEMP_R')

materfl = AFFE_MATERIAU(AFFE=_F(MATER=(Aluminum, ),
                                TOUT='OUI'),
                        AFFE_VARC=_F(CHAM_GD=tempS,
                                     NOM_VARC='TEMP',
                                     TOUT='OUI',
                                     VALE_REF=20.0),
                        MODELE=model)

result = MECA_STATIQUE(CHAM_MATER=materfl,
                       EXCIT=_F(CHARGE=mecabc),
                       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()
```

```python:circle2.commの内容
DEBUT(LANG='FR')

mesh = LIRE_MAILLAGE(FORMAT='MED',
                     UNITE=20)

model = AFFE_MODELE(AFFE=_F(MODELISATION=('3D', ),
                            PHENOMENE='MECANIQUE',
                            TOUT='OUI'),
                    MAILLAGE=mesh)

steel = DEFI_MATERIAU(ELAS=_F(ALPHA=1.18e-05,
                              E=212000.0,
                              NU=0.293))

Aluminum = DEFI_MATERIAU(ELAS=_F(ALPHA=2.3e-05,
                                 E=70600.0,
                                 NU=0.345))

mecabc = AFFE_CHAR_MECA(DDL_IMPO=(_F(DX=0.0,
                                     GROUP_MA=('fixX', )),
                                  _F(DY=0.0,
                                     GROUP_MA=('fixY', )),
                                  _F(DX=0.0,
                                     DY=0.0,
                                     DZ=0.0,
                                     GROUP_NO=('CPoint', )),
                                  _F(DX=0.0,
                                     DY=0.0,
                                     GROUP_MA=('CLine', ))),
                        MODELE=model)

tempS = CREA_CHAMP(AFFE=_F(NOM_CMP=('TEMP', ),
                           TOUT='OUI',
                           VALE=(120.0, )),
                   MODELE=model,
                   OPERATION='AFFE',
                   TYPE_CHAM='NOEU_TEMP_R')

materfl = AFFE_MATERIAU(AFFE=(_F(GROUP_MA=('Lower', ),
                                 MATER=(steel, )),
                              _F(GROUP_MA=('Upper', ),
                                 MATER=(Aluminum, ))),
                        AFFE_VARC=_F(CHAM_GD=tempS,
                                     NOM_VARC='TEMP',
                                     TOUT='OUI',
                                     VALE_REF=20.0),
                        MODELE=model)

result = MECA_STATIQUE(CHAM_MATER=materfl,
                       EXCIT=_F(CHARGE=mecabc),
                       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()
```

[^1]: [Salome-Meca 2020の使いかたの基本（SSLV）「解析条件の設定」](https://qiita.com/Jun_Tatsuno/items/f87220bbeb2ce0aaf3a7)