機械回路の記号解析 その2
1K Views
February 10, 25
スライド概要
はじめに
1. モビリティ法
1.1 モビリティとは
1.2 機械要素のモビリティ
1.3 機械回路と電気回路の対応関係
1.4 モビリティの直列、並列接続
1.5 直列質量について
２．モビリティの機械回路例
2.1 機械的蓄音機のモデル
2.2 電気回路LTspiceでの解析
2.3 機械回路のモビリティ表現
3. サイズモ系例
3.1 計測原理
3.2 変位計測
3.3 加速度計測
４．動吸振器例
4.1 動吸振器の運動方程式
4.2 電気回路LTspiceでの解析
4.3 動吸振器(AC解析：定点理論)
4.4 動吸振器（最適化例）
５．クランク軸のねじり振動
5.1 運動方程式（階差方程式）
5.2 階差方程式の解
5.3 クランク軸のねじり振動(例題)
5.4 クランク軸のねじり振動(例題)の解
5.5 電気回路LTspiceでの解析
BODE線図
過渡応答
6 連続体のモビリテイ
6.1 弾性梁のモビリテイ(方程式)
6.2 弾性梁のモビリテイ(静的変位)
6.3 弾性梁のモビリテイ(動的変位)
6.4 弾性梁のモビリテイ(共振)
6.5 弾性梁のモビリテイ(疑似集中定数系)
6.6 連続体を含む回路の計算（例）
つづく
補足：中田孝氏について
これまでに主に，ロボティクス・メカトロニクス研究，特にロボットハンドと触覚センシングの研究を行ってきました。現在は、機械系の学部生向けのメカトロニクス講義資料、そしてロボティクス研究者向けの触覚技術のサーベイ資料の作成などをしております。最近自作センサの解説を動画で始めました。https://researchmap.jp/read0072509 電気通信大学 名誉教授
関連スライド
各ページのテキスト
2025.2.10 機械回路の記号解析 その2 電気回路アナロジーによる解析手法 機械系のためのメカトロニクス サイズモ系 動吸振器 航空エンジンとプロペラの振動回路 下 条 誠 電気通信大学名誉教授 https://researchmap.jp/read0072509/ クランク軸系の振動回路 The University of Electro-Communications https://www.docswell.com/user/m_shimojo Department of Mechanical Engineering and Intelligent System
内 容 はじめに 1. モビリティ法 1.1 モビリティとは 1.2 機械要素のモビリティ 1.3 機械回路と電気回路の対応関係 1.4 モビリティの直列、並列接続 1.5 直列質量について ２．モビリティの機械回路例 2.1 機械的蓄音機のモデル 2.2 電気回路LTspiceでの解析 2.3 機械回路のモビリティ表現 3. サイズモ系例 3.1 計測原理 3.2 変位計測 3.3 加速度計測 ４．動吸振器例 4.1 動吸振器の運動方程式 4.2 電気回路LTspiceでの解析 4.3 動吸振器(AC解析：定点理論) 4.4 動吸振器（最適化例） ５．クランク軸のねじり振動 5.1 運動方程式（階差方程式） 5.2 階差方程式の解 5.3 クランク軸のねじり振動(例題) 5.4 クランク軸のねじり振動(例題)の解 5.5 電気回路LTspiceでの解析 ① BODE線図 ② 過渡応答 6 連続体のモビリテイ 6.1 弾性梁のモビリテイ(方程式) 6.2 弾性梁のモビリテイ(静的変位) 6.3 弾性梁のモビリテイ(動的変位) 6.4 弾性梁のモビリテイ(共振) 6.5 弾性梁のモビリテイ(疑似集中定数系) 6.6 連続体を含む回路の計算（例） つづく 補足：中田孝氏について 2
はじめに 3 機械の振動系と電気回路が同じように取り扱えると、知ったのは中田孝氏の書籍「工 学解析」でした。同書は、Lagrange運動方程式、Fourie解析、変分法、テンソル解 析など物理学者の数学的道具を機械工学を学ぶ学生に分かり易く解説したものです。 その中に機械回路の記号解析がありました。これがアナロジーというものに関心を 持った始めでした。 機械系と電気系のアナロジーについて学ぶに従い、「物理量の捉え方」の一つとして 流通量と位差量があること、これにより機械と電気系の物理変数の対応関係ができる こと、さらに流体系、熱回路、音響工学など広い分野の物理現象の解析にアナロジー の考え方が使えることに興味を覚えたものです。また双対性の概念は、各種物理現象 への新たな視点を得るきっかけとなりました。 さて、機械系メカトロニクス講義の続きとして、新たに機械回路をはじめます。これ は機械振動系を電気回路理論で理解し、回路シミュレータにより解析を行います。電 気電子技術に馴染みがない機械系学生に少しでもエレクトロニクスへの興味，理解を 深める機会になれば幸いです。今回は、電気回路でのインピーダンスに対応するモビ リティという概念を用いて機械回路の解析を行います。加速度/変位センサなどのサ イズモ系や、クランク軸のねじり振動の解析、そして連続形のモビリティについ解説 します。 参考文献 中田孝:工学解析(技術者のための数学手法),オーム社,1972. pythonのprogram： https://github.com/m4881shimojo/TTanaka_Eng/
１. モビリティ法 前回は、機械系の，”速度と力”を，電気系の”電圧と電流”に対応させた 力－電流アナロジーについて解説しました。 https://www.docswell.com/s/m_shimojo/57REQ7-2024-11-05-140854 今回は、電気回路でのインピーダンスに相当する概念をモビリティとしてい ます。また機械系では位置と力でシステムの状態を表現するのが普通です。 そこで，位置の時間微分を用いたモビリティ法について解説します 機械回路 電気回路 λ: mobility Z: impedance よく使うパラメータ 重要なパラメータ 変位と力 電圧と電流 このため微分演算子ｓを導入 4
1.1 モビリティ (mobility)とは 5 機械要素を通過する力と速度差との関係 機械系 v(t) 速度 f(t) f3 k f2 d f1 電気系 i3 電圧 e(t) L i2 R i i1 C m 𝑓 = 𝑓1 + 𝑓2 + 𝑓3 𝑖 = 𝑖1 + 𝑖2 + 𝑖3 𝑑𝑣 𝑓=𝑚 + 𝑑𝑣 + 𝑘 න 𝑣 𝑑𝑡 𝑑𝑡 𝑑𝑒 1 1 𝑖=𝐶 + 𝑒 + න 𝑒 𝑑𝑡 𝑑𝑡 𝑅 𝐿 モビリティ－： 𝑣 𝜆= 𝑓 インピーダンス： 𝑍= 𝐸 𝐼
1.1 定義： モビリティ (mobility)とは 𝑣1,2 = 𝑣1 − 𝑣2 = 𝜆𝑓 モビリティ(mobility) 機械回路 𝒇 𝒇 λ 機械要素を通過する 力と速度差との関係 𝒗𝟏,𝟐 𝑣1 速度差 𝑣2 通過する力 λ: f: v: mobility force velocity 1. 力ー電流アナロジーでは、機械系の”速度と力”を電気系の”電圧と電流”に対応させています。 モビリティ法でも同様な対応付けです 2. しかし、機械系では変位と力でシステムの状態を通常表現するため、位置の代わりにその時間 微分である速度を使う点なじみにくいかも知れません 3. そこで、微分演算子“ｓ”を導入して、変位による表現も可能とします 演算子： 𝑑 =𝑠 𝑑𝑡 𝑠 −1 = න 𝑑𝑡 6
1.2 機械要素のモビリティ ば ね 𝑣2 質 量 𝑣1 𝑎 = 1Τ𝑑 𝑣2 𝑣1 𝑓 𝑓 機 械 回 路 電 気 回 路 抵 抗 𝑐 = 1Τ𝑘 𝑣1 7 𝑣1,0 𝑉1 𝑣1,2 𝑣1,2 𝜆 = 𝑐𝑠 𝜆=𝑎 I 𝑉2 𝑉1,2 𝐿 𝑉1 m 𝑣0 𝜆 = 𝑚𝑠 −1 I 𝑉2 𝑉1,2 𝑣1 𝑚 R 𝑧 = 𝑗𝜔𝐿 𝑧=𝑅 λ: mobility z: impedance 𝑉1 𝑉1 𝑉1,0 I 𝐶 𝑧 = 𝑗𝜔𝐶 −1 演算子： 𝑑 =𝑠 𝑑𝑡
モビリティ (ダンパ) 8 ⚫ ダンパに①から②の方向に力ｆが加わると、 ⚫ ①－②間は圧縮される ① 𝑓 𝑎 ② 𝑓 = 𝑑 𝑣1 − 𝑣2 = 𝑑𝑣12 𝑣1,2 𝑣1 ⚫ その時の①－②間の速度差と力の関係は次となる 𝑣2 𝑣1,2 = 𝑎𝑓 d: 粘性抵抗係数 𝑁 𝑠/𝑚 モビリティの形式にすると 𝑣1,2 = 𝑎𝑓 𝑎: 𝑚𝑠 −1 N −1 𝜆 =𝑎 𝑎 = 1Τ𝑑 𝑚𝑠 −1 N −1 ダンパ (damper)とは、機械構造や建築物の振動を 減衰する装置です。減衰力の発生には、非圧縮性の 粘性流体抵抗などを利用します。ショックアブソー バー（shock absorber）とも呼ばれます。 Q: a=0とはどのようなことでしょう？
粘性抵抗と粘性係数の違い 9 By ChatGPT 粘性抵抗と粘性係数は、どちらも流体や機械の摩擦に関係する概念ですが、間違いやすいのでここに示します 粘性抵抗（Viscous Resistance） 粘性係数（Viscosity Coefficient, Dynamic Viscosity） 意味 流体や機械系の動きに対して生じる抵抗 力のこと 物質（特に流体）の粘り気の強さを示す 物理量 性質 物体の運動速度に比例して大きくなる 流体の内部摩擦を示し、温度によって変 化する 𝑑𝑢 𝜏=𝜂 𝑑𝑦 𝐹 = 𝑑𝑣 式 𝐹: 粘性抵抗 𝑁 𝑑: 粘性抵抗係数 𝑁 𝑠/𝑚 𝜏: せん断応力 𝑁Τ𝑚2 𝜂: 粘性係数 𝑃𝑎 ∙ 𝑠 = 𝑁𝑠Τ𝑚2 𝑑𝑢Τ𝑑𝑦 : 速度勾配 𝑣: 速度[𝑚Τ𝑠] 例 • • 油や水の中で動く物体が受ける抵抗 車のダンパー（ショックアブソーバー） の減衰 • • 水や油の流れやすさ 空気の粘性による摩擦
モビリティ (ばね) 10 ⚫ ばねに①から②の方向に力ｆが加わると ⚫ ①－②間は圧縮される ① 𝑐 ② 𝑓 𝑓 = 𝑘 𝑥1 − 𝑥2 𝑣1,2 𝑣1 𝑣2 𝑣1,2 = 𝑐𝑠 𝑓 𝜆 = 𝑐𝑠 𝑐 = 1Τ𝑘 𝑚/𝑁 𝑑 演算子： =𝑠 𝑑𝑡 ⚫ その時の①－②間の速度差と力の関係は次となる k: ばね剛性 N/𝑚 逆数での表現も一般的です 𝑥1 − 𝑥2 = 𝑐𝑓 c: コンプライアンス 𝑚/𝑁 𝑐 = 1Τ𝑘 ばねでは、速度よりも変位が一般的であるため、 変位で記述すると 𝑣1,2 = 𝑣1 − 𝑣2 = 𝑑 𝑥1 − 𝑥2 𝑑𝑡 𝑑 𝑣1,2 = 𝑐 𝑓 = 𝑐𝑠 𝑓 𝑑𝑡
モビリティ (質量) ⚫ 質量の入力端①の速度v1と出力端の速度v2 はどのよう に考えればよいでしょう ① ⚫ 質量と力の関係は f = ma となり、加速度a は慣性系 からの加速度となります。これは通常地球を原点とした 加速度と考えています 𝑓 m ⚫ すなわち、入力端①の速度をv1 とすると、出力端2は慣 性系である地球であり、v2 = 0 と考えてよいでしょう 𝑣1,0 𝑣1 𝑣2 = 0 𝑣1,0 = 𝑚𝑠 −1 𝑓 𝜆= 𝑚𝑠 −1 𝑚 𝑘𝑔 演算子： 𝑑 =𝑠 𝑑𝑡 𝑓 = ma = m 𝑑 𝑑 𝑣1 − 𝑣2 = m 𝑣 −0 𝑑𝑡 𝑑𝑡 1 𝑣1,0 = 𝑣1 − 0 = 1 1 න 𝑓 𝑑𝑡 = 𝑓 𝑚 𝑚𝑠 1 𝑣1,0 = 𝑓 = 𝑚𝑠 −1 𝑓 𝑚𝑠 11
1.3 機械回路と電気回路の対応関係 機械回路 電気回路 流通量 f: 要素を通る力 i: 要素を通る電流 位差量 v: 要素両端の速度差 v: 要素両端の電位差 ダンパと抵抗 𝜆=𝑎 𝑎 = 1Τ𝑑 d: 粘性抵抗係数 𝑁 𝑠/𝑚 バネとコイル 質量とコンデンサ 𝜆 = 𝑐𝑠 𝑐 = 1Τ𝑘 k:バネ定数[m/N] ｃ: コンプライアンス[N/m] 𝜆 = 𝑚𝑠 −1 m: 質量[kg] 𝑧=𝑅 R: レジスタンス[Ω] 𝑧 = 𝑗𝜔𝐿 L: ヘンリ[H] 𝑧 = 𝑗𝜔𝐶 −1 C: ファラッド[F] 𝑠 ⇒ 𝑗𝜔 12
1.4 モビリティの直列、並列接続 モビリティを直列または並列接続したときの合成モビリティについて示します （１）直列 (serial) （２）並列 (parallel) 𝑛 𝑛 𝜆 = 𝜆𝑖,𝑖+1 𝜆−1 = 𝜆−1 𝑖,𝑖+1 𝑖=1 𝑖=1 電気回路の抵抗の接続方法、直列と並列の場合の合成抵抗の算出方法と同じ 13
1.5 モビリティ (直列質量) 14 これまで質量は、入力端と出力端のもつ２端子素子と考えることはできないとしてきました。 しかし、ある機構を用いて要素を通過する力と要素間の速度差の2端子素子を実現できます 通過する力 速度差 𝑣1 𝑓 𝑣2 𝑣1,2 ① ② M=4m V P 𝑓 l M F 質量を無視できる2個の梃があり、左梃の中央Pを ピンで支え、質量Mに接続する。質量Mの速度は、 𝑉 = 𝑣1 − 𝑣2 Τ2 支点Pにかかる質点による反作用力Fは、 𝐹=𝑀 l 𝑑𝑉 1 𝑑 1 𝑑 = 𝑀 𝑣1 − 𝑣2 = 𝑀 𝑣1,2 𝑑𝑡 2 𝑑𝑡 2 𝑑𝑡 この反作用力につり合う入力は、 𝑣1,2 = 𝑚𝑠 −1 𝑓 このように直列質量が構成できます。但し、 機構の質量は4mとなっています。直列接 続すると入出力端からの見かけの質量は減 少しますが、機構の質量が軽くなるわけで はありません 𝑓 = 𝐹 Τ2 となる。よって、 𝑣1,2 = 2 2 4 𝐹= 2𝑓 = 𝑓 ≡ 𝑚𝑠 −1 𝑓 𝑀𝑠 𝑀𝑠 𝑀𝑠 この様に、①－②間に力が通過して、かつ速度差 がある2端子機構となります。これが直列電気容量 に類似する直列質量となります
1.5 モビリティ (直列質量) ⚫ 図は直列質量を直列接続したもので す。伝達する力fを変えない場合、速 度差v1,4はどうなるかを考察します 𝑣1,4 =4V 𝑣1,2 =2V 𝑣1 ① M=4m V 𝑓 P M l 𝑣2 ② 𝑣3 ③ 𝑓 𝑓 𝑣3,4 = 2𝑉 𝑣4 ④ M=4m V P l M l l 直列質量の直列接続 𝜆 = 1Τ2 𝑚𝑠 −1 これにより、直列接続の場合、見かけの質量は単 𝑓 ⚫ まず P点の速度Vとします。すると 伝達力はfとなり、そして速度差は、 次の様になります 𝑣1,2 = 𝑣3,4 = 2𝑉, デンサCの直列接続のときと同様な結果となりま 𝑣1,4 = 4𝑉 よって直列質量接続のモビリティーλは 𝜆= 𝑣1,4 𝑣1,2 + 𝑣3,4 = 𝑓 𝑓 また前頁で示したように 2 4 𝑣1,2 = 2𝑓 = 𝑓 𝑀𝑠 𝑀𝑠 以上より、λは次のように求まります 独のときと比べて1/２になります。これはコン す 15 𝜆= 𝑣1,2 + 𝑣3,4 8 2 = = 𝑓 𝑀𝑠 𝑚𝑠
見かけの質量 16 入力点からみた見かけの質量について計算してみよう 入力点 f 𝑀 = 4𝑚 1. 入力点にfを加えると、 P点では2fの力 が質量Mに作用する P M 2. また入力点の変位に対して、 P点での 変位（加速度）は1/2となる 2f 3. 即ち、 P点では2倍の力を加えたのに対 して1/2の変位の変化しかない 𝑓 P ただし、直列質量がこのような重力場 に置かれたら重さは4mですね 重 力 M M=4m 𝑎:加速度 𝑓 ② l 𝑔 l P点： 4. これは質量Mに対して、入力点からは 1/4の質量mと見かけ上なる ① 入力点： 𝑓 =m 𝑎 2𝑓 4𝑓 1 = 𝑎 = 4𝑚 𝑎 2
見かけの質量 17 ⚫ コンデンサの直列接続では合成容量は小さくなります。これを質量に当てはめて 考えると、質量を上手くつなぎ合わせれば軽くなることに相当します ⚫ このようなことができれば貨物輸送は楽になり、宇宙への進出も簡単でしょう。 電気回路では容量の直列接続で合成容量は減少します。ただし、“機械回路で直 列接続するとトータルの質量が軽くなる”は違います 直列質量装置の質量は M(4m)である 2mg 4mg P ① 4mg M M=4m リンク機構、ケースなどの 質量は無視できるとする l 𝑔 l M 重 力 ② M M
2. モビリティの機械回路例 18 簡単な振動回路の例 モビリティの機械回路例への適用として機械的蓄音機のモデルを示します。 （この章は飛ばしても良いです）
2.1 機械的蓄音機のモデル 19 図にはこの機械的蓄音機のモデルを示します。これは質量とダンパをバネを介し てクランクで駆動する回路です。レコード盤にある溝の凹凸がXecosωt とすると、 その時間微分ve が速度となります。この振動が、サウンドボックスとホーン内の 空気を動かし、ホーンより速度vの音波になって聞こえるモデルを表しています。 この時のfは針にかかる力に相当します
2.1 機械的蓄音機のモデル 20 レコード盤に溝が切ってあり、音声信号はその 溝のわずかな凹凸で記録してあります。 蓄音機の針がレコード盤の溝に触れるとレコー ド盤の回転に従って針が振動します。その振れ をサウウンドボックス、ホーンを通して、音を 再生するものです ホーン 弾性と質量と粘性でモデル化 エアチェンバー スタイラスバー トーンアーム レコード盤 溝 サウンド ボックス 針 レコード溝 機械的蓄音機 http://soundboxes.g2.xrea.com/audio2.html https://www.soundhouse.co.jp/contents/column/index?post=2257
2.2 LTspiceでの解析 21 LTspice 電流を計測(→力) 電圧を計測(→速度) ①電圧を計測(→速度) Cの値をリストの様に順次変化させる サウンドボックスとホーン内 の空気が振動してホーンより でる音波の速度ｖ ②電流を計測(→力) Lの値をリストの様に順次変化させる 青字はコメントアウトされていることを示す AC解析でBODE線図を描くコマンド レコード盤の溝をなぞる 針の受ける反力f
2.2 ①電圧を計測 20dB LTspiceでの解析 サウンドボックスとホーン内の空気が振動してホーン よりでる音波の速度ｖ 20dB 15m 10m 7m 5m 質量をパラメータ 弾性をパラメータ 3m 2m 1m 0.5m 0dB 22 0.01m 0.02m 0.03m 0.05m 0.07m 0.1m 0dB 0.3m 0.1m 0.5m 1m 20dB -20dB ホーンよりでる 音波速度v ホーンよりでる 音波速度v 1𝑚 ⟶ 1 × 10−3 𝑘𝑔 1𝑚 ⟶ 1 × 10−3 𝑚 ∙ 𝑁 −1 10Hz m：10-3 100Hz 1kHz 10kHz 10Hz 100Hz 1kHz 10kHz
2.2 ②電流を計測 LTspiceでの解析 レコード盤の溝をなぞる針の受ける反力 f 35dB 弾性をパラメータ 電流を計測(→力) 0.01m 0.02m 0.03m 0.05m 0.07m 0.1m 0dB 0.3m 0.5m レコード盤をなぞる 針にかかる反力 35dB 10Hz 1m Lの値(コンプライアンス)をリストの様に順次変化させる 1𝑚 ⟶ 1 × 10−3 𝑚 ∙ 𝑁 −1 100Hz 1kHz 10kHz 23
2.2 過渡応答と状態図 LTspiceでの解析 24 サウンドボックスとホーン内の空気が振動してホーン よりでる音波の速度ｖと振幅 速度 過渡応答図 インパルス応答 粘弾性抵抗 R1=2 積分器：速度を積分 して変位（振幅) impulse波形 粘弾性抵抗 速度 状態図 過渡応答解析コマンド 振幅
2.2 過渡応答と状態図 LTspiceでの解析 25 サウンドボックスとホーン内の空気が振動してホーン よりでる音波の速度ｖと振幅 粘弾性抵抗を変えた場合の過渡応答を示します 速度 過渡応答図 過渡応答図 R1=10 R1=1 状態図 状態図 速度 速度 速度
2.3 機械回路のモビリティ表現 26 機械回路をモビリティを用いて表現したものです。この機械回路のモビリティは， 𝜆 = 𝜆1 + 1 1 𝑎 = 𝑐𝑠 + = 𝑐𝑠 + 1 Τ 𝜆2 + 1 Τ 𝜆3 𝑚𝑠 + 𝑎−1 𝑎𝑚𝑠 + 1 蓄音機の針に発生する力fは、 𝑓= ホーンより速度vは 𝜐 = 𝑣𝑒 − 𝜆1 𝑓 = 1 − (10.10) (10.11) 𝑣𝑒 1 = 𝑎 𝜆 𝑐𝑠 + 𝑎𝑚𝑠 + 1 𝜆1 1 𝜆1 + 1 Τ 𝜆2 + 1 Τ 𝜆3 𝑣𝑒 = 1 1 1 1 𝑠2 + 𝑠+ 𝑓= 𝑠+ 𝑣 𝑎𝑚 𝑚𝑐 𝑐 𝑎𝑚 𝑒 1 1 1 𝑠 + 𝑠+ v = 𝑣𝑒 𝑎𝑚 𝑚𝑐 𝑐 2 (10.10) 1 𝑣 1 + 𝜆1 1Τ𝜆2 + 1Τ𝜆3 𝑒 (10.11) 𝑑 2 𝑓 𝑅 𝑑𝑓 𝑘 𝑑𝑣𝑒 𝑅 + + 𝑓 = 𝑘 + 𝑣𝑒 𝑑𝑡 2 𝑚 𝑑𝑡 𝑚 𝑑𝑡 𝑚 𝑑2 𝑣 𝑑𝑣 + 𝑅 + 𝑘𝑣 = 𝑘𝑣𝑒 𝑑𝑡 2 𝑑𝑡
2.3 機械回路のモビリティ表現 電気回路として表現したものです。この場合駆動電圧は、 となり，電圧（−ωXe）の正弦波交流となります。回路のインピーダンスは、 𝑉𝑒 = 𝑑 𝑋𝑒 cos 𝜔𝑡 = −𝜔𝑋𝑒 sin 𝜔𝑡 = 𝑅𝑒𝑎𝑙 −𝜔𝑋𝑒 𝑒 𝑗𝜔𝑡 𝑑𝑡 1 𝑅 1 − 𝜔2 𝐿𝐶 + 𝑗𝜔𝐿 𝑍 = 𝑗𝜔𝐿 + = 𝑗𝜔𝐶 + 𝑅 −1 1 + 𝑗𝜔𝐶𝑅 回路に流れる電流I は、 𝐼= 𝑉𝑒 𝑍 27
3. サイズモ系（seismic system）例 28 サイズモ系（seismic system）は、振動/変位する本体と、本体内部のばね要 素や減衰要素を介して支えられている質量要素からなる振動系です。 変位計や加速度計の原理として応用されています。 セ ン サ 本 体 m センサ出力 y 本体の振動/変位 を検出 k x d Ground GROUND 書籍（工学解析）では扱わず
3.1 計測原理（変位と加速度計測） サイズモ系 「変位」or 「加速度」センサとして使用します セ ン サ 本 体 m x+y y k x 29 d Ground （１）変位計測 センサ内部での ｍの変位 運動方程式： 𝑚 𝑥ሷ + 𝑦ሷ + 𝑐 𝑦ሶ + 𝑘𝑦 = 0 伝達関数： 𝑌(𝑠) 𝑠2 =− 2 𝑋(𝑠) 𝑠 + 2𝜁𝜔0 𝑠 + 𝜔02 𝜔0 : 固有振動数 𝜁： ダンピング比 𝜔0 = 𝑘 𝑑 ,𝜁 = 𝑚 2 𝑚𝑘 （２）加速度計測 𝑌(𝑗𝜔) ≅1 𝑋(𝑗𝜔) センサ本体の 変位 センサ内部での ｍの変位 𝑌(𝑗𝜔) 1 ≅ 𝜔 2 𝑋(𝑗𝜔) 𝜔02 センサが受けた 加速度 𝜔0 ≪ 𝜔 応用：地震計など （一定） 𝜔0 ≫ 𝜔 応用：加速度計など ✓ センサ内部での質 量ｍの変位とセン サが受けた加速度 の比は一定となる ✓ 質量ｍの変位ｙを 計測することで加 速度が計測できる
3.1 計測原理 自由振動： 𝑑 2 𝑥(𝑡) 𝑑𝑥(𝑡) 𝑚 +𝑑 + 𝑘𝑥(𝑡) = 0 𝑑𝑡 2 𝑑𝑡 簡略表現： 𝑑 2 𝑥(𝑡) 𝑑𝑥(𝑡) + 2𝜁𝜔 + 𝜔02 𝑥(𝑡) = 0 0 2 𝑑𝑡 𝑑𝑡 固有角振動数： 𝜔0 = 𝑘Τ𝑚 (natural angular frequency) 臨界粘性減衰係数： (critical viscous damping constant) 減衰比： (damping ratio) 振動peak比： 30 𝑑𝑐 = 2 𝑚𝑘 𝑑 𝑑 𝜁= = 𝑑𝑐 2 𝑚𝑘 𝑒 −𝜁𝜔0 Δ𝑡 セ ン サ 本 体
3.1 計測原理（変位と加速度の計測） 𝑠2 （１） 𝑌(𝑠) =− 2 変 𝑋(𝑠) 𝑠 + 2𝜁𝜔0 𝑠 + 𝜔02 位 計 𝑌(𝑗𝜔) 𝜔2 測 𝑋(𝑗𝜔) = 𝜔2 − 𝜔 2 + 2𝑗𝜁𝜔 𝜔 0 0 𝑌(𝑗𝜔) = 𝑋(𝑗𝜔) 31 セ ン サ 本 体 センサ出力 ：伝達関数 本体の振動/ 変位を検出 ：s=jω 周波数領域へ変換 𝜔2 𝜔02 − 𝜔 2 2 + 2𝑗𝜁𝜔0 𝜔 2 = 1 𝜔02 Τ𝜔 2 − 1 2 + 2𝑗𝜁 𝜔0 Τ𝜔 2 ≅𝟏 𝜔0 ≪ 𝜔 （y の値が、センサ本体の変位x となる） 𝑌(𝑠) 1 （２） 𝑌(𝑠) = 2 =− 2 ：伝達関数 加 ሷ 𝑠 𝑋(𝑠) 𝑠 + 2𝜁𝜔0 𝑠 + 𝜔02 𝑋(𝑠) 速 𝑌(𝑗𝜔) 1 度 = ：s=jω 周波数領域へ変換 計 𝜔 2 𝑋(𝑗𝜔) 𝜔02 − 𝜔 2 + 2𝑗𝜁𝜔0 𝜔 測 𝑌(𝑗𝜔) 1 1 𝟏 = = ≅ 𝜔 2 𝑋(𝑗𝜔) 𝝎𝟐𝟎 2 2 2 2 2 2 2 2 2 𝜔0 − 𝜔 + 2𝑗𝜁𝜔0 𝜔 𝜔0 1 − 𝜔 Τ𝜔0 + 2𝑗𝜁 𝜔Τ𝜔0 𝜔0 ≫ 𝜔
3.2 （１）変位センサ ( ω0<<ω ) 変位計測 セ ン サ 本 体 センサ本体の変位ｘはセンサ出力 ｙから計測できる 𝑌(𝑠) 𝑠2 =− 2 𝑋(𝑠) 𝑠 + 2𝜁𝜔0 𝑠 + 𝜔02 32 センサ内部で のｍの変位 センサ本体の 変位 𝑌(𝑗𝜔) ≅1 𝑋(𝑗𝜔) 𝜔0 ≪ 𝜔 変位センサになっている ω0=1, ω=200π x 積分 y 積分 C=0.1F, R=0.5Ω, L=10H 伝達関数からも分かるように符号は反転 (x+y)-x
3.2 （１）変位センサ 変位計測 33 ( ω0<<ω ) 前頁変位センサ例では、機械回路mdkの値は適当に選んだものです。 但し、減衰率は臨界減衰としています 機械系のmdkから電気回路のCRLへの対応は次の通りです C=m, R=1/d, L=1/k 機械回路 電気回路 前頁加速度 センサ例 m=0.1kg, d=2N s/m, k=0.1N/m C=0.1F, R=0.5Ω, L=10H 固有角振動数： 臨界粘性減衰係数： 減衰比： 𝜔0 = ( m=100g, d=2N s/m, k=0.0001N/mm ) 𝑘Τ𝑚 = 0. 1Τ0.1 = 1rad/s 𝑑𝑐 = 2 𝑚𝑘 = 2 0.1 ∙ 10 = 2N s/m 𝜁= 𝑑 2 = = 1.0 𝑑𝑐 2 臨界減衰 セ ン サ 本 体
3.2 変位計測（ダンパ係数の違い） （１）変位センサ ( ω0<<ω ) セ ン サ 本 体 減衰率は臨界減衰より小さくして いくとy出力が変動する センサ内部で のｍの変位 34 𝑌(𝑗𝜔) ≅1 𝑋(𝑗𝜔) 𝜔0 ≪ 𝜔 センサ本体の 変位 変位センサになっている x ω0=1, ω=200π 積分 100Ω (d=0.01N s/m) 10Ω (d=0.1N s/m) y 積分 0.5Ω (d=2N s/m) (x+y)-x
3.2 変位計測（ダンパ係数の違い） （１）変位センサ ( ω0<<ω ) セ ン サ 本 体 R (粘弾性抵抗d)をパラメータと してBODE線図を描く センサ内部で のｍの変位 センサ本体の 変位 35 𝑌(𝑗𝜔) ≅1 𝑋(𝑗𝜔) 𝜔0 ≪ 𝜔 変位センサになっている ω0=1, ω=200π 100Ω (d=0.01N s/m) 10Ω (d=0.1N s/m) 0.5Ω (d=2N s/m) 積分 臨界減衰 積分 BODE線図 0.1Hz C=0.1F, L=10H 1k (x+y)-x
3.3 加速度計測 （２）加速度センサ ( ω0>>ω ) セ ン サ 本 体 センサ本体の加速度はセンサ出力 ｙから計測できる センサ内部で のｍの変位 センサが受け た加速度 36 𝑌(𝑗𝜔) 1 ≅ 2（一定） 𝜔 2 𝑋(𝑗𝜔) 𝜔0 𝜔0 ≫ 𝜔 加速度センサになっている 入力速度 (dx) ω0=10k, ω=10π 積分 積分 cos 106・{(x+y)-x} y出力 (y) 𝑦∝ 𝑑 (𝑑𝑥) 𝑑𝑡 臨界減衰 sin L=0.1mH, C=0.1mF, R=2Ω 200ms
3.3 加速度計測 37 （２）加速度センサ ( ω0>>ω ) 前頁加速度センサ例では、機械回路mdkの値は適当に選んだものです。 但し、減衰率は臨界減衰としています 機械系のmdkから電気回路のCRLへの対応は次の通りです C=m, R=1/d, L=1/k 前頁加速度 センサ例 電気回路 機械回路 C=0.1mF, R=2Ω, L=0.1mH m=10-4kg, d=0.5N s/m, k=104N/m ( m=0.1g, d=0.5N s/m, k=10N/mm ) 固有角振動数： 臨界粘性減衰係数： 減衰比： 𝜔0 = 𝑘Τ𝑚 = 104Τ10−4 = 104rad/s 𝑑𝑐 = 2 𝑚𝑘 = 2 104 ∙ 10−4 = 2N s/m 𝑑 1 𝜁= = = 0.5 𝑑𝑐 2 セ ン サ 本 体
3.3 加速度計測 （２）加速度センサ ( ω0>>ω ) 減衰率は臨界減衰より小さく セ ン サ 本 体 センサ内部で のｍの変位 センサが受けた 加速度 していくと振動が残る 振動が残る 38 𝑌(𝑗𝜔) 1 ≅ （一定） 𝜔 2 𝑋(𝑗𝜔) 𝜔02 𝜔0 ≫ 𝜔 ω0=10k, ω=10π ω0=10k, ω=10π R=10Ω R=100Ω d=0.1N s/m d=0.01N s/m 振動が残る
3.3 加速度計測 （２）加速度センサ ( ω0>>ω ) セ ン サ 本 体 R (粘弾性抵抗d)をパラメータ としてBODE線図を描く センサ内部で のｍの変位 センサが受けた 加速度 39 𝑌(𝑗𝜔) 1 ≅ （一定） 𝜔 2 𝑋(𝑗𝜔) 𝜔02 𝜔0 ≫ 𝜔 加速度センサになっている 100Ω (d=0.01N s/m) 10Ω (d=0.1N s/m) 0.5Ω (d=2N s/m) 臨界減衰 BODE線図 1Hz L=0.1mH, C=0.1mF 10k ω0=10k
4. 動吸振器（dynamic vibration absorber） 40 対象物に補助的な質量体を取り付け、この質量体に対象物の振動を吸収さ せて代わりに振動させることで対象物の振動抑制を図るのが動吸振器です。 建造物における設計思想の「制振」に分類される装置に該当します https://ja.wikipedia.org/wiki/動吸振器 補助的な 質量体 b m k2 f x2 2 d a m1 x1 k1 書籍（工学解析）では扱わず
4.1 動吸振器の運動方程式 41 減衰付動吸振器の運動方程式 b m k2 f x2 2 𝑑2 𝑥2 (𝑡) 𝑑𝑥1 𝑡 𝑑𝑥2 (𝑡) 𝑚2 + 𝑑 − + − 𝑘2 𝑥1 𝑡 + 𝑘2 𝑥2 (𝑡) = 0 𝑑𝑡 2 𝑑𝑡 𝑑𝑡 d a m1 𝑑2 𝑥1 (𝑡) 𝑑𝑥1 (𝑡) 𝑑𝑥2 (𝑡) 𝑚1 + 𝑑 − + 𝑘1 + 𝑘2 𝑥1 𝑡 − 𝑘2 𝑥2 (𝑡) = 𝑓(𝑡) 𝑑𝑡 2 𝑑𝑡 𝑑𝑡 x1 k1 主系m1の振動を 従系m2が振動す ることで吸収する 行列形式で表示 𝑚1 0 0 𝑥ሷ 1 𝑑 + 𝑚2 𝑥ሷ 2 −𝑑 𝑘 + 𝑘2 −𝑑 𝑥ሶ 1 + 1 −𝑘2 𝑑 𝑥ሶ 2 𝜇 = 𝑚2 Τ𝑚1 ：主系と従系の質量の比 𝑝1 = 𝑘1 Τ𝑚1 𝑝2 = 𝑘2 Τ𝑚2 −𝑘2 𝑥1 𝑓 = 𝑘2 𝑥2 0 ：主系の固有振動数 ：従系の固有振動数 ：主振動系と副振動系の固有角振動数比 𝜅 = 𝑝2 Τ𝑝1 (カッパー) ：主振動系の固有角振動数を用いた減衰比 𝜁 = 𝑑Τ2𝑚2 𝑝1 https://ja.wikipedia.org/wiki/動吸振器 ：副振動系の減衰比 https://www.youtube.com/watch?v=KfK_PjZ4MQ0 分かり易かった
4.2 電気回路LTspiceでの解析 42 機械モデルを電気回路に変換する方法 b 補助的な 質量体 b m f x2 2 k2 C2 R d a 周期的外力 m1 L2 a x1 f C=m, R=1/d, L=1/k k1 1. 質量m をキャパシタンスC に置き換え配置を行う 2. キャパシタンスC の一方の端子はグランドに接続する 3. キャパシタンスC の他の端子はノード記号を付ける（a,b,c... など） 4. C1 L1 バネ要素k，粘性抵抗要素dをそれぞれ，インダクタンスL， 抵抗Rに変更する 5. バネ要素，粘性抵抗要素は，上記ノード間，ノードとグラン ド間を接続しているから，それらを結ぶように配置する 6. 外力(もしくは強制変位速度) は接続するノードに，力（もし くは速度）に対応する電流源(もしくは電圧源) を接続する R a L2 f L1 C1 C2 b
4.2 電気回路LTspiceでの解析 43 周波数特性の解析 機械モデルを変換した電気回路 LTspice パラメータ m1の変位x1 C1=C2=L1=L2=1 R=0.2, 0.4,2,4,20,40 m1=m2=k1=k2=1 d=5, 2.5, 0.5,0.25,0.05,0.025 積分 m2の変位x2 パラメータの選定は一例 （単純な場合を想定） 積分 振動系の減衰比 𝜁 = 𝑑Τ2𝑚2 𝑝 , d = 1Τ𝑅 振動系の減衰比をパラメータとする 𝑚1 = 𝑚2 = 𝑘1 = 𝑘2 = 𝑑 = 1
4.3 動吸振器(AC解析：定点理論) 15V m1応答(x1) １次 固有振動 質量比μと振動数比κを 決めると、どんなζでも P点Q点を通る （ 定点理論） パラメータ ２次 固有振動 P 44 C1=C2=L1=L2=1 R=0.2, 0.4,2,4,20,40 Q m1=m2=k1=k2=1 d=5, 2.5, 0.5,0.25,0.05,0.025 0V 15V m2応答(x2) 𝜁 = 𝑑Τ2𝑚2 𝑝 主振動系の固有角振動数を用い た減衰比 𝑝1 = 𝑝1 = 𝑘1 Τ𝑚1 1Τ1 = 1 𝑟𝑎𝑑 Τ𝑠 = 0.159𝐻𝑧 図の縦軸は変位を示す 0V 40mHz 160mHz 220mHz 400mHz https://deltapower.site/2020/07/20/vibrati on-dynamic-damper-fixed-point-theory/
4.4 動吸振器（最適化例） 45 最適化例：2つの共振ピーク の高さが同じようにする 𝜁 = 𝑑Τ2𝑚2 𝑝1 ⇒ 𝜁𝑜𝑝𝑡 = 1 2 1+𝜇 3𝜇 2 1+𝜇 𝜅 = 𝑝2 Τ𝑝1 ⇒ 1 副振動系 𝜅𝑜𝑝𝑡 = の減衰比 1+𝜇 𝜁𝑜𝑝𝑡 = 0.2165, ⇒ 𝑑 = 0.4330127 𝜅𝑜𝑝𝑡 = 0.5, ⇒ 𝑘2 = 0.25 𝜁 = 𝑑Τ2𝑚2 𝑝1 𝑚1 = 𝑚2 = 𝑘1 = 1 𝑘2 = 0.25 𝑑 = 0.4330127 𝑝1 = 𝑘1 Τ𝑚1 𝑝2 = 𝑘2 Τ𝑚2 𝜇 = 𝑚2 Τ𝑚1
4.4 動吸振器（最適化例） m1応答(x1) 4.925 46 最適化例：2つの共振ピーク の高さが同じようにする 𝜁𝑜𝑝𝑡 = 0.2165 ⇒ 𝑑 = 0.4330127 𝜁𝑜𝑝𝑡 = 0.2165 0.4341 2.4626 𝜅𝑜𝑝𝑡 = 0.5 ⇒ 𝑘2 = 0.25 m2応答(x2) 振動系の減衰比 𝜁 = 𝑑Τ2𝑚2 𝑝 , d = 1Τ𝑅 https://www.youtube.com/watch?v=KfK_PjZ4MQ0 分かり易かった
5. クランク軸のねじり振動 この図はエンジンのフライホイールとクランク軸とプロペラ軸を 表しています。この振動系の固有各振動数を求めてみます Torque 1 I 2 k I 3 k I 4 k I 5 k I 6 k クランク軸（crankshaft） pythonで記述したprogramは下記GitHubにあります https://github.com/m4881shimojo/TTanaka_Eng/ I kd Id 47
5. 概 クランク軸のねじり振動 48 要 1. 解析的手法 step1: 階差方程式を立てる step2: 固有方程式と固有を得る step3: 境界条件を満たす解を求める step4: 例題での計算例を示す 2. シミュレーター 電気回路LTspiceでの解析 境界条件 （左） 境界条件 （右）
5.1 運動方程式 49 この図はエンジンのフライホイールとクランク軸とプロペラ軸を表しています 1. 慣性モーメントは、ピストン機構バランス重りを考慮した平均値をI、ねじり コンプライアンスをcとしています 2. シリンダ数をnとして、境界条件を決める両端とx番目と(x+1)番目の軸部分 を示しています 3. クランクのねじれモビリティーをλc、慣性モビリティーをλmとしています。 またx番目と(x+1)番目の軸を捩じるトルクをτxとしています 1 Ig cg I 1 𝜆𝑔 2 c 𝜆𝑐 I x 2 x 𝜏2 𝜏𝑔 I c 左端 c n I 𝜆𝑐 Id cd n 𝜆𝑑 𝜏𝑥+1 𝜏𝑥 𝜆𝑚 I c 𝜆𝑐 x+1 𝜏𝑥−1 𝜏1 𝜆𝑚 x+1 𝜆𝑚 𝜏𝑛−1 𝜆𝑚 中間 𝜆𝑚 𝜏𝑑 𝜆𝑚 右端
5.1 階差方程式 運動方程式（階差方程式） 50 x x番目の閉回路での角速度差を求める step1: x点の角速度： xと(x+1)の間の角速度差： c 𝜆𝑚 𝜏𝑥 − 𝜏𝑥−1 I 𝜏𝑥−1 𝜆𝑐 𝜏 𝑥 x 𝜏𝑥−1 𝜆𝑚 𝜆𝑚 𝜏𝑥+1 − 𝜏𝑥 − 𝜆𝑚 𝜏𝑥 − 𝜏𝑥−1 = 𝜆𝑐 𝜏𝑥 整理すると 線形階差方程式： 𝜏𝑥−1 − 2 + 𝜆𝑐 𝜏 + 𝜏𝑥+1 = 0 𝜆𝑚 𝑥 (10.46) I 𝜏𝑥 (x+1)点の角速度： 𝜆𝑚 𝜏𝑥+1 − 𝜏𝑥 これから x+1 𝜆𝑐 𝜏𝑥+1 x+1 𝜏𝑥 𝜏𝑥+1 𝜆𝑚 慣性モーメント のモビリティ ∆𝜃ሶ = I𝑠 −1 𝜏 𝜆𝑚 = I𝑠 −1 ここで上式の解を次のようにします。A,Bは境界条件で定める定数、 αは階差方程式を満足する固有値（定数）です 𝜏𝑥 = 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 （クランク軸への入力を正弦波的としている？） (10.47) ねじれコンプライア ンストのモビリティ ∆𝜃ሶ = cs 𝜏 𝜆𝑐 = cs 正弦波的解析の場合：s →jω
5.1 運動方程式（階差方程式） 51 step2: 線形階差方程式： τxに想定する解： 𝜏𝑥−1 − 2 + → 𝜆𝑐 𝜏 + 𝜏𝑥+1 = 0 𝜆𝑚 𝑥 𝜏𝑥 = 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 𝐴 cos 𝑥 − 1 𝛼 + 𝐵 sin 𝑥 − 1 𝛼 − 2 + 𝜆𝑐 𝜆𝑚 これを整理すると 𝑐𝑜𝑠𝛼 − 1 + (10.46) を代入します 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 + 𝐴 cos 𝑥 + 1 𝛼 + 𝐵 sin 𝑥 + 1 𝛼 = 0 1 𝜆𝑐 2 𝜆𝑚 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 = 0 常にこの式が成り立つためには、 次の式が成立する必要があります 固有方程式 (eigen equation) 𝑐𝑜𝑠𝛼 = 1 + 1 𝜆𝑐 2 𝜆𝑚 α：固有値 これが線形階差方程式の固有方程式で、また根αは固有値と呼ばれます
5.1 step3: 運動方程式（左:境界条件） 52 固有方程式の解で境界条件を満たす解を求める 1 𝜆𝑔 境界条件（BC)： 𝜆𝑐 2 1 𝜏2 𝜆𝑔 Mobility： 𝜆𝑔 = 𝑗𝜔𝑐𝑔 + 𝑗𝜔𝐼 𝜏𝑔 𝑔 𝜏𝑔 の閉回路： 𝜆𝑔 𝜏𝑔 + 𝜆𝑚 𝜏𝑔 − 𝜏1 = 0 𝜏1 左端境界条件 𝜏1 の閉回路： 𝜆𝑚 𝜏1 − 𝜏𝑔 + 𝜆𝑐 𝜏1 + 𝜆𝑚 𝜏1 − 𝜏2 = 0 両式からτgを消去： − 2+ 𝜆𝑐 𝜆𝑚 + 𝜏 + 𝜏2 = 0 𝜆𝑚 𝜆𝑚 + 𝜆𝑔 1 固有方程式 2𝑐𝑜𝑠𝛼 = 2 + − 2𝑐𝑜𝑠𝛼 + 解を右記とすると 左端の 境界条件 𝜆𝑚 𝜏 + 𝜏2 = 0 𝜆𝑚 + 𝜆𝑔 1 𝜏𝑥 = 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 𝐴 𝜆𝑐 = 𝑗𝜔𝑐 𝜆𝑐 𝜆𝑚 𝜆𝑚 = 𝑗𝜔𝐼 −1 これを上式に代入してまとめる 𝜆𝑚 𝜆𝑚 cos 𝛼 − 1 + 𝐵 sin 𝛼 = 0 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 + 𝜆𝑔 (10.51)
5.1 step3: 運動方程式（右:境界条件） 固有方程式の解で境界条件を満たす解を求める 53 𝜆𝑐 n 𝜆𝑑 境界条件（BC)： 𝜆𝑑 Mobility： 𝜆𝑑 = 𝑗𝜔𝑐𝑑 + 1 𝑗𝜔𝐼𝑑 𝜏𝑛−1 𝜏d の閉回路： 𝜆𝑑 𝜏𝑑 + 𝜆𝑚 𝜏𝑑 − 𝜏𝑛−1 = 0 𝜏𝑑 右端境界条件 𝜏1 の閉回路： 𝜆𝑚 𝜏𝑛−1 − 𝜏𝑛−2 + 𝜆𝑐 𝜏𝑛−1 + 𝜆𝑚 𝜏𝑛−1 − 𝜏𝑑 = 0 境界条件からτdを消去 固有方程式を用いて整理 解を右記とすると 右端の 境界条件 𝐴 cos 𝑛𝛼 − 𝜆𝑚 − 2𝑐𝑜𝑠𝛼 𝜏𝑛−1 + 𝜏𝑛 = 0 𝜆𝑚 + 𝜆𝑑 𝜏𝑥 = 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 これを上式に代入してまとめる 𝜆𝑚 𝜆𝑚 cos 𝑛 − 1 𝛼 + 𝐵 sin 𝑛𝛼 − sin 𝑛 − 1 𝛼 = 0 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 + 𝜆𝑔 (10.52)
5.2 step3: 階差方程式の解 54 固有方程式の解で右＆左の境界条件を満たす解(10.51),(10.52)を求めた。両式には、 A,B項以外の定数項はないから、A,Bの係数からなる以下の行列式はゼロとなる 𝜆𝑚 cos 𝛼 − 1 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 cos 𝑛𝛼 − cos 𝑛 − 1 𝛼 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 sin 𝛼 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 sin 𝑛𝛼 − sin 𝑛 − 1 𝛼 𝜆𝑚 + 𝜆𝑔 =0 これを解くと（境界条件を満たす解は次の様になる） 𝜆𝑚 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 𝜆𝑚 𝜆𝑚 sin 𝑛 − 2 𝛼 − + sin 𝑛 − 1 𝛼 + sin 𝑛𝛼 = 0 𝜆𝑚 + 𝜆𝑑 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 + 𝜆𝑑 ただし 1 𝜆𝑐 𝑐𝐼 𝑐𝑜𝑠𝛼 = 1 + ≡ 1 − 𝜔2 2 𝜆𝑚 2 (10.53) 𝜆𝑐 = 𝑗𝜔𝑐 固有方程式 𝜆𝑚 = 𝑗𝜔𝐼 −1 これより固有振動数を求めるには原理的にはλはすべてωの関数であるから、上式 (10.53)(10.49)よりαを消去して、ωだけの固有方程式とし F 𝜔 =0 の根、ω1, ω2…を求める。これが全クランク軸のねじれ固有角周波数となる (10.49)
階差方程式と差分方程式 55 By ChatGPT 階差方程式と差分方程式の違いは？ ⚫ 「階差方程式」は、特に隣り合う項の差（階差）を扱うシンプルな形の方程式。 ⚫ 「差分方程式」は、より広い概念であり、高次の項を含む一般的な形式の方程 式も指す。 ⚫ つまり、階差方程式 ⊂ 差分方程式 という関係になります。 項目 階差方程式 差分方程式 定義 隣接する項の差（階差）を利用した 方程式 広く数列の項どうしの関係を表す 離散方程式 例 𝑥𝑛+1 − 𝑥𝑛 = 3 𝑥𝑛+2 + 𝑝𝑥𝑛+1 + 𝑞𝑥𝑛 = 0 特徴 1次の関係が多い 高次の項を含む場合もある 関係性 差分方程式の一部 階差方程式を含む、より一般的な 概念
5.3 クランク軸のねじり振動(例題) step4: 例題 ６シリンダ・エンジンにプロペラ軸を接続 d 1 2 3 4 5 6 左端境界： 𝐼𝑔 = 0, 𝜆𝑔 = ∞ 右端境界： 𝜆𝑑 = 𝑗𝜔𝑐𝑑 + 𝑗𝜔𝐼𝑑 −1 端以外： 𝜆𝑐 = 𝑗𝜔𝑐, 𝜆𝑚 = 𝑗𝜔𝐼 −1 解(10.53)式に代入して 整理すると sin 𝑛𝛼 − c I c I c I c I c I cd I 𝜆𝑚 sin 𝑛 − 1 𝛼 = 0 𝜆𝑚 + 𝜆𝑑 αを直接解く代わりに新しい変数（θ, μ）を導入 𝜇= 56 Id (10.54) 𝜃 = 𝑛−1 𝛼 (10.55) 𝜆𝑚 1 1 1 + = = 𝜆𝑚 + 𝜆𝑑 1 + 𝜆𝑑 Τ𝜆𝑚 1 + 𝐼 Τ𝐼𝑑 − 𝜔 2 𝑐𝑑 𝐼 1 + 𝐼 Τ𝐼𝑑 − 2 𝑐𝑑 Τ𝑐 1 − cos 𝛼 (10.56) (10.54)に(10.55),(10.56) を代入すると sin 𝜃 + 𝛼 − 𝜇 sin 𝜃 = 0 分解 cos 𝛼 − 𝜇 sin 𝜃 + sin 𝛼 cos 𝜃 = 0 cos 𝛼 − 𝜇 2 + sin 𝛼 2 sin 𝜃 + tan−1 sin 𝛼 =0 cos 𝛼 − 𝜇 (a)
5.3 step4: クランク軸のねじり振動(例題の解) 57 前頁(a)がゼロとなる条件：μ､α は実数であり、√ の所がゼロとなることはない よって(a)がゼロとなるためには、次式が成り立つ sin 𝜃 + tan−1 sin 𝛼 =0 cos 𝛼 − 𝜇 (10.57) θ、μを元に戻す sin 𝑛 − 1 𝛼 + tan sin 𝛼 −1 cos 𝛼 − 例題の解 1 + 𝐼 Τ𝐼𝑑 1 − 2 𝑐𝑑 Τ𝑐 1 − cos 𝛼 ｛ ｝内をf(α)とすれば =0 sin 𝛼 𝑓 𝛼 = 𝑛 − 1 𝛼 + tan−1 cos 𝛼 − (10.58)の解は、f(α)=mπ/2 1 + 𝐼 Τ𝐼𝑑 1 − 2 𝑐𝑑 Τ𝑐 1 − cos 𝛼 (m整数)となるαを求めればよい (10.58)
5.3 step4: クランク軸のねじり振動(例題) 58 【計算例】 書籍では重力単位系であったのでSI単位系に変換します 重力単位系： I=300 cm･kg･s2, Id=20 cm･kg･s2, c=6×10-9 rad･cm-1･kg-1, cd=8×10-8 rad･cm-1･kg-1 国際単位系: (SI) I=７7 kg･m2 ,Id=0.2xg kg･m2 ,c=6×10-7×(1/g)rad･kg-1･m-2･s2, cd= 8×10-6×(1/g)rad･kg-1･m-2･s g=9.8066 [m･s-2] 例題の解(10.58)に値を代入 sin 𝛼 𝑓 𝛼 = 6 − 1 𝛼 + tan−1 1 cos 𝛼 − 1 + 3𝑔Τ0.2𝑔 − 2 8 × 10−6 𝑔Τ6 × 10−7 𝑔 1 − cos 𝛼 sin 𝛼 𝑓 𝛼 = 5𝛼 + tan−1 cos 𝛼 − 1 16 − 26.667 1 − cos 𝛼 (10.58)の解は、f(α)=mπ (10.59) αは固有方程式(10.49)の解 (m整数)となるαを求めればよい また解αから次の様にクランク－プロペラ軸系の共振周波数が求められます 𝑓𝑝 = 𝜔𝑝 1 2 = 1 − cos 𝛼 2𝜋 2𝜋 𝑐𝐼 (10.60)
5.4 クランク軸のねじり振動(例題)の解 step4: f(α) cos 𝛼 − 1 16 − 26.667 1 − cos 𝛼 クランク－プロペラ軸系の共振周波数 900 6 720 5 540 4 360 2 59 f(α)= π, 2π, 3π, 4π, 5π, 6πとなる αの値を図から求める sin 𝛼 𝑓 𝛼 = 5𝛼 + tan−1 f(α)= mπ p αp(degree) fp(Hz) 1 29.6 60.61 2 57.5 114.12 3 (67.6) 131.98 4 90.8 168.93 5 120.3 205.78 6 150.2 229.28 3 180 1 f(α)=mπ (degree) 𝑓𝑝 = 𝜔𝑝 1 2 = 1 − cos 𝛼 2𝜋 2𝜋 𝑐𝐼 p310v3.py https://github.com/m4881shimojo/TTanaka_Eng/tree/main/crankshaft
クランク軸のねじり振動(例題)の解 f(α)= mπ 拡大図 p310v3.py program で表示範 囲とstep を調整 360 180 29.6 57.5 67.6 720 540 900 120.3 90.8 150.2 60
5.4 step4: クランク軸のねじり振動(例題)の解 f(α)= mπ 61 つぎに各共振周波数におけるねじれトルクの分布状態を求めてみる 1 クランク軸右端の境界条件： 𝐴 cos 𝑛𝛼 − 𝜏1 3 𝜏2 4 𝜏3 5 𝜏4 と置く 𝜏𝑥 = B 𝜏𝑥 = 𝐴 cos 𝛼𝑥 + 𝐵 sin 𝛼𝑥 d 6 𝜏5 𝜆𝑚 𝜆𝑚 cos 𝑛 − 1 𝛼 + 𝐵 sin 𝑛𝛼 − sin 𝑛 − 1 𝛼 = 0 𝜆𝑚 + 𝜆𝑔 𝜆𝑚 + 𝜆𝑔 クランク軸のねじりトルク： 𝐴 =𝜙 𝛼 𝐵 2 𝜏6 (10.52) (10.47) 𝜏𝑥 ⟹ 𝜏1 , 𝜏2 , ⋯ , 𝜏6 𝐴 cos 𝛼𝑥 + sin 𝛼𝑥 = B 𝜙 𝛼 cos 𝛼𝑥 + sin 𝛼𝑥 𝐵 sin 𝑛 − 1 𝛼 1 + 𝐼 Τ𝐼𝑑 − 2 𝑐𝑑 Τ𝑐 1 − cos 𝛼 𝜏𝑥 = B sin 𝛼 𝑥 − cos 𝑛 − 1 𝛼 cos 𝛼 − 1 + 𝐼 Τ𝐼𝑑 − 2 𝑐𝑑 Τ𝑐 1 − cos 𝛼 sin 𝛼 − (10.61)
5.4 step4: クランク軸のねじり振動(例題)の解 1 つぎに各共振周波数におけるねじれ トルクの分布状態 2 3 𝜏1 𝜏2 4 𝜏3 5 𝜏4 62 d 6 𝜏5 𝜏6 下の図は第五モード(P5, 205.78Hz)でのクランク軸にかかるトルク状態の分布状態を示します P=5, f5= 205.78Hz 1.5 1 0.5 0 1 2 3 4 5 6 𝜏1 𝜏2 𝜏3 𝜏4 𝜏5 𝜏6 -0.5 -1 -1.5 p311v2.py 次頁には、P1～P5の共振振動モードでのクランク軸にかかるトルク状態の分布状態を示します https://github.com/m4881shimojo/TTanaka_Eng/tree/main/crankshaft
5.4 step4: クランク軸のねじり振動(例題) P=1, f1=60.60Hz P=2, f1=114.11Hz 1.5 1.5 0.5 0.5 -0.5 1 2 3 4 5 6 -1.5 -0.5 1 1.5 0.5 0.5 1 2 3 4 5 6 -1.5 -0.5 1 4 5 6 2 3 4 5 6 5 6 -1.5 P=5, f5= 205.78Hz P=6, f6=229.28Hz 1.5 1.5 0.5 0.5 -1.5 3 P=4, f4=168.93z 1.5 -0.5 2 -1.5 P=3, f1=131.98Hz -0.5 63 1 2 3 4 5 6 -0.5 -1.5 1 2 3 4
5.5 電気回路LTspiceでの解析 クランク軸のねじり振動(例題) クランク軸（crankshaft） 電気回路モデルへ変換 LTspiceで解析 64
5.5 電気回路LTspiceでの解析(BODE線図) LTspice BODE線図 65 クランクープロペラ系の振動系 の電気回路 トルク（電流）周波数応答結果 （BODE線図) AC解析 機械系のパラメータ I=3/g= 0.3059kg･m2 Id=0.2/g= 0.02039 kg･m2 c= 6×10-7×(g)rad･kg-1･m-2･s2 = 5.884×10-6 rad･kg-1･m-2･s2 cd= 8×10-6×(g)rad･kg-1･m-2･s2 = 7.845×10-5 rad･kg-1･m-2･s2
5.5 電気回路LTspiceでの解析(BODE線図) 66 電気回路シミュレータ結果 BODE線図 LTspice トルク（電流）周波数応答結果 （BODE線図) トルク(電流) 20dB τ1 τ1 τ2 τ6 ねじり振動の共振周波数 10dB τ6 0dB -10dB 10Hz τ2 60Hz 100Hz 200Hz 1kHz p αp(°) fp(Hz) 1 29.6 60.61 2 57.5 114.12 3 (67.6) 131.98 4 90.8 168.93 5 120.3 205.78 6 150.2 229.28 右図はτ1,τ2,τ6(torque) を示します
5.5 LTspice 電気回路LTspiceでの解析(過渡応答) Transient線図 67 クランクープロペラ系 振動系の電気回路 トルク入力での 伝達トルク応答結果 （ステップ応答) （トルクは電流となる） 過渡解析 機械系のパラメータ I=3/g= 0.3059kg･m2 Id=0.2/g= 0.02039 kg･m2 c= 6×10-7×(g)rad･kg-1･m-2･s2 = 5.884×10-6 rad･kg-1･m-2･s2 cd= 8×10-6×(g)rad･kg-1･m-2･s2 = 7.845×10-5 rad･kg-1･m-2･s2
5.5 電気回路LTspiceでの解析(過渡応答) 68 電気回路シミュレータ結果 ステップ応答 トルク入力での step input 伝達トルク応答結果 （ステップ応答) τ1 （トルクは電流となる） τ2 τ3 τ4 τ3 τ6 20ms 100ms 200ms （コイル素子に対してステップ状の 電流を流すことは望ましくない）
5.5 LTspice 電気回路LTspiceでの解析(過渡応答) Transient線図 69 クランクープロペラ系の振動系 の電気回路 トルク入力での 伝達角速度/トルク応答結果 （sin波入力) （トルクは電流となる） 過渡解析 機械系のパラメータ I=3/g= 0.3059kg･m2 Id=0.2/g= 0.02039 kg･m2 c= 6×10-7×(g)rad･kg-1･m-2･s2 = 5.884×10-6 rad･kg-1･m-2･s2 cd= 8×10-6×(g)rad･kg-1･m-2･s2 = 7.845×10-5 rad･kg-1･m-2･s2
5.5 電気回路LTspiceでの解析(過渡応答) クランクープロペラ系 振動系の電気回路 80mV 伝達角速度 60.6Hz入力 𝜃ሶ 6 70 𝜃ሶ 7 トルク入力での 伝達角速度/トルク応答結果 𝜃ሶ1 𝜃ሶ 2 𝜃ሶ 3 𝜃ሶ4 （sin波入力:60.6Hz) 𝜃ሶ 5 0mV 𝜃ሶ1 𝜃ሶ 2 𝜃ሶ 3 𝜃ሶ4 𝜃ሶ 5 𝜃ሶ 6 トルク ねじり振動の共振周波数 τ1 0ms τ2 τ3 τ4 τ5 τ6 10ms 20ms p αp(°) fp(Hz) 1 29.6 60.61 2 57.5 114.12 3 (67.6) 131.98 4 90.8 168.93 5 120.3 205.78 6 150.2 229.28 𝜃ሶ 7
5.5 電気回路LTspiceでの解析(過渡応答) クランクープロペラ系 振動系の電気回路 250mV 114.1Hz ሶ 𝜃ሶ1 𝜃ሶ 2 𝜃ሶ 3 𝜃ሶ4 𝜃ሶ 5 𝜃6 71 伝達角速度 トルク入力での 𝜃ሶ 7 伝達角速度/角度応答結果 （sin波入力:114.1Hz) 0mV 𝜃ሶ1 𝜃ሶ 2 𝜃ሶ 3 𝜃ሶ4 𝜃ሶ 5 𝜃ሶ 6 回転角度 ねじり振動の共振周波数 𝜃7 𝜃1 20ms 40ms p αp(°) fp(Hz) 1 29.6 60.61 2 57.5 114.12 3 (67.6) 131.98 4 90.8 168.93 5 120.3 205.78 6 150.2 229.28 𝜃ሶ 7
6. 連続体のモビリテイ ここでは弾性梁の例について電気回路とのアナロジーに基づき 解析を行ってみます この弾性体梁のモビリティを求め、振動状態を解析します 細長い梁の左端を手にもって振動を与えてい る様子です。写真は19世紀の終わりに Marey（フランス）によってストロボ写真撮 影された世界最初のものだそうです。 LIFE, Jan.23, (1967) 72
6.1 弾性梁のモビリテイ(方程式) b 梁の動たわみ方程式 4 𝑦 2 𝜕 𝑌 𝜌𝑎 𝜕 𝑌 𝐸𝐼 4 = −𝜌𝑎 − 𝜕𝑥 𝑔 𝜕𝑡 2 剪断力の差 重量 pin支持 慣性力 はりの変位Yを静的ysと動的ydの部分に分ける ∆𝑥 𝐼, 𝐸, 𝜌, 𝑎 𝑜 ഥ 𝑗𝜔𝑡 Θ = 𝜃 + Θ𝑒 小片Δxの両端の剪断力の差は、小片に 作用する重量と慣性力に等しい （ここでは重力単位系を用いている） 梁の変位： 73 𝑥 h 断面 𝑌 𝑙 片持ち梁機構 分布質量の弾性梁の運動 𝐸: ヤング率 𝐼: 断面二次モーメント 𝜌: 単位長さの重量 𝑎: 断面積 𝑎 = 𝑏ℎ 𝐼 = 𝑏ℎ3 Τ12 𝑌 = 𝑦𝑠 + 𝑦𝑑 𝑑4 𝑦𝑠 1. 静的変位： 𝐸𝐼 = −𝜌𝑎 𝑑𝑥 4 𝑦𝑠 4 2 2. 動的変位： 𝐸𝐼 𝜕 𝑦𝑑 + 𝜌𝑎 𝜕 𝑦𝑑 = 0 𝑦𝑑 𝜕𝑥 4 𝑔 𝜕𝑡 2 重力単位系：質量の代わりに重量 （力）を基本単位として用いる
6.2 弾性梁のモビリテイ (静的変位) 74 𝑑4 𝑦𝑠 1. 静的変位： 𝐸𝐼 4 = −𝜌𝑎 𝑑𝑥 𝑥=𝑙 𝑥=0 1. 𝑦𝑠 = 0 (変位ゼロ) 2. 𝑑𝑦𝑠 =𝜃 𝑑𝑥 𝑑 2 𝑦𝑠 3. =0 𝑑𝑥 2 (曲げモーメントゼロ) 3 4. 𝑑 𝑦𝑠 =0 𝑑𝑥 3 (剪断力ゼロ) (角度θ) 2 解： 𝑦𝑠 = − 𝜌𝑎 4 𝜌𝑎𝑙 3 𝜌𝑎𝑙 2 𝑥 + 𝑥 − 𝑥 + 𝜃𝑥 24𝐸𝐼 6𝐸𝐼 4𝐸𝐼 (上記境界条件を満足する解は簡単に求まる） x=0における静的モーメント 𝑑2 𝑦𝑠 𝜌𝑎𝑙2 𝑀𝑠 = −𝐸𝐼 = 𝑑𝑥 2 𝑥=0 2 𝐸: ヤング率 𝐼: 断面二次モーメント 𝜌: 単位長さの重量 𝑎: 断面積 𝑎 = 𝑏ℎ 𝐼 = 𝑏ℎ3 Τ12 𝑔 = 980.0 𝑐𝑚 𝑠 −2 𝑙 = 400 𝑐𝑚, ℎ = 1 𝑐𝑚 𝐸 = 2 × 106 𝑘𝑔𝑓 𝑐𝑚−2 𝜌 = 7.8 × 10−3 𝑘𝑔 𝑐𝑚−3 (書籍では, h=0.05cm) 重力単位系：質量の代わりに重量（力） を基本単位として用いる
6.3 2. 動的変位： 弾性梁のモビリテイ(動的変位) 75 𝜕 4 𝑦𝑑 𝜌𝑎 𝜕 2 𝑦𝑑 𝐸𝐼 + =0 𝜕𝑥 4 𝑔 𝜕𝑡 2 弾性梁の動的変位は位置ｘと時間ｔ の偏微分方程式となる 今回、正弦波的たわみydについては 𝑦𝑑 = 𝑦ത𝑑 𝑒 𝑗𝜔𝑡 として、たわみをｘのみの関数 とすれば 𝑑4 𝑦ത𝑑 𝜌𝑎𝜔2 − 𝑦ത = 0 𝑑𝑥 4 𝑔𝐸𝐼 𝑑 𝐸: ヤング率 𝐼: 断面二次モーメント 𝜌: 単位長さの重量 𝑎: 断面積 𝑎 = 𝑏ℎ 𝐼 = 𝑏ℎ3 Τ12 𝑔 = 980.0 𝑐𝑚 𝑠 −2 𝑙 = 400 𝑐𝑚, ℎ = 1 𝑐𝑚 𝐸 = 2 × 106 𝑘𝑔𝑓 𝑐𝑚−2 𝜌 = 7.8 × 10−3 𝑘𝑔 𝑐𝑚−3 常微分方程式となる 重力単位系：質量の代わりに重量（力） を基本単位として用いる
6.3 2. 動的変位： 弾性梁のモビリテイ(動的変位) 76 𝑑 4 𝑦ത𝑑 𝜌𝑎𝜔2 − 𝑦ത = 0 𝑑𝑥 4 𝑔𝐸𝐼 𝑑 無重力状態に置かれた梁の一端に入力として ഥ 𝑒 𝑗𝜔𝑡 Θ=Θ 上記の正弦波的角度変位を与えると境界条件は 𝑥=0 1. 𝑦ത𝑑 = 0 (変位ゼロ) 𝑑𝑦ത𝑑 ഥ 2. =Θ 𝑑𝑥 (角度θ) 𝑥=𝑙 𝑑2 𝑦ത𝑑 3. =0 𝑑𝑥 2 (曲げモーメントゼロ) 𝐸: ヤング率 𝐼: 断面二次モーメント 𝜌: 単位長さの重量 𝑎: 断面積 𝑎 = 𝑏ℎ 𝐼 = 𝑏ℎ3 Τ12 𝑑3 𝑦ത𝑑 4. =0 𝑑𝑥 3 (剪断力ゼロ) 動的変位の常微分方程式の解は次の様になる。A,B,C,Dを境界条件から求めることで弾性梁の 正弦波的振動モードが求まることになる 𝜌𝑎𝜔2 4 𝜇 = 𝑦ത𝑑 = 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 + 𝐶 cos 𝜇𝑥 + 𝐷 sin 𝜇𝑥 𝑔𝐸𝐼 双曲線関数、三角関数は上記動的変位の常微分方程式の解であることは明らか。よってこれら関数の一次結合も解である
6.3 弾性梁のモビリテイ(動的変位) 77 弾性梁の正弦波的振動モード 解： 𝑦𝑑 𝑥, 𝑡 ഥ 𝑒 𝑗𝜔𝑡 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 + 𝐶 cos 𝜇𝑥 + 𝐷 sin 𝜇𝑥 =Θ 𝜌𝑎𝜔2 𝜇 = 𝑔𝐸𝐼 4 境界条件を満足するA,B,C,Dは以下のようになる 𝐴= 1 cosh 𝜇𝑙 sin 𝜇𝑙 − sinh 𝜇𝑙 𝑐𝑜𝑠 𝜇𝑙 2𝜇 (1 + cosh 𝜇𝑙 cos 𝜇𝑙) 𝐵= 1 cosh 𝜇𝑙 cos 𝜇𝑙 − sinh 𝜇𝑙 𝑠𝑖𝑛 𝜇𝑙 + 1 − 2𝑠𝑖𝑛2 𝜇𝑙 2𝜇 (1 + cosh 𝜇𝑙 cos 𝜇𝑙) C = −A ഥ Θ 𝐷 = −B 𝜇 境界条件を満たす A,B,C,D項 𝐸: ヤング率 𝐼: 断面二次モーメント 𝜌: 単位長さの重量 𝑎: 断面積 A,B,C,D項の求め方は次頁以降に示します 重力単位系：質量の代わりに重量（力） を基本単位として用いる
A,B,C,D項の求め方 １/4 𝑑4 𝑦ത𝑑 𝜌𝑎𝜔2 − 𝑦ത = 0 𝑑𝑥 4 𝑔𝐸𝐼 𝑑 1. 𝑦ത𝑑 = 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 + 𝐶 cos 𝜇𝑥 + 𝐷 sin 𝜇𝑥 これら双曲線関数、三角関数は左記の常微分方程式 の解になるのは明らかです このA,B,C,Dを次の境界条件から求める 𝑥=0 1. 𝑦ത𝑑 = 0 (変位ゼロ) 𝑑𝑦ത𝑑 ഥ 2. =Θ 𝑑𝑥 (角度θ) 2. 𝜌𝑎𝜔2 𝜇 = 𝑔𝐸𝐼 4 𝑥=𝑙 𝑑2 𝑦ത𝑑 3. =0 𝑑𝑥 2 (曲げモーメントゼロ) 𝑑3 𝑦ത𝑑 4. =0 𝑑𝑥 3 ഥ 𝑗𝜔𝑡 Θ = 𝜃 + Θ𝑒 (剪断力ゼロ) 𝑑𝑦ത𝑑 = 𝜇 𝐴 sinh 𝜇𝑥 + 𝐵 cosh 𝜇𝑥 − 𝐶 sin 𝜇𝑥 + 𝐷 cos 𝜇𝑥 𝑑𝑥 𝑑2 𝑦ത𝑑 3. = 𝜇 2 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 − 𝐶 cos 𝜇𝑥 − 𝐷 sin 𝜇𝑥 2 𝑑𝑥 𝑑3 𝑦ത𝑑 4. = 𝜇 3 𝐴 sinh 𝜇𝑥 + 𝐵 cosh 𝜇𝑥 + 𝐶 sin 𝜇𝑥 − 𝐷 cos 𝜇𝑥 3 𝑑𝑥 78
A,B,C,D項の求め方 ２/4 A,B,C,D項の求め方 𝑥=0 1. 𝑦ത𝑑 = 0 1. 𝑦ത𝑑 = 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 + 𝐶 cos 𝜇𝑥 + 𝐷 sin 𝜇𝑥 = A + C = 0 ∴ C = −A （変位ゼロ） 𝑑𝑦ത𝑑 = 𝜇 𝐴 sinh 𝜇𝑥 + 𝐵 cosh 𝜇𝑥 − 𝐶 sin 𝜇𝑥 + 𝐷 cos 𝜇𝑥 𝑑𝑥 𝑑𝑦ത𝑑 ഥ 2. =Θ 𝑑𝑥 2. （角度θ） 𝑑𝑦ത𝑑 ഥ =𝜇 𝐵+𝐷 =Θ 𝑑𝑥 ഥ Θ ∴𝐷 = −B 𝜇 ഥ 𝑗𝜔𝑡 Θ = 𝜃 + Θ𝑒 79
A,B,C,D項の求め方 ３/4 𝑥=𝑙 𝑑 2 𝑦ത𝑑 3. =0 𝑑𝑥 2 (曲げモーメントゼロ) 𝑑2 𝑦ത𝑑 3. = 𝜇 2 𝐴 cosh 𝜇𝑥 + 𝐵 sinh 𝜇𝑥 − 𝐶 cos 𝜇𝑥 − 𝐷 sin 𝜇𝑥 2 𝑑𝑥 𝑑2 𝑦ത𝑑 ഥ Τ𝜇 − B sin 𝜇𝑙 = 0 = 𝜇 2 𝐴 cosh 𝜇𝑙 + 𝐵 sinh 𝜇𝑙 + 𝐴 cos 𝜇𝑙 − Θ 2 𝑑𝑥 ഥ Τ𝜇 − B sin 𝜇𝑙 = 0 𝐴 cosh 𝜇𝑙 + 𝐵 sinh 𝜇𝑙 + 𝐴 cos 𝜇𝑙 − Θ ഥ Τ𝜇 sin 𝜇𝑙 𝐴 cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 + 𝐵 sinh 𝜇𝑙 + sin 𝜇𝑙 = Θ 𝑑 3 𝑦ത𝑑 4. =0 𝑑𝑥 3 𝑑3 𝑦ത𝑑 4. = 𝜇 3 𝐴 sinh 𝜇𝑥 + 𝐵 cosh 𝜇𝑥 + 𝐶 sin 𝜇𝑥 − 𝐷 cos 𝜇𝑥 3 𝑑𝑥 (剪断力ゼロ) 𝑑3 𝑦ത𝑑 = 𝜇 3 𝐴 sinh 𝜇𝑙 + 𝐵 cosh 𝜇𝑙 + 𝐶 sin 𝜇𝑙 − 𝐷 cos 𝜇𝑙 = 0 3 𝑑𝑥 ഥ Τ𝜇 − B cos 𝜇𝑙 = 0 𝐴 sinh 𝜇𝑙 + 𝐵 cosh 𝜇𝑙 − 𝐴 sin 𝜇𝑙 − Θ ഥ Τ𝜇 cos 𝜇𝑙 𝐴 sinh 𝜇𝑙 − sin 𝜇𝑙 + 𝐵 cosh 𝜇𝑙 + cos 𝜇𝑙 = Θ これらから前頁の結果を用いC,D項を消去すると、A,B項のみの連立方程式となる 80
A,B,C,D項の求め方 ４/4 A,B項は次の連立方程式となる cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 sinh 𝜇𝑙 − sin 𝜇𝑙 ഥ Τ𝜇 sin 𝜇𝑙 sinh 𝜇𝑙 + sin 𝜇𝑙 𝐴 Θ = ഥ Τ𝜇 cos 𝜇𝑙 cosh 𝜇𝑙 + cos 𝜇𝑙 𝐵 Θ cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 𝐴 = sinh 𝜇𝑙 − sin 𝜇𝑙 𝐵 ഥ Τ𝜇 sin 𝜇𝑙 sinh 𝜇𝑙 + sin 𝜇𝑙 −1 Θ ഥ Τ𝜇 cos 𝜇𝑙 cosh 𝜇𝑙 + cos 𝜇𝑙 Θ 1 cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 𝐴 = 𝐵 𝑊 −sinh 𝜇𝑙 + sin 𝜇𝑙 − sinh 𝜇𝑙 − sin 𝜇𝑙 cosh 𝜇𝑙 + cos 𝜇𝑙 ഥ Τ𝜇 sin 𝜇𝑙 Θ ഥ Τ𝜇 cos 𝜇𝑙 Θ A,B項の連立方程式 前頁までの結果からC,D項を 消去できるから、A,B項のみ の連立方程式となる 逆行列をかける 逆行列計算 𝑊 = cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 cosh 𝜇𝑙 + cos 𝜇𝑙 − sinh 𝜇𝑙 + sin 𝜇𝑙 sinh 𝜇𝑙 − sin 𝜇𝑙 = 𝑐𝑜𝑠ℎ2 𝜇𝑙 − 𝑠𝑖𝑛ℎ2 𝜇𝑙 + 𝑐𝑜𝑠 2 𝜇𝑙 + 𝑠𝑖𝑛2 𝜇𝑙 + 2 cosh 𝜇𝑙 cos 𝜇𝑙 = 1 + 1 + 2 cosh 𝜇𝑙 cos 𝜇𝑙 𝑊 = 2 (1 + cosh 𝜇𝑙 cos 𝜇𝑙) A,B項は次のようになる ഥ Τ𝜇 sin 𝜇𝑙 − sinh 𝜇𝑙 + sin 𝜇𝑙 Θ ഥ Τ𝜇 cos 𝜇𝑙 𝐴 = cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 Θ 式を整理 ഥ Τ𝜇𝑊 𝐴 = cosh 𝜇𝑙 sin 𝜇𝑙 − sinh 𝜇𝑙 𝑐𝑜𝑠 𝜇𝑙 Θ ഥ Τ𝜇 sin 𝜇𝑙 + cosh 𝜇𝑙 + 𝑐𝑜𝑠 𝜇𝑙 Θ ഥ Τ𝜇 cos 𝜇𝑙 𝐵 = − sinh 𝜇𝑙 + sin 𝜇𝑙 Θ 式を整理 ഥ Τ𝜇𝑊 𝐵 = cosh 𝜇𝑙 cos 𝜇𝑙 − sinh 𝜇𝑙 𝑠𝑖𝑛 𝜇𝑙 + 1 − 2𝑠𝑖𝑛2 𝜇𝑙 Θ 81
6.3 弾性梁のモビリテイ(変位計算例) 82 弾性梁の 変位計算例 静的＋動的変位 入力周波数:2Hz 条件パラメータ 𝑔 = 980.0 𝑐𝑚 𝑠 −2 𝑙 = 400 𝑐𝑚, ℎ = 1 𝑐𝑚 𝐸 = 2 × 106 𝑘𝑔𝑓 𝑐𝑚−2 𝜌 = 7.8 × 10−3 𝑘𝑔 𝑐𝑚−3 𝜃 =0.1 rad, Θ =0.2 rad 重力単位系 動的変位 1周期を10 (tnum)に 分割して表示している tnum=10 https://github.com/m4881shimojo/TTanaka_Eng/tree/main/beamMotion beamv2.py
弾性梁の 動的変位（励起振動数を変化させた） 0.2Hz 2Hz 図の上段は梁全体の変位を 下段は振動成分のみを表示 4Hz 梁全体の変位 振動成分 7Hz 10Hz 12Hz 梁全体の変位 振動成分 tnum=6 https://github.com/m4881shimojo/TTanaka_Eng/tree/main/beamMotion beamv2.py 83
6.4 弾性梁のモビリテイ(共振) 84 共振周波数： ഥ 𝑗𝜔𝑡 とすると O点での反作用モーメントは 𝑀 = 𝑀𝑒 𝑑2 𝑦ത𝑑 sinh 𝜇𝑙 cos 𝜇𝑙 − cosh 𝜇𝑙 sin 𝜇𝑙 ഥ ഥ = −𝐸𝐼 𝑀 = 𝜇𝐸𝐼 Θ 𝑑𝑥 2 𝑥=0 1 + cosh 𝜇𝑙 cos 𝜇𝑙 よって（回転の）モビリティーは 𝜆 𝑗𝜔 = 𝑖 𝜔 𝜆 𝑗𝜔 = 4 1 + cosh 𝜇𝑙 cos 𝜇𝑙 sinh 𝜇𝑙 cos 𝜇𝑙 − cosh 𝜇𝑙 sin 𝜇𝑙 𝜌 𝑎𝐸 3 𝐼 3 𝑔 ഥ 𝑗𝜔𝑡 Θ = 𝜃 + Θ𝑒 𝜕𝑦ത𝑑 Τ𝜕𝑥 ]𝑥=0 𝑀 𝜌𝑎𝜔2 𝜇 = 𝑔𝐸𝐼 4 となる 回転モビリティの極致を考えてみると、|λ(jwt)|＝０、|λ(jwt)|=∞の場合である 1. 𝜆 𝑗𝜔 =0 ：わずかな入力角変位で梁が猛烈に振れる 共振 2. 𝜆 𝑗𝜔 = ∞ ：梁を動かすのに手ごたえがない 𝜆 𝑗𝜔 = 0 ⟼ cosh 𝜇𝑙 cos 𝜇𝑙 + 1 =0 𝜆 𝑗𝜔 = ∞ ⟼ tanh 𝜇𝑙 − tan 𝜇𝑙 = 0 共振周波数が求められる （回転モビリティー分母がゼロ）
6.4 弾性梁のモビリテイ(共振) 85 共振周波数： 共振周波数 λ(jω)=0 cosh 𝜇𝑙 cos 𝜇𝑙 = −1 𝜇𝑙 = 𝑙 4 𝜌𝑎 𝜔 𝑔𝐸𝐼 λ(jω)=∞ tanh 𝜇𝑙 = tan 𝜇𝑙 μl frequency μl frequency μ1l=1.8751 0.5061 μ’1l=3.9266 2.2193 μ2l=4.6941 3.1717 μ’2l=7.0686 7.1921 μ3l=7.8548 8.8810 μ’3l=10.2102 15.0058 μ4l=10.9955 17.4028 μ’4l=13.3518 25.6608 μ5l=14.1372 28.7685 μ’5l=16.4934 39.1571 μ6l=17.2788 42.9752 μ’6l=19.6350 55.4948 1. 𝜆 𝑗𝜔 =0 ：わずかな入力角変位で梁が猛烈に振れる 2. 𝜆 𝑗𝜔 = ∞ ：梁を動かすのに手ごたえがない 共振周波数はμlから各種パラメータの値を用いて計算する 上記の値の場合：l=400cm,h=1cm,E= 2000000.0kgf/cm-2,ro= 0.0078kgcm-3,g=980.0cms-2,a=1cm2
非線形方程式の解を求める 状態 86 状態となる条件 𝜆 𝑗𝜔 = 0 cosh 𝜇𝑙 cos 𝜇𝑙 + 1 =0 𝜆 𝑗𝜔 = ∞ tanh 𝜇𝑙 − tan 𝜇𝑙 = 0 1. 関数をplotする 𝜇𝑙 = 𝑙 4 𝜌𝑎 𝜔 𝑔𝐸𝐼 次の方法で見つける ➢ plotの結果からゼロとなる区間[a,b]を求める 2. はさみうち法を用いて解を求める ➢ 区間[a,b]を指定してはさみうち法（数値解析）により解を求める https://github.com/m4881shimojo/TTanaka_Eng/tree/main/function_ans （programです） はさみうち法 関数をplotする ふろくに説明あります
弾性梁の 動的変位（共振例） 10000cm 3.17Hz 図の上段は梁全体の変位を 下段は振動成分のみを表示 20000cm 共振状態 共振状態 8.88Hz μ2l=4.6941 μ3l=7.8548 前頁で示した「 2Hz,4Hz,7Hz,10Hz 」と比べてみると変位が異常 https://github.com/m4881shimojo/TTanaka_Eng/ beamv2.py 87
弾性梁のモビリテイ 88 写真は19世紀の終わりに Marey（フランス）によって ストロボ写真撮影された世界 最初のものだそうです。 LIFE, Jan.23, (1967) 4Hz 縦軸の単位 長さは？ 1. この解説では、左図にある式は使用してい ません 2. また、左図でパラメータとしてh=0.05cm とありますが、これも採用していません。 これだと変形量が大きくなりすぎてしまい ます。私はh=1cmを使用しました 3. 私が用いた計算プログラムを下記GitHubに 置きます。計算式、パラメータなど正しい かは疑問ですが、参考になれば幸いです 0.5Hz p.339 中田孝:工学解析(技術者のための数学手法),オーム社,1972.
6.5 弾性梁のモビリテイ(疑似集中定数系) 89 弾性梁の回転モビリティ 𝑖 𝜔 𝜆 𝑗𝜔 = 4 1 + cosh 𝜇𝑙 cos 𝜇𝑙 sinh 𝜇𝑙 cos 𝜇𝑙 − cosh 𝜇𝑙 sin 𝜇𝑙 𝜌 𝑎𝐸 3 𝐼 3 𝑔 2 𝜌𝑎𝜔 𝜇4 = 𝑔𝐸𝐼 cosh(x) cos (x), sinh(x) sin (x)を級数展開 再構成+まとめる 弾性梁のモビリテイ 第一次固有振動数以下 ならば近似可能とした 0.4772 1 𝑙 𝜌 Τ𝑔 𝑙 3 3 𝜆 𝑗𝑓 = 𝑗 − + 1.6614 𝑓 − 60 2 5 𝑓 𝜌Τ𝑔 𝑏ℎ𝑙 3 𝑓 𝐸𝐼 𝐸 𝑏ℎ 1 1 1 = + 𝜆 𝜆1 𝜆2 + 𝜆3 梁の慣性モーメント＆コンプライアンスｃを求める 上記近似式に代入＆式を整理 𝜆 𝑗𝑓 = 𝑗 − 1 1 𝑙 + 0.528𝜋 𝑓 2𝜋𝐽 𝑓 𝐸𝐼 弾性梁の疑似集中定数系近似 𝜆 Θ 𝜆1 = 𝑗𝜔𝐽1 −1 𝜆2 = 𝑗𝜔𝐽2 −1 𝐽1 𝑀 𝐽2 𝜆3 = 𝑗𝜔𝑐 疑似集中定数系
6.5 弾性梁のモビリテイ(疑似集中定数系) 90 弾性梁の疑似集中定数系近似の導出についてこれから示します すでに導出した回転 モビリティ 𝑖 𝜔 𝜆 𝑗𝜔 = 4 1 + cosh 𝜇𝑙 cos 𝜇𝑙 sinh 𝜇𝑙 cos 𝜇𝑙 − cosh 𝜇𝑙 sin 𝜇𝑙 𝜌 3 𝐼3 𝑎𝐸 𝑔 𝜌𝑎𝜔2 𝜇 = 𝑔𝐸𝐼 4 cosh(x) cos (x)などを級数展開再構成+整理 3 2 3! 22 1− + 近似式： 𝜆 𝑗𝜔 = −𝑗𝜔 𝐸𝐼𝜇 4 𝑙3 4! 7! 𝜔 = 2𝜋𝑓, 𝐼 = 𝑏ℎ3 , 𝑎 = 𝑏ℎ の関係を代入すると 3 3 4 2 3! 2 3! 2 𝜇𝑙 4 + + + 8! 4! 7! 11! 0.4772 1 𝑙 𝜌 Τ𝑔 𝑙 3 3 𝜆 𝑗𝑓 = 𝑗 − + 1.6614 𝑓 − 60 2 5 𝑓 𝜌Τ𝑔 𝑏ℎ𝑙3 𝑓 𝐸𝐼 𝐸 𝑏ℎ 𝜇𝑙 8 ⋯ (10.139) 近似式の第一項は静止項 ここで梁の長さ l の慣性モーメントＪを求めてみる 𝑙 この梁の回転 モビリティ 𝜆 = 𝑠𝐽 −1 = 𝑗𝜔 𝐽 −1 → 𝜆 = 𝑗2𝜋𝑓 𝐽 −1 剛体梁と仮定し たモビリティ 3 1 0.4772 1 𝜆𝑟 = −𝑗 ∙ = −𝑗 ∙ 𝜌Τ𝑔 𝑏ℎ𝑙3 2𝜋𝑓 𝜌Τ𝑔 𝑏ℎ𝑙3 𝑓 1 1 𝐽 = 𝑀𝑙 2 = 𝜌Τ𝑔 𝑏ℎ𝑙 3 3 3
6.5 弾性梁のモビリテイ(疑似集中定数系) 91 近似式： 0.4772 1 𝑙 𝜌 Τ𝑔 𝑙 3 3 𝜆 𝑗𝑓 = 𝑗 − + 1.6614 𝑓 − 60 2 5 𝑓 𝜌Τ𝑔 𝑏ℎ𝑙3 𝑓 𝐸𝐼 𝐸 𝑏ℎ 剛体梁の慣性 モーメント 振動３次項 振動１次項 質量分布梁機構の疑似集中定数回路へ 入力端に未知の慣性モーメントJ1が集中し、次に質量を無視した 梁全体コンプライアンスに相当するｃのねじりバネが続き、その 先端に未知な慣性モーメントJ2が結合されたものとと考えられる 𝐽2 𝐽1 + 𝐽2 2 𝑐2𝜋𝑓 (10.141) 𝜆2 = 𝑗𝜔𝐽2 −1 𝐽1 𝑀 1 1 1 1 = + = 𝑗𝜔𝐽1 + 1 𝜆 𝜆1 𝜆2 + 𝜆3 + 𝑗𝜔𝑐 𝑗𝜔𝐽2 1 𝜆=𝑗 − + 2𝜋 𝐽1 + 𝐽2 𝑓 𝜆 Θ 𝜆1 = 𝑗𝜔𝐽1 −1 𝐽2 𝜆3 = 𝑗𝜔𝑐 𝐽1 = 0.274𝐽 𝐽2 = 0.726𝐽 1 𝐽 = 𝜌Τ𝑔 𝑏ℎ𝑙3 3 𝑐 = 0.5 𝑙Τ𝐸𝐼
6.5 弾性梁のモビリテイ(疑似集中定数系) 梁全体コンプライ アンスｃを求め、 (10.141) に代入 𝜆=𝑗 − 𝜃 𝑃𝑙2 𝑙 𝑐= = ൘ 𝑃𝑙 = 𝑀 2𝐸𝐼 2𝐸𝐼 1 𝐽2 + 2𝜋 𝐽1 + 𝐽2 𝑓 𝐽1 + 𝐽2 2 𝑐= 𝑙 2𝐸𝐼 𝑀 𝜋𝑙 𝑓 𝐸𝐼 𝐸𝐼 𝑃 𝜃 (10.143) 92 𝑙 梁の曲げ Τ𝑔 𝑙 3 3 0.4772 1 𝑙 𝜌 近似式：𝜆 𝑗𝑓 = 𝑗 − + 1.6614 𝑓 − 60 2 5 𝑓 𝜌Τ𝑔 𝑏ℎ𝑙3 𝑓 𝐸𝐼 𝐸 𝑏ℎ 3 𝜌Τ𝑔 𝑏ℎ𝑙 1 𝐽1 + 𝐽2 = ∙ 0.4772 2𝜋 1 𝐽1 + 𝐽2 = 𝜌Τ𝑔 𝑏ℎ𝑙3 = 𝐽 3 疑似集中定数系： 𝐽2 𝐽1 + 𝐽2 𝐽2 𝐽1 + 𝐽2 2 考慮外 𝜋𝑙 𝑙 𝑓 = 1.6614 𝑓 𝐸𝐼 𝐸𝐼 2 = 1 1 𝑙 𝜆 𝑗𝑓 = 𝑗 − + 0.528𝜋 𝑓 2𝜋𝐽 𝑓 𝐸𝐼 1.6614 = 0.7262 𝜋 第一次固有振動数以下ならば 近似可能とした
6.6 連続体を含む回路の計算（例） 図は航空発動機とプロペラの結合を示したものでピス トン、クランクは一括してλ0の集中慣性モビリティ、 またプロペラ軸はλ1の集中コンプライアンスモビリ ティ、プロラハブは集中慣性モビリティλ2とし、その 先に分布質量梁と仮定したプロペラの連続体モビリ ティーλpが接続されたものとします 𝜆0 = 𝑗𝜔𝐽0 −1 𝜆𝑝 𝐽0 𝜆2 = 𝑗𝜔𝐽2 −1 𝜆1 = 𝑗𝜔𝑐1 𝜆𝑝 𝜆𝑝 𝜆1 𝜆0 𝜆2𝑝 𝑐1 左端入力端からのモビリティは次の様になります ここに数式を入力します。 1 1 1 = + 𝜆 𝜆0 𝜆1 + 𝜆2𝑝 1 1 = + 1 𝜆0 𝜆 + 1 1 1 1 + + 𝜆2 𝜆𝑝 𝜆𝑝 93 𝜆𝑝 並列 𝜆2 航空発動機とプロペラの振動回路 並列回路のモビリティとなるため、その逆数の和を計算する。まず系は発動機系のλ0 と、プロペラ軸λ1&プロペラλpが並列回路となる。そしてまた後者はプロペラハブλ2と プロペラλpの並列となる。
6.6 連続体を含む回路の計算（例） 94 よってモビリティは次の様になる 1 1 𝜆 𝜔 = = 1 1 1 + 𝑗𝜔𝐽 + 0 1 1 𝜆0 𝜆1 + 1 𝑗𝜔𝑐 + 1 1 1 2 + + 𝑗𝜔𝐽1 + 𝜆2 𝜆𝑝 𝜆𝑝 𝜆𝑝 𝜔 ここでプロペラのモビリティーが弾性梁で代用できるとすると −𝑗 𝜆 𝑗𝜔 = 𝜔𝐽0 − 1 𝜔𝑐1 − 𝜔 𝑓 𝜔 = 4 1 𝜔𝐽1 − 2 𝑓(𝜔) 1 + cosh 𝜇𝑙 cos 𝜇𝑙 sinh 𝜇𝑙 cos 𝜇𝑙 − cosh 𝜇𝑙 sin 𝜇𝑙 𝜌 3 𝐼3 𝑎𝐸 𝑔 𝜇4 = (10.149) 𝜌𝑎𝜔2 𝑔𝐸𝐼 (10.149)の分母は実数であるから角周波数ωを変数として計算できます。ここで |λ(jω)|=∞となるωよりプロペラ発動機系の数次にわたる共振周波数が求められます。 pythonで計算させるのは簡単と思います。今回具体的数値が不明なためパスしました。興味がある方は、 すでに似たような数値計算・PLOTを行っていますので参考にしてください。
つづく ⚫ やっとpart2です。今回は弾性梁解析に時間がかかりました。計算結果が実際と大き く違ったためです。これを書籍での板厚を0.5mmから10mmに変更してみました。計算 結果はそれらしい値になりましたが、しかし、どうなんでしょうか。 ⚫ 書籍「工学解析」は重力単位系を使っています。当時、機械系では重力単位系が主で した。また学部の材料力学では英語の教科書でしたので単位がヤードポンド法 (1foot=12inch、1pound=16ounceなど長さや重さの単位)で混乱した思い出があります。今回も単位 系の問題で私がミスをしているのかもしれません。 ⚫ LTspiceのプログラムはスライドにある (.asc)そのものです。簡単に作成できますの で興味のある方は試してください。(LTspiceを使える前提） ⚫ ただし、今回はpythonを使った計算がいくつかあります。このプログラムに関しては GitHubに置いておきます。パラメータなどを変化させて挙動の変化を確認することも 可能です。 https://github.com/m4881shimojo/TTanaka_Eng/ ⚫ なおサイズモ系や動吸振器については工学解析の書籍では扱っていません。よって、 もし内容に間違がある場合は私の間違いです。 ⚫ 最近体調などの問題で根を詰めての仕事が難しくなりました。次回は少し時間がかか るかもしれません。 ⚫ アナロジーは形式的なもので、物理的意味はない、という意見もあります。どうなの でしょうか、私も知りたいところです 95
参考文献 96 中田孝:工学解析(技術者のための数学手法),オーム社,1972、総頁572.
補足：中田孝氏について 中田 孝（1908- 2000）：日本の機械工学者。東京工業大学名誉教授。 元日本機械学会会長。元計測自動制御学会会長。元日本学士院会員。日本 学士院賞受賞。 https://ja.wikipedia.org/wiki/中田孝 中田孝先生は、大学院時代に私が所属した研究室の前教授でした。著書「転位歯車」は 修論では大変お世話になりました。しかし、いろいろな逸話などは聞くことはありまし たが、中田先生にはお会いしたことは残念ながらありませんでした。先日、同期で研究 室を引き継いだ北條さん（東工大名誉教授）に会った時、その話をしたところ、中田先 生の講義ビデオを送って頂き、初めてその肉声に触れることができました。感謝します。 そのビデオを見ながら、中田先生は、歯車の研究で日本学士院賞を45歳で受賞したのを はじめ、自動制御、NC工作機などの分野も手掛けた非常にエンジニアリングセンスの高 い方であり、また数学や電気電子工学に対する造詣も大変深い方だと認識を新たにしま した。またビデオの中で“歯型理論にはあまり手を付けたくなかったが電気回路理論の 応用に興味を持っていたので研究を進めた”とのあたりは、非常に親近感を感じました。 中田孝先生の紹介記事です【日本の工作機械を築いた人々】 大河出版「応用機械工学」1988年4~5月号掲載、 https://www.sme-japan.org/library5.pdf 中田先生のNC工作機械を開発したのときの開発談が語られています。大変面白い内容ですのでお勧めです。 97
補足（私の卒論修論） 98 【卒研】私は卒業研究では電通大の梶谷誠先生の研究室で歯車を加工するホブ盤の制御回路を作成しました。但し、装置は東 工大石川研にありました。開発する制御回路は、ホブ盤主軸の回転角度を計測して、切削工具を駆動する電気油圧パルスモー タを制御するものです。この電気油圧パルスモータは、電気パルスモータの出力をプランジャー形油圧モータでトルク増幅す るという仕組みでした。初めてホブ盤を見たときの印象は巨大な鉄の塊で、制御ミスでパワフルな油圧モータなどを暴走させ たらどうなるかと、恐れを抱いたものです。 試作した回路は、単純な論理素子（TTL)を組合わせた回路です。当時は論理素子の黎明期で、詳しい書籍など無い時に初 歩の初歩から学び始め、TI社の分厚いハンドブックで論理素子の細かな仕様をチェックしながら、連日連夜、回路を試作して はオシロスコープなどを使い動作を確認して、回路を設計していました。本当に充実した楽しい時間でした。ご指導頂いた梶 谷先生には心から感謝しています。 しかし、実機実験では散々で、必死で改良を重ね何とか動作させました。やっと加工できた歯車を見たときの感激は今も忘 れません。しかし、歯すじ誤差に問題がありました。その原因は切削加工力の増加に伴うモータ回転遅れです。これは電気油 圧パルスモータの最大の長所でもある開ループ制御の限界だったのでしょう。(後に中田先生が電気油圧パルスモータの開発 にも深く関与していたと知りました。奇妙な縁を感じます） 【就職】卒業後はX社に就職しました。理由はX社の子会社であるF社でのロボット開発に関心があったためです。当時F社へ はX社に就職する必要がありました。そして、入社時実習は希望先のF社でした。ところが、本配属では全く違うX社部署だっ たため、すぐに退職して東工大院の石川研(中田研の後継研究室)へ進学することにしました。人生にifはつきものですが、F 社に入っていたら全く違う人生を送ったようにも思います。 【修論】修論では習得した論理回路の設計技術を用いて、自動歯車精度測定機の開発に係り、歯形計算回路と測定機の制御回 路の開発を行いました。また修士2年の中頃、当時の先進機器であるミニコン(HP-21MX)が納入されたので、データ処理 （Fourier解析など）に利用することにしました。それからは、ミニコンと歯車測定機とのinterface回路設計から始めて、各 種プログラム作りと正月もない日々でした。またなんと、ミニコンは紙テープベースの開発システムで、コンパイラもライブ ラリもリンカも、全て毎回毎回、そのsoftwareを紙テープから読込ませる必要があり、コンパイル結果なども紙テープで TTY(110 baud rates)出力という、本当に気が狂うようなシステムでした(これと比べるとDOSは凄い)。実は、修論にはミ ニコンを利用する必要は全くなかったのですが、好奇心から無理にお願いして利用したのが真実です。これらの成果は機械学 会誌に論文として掲載できました。そして修了後は、国の研究機関に就職することに致しました。 大学時代に、昼夜を忘れて一心不乱に学んだ論理回路設計、計算機の経験は私の一つの財産になったように思います。 (DOS: Disk Operating System)
ふろく 非線形方程式の解を求めます f(x)=0 での“ｘ”を数値解法で求めます https://github.com/m4881shimojo/TTanaka_Eng/tree/main/function_ans はさみうち法 関数をplotする 99
非線形方程式の解を求める 100 次の2stepで見つける 1. 関数をplotする ➢ plotの結果からゼロとなる区間[a,b]を求める 2. はさみうち法を用いて解を求める ➢ 区間[a,b]を指定してはさみうち法（数値解析）により解を求める 例） 非線形関数f(x)： cosh 𝜇𝑙 cos 𝜇𝑙 + 1 = 0 f(x) PLOT a 0 解 1. 関数をplotする x c b 次を繰返す f(c)<0 f(c)>0 → [c,b]間に解がある → [a,c]間に解がある 2. はさみうち法を用いて解を求める
非線形方程式の解を求める 1. 関数をplotする 101 𝜆 𝑗𝜔 = 0 ⟼ cosh 𝜇𝑙 cos 𝜇𝑙 = −1 r_val=np.cosh(x)*np.cos(x)+1 関数のPLOT波形 https://github.com/m4881shimojo/TTanaka_Eng/tree/main/function_ans fnc_plot_inf_zero.py
非線形方程式の解を求める 𝜆 𝑗𝜔 = 0 1. 関数をplotする 𝜆 𝑗𝜔 = 0 ⟼ cosh 𝜇𝑙 cos 𝜇𝑙 = −1 f(x)=np.cosh(x)*np.cos(x)+1 関数がゼロになる区間[ a,b ] ① [1.5，2.5] ② [ 4,5 ] ① ③ [ 7,8 ] ④ [ 10.5,11.5 ] f(x)=0 ⑤ [ 13.5,14.5 ] ⑥ [ 16.5,17.5 ] ⑦ [ 20,21 ] ⑧ ・・・ https://github.com/m4881shimojo/TTanaka_Eng/tree/main/function_ans 102 fnc_plot_zero.py
非線形方程式の解を求める 𝜆 𝑗𝜔 = 0 2. はさみうち法を用いて解を求める 103 𝜆 𝑗𝜔 = 0 ⟼ cosh 𝜇𝑙 cos 𝜇𝑙 = −1 区間[a,b]を指定してはさみうち法 （数値解析）により解を求める f(x) a x c 0 b 新たな区間とする ① f(a)とf(b)の値が異符号となる区間 [a,b] を見つける ② 線分とｘ軸の交点を見つける ③ f(c)とf(a)が同符号の場合、[c,b]を 新たな区間[a,b]とし、区間を狭め る ④ 異なる符号の場合、[a,c]を新たな 区間[a,b]とし、区間を狭める [a,b] ⑤ f(c)が望む精度になるまで繰返す https://github.com/m4881shimojo/TTanaka_Eng/tree/main/function_ans fnc_hasamiuchi_zero.py
104