[[授業]]



*第1回(2021.8.24) [#m84551e8]

- 準備
-- 参考資料:&ref(3年学生実験2018/B2.pdf); の 3.1節、3.2節、付録A, B
-- 複素振幅による調和振動の表現 ... 3.1節、式(2)
-- 圧力と粒子(体積)速度に関する基本の式 ... 付録A、式(27)と(31)
- 目的:電力フィードバック型熱音響発電機の安定性解析用mファイルを理解する
&ref(make_graph2_fixed_add_label_moved.m);
-- [[篠田他2018>https://www.jstage.jst.go.jp/article/jasj/74/6/74_313/_pdf/-char/ja]]
-- [[小林、梅本他2017>https://www.jstage.jst.go.jp/article/jasj/73/1/73_12/_pdf/-char/ja]]
-- dl_core, dl_pipe, l_sens
-- コア部の周波数応答の算出 &ref(load_data_Gstack_AB_fixed.m); ... [小林、梅本他](2)式
   tmp = [A2_spk1_g.resp(k), A2_spk2_g.resp(k); B1_spk1_g.resp(k), B1_spk2_g.resp(k)] / [A1_spk1_g.resp(k), A1_spk2_g.resp(k); B2_spk1_g.resp(k), B2_spk2_g.resp(k)];
--- &ref(load_data_raw_AB_fixed.m); ... [小林、梅本他](3)式
 A1_jw = ( exp(-i*k*dl).*ps1_jw - exp(-i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 B1_jw = (-exp( i*k*dl).*ps1_jw + exp( i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 A2_jw = (-exp( i*k*dl).*ps2_jw + exp( i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));
 B2_jw = ( exp(-i*k*dl).*ps2_jw - exp(-i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));
-- リニアモータ部の周波数応答の算出 &ref(load_data_H.m);
-- その他、管路部・回路部の周波数応答の算出 
-- ナイキストの安定判別に基づく発振周波数などの解析 &ref(nyquist_analysis_spline_frange_fixed.m);

#ref(2021.08.24-1.jpg,left,noimg,板書1)
#ref(2021.08.24-2.jpg,left,noimg,板書2)
#ref(2021.08.24-3.jpg,left,noimg,板書3)
#ref(2021.08.24-4.jpg,left,noimg,板書4)


*第2回(2021.8.31) [#y90d142b]
- 前回のつづき
-- コア部の周波数応答の算出
--- &ref(calibrate_fixed.m); ... [小林、梅本他2017]3.1節、(7)式、[篠田他2018](4)式
 function [tmp_g] = calibrate_fixed(tmp_g, G1_g, G2_g, G1m_g, G2m_g);
 tmp_g(1,1) = tmp_g(1,1) * G1m_g(1,1) / G1_g(1,1);
 tmp_g(2,2) = tmp_g(2,2) * G1m_g(2,2) / G1_g(2,2);
 tmp_g(1,2) = tmp_g(1,2) * G2m_g(1,2) / G2_g(1,2);
 tmp_g(2,1) = tmp_g(2,1) * G2m_g(2,1) / G2_g(2,1);
--- 5段接続
 Gcore1_g = G_g; % コア1段分
 for m = 2:5
     G_g = lft(G_g, Gcore1_g, 1, 1); % m stage
 
     Ltube = Ltube - Lcore1; % リニアモータ間距離 - コア1段分の長さ
 end
-- リニアモータ部の周波数応答の算出
--- &ref(load_data_H.m); ... [小林、梅本他2017](12)式相当(入力 vs, i を入れ替えて定義しているので注意)
 tmp = [Ain_spk1_g.resp(k), Ain_spk2_g.resp(k); Bout_spk1_g.resp(k), Bout_spk2_g.resp(k)] /[vs_spk1_g.resp(k), vs_spk2_g.resp(k); i_spk1_g.resp(k), i_spk2_g.resp(k)];
--- &ref(load_data_raw_ABvi.m); ... 進行波圧力成分の複素振幅 A, B はコア部の左側 A1, B1 と同じ
 ps_jw = frdata(:,4).*exp(i*pi/180.0*frdata(:,5))*nagano_keiki_coeff;
 pc_jw = frdata(:,2).*exp(i*pi/180.0*frdata(:,3))*nagano_keiki_coeff;
 vi_jw = frdata(:,6).*exp(i*pi/180.0*frdata(:,7));
 vs_jw = frdata(:,8).*exp(i*pi/180.0*frdata(:,9));
 i_jw = (vs_jw - vi_jw)/R;
 (略)
 Ain_jw = ( exp(-i*k*dl).*ps_jw - exp(-i*k*(ls+dl)).*pc_jw)./( 2*i*sin(k*ls));
 Bout_jw = (-exp( i*k*dl).*ps_jw + exp( i*k*(ls+dl)).*pc_jw)./( 2*i*sin(k*ls));
--- 補正 ... [篠田他2018](10)式
 F_g = G2m_g(2,1) / G2_g(2,1); % 補正係数(2015.1.24 3頁参照)
 H_g = mdiag(1, F_g)*H_g; % 補正。
 H_g = H_g * [0, 1; 1, 0]; % スター積利用のため列を入れ替え(これでH_gの右から [i; v]がかかるようになる)
--- 左右のリニアモータの伝達マトリクス ... [篠田他2018](11)式
 H1_g = H_g; % 入力入れ替え済
 H2p_g = [0, 1; 1, 0]*H1_g*[1, 0; 0, -1]; % 梅本論文(20)式の下の式。
-- 管路部の周波数応答の算出 ... [篠田他2018](12)式
 tmp = frd(exp(-1i*k*Ltube), w); % (20)式
 Gtube_g = [0, tmp; tmp, 0];
-- 回路部の周波数応答の算出
 Z1 = Rinf;
 Z2 = Rinf;
 Z3 = R;
 Tcirc = [1 + Z3/Z2, -(Z1 + Z2 + Z3)/(Z1*Z2); -Z3, 1 + Z3/Z1]; % (19)式
 Gcirc = chain_inv(Tcirc); % その下の式


#ref(2021.08.31-1.jpg,left,noimg,板書1)
#ref(2021.08.31-2.jpg,left,noimg,板書2)

*第3回(2021.9.16) [#e30254ae]
- 目的(つづき):電力フィードバック型熱音響発電機の安定性解析用mファイルを理解する
-- 管路部の周波数応答の算出 ... [篠田他2018](12)式
 tmp = frd(exp(-1i*k*Ltube), w); % (20)式
 Gtube_g = [0, tmp; tmp, 0];
-- 回路部の周波数応答の算出 ... [篠田他2018](14)式
 Z1 = Rinf;
 Z2 = Rinf;
 Z3 = R;
 Tcirc = [1 + Z3/Z2, -(Z1 + Z2 + Z3)/(Z1*Z2); -Z3, 1 + Z3/Z1]; % (19)式
 Gcirc = chain_inv(Tcirc); % その下の式
-- コア部から見た管路の周波数応答の算出 ... [篠田他2018](15)式
-- コア部から見た管路の周波数応答の算出 ... [篠田他2018](15)式 &ref(chain_inv.m);
 Gother_g = lft(lft(chain_inv(H2p_g), Gcirc,  1, 1), lft(chain_inv(inv(H1_g)), Gtube_g, 1, 1), 1, 1); % (22)式
-- ナイキストの安定判別に基づく発振周波数などの解析 &ref(nyquist_analysis_spline_frange_fixed.m);
 result{num} = nyquist_analysis_spline_frange_fixed(G_g*[0, 1; 1, 0]*Gother_g*[0, 1; 1, 0], 20, fmin(num), fmax(num))
&br;
 function [s] = nyquist_analysis_spline_frange_fixed(L_g, m, fmin, fmax);
 phi0_g = (1 - L_g(1,1))*(1 - L_g(2,2)) - L_g(1,2)*L_g(2,1);
- 目的2:相関法、フーリエ級数におけるフーリエ係数の導出(2018年度第6回(2018.7.3)参照)

//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
&color(black,red){&size(20){!!!!!!!!!!!!!!!!!!!!!!!!!!!!!以下は過去の情報です!!!!!!!!!!!!!!!!!!!!!!!!!!};};
#ref(2021.09.16-1.jpg,left,noimg,板書1)
#ref(2021.09.16-2.jpg,left,noimg,板書2)
#ref(2021.09.16-3.jpg,left,noimg,板書3)
#ref(2021.09.16-4.jpg,left,noimg,板書4)

-- 一般の場合:Z(L) = α ρ0 c0
- 2センサ法による進行波圧力成分の計算(去年の第3回も参照) ... https://doi.org/10.5687/iscie.28.392 (22)式の下の式; &ref(load_data_raw_AB_fixed_rows_changed_PCB.m);
 pc2_jw = frdata(:,4).*exp(i*pi/180.0*frdata(:,5))/pcb_coeff_12638; % 
 ps1_jw = frdata(:,8).*exp(i*pi/180.0*frdata(:,9))/pcb_coeff_12992; % dummy
 ps2_jw = frdata(:,2).*exp(i*pi/180.0*frdata(:,3))/pcb_coeff_12637; % 
 pc1_jw = frdata(:,6).*exp(i*pi/180.0*frdata(:,7))/pcb_coeff_12991; % dummy
 
 A1_jw = ( exp(-i*k*dl).*ps1_jw - exp(-i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 B1_jw = (-exp( i*k*dl).*ps1_jw + exp( i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 A2_jw = (-exp( i*k*dl).*ps2_jw + exp( i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));
 B2_jw = ( exp(-i*k*dl).*ps2_jw - exp(-i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));

#ref(2019.06.27-1.jpg,left,noimg,板書1)
#ref(2019.06.27-2.jpg,left,noimg,板書2)
#ref(2019.06.27-3.jpg,left,noimg,板書3)
#ref(2019.06.27-4.jpg,left,noimg,板書4)

*第3回(2019.7.4) [#cd6371f6]
*第4回(2021.9.27) [#rae588c5]

- 伝達マトリクスと因果的システム表現
- 例:長さlで断面が一様な管
- 伝達マトリクスに対応する伝達関数行列は非因果的
- 因果的システム表現(scattering matrix)
- 目的2:相関法、フーリエ級数におけるフーリエ係数の導出(2018年度第6回(2018.7.3)参照)
- 伝達マトリクスと因果的システム表現 ... [[小林、山田2017>https://www.jstage.jst.go.jp/article/iscie/28/9/28_392/_pdf/-char/ja]]
-- 例:長さlで断面が一様な管
-- 伝達マトリクスに対応する伝達関数行列は非因果的
-- 因果的システム表現(scattering matrix)

#ref(2021.09.27-1.jpg,left,noimg,板書1)
#ref(2021.09.27-2.jpg,left,noimg,板書2)
#ref(2021.09.27-3.jpg,left,noimg,板書3)
#ref(2021.09.27-4.jpg,left,noimg,板書4)


//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
&color(black,red){&size(20){!!!!!!!!!!!!!!!!!!!!!!!!!!!!!以下は過去の情報です!!!!!!!!!!!!!!!!!!!!!!!!!!};};

- ナイキストの安定判別とシミュレーション
-- Matlabソース &ref(simu.m);
-- Simulinkモデル &ref(simple.mdl);
-- その他、必要なもの ... nyquistplot.m, setoptions.m
//&ref(nyquistplot.m); &ref(setoptions.m);

#ref(2019.07.04-1.jpg,left,noimg,板書1)
#ref(2019.07.04-2.jpg,left,noimg,板書2)
#ref(2019.07.04-3.jpg,left,noimg,板書3)
#ref(2019.07.04-4.jpg,left,noimg,板書4)

*第4回(2019.7.18) [#s24b85d5]

- 行列表現によるナイキスト安定判別:因果的システム表現を用いる場合と、伝達マトリクスを用いる場合
- 例題(Matlab)
-- 定在波管の場合(前回の例と同一)&ref(simu_add.m);
--- 因果的システム表現による解析結果は、前回と同一になる
--- 伝達マトリクスは、分割されたシステムに対して定義できない(解析不可)
-- ループ管の場合 &ref(simu_loop.m); &ref(chain.m);
--- 伝達マトリクスを因果的システムと誤って使用すると、意味のないナイキスト軌跡が得られる
--- 因果的システム表現による解析結果 ... 発振周波数が単純計算した共振周波数に合わないかも、■■■要検討■■■

#ref(2019.07.18-1.jpg,left,noimg,板書1)
#ref(2019.07.18-2.jpg,left,noimg,板書2)
#ref(2019.07.18-3.jpg,left,noimg,板書3)

*第5回(2019.7.25) [#zf8bee59]

- 第4回の補足:ナイキストの安定判別が単純計算したループ管の共振周波数に合わない理由 &ref(simu_loop_fixed.m);
- 第3回の補足:「反射係数」が共振周波数に与える効果(複素数の場合)...電力フィードバックの共振周波数のヒント &ref(simu_comp_r.m);
- 段付の管を含む共振周波数の計算 &ref(check_freqresp_calibration.m); &ref(freqresp_calibration.zip);

#ref(2019.07.25-1.jpg,left,noimg,板書1)
#ref(2019.07.25-2.jpg,left,noimg,板書2)



- 前回のつづき
-- 粒子(体積)速度に関する基本の式(31)
-- 特性音響インピーダンス
- 周波数応答実験結果とモデルの比較
-- 周波数応答実験の結果(3年学生実験の装置): &ref(spk1.dat); &ref(spk2.dat);
-- mファイル: &ref(freqresp_distx_p.m); &ref(p_tilde_down.m); &ref(p_tilde_up.m); 

#ref(2018.04.27-1.jpg,left,noimg,板書1)
#ref(2018.04.27-2.jpg,left,noimg,板書2)

*第3回(2018.5.9) [#g7d6b9f3]

- 管内の圧力分布(前回の続き):&ref(freqresp_distx_p_fixed.m);
-- 周波数が高くなると、(圧力の)腹と節の数が増える
-- マイクの位置が節に重なる周波数で、ゲイン特性が下向きのピークを示す
-- wからz,yまでの周波数応答における違い:リファレンスマイクは何度も節を横切るが、エラーセンサは周波数が十分高くなって初めて節を横切る
-- ...
- むだ時間要素(導波管、2センサ法、指向性センサ、etc. の基本)の伝達関数と周波数応答
- 2センサ法による進行波圧力成分の計算 &ref(load_data_raw_AB_fixed_rows_changed_PCB.m);
 pc2_jw = frdata(:,4).*exp(i*pi/180.0*frdata(:,5))/pcb_coeff_12638; % 
 ps1_jw = frdata(:,8).*exp(i*pi/180.0*frdata(:,9))/pcb_coeff_12992; % dummy
 ps2_jw = frdata(:,2).*exp(i*pi/180.0*frdata(:,3))/pcb_coeff_12637; % 
 pc1_jw = frdata(:,6).*exp(i*pi/180.0*frdata(:,7))/pcb_coeff_12991; % dummy
 
 A1_jw = ( exp(-i*k*dl).*ps1_jw - exp(-i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 B1_jw = (-exp( i*k*dl).*ps1_jw + exp( i*k*(ls+dl)).*pc1_jw)./( 2*i*sin(k*ls));
 A2_jw = (-exp( i*k*dl).*ps2_jw + exp( i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));
 B2_jw = ( exp(-i*k*dl).*ps2_jw - exp(-i*k*(ls+dl)).*pc2_jw)./( 2*i*sin(k*ls));

#ref(2018.05.09-1.jpg,left,noimg,板書1)
#ref(2018.05.09-2.jpg,left,noimg,板書2)

*第4回(2018.5.14) [#c812f624]

- 前回のつづき:2センサ法による進行波圧力成分の計算 ... https://doi.org/10.5687/iscie.28.392 (22)式の下の式

#ref(2018.05.14-1.jpg,left,noimg,板書1)
#ref(2018.05.14-2.jpg,left,noimg,板書2)
#ref(2018.05.14-3.jpg,left,noimg,板書3)
#ref(2018.05.14-4.jpg,left,noimg,板書4)
#ref(2018.05.14-5.jpg,left,noimg,板書5)

*第5回(2018.%%5.31%%6.11) [#o9772802]

- 2センサ法の限界(特定の周波数では算出できない)
- 指向性センサの周波数帯域&br;
/home/baba/spkcont/&ref(spkcont_app.c); (参考:&ref(spkcont_module.c); &ref(spkcont.h);):
 #ifdef PCB
     As1 = sqrt(msg->as1*msg->as1 + msg->bs1*msg->bs1) / PCB_COEFF_12638;
     Ac1 = sqrt(msg->ac1*msg->ac1 + msg->bc1*msg->bc1) / PCB_COEFF_12637;
     As2 = sqrt(msg->as2*msg->as2 + msg->bs2*msg->bs2) / PCB_COEFF_12992;
     Ac2 = sqrt(msg->ac2*msg->ac2 + msg->bc2*msg->bc2) / PCB_COEFF_12991;
 #else // NAGANO_KEIKI
     As1 = sqrt(msg->as1*msg->as1 + msg->bs1*msg->bs1) * NAGANO_KEIKI_COEFF;
     Ac1 = sqrt(msg->ac1*msg->ac1 + msg->bc1*msg->bc1) * NAGANO_KEIKI_COEFF;
     As2 = sqrt(msg->as2*msg->as2 + msg->bs2*msg->bs2) * NAGANO_KEIKI_COEFF;
     Ac2 = sqrt(msg->ac2*msg->ac2 + msg->bc2*msg->bc2) * NAGANO_KEIKI_COEFF;
 #endif
                
 // Phase
 Ps1 = atan(msg->bs1/msg->as1) + M_PI * (msg->as1 < 0);
 Pc1 = atan(msg->bc1/msg->ac1) + M_PI * (msg->ac1 < 0);
 Ps2 = atan(msg->bs2/msg->as2) + M_PI * (msg->as2 < 0);
 Pc2 = atan(msg->bc2/msg->ac2) + M_PI * (msg->ac2 < 0);
 
 A1 = sqrt(As1*As1 - 2*As1*Ac1*cos(k*L + Ps1 - Pc1) + Ac1*Ac1) / (2*fabs(sin(k*L)));
 B1 = sqrt(As1*As1 - 2*As1*Ac1*cos(k*L + Pc1 - Ps1) + Ac1*Ac1) / (2*fabs(sin(k*L)));
 A2 = sqrt(Ac2*Ac2 - 2*Ac2*As2*cos(k*L + Pc2 - Ps2) + As2*As2) / (2*fabs(sin(k*L)));
 B2 = sqrt(Ac2*Ac2 - 2*Ac2*As2*cos(k*L + Ps2 - Pc2) + As2*As2) / (2*fabs(sin(k*L)));

#ref(2018.06.11-1.jpg,left,noimg,板書1)
#ref(2018.06.11-2.jpg,left,noimg,板書2)
#ref(2018.06.11-3.jpg,left,noimg,板書3)

*第6回(2018.7.3) [#j181ed3d]
- 前回のつづき
- /home/baba/spkcont/&ref(spkcont_module.c); の一部
 if(theta > 2*PERIODS*M_PI){
         theta -= 2*PERIODS*M_PI;
         msg->as1 = 2.* as1 / count;
         msg->bs1 = 2.* bs1 / count;
         msg->ac1 = 2.* ac1 / count;
         msg->bc1 = 2.* bc1 / count;
         msg->as2 = 2.* as2 / count;
         msg->bs2 = 2.* bs2 / count;
         msg->ac2 = 2.* ac2 / count;
         msg->bc2 = 2.* bc2 / count;
 
         as1 = 0;
         bs1 = 0;
         ac1 = 0;
         bc1 = 0;
         as2 = 0;
         bs2 = 0;
         ac2 = 0;
         bc2 = 0;
         msg->count = count;
         count = 0;
 }
 (略)
 ad_conv(&y);
 
 ps1 = y[1]; // / PCB_COEFF_12638;
 pc1 = y[2]; // / PCB_COEFF_12637;
 //              pc1 = y[0]; // / PCB_COEFF_12637;
 pc2 = y[2]; // / PCB_COEFF_12991;
 ps2 = y[3]; // / PCB_COEFF_12992;
 
 as1 += ps1 * cos(theta);
 bs1 -= ps1 * sin(theta);
 ac1 += pc1 * cos(theta);
 bc1 -= pc1 * sin(theta);
 as2 += ps2 * cos(theta);
 bs2 -= ps2 * sin(theta);
 ac2 += pc2 * cos(theta);
 bc2 -= pc2 * sin(theta);
-- 相関法、フーリエ級数におけるフーリエ係数の算出
-- A/D変換器の ch0, ch1, ch2, ch3 の電圧がそれぞれ y[0], y[1], y[2], y[3] に入る。ただし、馬場君の実験ではch0, ch3 に結線していないので、y[0], y[3] にはゴミが入る。

#ref(2018.07.03-1.jpg,left,noimg,板書1)
#ref(2018.07.03-2.jpg,left,noimg,板書2)
#ref(2018.07.03-3.jpg,left,noimg,板書3)

//■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
&color(black,red){&size(20){!!!!!!!!!!!!!!!!!!!!!!!!!!!!!以下は過去の情報です!!!!!!!!!!!!!!!!!!!!!!!!!!};};

- 長さlの管の伝達マトリクス(入出力関係が非因果的)
- 進行波圧力振幅を入出力とするシステム(因果的)



#ref(2015.04.30-1.jpg,left,noimg,板書1)
#ref(2015.04.30-2.jpg,left,noimg,板書2)
#ref(2015.04.30-3.jpg,left,noimg,板書3)
#ref(2015.04.30-4.jpg,left,noimg,板書4)

*第3回(2015.5.29) [#c963c084]

- ナイキストの安定判別(開ループ系が安定な場合)
- むだ時間系の安定判別
- MIMO版 ナイキストの安定判別
- ループ管の安定判別

#ref(2015.05.29-1.jpg,left,noimg,板書1)
#ref(2015.05.29-2.jpg,left,noimg,板書2)
#ref(2015.05.29-3.jpg,left,noimg,板書3)

*第4回(2015.8.26) [#fe774753]

- 2センサ法
- 指向性センサ

#ref(2015.08.26-1.jpg,left,noimg,板書1)
#ref(2015.08.26-2.jpg,left,noimg,板書2)
#ref(2015.08.26-3.jpg,left,noimg,板書3)
#ref(2015.08.26-4.jpg,left,noimg,板書4)
#ref(2015.08.26-5.jpg,left,noimg,板書5)

第5回(2015.8.27)

- リニアモータの電気-音響特性に基づく管路の計算手順
- CHAIN とその逆変換

#ref(2015.08.27-1.jpg,left,noimg,板書1)
#ref(2015.08.27-2.jpg,left,noimg,板書2)
#ref(2015.08.27-3.jpg,left,noimg,板書3)


トップ   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS