AI摘要：本文讨论了环形活塞声场的计算和仿真。环形活塞声场可以视为大圆形活塞减去小圆形活塞，通过积分计算得到声场表达式。文章还介绍了指向性函数的计算方法，并提供了Python代码进行仿真，绘制了指向性图。最后，文章展示了仿真结果的图像。

Powered by AISummary.

超声——环形活塞声场

远场声压

环形活塞场可以看做是大圆形活塞减去小圆形活塞。

根据圆形活塞声场的推导远场声压公式，有：

$p=j\frac{\rho_0C_0k}{2\pi r} u_a e^{j(wt-kr)} \int \sigma d\sigma \int e^{jk\sigma sin\theta cos\varphi}d\varphi$，其中σ的积分范围为0~a，φ的范围为0~2π。

假设外径是a2，内径是a1，则σ的积分范围为a1~a2。则可以看成是0~a2的积分减去0~a1的积分。

$p_{a1}= j\frac{\rho_0wa_1^2}{2 r} u_a e^{j(wt-kr)}\frac{2J_1(ka_1\sin\theta)}{ka_1\sin\theta}=j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{a_1J_1(ka_1\sin\theta)}{k\sin\theta}$

$p_{a2}= j\frac{\rho_0wa_2^2}{2 r} u_a e^{j(wt-kr)}\frac{2J_1(ka_2\sin\theta)}{ka_2\sin\theta}=j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{a_2J_1(ka_2\sin\theta)}{k\sin\theta}$

两者差为：

$p=j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{[a_2J_1(ka_2\sin\theta)-a_1J_1(ka_1\sin\theta)]}{k\sin\theta}$

指向性

$D(\theta)=\frac{P_a(\theta)}{P_a(\theta)|_{\theta=0}}=\frac{\frac{[a_2J_1(ka_2\sin\theta)-a_1J_1(ka_1\sin\theta)]}{k\sin\theta}}{\frac{[a_2J_1(ka_2\sin\theta)-a_1J_1(ka_1\sin\theta)]}{k\sin\theta}|_{\theta=0}}$

θ=0时，$\frac{J_1(0)}{0}=\frac{1}{2}$

$P_a(\theta)|_{\theta=0}=P_{a_2}(\theta)|_{\theta=0}-P_{a_1}(\theta)|_{\theta=0}=j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{a_2^2}{2}-j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{a_1^2}{2}=j\frac{\rho_0w}{r} u_a e^{j(wt-kr)}\frac{a_2^2-a_1^2}{2}$

则上式为：

$D(\theta)=|\frac{\frac{2[a_2J_1(ka_2\sin\theta)-a_1J_1(ka_1\sin\theta)]}{k\sin\theta}}{a_2^2-a_1^2}|=|\frac{2a_2J_1(ka_2\sin\theta)-2a_1J_1(ka_1\sin\theta)}{k\sin\theta(a_2^2-a_1^2)}|$

仿真

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j1  # 导入一阶贝塞尔函数

# %%
def calculate_D(k, a1, a2, start=-np.pi / 2, stop=np.pi / 2, num=1000):
    """
    计算指向性函数
    参数：
        k: 波数
        a1: 半径1
        a2: 半径2
        start: 起始角度（默认 -π/2）
        stop: 终止角度（默认 π/2）
        num: 采样点数
    返回：
        theta数组和D数组
    """
    theta = np.linspace(start, stop, num)
    numerator = 2 * a2 * j1(k * a2 * np.sin(theta)) - 2 * a1 * j1(k * a1 * np.sin(theta))
    denominator = k * np.sin(theta) * (a2 ** 2 - a1 ** 2)
    # 处理分母为0的情况（使用掩码）
    mask = np.abs(denominator) < 1e-10
    D = np.zeros_like(theta)
    D[~mask] = np.abs(numerator[~mask] / denominator[~mask])
    D[mask] = 0  # 这里简单设为0，实际可根据极限情况调整
    return theta, D

# %%
# 创建极坐标绘图画布
fig = plt.figure()
ax = fig.add_subplot(111, projection='polar')  # 通过projection='polar'创建极坐标系

# 计算并处理数据
theta, D = calculate_D(k=3, a1=0.5, a2=1.0)  # 修正参数传递方式
D = np.abs(D)                       # 取绝对值确保指向性值为正


# 绘制极坐标图
ax.plot(theta, D)                    # 在极坐标系中绘制角度-指向性曲线
ax.set_xlabel(r'$\theta$')          # 设置角度轴标签，使用LaTeX数学符号
ax.set_title(r'$ka=3$', loc='left')  # 左对齐标题，显示当前ka乘积值
ax.grid(True)                       # 开启极坐标网格
plt.show()                           # 显示完整图形

# %%
# 绘图
fig, axes = plt.subplots(2, 2, figsize=(12, 12),  # 新增figsize参数调整整体尺寸
                        subplot_kw={'projection': 'polar'})
# 固定 a1 和 a2 的值
a1 = 0.5
a2 = 1.0
k_values = [1, 3, 4, 10]

for ax, k in zip(axes.ravel(), k_values):
    # 修正参数传递，传入 k、a1 和 a2
    theta, D = calculate_D(k=k, a1=a1, a2=a2)
    D = np.abs(D)
    ax.plot(theta, D)
    ax.set_title(f'k={k}', pad=15, loc='left')  # 修改标题显示 k 值
    ax.grid(True)

plt.tight_layout()
plt.show()



# %%
# 新增笛卡尔坐标系绘图
plt.figure(figsize=(8, 6))  # 新建画布
# 定义 ka_values
ka_values = [1, 3, 4, 10]
# 固定 a1 和 a2 的值
a1 = 0.5
a2 = 1.0

for ka in ka_values:
    # 计算 k 值，假设 a = 1
    k = ka
    # 修正参数传递，传入 k、a1 和 a2
    theta, D = calculate_D(k=k, a1=a1, a2=a2, start=0.0001, stop=np.pi*2, num=1000)
    x = np.abs(ka * np.sin(theta))
    plt.plot(x, D, label=f'ka={ka}')
   
    # 新增过零点检测
    zero_crossings = np.where(np.diff(np.sign(D)))[0]  # 检测符号变化点
    if len(zero_crossings) > 0:
        x_zeros = x[zero_crossings] + np.diff(x)[zero_crossings] * (-D[zero_crossings]/(D[zero_crossings+1] - D[zero_crossings]))
        plt.scatter(x_zeros, np.zeros_like(x_zeros), color='red', marker='x', s=50, zorder=5, label='Zero Crossing')
        
        # 新增数值标注
        for xz in x_zeros:
            plt.text(xz, -0.05, f'{xz:.2f}',  # 向下偏移0.05避免重叠
                    ha='center', va='top', 
                    rotation=45, fontsize=8,
                    color='red')

# 添加图例去重处理
handles, labels = plt.gca().get_legend_handles_labels()
by_label = dict(zip(labels, handles))
plt.xlabel(r'kasin($\theta$)')
plt.ylabel(r'$\frac{2J_1(kasin(\theta))}{kasin(\theta)}$', rotation=0)
plt.legend(by_label.values(), by_label.keys())
plt.grid(True)  
plt.tight_layout()
plt.show()

指向性图如下所示。

近场声压

参考圆形活塞声场的处理，不是一般性的，研究的是中心轴上的声场，即x轴。观察点P在x轴上，不再是xoz平面上。则θ=0。

近场时，r不是远大于a，使用z表示r。

$p=j\frac{\rho_0C_0k}{2\pi h}\int dQ_0e^{j(wt-kh)}=j\frac{\rho_0C_0k}{2\pi h} \int u_a dS e^{j(wt-kh)}=j\frac{\rho_0C_0k}{2\pi h} u_a e^{j(wt-kh)} \int dS $

这里环元面积变为外径圆和内径圆的差：$p=j\frac{\rho_0C_0k}{2\pi h} u_a e^{j(wt-kh)}\int_{a_1}^{a_2} 2\pi\sigma d\sigma$

那么环元到观察点的距离分别为：（z为轴向的圆点到观察点的距离）

$h_1=\sqrt{z^2+{a_1}^2}$

$h_2=\sqrt{z^2+{a_2}^2}$

由于$dS=2\pi\sigma d\sigma=2\pi h\partial h$，σ=a1时，h=h1；σ=a2时，h=h2。

故p为：

$p=j\frac{\rho_0C_0k}{2\pi h} u_a e^{j(wt-kh)}\int_0^{a} 2\pi\sigma d\sigma=j\rho_0C_0ku_ae^{jwt}\int_{h_1}^{h_2}e^{-jkh}dh$

$p=j\rho_0C_0ku_ae^{jwt}[\frac{1}{-jk}e^{-jkh}|_z^R]=-\rho_0C_0u_ae^{jwt}[e^{-jkh_2}-e^{-jkh_1}]=-\rho_0C_0u_ae^{jwt}[e^{-jk \sqrt{z^2+{a_2}^2}}-e^{-jk \sqrt{z^2+{a_1}^2}}]$

仿真

# 定义常量（水介质）
rho0 = 998  # 水密度（kg/m³）
C0 = 1480  # 水中声速（m/s）
ua = 1  # 速度幅值（假设值）
f = 10e3  # 频率（Hz）
w = 2 * np.pi * f  # 角频率
k = 2 * np.pi * f / C0  # 波数

# 半径假设值
a1 = 0.5
a2 = 1

def near_field_pressure(z, k, a1, a2):
    """
    根据给定的近场声压公式计算声压
    :param z: 轴向距离
    :param k: 波数
    :param a1: 半径1
    :param a2: 半径2
    :return: 声压值
    """
    h1 = np.sqrt(z ** 2 + a1 ** 2)
    h2 = np.sqrt(z ** 2 + a2 ** 2)
    p = -rho0 * C0 * ua * np.exp(1j * w * 0) * (np.exp(-1j * k * h2) - np.exp(-1j * k * h1))
    return p

# 轴向距离范围（调整为更适合水中高频的尺度）
z_values = np.linspace(0.1, 50, 1000)  # 从0.1mm到10mm

# 计算声压
pressure_values = np.array([near_field_pressure(z, k, a1, a2) for z in z_values])
pressure_magnitude = np.abs(pressure_values)

# 找到最后一个极大值的索引
diff_signal = np.diff(np.sign(np.diff(pressure_magnitude)))
maxima_indices = np.where(diff_signal == -2)[0] + 1  # 极大值点索引
if len(maxima_indices) > 0:
    last_max_index = maxima_indices[-1]
    last_max_x = z_values[last_max_index]
    last_max_y = pressure_magnitude[last_max_index]
    print(f"最后一个极大值的x坐标: {last_max_x}")
    print(f"最后一个极大值的y坐标: {last_max_y}")

# 绘图
plt.figure(figsize=(10, 6))  # 调整图表大小
plt.plot(z_values, pressure_magnitude)

# 标记最后一个极大值点
if len(maxima_indices) > 0:
    # 绘制极大值点
    plt.scatter(last_max_x, last_max_y, color='r', marker='o', s=80)
    
    # 在右侧添加文本注释
    plt.annotate(
        f'({last_max_x:.2f}m, {last_max_y/1e6:.2f}MPa)',
        xy=(last_max_x, last_max_y),
        xytext=(15, 0),  # 文本相对于点的偏移量（向右15个点）
        textcoords='offset points',
        fontsize=10,
        color='red',
        ha='left',  # 水平对齐方式：左对齐
        va='center',  # 垂直对齐方式：居中
        bbox=dict(boxstyle='round,pad=0.3', fc='yellow', alpha=0.7)  # 添加黄色半透明背景
    )

plt.xlabel('Axial Distance z (m)')
plt.ylabel('Sound Pressure Magnitude (MPa)')
plt.title('Near - field Sound Pressure Distribution')
plt.grid(True)
plt.tight_layout()  # 确保标签不被截断
plt.show()