静解析CAE屋さんのための動解析FEM(有限要素法)入門

CAEで静解析応力解析しかやったことがない方のための動解析講座。怪しいところ有ったらご連絡ください(^^;)

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
  1. --/--/--(--) --:--:--|
  2. スポンサー広告

モーダル法過渡応答の原理の基本その2

前回

(ある時間の変形)=A×(1次の固有モード形状)+B×(2次の固有モード形状)+C×(3時の固有モード形状)+…(無限に続く)

と書きましたが、

以下自問自答コーナー(^^)

「前、固有モード形状は形状のみ(変位の相対値)意味があると言ったのに、それでは足し算とか掛け算とかできないのでないのですか?

「そのとおりです。

http://freecaetester.blog62.fc2.com/blog-entry-629.html

にもかいてあります。

よって、ある値を基準に「変位」を決めています。

基準をとる方法は、実はいろいろあります。

たとえば変形の一番大きいところを1mmとして、残りはその一番大きいところとの比から計算する、ということもできます。

実際は、

モード質量を正規化(「1」)

になるように決める事が多いです。。

モード質量が1とは何かは、また機会があれば述べましょう(^^)。」


では次の質問

「そもそも本当に固有モードを足し合わせれば振動は表現できるか?」

「できます。直感的に理解する方法としては、現実の構造物の固有モードは無限にあるからです。

無限に形状があるということは、いつかは表せるということですよね(^^)。」


「でも有限要素法では固有値の数は有限だし、さらに実際計算するときは荷重振動数の3倍で十分とか20個で十分とか言って、全部計算しないではないですか?」

「そこがまさにポイントです!。まず、一般的に低い振動数のほうが、高い振動数より影響が小さいのです。なぜなら、高い振動数は同じ振幅だと速度が速くなるので減衰しやすい。よって減衰してしまうか、そもそも振幅が小さいから無視できると場合が多い、と言うわけです。」


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
スポンサーサイト
  1. 2012/06/03(日) 22:56:15|
  2. 未分類
  3. | トラックバック:0
  4. | コメント:0

果たしてモード形状の足し合わせでどんな形状も表現できるのか

前回まで、各固有モード形状(正規化したもの)を何倍かしたものを足し合わせると、過渡解析のある時間の結果が得られると書きました。

しかし、モード形状の足し合わせだけで本当に表すことができるのでしょうか…

理論云々も重要ですが、それは専門書に譲って、実際のモデルで確かめてみましょう。

簡単のため、梁要素で説明しようと思ったのですが、

CalculiXの梁要素は少し特殊で扱いにくいので、Femap with NX Nastranを使います。

モデルはこんな感じ…

ds465

節点番号は左から1,2,3,4とします。


片端は完全固定、あとの節点はy方向のみ動くとします。

ちなみにNastranバルクデータを作成してから、それをFemapに読み込ませています。

バルクデータは以下のとおりです。

sol 103
TIME 10
CEND
SPC=1
METHOD=10
DISP=ALL
BEGIN BULK
EIGARL,10,,,5
MAT1,1,210E9,,0.3,7860
PBAR,100,1,2.5e-5,5.20833e-11,5.20833e-11
CBAR,1,100,1,2,0,1,0
CBAR,2,100,2,3,0,1,0
CBAR,3,100,3,4,0,1,0
GRID,1,,0,0,0
GRID,2,,0.25,0,0
GRID,3,,0.50,0,0
GRID,4,,0.75,0,0

SPC1,1,123456,1
SPC1,1,13456,2,3,4

CEND

もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/04(月) 22:56:20|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

3自由度モデルの固有モード形状を計算する

前回のモデルを固有値解析して見ます。

モード形状は

モード1(34.1Hz)
ds466


正規化されたy方向の変位は、左側の節点から

{0, -1.84, -3.19, -3.68}

になります。

モード2(93.1Hz)
ds467


正規化された変位は

{0.0,3.68,0.0,-3.68}


モード3127.1Hz
ds468


{0.0,-1.84,3.19,-3.68}


入力データでは、5個固有値を求めるように入力しているのですが、

モデルの計算する自由度が3つ(節点2,3,4のy方向)だけなので

固有値も3つしかありません。

この4節点のモデルはこの3つの形状を何倍かすれば使えばどんな形状で表現することができるのです。。

もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/05(火) 23:09:55|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

固有ベクトルとは何ぞや

前回、求めた固有モード形状を使って説明します。

数学的にベクトルを使って表現します。

ベクトルといっても、

変位を縦に並べただけです

ds469

拘束点は変位が0なので省略して

ds470


モード1をベクトル方式であらわすと
ds471

モード2は
ds472

モード3は
ds473


これらを固有ベクトルと言います。。。

この固有ベクトルがモーダル法過渡解析を解く鍵なのです。

前に各固有モード形状(正規化したもの)を何倍かしたものを足し合わせると、

過渡解析のある時間の結果が得られると書きましたが、

言い換えると、

任意の時間のどんな変形状態でも

この3つの固有ベクトルをそれぞれ何倍かした和によって

表すことができる

数式で書くと

ds474


このa,b,c(定数)が求められるということです。


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/06(水) 23:18:04|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

固有ベクトルにかかる係数を求める

では早速求めてみましょう。

例えば、何秒後かに

ds475

図で書くと

ds476

のような変形をしたとしましょう。



この場合は、ベクトル式で表すと

ds477


となります。

これは

0=-1.84a+3.68b-1.84c
1=-3.19a+3.19c
1=-3.68a-3.68b-3.68c

の連立方程式を解けば求まります。

a=-0.247
b=-0.090
c=0.066

となり、表すことができましたね…

って、

式3つに未知数3つなんだから求まるのは当たり前です…

これは作成したモデルのすべての固有モード形状を使えば、

そのモデルの節点が移動したどんな変形でも表すことができる(=正しい解になる)

ことをあらわしています。


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/07(木) 23:24:43|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

すべての固有モード形状を使わなくても計算できるとはどういうことか

前回、どんな変形もモード形状を何倍かしたものの和であらわすことを説明しました。

この例では全部で3つしか固有値がないので、すべてのモード形状の結果(モードベクトル)を使うことができましたが、

でも現実的にはこのような簡単なモデルを使用することはありません。。。

固有値は節点の数以上あるので、すべての固有モード形状を使うことは不可能です。

でも、モーダル過渡応答解析を行ったときには、すべてのモード形状を使わずに計算していましたよね。。。

1次モードから多くても20次モードまであれば十分などと書きました。

どういうことでしょうか。。。

答えを言うと、「計算が合うモデル」のみ計算していたのです。

たとえば

ベクトルであらわすと

ds480

というような形状が現れたとしましょう。

図で書くと

ds481

な感じでしょうか。。。

ds482


a,b,cを求めると、

a=1.5
b=1.2
c=0.0

係数cが0になりました。

係数cはモード3の形状にかかる係数でしたので、

つまりこの形状だとモード3の形状は計算に使っていない



省略できると言うことです!

仮想ブログ読者:「ちょっと待ってください。今回たまたまcが0になる形状だったのですが、一般的にはそうではないでしょう?」

わかりました。一般的にどうかは次回議論してみましょう。


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/12(火) 00:03:47|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

モードが打値切れる理由は…

前回、c(モード3の係数)=0の場合は、モード3は省略できると書きました。

この場合はモード1とモード2のみのモード形状があれば計算ができるということになります。

しかし、こんな都合のよい場合だけではないとおもいますよね。

たとえば、結果が

ds483

の場合、

ds484

a,b,cの値は

ds485

から、
a=-0.4944
b=0.0905
c=0.13248

でcの結果を無視できません。。。


ここで、少し以前確認した

モーダル法過渡応答を行うときは

「対象とする振動数の3倍程度の振動数が含まれるように固有モードの数を指定する」

ということを思い出してみましょう。

今回の固有振動数は

モード1:34.1Hz
モード2:93.1Hz
モード3 :127.1Hz

です。

このことから

もし、加振振動数が30~40Hzまでであれば、モード1とモード2で十分

つまりcは求めなくてもよい。

加振振動数が50Hzを越えてくればモード3も必要

つまりcを求めないと精度が悪くなる

ということになります。

うーん、こうなったら実際30Hzや50Hzの荷重をこのモデルにかけて

a,b,cの値を調べてみることにしましょう(^^)


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/12(火) 23:49:33|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

モード形状(固有ベクトル)にかかる係数の計算方法

さて今までは、変形をたとえば

ds483

のように仮定して、a,b,cの値を

ds485a


のように式を立てて逆算していました。

でもこれはあくまでも逆算であって、

実際シミュレーションする際は実際の変形がわかりません。、

このような方法でa,b,cはもちろん求めることはできません。


ではどうやって計算しているかという話は、もう少し後回しにしますが、

動解析の計算をしているので、動的な力の釣り合いを計算していることは予測できると思います。

具体的には、慣性力やら外力やらが変形する力とつりあっている

[M]{a}+[K]{u}={f}

http://freecaetester.blog62.fc2.com/blog-entry-589.html

という式が基本になります。


また減衰を考慮する場合は減衰力もこの式にいれなければなりません。。。

減衰力が速度に比例する、としてあげると

上の式は

[M]{a}+[B]{v}+[K]{u}={f}

という式になり、実際FEMソフトはこの方程式を解いています。

(ここのaは加速度ベクトルをあらわすaで、係数のaとは違います)

[B]は減衰マトリックス、{v}は速度ベクトルなのです。

この式の意味など詳しく追っていると、また長くなるので、

それはまた別の機会に書くことにして、とりあえず話を進めます(^^)。


[M]{a}+[B]{v}+[K]{u}={f}

何らかの処理

モードにかかる係数が求まる

モードにその係数をかけて変位が求まる


この何らかの処理についてはまた後でお話しすることにして、

荷重が30Hz,50Hzのときの

Nastranで作成したデータのモードにかかる係数を見てみたいと思います。


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/13(水) 23:05:00|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

Nastranでモード形状にかかる係数を求める場合の入力データ

それでは、Nastranでモーダル過渡応答解析のデータを作成します。

まず節点番号4に与える荷重の振動数が30Hzの場合、

Femapで作成して、Nastran入力ファイルにエクスポートすると、以下のようになります。


INIT MASTER(S)
NASTRAN SYSTEM(442)=-1, SYSTEM(319)=1
ID F0Hz,Femap
SOL SEMTRAN
TIME 10000
CEND
SDISP=ALL
TITLE = NX Nastran 時刻歴 解析セット
ECHO = NONE
DISPLACEMENT(SORT1,PRINT) = ALL
SPC = 1
DLOAD = 1
METHOD = 1
SDAMPING = 5
TSTEP = 1
MEFFMASS(PRINT,GRID=0,PARTFAC,MEFFM,MEFFW,FRACSUM,SUMMARY) = YES
BEGIN BULK
$ ***************************************************************************
$ Written by : Femap with NX Nastran
$ Version : 10.2.0
$ Translator : NX Nastran
$ From Model : E:\work\dynamic_seminar\modal_dynamic\NASTRAN\30Hz.modfem
$ Date : Thu Jun 14 23:03:08 2012
$ Output To : E:\work\dynamic_seminar\modal_dynamic\NASTRAN
$ ***************************************************************************
$
PARAM,POST,-1
PARAM,OGEOM,NO
PARAM,AUTOSPC,YES
PARAM,GRDPNT,0
CORD2C 1 0 0. 0. 0. 0. 0. 1.+FEMAPC1
+FEMAPC1 1. 0. 1.
CORD2S 2 0 0. 0. 0. 0. 0. 1.+FEMAPC2
+FEMAPC2 1. 0. 1.
EIGRL 1 10 0 MASS
$ Femap with NX Nastran Load Set 1 : Untitled
$ Femap with NX Nastran Function 5 : 臨界減衰比(ζ) 関数
TABDMP1 5 CRIT +
+ 0. .05 5000. .05 10000. .05ENDT
$ Femap with NX Nastran Function 3 : vs. 時間 関数
TABLED2 3 0. +
+ 0. 0. .001.0941083 .002.1873813 .003.2789911+
+ .004.3681246 .005.4539905 .006.5358268 .007.6129071+
+ .008.6845471 .009.7501111 .01 .809017 .011 .860742+
+ .012.9048271 .013.9408808 .014.9685832 .015.9876883+
+ .016.9980267 .017.9995066 .018.9921147 .019.9759168+
+ .02.9510565 .021.9177546 .022.8763067 .023.8270806+
+ .024.7705132 .025.7071068 .026 .637424 .027.5620834+
+ .028.4817537 .029.3971479 .03 .309017 .031.2181432+
+ .032.1253332 .033.0314108 .034-.062791 .035-.156434+
+ .036 -.24869 .037-.338738 .038-.425779 .039-.509041+
+ .04-.587785 .041-.661312 .042-.728969 .043-.790155+
+ .044-.844328 .045-.891007 .046-.929776 .047-.960294+
+ .048-.982287 .049-.995562 .05 -1. .051-.995562+
+ .052-.982287 .053-.960294 .054-.929776 .055-.891007+
+ .056-.844328 .057-.790155 .058-.728969 .059-.661312+
+ .06-.587785 .061-.509041 .062-.425779 .063-.338738+
+ .064 -.24869 .065-.156434 .066-.062791 .067.0314108+
+ .068.1253332 .069.2181432 .07 .309017 .071.3971479+
+ .072.4817537 .073.5620834 .074 .637424 .075.7071068+
+ .076.7705132 .077.8270806 .078.8763067 .079.9177546+
+ .08.9510565 .081.9759168 .082.9921147 .083.9995066+
+ .084.9980267 .085.9876883 .086.9685832 .087.9408808+
+ .088.9048271 .089 .860742 .09 .809017 .091.7501111+
+ .092.6845471 .093.6129071 .094.5358268 .095.4539905+
+ .096.3681246 .097.2789911 .098.1873813 .099.0941083+
+ .1-8.51-15ENDT
TLOAD1 101 102 LOAD 3
FORCE 102 4 0 1. 0. 1. 0.
DLOAD 1 1. 1. 101
TSTEP 1 100 .001 1
$ Femap with NX Nastran Constraint Set 1 : NASTRAN SPC 1
SPC1 1 123456 1
SPC1 1 13456 2
SPC1 1 13456 3
SPC1 1 13456 4
$ Femap with NX Nastran Property 100 : バー プロパティ
PBAR 100 1 2.5-55.208-115.208-11 0. 0. +PR 2S
+PR 2S 0. 0. 0. 0. 0. 0. 0. 0.+PA 2S
+PA 2S 0.
$ Femap with NX Nastran Material 1 : 等方性 マテリアル
MAT1 1 2.1+118.077+10 .3 7860. 0. 0.
GRID 1 0 0. 0. 0. 0
GRID 2 0 .25 0. 0. 0
GRID 3 0 .5 0. 0. 0
GRID 4 0 .75 0. 0. 0
CBAR 1 100 1 2 0. 1. 0.
CBAR 2 100 2 3 0. 1. 0.
CBAR 3 100 3 4 0. 1. 0.
ENDDATA


データの詳細はNastranのマニュアル見てください(^^;)

モード形状(固有ベクトル)にかかる係数の出力は

SDISP=ALL

で設定しています。

減衰は臨界減衰値の0.05を入力しています。

FemapはTload2を使えないようなので、表形式の荷重の入力(TABLED2)になっています。


もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/15(金) 00:06:39|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

Nastran f06ファイルの中のモード形状にかかる係数の結果が書いてある場所

では結果を見てみましょう

Nastranユーザーにはおなじみのf06ファイルです。

f06ファイルを開いてみていくと以下のようなところがあるはずです。

ds486

まずTime=0.0の

DISPLACEMENT VECTOR (SOLUTION SET)

と書いてあってその後に、その後に節点番号1の変位らしきものが書いてあります。。。

続けて
Time=1.0E-3

のものが同様に書いてあります。


実はこれが、各時間におけるモード形状にかかる係数、つまり

ds487

であらわしたときの、係数a、b、cなのです。

えっ、結果の表記ではPoint IDが1のx方向変位、y方向変位、z方向変位ではないかって…

どうも表記とは関係なく、モードの低い順から係数を出しているだけのようです。
(フォーマットは借り物…)



もしよかったらクリックしてください。順位が上がります。。。

にほんブログ村 科学ブログ 技術・工学へ
にほんブログ村

↓tsunodakoのつぶやき
tsunodakoをフォローしましょう
  1. 2012/06/17(日) 23:07:43|
  2. 解析講座
  3. | トラックバック:0
  4. | コメント:0

カレンダー

プルダウン 降順 昇順 年別

05月 | 2012年06月 | 07月
- - - - - 1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30


プロフィール

つのだこたろう

Author:つのだこたろう
2012年8月以前の記事は管理人別ブログ http://freecaetester.blog62.fc2.com/と重複します。

最新記事

最新コメント

最新トラックバック

月別アーカイブ

カテゴリ

未分類 (10)
振動解析 Part2 (44)
解析講座 (329)
ひとりごと (1)

検索フォーム

RSSリンクの表示

リンク

このブログをリンクに追加する

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。