地震計の物理学3-躍動する舞台-

前へ　地震計の物理学2-道化師のダンス-

本項

次へ Coming Soon

机を揺らしてみる

机を揺らすと簡単に言ってしまったが、想像するだけでもそれなりの重労働だ。その上、与えられた関数に沿うように精密に揺らさなければいけない。とてもじゃないが人間技ではない。だが心配する必要はなかった。弊社は地震計メーカーだ。なんと自由自在に勝手に揺れてくれる超高性能な揺れる机、振動台がある。出荷される地震計の検査はもちろん、開発中の新製品の性能試験にも使われる超重要設備だ。出来るだけ外から余計な振動が伝わらないように、室内にあるにも関わらず、振動台を設置してある基礎が社屋の基礎と分離されている念の入れようである。こいつを使って揺らしていこう。もちろん頭の中で、だ。そう、ただの宣伝である。

この章から読み始めた人は地震計の物理学1を流し読みしてきてほしい。メーカーがしっかり技術的な解説をしていると思って読むと求めているものとは違う可能性がある。

どう揺らす？

基本的には前章まででやった運動方程式に振動する力の項を加えてあげればよい。問題はその力の形だ。大発見をした時に友人の肩を掴みガクガク揺さぶるように線形のジグザグな力を掛けるか、地震波の様に様々な周期振幅で複雑に揺らすか、もっとスマートにフラフープを回す様に揺らすのか。どんな力を考えるのが適切なのか、その力をどう表現するのか、とても難しい問題だが、幸いにも我々にはジョセフ・フーリエ先生がついている。彼によれば関数は三角関数の足し算で表現出来る。フーリエ変換である。波形解析時にも大活躍してくれる。三角関数の足し算で表せるなら話は簡単だ。と言うことで取り敢えずシンプルな正弦波を考えることにしよう。世の中、簡単に割り切れるものばかりではないがシンプルな事はいいことだ。

厳密に言えばそもそも本当にフーリエ変換出来るのかとか連続なフーリエ変換と離散的な級数展開の差とか、考える事は山ほどあるが、どうせ地震計はデジタルなのだ。時間的にも空間的にも離散的にしか取得できない。シンプルな足し算で良いことにしておこう。

余談1 地震計の分解能

ちなみに、特別なものを除いて大抵の地震計の分解能は24bitである。ノイズ等の影響も考えると実行分解能は21bit程度になる。今まで2Gまでしか測れなかったが4Gまで測れる高性能な地震計を導入しました！と言う文句を時々見るが最大値が2倍でも分解能は24bitのままだ。つまり精度が半分になってしまう。もちろん4Gまで測れるセンサーが高性能でないわけではない。しかし、精々2G程度しか揺れない様な場所で4Gのセンサーを導入してしまうと無駄に精度が犠牲になってしまう。観測したい場所に適した性能を選ぶことが重要なのだ。

余談2 サンプリング周波数

話が逸れたついでに時間の分解能の話もしておこう。大体の観測では100Hzでサンプリングしている様である。フーリエ先生の知恵を拝借すると、これは50Hz程度までの波を分析に使えることを意味する。センサーの周波数感度やデータが巨大になってしまう問題もあるため、より高速にサンプリングしたからといって一概にいいデータが取れるとは言えないが、観測したいものに応じて適切に設定することでよりよい分析ができるだろう。データ容量に問題がないなら取り敢えず1kHzで観測しておいて、分析時の必要に応じて間引いて100Hzに落としてもよい。なお、弊社のAccuSEIS Centoは常時1kHzでサンプリングし、収録時に設定された周波数まで間引いている。高速サンプリングの御用命がございましたらご検討いただけると幸いです。

余談3 時刻精度

ここまで来たらとことんそれておこう。時刻精度の話もしておきたい。大抵の地震計はGPS衛星が発信している時刻情報で時刻を合わせているため非常に正確で、離れた場所のデータであっても簡単に波形を並べることが出来る。時刻精度が信頼できない場合は、近距離地点であれば波の立ち上がりを揃え、そうでなければ職人技になるそうだ。
GPSに同期しているから非常に正確と述べたが、実はその正確さにも色々と差がある。まずは同期頻度だ。同期と同期の間は収録機の時計で合わせることになる。当然、その頻度が高いほど正確になるが、CPUの処理性能や消費電力等の問題から頻度を10分や30分に1度と下げて運用することがある。海底地震計などは電池や電波受信の問題で、沈める前に合わせたあとは時刻修正が出来ない。そのため省エネで精度の高い時計をつけてあげる必要がある。
また、波形データの時刻をどの時点でつけるかによっても差が出てくる。1番低コストで簡単なのは収録機にデータが格納された段階の時刻を記録することだ。しかし、これではデータ処理やAD変換によるタイムラグが生じてしまう。特にデータ処理に左右されてしまうとタイムラグにその幅のばらつきまでのってしまい、折角GPSに合わせていても精度が落ちてしまう。1番良いのはAD変換をかけるタイミングを時刻管理することだ。その点、弊社のAccuSEIS Cento は可能な限りGPSに常時同期し、AD変換のタイミングを10μ秒(同一ユニット内は1μ秒)の精度で合わせているため非常に正確である。これが言いたかった！

強制振動と共振

世の中、広告に溢れている。そのあり方は度々議論にあがり、よりよい広告が模索されているが、いずれにせよ生活と広告は切っても切り離せないものだ。この記事もやっとCMが終わったので本題に入ることにする。

今まで扱ってきた

x''+px'+qx=0

の形でかける微分方程式は斉次方程式と呼ばれる。これから扱うのは、これにさらに外力r(t)が加わった非斉次方程式

x''+px'+qx=r(t)

である。ここでr(t)は非斉次項と呼ばれ、非斉次方程式の一般解は対応する斉次方程式の一般解に、非斉次項に対応する特解を加えたものになる。これには線形性があるため、非斉次項が複数の関数の足し算で書かれていた場合、特解はそれぞれの非斉次項に対応する特解を足し合わせた物で良い。ということは、フーリエ級数を考えると代表する三角関数を一つ計算してあげるだけで大体の振動現象を考えることが出来る。なんだこれは。あまりに都合良く出来すぎていやしないか。なお詳細は省略するが、特解の線形性は非斉次項と対応する解をr_1(t):x_1(t),r_2(t):x_1(t),R(t)=r_1(t)+r_2(t):X(t)=x_1(t)+x_2(t)とでも置いて複数用意して方程式に代入してあげれば簡単に確認することができる。気が向いたらやってみてほしい。

強制振動

さて、まずはシンプルに減衰のない単振動を考えよう。バネを繋いでいる土台、つまり机をX(t)=A_0\sin{\Omega t}\:(\Omega\neq\omega)で揺らす。x(t)は机に対する相対変位だと考えよう。この時、錘は外力によって振動させられるので強制振動と呼ぶ。フーリエ級数展開を意識するなら本当はX(t)=A_0\sin{\Omega t}+B_0\cos{\Omega t}とでもしたほうが良いのだろうし、なんならラプラス変換まで導入すればもっと一般に議論できるのだが、記事の趣旨としてラプラス変換はさすがに踏み込みすぎだし、\sinだけ考えても結論は変わらないのでシンプルな形でいいことにしておこう。気が向いたら\cosも同様に計算して足してみてほしい。さて、 \frac{d^2X(t)}{dt^2}=-\Omega^2A_0\sin{\Omega t} より、運動方程式は

\begin{aligned} &\frac{d^2(x(t)+X(t))}{dt^2}\\ &=\frac{d^2x(t)}{dt^2}+\frac{d^2X(t)}{dt^2}\\ &=\frac{d^2x(t)}{dt^2}-\Omega^2A_0\sin{\Omega t}\\ &=-\frac{k}{m}x(t)\\ &\therefore \frac{d^2x}{dt^2}+\omega^2x= \Omega^2A_0\sin{\Omega t}\quad\Bigl(\omega^2=\frac{k}{m}\Bigr)\\ \end{aligned}

となり、非斉次方程式が求められた。一応、斉次方程式の一般解を確認しておこう。

\begin{aligned} x(t)&=c_1e^{i\omega t}+c_2e^{-i\omega t}\\ &=c_1(\cos{\omega t}+i\sin{\omega t})+c_2(\cos{\omega t}-i\sin{\omega t})\\ &=a\cos{\omega t}+b\sin{\omega t}\\ &\bigl(a=(c_1+c_2),b=i(c_1+c_2)\bigr) \end{aligned}

であった。特解は地道に積分して求める方法もあるのだが、この場合、当たりをつけて代入して確認する方が簡単なのでそうする。x=c\sin{\Omega t}+d\cos{\Omega t}とおいて運動方程式に代入すると

\begin{aligned} -\Omega^2&(c\sin{\Omega t}+d\cos{\Omega t})+\omega^2(c\sin{\Omega t}+d\cos{\Omega t})\\ =&(\omega^2-\Omega^2)c\sin{\Omega t}+ (\omega^2-\Omega^2)d\cos{\Omega t}\\ =&\Omega^2A_0\sin{\Omega t} \end{aligned}

よりc=\Omega^2A_0/(\omega^2-\Omega^2),d=0となる。したがって一般解は

\begin{aligned} x(t)=a\cos{\omega t}+b\sin{\omega t}+\frac{\Omega^2A_0}{\omega^2-\Omega^2}\sin{\Omega t} \end{aligned}

となる。さて、初期条件を与えてどんな運動をするか見ていく。とりあえずx(0)=x_0,x'(0)=v_0としておこう。

\begin{aligned} x(0)&=a\cos{(\omega・0)}+b\sin{(\omega・0)}+\frac{\Omega^2A_0}{\omega^2-\Omega^2}\sin{(\Omega・0)}\\ &=a\\ &=x_0\\ \\ x'(0)&=-\omega a\sin{(\omega・0)}+\omega b\cos{(\omega・0)}+\Omega\frac{\Omega^2A_0}{\omega^2-\Omega^2}\cos{(\Omega・0)}\\ &=\omega b+\frac{\Omega^3A_0}{\omega^2-\Omega^2}\\ &=v_0\\ \\ b&=\frac{1}{\omega}・\Bigl(v_0-\frac{\Omega^3A_0}{\omega^2-\Omega^2}\Bigr)\\ \end{aligned}

より

\begin{aligned} x(t)=x_0\cos{\omega t}+\frac{1}{\omega}・\Bigl(v_0-\frac{\Omega^3A_0}{\omega^2-\Omega^2}\Bigr)\sin{\omega t}+\frac{\Omega^2A_0}{\omega^2-\Omega^2}\sin{\Omega t} \end{aligned}

x_0=0,\;v_0=0とすると

\begin{aligned} x(t)&=\frac{\Omega^2A_0}{\omega^2-\Omega^2}\Bigl(\sin{\Omega t}-\frac{\Omega}{\omega}\sin{\omega t}\Bigr)\\ &=\frac{A_0}{\frac{\omega^2}{\Omega^2}-1}\Bigl(\sin{\Omega t}-\frac{\Omega}{\omega}\sin{\omega t}\Bigr) \end{aligned}

となる。

\Omega\to\omegaと、机が揺れる周波数を固有周波数に近づけていくと振幅が増大していく。これを共振と呼ぶ。この時、うなりと呼ばれる現象が見られる。ギターなど楽器の調律をするときに利用するあれだ。基準音になる音叉を鳴らし、弦を弾いてウォンウォン鳴るうなりを聴く。調律をして音が一致していくとうなりの周期が伸びていき、一致するとうなりが消え、通り過ぎるとまたうなり出すので良い感じのところで合わせるのだ。

https://www.desmos.com/calculator/44lb8ir9ya

 

また、初期条件をx_0=0,\;v_0=\frac{\Omega^3A_0}{\omega^2-\Omega^2}と都合よく設定してあげると斉次解の部分が消え特解だけになる。

\begin{aligned} x(t)=\frac{\Omega^2A_0}{\omega^2-\Omega^2}\sin{\Omega t} \end{aligned}

入力された振幅A_0、角周波数\Omegaの振動が周波数はそのままに振幅が\frac{\Omega^2}{\omega^2-\Omega^2}だけ増幅されて振動する。下図は入力角周波数に対する応答振幅のグラフである。固有角周波数\omegaに近いほど振幅が大きくなるだけでなく、\omegaを境に位相がずれて振幅がひっくり返っている。また、固有周波数に対して入力周波数が十分に小さいとき、錘はほぼ動かず静止し、逆に十分早く振ってあげると入力振幅に近づいていく。

実際にはそんな都合のいい初期条件を与える事は困難だが、前回やった減衰のある場合を思い出してほしい。十分に時間が経てば斉次解は減衰し無視できた。つまり減衰のある場合は特解こそが系の挙動を決める重要な項になるのだ。早速減衰のある場合を見ていきたいのだが、その前に今回除外した\Omega=\omegaの場合を見ておこう。

共鳴振動

机を揺らす周波数を固有周波数に近づけると振幅が極端に増大することがわかった。それでは完全に一致させた時どうなるか見ていこう。運動方程式はさっきと同様に

\frac{d^2x}{dt^2}+\omega^2x= \omega^2A_0\sin{\omega t}

とする。非斉次項の角周波数が\omegaであることに注意して欲しい。
特解も同様にx=c\sin{\omega t}+d\cos{\omega t}としたいところだが、これは斉次解になっており特解にならない。そこで時間tをかけてx=t(c\sin{\omega t}+d\cos{\omega t})としてみると

\begin{aligned} x'(t)&=(c\sin{\omega t}+d\cos{\omega t})+\omega t(c\cos{\omega t}-d\sin{\omega t})\\ &=-(\omega dt-c)\sin{\omega t}+(\omega ct+d)\cos{\omega t}\\ \\ x''(t)&=-\omega d\sin{\omega t}-\omega(\omega dt-c)\cos{\omega t}\\ &\qquad+\omega c\cos{\omega t}-\omega(\omega ct+d)\sin{\omega t}\\ &=-2\omega(d\sin{\omega t}-c\cos{\omega t})\\ &\qquad-\omega^2t(c\sin{\omega t}+d\cos{\omega t}) \end{aligned}

より、運動方程式は

-2\omega(d\sin{\omega t}-c\cos{\omega t})=\omega^2A_0\sin{\omega t}

となるのでd=-\frac{\omega A_0}{2},c=0を得る。したがって一般解は

x(t)=(a-\frac{\omega A_0}{2}t)\cos{\omega t}+b\sin{\omega t}

である。時間tをかけると特解になると言うのは少し狐につままれた様な感じもするが、矛盾なく求まるのだから仕方ない。初期条件x(0)=x_0,x'(0)=v_0を与えると、

\begin{aligned} x(0)&=(a-\frac{\omega A_0}{2}・0)\cos{(\omega・0)}+b\sin{(\omega・0)}\\ &=a\\ &=x_0\\ \\ x'(0)&=-\frac{\omega A_0}{2}\cos{(\omega・0)}-\omega(a-\frac{\omega A_0}{2}・0)\sin{(\omega・0)}+\omega b\cos{(\omega・0)}\\ &=-\frac{\omega A_0}{2}+\omega b\\ &=v_0 \end{aligned}

よりa=x_0,b=\frac{2v_0+A_0\omega}{2\omega}なので解は

\begin{aligned} x(t)=\Bigl(x_0-\frac{\omega A_0}{2}t\Bigr)\cos{(\omega t)}+\frac{2v_0+A_0\omega}{2\omega}\sin{(\omega t)} \end{aligned}

振幅に時間tに比例する項を含んでいるため、時間経過とともに振幅が増大していく。下図はx_0=0,v_0=0のときのグラフである。

 

減衰のある強制振動

減衰がある場合も同様に机を振動させよう。非斉次項をr=\omega^2A_0\sin{\Omega t}としたいところだが、減衰振動の斉次解で\Omegaを使ってしまっている。こう言うパラメータに適切な文字を当てていくのは実に難しい。手計算でノートに落書きしているときは取り敢えずなんでもいいが、人に読んでもらおうとなるとわかりやすい方がいい。今、非斉次項はフーリエ展開された地震動を想定しているので序数のnをつけて周波数を\omega_nと置いてあげることにしよう。合わせてA_0もA_nとする。この場合、遡って単振動の部分も書き換えるのが筋なのであろうが、話のネタにもなったわけだしあえてそのままにしよう。そのくらいのノリの記事である。

さて、考える運動方程式は

\begin{aligned} \frac{d^2x}{dt^2}+\nu\frac{dx}{dt}+\omega^2x=\omega_n^2A_n\sin{\omega_n t} \end{aligned}

である。対応する斉次方程式の解を確認しておこう。

\begin{aligned} x(t)=&\left\{ \begin{array}{} (a\cos{\Omega t}+b\sin{\Omega t})e^{-\frac{\nu}{2}t} &(\nu^2-4\omega^2<0)\\ (at+b)e^{-\frac{\nu}{2}t}&(\nu^2-4\omega^2=0)\\ (ae^{\gamma t}+be^{-\gamma t})e^{-\frac{\nu}{2}t} &(\nu^2-4\omega^2>0)\\ \end{array} \right.\\ \\ &\Biggl(\:\Omega=\sqrt{\omega^2-\frac{\nu^2}{4}}\:,\:\gamma=\sqrt{\frac{\nu^2}{4}-\omega^2}\:\Biggr) \end{aligned}

今回も特解を求めた後、それぞれの場合について挙動を検証していくことになる。少々大変そうだが、幸いにも特解は共通の一つで済む。単振動の時の様に\omega=\omega_nを分けて考える必要がないのも救いである。今回もx(t)=c\sin{\omega_nt}+d\cos{\omega_nt}として運動方程式に代入し整理することで定数c,dを求めていこう。

\begin{aligned} x(t)=&c\sin{\omega_nt}+d\cos{\omega_nt}\\ x'(t)=&-\omega_nd\sin{\omega_nt}+\omega_nc\cos{\omega_nt}\\ x''(t)=&-\omega_n^2c\sin{\omega_nt}-\omega_n^2d\cos{\omega_nt} \end{aligned}

より

\begin{aligned} &(-\omega_n^2c-\omega_n\nu d+\omega^2c)\sin{\omega_nt}+(-\omega_n^2d-\omega_n\nu c+\omega^2d)\cos{\omega_nt}\\ &\qquad=\omega_n^2A_n\sin{\omega_nt} \end{aligned}

したがって

\left\{ \begin{array}{} (\omega^2-\omega_n^2)c-\omega_n\nu d=&\omega_n^2A_n\\ (\omega^2-\omega_n^2)d+\omega_n\nu c=&0 \end{array} \right.

となる。見るからに厄介だが、頑張って解いていこう。

\begin{aligned} c=-\frac{\omega^2-\omega_n^2}{\omega_n\nu}d \end{aligned}

より

\begin{aligned} &-\frac{(\omega^2-\omega_n^2)^2}{\omega_n\nu}d -\frac{(\omega_n\nu)^2}{\omega_n\nu} d=\omega_n^2A_n\\ \Leftrightarrow&\quad \bigl\{(\omega^2-\omega_n^2)^2+\omega_n^2\nu^2\bigr\}d=-\omega_n^2A_n・\omega_n\nu\\ \therefore&\quad d=-\frac{\omega_n^2A_n・\omega_n\nu}{(\omega^2-\omega_n^2)^2+\omega_n^2\nu^2} \end{aligned}

また

\begin{aligned} c=\frac{\omega_n^2A_n・(\omega^2-\omega_n^2)}{(\omega^2-\omega_n^2)^2+\omega_n^2\nu^2} \end{aligned}

よって特解はB^2=(\omega^2-\omega_n^2)^2+\omega_n^2\nu^2として

\begin{aligned} x(t)&=\frac{\omega_n^2A_n}{B^2} \bigl\{(\omega^2-\omega_n^2)\sin{\omega_nt}-\omega_n\nu\cos{\omega_nt}\bigr\}\\ &=\frac{\omega_n^2A_0}{B} \Biggl(\frac{\omega^2-\omega_n^2}{B}\sin{\omega_nt} -\frac{\omega_n\nu}{B}\cos{\omega_nt}\Biggr)\\ &=\frac{\omega_n^2A_n}{B}(\sin{\omega_nt}\cos{\delta}-\cos{\omega_nt}\sin{\delta})\qquad\Bigl(\tan{\delta}=\frac{\omega_n\nu}{\omega^2-\omega_n^2}\Bigr)\\ &=A\sin{(\omega_nt-\delta)}\qquad\Bigl(A=\frac{\omega_n^2A_n}{B}\Bigr) \end{aligned}

なんとか求められた。計算はややこしかったが特解はシンプルな正弦波の形になっている。一応\omega=\omega_nの時を考えると、

\begin{aligned} x(t)=-\frac{\omega}{\nu}A_n\cos{\omega t} \end{aligned}

と幾分簡単になる。減衰し0に近づく斉次解に特解を足すことを考えると、初期条件の影響を受ける斉次解は時間経過につれ減衰し十分に小さくなり特解の示す単振動に落ち着くことがわかる。つまり、電源を入れ観測を始めてから減衰項\nuで決まる十分な時間が経過した後は特解にのみ注目してあげれば良いと言える。一応初期条件を与え一般解をみていこう。

・\nu^2-4\omega^2<0の時、一般解は

x(t)=(a\cos{\Omega t}+b\sin{\Omega t})e^{-\frac{\nu}{2}t}+A\sin{(\omega_nt-\delta)}

初期条件x(0)=x_0,x'(0)=v_0を与えると

\begin{aligned} x(0)&=(a\cos{(\Omega・0)}+b\sin{(\Omega・0)})e^{-\frac{\nu}{2}・0}+A\sin{(\omega_n・0-\delta)}\\ &=a-\frac{\omega_n\nu A}{B}\\ &=x_0\\ \\ x'(0)&=\Bigl\{(-\frac{\nu}{2}a+\Omega b)\cos{(\Omega・0)} +(-\frac{\nu}{2}b-\Omega a)\sin{(\Omega・0)}\Bigr\}e^{-\frac{\nu}{2}・0}\\ &\qquad\qquad\qquad+\omega_nA\cos{(\omega_n・0-\delta)}\\ &=-\frac{\nu}{2}a+\Omega b+\frac{(\omega^2-\omega_n^2)\omega_nA}{B}\\ &=v_0 \end{aligned}

より

\begin{aligned} a&=x_0+\frac{\omega_n\nu A}{B}\\ \\ b&=\frac{\nu x_0}{2\Omega}+\frac{\omega_n\nu^2A}{2B\Omega}-\frac{(\omega^2-\omega_n^2)\omega_nA}{B \Omega}+\frac{v_0}{\Omega}\\ &=\frac{\omega_nA}{\Omega B}\bigl\{\frac{\nu^2}{2}-(\omega^2-\omega_n^2)+\frac{\nu Bx_0+2Bv_0}{2\omega_nA}\bigr\} \end{aligned}

よって解は

\begin{aligned} x(t)&=\Bigl[(x_0+\frac{\omega_n\nu A}{B})\cos{\Omega t}+\frac{\omega_nA}{\Omega B}\bigl\{\frac{\nu^2}{2}-(\omega^2-\omega_n^2)+\frac{\nu Bx_0+2Bv_0}{2\omega_nA}\bigr\}\sin{\Omega t}\Bigr]e^{-\frac{\nu}{2}t}\\ &\qquad\qquad\qquad+A\sin{(\omega_nt-\delta)} \end{aligned}

 

・\nu^2-4\omega^2=0の時、一般解は

\begin{aligned} x(t)=(at+b)e^{-\frac{\nu}{2}t} +A\sin{(\omega_nt-\delta)} \end{aligned}

初期条件x(0)=x_0,x'(0)=v_0を与えると

\begin{aligned} x(0)&=(a・0+b)e^{-\frac{\nu}{2}・0} +A\sin{(\omega_n・0-\delta)}\\ &=b-\frac{\omega_n\nu A}{B}\\ &=x_0\\ \\ x'(0)&=ae^{-\frac{\nu}{2}・0} -\frac{\nu}{2}(a・0+b)e^{-\frac{\nu}{2}・0} +\omega_nA\cos{(\omega_n・0-\delta)}\\ &=a-\frac{\nu}{2}b+\frac{(\omega^2-\omega_n^2)\omega_nA}{B}\\ &=v_0 \end{aligned}

より

\begin{aligned} b&=x_0+\frac{\omega_n\nu A}{B}\\ &=\frac{\omega_nA}{B}\bigl(\nu +\frac{Bx_0}{\omega_nA}\\ \\ a&=\frac{\nu x_0}{2}+\frac{\omega_n\nu^2A}{2B}-\frac{(\omega^2-\omega_n^2)\omega_nA}{B}+v_0\\ &=\frac{\omega_nA}{B}\bigl\{\frac{\nu^2}{2}-(\omega^2-\omega_n^2)+\frac{B\nu x_0+2Bv_0}{2\omega_nA}\bigr\} \end{aligned}

よって解は

\begin{aligned} x(t)= \frac{\omega_nA}{B}\Bigl[\bigl\{\frac{\nu^2}{2}-(\omega^2-\omega_n^2)+\frac{B\nu x_0+2Bv_0}{2\omega_nA}\bigr\}t+(\nu+\frac{Bx_0}{\omega_nA})\Bigr]e^{-\frac{\nu}{2}t}+A\sin{(\omega_nt-\delta)} \end{aligned}

 

・\nu^2-4\omega^2>0の時、一般解は

\begin{aligned} x(t)=(ae^{\gamma t}+be^{-\gamma t})e^{-\frac{\nu}{2}t}+A\sin{(\omega_nt-\delta)} \end{aligned}

初期条件x(0)=x_0,x'(0)=v_0を与えると

\begin{aligned} x(0)&=(ae^{\gamma・0}+be^{-\gamma・0})e^{-\frac{\nu}{2}・0}+A\sin{(\omega_n・0-\delta)}\\ &=a+b-\frac{\omega_n\nu A}{B}\\ &=x_0\\ \\ x'(0)&=\gamma(ae^{\gamma・0}-be^{-\gamma・0})e^{-\frac{\nu}{2}・0} -\frac{\nu}{2}(ae^{\gamma・0}+be^{-\gamma・0})e^{-\frac{\nu}{2}・0} +\omega_nA\cos{(\omega_n・0-\delta)}\\ &=\gamma(a-b)-\frac{\nu}{2}(a+b) +\frac{(\omega^2-\omega_n^2)\omega_nA}{B}\\ &=v_0 \end{aligned}

より

\begin{aligned} a+b&=x_0+\frac{\omega_n\nu A}{B}\\ &=\frac{\omega_nA}{B\gamma}(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)\\ \\ a-b&=\frac{\nu}{2\gamma}\frac{\omega_nA}{B\gamma}(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)-\frac{(\omega^2-\omega_n^2)\omega_nA}{B\gamma}+v_0\\ &=\frac{\omega_nA}{B\gamma}\Bigl\{\frac{\nu}{2\gamma}(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)-(\omega^2-\omega_n^2)+\frac{B\gamma v_0}{\omega_nA}\Bigr\}\\ \\ a&=\frac{\omega_nA}{2B\gamma}\Bigl\{(1+\frac{\nu}{2\gamma})(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)-(\omega^2-\omega_n^2)+\frac{B\gamma v_0}{\omega_nA}\Bigr\}\\ \\ b&=\frac{\omega_nA}{2B\gamma}\Bigl\{(1-\frac{\nu}{2\gamma})(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)+(\omega^2-\omega_n^2)-\frac{B\gamma v_0}{\omega_nA}\Bigr\}\\ \end{aligned}

よって解は

\begin{aligned} x(t)&= \frac{\omega_nA}{2B\gamma}\Bigl[\Bigl\{(1+\frac{\nu}{2\gamma})(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)-(\omega^2-\omega_n^2)+\frac{B\gamma v_0}{\omega_nA}\Bigr\}e^{\gamma t}\\ &\quad+\Bigl\{(1-\frac{\nu}{2\gamma})(\frac{B\gamma x_0}{\omega_nA}+\gamma\nu)+(\omega^2-\omega_n^2)-\frac{B\gamma v_0}{\omega_nA}\Bigr\}e^{-\gamma t}\Bigr]e^{-\frac{\nu}{2}t}\\ &\qquad\qquad+A\sin{(\omega_nt-\delta)} \end{aligned}

ガシガシと力技で計算してみたわけだが、どの場合も大きな差はなくe^{-\frac{\nu}{2}t}ですぐに減衰し、特解部分のA\sin{(\omega_nt-\delta)}に落ち着くことがわかる。下図はx_0=0,v_0=0のときのグラフである。実践が一般解、点線が特解、赤線が振幅を1/10にした入力振動だ。desmosなのでパラメータをいじって遊んでみてほしい。

 

https://www.desmos.com/calculator/yzorhmcuzw

口数が減ったことから察した読者もおられるかも知れないが、流石にこれだけの式を\TeX表記で書いていくのは疲れる。地震計としてはほとんど寄与しない斉次解部分を含む一般解の初期値解ともなると尚更だ。計算間違いや打ち損じでエラーが出たときなどどこをどう直していいやら途方に暮れるばかりである。日々ソースコードとにらめっこしてデバッグしているプログラマーには尊敬の念が絶えない。設計課の方を向き拝んでおく。今度爪の垢でももらってこよう。書くのも大変だが読むのも大変だ。ここまで読んでくれたことを嬉しく思うし、何かの助けになることを願ってやまない。なお、例によってアニメーションは単振動のグラフでY(x)=の後を書き換えることでみることができる。

初期値はほとんど寄与しないと言ったが減衰が小さく固有周期の長いセンサーの場合は少し問題になる。現場で設置してからちゃんと動いているかテスト確認できる様になるまでの待ち時間が長くなってしまうのだ。長時間待ってちょっとまずいぞと再調整する羽目になったら、安定するまで再度長時間待たなければならない。検出器不具合の現地調査ともなればどれほど待たなければならないのか。地震観測は作る人、測る人だけでなく、現場で設置調整する工事部隊の多大な貢献があるのだ。設置を疎かにしては正確な地震観測は出来ない。

さて、次回は今回得られた速度抵抗のある強制振動の性質を実際の地震計と照らし合わせながら色々とみていくことにしよう。いよいよ大詰めである。

前へ　地震計の物理学2-道化師のダンス-

本項

次へ Coming Soon