Salome-Meca2019による丸棒のねじり解析

#解析対象の概要
※ ここでは，数値的解析精度検討の目的で指標にできる数量を現実問題をカバーするよう5桁で示しているが，実用上の目的に応じて適切に読み替えていただきたい．

+ 材料特性

ヤング率 *E* = 200 GPa
ポアソン比 *&nu;* = 0.3
とすると，
横弾性係数 *G* = *E*/{2(1+*&nu;*)} = 76.923 GPa

+ 丸棒軸の寸法

直径 *d* = 100 mm
軸長さ *L* = 500 mm
とすると，
断面二次極モーメント *I<sub>p</sub>* = *&pi;d*<sup>4</sup>/32 = 9.8175×10<sup>6</sup> mm<sup>4</sup>

+ 荷重条件
上記の丸棒の下端を固定し，上端に *T* = 1.0×10<sup>6</sup> N・mmのトルクを与える．


# 微小変形弾性理論解

変形が微小で，軸方向直線が変形後も直線形状を維持すると仮定すると，せん断ひずみ *&gamma;* に対するフックの法則と下図の幾何学的関係から，

端部のねじれ角 *&phi;* = *TL*/*GI<sub>p</sub>* = 6.6208×10<sup>-4</sup> rad (3.7934×10<sup>-2</sup>°)
上端軸表面の接線方向変位 *u* = (*d*/2)×*&phi;* = 3.3104×10<sup>-2</sup> mm
表面のねじり応力 *&tau;* = 16 *T*/*πd*<sup>3</sup> = 5.0929 MPa
純粋せん断応力下でのミーゼスの相当応力 *&sigma;* は，*&sigma;* = 3<sup>1/2</sup> *&tau;* = 8.9944 MPa
![roundbarmodel.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/608825/89b98560-b89f-cf2c-a35b-4539e7c86abe.jpeg)

#モデル化

+ 主に六面体ソリッド要素となるよう規則的に丸棒モデルを要素分割し，二次要素化する

参考：SALOME-Mecaの使用法解説：1.1 基本 メッシュ作成 
(http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?SALOME-Meca%A4%CE%BB%C8%CD%D1%CB%A1%B2%F2%C0%E2)

各種設定は，角柱のねじり解析に関する以下の記事を参考にバージョン違いを考慮して修正
岐阜高専 Salomeの使い方 トルク解析
 (http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?plugin=attach&refer=%C2%E8%A3%B1%A3%B4%B2%F3%CA%D9%B6%AF%B2%F1%A1%A7H240602&openfile=%A5%C8%A5%EB%A5%AF%B2%F2%C0%CF.pdf)

+ 設定したグループは下図（メッシュ上での位置とオブジェクトブラウザ）
+ 下端の面グループ fix を固定（DDL_IMPO）
+ 上端より100 mmだけ*z*軸方向の位置に独立節点 Moment を設定し，上端面の節点グループ load に含める
+ 節点グループload に対して剛体リンク（LIAISON_SOLIDE）を設定
+ 節点Momentに離散要素（0D Element）を設定し，Momentに*z*軸回りのモーメント *T* を与える（FORCE_NODALE, MZ=1000000.0）
+ 上端面で*x*軸に沿うエッジにグループlineXを設定し，応力分布をテキスト出力する（IMPR_RESU, FORMAT='RESULTAT'）
+ *x*= 50 mmとなる*z*軸方向の表面線にエッジグループlineZを設定し，変位分布をテキスト出力する（IMPR_RESU, FORMAT='RESULTAT'）
+ テキスト出力の際，節点座標をプリント（IMPR_COOR='OUI'）
+ その他の物理量はバイナリファイルに出力（IMPR_RESU, FORMAT='MED'）

![groups.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/608825/428611b9-1713-5378-444f-6f8e80be345b.jpeg)


#結果
変形図と変位のマグニチュードのカラーコンター図
上端表面変位のマグニチュードが軸端部表面接線方向変位の理論解  *u* = (*d*/2)×*&phi;* = 3.3104×10<sup>-2</sup> mm に近い．
※ ソリッド要素の変位ベクトルはデカルト座標系で与えられるため，変位は接線方向に向き，変形後の外周が円の外側にはみ出るように見える．
![TortionOfShaft0503.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/608825/5b4f26f9-6fb7-c44f-cf5a-4dd2a260bada.png)
#コマンドファイル

```
DEBUT(LANG='EN')

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

mesh0 = CREA_MAILLAGE(CREA_POI1=_F(GROUP_NO=('Moment', )),
                      MAILLAGE=mesh)

model = AFFE_MODELE(AFFE=(_F(MODELISATION=('3D', ),
                             PHENOMENE='MECANIQUE',
                             TOUT='OUI'),
                          _F(GROUP_MA=('Moment', ),
                             MODELISATION=('DIS_TR', ),
                             PHENOMENE='MECANIQUE')),
                    MAILLAGE=mesh0)

elemprop = AFFE_CARA_ELEM(DISCRET=_F(CARA='K_TR_D_N',
                                     GROUP_MA=('Moment', ),
                                     SYME='OUI',
                                     VALE=(1.0, 1.0, 1.0, 1.0, 1.0, 1.0)),
                          MODELE=model)

mater = DEFI_MATERIAU(ELAS=_F(E=200000.0,
                              NU=0.3))

fieldmat = AFFE_MATERIAU(AFFE=_F(MATER=(mater, ),
                                 TOUT='OUI'),
                         MODELE=model)

load = AFFE_CHAR_MECA(DDL_IMPO=_F(DX=0.0,
                                  DY=0.0,
                                  DZ=0.0,
                                  GROUP_MA=('fix', )),
                      FORCE_NODALE=_F(GROUP_NO=('Moment', ),
                                      MZ=1000000.0),
                      LIAISON_SOLIDE=_F(GROUP_NO=('load', )),
                      MODELE=model)

reslin = MECA_STATIQUE(CARA_ELEM=elemprop,
                       CHAM_MATER=fieldmat,
                       EXCIT=_F(CHARGE=load),
                       MODELE=model)

reslin = CALC_CHAMP(reuse=reslin,
                    CONTRAINTE=('SIGM_ELNO', 'SIGM_NOEU'),
                    CRITERES=('SIEQ_ELNO', 'SIEQ_NOEU'),
                    FORCE=('REAC_NODA', ),
                    RESULTAT=reslin)

IMPR_RESU(FORMAT='MED',
          RESU=_F(MAILLAGE=mesh0,
                  NOM_CHAM=('DEPL', 'REAC_NODA', 'SIGM_ELNO', 'SIGM_NOEU', 'SIEQ_ELNO', 'SIEQ_NOEU'),
                  RESULTAT=reslin),
          UNITE=8)

IMPR_RESU(FORMAT='RESULTAT',
          RESU=(_F(GROUP_NO=('Moment', ),
                   NOM_CHAM=('DEPL', ),
                   RESULTAT=reslin),
                _F(GROUP_MA=('lineX', ),
                   IMPR_COOR='OUI',
                   NOM_CHAM=('SIGM_NOEU', 'SIEQ_NOEU'),
                   RESULTAT=reslin),
                _F(GROUP_MA=('lineZ', ),
                   IMPR_COOR='OUI',
                   NOM_CHAM=('DEPL', ),
                   RESULTAT=reslin)),
          UNITE=2)

FIN()
```

# 評価
+ 数値解
上端面をリンクした節点（離散要素Moment）のねじり角 *&phi;* = 6.6215×10<sup>-4</sup> rad
座標(50, 0, 500)の節点（lineXの端部節点）上のねじり応力 *&tau;* = 5.1510×10 MPa

+ 理論解
端部のねじれ角 *&phi;* = *TL*/*GI<sub>p</sub>* = 6.6208×10<sup>-4</sup> rad (3.7934×10<sup>-2</sup>°)
表面のねじり応力 *&tau;* = 16 *T*/*πd*<sup>3</sup> = 5.0929 MPa

+ 理論解の仮定の確認
軸表面の直線（lineZ）上の接線方向変位（*y*方向変位）が軸方向距離に比例し，lineZが直線を維持することを確認．

![LineZ2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/608825/879d85c9-7360-5d6d-baeb-c585dbff0f78.jpeg)

本研究の一部は科学研究費補助金（18K02963，代表：藤岡照高）の助成を受けて実施した．