第030篇：流-固-声耦合分析

摘要

流-固-声耦合是描述流体流动、固体结构振动与声波传播之间相互作用的复杂多物理场过程。本主题系统阐述流-固-声耦合的基本理论、数学模型和数值求解方法。首先介绍声学基础理论，包括声波方程、声压与质点速度的关系；然后建立流-固-声耦合控制方程，包括流体域的Navier-Stokes方程、固体域的结构动力学方程和声域的Helmholtz方程；接着讨论流-固界面、固-声界面的耦合条件；最后通过六个典型工程案例展示流-固-声耦合分析的实际应用，包括管道流噪声、风机叶片气动噪声、潜艇声隐身、汽车风噪、消声器设计和建筑声学等。每个案例均配有详细的Python仿真代码和可视化结果，帮助读者深入理解流-固-声耦合问题的建模与求解技术。

关键词

流-固-声耦合，气动声学，结构振动，声辐射，声传输，有限元方法，边界元方法，工程应用

---

1. 声学基础理论

1.1 声波方程

声波是弹性介质中传播的机械波，其传播满足波动方程。对于理想流体介质，声波方程为：

$${\nabla }^{2}p-\frac{1}{c^2}\frac{{\partial }^{2}p}{\partial t^2}=0$$

其中，$p$为声压，$c$为声速，$c=\sqrt{\gamma RT}$，对于空气$c\approx 340$ m/s。

对于简谐声波（时间因子$e^{j\omega t}$），波动方程变为Helmholtz方程：

$${\nabla }^{2}p+k^{2}p=0$$

其中，$k=\omega /c=2\pi /\lambda$为波数，$\omega$为角频率，$\lambda$为波长。

1.2 声场基本参数

声压（Sound Pressure）：声波引起的压强变化，单位Pa。人耳可听声压范围：$2\times 10^{-5}$ Pa（听阈）~ 20 Pa（痛阈）。

声强（Sound Intensity）：单位面积上的声功率：

$$I=\frac{p_{rms}^2}{\rho c}$$

其中，$p_{rms}$为声压有效值，$\rho$为介质密度。

声功率（Sound Power）：声源辐射的总声能流：

$$W={\int }_{S}I\cdot dS$$

声压级（Sound Pressure Level）：

$$L_p=20{\log }_{10}\left(\frac{p_{rms}}{p_{ref}}\right)$$

其中，$p_{ref}=2\times 10^{-5}$ Pa为参考声压。

1.3 声波的反射、透射与吸收

当声波遇到不同介质的界面时，会发生反射、透射和吸收。

反射系数：

$$R={\left(\frac{Z_2-Z_1}{Z_2+Z_1}\right)}^2$$

透射系数：

$$T=\frac{4Z_1{Z}_2}{(Z_2+Z_1{)}^2}$$

其中，$Z=\rho c$为特征阻抗。

吸声系数：

$$\alpha =1-R$$

---

2. 流-固-声耦合控制方程

2.1 流体域控制方程

流体域采用可压缩Navier-Stokes方程描述：

连续性方程：

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$

动量方程：

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)=-\nabla p+\mu {\nabla }^2\mathbf{u}+\frac{1}{3}\mu \nabla (\nabla \cdot \mathbf{u})$$

对于声学问题，通常采用线性化欧拉方程：

$${\rho }_0\frac{\partial \mathbf{u}}{\partial t}=-\nabla p$$

$$\frac{\partial p}{\partial t}=-{\rho }_0{c}^2\nabla \cdot \mathbf{u}$$

2.2 固体域控制方程

固体域采用弹性动力学方程描述：

$${\rho }_s\frac{{\partial }^2\mathbf{d}}{\partial t^2}=\nabla \cdot \mathbf{\sigma }+{\mathbf{f}}_b$$

其中，$\mathbf{d}$为位移矢量，$\mathbf{\sigma }$为应力张量，${\mathbf{f}}_b$为体积力。

应力-应变关系（线弹性）：

$$\mathbf{\sigma }=\mathbf{C}:\mathbf{\epsilon }$$

其中，$\mathbf{C}$为弹性张量，$\mathbf{\epsilon }$为应变张量。

2.3 声域控制方程

声域采用Helmholtz方程描述：

$${\nabla }^{2}p+k^{2}p=-Q$$

其中，$Q$为声源项。

对于结构振动辐射的声场，声源项与结构表面振动速度相关：

$$Q=j\omega {\rho }_0{v}_n\delta (\mathbf{r}-{\mathbf{r}}_s)$$

其中，$v_n$为结构表面法向振动速度。

2.4 耦合界面条件

流-固界面：

• 运动学条件（速度连续）：
  $${\mathbf{u}}_f={\mathbf{u}}_s=\frac{\partial \mathbf{d}}{\partial t}$$

• 动力学条件（应力平衡）：
  $${\mathbf{\sigma }}_f\cdot \mathbf{n}={\mathbf{\sigma }}_s\cdot \mathbf{n}$$

固-声界面：

• 速度连续：
  $$v_n=\frac{\partial d_n}{\partial t}=-\frac{1}{j\omega {\rho }_0}\frac{\partial p}{\partial n}$$

• 压力平衡：
  $$p={\mathbf{\sigma }}_s\cdot \mathbf{n}$$

---

3. 数值求解方法

3.1 有限元方法（FEM）

有限元方法是求解流-固-声耦合问题的主要数值方法。

声场有限元离散：

对Helmholtz方程应用Galerkin方法，得到有限元方程：

$$(\mathbf{K}-k^2\mathbf{M})\mathbf{p}=\mathbf{F}$$

其中：

• $\mathbf{K}$为刚度矩阵：$K_{ij}={\int }_{\Omega }\nabla N_i\cdot \nabla N_{j}d\Omega$
• $\mathbf{M}$为质量矩阵：$M_{ij}={\int }_{\Omega }{N}_i{N}_{j}d\Omega$
• $\mathbf{F}$为载荷向量：$F_i={\int }_{\Omega }Q{N}_{i}d\Omega$

结构有限元离散：

$$({\mathbf{K}}_s-{\omega }^2{\mathbf{M}}_s)\mathbf{d}={\mathbf{F}}_s$$

3.2 边界元方法（BEM）

边界元方法特别适合求解无限域声辐射问题。

Helmholtz边界积分方程：

$$c(\mathbf{r})p(\mathbf{r})={\int }_S\left[\frac{\partial G(\mathbf{r},{\mathbf{r}}^{\prime })}{\partial n^{\prime }}p({\mathbf{r}}^{\prime })-G(\mathbf{r},{\mathbf{r}}^{\prime })\frac{\partial p({\mathbf{r}}^{\prime })}{\partial n^{\prime }}\right]d{S}^{\prime }$$

其中，$G(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{e^{-jk|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|}$为自由场Green函数。

3.3 耦合求解策略

直接耦合法：

将流体、固体、声场方程组装成统一的代数方程组：

$$\left[\begin{array}{ccc}{\mathbf{K}}_{ff}& {\mathbf{K}}_{fs}& 0\\ {\mathbf{K}}_{sf}& {\mathbf{K}}_{ss}& {\mathbf{K}}_{sa}\\ 0& {\mathbf{K}}_{as}& {\mathbf{K}}_{aa}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_f\\ {\mathbf{d}}_s\\ {\mathbf{p}}_a\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_f\\ {\mathbf{F}}_s\\ {\mathbf{F}}_a\end{array}\right]$$

迭代耦合法：

1. 求解流体域，获得流固界面压力
2. 将流体压力作为载荷求解固体响应
3. 将固体振动速度作为边界条件求解声场
4. 检查收敛性，若不收敛返回步骤1

---

4. 工程应用案例

4.1 案例1：管道流噪声分析

问题描述：流体在管道中流动时，由于湍流和边界层分离会产生噪声。分析管道流噪声的产生机理和传播特性。

管道参数：

• 管道直径：$D=0.1$ m
• 管道长度：$L=2$ m
• 流速：$U=10$ m/s
• 流体密度：$\rho =1.225$ kg/m³
• 声速：$c=340$ m/s

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例1：管道流噪声分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题030'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("="*60)
print("案例1：管道流噪声分析")
print("="*60)

# 管道参数
D_pipe = 0.1  # 管道直径 (m)
L_pipe = 2.0  # 管道长度 (m)
U_flow = 10.0  # 流速 (m/s)
rho_air = 1.225  # 空气密度 (kg/m³)
c_sound = 340.0  # 声速 (m/s)
nu = 1.5e-5  # 运动粘度 (m²/s)

# 计算雷诺数
Re = U_flow * D_pipe / nu
print(f"雷诺数: {Re:.0f}")

# 计算马赫数
Ma = U_flow / c_sound
print(f"马赫数: {Ma:.3f}")

# 湍流边界层噪声（简化模型）
# 使用Lighthill声类比理论

# 频率范围
f = np.linspace(50, 5000, 100)  # Hz
omega = 2 * np.pi * f  # rad/s
k = omega / c_sound  # 波数

# 湍流边界层压力脉动（简化模型）
# 使用Corcos模型
def turbulence_spectrum(f, U, delta):
    """湍流压力脉动谱"""
    f_delta = f * delta / U
    # 简化的湍流谱模型
    phi = 1e-4 / (1 + (f_delta/0.5)**2)
    return phi

# 边界层厚度（简化估算）
delta_bl = 0.37 * L_pipe / (Re**0.2)
print(f"边界层厚度: {delta_bl*1e3:.2f} mm")

# 声功率计算
# 使用Curle方程估算声辐射
S_pipe = np.pi * D_pipe * L_pipe  # 管道表面积

# 声功率谱
W_acoustic = np.zeros_like(f)
for i, fi in enumerate(f):
    phi_pp = turbulence_spectrum(fi, U_flow, delta_bl)
    # 简化的声功率计算
    W_acoustic[i] = rho_air * c_sound**3 * phi_pp * S_pipe * (Ma**5)

# 总声功率
W_total = np.trapz(W_acoustic, f)
print(f"总声功率: {W_total*1e6:.2f} μW")

# 声功率级
W_ref = 1e-12  # 参考声功率
L_W = 10 * np.log10(W_total / W_ref)
print(f"声功率级: {L_W:.1f} dB")

# 声压级计算（距离管道1m处）
r_observer = 1.0  # m
I_acoustic = W_total / (4 * np.pi * r_observer**2)
p_rms = np.sqrt(I_acoustic * rho_air * c_sound)
p_ref = 2e-5  # Pa
L_p = 20 * np.log10(p_rms / p_ref)
print(f"1m处声压级: {L_p:.1f} dB")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 声功率谱
ax1 = axes[0, 0]
ax1.semilogx(f, 10*np.log10(W_acoustic/1e-12 + 1e-20), 'b-', linewidth=2)
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Sound Power Level (dB)', fontsize=11)
ax1.set_title('Pipe Flow Noise Spectrum', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_xlim([50, 5000])

# 2. 管道几何示意图
ax2 = axes[0, 1]
pipe_rect = Rectangle((0, -D_pipe/2), L_pipe, D_pipe, 
                       linewidth=2, edgecolor='blue', facecolor='lightblue', alpha=0.5)
ax2.add_patch(pipe_rect)
# 流动箭头
for x in np.linspace(0.2, 1.8, 5):
    ax2.arrow(x, 0, 0.3, 0, head_width=0.01, head_length=0.05, fc='red', ec='red')
ax2.set_xlim([-0.2, 2.2])
ax2.set_ylim([-0.1, 0.1])
ax2.set_xlabel('Length (m)', fontsize=11)
ax2.set_ylabel('Radial (m)', fontsize=11)
ax2.set_title('Pipe Geometry and Flow', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')

# 3. 边界层速度分布
ax3 = axes[1, 0]
y_wall = np.linspace(0, delta_bl*3, 100)
# 对数律速度分布
u_plus = np.log(y_wall * U_flow / nu) / 0.41 + 5.0
u_plus[y_wall == 0] = 0
u_velocity = u_plus * nu / delta_bl
u_velocity = np.clip(u_velocity, 0, U_flow)
ax3.plot(u_velocity, y_wall*1e3, 'g-', linewidth=2)
ax3.set_xlabel('Velocity (m/s)', fontsize=11)
ax3.set_ylabel('Distance from Wall (mm)', fontsize=11)
ax3.set_title('Boundary Layer Velocity Profile', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 4. 声压级随距离衰减
ax4 = axes[1, 1]
r_distances = np.linspace(0.5, 10, 50)
L_p_distance = L_p - 20 * np.log10(r_distances / r_observer)
ax4.plot(r_distances, L_p_distance, 'r-', linewidth=2)
ax4.set_xlabel('Distance (m)', fontsize=11)
ax4.set_ylabel('Sound Pressure Level (dB)', fontsize=11)
ax4.set_title('SPL Attenuation with Distance', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_pipe_flow_noise.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例1完成，图片已保存")

# 噪声评价
if L_p > 85:
    print(f"\n⚠ 噪声超标: {L_p:.1f} dB")
    print("建议降噪措施:")
    print("  1. 增加管道壁厚")
    print("  2. 安装消声器")
    print("  3. 包裹吸声材料")
else:
    print(f"\n✓ 噪声水平可接受: {L_p:.1f} dB")

4.2 案例2：风机叶片气动噪声分析

问题描述：风力机叶片在旋转过程中，由于来流湍流、叶片-塔架干扰和叶尖涡脱落产生气动噪声。分析风机叶片的气动噪声特性。

风机参数：

• 叶片长度：$R=40$ m
• 转速：$\Omega =20$ rpm
• 来流风速：$V_0=12$ m/s
• 叶片数：$B=3$

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例2：风机叶片气动噪声分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Wedge
import warnings
warnings.filterwarnings('ignore')

print("\n" + "="*60)
print("案例2：风机叶片气动噪声分析")
print("="*60)

# 风机参数
R_blade = 40.0  # 叶片长度 (m)
Omega_rpm = 20.0  # 转速 (rpm)
Omega = Omega_rpm * 2 * np.pi / 60  # rad/s
V0 = 12.0  # 来流风速 (m/s)
B_blades = 3  # 叶片数
c_sound = 340.0  # 声速 (m/s)

# 计算叶尖速度
V_tip = Omega * R_blade
print(f"叶尖速度: {V_tip:.1f} m/s")

# 计算叶尖马赫数
Ma_tip = V_tip / c_sound
print(f"叶尖马赫数: {Ma_tip:.3f}")

# 叶片通过频率（BPF）
BPF = B_blades * Omega_rpm / 60
print(f"叶片通过频率: {BPF:.1f} Hz")

# 频率范围
f = np.linspace(20, 2000, 200)  # Hz

# 噪声源模型
# 1. 来流湍流噪声（低频）
def inflow_turbulence_noise(f, V0, R, B):
    """来流湍流噪声"""
    St = f * R / V0  # Strouhal数
    # 简化的湍流噪声模型
    SPL = 50 - 10 * np.log10(1 + (St/0.1)**2)
    return SPL

# 2. 叶尖涡噪声（中频）
def tip_vortex_noise(f, V_tip, R):
    """叶尖涡噪声"""
    St_tip = f * 0.1 * R / V_tip  # 基于叶尖厚度的Strouhal数
    SPL = 60 - 20 * np.log10(1 + (St_tip - 0.2)**2 / 0.01)
    return SPL

# 3. 叶片-塔架干扰噪声（BPF及其谐波）
def blade_tower_interaction(f, BPF, B):
    """叶片-塔架干扰噪声"""
    # 在BPF及其谐波处出现峰值
    harmonics = np.arange(1, 6) * BPF
    SPL = 40 * np.ones_like(f)
    for h in harmonics:
        SPL += 20 * np.exp(-(f - h)**2 / (2 * 10**2))
    return SPL

# 总噪声
SPL_inflow = inflow_turbulence_noise(f, V0, R_blade, B_blades)
SPL_tip = tip_vortex_noise(f, V_tip, R_blade)
SPL_bti = blade_tower_interaction(f, BPF, B_blades)

# 声压级叠加
p_inflow = 10**(SPL_inflow/20)
p_tip = 10**(SPL_tip/20)
p_bti = 10**(SPL_bti/20)
p_total = np.sqrt(p_inflow**2 + p_tip**2 + p_bti**2)
SPL_total = 20 * np.log10(p_total)

# 总声压级（A计权）
def A_weighting(f):
    """A计权滤波器"""
    c1 = 12194**2
    c2 = 20.6**2
    c3 = 107.7**2
    c4 = 737.9**2
    f2 = f**2
    num = c1 * f2**2
    den = (f2 + c2) * np.sqrt((f2 + c3) * (f2 + c4)) * (f2 + c1)
    RA = num / den
    A_weight = 20 * np.log10(RA) + 2.0
    return A_weight

SPL_A = SPL_total + A_weighting(f)
L_A_total = 10 * np.log10(np.trapz(10**(SPL_A/10), f) / (f[-1] - f[0]))
print(f"总A计权声压级: {L_A_total:.1f} dB(A)")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 噪声频谱
ax1 = axes[0, 0]
ax1.semilogx(f, SPL_inflow, 'b--', linewidth=1.5, label='Inflow Turbulence')
ax1.semilogx(f, SPL_tip, 'g--', linewidth=1.5, label='Tip Vortex')
ax1.semilogx(f, SPL_bti, 'r--', linewidth=1.5, label='Blade-Tower Interaction')
ax1.semilogx(f, SPL_total, 'k-', linewidth=2, label='Total')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Sound Pressure Level (dB)', fontsize=11)
ax1.set_title('Wind Turbine Noise Spectrum', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim([20, 2000])

# 2. 风机俯视图
ax2 = axes[0, 1]
theta_blades = np.linspace(0, 2*np.pi, B_blades, endpoint=False)
for theta in theta_blades:
    ax2.arrow(0, 0, R_blade*np.cos(theta), R_blade*np.sin(theta), 
              head_width=2, head_length=3, fc='blue', ec='blue', linewidth=2)
# 塔架
tower = Circle((0, 0), 2, color='gray', alpha=0.5)
ax2.add_patch(tower)
ax2.set_xlim([-50, 50])
ax2.set_ylim([-50, 50])
ax2.set_aspect('equal')
ax2.set_xlabel('X (m)', fontsize=11)
ax2.set_ylabel('Y (m)', fontsize=11)
ax2.set_title('Wind Turbine Top View', fontsize=12, fontweight='bold')
# 来流方向
ax2.arrow(-45, -45, 10, 0, head_width=2, head_length=2, fc='red', ec='red')
ax2.text(-35, -48, 'Wind', fontsize=10, color='red')

# 3. 噪声指向性
ax3 = axes[1, 0]
theta = np.linspace(0, 2*np.pi, 360)
# 简化的指向性模型（偶极子特性）
D_theta = np.sin(theta)**2
ax3.polar(theta, D_theta, 'b-', linewidth=2)
ax3.set_title('Noise Directivity Pattern', fontsize=12, fontweight='bold', pad=20)

# 4. 噪声随距离衰减
ax4 = axes[1, 1]
distances = np.linspace(50, 500, 100)
# 球面扩散 + 大气吸收（简化）
alpha_atm = 0.005  # dB/m
L_p_distance = L_A_total - 20*np.log10(distances/50) - alpha_atm*(distances-50)
ax4.plot(distances, L_p_distance, 'r-', linewidth=2)
ax4.axhline(y=45, color='g', linestyle='--', label='Night Limit (45 dB)')
ax4.axhline(y=50, color='orange', linestyle='--', label='Day Limit (50 dB)')
ax4.set_xlabel('Distance (m)', fontsize=11)
ax4.set_ylabel('Sound Pressure Level (dB(A))', fontsize=11)
ax4.set_title('Noise Propagation', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_wind_turbine_noise.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2完成，图片已保存")

# 噪声合规性检查
if L_A_total > 45:
    print(f"\n⚠ 夜间噪声超标: {L_A_total:.1f} dB(A)")
    print("建议措施:")
    print("  1. 降低转速")
    print("  2. 优化叶片翼型")
    print("  3. 增加轮毂高度")
else:
    print(f"\n✓ 噪声水平符合标准: {L_A_total:.1f} dB(A)")

4.3 案例3：潜艇声隐身分析

问题描述：潜艇的声隐身性能是其在现代海战中生存的关键。分析潜艇壳体振动和声辐射特性，评估声隐身性能。

潜艇参数：

• 艇长：$L=100$ m
• 直径：$D=10$ m
• 壳体厚度：$t=0.05$ m
• 材料：钢，$E=210$ GPa，$\rho =7850$ kg/m³

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例3：潜艇声隐身分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import warnings
warnings.filterwarnings('ignore')

print("\n" + "="*60)
print("案例3：潜艇声隐身分析")
print("="*60)

# 潜艇参数
L_sub = 100.0  # 艇长 (m)
D_sub = 10.0   # 直径 (m)
t_shell = 0.05  # 壳体厚度 (m)
E_steel = 210e9  # 弹性模量 (Pa)
rho_steel = 7850  # 钢密度 (kg/m³)
nu = 0.3  # 泊松比

# 海水参数
rho_water = 1025  # 海水密度 (kg/m³)
c_water = 1500    # 海水中声速 (m/s)

# 频率范围
f = np.linspace(10, 1000, 200)  # Hz
omega = 2 * np.pi * f
k_water = omega / c_water

# 壳体振动模态（简化圆柱壳模型）
# 环向模态频率
n_modes = np.arange(0, 10)  # 环向波数
f_ring = np.zeros_like(n_modes, dtype=float)

for i, n in enumerate(n_modes):
    if n == 0:
        # 呼吸模态
        f_ring[i] = 1/(2*np.pi*D_sub/2) * np.sqrt(E_steel / (rho_steel * (1 - nu**2)))
    else:
        # 环向弯曲模态
        beta = t_shell / (D_sub/2)
        f_ring[i] = n * (n**2 - 1) / (2*np.pi*(D_sub/2)**2) * np.sqrt(E_steel*t_shell**2 / (12*rho_steel*(1-nu**2)))

print("环向模态频率:")
for i, n in enumerate(n_modes[:5]):
    print(f"  n={n}: {f_ring[i]:.1f} Hz")

# 声辐射效率（简化模型）
def radiation_efficiency(f, D, c):
    """圆柱壳声辐射效率"""
    k = 2 * np.pi * f / c
    ka = k * D/2
    # 简化的辐射效率模型
    sigma = np.zeros_like(f)
    for i, kai in enumerate(ka):
        if kai < 1:
            sigma[i] = kai**2 / 2
        else:
            sigma[i] = 1 / np.sqrt(kai)
    return sigma

sigma_rad = radiation_efficiency(f, D_sub, c_water)

# 壳体振动响应（简化的点激励响应）
F_excitation = 1000  # 激励力 (N)
# 假设在第一个环向模态频率处共振
f_resonant = f_ring[2]  # n=2模态
Q_factor = 50  # 品质因数

vibration_velocity = np.zeros_like(f)
for i, fi in enumerate(f):
    # 简化的频率响应
    H = 1 / (1 + 1j * Q_factor * (fi/f_resonant - f_resonant/fi))
    vibration_velocity[i] = np.abs(H) * F_excitation / (rho_steel * t_shell * 2*np.pi*fi)

vibration_velocity = np.clip(vibration_velocity, 0, 0.1)

# 声辐射功率
S_sub = np.pi * D_sub * L_sub  # 湿表面积
W_radiated = 0.5 * rho_water * c_water * S_sub * sigma_rad * vibration_velocity**2

# 声源级
SL = 10 * np.log10(W_radiated / 1e-12) + 170.8  # dB re 1 μPa @ 1m

# 目标强度（TS）
TS = 10 * np.log10(S_sub / (4*np.pi))
print(f"\n目标强度: {TS:.1f} dB")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 壳体模态频率
ax1 = axes[0, 0]
ax1.bar(n_modes[:6], f_ring[:6], color='steelblue', edgecolor='black')
ax1.set_xlabel('Circumferential Mode Number n', fontsize=11)
ax1.set_ylabel('Natural Frequency (Hz)', fontsize=11)
ax1.set_title('Submarine Hull Modal Frequencies', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3, axis='y')

# 2. 辐射效率
ax2 = axes[0, 1]
ax2.semilogx(f, sigma_rad, 'b-', linewidth=2)
ax2.set_xlabel('Frequency (Hz)', fontsize=11)
ax2.set_ylabel('Radiation Efficiency', fontsize=11)
ax2.set_title('Acoustic Radiation Efficiency', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xlim([10, 1000])

# 3. 声源级频谱
ax3 = axes[1, 0]
ax3.semilogx(f, SL, 'r-', linewidth=2)
ax3.set_xlabel('Frequency (Hz)', fontsize=11)
ax3.set_ylabel('Source Level (dB re 1 μPa @ 1m)', fontsize=11)
ax3.set_title('Submarine Acoustic Source Level', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim([10, 1000])

# 4. 潜艇几何示意图
ax4 = axes[1, 1]
# 艇体
hull = Ellipse((0, 0), L_sub/5, D_sub/5, angle=0, 
               facecolor='gray', edgecolor='black', linewidth=2)
ax4.add_patch(hull)
# 指挥塔
conning_tower = Ellipse((15, 8), 20, 6, angle=0,
                        facecolor='darkgray', edgecolor='black', linewidth=2)
ax4.add_patch(conning_tower)
ax4.set_xlim([-60, 60])
ax4.set_ylim([-20, 20])
ax4.set_aspect('equal')
ax4.set_xlabel('Length (m)', fontsize=11)
ax4.set_ylabel('Beam (m)', fontsize=11)
ax4.set_title('Submarine Geometry', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_submarine_acoustics.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3完成，图片已保存")

# 声隐身评估
SL_max = np.max(SL)
if SL_max > 140:
    print(f"\n⚠ 声隐身性能差: 最大声源级 {SL_max:.1f} dB")
    print("建议降噪措施:")
    print("  1. 增加壳体厚度")
    print("  2. 敷设消声瓦")
    print("  3. 优化壳体结构")
else:
    print(f"\n✓ 声隐身性能良好: 最大声源级 {SL_max:.1f} dB")

4.4 案例4：汽车风噪分析

问题描述：汽车高速行驶时，风噪是主要的噪声源之一。分析汽车外形对风噪的影响，优化车身设计。

汽车参数：

• 车长：$L=4.5$ m
• 车宽：$W=1.8$ m
• 车高：$H=1.4$ m
• 车速：$V=120$ km/h

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例4：汽车风噪分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch
import warnings
warnings.filterwarnings('ignore')

print("\n" + "="*60)
print("案例4：汽车风噪分析")
print("="*60)

# 汽车参数
L_car = 4.5  # 车长 (m)
W_car = 1.8  # 车宽 (m)
H_car = 1.4  # 车高 (m)
V_car = 120 / 3.6  # 车速 (m/s)
rho_air = 1.225  # 空气密度 (kg/m³)
c_sound = 340  # 声速 (m/s)
nu = 1.5e-5  # 运动粘度 (m²/s)

# 计算雷诺数
Re = V_car * L_car / nu
print(f"雷诺数: {Re:.2e}")

# 计算马赫数
Ma = V_car / c_sound
print(f"马赫数: {Ma:.3f}")

# 风噪源位置
# 1. A柱（前挡风玻璃与侧窗连接处）
# 2. 后视镜
# 3. 车顶
# 4. 车尾涡流区

# 风噪频谱（简化模型）
f = np.linspace(100, 5000, 200)  # Hz

# 各噪声源贡献
# A柱噪声（宽频）
SPL_A_pillar = 70 + 20*np.log10(V_car/30) - 10*np.log10(1 + (f/1000)**2)

# 后视镜噪声（宽频）
SPL_mirror = 65 + 20*np.log10(V_car/30) - 8*np.log10(1 + (f/800)**2)

# 车顶噪声（低频为主）
SPL_roof = 60 + 20*np.log10(V_car/30) - 15*np.log10(1 + (f/500)**2)

# 车尾涡流噪声（低频）
SPL_wake = 68 + 20*np.log10(V_car/30) - 5*np.log10(1 + (f/200)**2)

# 总风噪（能量叠加）
p_A = 10**(SPL_A_pillar/20)
p_M = 10**(SPL_mirror/20)
p_R = 10**(SPL_roof/20)
p_W = 10**(SPL_wake/20)
p_total = np.sqrt(p_A**2 + p_M**2 + p_R**2 + p_W**2)
SPL_total = 20 * np.log10(p_total)

# 车内噪声（考虑车身隔声）
TL = 25 + 5*np.log10(f/1000)  # 传递损失 (dB)
SPL_interior = SPL_total - TL

# 总声压级
L_A_weighted = SPL_interior + A_weighting(f)
L_A_total = 10 * np.log10(np.trapz(10**(L_A_weighted/10), f) / (f[-1] - f[0]))
print(f"车内总A计权声压级: {L_A_total:.1f} dB(A)")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 风噪频谱
ax1 = axes[0, 0]
ax1.semilogx(f, SPL_A_pillar, 'b--', linewidth=1.5, label='A-Pillar')
ax1.semilogx(f, SPL_mirror, 'g--', linewidth=1.5, label='Side Mirror')
ax1.semilogx(f, SPL_roof, 'r--', linewidth=1.5, label='Roof')
ax1.semilogx(f, SPL_wake, 'm--', linewidth=1.5, label='Wake')
ax1.semilogx(f, SPL_total, 'k-', linewidth=2, label='Total Exterior')
ax1.semilogx(f, SPL_interior, 'orange', linewidth=2, label='Interior')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Sound Pressure Level (dB)', fontsize=11)
ax1.set_title('Automotive Wind Noise Spectrum', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim([100, 5000])

# 2. 汽车侧视图
ax2 = axes[0, 1]
# 车身
body = FancyBboxPatch((0, 0.3), L_car, H_car-0.3, 
                      boxstyle="round,pad=0.1", 
                      facecolor='lightblue', edgecolor='blue', linewidth=2)
ax2.add_patch(body)
# 车轮
wheel1 = plt.Circle((0.8, 0.3), 0.3, color='black')
wheel2 = plt.Circle((3.7, 0.3), 0.3, color='black')
ax2.add_patch(wheel1)
ax2.add_patch(wheel2)
# 后视镜
ax2.plot([0.8, 0.6], [1.0, 1.2], 'k-', linewidth=3)
ax2.add_patch(plt.Circle((0.6, 1.2), 0.15, color='gray'))
# 噪声源标记
ax2.plot([0.5], [1.3], 'ro', markersize=10, label='A-Pillar')
ax2.plot([0.6], [1.2], 'go', markersize=10, label='Mirror')
ax2.plot([2.25], [1.4], 'mo', markersize=10, label='Roof')
ax2.plot([4.3], [0.8], 'ko', markersize=10, label='Wake')
ax2.set_xlim([-0.5, 5.5])
ax2.set_ylim([0, 2])
ax2.set_aspect('equal')
ax2.set_xlabel('Length (m)', fontsize=11)
ax2.set_ylabel('Height (m)', fontsize=11)
ax2.set_title('Vehicle Side View with Noise Sources', fontsize=12, fontweight='bold')
ax2.legend(loc='upper right')

# 3. 噪声随车速变化
ax3 = axes[1, 0]
V_range = np.linspace(60, 180, 50)  # km/h
# 风噪与速度的6次方成正比（简化）
SPL_vs_speed = 60 + 60*np.log10(V_range/60)
ax3.plot(V_range, SPL_vs_speed, 'b-', linewidth=2)
ax3.axhline(y=70, color='r', linestyle='--', label='Comfort Limit')
ax3.set_xlabel('Vehicle Speed (km/h)', fontsize=11)
ax3.set_ylabel('Interior Noise Level (dB(A))', fontsize=11)
ax3.set_title('Interior Noise vs Speed', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. 车身压力分布（简化）
ax4 = axes[1, 1]
x_car = np.linspace(0, L_car, 100)
# 简化的压力系数分布
Cp = np.zeros_like(x_car)
for i, xi in enumerate(x_car):
    if xi < 0.1*L_car:  # 前部高压区
        Cp[i] = 0.8
    elif xi < 0.3*L_car:  # 前部过渡区
        Cp[i] = 0.8 - 1.2*(xi - 0.1*L_car)/(0.2*L_car)
    elif xi < 0.7*L_car:  # 中部低压区
        Cp[i] = -0.2
    else:  # 后部恢复区
        Cp[i] = -0.2 + 0.3*(xi - 0.7*L_car)/(0.3*L_car)

colors = plt.cm.coolwarm((Cp + 0.5))
ax4.scatter(x_car, np.zeros_like(x_car), c=Cp, cmap='coolwarm', s=50)
ax4.plot(x_car, Cp, 'k-', linewidth=2)
ax4.set_xlabel('Position along Vehicle (m)', fontsize=11)
ax4.set_ylabel('Pressure Coefficient Cp', fontsize=11)
ax4.set_title('Surface Pressure Distribution', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.set_ylim([-0.5, 1.0])

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_automotive_wind_noise.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4完成，图片已保存")

# 风噪评价
if L_A_total > 70:
    print(f"\n⚠ 车内噪声偏高: {L_A_total:.1f} dB(A)")
    print("建议优化措施:")
    print("  1. 优化A柱外形")
    print("  2. 改进后视镜设计")
    print("  3. 增强车身密封")
else:
    print(f"\n✓ 车内噪声水平良好: {L_A_total:.1f} dB(A)")

4.5 案例5：消声器声学性能分析

问题描述：消声器是控制排气噪声的关键部件。分析扩张室消声器的声学性能，优化消声设计。

消声器参数：

• 入口管直径：$d_1=0.05$ m
• 扩张室直径：$D=0.15$ m
• 扩张室长度：$L=0.3$ m
• 出口管直径：$d_2=0.05$ m

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例5：消声器声学性能分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

print("\n" + "="*60)
print("案例5：消声器声学性能分析")
print("="*60)

# 消声器参数
d_inlet = 0.05   # 入口管直径 (m)
D_chamber = 0.15  # 扩张室直径 (m)
L_chamber = 0.3   # 扩张室长度 (m)
d_outlet = 0.05   # 出口管直径 (m)
c_sound = 340     # 声速 (m/s)

# 计算扩张比
m_expansion = (D_chamber/d_inlet)**2
print(f"扩张比: {m_expansion:.1f}")

# 频率范围
f = np.linspace(50, 2000, 500)  # Hz
k = 2 * np.pi * f / c_sound

# 扩张室消声器传递损失（理论公式）
def expansion_chamber_TL(f, L, D, d, c):
    """计算扩张室消声器传递损失"""
    k = 2 * np.pi * f / c
    m = (D/d)**2  # 扩张比
    
    # 传递损失公式
    TL = 10 * np.log10(1 + 0.25*(m - 1/m)**2 * np.sin(k*L)**2)
    return TL

TL = expansion_chamber_TL(f, L_chamber, D_chamber, d_inlet, c_sound)

# 最大传递损失频率
f_max_TL = c_sound / (4 * L_chamber)
print(f"最大传递损失频率: {f_max_TL:.1f} Hz")

# 通过频率（传递损失为零）
f_pass = c_sound / (2 * L_chamber)
print(f"通过频率: {f_pass:.1f} Hz")

# 插入损失（考虑源阻抗和负载阻抗）
# 简化的插入损失计算
IL = TL - 1  # dB

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 传递损失频谱
ax1 = axes[0, 0]
ax1.plot(f, TL, 'b-', linewidth=2)
ax1.axvline(x=f_max_TL, color='r', linestyle='--', label=f'f_max = {f_max_TL:.0f} Hz')
ax1.axvline(x=f_pass, color='g', linestyle='--', label=f'f_pass = {f_pass:.0f} Hz')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Transmission Loss (dB)', fontsize=11)
ax1.set_title('Muffler Transmission Loss', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim([50, 2000])

# 2. 消声器几何示意图
ax2 = axes[0, 1]
# 入口管
inlet = Rectangle((0, 0.04), 0.2, 0.02, facecolor='gray', edgecolor='black')
ax2.add_patch(inlet)
# 扩张室
chamber = Rectangle((0.2, 0), L_chamber, D_chamber, 
                    facecolor='lightblue', edgecolor='blue', linewidth=2)
ax2.add_patch(chamber)
# 出口管
outlet = Rectangle((0.2+L_chamber, 0.04), 0.2, 0.02, facecolor='gray', edgecolor='black')
ax2.add_patch(outlet)
# 标注
ax2.annotate('Inlet', xy=(0.1, 0.08), fontsize=10, ha='center')
ax2.annotate('Expansion Chamber', xy=(0.2+L_chamber/2, 0.18), fontsize=10, ha='center')
ax2.annotate('Outlet', xy=(0.4+L_chamber, 0.08), fontsize=10, ha='center')
ax2.set_xlim([-0.1, 0.7])
ax2.set_ylim([-0.05, 0.25])
ax2.set_aspect('equal')
ax2.set_xlabel('Length (m)', fontsize=11)
ax2.set_title('Expansion Chamber Muffler', fontsize=12, fontweight='bold')
ax2.axis('off')

# 3. 不同扩张比的传递损失比较
ax3 = axes[1, 0]
m_values = [4, 9, 16, 25]
colors = ['blue', 'green', 'red', 'purple']
for m, color in zip(m_values, colors):
    D = d_inlet * np.sqrt(m)
    TL_m = expansion_chamber_TL(f, L_chamber, D, d_inlet, c_sound)
    ax3.plot(f, TL_m, color=color, linewidth=2, label=f'm = {m}')
ax3.set_xlabel('Frequency (Hz)', fontsize=11)
ax3.set_ylabel('Transmission Loss (dB)', fontsize=11)
ax3.set_title('Effect of Expansion Ratio', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim([50, 2000])

# 4. 不同长度的传递损失比较
ax4 = axes[1, 1]
L_values = [0.15, 0.3, 0.45, 0.6]
for L, color in zip(L_values, colors):
    TL_L = expansion_chamber_TL(f, L, D_chamber, d_inlet, c_sound)
    ax4.plot(f, TL_L, color=color, linewidth=2, label=f'L = {L*100:.0f} cm')
ax4.set_xlabel('Frequency (Hz)', fontsize=11)
ax4.set_ylabel('Transmission Loss (dB)', fontsize=11)
ax4.set_title('Effect of Chamber Length', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)
ax4.set_xlim([50, 2000])

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_muffler_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5完成，图片已保存")

# 消声器性能评价
TL_avg = np.mean(TL)
if TL_avg < 10:
    print(f"\n⚠ 消声效果不足: 平均传递损失 {TL_avg:.1f} dB")
    print("建议改进措施:")
    print("  1. 增加扩张比")
    print("  2. 增加扩张室长度")
    print("  3. 采用多腔设计")
else:
    print(f"\n✓ 消声效果良好: 平均传递损失 {TL_avg:.1f} dB")

4.6 案例6：建筑声学分析

问题描述：室内声学环境对人们的舒适度和工作效率有重要影响。分析房间声学参数，优化室内声学设计。

房间参数：

• 房间尺寸：$10\times 8\times 3$ m³
• 用途：会议室
• 容积：$V=240$ m³
• 总表面积：$S=268$ m²

Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例6：建筑声学分析
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

print("\n" + "="*60)
print("案例6：建筑声学分析")
print("="*60)

# 房间参数
L_room = 10.0  # 长度 (m)
W_room = 8.0   # 宽度 (m)
H_room = 3.0   # 高度 (m)
V_room = L_room * W_room * H_room  # 容积 (m³)
S_room = 2*(L_room*W_room + L_room*H_room + W_room*H_room)  # 总表面积 (m²)

print(f"房间容积: {V_room:.1f} m³")
print(f"总表面积: {S_room:.1f} m²")

# 材料吸声系数（不同频率）
freq = np.array([125, 250, 500, 1000, 2000, 4000])  # Hz

# 地面（地毯）
alpha_floor = np.array([0.05, 0.10, 0.20, 0.30, 0.40, 0.50])
# 墙面（石膏板）
alpha_wall = np.array([0.30, 0.25, 0.20, 0.15, 0.10, 0.10])
# 天花板（矿棉板）
alpha_ceiling = np.array([0.20, 0.40, 0.60, 0.80, 0.85, 0.80])
# 窗户（玻璃）
alpha_window = np.array([0.35, 0.25, 0.18, 0.12, 0.07, 0.05])

# 各表面面积
S_floor = L_room * W_room
S_ceiling = S_floor
S_wall_total = 2*(L_room + W_room) * H_room
S_window = 4.0  # 窗户面积 (m²)
S_wall = S_wall_total - S_window

# 计算各频率的总吸声量
A_floor = S_floor * alpha_floor
A_ceiling = S_ceiling * alpha_ceiling
A_wall = S_wall * alpha_wall
A_window = S_window * alpha_window

A_total = A_floor + A_ceiling + A_wall + A_window

# 平均吸声系数
alpha_avg = A_total / S_room

# 赛宾混响时间
T60_sabine = 0.161 * V_room / A_total

# 艾润混响时间（考虑空气吸收）
# 空气吸收系数（简化）
m_air = np.array([0.0, 0.0, 0.0, 0.0, 0.01, 0.03])  # 1/m
T60_eyring = 0.161 * V_room / (-S_room * np.log(1 - alpha_avg) + 4*m_air*V_room)

print("\n混响时间:")
for i, f in enumerate(freq):
    print(f"  {f} Hz: T60 = {T60_sabine[i]:.2f} s (Sabine), {T60_eyring[i]:.2f} s (Eyring)")

# 最佳混响时间（会议室）
T60_optimal = 0.5  # s
print(f"\n会议室最佳混响时间: {T60_optimal:.1f} s")

# 语言传输指数（STI）估算（简化）
# STI与混响时间和信噪比相关
SNR = 25  # dB
STI = 1 - np.exp(-0.32 * SNR/30) * np.exp(-0.25 * T60_sabine[3]/2)
print(f"估算语言传输指数 (STI): {STI[3]:.2f}")

# 清晰度评价
if STI[3] > 0.75:
    speech_quality = "优秀"
elif STI[3] > 0.6:
    speech_quality = "良好"
elif STI[3] > 0.45:
    speech_quality = "一般"
else:
    speech_quality = "差"
print(f"语言清晰度: {speech_quality}")

# 声压级分布（简化的扩散场模型）
# 声源功率
W_source = 0.001  # W (正常说话)
# 房间常数
R_room = A_total / (1 - alpha_avg)
# 距离声源不同位置的声压级
r_distances = np.linspace(1, 8, 100)
SPL_direct = 10*np.log10(W_source / (4*np.pi*r_distances**2)) + 120  # 直达声
SPL_reverb = 10*np.log10(4*W_source/R_room) + 120  # 混响声
# 总声压级
SPL_total_room = 10*np.log10(10**(SPL_direct/10) + 10**(SPL_reverb/10))

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 1. 混响时间频谱
ax1 = axes[0, 0]
ax1.semilogx(freq, T60_sabine, 'bo-', linewidth=2, markersize=8, label='Sabine')
ax1.semilogx(freq, T60_eyring, 'rs-', linewidth=2, markersize=8, label='Eyring')
ax1.axhline(y=T60_optimal, color='g', linestyle='--', linewidth=2, label=f'Optimal ({T60_optimal}s)')
ax1.set_xlabel('Frequency (Hz)', fontsize=11)
ax1.set_ylabel('Reverberation Time T60 (s)', fontsize=11)
ax1.set_title('Room Reverberation Time', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim([100, 4000])

# 2. 材料吸声系数
ax2 = axes[0, 1]
ax2.plot(freq, alpha_floor, 'b-o', linewidth=2, markersize=6, label='Floor (Carpet)')
ax2.plot(freq, alpha_wall, 'r-s', linewidth=2, markersize=6, label='Wall (Gypsum)')
ax2.plot(freq, alpha_ceiling, 'g-^', linewidth=2, markersize=6, label='Ceiling (Mineral Wool)')
ax2.plot(freq, alpha_window, 'm-d', linewidth=2, markersize=6, label='Window (Glass)')
ax2.set_xlabel('Frequency (Hz)', fontsize=11)
ax2.set_ylabel('Absorption Coefficient', fontsize=11)
ax2.set_title('Material Absorption Coefficients', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim([100, 4000])
ax2.set_ylim([0, 1])

# 3. 声压级随距离变化
ax3 = axes[1, 0]
ax3.plot(r_distances, SPL_direct, 'b--', linewidth=2, label='Direct Sound')
ax3.plot(r_distances, SPL_reverb[3]*np.ones_like(r_distances), 'r--', linewidth=2, label='Reverberant Sound')
ax3.plot(r_distances, SPL_total_room[3], 'g-', linewidth=2, label='Total SPL')
ax3.set_xlabel('Distance from Source (m)', fontsize=11)
ax3.set_ylabel('Sound Pressure Level (dB)', fontsize=11)
ax3.set_title('SPL Distribution in Room', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. 房间平面示意图
ax4 = axes[1, 1]
room_outline = Rectangle((0, 0), L_room, W_room, 
                         fill=False, edgecolor='black', linewidth=2)
ax4.add_patch(room_outline)
# 声源位置
ax4.plot([L_room/2], [W_room/2], 'ro', markersize=15, label='Sound Source')
# 听众位置
listener_positions = [(3, 3), (5, 4), (7, 3), (4, 6), (6, 6)]
for i, pos in enumerate(listener_positions):
    ax4.plot([pos[0]], [pos[1]], 'bs', markersize=10)
ax4.set_xlim([-1, 11])
ax4.set_ylim([-1, 9])
ax4.set_aspect('equal')
ax4.set_xlabel('Length (m)', fontsize=11)
ax4.set_ylabel('Width (m)', fontsize=11)
ax4.set_title('Room Layout', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_room_acoustics.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6完成，图片已保存")

# 声学设计评价
T60_avg = np.mean(T60_sabine)
if T60_avg > 1.0:
    print(f"\n⚠ 混响时间过长: 平均 {T60_avg:.2f} s")
    print("建议改进措施:")
    print("  1. 增加吸声材料")
    print("  2. 使用吸声天花板")
    print("  3. 添加吸声屏风")
elif T60_avg < 0.3:
    print(f"\n⚠ 混响时间过短: 平均 {T60_avg:.2f} s")
    print("建议改进措施:")
    print("  1. 减少吸声材料")
    print("  2. 使用反射性表面")
else:
    print(f"\n✓ 混响时间合适: 平均 {T60_avg:.2f} s")

print("\n" + "="*60)
print("所有案例仿真完成！")
print("="*60)

print("\n生成的图片文件:")
print("  - case1_pipe_flow_noise.png")
print("  - case2_wind_turbine_noise.png")
print("  - case3_submarine_acoustics.png")
print("  - case4_automotive_wind_noise.png")
print("  - case5_muffler_analysis.png")
print("  - case6_room_acoustics.png")