火力発電所の現場で役立てた数値解析手法１．火力発電所での技術検討の実例電気が届くまでの概略の経路送電線発電所タービン発電機変圧器変電所配電線一般家庭工場変圧器 1 送電線の地絡故障通常時ＩＥＸＧＸＴｒＩＸＬｉｎｅＸＬｏａｄＲＬｏａｄＸＬｉｎｅＸＬｏａｄＲＬｏａｄ事故時Ｉ’ ＥＸＧＸＴｒＩ’ 2 高速度再閉路時間発電機負荷雷による地絡事故の発生地絡保護装置による検出地絡事故の両側にある遮断器の自動開放事故電流の消滅発電機負荷自動開放自動開放消滅地絡事故の両側にある遮断器の自動投入発電機負荷地絡事故からの復旧自動投入自動投入 3 高速度再閉路が機器に与える影響：タービンが発電機を回すトルク：発電機が送り出す電流により押し返されるトルク（反作用）タービン発電機地絡事故が発生し、高速度再閉路をすると・・・電流により押し返されるトルクが大きく変動し、回転軸に大きな負担がかかる。この現象の解析は電気と機械の複合問題 4 質点の運動の基本式質点Ｆ［Ｎ］力剛体の回転慣性モーメントＩ［ｋｇ･ｍ２］質量ｍ［ｋｇ］速度ｖ［ｍ／ｓ］角速度加速度ａ［ｍ／ｓ２］角加速度運動方程式仕事率 F  ma Ｐ［Ｗ］＝ＦｖＴ［Ｎ・ｍ］トルク ω［ｒａｄ／ｓ］ d ［ｒａｄ／ｓ２］ dt 運動方程式仕事率 d T I dt Ｐ［Ｗ］＝Ｔω 5 剛体の回転運動の基本式質点Ｆ［Ｎ］力剛体の回転慣性モーメントＩ［ｋｇ･ｍ２］質量ｍ［ｋｇ］速度ｖ［ｍ／ｓ］角速度加速度ａ［ｍ／ｓ２］角加速度運動方程式仕事率 F  ma Ｐ［Ｗ］＝ＦｖＴ［Ｎ・ｍ］トルク ω［ｒａｄ／ｓ］ d ［ｒａｄ／ｓ２］ dt 運動方程式仕事率 d T I dt Ｐ［Ｗ］＝Ｔω ここにトルクが出てくるが、これがそのまま軸にかかるトルクというわけではない 6 タービンと発電機をつなぐ軸にかかるトルク：タービンが発電機を回すトルク：発電機が送り出す電流により押し返されるトルク（反作用）タービン発電機発電機の回転運動の式とタービンの回転運動の式を個別に立てる。軸のねじれにより伝達されるトルクは作用・反作用の法則で逆向き・同じ大きさ 7 発電機およびタービンにはたらくトルクﾀｰﾋﾞﾝ内部を流れる蒸気によるトルクタービンローター Tt 運転時の回転方向（正方向）発電機を流れる電流によるトルク発電機ローター  Ta  Tg Ta 軸のねじれにより伝達されるトルク（大きさ同じで逆向き） 8 発電機の運動方程式回転運動の基本式運転時の回転方向（正方向）発電機を流れる電流によるトルク発電機ローター  Tg Ta d T I dt Ta  Tg  I g d g dt Pg d g Ta   Ig g dt 軸のねじれにより伝達されるトルク 9 タービンの運動方程式回転運動の基本式ﾀｰﾋﾞﾝ内部を流れる蒸気によるトルク d T I dt dt Tt  Ta  I t dt dt  Ta  I t t dt Pt タービンローター Tt 運転時の回転方向（正方向）  Ta 軸のねじれにより伝達されるトルク 10 ねじれをどう考えるかタービンローター発電機ローター仮想の目印線目仮想の目印線仮想の目印線 g g  d g t dt 一定回転速度なら d t t  dt g  t  2f 11 軸のねじれによるトルクはタービンローター Tt 運転時の回転方向（正方向）発電機ローター  Ta  Tg Ta 仮想の目印線 g 軸のねじれによるトルク伝達は仮想の目印線 t Ta  k (t  g ) 12 ここまでの方程式を列挙発電機の回転の運動方程式タービンの回転の運動方程式発電機の角速度は回転角の微分 Ta  Pg g  Ig d g dt dt  Ta  I t t dt d g g  dt Pt タービンの角速度は回転角の微分 d t t  dt 軸のねじれによるトルク伝達 Ta  k (t  g ) 変数７つに式５つ・・・短時間だから Pt 一定と考えれば残るは１つ 13 発電機の出力は EV Pg  sin  Z E 発電機の内部誘起電圧 V 発電機の端子電圧 Z 発電機の同期インピーダンス  発電機の「(内部)相差角」  ってなんだろう？ 14 観測者が定格回転数で回っていると考えるタービンローターそして、出力ゼロの時の仮想の目印線を基準にして回転角を表現する。発電機ローター仮想の目印線目  g から定格回転分を差し引いた分が相差角だと考えることができる  t についても同様の変換をして    g  2ft    t  2ft 15 これですべての方程式がそろった発電機の回転の運動方程式タービンの回転の運動方程式発電機の角速度は回転角の微分 Ta  Pg g  Ig d g dt dt  Ta  I t t dt d g g  dt Pt タービンの角速度は回転角の微分 d t t  dt 軸のねじれによるトルク伝達 Ta  k (t  g ) 発電機の出力式変数変換の２式 EV Pg  sin  Z    g  2ft Pt は一定値    t  2ft 16 連立微分方程式を整理すると EV sin  d 2 Z I g 2  k (   )   d dt + 2f dt Pt d 2 I t 2 + k (   )  d dt + 2f dt 解けそうにない？ 17 実際に解析したモデル HIP θ1 (φ1) I1 LP-A θ2 (φ2) I2 k12 I3 k23 P２ P１ d  Ge θ4 (φ4) LP-B θ3 (φ3) k34 P４ P３ 2 高中圧ﾀｰﾋﾞﾝ I 1 1 dt 2 d  + k 1 ( 1  ) 2 2 低圧Ａﾀｰﾋﾞﾝ I 2 2 dt 2 d  I4 + k 2 ( 2  3 P1  0  + )  k 1 ( 1 2 低圧Ｂﾀｰﾋﾞﾝ I 3 3 dt 2 d  + k 3 ( 3 2 発電機 I 4 4 dt 2 もっと解けそうにない？  k 3 ( 3  4 )  k 2 ( 2 1  ) P2  + P3   3)   0 +  P 4  4)   0 +  2  0 2 3 4 刻めば解ける！ 18 ２．副産物としてできた単純明快な数値解析手法最も基本的なオイラー法による数値解析微分の定義式（速度と位置の例） x(t + t )  x(t ) v(t )  lim t 0 t t v(t )  が十分小さいときの近似式変形すると x(t + t )  x(t ) t x(t + t )  x(t ) + v(t )t ある瞬間における速度が、その後の微小時間継続するとみなして、微小時間経過後の位置を求めるという考え方 19 オイラー法とは・・・直線近似 x(t + t )  x(t ) + v(t )t x(t ) ：位置この点における接線(速度)が曲線を近似していると考える t t＋t t ：時刻 20 オイラー法の計算例１分きざみで速度計の読みがわかっているとき時刻０分の時の速度を０～１分の速度だとみなす方法時刻（分） 0 1 2 3 4 5 6 7 8 9 10 速度計（m/分） 0 300 550 610 580 590 680 450 490 570 540 時間（分） 0～1 1～2 2～3 3～4 4～5 5～6 6～7 7～8 8～9 9～10 代表速度(m/分) 0 300 550 610 580 590 680 450 490 570 時間内の移動(m) 0 300 550 610 580 590 680 450 490 570 （＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）時刻（分） 0 1 2 3 4 5 6 7 8 9 10 位置(m) 0 0 300 850 1460 2040 2630 3310 3760 4250 4820 区間の速度を初期速度で代表させる「前進オイラー法」とも呼ばれる 21 別のオイラー法時刻（分） 0 1 2 3 4 5 6 7 8 9 10 速度計（m/分） 0 300 550 610 580 590 680 450 490 570 540 時間（分） 0～1 1～2 2～3 3～4 4～5 5～6 6～7 7～8 8～9 9～10 代表速度(m/分) 300 550 610 580 590 680 450 490 570 540 時間内の移動(m) 300 550 610 580 590 680 450 490 570 540 （＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）（＋）時刻（分） 0 1 2 3 4 5 6 7 8 9 10 位置(m) 0 300 850 1460 2040 2630 3310 3760 4250 4820 5360 区間の速度を終期速度で代表させる「後退オイラー法」と呼ばれる 22 運動方程式でのオイラー法の活用位置と速度の関係式 x(t + t )  x(t ) + v(t )  t 速度と加速度の関係式 v(t + t )  v(t ) + a(t )  t 運動方程式では加速度が関数として表されるので、上の２式を使って計算を進めることができる。では、一番単純な単振動で計算してみよう。 23 単振動をオイラー法で計算自然長の状態おもり質量 m x0 x ・おもりの質量：ｍ＝１［ｋｇ］・ｔ＝０における位置：ｘ(0)＝０.1［ｍ］、速度：ｖ(0)＝０［ｍ／ｓ］・ばね定数：ｋ＝１［Ｎ／ｍ］(ばねの質量は無視できる) ・ばねの力＝－ｘなので運動方程式はａ(t)＝－ｘ(t) 24 表計算ソフトに式を入れてみるａ（ｔ）＝－ｘ（ｔ） v（ｔ+Δｔ）＝ v（ｔ）＋ａ（ｔ）・Δｔ x（ｔ+Δｔ）＝ x（ｔ）＋v（ｔ）・Δｔ時刻をΔｔ進める 25 オイラー法で計算・グラフ作成ができた！？ 26 しかし、よく見ると振動が発散していく・・・刻み時間 0.01秒刻み時間 0.005秒オイラ－法では時間の経過と共に，振動が発散細かくきざめば精度は向上するが・・・やはり振幅に拡大傾向 27 オイラー法を改善してみる速度と加速度の関係式 v(t + t )  v(t ) + a(t )  t 位置と速度の関係式 x(t + t )  x(t ) + v(t )  t 速度を平均速度に v(t ) + v(t + t ) x(t + t )  x(t ) +  t 2 結果は・・・振幅の増大は改善されたが依然として発散する傾向 28 さらに改善してみる位置と速度の関係式 x(t + t )  x(t ) + v(t )  t 速度を平均速度に v(t ) + v(t + t ) x(t + t )  x(t ) +  t 2 平均速度をこのように考える 0.5v(t ) + 0.5v(t + t ) 初期速度と終期速度の比率が０．５：０．５ 29 初期速度と終期速度の最適比修正オイラー法をさらに修正するための視点として、ｘ（ｔ+Δｔ）算出に用いるｖ（ｔ）とｖ（ｔ+Δｔ）の比率に着目オイラー法修正オイラー法さらなる改善区間初期速度ｖ（ｔ）の比率区間終期速度ｖ（ｔ+Δｔ）の比率評価１０精度低い０．５０．５やや改善？？？最適比は０.４：０.６？ｏｒ０.３：０.７？ｏｒ・・・・・？？？ 30 思い切って０：１で試してみると振幅の拡大傾向がなくなった区間速度を終期速度で代表させるのがベスト 31 求めた最適比による式をながめる速度と加速度の関係式位置と速度の関係式 v(t + t )  v(t ) + a(t )  t x(t + t )  x(t ) + v(t + t )t この式は後退オイラー法の式 32 前進オイラー・後退オイラーの組み合わせ計算ステップ加速度から速度への計算速度から位置への計算計算方法 v(t + t )  v(t ) + a(t )  t x(t + t )  x(t ) + v(t + t )t オイラー法の種類前進オイラー後退オイラー｢ミックスオイラー（Ｍｉｘｅｄ-Ｅｕｌｅｒ）法」と呼ぼう 33 文献調査結果【各手法の代表的な解説内容】非減衰自由振動方程式：ｘ”＋ ( 2π) T 時間刻み幅 h=0.02T 2 ｘ = 0 に対する時刻歴応答波形（参考）破線は厳密解 h=0.05T ○オイラー法：『大きく発散。……不安定であり，精度も悪いため，実際に用いられることは少ない』 ○ホイン法（修正オイラー法）：『発散が現れる。オイラー法より精度は改善され２次の精度を有するが……あまり実用されていない』 ○後方(後退）オイラー法：『安定ではあるが，精度が悪く，実用には不向き』 h=0.1T 発散する事実がべられている程度述単純に計算を進めることができず繰り返し計算が必要に。しかも解析結果は減衰する。 34 ルンゲクッタ４次法とミックスオイラー法との比較（ばねの単振動）ＲＫ法（４次）ｔ x( t ) 0 0. 01 0. 02 0. 03 0. 04 0. 05 0. 06 0. 07 0. 08 0. 09 0. 1 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. v( t ) 1 998027 992115 982287 968583 951057 929777 904827 876307 844328 809017 0 -0. 39452 -0. 78749 -1. 17735 -1. 56256 -1. 94161 -2. 31299 -2. 67525 -3. 02695 -3. 3667 -3. 69316 kv1 kx1 -0. 394784176 -0. 39400516 -0. -0. 391671185 -0. -0. 387791464 -0. -0. 382381308 -0. -0. 375462067 -0. -0. 36706105 -0. -0. 357211411 -0. -0. 345952022 -0. -0. 333327319 -0. -0. 319387125 -0. 0 003945244 007874918 011773514 015625644 019416108 023129945 026752499 030269472 033666986 036931632 kv2 kx2 kv3 kx3 kv4 kx4 -0. 394784176 -0. 001973921 -0. 39439454 -0. 001973921 -0. 394004903 -0. 003943945 -0. 3932264 -0. 00591527 -0. 392837532 -0. 005911376 -0. 391671442 -0. 00787362 -0. 390116739 -0. 009833274 -0. 389730175 -0. 009825502 -0. 387792233 -0. 01177222 -0. 385467466 -0. 013712471 -0. 385084731 -0. 013700851 -0. 382382585 -0. 015624361 -0. 379296929 -0. 017537551 -0. 378919534 -0. 017522129 -0. 375463848 -0. 01941484 -0. 371629481 -0. 021293418 -0. 371258915 -0. 021274255 -0. 367063328 -0. 023128697 -0. 362495382 -0. 02496525 -0. 362133107 -0. 024942422 -0. 357214177 -0. 026751276 -0. 351930679 -0. 028538556 -0. 351578126 -0. 028512152 -0. 345955264 -0. 03026828 -0. 339977068 -0. 031999232 -0. 339635627 -0. 031969358 -0. 333331025 -0. 033665829 -0. 326681722 -0. 035333623 -0. 326352741 -0. 035300395 -0. 319391281 -0. 036930514 -0. 312097113 -0. 038528568 -0. 311781891 -0. 038492118 -0. 304191046 -0. 040049451 計算の手数が約１／４に圧縮，精度は同等 x[m] ミックスオイラー法ｔ x( t ) 0 0. 01 0. 02 0. 03 0. 04 0. 05 0. 06 0. 07 0. 08 0. 09 0. 1 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1 996052 988172 976391 960755 941326 918181 891411 861122 827434 790479 v( t ) 0 -0. 39478 -0. 78801 -1. 17812 -1. 56359 -1. 94288 -2. 3145 -2. 67698 -3. 0289 -3. 36886 -3. 69551 a( t ) -39. -39. -39. -38. -37. -37. -36. -35. -33. -32. -31. 4784 3226 0115 5464 9291 1621 2483 1915 9957 6658 2068 20 15 10 5 0 0 2 4 6 8 -10 -15 10 12 14 16 ｔ[s] -5 ・・・：オイラー法：ルンゲクッタ（４次）法：ミックスオイラー法 -20 35 ミックスオイラー法の限界振幅は発散しないが、位相のずれが誤差として発生する 36 手法の限界は事実上問題あり？ ○ 「振動」のような実際の問題を取り扱う場合、「位相」は「振幅」ほど重要ではない。 ○ 刻み時間を細かくすることによって改善が可能また、ﾙﾝｹﾞｸｯﾀ法との併用も可能ミックスオイラー法ルンゲクッタ法特徴活用方法併用の効果単純精度やや劣概略検討スクリーニングダブルチェックによ複雑精度良スクリーニング後の詳細検討る解析結果の確認が可能 37 何故か？（概念図で理解）【オイラー法】 t=0 加速度 a t=Δt t=2･Δt 【修正オイラー法】 t=0 t=Δt 加速度 a 速度速度 v v 変位変位 x a → v → x の段階的近似において，時 t=2･Δt x オイラー法よりは改善刻の適用は全体的に遅れ気味【ミックスオイラー法】 t=0 t=Δt 加速度 a （参考）後退オイラー法 t=2･Δt 加速度速度 v v 変位変位全体的に補正しあう効果が期待できる t=Δt t=2･Δt a 速度 x a → v 時刻の適用は遅れ v → x 時刻の適用は進み t=0 x a → v → x の段階的近似において，時刻の適用は全体的に進み気味 38 テイラー展開的な理解① 【テイラー展開】 1 1 2 n x (t + Δt）＝x (t ) + Δt  x (t ) + ( Δt )  x (t ) + ･･･ + ( Δt )  x n (t ) + ･･･ 2 n! 第１項第２項第３項第（ｎ+１）項【オイラー法】 v(t + Δt）＝v(t ) + Δt  a(t ) x( t + Δt）＝x( t ) + Δt  v( t ) テイラー展開の第２項までに相当 39 テイラー展開的な理解② 【テイラー展開】 1 1 2 n x (t + Δt）＝x (t ) + Δt  x (t ) + ( Δt )  x (t ) + ･･･ + ( Δt )  x n (t ) + ･･･ 2 n! 第１項第２項第３項第（ｎ+１）項【修正オイラー法】 v( t + Δt）＝v( t ) + Δt  a( t ) x( t + Δt）＝x( t ) + Δt 2    v( t ) + v( t + Δt ) 1 2 x( t + Δt）＝x( t ) + Δt  v( t ) + ( Δt )  a( t ) 2 テイラー展開の第３項までに相当 40 テイラー展開的な理解③ 【テイラー展開】 1 1 2 n x (t + Δt）＝x (t ) + Δt  x (t ) + ( Δt )  x (t ) + ･･･ + ( Δt )  x n (t ) + ･･･ 2 n! 第１項第２項第３項第（ｎ+１）項【ミックスオイラー法】 v( t + Δt）＝v( t ) + Δt  a( t ) x( t + Δt）＝x( t ) + Δt  v( t + Δt ) x( t + Δt）＝x( t ) + Δt  v( t ) + ( Δt )  a( t ) 2 テイラー展開第４項以下の合計を第３項に等しいとおいた式に相当 41 エピローグタービン発電機の軸ねじれの解析結果（ミックスオイラー法を活用）主要計算シート諸元入力シート出力（グラフ）算出したトルクを強度評価に活用 42 ミックスオイラー法を役立てる解析が必要となる問題の認識背景となる知識モデル化（方程式化）解析手法の活用解析ツール解析作業結果の妥当性チェックネックだった解析作業の難解さがﾐｯｸｽｵｲﾗｰ法により解消!! モデル化への意欲が増す解析結果の評価問題の解決へ 43 エピソード･･･発電機の内部相差角を測定するアルミ箔発電機レーザー装置（検出器つき）発電機端子電圧レーザー光線パルス時刻時刻ずれ発電機無負荷運転時ずれ発電機負荷運転時おわり 44