振動工学にトライ2（固体連続体、平板振動）

4　今回は後半の一般解完成までを挙げる。

モード形状、線形重ね合わせ、フーリエ級数の直交関係を利用した初期値係数の求解、そして最後に一般解完成まで辿る。ここまで記載されているテキストは見当たらない。

-------------------------------------------------------------------------------

---------------------------------------------------------------------------------

5　平板振動の初期条件が分かっていれば、関数で入力すれば結果が得られるはずだ。そのような例題をまだ探せていない。見つかったら紹介しよう。

気が付いた方がいるかもしれないが、もう一つの境界条件が構成式を満足する事を確認していない。そこで次回は、モーメントの境界条件が構成式を満足しているかどうかを確認する。

6　参考図書

今回に平板振動で参考した本を挙げてみる。

（1）音響の本ではあるが、平板振動に安田テキストが役に立った。

機械音響学

（2）JSMEのテキストにも平板振動が記載されています。

振動学 (JSMEテキストシリーズ 6-1)

（3）フーリエ級数や直交関係については馬場テキストが解りやすかった。

偏微分方程式キャンパス・ゼミ 改訂6

振動工学にトライ2（固体連続体、平板振動）

1　はじめに

連続体の振動の2回目は平板を取り上げ、複数のアプリケーションで評価してみたい。

これまで通り、出来る限りシンプルに進めたいので3Dを目指すが、計算量の関係で2Dを採用せざるを得ない場面が出てくるかもしれない。

2　解説

いつもながら、このブログの趣旨から外れ無いよう、構成式の導出は極力避ける。詳しくは手元にある材料力学や振動工学の教科書を参照して頂きたい。

今回は薄板（キルヒホッフ・ラブ理論）の曲げ振動（横振動）の計算例を示す。計算例と強く連動しているので、構成式のみを天下り的に利用している。

偏微分方程式（PDE）の構成式を、変数分離法により複数の常微分方程式（ODE）に変えて、拘束（境界条件）から固有値を求め、固有角振動数ω(rad/s)と固有周波数f(Hz)まで計算する。例によってMapleによる解説と具体的な計算例を下に掲載する。

当然ながら方持ち梁より難しくなって、数式が非常に多くなった。その分丁寧に中味を追った記述を志したので勘弁願いたい。Mapleコマンドも随所で補足している。

---------------------------------------------------------------------------------

----------------------------------------------------------------------------------

3　構造体で固有振動数が重要視されるのは低次のもので、実際にはn=3程度までであろう。このタイプは構成式導出までに以下の仮定が設けられている。

・板厚は薄く一定で、撓みは小さいので、平面応力問題とする。

・板の中心面に垂直な平面は、中央面に対して垂直を保って撓む（キルヒホッフ仮定）。

・板厚方向のせん断による変形は無視。

・材質は、弾性材料で、任意の断面において均質。

厚みのある板を扱う場合は、せん断変形も考慮したミンドリン・レイスナー理論を考える必要がある。

次回は、固有関数（モード形状）を解説する。

振動工学にトライ（固体連続体）

4　過度応答の解説

Mapleによる解説の後半を掲載する。式番号が前のものもあるので、適宜前回の解説を参照されたい。ここで肝心な事は、数学的には、固有関数Yn(x)が直交性を持つ事を利用して積分定数C1nとC2nが求められる事だ。またMapleにおいて、離散計算データでも、配列にすることでプロット可能な事だ。

5　参考図書

梁の振動が記載されている本を色々挙げてみる。

（1）振動全般は藤田テキストが役に立った。

振動工学 新装版:振動の基礎から実用解析入門まで

（2）麻生テキストには梁の様々な振動形態が紹介されている。個人的には最も解析的な本と感じた。

機械力学 (新機械工学シリーズ)

（3）建築畑からは、宮本テキストを挙げる。機械系と若干言葉が異なる所も趣がある。

構造工学（第４版）

（4）JSMEテキストにのみ、連続体としての平板の曲げ振動が載っていた。流石JSMEだ。

振動学 (JSMEテキストシリーズ 6-1)

（5）最後に、FreeCADの参考書も紹介する。

基礎からのFreeCAD[四訂版] (I/O BOOKS)

振動工学にトライ（固体連続体）

3　FreeCADによるCAE結果

片持ち梁の曲げ振動の計算例として、FreeCADに搭載されているFEMツールで実行したものを紹介する。

前回、「次回は初期条件を利用する過度応答を解説する」としていたが、数式で目が疲れるのでMapleの解説は後回しにした。

FreeCADはオープンソースの3DCADで誰でも自由にダウンロードして使用できる。日本語参考書も多数あり、画面もほとんど日本語にローカライズされているので非常に使い勝手が良い。解析に必要なメッシュ機能や応力解析や振動解析、結果を描画する機能も有するので、オールインワンと言って良い。将来、市販のハイエンド専門ツール（SolidWorks、ANSYSなど）を使う際でも参考になるはずだ。

最初に、全ての作業が終了した後の画面全体を見て欲しい。左側（コンボビュー）のモデルタブをクリックすると、処理が下のようにツリー状に見える。

・Body（CAD図）

・Analysis（解析ファイル）

・MaterialSolid（材料設定）

・ConstraintFixed（課した境界条件）

・FEMmeshNetgen（メッシュ作成）

・SolverCcxTools（計算ソルバーCalculiX起動）

・CCX_EigenMode_1_Results（固有モード1の計算結果）

・Pipeline_CCX_EigenMode_1_Results（固有モード1の計算結果描画）

・ccx_dat_file（CalculiX計算結果ファイル）

右側（3Dビュー）には作成したCAD図が見える。片持ち梁なので、拘束かけた箇所が赤いスタンドで固定されている事が判る。

図1　作業終了後の画面全体

梁の形状寸法は前回Maple解説に示した通り。正方形断面（1cm x 1cm）で長さ1mである。

次にFreeCADに内蔵されている標準のメッシャーNetgenでメッシュ切りした条件を示す。

表1　FEMメッシュデータ

ノード数が40000を超えているので、十分なメッシュである事が解る。梁の断面が長さに比べて小さいので、メッシュ後の絵を見ながら、気を遣って2次精度でメッシュ切りを行った。細長い場合にはこのような慎重さが重要である。

次にソルバーCalculixによるFEM計算に必要な材料データを示す。

表2　FEM材料データ

以下にモード1、3、5ならびに7のx方向変位（上下）を示す。モード次数が偶数の箇所はy方向（奥行き）への曲げ振動を表しているようだ。値が同じなので掲載しなかった。

梁が曲がっている図は大袈裟にスケールアップしているので、傾向だけを掴んでもらえれば良い。モード次数が高くなるにつれて節の数が増えていく。

図2(a)　モード1

図2(b)　モード3

図2(c)　モード5

図2(d)　モード7

最後にccx_dat_fileの一部を抜粋して載せる。

表3　ccx_dat_file

最上段のパートは、左からモード次数n、固有値、固有角振動数ωn、固有周波数fnである。一番右端は虚数部なのでゼロになっている。前回のMaple解説の式(70)から(77)と合致している事が判る。

上から2段目は、PARTICIPATION FACTOR（刺激係数）で、値が大きいほど励振方向による振れやすさを表している。モード次数1と2を見ると、x成分とy成分の値が反対になっている。また回転も見られる。

上から3段目は、EFFECTIVE MODAL MASS（有効質量）で、値の大きな振動モードほど、元の振動に寄与している。

次回はMaple解説の過度応答に進む。

振動工学にトライ（固体連続体）

1　はじめに

今回から連続体の振動の問題を取り上げ、複数のアプリケーションで評価してみたい。

これまで通り、出来る限りシンプルに進めたいので3Dを目指したいが、計算量の関係で2D（x,t）を採用せざるを得ない箇所が出てくる可能性が高い。

2　解説

これまで通り、このブログの趣旨から外れるので、構成式の導出は極力避けている。詳しくは手元にある材料力学や振動工学の教科書を参照して頂きたい。

1回目は片持ち梁の曲げ振動の計算例を示す。今回は計算例と強く連動しているので、構成式のみを天下り的に掲載した。

偏微分方程式（PDE）の構成式を、変数分離法により複数の常微分方程式（ODE）に変えて、梁の境界条件から固有値を求め、固有角振動数と固有周波数までを求める。例によってMapleによる解説と具体的な計算例を下に掲載する。緑色の文章は参考まで載せた箇所なので、識別出来るように緑色にした。

---------------------------------------------------------------------------------

-----------------------------------------------------------------------------------

構造体で固有振動数が重要視されるのは低次のもので、実際にはn=3程度までであろう。梁の構成式も単純梁（オイラー・ベルヌーイ理論）を用いた。このタイプは構成式導出までに以下の仮定が設けられている。

・梁の中心軸は引っ張りまたは圧縮を受けない。

・せん断による変形は無視（変形前の梁の中心軸に垂直な断面は平面を保ち、変形後も中心軸に垂直）

・回転慣性の影響を無視。

・材質は、弾性材料で、任意の断面において均質。

これ以外でも高次になるほど、節になる箇所が増えるので、節間の距離が減少していって、梁が波打つような形となり厳密性が損なわれていく。そしてモード次数が増加するにつれて回転の運動エネルギーを無視できなくなる、などの注意点がある。

梁の断面積が長さに比べて大きい場合は、より回転運動やせん断変形も無視できなくなるので、これらの影響を考慮したい場合はチモシェンコ梁を考える必要にがある。

次回は初期条件を利用する過度応答を解説する。

流体力学にトライ（円柱流れ）

層流のなかでも、円管を流れるハーゲンポアズイユ流れの次に、流れに逆らうようにおいた物体に働く影響、今回は円柱に働くケースだ。

3　FlexPDEによるFEM解（数値計算）を11月17日に示しているので、計算に用いたコードについて補足する。まず説明する前にアップしよう。

（1）コード

--------------------------------------------------------------------------------

-----------------------------------------------------------------------------

(2)　コードの説明

このブログの趣旨（プログラミングは最小限にしてアプリケーションを有効利用する）から外れるかもしれないが、数値計算の部分をコーディングしなくてもいいと考えればどうだろう。

TITLE文：これまでも記載した通りフリーに書いて構わない。

SELECT文：今回はstagesコマンドを使って、流入口からのx方向流入速度vx0（一様流れ、ベンチマークでは放物線形状）を少しずつ増やして行った。乱流に遷移しない様に慎重に、試行錯誤を重ねた結果だ。

COORDINATES文：今回x-y二次元デカルト座標を用いるので、省略可能である。

VARIABLES文：vx、vyおよびpが未知数である。

DEFINITIONS 文：構成式中の定数や定義を指定する。本問題ではモデル形状（ダクト幅や長さ及び円柱半径r＝0.05mなど）と粘度や密度を指定した。またvx0の値を1/100000m/sから0.2m/sまで段階的に分け、stagesコマンドと連携させた。MReは修正レイノルズ数と呼ばれているもので、代表速度は形状内の最大流速値を拾うようになっていて、一定ではない。v＝vector(vx, vy)はそのままベクトルにしたもの。vm=magnitude(v)はベクトルvの大きさ。nx=normal(vector(1, 0))は単位ベクトルの法線方向成分を意味している。

NSEQをそのまま表示するとコードが長くなるので、繰り返し計算する箇所を記号化したものである（16aと16b参照）。

8月29日にアップした原稿の(7)式のz座標を削除したものを記載する。

またdel2(vx)はFlexPDE特有のラプラス演算子で以下を表す。

natp式については、最後に説明する。

EQUATIONS文：構成方程式を指定するエリアなので、未知数vxとvyおよびpの3本である。圧力の式は8月29日にアップした式(13)そのままで、dens_termという記号を用いている。

BOUNDARIES文：境界条件を指定するエリアで、valueはディリクレ型、naturalはノイマン型で境界値を与える。境界条件を指定し易くするため形状全体をdomain、ダクト外形をouter、円柱をcylinderで区別した。

pについては、流入口と壁側さらに円柱表面にnatural(p)=natp式を指定。出口側にのみvalueで大気圧（0）を指定している。

vxについては、ダクト流入口にvalue(vx)=vx0（最終的にベンチマーク指定0.2m/sに達する様にstagesコマンドで制御）。そしてダクト壁や円柱表面は滑りなしとしてvalue(vx)=0ゼロを指定。ダクト出口は何かしらの影響があるかもしれないのでnatural(vx)=0とした。

vyについては、ダクト壁や円柱表面は滑りなしとしてvalue(vy)=0を指定。ダクト流入口や出口については、何かしらの影響があるかもしれないのでnatural(vy)=0とした（逆にゼロ指定はしない方がいい）。

PLOTS文：今まで見てきたとおり様々なプロットができるエリアだ。vectorコマンドの後ろにオプションzoomを記載すると、拡大したい箇所を見る事が出来る便利な機能だ。この機能で双子渦がはっきり再現されている。この機能はcontour(p)でも利用できる。

END文：最後はENDでプログラムを閉じる。

(3)　natp式の説明

圧力pが既知であればvalue文で指定出来る。しかし流れが存在する場合は他の境界条件を何とか持ってくるしかない。スカラー場pの勾配∇pは外向きの法線ベクトル場と見なす事が出来る。これを使ってみよう（イメージとして図10参照）。

図10　勾配の幾何学的表現

natp式の原型は8月29日にアップした式(6mod)である。圧力項を左辺に持ってくると、

これを利用して、以下の式を考える。

ここで、nは境界の各点において外向きに立てた単位法線ベクトル。

式(19)に作用させると、

さらに具体的に展開していく。

ここで、nx、nyはそれぞれnと座標軸（x-y）とがなす方向余弦を表す。

ここで、先に挙げた式(16aと16b)を用意した理由が明らかになったはずだ。

4　最後に

本来ならFlexPDEと他のアプリケーションとの計算結果比較を目的としていたが、以下のような理由で断念した。Featflowベンチマーク結果との比較に終わった事を残念に思う。

・流体解析ツールはLinux系オープンソースが多い。Linux系のOSでは、インストールそのものに手間が掛かる。

・windowsやmacで、フリーライセンスで使えるものがほとんどない。

・このブログの性格上、プログラミング言語系（Fortrun、Python、C、Pascal）は避けたい。

5　参考図書

（1）境界条件の与え方は河村テキストが役に立った。内容は簡単ではないが、流体解析を目指すなら一読を勧める。

流体解析 2 (応用数値計算ライブラリ)

（2）竹内テキストは、FEMの数学的なバックグラウンドとして役立った3。

計算力学(第2版):有限要素法の基礎

流体力学にトライ（円柱流れ）

今回から円柱流れを考えてみる。

これまでのパターンではMapleで作成した解説を先に公表していたが、今回は他アプリケーションの結果をベンチマークとして、FlexPDEのFEM結果を比較してみたい。

1　ベンチマークとして利用したアプリケーション

FeatflowのHPを開いてBenchmarksをクリックすると、Single Phase Navier-Stokes Equationsのタイトルで円柱流れの結果が4件見つけられるはず。一番初歩的な「RE20, Laminar」のケースを今回は採用した。計算条件は以下の通り

（1）計算モデル形状

X-Y方向の2Dモデル。幅0.42m（Y）、長さL 2.2m （X）のダクトに流体を流す（深さ方向Zは1mとする）。ダクト流入口から0.2mの位置に、流れに抵抗するように半径r 50mmの円柱をZ方向に立てる。

（2）流体条件

密度ρ 1.0 kg/m3の気体。温度は規定なし。動粘度ν 0.001 m2/s（ρ=1.0なので、粘度μは0.001 Pa・sとなる）。ダクト入口から一様流速u 0.2 m/sで流し込む。

（3）境界条件

vx=0：ダクト壁、円柱表面

vy=0：ダクト壁、円柱表面

p=0：ダクト出口

2　FlexPDEの結果（FEM計算）

（1）モデルのグリッド（自動グリッドでノード105、セル165）...図1

流体は左側から流入する。

（2）修正レイノルズ数（MRe）...図2

カラーバーを見た通り全体的に約32.17なので、まだ層流状態である。ちなみにFeatflowベンチマークではRe20。

（3）流速コンター（vm、つまりsqrt(vx^2+vy^2)）...図3

円柱と壁の間が流速最大0.36m/sで、円柱後部は流速最低0m/sである事が判る。またダクト出口では0.28m/s最大で放物形に分布していることが図8で判明する。

先のベンチマークの流速コンター図で読み取れた値を下に表記する。概略近い値を示している。

　　　　　　　　　　　　　　　　FlexPDE  　　Featflow

　　　　円柱ー壁間最大部（m/s）　　0.36　　　　0.4

　　　　円柱後部（m/s）　　　　　　　0　　　　　0

　　　　ダクト出口最大部（m/s）　　0.28　　　　0.28

（4）流速ベクトル（vx,vy）の拡大...図4

円柱後部で渦を巻いている、つまり双子渦の発生がはっきり判る。流体力学のテキストではレイノルズ数 6<Re<40の範囲で双子渦が見られると記載ある。

（5）圧力コンター（p）...図5

ダクト流入口の圧力は約0.07Pa、円柱先端部は0.13Paで最大圧力を示す。ダクト出口がゼロ（大気圧）になっているのは境界条件で指定しているからだ。

先のベンチマークの圧力コンター図で読み取れた値を下に表記する。概略近い値を示している。

　　　　　　　　　　　　　　FlexPDE  　Featflow

　　　　ダクト流入口（Pa）　　0.07　　　　0.074

　　　　円柱先端部（Pa）　　　0.13　　　　0.12

（6）圧力コンター（p）の拡大...図6

円柱左側に圧力の最大点0.13Paが、円柱最大径を過ぎた直後に圧力の最低点-0.01Paが見える。先のベンチマークの圧力コンター図も同様な絵である。

そのベンチマークには、円柱前後の差圧が0.1Pa実測されたことも記載されている。

ともかく圧力による抗力発生が容易に予想される絵だ。

（7）入口流速（vx）...図7

境界条件で一様流速0.2m/sを指示した通りになっている（ム？ギザギザじゃないの？と訝った方がいるかもしれないが、無料使用の範囲では精度がこの程度である）。

（8）出口流速（vx）...図8

0.28m/sを頂点に放物形の分布を示している事が判る。

（9）質量保存則（div_vのコンター）...図9

連続の式を計算した。カラーバーから、円柱付近を除き全体的にゼロと見て取れる。うまく発散して質量保存則が成立している証拠だ。

計算に使用したFlexPDEのコードとその説明は次回とする。