非線形音響研究会資料(11-5)
音響流駆動力と音響放射圧による超音波空気ポンプ内部の有限要素解析
和田有司，小山大介，中村健太郎（東工大・精研）Abstract : An ultrasonic air pump utilizing acoustic streaming is discussed and its efficient simulation method using finite element analysis (FEA) is suggested in this report. The pump induces air streaming in a thin air layer between a bending transducer and a reflector by exciting an intense sound field. The absolute sound pressure and streaming velocity in the air layer are calculated through following three steps: 1. linear acoustic analysis, 2. calculation of acoustics streaming driving force and radiation pressure, and 3. static incompressible flow analysis. The calculation results are compared with experimental results for the sound pressure and flow velocity. As a result, the sound pressure is calculated in the derivation of the error of 8.5%, and the measurement and analysis results for the streaming distribution are in good agreement.
keywords : acoustic streaming, air pump, bending vibration, driving force, FEA, radiation pressure
1
はじめに
近年，携帯端末などの小型薄型化及び高出力化 が進み，従来のファンやポンプでは不可能な狭空間 への各種気体供給の需要が存在する。我々は，超音 波領域で駆動するたわみ振動板と反射板にはさま れた薄い空気層に強い音場を励振する事で音響流 を発生させる空気ポンプを提案し検討しているが [1, 2, 3, 4]，これには機械的可動部がなく，小型化・ 薄型化が可能である。この超音波空気ポンプの最適 化には，装置内部の適切な音場および音響流の適切 なシミュレーションが必要である。 音響流の計算手法については古くから多くの検討 がなされており，二次微小項を考慮することで音響 流駆動力を算出し，解析的に流れを解く方法として， W. L. Nyborgの解説[5, 6]，およびO. V. Rudenko, S. I. Soluyanの著書[7]がよく知られている。しか しながら，解析的なアプローチでは極めて理想的な 条件のもとでしか流体方程式を解くことができず， 三次元流れを検討することは不可能である。 一方，近年の計算機性能の向上により圧縮性流 体方程式を直接解くことによって，音響流の計算を 行った報告もされている[4, 8]。とりわけ，共鳴管 内の音響流については多くの数値計算の報告がなさ れている[9, 10, 11, 12, 13, 14]。しかしながら，超 音波の周期に対して音響流の流れの発生や遷移にか かる時間のオーダーの差により，計算領域を二次元 に制限したり，粗いメッシュで計算を行ったりして もまだまだ計算負荷は大きく，実際の三次元デバイ スへの適用は困難なものとなっている。 本報告では有限要素解析(FEA)を用いて，Nyborg の音響流駆動力と音響放射圧を考慮することで音響 解析と流体解析を分離し，低計算負荷でポンプ能力 を判断する上で必要な音響流を計算する方法を提案 し，計算結果と測定結果を比較することで計算方法 の検証を行った。2
音響流と駆動力
Nyborg[5]は，以下のようにして音響流の導出が 可能であることを示した。はじめに圧力p，密度ρ， 流速uをそれぞれ以下のように展開する。 p− p0 = p1+ p2+ · · · , hp − p0i = p2 (1) ρ − ρ0 = ρ1+ ρ2+ · · · , hρ − ρ0i = ρ2 (2) u = u1+ u2+ · · · , hui = u2 (3) ここでh· · ·iは音波の周期より十分大きく，流れの 遷移時間より小さい間で時間平均とる操作を示し， p0, ρ0は定数である。添字 0は定数(大気圧 p0 = 1.0 × 105Pa，空気の密度ρ0= 1.3 kg/m3)，1は時間 的に振動する”音”の成分であり，この時間平均は零 となる。二次高調波もこの項に含まれる。2は時間 的には直流であるが，空間的に分布を持つ定常成分 であり，二次の微小量である。 圧縮性流体方程式に代入，3次以下の微小項を無 視すると，音響流U= u2+ hρ1u1i/ρ0は非圧縮性流 体方程式として扱うことができる。 ∂U ∂t + (U · ∇)U − η ρ0∇ 2U = −∇p2 ρ0 + F (4) ここでηは空気の粘性定数である。上式右辺に出現Fig. 1: Block diagram for the analysis.
するFが音響流駆動力であり，
F= −h(u1· ∇)u1i − hu1div u1i (5)
と表される。しかしながら，同様に右辺に現れるp2 の適切な値は文献[5]の時点では検討されておらず， 減衰進行波から発生する直進流を解く上では，人為 的に定常圧に勾配をもたせない限りこの値はほぼ零 とみなして差し支えなかった。 その後，Nyborg[15]は音波による小物体の浮揚 について，音響放射圧p2が，音のポテンシャルエ ネルギーと運動エネルギーの差から以下のように求 めることができるとしている。 p2= hepi − heki = |p1| 2 4ρ0c2 − ρ0|u1|2 4 (6) 以上のことから，音場からの流体に対する力を考 える場合，これら二つの影響を足しあわせた総合的 な力を考える必要がある。この総合的な力をFe f f とし，音響流が定常流∂U/∂t = 0であるとすると， 音響流の算出に必要な流体方程式は式(7, 8)のよう に表すことができる。 (U· ∇)U − η ρ0∇ 2U= F e f f (7)
Fe f f = −h(u1·∇)u1i−hu1div u1i−
∇hepi ρ0 + ∇heki ρ0 (8)
3
計算手順
Fig.1に計算手順のブロック図を示す。初めに，有 限要素法を用いて圧電-構造-音響連成解析を行い， 装置全体の粘性減衰を含む線形音場のシミュレー ションを行う。シミュレーション結果から，有限要 素格子の各節点上の音圧p1および粒子速度u1の結 果を複素数で得る。 次に，粒子速度u1，ポテンシャルエネルギー，お よび運動エネルギーについて有限要素格子上で数値Fig. 2: Basic structure of the device.
Fig. 3: Bending vibration mode on the transducer at the resonant frequency of 26.2 kHz.
微分を行う。微分の結果から，構成する各要素を計 算し足し合わせることで，各節点におけるFe f f を 算出する。 最後に，音場解析に用いたものと同一の格子上 で，有限要素法による非圧縮定常流体解析を行う。 先程求めた，Fe f f を体積力として各節点に入力す ることで流体を駆動する。シミュレーションにより 得られた流れ場が音響流に相当する。
4
装置の構成
Fig. 2に超音波空気ポンプの形状を示す。PZTの 矩形板(長さ10 mm，幅 30 mm，厚さ 0.4 mm)を 接着した長さ20 mm，幅30 mmのアルミニウム製 振動板に対し，振動板と同形状のアクリル製反射板 を音波の波長より十分小さい空気層を介して配置す る。振動板右側のギャップを1 mmに固定して反射 板を4度傾けることで，進行波音場をつくりだし一 方向の音響流を誘起する。なお，座標系について振 動板中央から高さ0.5 mmの空気層中に原点にとり， 長さ方向にx軸，幅方向にy軸，厚さ方向にz軸を とる。今回は，計算結果の検討のために進行波が発 生しにくい0度の場合についても結果を示す。 Fig. 3に示すように，振動板に長さ方向(x軸方 向)に沿って1次たわみ振動を励振する。強力な音Table 1: Material properties used in calculations. Air
Sound speed [m/s] 340 Density [kg/m3] 1.29 Viscosity coefficient [×10−6Pa sec] 1.62
PZT
Young’s modulus [GPa] 74.7 Density [kg/m3] 7700
Poisson rate 0.290 Electromechanical coupling factor 0.710 Relative permittivity 1450 Q-farctor 200 Aluminum
Young’s modulus [GPa] 70.3 Density [kg/m3] 2698 Poisson rate 0.350 Q-farctor 200
Table 2: Specifications of Computer, FEA, and the re-sources spent for the calculations.
CPU AMD Phenom 9550 Quad-Core Processor 2.20 GHz Memory 8.0 Gbyte Operation System CentOS Linux 5.4
Nodes 25,662
Elements 23,142 Shared-Memory Parallel 2 process Memory used 510.8 MByte Calculation time for sound 357.4 sec Calculation time for fluid 198.3 sec
場を形成するために，振動板振動分布と空気層音圧 分布の相関が大きい長さ方向(x方向)に腹が3つ， 幅方向(y方向)に腹が1つの(nx, ny)= (3, 1)モード の音場を励振する。(3,1)モードの場合 f = 26.2 kHz であり，この周波数で1次たわみ振動するよう振動 板の厚みを2 mmとした。測定，計算ともに振動板 の最大振動速度振幅が0.5 m/sとなるように印加電 圧の調整を行った。
5
計算条件
(
有限要素モデル
)
Fig.4にそれぞれ音場と流体の有限要素モデルを 示す。装置全体を32x38x15 mm3の直方体空気塊で 覆い，音場および流れ場の計算を行う。音場と流 体では境界条件の種類が異なるが有限要素格子自Fig. 4: FEA model.
体は同一である。有限要素法解析にはANSYS11.0
（ANSYS Inc.)を用いた。Table 1, 2にはそれぞれ解 析に用いた材料定数，および計算機のスペックと計 算時間を示す。 計算格子のメッシュサイズは1 mm，反射板と振 動板に挟まれる空気ギャップの厚み方向に対しては 0.5 mmとした。本来ならば，定在波により発生す る厚み方向渦の考慮のために，特に厚み方向に対し てはより細かいメッシュサイズを採用することが望 ましい。しかしながら，ANSYSライセンス上の要 素数制限や開領域の表現に必要な領域の確保の制約 に対して，ポンプの吸込・吐出能力に対してこの循 環渦の重要性はそこまで大きくないため，ポンプ能 力の評価のために1 mmのメッシュサイズで差し支 えないと判断した。 音場解析では周波数応答解析を用いて，PZT両 端に電圧境界条件を設定することにより駆動する。 空気塊外部には音波の吸収境界条件を，空気と振動 子・反射板との界面には粘性減衰[16]を設定した。 流体解析では，音場解析の結果から計算された駆 動力Fe f f を各接点に体積力として入力，非圧縮定 常流体解析を行った。空気塊外部には流れの吸収条 件として圧力が零となる条件を設定し，空気と振動 子・反射板との界面には粘着条件を設定した。
6
測定手法
6.1 音圧分布 空気層中の音場分布の測定には，光ファイバプ ローブを用いた[17]。Fig.5に測定系を示す。音圧 変化により生じる空気の屈折率変化をファイバ端Circulator Ref. Sig. Function Generator Transducer ASE light source PD
Lock−in Amp Amp.
Mesurement object
Single−mode fiber
Fig. 5: Measurement setup for sound pressure.
Fig. 6: Measurement setupt for PIV.
面における反射率変調として検出し，音圧の絶対 値を測定を算出する。測定は−15 < x < 15 mm， −15 < y < 15 mmの領域を1 mmごと測定した。 6.2 流速分布 流れ場の測定には粒子画像流速測定法(PIV)を用 いた。Fig.6に測定系を示す。トレーサーして石松 子を振動板上に散布したうえで装置を駆動して空 気流れを起こし，LED照明のもと高速度カメラで 石松子の動きを撮影した。振動板の振動の影響の 写り込みを軽減するため，撮影周期は駆動周波数の ちょうど1/40となる655 fpsとした。露光時間は10 µs，−13.7 < x < 20.1 mm, −9.6 < y < 11.4の領域 を512×328ピクセルの画像として1500フレーム分 (1.16秒)撮影行った。 PIV の設定としては，画像を 20×20 ピクセル (1.50×1.50 mm)の格子に切りその領域のパターンを トレースし，2画像間のパターンの変位と撮影周期 から流速を算出し，ノイズの影響を軽減するため全 フレームでの平均を取ることで，流速を算出した。
Fig. 7: Sound pressure distribution.
7
測定結果と計算結果の比較
7.1 音圧分布 Fig.7に，音圧振幅分布の計算結果(上段)，x軸上 における音圧振幅(中段)・位相(下段)の計算・測定 結果を示す。反射板の傾きが0度の時は振幅1500 Pa，4度の時は振幅800 Pa程度の(3,1)モードの音 場が発生していることが確認できる。位相分布から 判断して傾き0度の時はほとんど完全な定在波であ るのに対して，傾き4度の時では左から右への進行 波成分が発生しているといえる。 測定結果と計算結果の誤差の標準偏差はそれぞれ 140.2 Pa(8.7%)，68.5 Pa(8.6%)であった。誤差の原 因としては測定に用いた光ファイバプローブに起因 する誤差が100 Pa程度生じる事を考慮にいれれば， 空気ポンプの設計に必要な計算精度としては十分で あると考えられる。 7.2 駆動力分布 Fig.8 に計算結果から算出した，音響流駆動力 F(上 段)，音 響 流 駆 動 力 F，音 響 放 射 圧 の 勾 配 -gradp2/ρ0，これらの和である入力Fe f f のx方向成 分のx軸上における分布(中段)，Fe f f の分布(下段) を示す。傾き0度・4度の場合共に，Fの分布ではFig. 8: Driving force distribution. ポンプ内部において多くの対向するベクトルが見ら れるのに対して，図の中段で示されるようにその成 分の多くが音響放射圧の勾配と相殺され，下段に示 されるように対向するベクトルの多くが失われる。 傾き0度と4度の違いとしては，x = −5 mm付 近に見られる左向きベクトルの存在で，これが左向 き流れを誘引する主な原因となっている。 7.3 流れ分布 Fig.9，10にそれぞれ傾き0度と4度の時の流れ 分布の計算結果とPIVによる測定結果を示す。 傾き0度の時は流体を駆動する力がほぼ対称であ るため，計算結果では力が拮抗し両方向の流れが， PIV測定結果では，計算では考慮していない循環渦 に起因する乱雑なベクトル分布が見られる。いずれ の結果よりも，傾き0度の場合はポンプ能力におい ては好ましくない結果が示されていると言える。 これに対して，傾き4度の時は，計算・測定結果 共に左向き流れを示すベクトルが多く見られ，計算 結果と測定結果はよく一致していると言える。
Fig. 9: Flow velocity distribution when the gradient of the reflector is 0 degree.
Fig. 10: Flow velocity distribution when the gradient of the reflector is 4 degrees.
8
まとめ
本報告では，有限要素法を用いた超音波空気ポン プ内部における音場および流れ場の計算手法の提案 を行い，測定結果との比較を行った。その結果，音 圧分布の計算結果は測定結果と誤差5.8%程度でよ く一致，流れ場の結果ではポンプ能力を議論するう えで差し支えない程度には，流れの傾向が一致して いる結果が得られた。参考文献
[1] H. Takei, D. Koyama, K. Nakamura, and S. Ueha: Jpn. J. Appl. Phys. 47 (2008) 4276.
[2] H. Takei, D. Koyama, K. Nakamura, and S. Ueha: Denshi Joho Tsushin Gakkai Ronbunshi A 91 (2008) 1152 [in Japanese].
[3] D. Koyama, Y. Wada, K. Nakamura, M. Nishikawa, T. Nakagawa and H. Kihara: IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 57 (2010) 253.
[4] Y. Wada, D. Koyama, K. Nakamura: Jpn. J. Appl. Phys. 49 07HE15 (2010).
[5] Wesley L. Nyborg : J. Acoust. Soc. Am., 25 (1953) pp. 68-75.
[6] Wesley L. Nyborg : J. Acoust. Soc. Am., 30 (1958) pp. 329-339.
[7] O. V. Rudenko and S. I. Soluyan ”THEO-RETICAL FOUNDATIONS OF NONLINEAR ACOUSTICS” (Consultants Bureau, New York, 1977)
[8] L. P. Cheng and S. Y. Zhang : Appl. Phys. Lett. 90, 244106 (2007).
[9] M. Kawahashi and M. Arakawa: JSME Int. J. 39B (1996) pp. 280-286.
[10] T. Yano: J. Acoust. Soc. Am. 106 (1999) L7-L12. [11] A. Alexeev and C. Gutfinger: Phys. Fluids 15
(2003) pp. 3397-3408.
[12] M. K. Aktas and B. Farouk: J. Acoust. Soc. Am. 116(2004) pp. 2822-2831.
[13] J. P. Boris, A. M. Landsberg, E. S. Oran, and J. H. Gardner, ”LCPFCT―A flux-corrected trans-port algorithm for solving generalized continuity equations,” Report No. NRL/MR/6410-93-7192, (1993).
[14] 矢野 猛:ながれ24(2005) pp. 371-380.
[15] Wesley L. Nyborg : J. Acoust. Soc. Am., 42 (1967) pp. 947-952.
[16] 早坂寿雄, 吉川昭吉郎著 ”音響振動論” 丸善
(1974)
[17] H. Takei, T. Hasegawa, K. Nakamura, and S. Ueha: Jpn. J. Appl. Phys. 46 (2007) 4555.
