シミュレーションのトビラ

シミュレーション技術に関する個人的メモ

EXCELでスピーカーシミュレーション(9)

9.ダブルバスレフ型エンクロージャ

 以前の記事で位相反転型(バスレフ型)エンクロージャの周波数特性をシミュレーションしました。世の中にはポートと容積をさらにくっつけて二段階にしたダブルバスレフ型というのもあります。また、さらにくっつけて三段階にしたトリプルバスレフ型を自作したなんて人もいるかもしれません。
 理論上は、ポートと容積をつなげるのは、機械振動としてはバネと重りを直列につなぐのと等価です。また電気回路としてはインダクタとコンデンサを直列に繋いだものを、並列に重ねていくのと等価です。計算するとわかりますが、ポートと容積を繋げた分だけ、複数の周波数で共振が起こります。ひとつめのポートの共振に加え、追加したポートによってさらに低い周波数でも共振させることで、シングルバスレフよりもさらに低音域に周波数特性を伸ばすことを期待する形式が、ダブルバスレフだったりトリプルバスレフだったりします。
 ですので、シングルバスレフの考え方を発展させれば、ダブルバスレフやトリプルバスレフも計算でき、理論上はいくつでも繋げられます。ただ、スピーカーキャビネットの大きさとしては、やってもダブルバスレフまでが限度でしょうね。
 ダブルバスレフの構造は下図のようになります。シングルバスレフの構造に対して、キャビネットが2つに仕切られ、その間をポートで繋いだ構成です。

ダブルバスレフ型エンクロージャの構造

9-1.ダブルバスレフの等価回路

 ダブルバスレフの等価回路は、下図のようになります。

ダブルバスレフ型の等価回路
 シングルバスレフの等価回路に対して、さらにインピーダンスが並列で加わった形となります。ドライバーユニットから背面の第1容積内に放射された音波は、第1ポートを介して第2容積内に伝わり、更に第2ポートを介して外部に放射されます。第1容積のバネと第2容積のバネが第1ポート内の空気の質量でつながる形となるため、第2容積は第1ポートに対して並列で加わるのです。

9-2.等価回路の解き方

 さて、この等価回路をキルヒホッフの第2法則を適用して解いてみましょう。この等価回路では粒子速度は3つの閉回路に流れる形となるため、等価回路から連立方程式を立てると、

\begin{cases} F &= z_1V_1+z_2(V_1-V_2) \\ 0 &= z_2(V_2-V_1) + z_3V_2 + z_4(V_2-V_3) \\ 0 &= z_4(V_3-V_2) + z_5V_3 \end{cases}

となります。
 これを整理すると
\begin{cases} F &= (z_1+z_2)V_1 \hspace{48pt} -z_2V_2 \\ 0 &= \hspace{20pt} -z_2V_1+(z_2+z_3+z_4)V_2 \hspace{28pt} -z_4V_3 \\ 0 &= \hspace{20pt} -z_4V_1 \hspace{82pt} +(z_4+z_5)V_3 \end{cases}

となります。今後の計算をしやすくするために、
 \begin{align} z_1+z_2&=A \\ z_2+z_3+z_4&=B \\ z_4+z_5&=C\end{align}

としてまとめておくと、
\begin{cases} F &= \hspace{9pt} AV_1-z_2V_2 \\ 0 &= -z_2V_1 \hspace{2pt} +BV_2 -z_4V_3 \\ 0 &= -z_4V_1 \hspace{34pt} +CV_3 \end{cases}

となります。
 さて、この方程式はV_1V_2V_3と3つの解がありますので、ここでは掃き出し法(ガウスの消去法)を使って解いてみます。

9-3.連立方程式の解

 まず、上記の連立方程式を拡大係数行列の形で表しますと、

\left( \begin{array}{ccc|c}
A & -z_2 & 0 & F \\
  -z_2 & B & -z_4 & 0 \\
0 & -z_4 & C & 0 
\end{array} \right)

のようになります。この行列に基本変形を繰り返して左側が単位行列になるようにします。
 基本変形とは、下記の3つの変形をいいます。
 ①各行または各列を定数倍する
 ②各行または各列を入れ替える
 ③各行または各列を定数倍したものを各行または各列に加える
 まず、1行目をAで割ります。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A} \\
 -z_2 & B & -z_4 & 0 \\
0 & -z_4 & C & 0 
\end{array} \right)

 次に、1行目にz_2を乗じたものを2行目に足します。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A}\\
0 & B-\frac{z_2^2}{A} & -z_4 & \frac{z_2}{A}F \\
0 & -z_4 & C & 0 
\end{array} \right)

 計算がややこしくなるので、B-\frac{z_2^2}{A}=Dとしましょう。
 次に、2行目をDで割ります。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A}\\
0 & 1 & -\frac{z_4}{D} & \frac{z_2}{AD}F \\
0 & -z_4 & C & 0 
\end{array} \right)

 2行目にz_4を乗じたものを3行目に足します。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A}\\
0 & 1 & -\frac{z_4}{D} & \frac{z_2}{AD}F \\
0 & 0 & C-\frac{z_4^2}{D} & \frac{z_2z_4}{AD}F
\end{array} \right)

 C-\frac{z_4^2}{D}=Eとします。3行目をEで割ります。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A} \\
0 & 1 & -\frac{z_4}{D} & \frac{z_2}{AD}F \\
0 & 0 & 1& \frac{z_2z_4}{ADE}F
\end{array} \right)

 これで左下半分が0になりました。すでにV_3が求まっているのが分かります。この手順を前進消去といいます。
 次に、3行目に\frac{z_4}{D}を乗じたものを2行目に足します。
\left( \begin{array}{ccc|c}
1 & -\frac{z_2}{A} & 0 & \frac{F}{A} \\
0 & 1 & 0 & \frac{z_2(DE+z_4^2)}{AD^2E}F \\
0 & 0 & 1& \frac{z_2z_4}{ADE}F
\end{array} \right)

 2行目に\frac{z_2}{A}を乗じたものを1行目に足します。
\left( \begin{array}{ccc|c}
1 & 0 & 0 & \frac{DE(AD+z_2^2)+z_2^2z_4^2}{A^2D^2E}F\\
0 & 1 & 0 & \frac{z_2(DE+z_4^2)}{AD^2E}F \\
0 & 0 & 1& \frac{z_2z_4}{ADE}F
\end{array} \right)

 これで右上半分も0になりました。この手順を後退代入といいます。左側が単位行列になりましたので連立方程式の解は、 
\begin{cases} V_1 &= \frac{DE(AD+z_2^2)+z_2^2z_4^2}{A^2D^2E}F \\ V_2 &= \frac{z_2(DE+z_4^2)}{AD^2E}F \\ V_3 &= \frac{z_2z_4}{ADE}F \end{cases}

となります。

9-4.ダブルバスレフの音圧

 ダブルバスレフ型エンクロージャから放射される音は、ドライバーユニットと第2ポートから放射される音を合わせたものになります。位相反転型エンクロージャの説明にありましたように、ドライバーユニットの粒子速度とポートの粒子速度は位相が反転しているので、トータルの音圧は、ドライバーユニットの体積速度からポートの体積速度を引くのでした。
 したがって、ダブルバスレフのトータルの体積速度は、VS=S_dV_1-S_{p2}V_3となり、音圧は、

 p= j \omega \rho_0 \dfrac {S_dV_1-S_{p2}V_3 }{2 \pi x} e^{-jkx}

となります。

9-5.音圧をEXCELでシミュレーションする

 前回までのファイルに続いて入力していきます。追加する定数は下図の通りです。正直適当ですが、第2ポートの共振周波数の方が第1ポートの共振周波数よりも低くなるように調整します。


 ちなみに、第1容積・第1ポートと第2容積・第2ポートの区別がつくように、以前入力したキャビネットとポートの符号を変えておくと分かりやすいと思います。また、前回の位相反転型エンクロージャのシミュレーションのときから少しだけ数値を変更しました。

 空洞のスチフネスは\frac{\gamma P_0}{V_c}S_d^2で表すことができました。第2キャビネット(第2容積)のスチフネスs_{c2}は、振動板面積S_dに当たる部分が、第1ポートの開口面積S_{p1}に代わりますので、
s_{c2}=\frac{\gamma P_0}{V_{c2}}S_{p1}^2

となります。
 それでは、各列にインピーダンスを入力していきます。ちなみにz_1z_2は前回入力したAC列とAD列がそのまま使えますので、z_3以降を入力します。
 まず、列BUに第1ポートの質量m_{p1}、列BVに第2キャビネットのスチフネスs_{c2}、列BWに第2ポートの質量m_{p2}と放射質量m_{pa}の和を入力します。
 セルBU2には、
=IMPRODUCT(I2,$B$19)

 セルBV2には、
=IMDIV($B$33,I2)

 セルBW2には、
=IMPRODUCT(I2,$B$38+$B$40)

と入力してオートフィルします。こんな感じです。

 そして、列BX~BZ列にz_3z_4z_5を入力します。
 セルBX2には、
=IMPRODUCT( $B$23,IMSUM( $B$20,BU2))

 セルBY2には、
=IMPRODUCT( $B$23,IMSUM( BV2,$B$34))

 セルBZ2には、
=IMPRODUCT( $B$23,$B$42,IMSUM( BW2,( $B$39+$B$41)))

と入力してオートフィルします。こんな感じです。

 各インピーダンスには変圧比が掛かることに注意します。z_3z_4には、振動板と第1ポートの変圧比{(\frac{S_d}{S_{p1}})}^2が掛かり、z_5には振動板と第1ポートの変圧比{(\frac{S_d}{S_{p1}})}^2にさらに第1ポートと第2ポートの変圧比{(\frac{S_{p1}}{S_{p2}})}^2が掛かります。
 計算の簡略化のために置き換えた部分についてもあらかじめ入力しておきます。
 セルCA2には、Aとして置き換えた
=IMSUM(AC2,AD2)

 セルCB2には、Bとして置き換えた
=IMSUM(AD2,BX2,BY2)

 セルCC2には、Cとして置き換えた
=IMSUM(BY2,BZ2)

 セルCD2には、Dとして置き換えた
=IMSUB( CB2,IMDIV( IMPOWER( AD2,2),CA2))

 セルCE2には、Eとして置き換えた
=IMSUB( CC2,IMDIV( IMPOWER( BY2,2),CD2))

を入力してオートフィルします。こんな感じです。

 V_1V_3の計算式は、直接入力してもいいのですが、入力式が複雑で間違いやすいので、列CF~列CJに中間の計算をしておくと計算ミスを防止できます。

 そして、セルCK2には、V1として
=IMPRODUCT( IMDIV( IMSUM( CF2,CG2),CH2),J3)

 セルCL2には、V2として
=IMPRODUCT( IMDIV( IMPRODUCT( AD2,CI2),CJ2),J3)

 セルCM2には、V3として
=IMPRODUCT( IMDIV( IMPRODUCT( AD2,BY2),IMPRODUCT( CA2,CD2,CE2)),J3)

と入力してオートフィルします。こんな感じです。

 あとはこれまでと同じ要領で、列CN~列CPには振動板の音圧、デシベル値、位相を入力し、列CQ~列CSには第2ポートの音圧、デシベル値、位相を入力し、列CT~列CVにはトータルの音圧、デシベル値、位相を入力します。群遅延は列CWにトータルのものだけ求めておきましょう。
 こんな感じになっていればOKです。

 各グラフはこのようになります。



 ダブルバスレフは第1ポートの共振周波数の少し上に大きなディップができるのが特徴です。そのため、理論上はトリプルバスレフにしてディップの無い特性にした方が良いのですが、実際問題トリプルバスレフにすると大きすぎるし、機械抵抗が増大してポートの共振が弱まるので、実際にはシングルで十分な気がします。

9-6.ダブルバスレフのインピーダンス特性を求める。

 ダブルバスレフ型エンクロージャの場合は、下記の式(30)のz+z_0に、各機械インピーダンスの合成インピーダンスを代入します。

  Z =\dfrac{(BL)^2}{z+z_0}+(Z+Z_0)  \tag{30}

 ダブルバスレフ型エンクロージャの等価回路より、各機械インピーダンスの合成インピーダンスは、
 z_t = z_1+\dfrac{\alpha z_2}{\alpha +z_2}

となります。ただし、
 \alpha = z_3 + \dfrac{z_4z_5}{z_4+z_5}

としてまとめています。
 このz_tを式(30)のz+z_0と置き換えて、
 Z =\dfrac{(BL)^2}{z_t}+R_c+j \omega L_c

となります。
 EXCELで計算してみましょう。
 列CXに\alphaを計算しておきます。
 セルCX2に、
=IMSUM( BX2,IMDIV( IMPRODUCT( BY2,BZ2),IMSUM( BY2,BZ2)))

と入力してオートフィルします。次に、列CYにz_tを計算します。
 セルCY2に、
=IMSUM( BK2,IMDIV( IMPRODUCT( CX2,BL2),IMSUM( CX2,BL2)))

と入力してオートフィルします。そして、列CZにZを計算します。
 セルCZ2に、
=IMSUM( IMDIV( ( $B$2*$B$3)^2,CY2),$B$5,IMPRODUCT( I2,$B$6))

と入力してオートフィルします。列DAには絶対値を計算します。こんな感じです。

 グラフはこのようになります。

 ダブルバスレフの場合には、ピークが3つ立ちます。ピーク間のディップが各ポートの共振周波数です。

 ということで番外編としてダブルバスレフ型エンクロージャの特性をシミュレーションしていました。これを応用すればトリプルバスレフなんかも計算できるのですが、連立方程式を解くのが鬼のようにめんどくさい割に、再三になりますが、実際に作ると極めて大きくなりますので、現実的ではないと思います。

前へ 8.インピーダンス特性を求める
次へ 10.シミュレーションに使う定数の求め方