振動解析における回転系座標についての検証を行う。

片持ち梁

X=100,Y=10,Z=20の片持ち梁を作成。画像奥側が固定面。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti1.png

等方性

細かい条件については、根本の卒論日誌:片持ち梁の振動解析(ソリッド要素)と同じ。

振動モード解析値(Hz)
水平一次20.8188
鉛直一次39.9531
水平二次124.4
ねじれ145.725
鉛直二次213.841

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti2.png

異方性

等方性から、材料条件のみを変更。

ihou = DEFI_MATERIAU(
 ELAS_ORTH=_F(
   E_L=6000.0, 
   E_N=240.0, 
   E_T=240.0, 
   G_LN=400.0, 
   G_LT=400.0, 
   G_TN=400.0, 
   NU_LN=0.016, 
   NU_LT=0.016, 
   NU_TN=0.016, 
   RHO=3.8e-07
 )
振動モード解析値(Hz)
水平一次19.7614
鉛直一次35.8196
ねじれ64.7791
水平二次104.701
鉛直二次148.822

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti3.png

異方性+回転+回転座標適用

梁を回転させて、かつ回転座標を適用させる。

材料設定について、上のELAS_ORTHをそのまま適用。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti7.png

回転座標については

AFFE_CARA_ELEM > MASSIF > ANGL_REP
入力欄を3つ作成し、それぞれがX,Y,Zに対応している。
それぞれの軸回りでいくら回転させるかを入力。
ex)今回の場合、Y軸を90度回転させるので、1番目に0、2番目に90、3番目に0を入力

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti5.png

振動モード解析値(Hz)
水平一次19.7614
鉛直一次35.8196
ねじれ64.7791
水平二次104.701
鉛直二次148.822

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/katamoti6.png

回転座標を適用する前と同じ結果になったため、ANGL_REPで回転座標を適用できるといえる。

結果

片持ち梁の場合は、ANGL_REPでの回転座標適用は有効だと考えられる。

単純梁

次に、同じくX=100,Y=10,Z=20で単純梁を作成。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun1.png

固定方法については、Salome-Meca演習_単純梁のジオメトリの作成と同じ。

等方性

材料設定は片持ち梁と同じ。

振動モード解析値(Hz)
水平一次128.488
鉛直一次157.321
ねじれ233.374
水平二次349.164
ねじれ447.339

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun2.png

固定箇所が回転しているように見えるが、確認したところ回転していなかった。

異方性

材料設定は片持ち梁と同じ。

振動モード解析値(Hz)
ねじれ75.3159
水平一次88.7695
鉛直一次110.323
ねじれ163.276
水平二次200.492

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun3.png

異方性+回転+回転座標適用

片持ち梁と同様に梁を回転させ、ANGL_REPを適用させる。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun4.png

振動モード解析値(Hz)
ねじれ75.7328
水平一次89.0501
鉛直一次110.629
ねじれ163.634
水平二次201.511

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun5.png

回転前と少しだけ振動数が違う。

回転前解析値(Hz)回転後解析値(Hz)
75.315975.7328
88.769589.0501
110.323110.629
163.276163.634
200.492201.511

ANGL_REPが正しく適用されていない可能性もあるため、比較用に回転前の梁(X軸方向が長軸)に回転座標を適用させてみる。

【比較用】異方性+回転座標で梁を回転させなかった場合

振動モード解析値(Hz)
水平一次30.0518
鉛直一次33.695
水平二次82.913
鉛直二次115.306
ねじれ120.604

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/tanjyun6.png

異方性を持たせて回転なしのときと、異方性+回転+回転座標適用のときとも違う結果になったため、 これを見る限り、回転座標は適用されているようだ。

結果

単純梁の場合は、ANGL_REPでの回転座標適用は有効だと考えられるが、 回転前と後で少しだけ振動数が違った。

考えられる原因としては、単純梁は線要素での固定であるため、片持ち梁のような面要素での固定よりも 固定部分の周りが動きやすいことから、

その違いが振動数となって現れた可能性は考えられる。

アーチ

アーチへの適用も検証する。 図の3点を使用してアーチを作成。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch1.png

等方性

材料設定は以下のとおり。

mater = DEFI_MATERIAU(
 ELAS=_F(
   E=7500.0, 
   NU=0.4, 
   RHO=3.8e-10
 )
振動モード解析値(Hz)
水平一次4.1721
鉛直二次9.44272
水平二次11.6283
対称一次20.8357
水平三次23.306
対称二次38.3033
水平四次39.7125
鉛直四次41.8928
水平五次60.2621
対称三次65.093

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch2.png

【比較用】異方性

異方性アーチの場合はELAS_ISTRを使用するのがいいらしい。(及川さん談)

ihou = DEFI_MATERIAU(
 ELAS_ISTR=_F(
   E_L=7500.0, 
   E_N=300.0, 
   G_LN=400.0, 
   NU_LN=0.016, 
   NU_LT=0.016, 
   RHO=3.8e-10
 )
振動モード解析値(Hz)
水平一次3.01718
鉛直二次5.94478
水平二次8.77324
対称一次13.5595
水平三次17.4116
対称二次22.7029
水平四次27.1765
鉛直四次30.9601
水平五次40.2001
対称三次46.4758

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch3.png

アーチの場合下記の回転座標適用をしないと、ただの直交異方性が適用されてしまうので、参考値とする。

異方性+回転座標適用

アーチの回転座標適用にはANGL_AXEとORIG_AXEコマンドを使用する。 (上のANGL_REPと同じ位置にある)

ANGL_AXEは回転角(図を参照)、ORIG_AXEはアーチの中心座標を指定する。

elemprop = AFFE_CARA_ELEM(identifier='2:1',
                         MASSIF=_F(ANGL_AXE=(0.0, 90.0),
                                   GROUP_MA=('Cut_1', ),
                                   ORIG_AXE=(0.0, 10000.0, 14500.0)),
                         MODELE=model)
振動モード解析値(Hz)
水平一次2.74942
鉛直二次5.93091
水平二次8.07844
対称一次13.4888
水平三次16.4706
対称二次22.104
水平四次27.0818
鉛直四次29.4049
水平五次39.9771
対称三次44.0591

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch4.png

上の異方性のみの材料設定よりも、固定端近くの剛性が高いように見える。

【比較用】回転座標の適用について

ORIG_AXEでアーチの中心座標を指定するが、ANGL_AXEで回転角を調整する前の中心座標か調整後の中心座標かわからなかったので、検証用として回転角を調整前の中心座標(X=14500,Z=0)を入力して解析する。

elemprop = AFFE_CARA_ELEM(identifier='2:1',
                         MASSIF=_F(ANGL_AXE=(0.0, 90.0),
                                   GROUP_MA=('Cut_1', ),
                                   ORIG_AXE=(14500.0, 10000.0, 0.0)),
                         MODELE=model)
振動モード解析値(Hz)
水平一次2.98803
鉛直二次7.66834
水平二次8.77053
水平三次16.3103
対称一次18.0064
水平四次27.5536
対称二次29.5963
鉛直四次35.3419
水平五次41.2204
対称三次54.2841

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch5.png

他の結果よりも梁が横に伸びたように見える。

ELAS_ORTHを使用した場合

ihou2 = DEFI_MATERIAU(identifier='5:1',
                     ELAS_ORTH=_F(E_L=7500.0,
                                  E_N=300.0,
                                  E_T=300.0,
                                  G_LN=400.0,
                                  G_LT=400.0,
                                  G_TN=400.0,
                                  NU_LN=0.016,
                                  NU_LT=0.016,
                                  NU_TN=0.016,
                                  RHO=3.8e-10))
振動モード解析値(Hz)
水平一次1.91278
鉛直二次5.51539
水平二次5.67227
鉛直三次11.627
対称一次12.6893
水平三次19.7471
21.9378
鉛直四次24.8051
水平五次29.5828
対称二次37.1475

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch6.png

ELAS_ISTRと全く違う値になった。 回転座標は適用させているので、及川さんのアーチ橋に異方性を適用させたら振動数が変わったのはこれが原因かも?

まとめ

等方性と異方性について比較すると、異方性の方が小さい値の振動数になる傾向があることがわかった。

また、アーチの異方性については、ELAS_ISTRを使用し、ANGL_AXEで回転角を、ORIG_AXEで回転の原点を指定する。

11/22追記

code_asterマニュアル「Operator DEFI_MATERIAU」P21に、以下の記載がある。

3.6 Keywords factor ELAS_ISTR, ELAS_ISTR_FO, ELAS_VISCO_ISTR
Definition of the constant elastic characteristics (ELAS_ISTR) or functions of the temperature in the
case of the transverse isotropy for the isoparametric solid elements (ELAS_ISTR_FO). It is possible to
define   constant   transverse   isotropic   characteristics   viscoelastic   (ELAS_VISCO_ISTR)  at  a  given
frequency, the modules to be informed are then complex.
By taking again the same notations as for the orthotropism, the transverse isotropy means here that
isotropy is in the plan  (L,T)  and that the orthotropism is thus carried by the direction  N  [R4.01.02].
One can draw the attention of  the reader to the fact that this convention differs from  a usual
convention   which   indicates   by   “longitudinal   direction”   the   direction   of   orthotropism   of   isotropic
transverse materials.

異方性のみ(回転なし)

振動モード解析値(Hz)
水平一次1.7197
鉛直二次5.24329
対称一次5.49752
鉛直三次10.7901
逆対称一次?12.6282
鉛直四次18.3927
水平四次22.1162
鉛直四次24.6472
鉛直五次27.604
対称二次36.8846

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch7.png

異方性+回転座標

上の結果にANGL_AXEとORIG_AXEを適用させる。 おそらくこれが正しい(マニュアル上は)アーチへの異方性の適用のさせ方。

振動モード解析値(Hz)
水平一次1.11295
鉛直二次3.01786
対称一次3.3733
鉛直三次6.01338
逆対称一次?6.84415
対称二次9.46839
鉛直四次10.0918
水平四次15.0331
鉛直五次15.4198
鉛直六次20.5621

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch8.png

【比較用】異方性(ELAS_ORTH)+回転座標

比較用に、ELAS_ORTHについても修正したものを使用して解析する。

ihou4 = DEFI_MATERIAU(
 ELAS_ORTH=_F(
   E_L=300.0, 
   E_N=7500.0, 
   E_T=300.0, 
   G_LN=400.0, 
   G_LT=400.0, 
   G_TN=400.0, 
   NU_LN=0.016, 
   NU_LT=0.016, 
   NU_TN=0.016, 
   RHO=3.8e-10
 )
振動モード解析値(Hz)
水平一次1.31775
対称一次3.43991
鉛直二次3.59629
逆対称一次?6.90519
鉛直三次7.11015
対称二次9.55552
鉛直四次12.1069
対称三次15.4075
鉛直五次18.4823
水平五次24.0486

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/arch9.png

まとめ

以下に等方性アーチの結果と、(マニュアル上はおそらく正しい)異方性アーチの結果についてまとめる。

振動モード解析値(Hz)
水平一次4.1721
鉛直二次9.44272
水平二次11.6283
対称一次20.8357
水平三次23.306
対称二次38.3033
水平四次39.7125
鉛直四次41.8928
水平五次60.2621
対称三次65.093
振動モード解析値(Hz)
水平一次1.11295
鉛直二次3.01786
対称一次3.3733
鉛直三次6.01338
逆対称一次?6.84415
対称二次9.46839
鉛直四次10.0918
水平四次15.0331
鉛直五次15.4198
鉛直六次20.5621

比較すると、全体的に1/3程度の振動数になっている。

斜めの梁

(10,10,100)でジオメトリを作成後、Y軸を軸として45°回転させる。

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/naname1.png

これを片持ち梁として解析→伸びを手計算と比較し、グローバル座標系/ローカル座標系について考察を行う。

解析

DX変位

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/naname3.png

DY変位

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/naname4.png

ほとんど変化はない。

DZ変位

http://www.str.ce.akita-u.ac.jp/~gotouhan/ishikuro/kaiten/naname5.png

計算結果との比較

解析値の合力は\( \frac{DX+DZ}{2}×\frac{1}{\sqrt{2}} \)で計算。

理論値解析値(DX伸び)解析値(DZ伸び)解析値(合力)相対誤差(%)
0.01333330.01872440.01869840.0132310-0.00076725

理論値と解析値の相対誤差はかなり小さいため、解析自体に問題はないと思われるが、変位、各種応力についてはX,Y,Z座標ごとにしか表示させることはできなかった。

【比較用】斜めじゃない梁

斜めにしない(Z軸に並行な)梁について比較する。条件は同じ。

理論値解析値相対誤差(%)
0.01333330.0132375-0.0007185

ひずみ、応力の比較

ひずみ、応力についても比較する。

EPSI_NOEU(XX)EPSI_NOEU(YY)EPSI_NOEU(ZZ)EPSI_NOEU(XY)EPSI_NOEU(XZ)EPSI_NOEU(YZ)
斜めの梁8E-05-0.0001066678E-05-1.14879E-180.000186667-4.524E-18
並行な梁-5.33333E-05-5.33333E-050.0001333331.43741E-18-3.79471E-19-7.54424E-19
SIGM_NOEU(XX)SIGM_NOEU(YY)SIGM_NOEU(ZZ)SIGM_NOEU(XY)SIGM_NOEU(XZ)SIGM_NOEU(YZ)
斜めの梁1-5.15282E-141-6.15422E-151-2.42357E-14
並行な梁1.88738E-156.10623E-1517.70044E-15-2.03288E-15-4.04156E-15

同じ結果ではないため、斜めの梁でも並行な梁でも同じ座標系が適用されているわけではないようだ。

ANGL_REPとANGL_EULER

ANGL_REPとANGL_EULERの場合でY軸を軸に45°回転を適用させて、座標系が変化しているか見る。

①elemprop = AFFE_CARA_ELEM(identifier='2:1',
                         MASSIF=_F(ANGL_REP=(0.0, 45.0, 0.0),
                                   GROUP_MA=('naname', )),
                         MODELE=model)
②elemprop = AFFE_CARA_ELEM(identifier='2:1',
                         MASSIF=_F(ANGL_EULER=(0.0, 45.0, 0.0),
                                   GROUP_MA=('naname', )),
                         MODELE=model)

①と②と上の結果を比較する。

EPSI_NOEU(XX)EPSI_NOEU(YY)EPSI_NOEU(ZZ)EPSI_NOEU(XY)EPSI_NOEU(XZ)EPSI_NOEU(YZ)
ANGL_REP8E-05-0.0001066678E-05-1.14879E-180.000186667-4.524E-18
ANGL_EULER8E-05-0.0001066678E-05-1.14879E-180.000186667-4.524E-18
SIGM_NOEU(XX)SIGM_NOEU(YY)SIGM_NOEU(ZZ)SIGM_NOEU(XY)SIGM_NOEU(XZ)SIGM_NOEU(YZ)
ANGL_REP1-5.15282E-141-6.15422E-151-2.42357E-14
ANGL_EULER1-5.15282E-141-6.15422E-151-2.42357E-14

変化がないので、座標系の適用方法が間違っているかもしれない。

メモ

salomeメモの「局所座標系」(近藤さん編集?)を確認すると、以下の記載がある。

                    MASSIF=(_F(GROUP_MA='moto',               ←回転なしの梁に適用
                                ANGL_EULER=(0.0,0.0,0.0,),     #(\psi,\theta,\phi)
                                ),
                             _F(GROUP_MA='alp330',       ←θx=330°傾けた梁に適用
                                ANGL_EULER=(0.0,330.0,270.0,), #(\psi,\theta,\phi)
                                ),
                             _F(GROUP_MA='bt30',               ←θy=30°傾けた梁に適用
                                ANGL_EULER=(90.0,30.0,270.0,), #(\psi,\theta,\phi)
                                ),
                             ),
                     );

270.0の数字がよくわからない。

3つ目の数字はθzの回転? だとしたら、オイラー角の適用のさせ方か考え方が間違っている可能性があるかも。

'21.12.02 近藤より

ここでは,\( E_{xx}=E_{yy}=\frac{E_{zz}}{25} \)としているので,\( \phi \)の値は何であっても結果は変わらないと思いますが,より一般化するために後藤さんの軸のとり方に合わせて載荷方向がy軸の正となるように回転させました(うろ覚え). ちなみにANGL_EULERは若干ややこしい(下記)ので,このようなサイトでチェックするのがいいかもです.

追記

まさかの近藤さんからコメントが・・・ありがとうございます。

もう一度マニュアルについて確認したところ、近藤さんのコメントの通りZ→X→Z'の順番の回転でした。 (ずっとオイラー角はx→y→zの順だと思ってました)

ななめの梁を異方性にして等方性と比較

まず、オイラー角の設定を以下に修正。

elemprop = AFFE_CARA_ELEM(identifier='2:1',
                         MASSIF=_F(ANGL_EULER=(90.0, 45.0, 0.0),
                                   GROUP_MA=('naname', )),
                         MODELE=model)

これで斜めの梁の中心にZ軸が通る(はず)

材料設定を2種類作る。

mater = DEFI_MATERIAU(identifier='3:1',
                     ELAS=_F(E=7500.0,  #等方性
                             NU=0.4))
mater = DEFI_MATERIAU(identifier='3:1',
                     ELAS_ORTH=_F(E_L=300.0, #異方性
                                  E_N=7500.0,
                                  E_T=300.0,
                                  G_LN=400.0,
                                  G_LT=400.0,
                                  G_TN=400.0,
                                  NU_LN=0.016,
                                  NU_LT=0.016,
                                  NU_TN=0.016))

E_Nを7500にして、ANGL_EULERの設定が正しければ材料の強軸方向ヤング率が等方性と異方性で同じなので、先端の変位、応力はほぼ同じになるはず。

EPSI_NOEU(XX)EPSI_NOEU(YY)EPSI_NOEU(ZZ)EPSI_NOEU(XY)EPSI_NOEU(XZ)EPSI_NOEU(YZ)
等方性8E-05-0.0001066678E-05-1.14879E-180.000186667-4.524E-18
異方性8E-05-0.0001066678E-052.28756E-140.000186667-2.56527E-13
SIGM_NOEU(XX)SIGM_NOEU(YY)SIGM_NOEU(ZZ)SIGM_NOEU(XY)SIGM_NOEU(XZ)SIGM_NOEU(YZ)
等方性1-5.15282E-141-6.15422E-151-2.42357E-14
異方性1-1.76453E-111-1.83004E-111-2.05221E-10

等方性、異方性でほぼ同じ結果が出たため、正しく適用できたようだ。

気になること

アーチ(変形)

アーチの異方性をsalome_meca2019のAFFE_CARA_ELEM > MASSIF > ANGL_AXEを用いて 回転対称性、横方向等方性の軸指定でやってみる。

local axis


トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-07-20 (木) 13:35:19