音響流の解析

By M.Sato

July 9, 2026

1. 概要

音響流（Acoustic Streaming）とは、強力な超音波が流体中を伝播する際、音波の減衰や境界との摩擦（レイノルズ応力）によって生じる定常的な流れの現象である。本来は振動にすぎない音波が、流体を駆動する効果を持つことが特徴である。

この現象は、マイクロ流体デバイス（Lab-on-a-Chip）内での細胞・粒子の非接触なミキシングや操作、高効率な局所冷却技術、超音波洗浄のメカニズム解明など、多岐にわたる分野で応用されている。

解析に際しては、圧縮性流体の非線形振動としてナビエ・ストークス方程式を解くことが必要であるが、効率的に解析を行うために、摂動法に基づき一次オーダーの音響場と二次オーダーの流れ場を分離して解析する手法などが用いられている。

ここでは、過去のOpenFOAMの解析例を参考にして、最新版のOpenFOAM(v2512)にソルバーを移植して解析を行った。なお今回実施したのは一次オーダーの音響場の計算である。二次オーダーの流れ場については今後機会があれば実施する予定である。

2. sonicLiquidFoam

OpenFOAMには、圧縮性の液体を対象としたsonicLiquidFoamというソルバーが存在する。今回はsonicLiquidFoamを改造して新しいソルバーを作成した。

sonicLiquidFoamでは、以下に示す連続の式と運動量方程式を解くことで流体の運動を計算している。

$$ \begin{align} & \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \notag \\ & \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot ( \mu \nabla \mathbf{u} ) \notag \end{align} $$

ここで、$\rho$は密度、$\mathbf{u}$は速度ベクトル、$p$は圧力、$\mu$は粘性係数である。

液体では、通常の理想気体の状態方程式（$p = \rho R T$）ではなく、液体の体積弾性率を考慮した状態方程式が使用される。

$$ \begin{align} \rho = \rho_0 + \psi (p-p_0) = \tilde{\rho} + \psi p \notag \end{align} $$

ここで$\psi$は流体の圧縮率であり、$\rho_0$は基準密度、$p_0$は基準圧力である。 $dp/d\rho=c^2$である（$c$は音速）ことから、$\psi=\frac{1}{c^2}$となる。なお、$\tilde{\rho}=\rho_0 - \psi p_0$である。

2.1. 計算アルゴリズム

OpenFOAMでは有限体積法を使って離散化を行うことになる。運動方程式は最終的に以下の連立方程式に変換される。

$$ \begin{align} A_P \mathbf{u}_P + \Sigma A_N \mathbf{u}_N &= -\nabla p \notag \end{align} $$

左辺は流速$\mathbf{u}$に関する行列になっているが、対角成分を$A$、非対角成分を$H(=H(\mathbf{u}))$として形式的に表すと次式が得られる。

$$ \begin{align} \mathbf{u} = \frac{H}{A} - \frac{1}{A} \nabla p = \tilde{\mathbf{u}} - \frac{1}{A} \nabla p \notag \end{align} $$

ここで、$\tilde{\mathbf{u}}=\frac{H}{A}$である。この関係式を連続の式に代入することで、圧力に関する方程式が得られる。

$$ \begin{align} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) &= \frac{\partial \psi p}{\partial t} + \nabla \cdot \left[ (\tilde{\rho} + \psi p) \left( \tilde{\mathbf{u}} - \frac{1}{A} \nabla p \right) \right] \notag \\ &= \frac{\partial \psi p}{\partial t} + \nabla \cdot (\tilde{\rho}\tilde{\mathbf{u}}) +\nabla \cdot (\psi \tilde{\mathbf{u}}p) - \nabla \cdot \left( \frac{\rho}{A} \nabla p \right) = 0 \notag \end{align} $$

これらの流速の方程式と圧力の方程式を連立して解くことで解が得られる。なおコーディングでは、$\psi\tilde{\mathbf{u}}$は、phidという変数に格納されている。

以下にsonicLiquidFoam.Cのソースコードの抜粋を示す。

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "compressibleCourantNo.H"

        solve(fvm::ddt(rho) + fvc::div(phi));

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            fvVectorMatrix UEqn  // 運動方程式の計算
            (
                fvm::ddt(rho, U)
              + fvm::div(phi, U)
              - fvm::laplacian(mu, U)
            );

            solve(UEqn == -fvc::grad(p));

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                volScalarField rAU("rAU", 1.0/UEqn.A());
                surfaceScalarField rhorAUf
                (
                    "rhorAUf",
                    fvc::interpolate(rho*rAU)
                );

                U = rAU*UEqn.H();

                surfaceScalarField phid
                (
                    "phid",
                    psi
                   *(
                       fvc::flux(U)
                     + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                    )
                );

                phi = (rhoO/psi)*phid;

                fvScalarMatrix pEqn  // 圧力方程式の計算
                (
                    fvm::ddt(psi, p)
                  + fvc::div(phi)
                  + fvm::div(phid, p)
                  - fvm::laplacian(rhorAUf, p)
                );

                pEqn.solve();

                phi += pEqn.flux();

                solve(fvm::ddt(rho) + fvc::div(phi));
                #include "compressibleContinuityErrs.H"

                U -= rAU*fvc::grad(p);
                U.correctBoundaryConditions();
            }
        }

        rho = rhoO + psi*p;

        runTime.write();

        runTime.printExecutionTime(Info);
    }

3. 摂動法を使った解析手法

超音波の流れ場は基本的にsonicLiquidFoamを使って計算可能であるが、周波数がMHzオーダの超音波の場合、時間刻み幅を振動周期よりも小さくとり、また同様に空間メッシュを波長よりも小さくとる必要があるため、計算時間が膨大になる。

💡 音響流の時間スケール（定常流れの発達）は音波の振動周期より桁違いに長い。摂動法で一次オーダーの音響場と二次オーダーの流れ場を分離することで、それぞれに適した時間スケールで効率的に解析できる。

ここでは摂動法に基づき一次オーダーの音響場と二次オーダーの流れ場を分離して解析する手法を説明する。摂動法では、$\rho,u,p$を微小量$\epsilon$を使って級数に展開する。

$$ \begin{align} \rho &= \rho_0 + \rho’_1 \epsilon + \rho’_2 \epsilon^2 + … \notag \\ \mathbf{u} &= \mathbf{u}_0 + \mathbf{u}’_1 \epsilon + \mathbf{u}’_2 \epsilon^2 + … \notag \\ p &= p_0 + p’_1 \epsilon + p’_2 \epsilon^2 + … \notag \end{align} $$

ここで、添字$0,1,2$はそれぞれ、基準値（音響場なし）、一次オーダー項（音響場）、二次オーダー項を表している。これらの関係式を基礎式に代入し、$\epsilon$のオーダーごとに整理する。

$\epsilon$の一次オーダー項を整理すると、次式が得られる。なおここでは$\rho’_1 \epsilon,$ $\mathbf{u}’_1 \epsilon,$ $\ p’_1 \epsilon$ を $\rho_1,\mathbf{u}_1,p_1$と書き直している。

$$ \begin{align} \frac{\partial \rho_1}{\partial t} &= -\rho_0 \nabla \cdot \mathbf{u}_1 \notag \\ \rho_0 \frac{\partial \mathbf{u}_1}{\partial t} &= -\nabla p_1 + \mu \nabla^2 \mathbf{u}_1 \notag \\ p_1 &= c^2 \rho_1 \notag \end{align} $$

同様に$\epsilon$の二次オーダー項を整理すると、次式が得られる。なおここでは$\rho’_2 \epsilon^2,$ $\mathbf{u}’_2 \epsilon^2,$ $\ p’_2 \epsilon^2$ を $\rho_2,\mathbf{u}_2,p_2$と書き直している。

$$ \begin{align} \frac{\partial \rho_2}{\partial t} &= -\rho_0 \nabla \cdot \mathbf{u}_2 - \rho_1 \nabla \cdot \mathbf{u}_1 - \mathbf{u}_1 \cdot \nabla \rho_1 \notag \\ \rho_0 \frac{\partial \mathbf{u}_2}{\partial t} &= -\nabla p_2 + \mu \nabla^2 \mathbf{u}_2 - \rho_1 \frac{\partial \mathbf{u}_1}{\partial t} - \rho_0 \mathbf{u}_1 \cdot \nabla \mathbf{u}_1 \notag \\ p_2 &= c^2 \rho_2 \notag \end{align} $$

$\rho_1,\mathbf{u}_1,p_1$は一次オーダー項の計算で求められているため既知である。音響に起因する流れ（音響流）はこれらの二次オーダー項の方程式を時間平均することで得られる。以下に時間平均した方程式を示す。$\langle \rangle$で囲んだ項は時間平均されたソース項を表す。

$$ \begin{align} \frac{\partial \rho_2}{\partial t} &= -\rho_0 \nabla \cdot \mathbf{u}_2 + \langle - \nabla \cdot (\rho_1 \mathbf{u}_1 ) \rangle \notag \\ \rho_0 \frac{\partial \mathbf{u}_2}{\partial t} &= -\nabla p_2 + \mu \nabla^2 \mathbf{u}_2 + \langle - \rho_1 \frac{\partial \mathbf{u}_1}{\partial t} - \rho_0 \mathbf{u}_1 \cdot \nabla \mathbf{u}_1 \rangle \notag \\ p_2 &= c^2 \rho_2 \notag \end{align} $$

3.1. 一次オーダー近似式の実装

OpenFOAMの標準ソルバーでは一次オーダー近似の音響場を計算することはできないため、 sonicLiquidFoamを改造してperturbFirstFoamを作成した。

このソルバーは連続式において密度変化を計算（すなわち圧縮性を考慮）し、一方、運動方程式では非圧縮流れを計算している。以下に基礎方程式を再記する。

$$ \begin{align} \frac{\partial \rho_1}{\partial t} &= -\rho_0 \nabla \cdot \mathbf{u}_1 \notag \\ \frac{\partial \mathbf{u}_1}{\partial t} &= -\nabla p_1 + \nu \nabla^2 \mathbf{u}_1 \notag \\ p_1 &= c^2 \rho_1 \notag \end{align} $$

ここで$\nu$は動粘性係数、$p_1$は基準密度で割ったkinematic pressure $[m^2/s^2]$である。この方程式系では密度は無次元値として扱われる。

以下にperturbFirstFoam.Cのソースコードの抜粋を示す。

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            fvVectorMatrix UEqn  // 運動方程式の計算
            (
                fvm::ddt(U)
              - fvm::laplacian(nu, U)
            );

            solve(UEqn == -fvc::grad(p));

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                volScalarField rAU("rAU", 1.0/UEqn.A());
                surfaceScalarField rhorAUf
                (
                    "rhorAUf",
                    fvc::interpolate(rho*rAU)
                );

                U = rAU*UEqn.H();

                surfaceScalarField phid
                (
                    "phid",
                    psi
                   *(
                       fvc::flux(U)
                    )
                );

                phi = (1.0/psi)*phid;

                fvScalarMatrix pEqn  // 圧力方程式の計算
                (
                    fvm::ddt(psi, p)
                  + fvc::div(phi)
                  - fvm::laplacian(rAU, p)
                );

                pEqn.solve();

                phi += pEqn.flux();

                #include "compressibleContinuityErrs.H"

                U -= rAU*fvc::grad(p);
                U.correctBoundaryConditions();
            }
        }

        rho = 1.0 + psi*p;  // 密度を更新

        runTime.write();

        runTime.printExecutionTime(Info);
    }

4. 解析例

解析例を参考にして解析を行った。なお今回は一次オーダーの音響場の計算のみを実施した。二次オーダーの流れ場の計算は今後に実施する予定である。

4.1. 解析モデル

元資料の解析モデルは全体モデルとして解析を行っているが、ここでは対称性を考慮して半分モデルとした。解析領域のサイズおよび境界条件は、図に示す通りである。左下の境界は音源として、圧力を時間変化するsin波として指定している。それ以外は固定壁とした。

圧力指定（音源）：$p=A\sin(\omega t)$（$A=2\times10^4$ [Pa]、$\omega=2\pi f$、$f=10^7$ [Hz]）

図1 — 解析モデル（解析領域と境界条件）

計算に使用したパラメータは以下の通りである。

パラメータ	値
基準圧力 $p_0$	$10^5$ [Pa]
基準密度 $\rho_0$	$1000$ [kg/m³]
粘度 $\mu$	$10^{-3}$ [Pa·s]
音速 $c$	$1440$ [m/s]
圧縮率 $\psi$	$4.82\times10^{-7}$ [s²/m²]
圧力振幅 $A$	$2\times10^4$ [Pa]
振動周波数 $f$	$10^7$ [Hz]