振動工学にトライ2（固体連続体、平板振動）
7 ここまで、もう一つの境界条件（モーメントゼロ）が構成式の一般解を満足する事を確認していない。そこで今回は、その境界条件が一般解を満足しているかどうかを確認する。
例によって、Mapleで作成した解説を下に掲載する。
------------------------------------------------------------------------------
-----------------------------------------------------------------------------
8 最後に
不思議な事に、平板振動のz方向撓みを求める構成式から一般解を導出するのに、境界条件の片方が不要だった。こんな事に遭遇したのは初めてだ。
また平板振動については、FEMツールが言うことを聞かないので、Mapleだけの評価になってしまった事が残念である。
振動問題はこれくらいにして、次は音響関連にトライしてみる。
振動工学にトライ2（固体連続体、平板振動）
4 今回は後半の一般解完成までを挙げる。
モード形状、線形重ね合わせ、フーリエ級数の直交関係を利用した初期値係数の求解、そして最後に一般解完成まで辿る。ここまで記載されているテキストは見当たらない。
-------------------------------------------------------------------------------
---------------------------------------------------------------------------------
5 平板振動の初期条件が分かっていれば、関数で入力すれば結果が得られるはずだ。そのような例題をまだ探せていない。見つかったら紹介しよう。
気が付いた方がいるかもしれないが、もう一つの境界条件が構成式を満足する事を確認していない。そこで次回は、モーメントの境界条件が構成式を満足しているかどうかを確認する。
6 参考図書
今回に平板振動で参考した本を挙げてみる。
（1）音響の本ではあるが、平板振動に安田テキストが役に立った。
（2）JSMEのテキストにも平板振動が記載されています。
（3）フーリエ級数や直交関係については馬場テキストが解りやすかった。
振動工学にトライ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だ。
（5）最後に、FreeCADの参考書も紹介する。
![基礎からのFreeCAD[四訂版] (I/O BOOKS)](https://m.media-amazon.com/images/I/4113L46G-QL._SL500_.jpg "基礎からの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）竹内テキストは、FEMの数学的なバックグラウンドとして役立った3。
