遠藤の卒論日誌 - *卒論日誌

目次

• 卒論日誌
• 研究内容
• ソリッド要素で螺旋折り円筒(2/1)
• ソリッド要素で線形解析、幾何学非線形有限変位微小変形(12/27)
• 土木学会東北支部原稿(1/20〜25)
• 再開(1/10)
• 中間発表まとめ中(12/16)
• メモ(12/9)
• 後藤メモ(2016/12/8)
• 進捗、有限変位、微小変形解析のためのCOQUE_3Dの使い方(12/5)
• 追加、変更分のcode
• 計算結果
• 今後の目標(11/11)、(11/18)、(11/28)、(12/2)
• モデル改良(11/7)
• 現在の状況(10/24)
• 現在の状況(10/7)
• LIAISON_UNIFを節点に適用
• シェル要素の中心応力の解析
• 螺旋折り円筒のモデリングのやり方(10/4)
• 螺旋折円筒の解析結果
• 現在の状況(9/12)
• 現在の状況(9/5)
• 現在の状況(8/23)
• 現在の状況(8/14)
• 現在の状況(8/7)
• 現在の状況(8/1)
• LIAISON_UNIFについて
• 反転螺旋折円筒のモデリング(途中)
• 次回までの目標
• 現在の状況(7/25)
• 現在の状況(7/21)
• salomeでシェル要素の3ヒンジアーチの解き方
• code-asterで非線形の解析をさせる方法
• 次の目標
• 目標(6/24)
• 線形解析と幾何学的非線形解析の違い(6/24)
• 当面の目標(6/13)
• 課題6 (6/27)
• 課題5 (6/20)
• 課題4 (6/13)
• 課題3 (6/6)
• 課題2 (5/23)
• 課題 (5/16)

卒論日誌

研究内容

 反転螺旋折り円筒のFEM解析（幾何学非線形を考慮した場合と、線形解析とでどれくらい差が出

 （反転）螺旋降り円筒の有限要素モデル化
 卒論タイトル、折り畳み可能な螺旋折り円筒の数値モデル

ソリッド要素で螺旋折り円筒(2/1)

ソリッド要素で線形解析、幾何学非線形有限変位微小変形(12/27)

 シェル要素でできなかったが立体要素なら可能(?)

 線形解析

 幾何学非線形有限変位微小変形

土木学会東北支部原稿(1/20〜25)

再開(1/10)

 不安定な変形のため計算が複雑になるのでどうしてもエラーがでる
 解決法としてRELATION='ELAS_HYPER',salomeバージョン10.5で使用できるらしい。
 使用コード(ELAS_HYPER不使用)

中間発表まとめ中(12/16)

メモ(12/9)

 エラスティカ問題、ニュートン反復

後藤メモ(2016/12/8)

進捗、有限変位、微小変形解析のためのCOQUE_3Dの使い方(12/5)

追加、変更分のcode

 メッシュは三角形二次要素にしてある

追加　MAILf=CREA_MAILLAGE(MAILLAGE=MAIL,
                   MODI_MAILLE=_F(TOUT='OUI',
                                  OPTION='TRIA6_7',),);


変更　MODE=AFFE_MODELE(MAILLAGE=MAILf,
           　　　      AFFE=_F(TOUT='OUI',
            　　　             PHENOMENE='MECANIQUE',
             　　　            MODELISATION='COQUE_3D',),);

変更　cara=AFFE_CARA_ELEM(MODELE=MODE,
        　　　            COQUE=_F(GROUP_MA=('Group_1','Group_4',),
          　　　                   EPAIS=0.01,
           　　　                  COQUE_NCOU=1,),);
〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜
追加　MODE2=AFFE_MODELE(MAILLAGE=MAIL,
             　　　    AFFE=_F(TOUT='OUI',
               　　　          PHENOMENE='MECANIQUE',
               　　　          MODELISATION='3D',),);

追加　PROJ=PROJ_CHAMP(RESULTAT=RESU,
             　　　   MODELE_1=MODE,
              　　　  MODELE_2=MODE2,);



変更　IMPR_RESU(FORMAT='MED',
       　　　   UNITE=80,
        　　　  RESU=_F(RESULTAT=PROJ,),);

計算結果

 右端の辺にz方向に変位を与えた時のx方向の変位、赤が最小で青が最大。

 まだ非線形解析なので正しい結果は出ていないが、DKT解析では出なかったx方向の変位が出ている。

今後の目標(11/11)、(11/18)、(11/28)、(12/2)

 ･折りたたみのアニメーション化
 ･モデル完成
 ･弧長増分法をつかってみる
 ･剛体に蝶番を入れいみる
 ･条件いろいろ変えてみる(荷重で計算、厚さ変更)
 ･非線形解析の大変形?
 ･解析方法を変えてみる(COQUE_3D)
 ･z方向の変位を与えて、x方向の変位を出す→COQUE_3Dを使い、有限(大)変位、微小変形

モデル改良(11/7)

 ヤング率の大きい面で円筒の上部を塞いでモデルの穴をなくす。
 ※DDL_IMPOとLIAISON_UNIFは同一面、線、点上に適用できない?グループ化に工夫が必要。

現在の状況(10/24)

 使用コードと使用モデルの各数値

DEBUT();

MA=DEFI_MATERIAU(ELAS=_F(E=70000.0,
                         NU=0.3,),);

MAIL=LIRE_MAILLAGE(FORMAT='MED',);

MODE=AFFE_MODELE(MAILLAGE=MAIL,
                 AFFE=_F(TOUT='OUI',
                         PHENOMENE='MECANIQUE',
                         MODELISATION='DKT',),);

cara=AFFE_CARA_ELEM(MODELE=MODE,
                    COQUE=_F(GROUP_MA=('Group_1','Group_2','Group_3','Group_4','Group_5','Group_6','Group_7','Group_8','Group_9','Group_10','Group_11','Group_12',),
                             EPAIS=0.01,),);

MATE=AFFE_MATERIAU(MAILLAGE=MAIL,
                   AFFE=_F(TOUT='OUI',
                           MATER=MA,),);

CHAR=AFFE_CHAR_MECA(MODELE=MODE,
                    DDL_IMPO=(
                        _F(GROUP_MA=('a','b','c','d','e','f',),
                           DZ=-60,),
                        _F(GROUP_NO=('61','62','63','64','65','66',),
                           DZ=-60,),
                        _F(GROUP_MA=('g',),
                           DX=0,
                           DY=0,
                           DZ=0,),
〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜
                        _F(GROUP_NO=('72',),
                           DX=0,
                           DY=0,
                           DZ=0,),),

                    LIAISON_UNIF=(
                        _F(GROUP_NO=('112','122',),
                           DDL=('DX','DY','DZ',),),
〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜
                        _F(GROUP_NO=('12110','12210',),
                           DDL=('DX','DY','DZ',),),
                           ),);

RESU=MECA_STATIQUE(MODELE=MODE,
                   CHAM_MATER=MATE,
                   CARA_ELEM=cara,
                   EXCIT=_F(CHARGE=CHAR,),);

RESU=CALC_ELEM(reuse =RESU,
               CHAM_MATER=MATE,
               RESULTAT=RESU,
               REPE_COQUE=_F(GROUP_MA=('Group_1','Group_2','Group_3','Group_4','Group_5','Group_6','Group_7','Group_8','Group_9','Group_10','Group_11','Group_12',),
                             NIVE_COUCHE='SUP',),
               OPTION=('SIGM_ELNO','SIEQ_ELNO',),);

RESU=CALC_NO(reuse =RESU,
             RESULTAT=RESU,
             OPTION=('SIGM_NOEU','SIEQ_NOEU',),);

IMPR_RESU(FORMAT='MED',
          UNITE=80,
          RESU=_F(MAILLAGE=MAIL,
                  RESULTAT=RESU,
                  NOM_CHAM=('SIGM_NOEU','SIEQ_NOEU','DEPL',),),);

FIN();

現在の状況(10/7)

LIAISON_UNIFを節点に適用

                    LIAISON_UNIF=(
                        _F(GROUP_NO=('11','12','13',),
                           DDL=('DX','DY','DZ',),),),

シェル要素の中心応力の解析

 SUPをMOYにすると中心部の解析になる。Group_1は解析する面グループ。

RESU=CALC_ELEM(reuse =RESU,
               CHAM_MATER=MATE,
               RESULTAT=RESU,
               REPE_COQUE=_F(GROUP_MA=('Group_1',),
                             NIVE_COUCHE='SUP',),
               OPTION=('SIGM_ELNO','SIEQ_ELNO',),);

螺旋折り円筒のモデリングのやり方(10/4)

 今まで去年の山内さんのモデリングのやり方を採用していたが、それはまちがいと判明。なのでこれは修正版。
 まず始めに下図の条件が螺旋折り円筒に必須。(n角形螺旋折り円筒)

 まず、平面に条件にあった図を作る。

 それを組み立てて(ここが難しい)

 完成、少しずれて完璧なものは作れなかった。

螺旋折円筒の解析結果

 間違っていたので削除。

現在の状況(9/12)

 螺旋折円筒モデル作成終了、code-asterも正常に動いた。

現在の状況(9/5)

 一辺10高さ10角度30の五角形螺旋折円筒モデル作成終了、あとは来週までにcode-asterを書いて計算、結果を考察する。

現在の状況(8/23)

現在の状況(8/14)

RESU=CALC_ELEM(reuse =RESU,
               CHAM_MATER=MATE,
               RESULTAT=RESU,
               REPE_COQUE=_F(GROUP_MA=('Group_32','Group_42',),
                             NIVE_COUCHE='SUP',),
               OPTION=('SIGM_ELNO','SIEQ_ELNO',),);

 蝶番の数は17個でやってみる。
 来週までの目標は、螺旋折り円筒を作ってみる。

現在の状況(8/7)

 先週の結果を踏まえて、蝶番をミシン目にしたらどうなるか調べようとしたが、LIAISON_UNIFを微小な直線に適用することが難しく、それまでしか出来なかった。

• 以下、エラー原因と改善法

 ･モデルは各々メッシュ作成しBuild Compoundでくっつける。
 ･メッシュを作成するとき2D:Automatic〜を使うとエラーが出るのでNetgen 1D-2Dを使う。Second Orderのチェックは外す。
 ･code-asterからMODI_MAILLAGE関連の文を消す。
 
 今月は勉強で忙しいので来週までの目標としてミシン目蝶番が使えるのか性能調査をする。

現在の状況(8/1)

LIAISON_UNIFについて

 ･隣り合うグループ同士に適用できない。
 ･点グループには適用できない。
 ･LIAISON_UNIFで選んだ物は完全に回転しなくなり、剛体になる。
 ･現状、これ以外で蝶番が作れるとは考えられない。

反転螺旋折円筒のモデリング(途中)

 前述の「LIAISON_UNIFについて」より普通の反転螺旋折円筒はすべての辺を蝶番にすることができないと分かったので、蝶番をミシン目にしてモデリングする。
 モデリング方法は去年の山内さんの卒論日誌を参考に作る。※間違いと判明(10/4のモデリングより)

次回までの目標

 ミシン目蝶番の応力と変位を調べる。結果に応じて、ミシン目蝶番を使って反転螺旋折円筒を作るか、他の方法を模索する。

現在の状況(7/25)

 非定常のものを解析するときは、力を加えるのではなく、変位を与えてでないと計算できない。
 これで、(7/21)まででわかったことと合わせれば、反転螺旋折円筒の非線形解析が出きるようになった（はず...）。

現在の状況(7/21)

salomeでシェル要素の3ヒンジアーチの解き方

• モデリング

 Rectangleで2つの平面を作り並べ、Compoundでひとつのモデルとし、全体面、固定する両端、2平面の隣り合うそれぞれの辺をGroupのCreateでグループ化する。

• mesh作成

 meshは、2D: Automatic Quadrangulationの1D: Local Lengthで好きな数値で切る。
 次に、できたMesh_1を右クリックしてcreate groups from geometryを選び、モデリングで作ったグループをmeshに関連づけさせる。

• code-aster

 まず、シェル要素を3Dで解析するには、解析方法として'DKT'を用いる。

MODE=AFFE_MODELE(MAILLAGE=MAIL,
                 AFFE=_F(TOUT='OUI',
                         PHENOMENE='MECANIQUE',
                         MODELISATION='DKT',),);

 これで2Dモデルを3D解析できる。(DKTの他にもCOQUE_3Dがあるが、これは変位や回転が大きい物に使う。)
 次に厚さを定義する。面全体グループの'Group_3'の厚さを0.001にするというコマンド。これを入れないと計算できない。

cara=AFFE_CARA_ELEM(MODELE=MODE,
                    COQUE=_F(GROUP_MA='Group_3',
                             EPAIS=0.001,),);

 この時、次のようにMECA_STATIQUEにCARA_ELEM=cara,も追加しないと計算失敗となる。

RESU=MECA_STATIQUE(MODELE=MODE,
                   CHAM_MATER=MATE,
                   CARA_ELEM=cara,
                   EXCIT=_F(CHARGE=CHAR,),);

 境界、荷重条件は、両端'Group_1'を固定し、2平面の隣り合うそれぞれの辺'Group_2'にFORCE_ARETEで線載荷する。

CHAR=AFFE_CHAR_MECA(MODELE=MODE,
                    DDL_IMPO=(
                        _F(GROUP_MA='Group_1',
                           DX=0,
                           DY=0,
                           DZ=0,
                           DRX=0,
                           DRY=0,
                           DRZ=0,),),

                    FORCE_ARETE=_F(GROUP_MA='Group_2',
                                   FZ=-1000.0,),
                    );

 そこにLIAISON_UNIFを追加して'Group_3'の2辺の変位だけを同じにする。

CHAR=AFFE_CHAR_MECA(MODELE=MODE,
                    DDL_IMPO=(
                        _F(GROUP_MA='Group_1',
                           DX=0,
                           DY=0,
                           DZ=0,
                           DRX=0,
                           DRY=0,
                           DRZ=0,),),

                    LIAISON_UNIF=(
                        _F(GROUP_MA='Group_3' ,
                           DDL=('DX','DY','DZ',),),
                          ),

                    FORCE_ARETE=_F(GROUP_MA='Group_3',
                                   FZ=-1000.0,),
                    );

 変位のみ同じで回転は継承していなかったので、蝶番になっているはず。計算結果も妥当なものになった。

code-asterで非線形の解析をさせる方法

 STAT_NON_LINEで非線形での解析ができる(らしい。エラーが出るので検証中、CARA_ELEM=caraのせい？)
 他に時間や降伏などの設定も必要。

sig400=DEFI_FONCTION(
                     NOM_PARA='EPSI',
                     VALE=(0.00114286,240000000,
                           0.0016,270000000,
                           0.0032,290000000,
                           0.0064,307000000,
                           0.0128,327000000,
                           0.02,341000000,
                           0.027,350000000,
                           0.04,360000000,
                           0.07,377000000,
                           0.1,390000000,
                           0.15,413000000,
                           0.2,430000000,
                           0.3,450000000,
                           ),
                     PROL_DROITE='LINEAIRE',
                     PROL_GAUCHE='LINEAIRE',);

MA=DEFI_MATERIAU(ELAS=_F(E=6.0,
                         NU=0.3,),
                 TRACTION=_F(SIGM=sig400,),);

〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜
inst=DEFI_LIST_REEL(DEBUT=0,
                    INTERVALLE=_F(JUSQU_A=10,
                                  NOMBRE=10,),);

fonc1=DEFI_FONCTION(NOM_PARA='INST',
                    VALE=(0,0,
                          10,0,),);

RESU=STAT_NON_LINE(MODELE=MODE,
                   CHAM_MATER=MATE,　　　　　　　　　　　　　 
　　　　　　　　　　　　　　　　　　　　　　　　　　←(「CARA_ELEM=cara,」を追加(7/25))
                   EXCIT=_F(CHARGE=CHAR,
                            FONC_MULT=fonc1,),
                   COMP_INCR=_F(RELATION='VMIS_ISOT_TRAC',),
                   INCREMENT=_F(LIST_INST=inst,),);

次の目標

 Asterで不静定を計算する。
 Asterで非線形での解析。

目標(6/24)

 板モデルにて蝶番を入れる

線形解析と幾何学的非線形解析の違い(6/24)

 ・線形解析、上下面対称の曲げ応力が発生。材料力学の公式と同じ計算結果になる。
 ・幾何学的非線形解析、曲げ応力+膜応力が発生。材料力学の公式と違う結果になる。(材料力学の公式と比べて、上下面での最大応力が小さくなる。)
 Asterで計算する時に STAT_NON_LINE と、DYNA_NON_LINE を使う。

当面の目標(6/13)

 salomeで要素と要素の間に蝶番入れられるか（境界条件、接続条件）
 3ヒンジアーチを参考にして考える

課題6 (6/27)

 l=100(mm)厚さ1(mm)縦横10(mm)のU型片持梁について先端にP=100(N)を与える、応力分布をsalomeを使いグラフ化。
 拘束面

 一要素隣

 梁の中央

課題5 (6/20)

 課題4から、側面をy方向に拘束する

 一要素隣の場合、より線形に近かった(梁の固定端では数値が正確に出ない？)

 梁の真ん中の場合

課題4 (6/13)

 B=10(mm) H=10(mm) L=200(mm) の固定端片持ち梁にて先端にP=100(N)が加えられるとき、　meshを立方体要素、length=2(mm) とした時salomeとgnuplotを使い、
 1.固定端での応力を縦軸(z)にしてグラフ化
 2.断面(z=0)を3次元プロットに一緒にプロットする
 3.直方体要素で、固定端から1要素となりの断面の応力分布は、どれくらい変わるか
 4.梁の真ん中の断面の応力分布はどうか（以上の3箇所を比較）
 5.手計算σ=MIyとどれくらい合うか。手計算で求まる面を3次元プロットにあわせてプロット。
 
 1.2一緒にプロット

 1.3.4の比較

 5
 固定端

 固定端の一要素隣

 梁の中央

課題3 (6/6)

 B=10(mm) H=10(mm) L=200(mm) の片持ち梁にて先端にP=100(N)が加えられるとき、　meshを立方体要素、length=2(mm) とした時固定端での応力分布をgnuplotを使いグラフ化。

課題2 (5/23)

  初等梁　v=Pl^3/3EI=0.0533333
 Timoshenko梁　v=Pl^3/3EI+Pl/kGA=0.0536533
 B=10(mm) H=10(mm) L=200(mm) 

四面体要素、線形

四面体要素、二次要素

立方体要素、線形

立方体要素、二次要素

課題 (5/16)

 初等梁　v=Pl^3/3EI=0.0533333
 Timoshenko梁　v=Pl^3/3EI+Pl/kGA=0.0536533
 B=10(mm) H=10(mm) L=200(mm)