LCR and Q
Q
とある方からLCの共振回路とQについてのお題が出された。某ページを見ていたら共振回路とQに言及した説明が有り、疑問が有ったそうな。色々調べてそんなのEXCELで簡単に計算できるよとうそぶいている人も見つけたとの事。
何でもEXCELで解決するのはどんなものだろう。その人には、
akae@pen:/tmp$ apt show ngspice Description: Spice circuit simulator NGSpice is a mixed-level/mixed-signal circuit simulator. Its code is based on three open source software packages: Spice3f5, Cider1b1 and Xspice.
こういうシュミレータを紹介した。
オイラーもやってみるかな。使う道具は何にする？ pythonは論外としてoctaveとかjuliaとかかな。最近lispに触っていないので、 思い出したように Maxima の使い方 とか、Maximaのまとめ を取り出してきて、ごにょごにょやってみる。
intro
分かり易いようにLCRの直列接続を考える。Lに直列に入っているRは、コイルの抵抗分を模している。Cの等価回路に出て来るリーク抵抗は、面倒なので無視（無限大を仮定）。
(%i22) z(f) := r + (2 * %pi * f*l) * %i - %i / (2 * %pi *f*c);
- %i
(%o22) z(f) := r + (2 %pi f l) %i + ---------
2 %pi f c
(%i24) r: 10;
(%o24) 10
(%i25) l: 2e-6;
(%o25) 2.0E-6
(%i26) c: 100e-12;
(%o26) 1.0E-10
(%i28) float( z(3e6) );
(%o28) 10.0 - 492.817365129907 %i
(%i29) float( z(13e6) );
(%o29) 40.93593868521128 %i + 10.0
(%i31) float( 1 / (2 * %pi * sqrt( c * l)));
(%o31) 1.125395395196382E+7
(%i42) plot2d( abs(z(f)), [f, 5e6, 20e6]);
(%o42) [/tmp/maxout2842.gnuplot_pipes]
インピーダンスは、周波数の関数だ。抵抗は10オーム、コイルは2uH、コンデンサーは100pFとしよう。
この条件で、3MHzと13MHzの時のインピーダンスを求めてみる。この範囲内で虚部が反転してるので、どこかに共振周波数が有るはず。アマチュア無線の国家試験に出て来る式で、共振周波数を求めると、11.23MHzとなった。グラフに書いてみよう。
計算結果は複素数になるんで、それをそのままプロットしようとすると、
plot2d: expression evaluates to non-numeric value everywhere in plotting range. plot2d: nothing to plot.
そんな無理な事は出来ませんと言われる。仕方なしに絶対値を取って表示させたよ。
%i12) plot2d([realpart(z(f)), imagpart(z(f))], [f, 5e6, 20e6]); (%o12) [/tmp/maxout96159.gnuplot_pipes]
調べた限りでは、一発で複素数のグラフを書く事は出来ないけど、上記のように、実部、虚部に分解して同時に表示する事は可能だ。
参考までに、インピーダンスの絶対値も載せたグラフを置いておく。 impedance.svg
ここまでは、導入部だな。
draw graph
maximaでは方程式を解いてくれるsolveが有る。でも、無から有は生まれない。正しい式を立てられてこそ、maximaが生きてくる。
さて問題は、ある周波数に共振するLCの組み合わせは無数にある事だ。共振回路の特性を表すものにQと呼ばれるものがある。LCの組み合わせによってQが変化するのかと言う問いだ。
Qは、共振回路の尖鋭度とも言えるもの。定義は、共振周波数fcから3dB(1/sqrt(2))下がった所の周波数flとfhによって決まる。(但し、fl<fc, fh>fc とする）
q = fc /(fh - fl)
平たく言うと、共振周波数をバンド幅で除した値だ（qは、クォリティーの頭文字）。 このままでは、このqな式は使い道が無い。
仮に1Vの高周波電源に、インピーダンスz(f)を接続したら、そこを流れる電流はオームの法則で、1/z(f) になるな。この電流の周波数変化を図にすれば、目測でｑの大小が分かるはず（ですよね）。
で、このインピーダンスの逆数は アドミッタンス（単位はシーメンス）って言うんだななんて事を思い出しながら、雑な関数を定義する。なお、インピーダンスは電流の流れにくさを表す指標。アドミッタンスは、電流の流れ易さを表しているって事だ。（アドミッタンスの実部はコンダクタンス、虚部はサセプタンスって名前がついてるけど、面倒なこったい。まて、名前が付いてて便利に使う、高い周波数領域に住んでる人がいるよ）
(%i31) g(l,c,f) := 1 / (r + (6.28*f*l - 1/(6.28*f*c))* %i );
1
(%o31) g(l, c, f) := ----------------------------
1
r + (6.28 f l - --------) %i
6.28 f c
勿論、6.28は、2 * %pi の積り。この%piが残った式で、値を計算すると、後でfloatしなくちゃいけなくなるんで、生活の知恵を出したよ。
(%i33) plot2d([abs(g(2e-6,100e-12,f)), abs(g(1e-6,200e-12,f))], [f, 5e6, 20e6]); (%o33) [/tmp/maxout732.gnuplot_pipes]
qcheck.svg に、グラフを載せた。青色の方は、2uH,100pFの組み合わせ。赤色は1uH,200pFの同調回路。抵抗は10Rとして共通。同調周波数では、1/10 = 100mA の電流が流れている事が確認出来る。
結論は、LCの組み合わせで、Lが大の方がQは高くなるって事だな。
solve
図的な解法は上で済んだけど、Qを計算で求められないかな？ 計算の都合上、今回はRを１オームと固定しておく。そうすれば、同調点で１Aの電流が流れる事になり、式が建てやすい。
Qの定義からfl,fhを表す値は、1/sqrt(2) = 0.7071… 丸めてしまって0.7としておこう。 maximaには方程式を解くsolveって関数が用意されてる。結果を数値で欲しいので、floatで数値化すればよい。ただfloatを使うとやたら細かい数字が出て来るので、桁数制御の出来る、bfloatを使う事にする。fpprecで3桁出力になるようにした。
(%i36) bfloat(solve(abs(g(2e-6,100e-12,f)) = 0.7, f)), fpprec: 3;
rat: replaced -0.7 by -7/10 = -0.7
rat: replaced -1.592356687898089E+9 by -250000000000/157 = -1.592356687898089E+9
rat: replaced 1.256E-5 by 157/12500000 = 1.256E-5
! 2 !
5.1b-10 !2.46b4 %i f + 1.96b9 f - 3.13b18 %i!
(%o36) [---------------------------------------------- = 1.43b0]
abs(f)
あれれ、これ答えになってないじゃん。苦肉の策としてabsを外してみるか。
(%i37) bfloat(solve( g(2e-6,100e-12,f) = 0.7, f)), fpprec: 3;
rat: replaced -0.7 by -7/10 = -0.7
rat: replaced -1.592356687898089E+9 by -250000000000/157 = -1.592356687898089E+9
rat: replaced 1.256E-5 by 157/12500000 = 1.256E-5
(%o37) [f = 9.1b-4 %i (1.24b10 %i - 1.88b7),
f = - 9.1b-4 %i (1.24b10 %i + 1.88b7)]
答えが２つ出て来て期待したんだけど、
(%i40) float(abs (9.1b-4 * %i * (1.24b10 * %i - 1.88b7))); (%o40) 1.128401296896028E+7 (%i41) float(abs (-9.1b-4 * %i * (1.24b10 * %i + 1.88b7))); (%o41) 1.128401296896028E+7
何故か、共振周波数が結果として得られた。何か勘違いしてるかな？
公式有ったぞ
ここまで我流でやってきたけど、行き詰まってしまったので、文献を当たってみる。
有名なので、わんさかと出てくるな。
(%i1) cq(l,c,r) := sqrt( l / c) / r;
l
sqrt(-)
c
(%o1) cq(l, c, r) := -------
r
(%i2) cq(2e-6,100e-12,1.0);
(%o2) 141.4213562373095
(%i3) cq(1e-6,200e-12,1.0);
(%o3) 70.71067811865476
直列共振回路のQを計算してみた。Calc Q だから略して関数名は CQね。別に狙った訳ではない。
で、CQにかかせないのがアンテナ。オイラーにとっちゃ悩みの種。伸縮自在の釣竿アンテナを作ろうと思って、釣具店でグラスロッドを探し回ったけで、手に入らず。
今はみんなカーボンに移行してしまったそうだ。カーボンなら導電性が有るから何とかならないか？ そんなのやめて輪だよ。
by hand
上でもう結論が出てしまったけど、そこに至るまで同調カーブから、値を読み取って、そこからQを計算するスクリプトを書いていたんだ。折角だから載せておく。
# Calc Q
def cq(l,c,r)
step = 0.01e6 # scan step
fc = 1 / (6.28 * Math.sqrt(l * c))
th = 0.7 / r # detect level
xl = 6.28 * l
xc = 6.28 * c
fh = fc + step
while (1 / Complex(r, xl * fh - 1 / (xc * fh))).abs > th do
fh += step
end
fl = fc - step
while (1 / Complex(r, xl * fl - 1 / (xc * fl))).abs > th do
fl -= step
end
fc / (fh - fl)
end
puts cq(2e-6, 100e-12, 1.0)
puts cq(1e-6, 200e-12, 1.0)
久しぶりにruby語です。やってる事は単純で、同調点から高い方にスキャンしながら、3dB下がった周波数(fh)を見つける。次は低い方に向かって見つける。後は定義通りに計算。
sakae@pen:/tmp$ ruby cq.rb 112.59662120804897 62.553678448916095
理論値と大分違うな。
sakae@pen:/tmp$ irb irb(main):001:0> require './cq' : => true irb(main):002:0> cq(4e-6, 50e-12, 10.0) => 26.808719335249755
rubyのreplであるirbから、新しい計算も出来る。
MLA
上でやった直列回路、Rに物凄く敏感だった(反比例ですから)。Rが小さければ、非常に大きな電流が流れるはず。 最近、アマチュア無線業界ではマグネチック・ループ・アンテナ（略してMLA）がもてはやされている。魅力なのは、かさばらない事。狭い敷地しかない日本で流行って当然。
電磁波は、電界、磁界、電界、磁界って連なりで遠くまで届く。強い磁界を発生出来れば強い電波になる。
磁界の強さはATに比例する。小学校の理科の時間に倣う。電磁石の強さは、流れる電流（A)と巻数（T）の積になるよってやつだ。
直列共振でガンガン電流を流しましょ。巻数が多いと作るの大変だからフラフープみたいに１巻でいいよって妥協。そのうちに、誰かが常温超伝導物質を発明してくれるでしょ。そしたら、電流は天井知らずで流せるよ。しめしめ。
まて、コンデンサーの事を忘れていないか。有名な式、Q=CVってのを思い出せ！ ここに出て来たQは、上に出てきたQとは全く別物。電荷量を表している。雑に言うと電流だな。
電流が無限大になったら、コンデンサーにかかる電圧も無限大っぽくなる。そんなコンデンサー（しかも同調を取らなきゃいけないのでバリコン）なんて有るのか？
夢は潰えますなあ。