4 過度応答の解説
Mapleによる解説の後半を掲載する。式番号が前のものもあるので、適宜前回の解説を参照されたい。ここで肝心な事は、数学的には、固有関数Yn(x)が直交性を持つ事を利用して積分定数C1nとC2nが求められる事だ。またMapleにおいて、離散計算データでも、配列にすることでプロット可能な事だ。
5 参考図書
梁の振動が記載されている本を色々挙げてみる。
（1）振動全般は藤田テキストが役に立った。
（2）麻生テキストには梁の様々な振動形態が紹介されている。個人的には最も解析的な本と感じた。
（3）建築畑からは、宮本テキストを挙げる。機械系と若干言葉が異なる所も趣がある。
（4）JSMEテキストにのみ、連続体としての平板の曲げ振動が載っていた。流石JSMEだ。
（5）最後に、FreeCADの参考書も紹介する。
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。
今回から円柱流れを考えてみる。
これまでのパターンでは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.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のコードとその説明は次回とする。
4 FlexPDEによるFEM計算結果
（1）計算モデル
円柱座標でかつ軸対象流れ。
内径d=50mm、長さL=0.8mの直円管に、密度ρ=900kg/m3、粘度μ=0.02Pa s、の油を体積流量Q=12 litter/min流し込む。
（2）結果
・モデルのグリッド（自動グリッドでノード24、セル26）...図1
長手方向がz方向で、管断面（緑色）は半割状態で計算した、よって管中心はR=0である。
・管路内の流速(m/s)とせん断応力(Pa) (vzとtauエレベーション図)...図2a,b,c
上から、管入口部、管中間部、管出口部を示す。オレンジ色がz方向流速vz、ブルーがせん断応力tauである。
図2a
境界条件として入口に平均流速umean=Q/ (πr^2)を採用した（Boundaries分でvalue(vz)=umean）ので、図2aのvzは管中心で約0.1m/sを示している。
vz、tauいずれも入口においては物理的な意味を正確には表せていない。
図2b
管中間部においては、vzは管中心で最大になり、約0.145m/s。
せん断応力は最大が管壁にあり、約0.24Paである事が判る。Maple解説では0.33Paなので、73%程度に見積もられている（理由はRe）。
図2c
出口部においては、管中心（R=0）で流速vzが最大で、約0.15m/s。Maple解説では約0.2m/sなので、これも25％小さく見積もられている。
せん断応力は最大が管壁にあり、約0.24Paである事が判る。
・最大流速(vmのコンター図)...図3
vmはベクトルv（vz、vr）の大きさ（マグニチュード、つまりsqrt(vz^2+vr^2)）を表す。管軸R=0ではvrの影響はずいぶん小さく、殆どvzが支配的と言って良い、0.15m/s。
・レイノルズ数（Reのコンター図）...図4
Maple解説では教科書通り代表長さにd、平均流速umeanを選び、Reは約229である。図3からFEMでは167.8なので、この値も小さく73％程度に見積もられている。値は全体的に同じと見て取れる。
・圧力分布（pのコンター図）...図5
出口は境界条件でvalue(p)=0を指定しているので、0Pa。入口は約19Paに計算されている。Maple解説では20.8Paなので、9割方合っていると言える。
全体図では判り辛いが、入口を入った直後の管壁にせん断応力の最大値0.37Paが見られる。管中間部は図2bと同じ。
・質量保存則（div_vのコンター図）...図7
円柱座標における連続の式を計算した。ムラは見えるがカラーバーから全体的にゼロと見て取れる。うまく発散して質量保存則が殆ど成立している事が判る。
・渦度（θ回りの回転のコンター図）
円柱座標での流速vの回転は3つ存在するが、軸対象を仮定すると1つしか残らない。それを計算したものである。傾向はせん断応力tauに近い。決して非回転流れとは言えないね、粘度を加味している結果が物語っている。
（3）コードの説明
ここまで来て、コードを見たいなと思ってもらうと嬉しいけど...。
------------------------------------------------------------------------------
---------------------------------------------------------------------------------
プログラムコードの補足を加えたい。
TITLEは概要や計算の背景などを記載しても構わない。参考図書も記録しておくと後で役立つ。
SELECTでは計算誤差の値をerrlimで指定できる。今回はアプリケーションの制限だろうが、2.5e-2以下の計算では乱流を呈してしまったので、ここでリミットをかけた。
COORDINATESは座標系を指定するが、円柱座標で計算するのでycylinderを指定した。VARIABLESでは変数を指定する。
DEFINITIONSでは構成式中の定数や定義などを指定する。本問題ではモデルの形状、粘度、密度を指定した。単位系は原則SI。
またレイノルズ数基本式、umean（入口流速）の基本式、流速のベクトル定義、速度ベクトルの発散の式（質量保存）、速度ベクトルの回転の式（渦度）、せん断応力tauの基本式を入れている。手元にある流体力学の教科書を参考にしてほしい
EQUATIONSでは、その名の通り構成方程式を指定する、今回はNSEQと圧力pの式を利用した。円柱座標で計算させるためにデカルト座標の構成式を使わず、修正したものになっている。
BOUDARIESは境界条件を指定するエリアで、肝の所だ。valueはディリクレ形境界条件で、入口でvz=umean、出口でp=0を指定。加えて管壁ではvz=0（すべりなし）とした。管軸中心においては、vr、vz、pをnatural（ノイマン）形境界条件とした。
PLOTSは様々なplotコマンドを入力するエリアだ。contourコマンドを多用しているが、具体的な数値が欲しい時は座標を指定してelevationコマンドを使った。
contourコマンドの後ろにpaintedオプションを使うと、記号コンターから色コンターに変えられる。全体を俯瞰したい時に使える。またreportオプションでは指定した計算値を図の下に出力してくれる。例えば図3ではglovalmax(vm)とReが出力されている。これでReの元になった代表流速値が判るわけだ。
最後はENDでプログラムを閉じる。
5 最後に
面倒な数値計算テクニックが無くても、この程度のコードでNSEQが解けてしまう、非圧縮性で層流の範囲ではあるが、現実流体の粘性が加味できるので良しとしよう。
FlexPDEによるFEM結果はMaple解説での値より、せん断応力、最大流速およびレイノルズ数などの数値が小さく計算されている。これはレイノルズ数の計算式が異なる事に原因が有りそうだ。
このように複数のアプリケーションで同一問題を解かせると違いが生じる事が往々にしてある。大事なことは、何故そうなのか、原因がどこにあるのかを突き詰めてみることだ。
6 参考図書
（1）NSEQなどの座標変換に関しては、岡本テキストが役に立った。
（2）巽テキストは王道です。
（3）FlexPDEのコードについてはGunnar Backstromのテキストが役に立った。
3 FEM計算結果
FEMシミュレーション結果はFlexPDEで示すが、計算に用いる構成式を少し補足しておく。
（1）ナビエ・ストークス方程式（以下NSEQ）
（2）圧力の方程式
（1）NSEQ
この方程式の導出は手もとの教科書を参照してほしい。質量保存則（オイラーの連続の式）に、運動量保存則（オイラーの運動方程式）とニュートンの粘性法則を加味して導かれている。ここではFEMに載せる形までを与える。
通常の教科書では、非圧縮性流体の場合は以下のように記述されている。ボールド体はベクトルである。
ここで、u(m/s)：流速、t(s)：時間，ρ(kg/m3)：流体密度、p(Pa)：圧力，μ(Pa s)：流体粘度、F(N/kg)：単位質量当たりの外力。
式(1)の左辺は慣性項、右辺第1項は圧力項、第2項は粘性項（拡散項という表現もする）、第3項は外力項という。
左辺のD/Dtはラグランジュ微分（実質微分ともいう）で、中味を書き下すと以下の様になる。
式(2)の右辺第1項を非定常項、右辺第2項を移流項（対流項という表現もする）。
ここで、・は内積を示す。少しは読み易くなっただろうか。
デカルト座標表示にする前に、問題を少しシンプルにする。
まず、定常を扱うことにする、よって非定常項は削除。
両辺にρを掛けて、
また、外力を無視する（つまり重力を考えない仮定を入れる）と、
この式(6)をデカルト座標表示に変える。冗長だが構造が理解し易くなる。
ここで、u(m/s)：x方向流速、v(m/s)：y方向流速、w(m/s)：z方向流速。
・円柱座標
次に、円管流れに良く用いられる円柱座標(r,θ,z)に式(7)を変換しておこう。これも数学的に導くのは他の教科書に譲って、結果のみ掲載したい。
ここで、ur、uθ、uz：それぞれr、θ、z方向速度(m/s)。
式(7)と式(8)は、座標を変えた3Dだが、計算機的には確実に式(8)は負荷が大きい。
・軸対称流れ
次に円管流れの中でも、軸対称での解析で十分な場合（回転流は無視）はθ項を簡単にできる。つまりuθ＝0,θによる微分もゼロとすると、
これで、だいぶシンプルになった。
（2）圧力の方程式
式(9)をみれば、未知数はur、uz、pの3つなので、圧力pの方程式も用意しなければ NSEQは解けない。まず座標に依存しない式(6)の両辺に発散を施すのだが、その前に若干修正しよう。
両辺に発散を施すと、
式(10)の左辺第3項に、以下の公式(11)を用いると、
連続の式 ∇・u = 0から、式(12)の左辺第3項はネグれるが、式(9)と式(12)のPDEが、「全ての場合において連続の式を満たす」と仮定するのは正しくない。実際は特別の場合にだけ成立する事が判っているからだ。問題を解いて行く過程で、最終的に発散は消えてしまうので、式(12)にある係数fvを掛けた∇・uを加算しても決して誤りではない。
ここで、
ここでC：定数、経験的に10000とする。L(m)：代表長さ。
式(13)に式(14)を代入して、最終的に円柱座標に変換しよう。座標変換手順は例に依って割愛するので、数学の参考書に譲る。
FEM計算の構成式は式(9)と式(15)である。だいぶ話が長くなったので、FlexPDEの結果は次回に持ち越す。
![基礎からのFreeCAD[四訂版] (I/O BOOKS)](https://m.media-amazon.com/images/I/4113L46G-QL._SL500_.jpg "基礎からのFreeCAD[四訂版] (I/O BOOKS)")
