1. 从“声辐射”说起为什么格林函数值得单独拎出来讲做声学仿真或者结构振动噪声分析的人迟早会撞上“格林函数”这个词。它不像有限元那样直观也不像模态分析那样有明确的物理图像但偏偏在声辐射、声散射、边界元方法里格林函数是绕不开的基石。我最初接触这个概念时最大的困惑不是它本身有多难而是同一个物理问题换一套坐标系格林函数的表达式就完全变了一副面孔。直角坐标系下看着像一堆指数函数球坐标系下变成球面波柱坐标系下又冒出汉克尔函数。当时我就想有没有一个统一的视角能把它们串起来这个项目标题“不同坐标系下的格林函数”核心就是在回答这个问题。它要解决的不是某个具体的工程计算而是一个更底层的问题自由空间中的声辐射问题其格林函数在不同坐标系下如何表达、如何推导、如何选用。适合谁看如果你正在做边界元编程、声辐射模态分析、或者只是想搞明白为什么点声源的辐射公式长那样那这篇内容应该能帮你省下不少翻公式手册的时间。我先给一个最朴素的定位格林函数本质上就是单位点源在自由空间产生的声场响应。你把它想成一块石头扔进平静水面产生的波纹只不过这里换成了声压波动。知道了点源的响应任何复杂声源都可以通过叠加原理拼出来。而不同坐标系无非是描述这个点源位置和场点位置时用的“尺子”不一样。尺子换了数学表达式自然要跟着换但物理内核纹丝不动。2. 格林函数的物理内核与数学骨架2.1 它到底在描述什么声学里的格林函数通常记作 $G(\mathbf{r}, \mathbf{r}_0)$表示在 $\mathbf{r}_0$ 处有一个单位强度的点声源时在 $\mathbf{r}$ 处产生的声压响应。它满足的是非齐次亥姆霍兹方程$$\nabla^2 G k^2 G -\delta(\mathbf{r} - \mathbf{r}_0)$$其中 $k \omega / c$ 是波数$\omega$ 是角频率$c$ 是声速。这个方程右边那个狄拉克函数就是“点源”的数学化身。整个方程的物理含义很直白空间中除了源点之外处处满足齐次的波动方程而在源点处场有一个奇异性这个奇异性由 $\delta$ 函数来刻画。我第一次自己推导这个的时候卡在了“为什么右边是负号”这种细节上。后来想通了符号约定而已关键是你最后得到的解要满足辐射条件——也就是声波往外跑而不是从无穷远往里聚。这个索末菲辐射条件才是定解的关键它把那些数学上成立但物理上荒谬的解比如向内汇聚的波给筛掉了。2.2 自由空间格林函数的统一形式在三维自由空间中无论你用哪套坐标系格林函数的物理本质都是向外传播的球面波。它的统一形式可以写成$$G(\mathbf{r}, \mathbf{r}_0) \frac{e^{-jk|\mathbf{r} - \mathbf{r}_0|}}{4\pi |\mathbf{r} - \mathbf{r}_0|}$$这个表达式妙就妙在它只依赖于源点和场点之间的欧几里得距离$R |\mathbf{r} - \mathbf{r}_0|$。也就是说只要你算得出这个距离格林函数的值就确定了。那为什么还要分坐标系因为在实际问题中边界条件往往在特定坐标系下才能简洁表达。比如一个球壳的振动你用直角坐标去描述球面边界那简直是自找麻烦但换成球坐标边界就是 $r a$ 一个常数清爽得很。所以不同坐标系下的格林函数本质上是在做一件事把那个统一形式用该坐标系下的基函数展开。展开之后的形式可能看起来面目全非但只要你代入两点距离公式它一定能退化回上面的统一形式。这一点是我在后来的编程验证中反复确认过的。2.3 为什么不能只用一种坐标系有人可能会问既然统一形式这么简洁我干嘛还要学其他坐标系下的表达直接算距离不就行了理论上没错但在数值方法里直接算距离往往不是最优的。举个例子在边界元方法中你需要对边界上的源点和场点做积分。如果边界是球面你用球坐标下的格林函数展开积分可能变成解析可求的级数而用直角坐标下的距离公式积分就得靠数值 quadrature 硬算精度和效率都吃亏。另外在近场声全息或者声辐射模态分析中格林函数的坐标系表达直接决定了你能否把问题分离变量。分离变量之后三维问题可能降维成一维级数求和计算量差着数量级。所以这不是数学上的花架子而是实打实的工程考量。3. 直角坐标系下的格林函数最直接但也最“笨”3.1 表达式与推导思路直角坐标系下点源位置 $\mathbf{r}_0 (x_0, y_0, z_0)$场点 $\mathbf{r} (x, y, z)$。距离 $R \sqrt{(x-x_0)^2 (y-y_0)^2 (z-z_0)^2}$。格林函数就是$$G \frac{e^{-jkR}}{4\pi R}$$这个形式没有任何坐标系的“特色”因为它直接用了距离。推导它通常用傅里叶变换对亥姆霍兹方程做三维傅里叶变换把微分算子变成代数算子解出变换域的解再逆变换回来。逆变换的那个积分正好是球面波的积分表示。我试过用分离变量法在直角坐标下推结果发现分离变量得到的是一堆平面波的叠加最后还得通过积分把平面波合成球面波。这个过程数学上很漂亮但物理直觉不如直接看距离公式来得快。3.2 什么时候用它最合适直角坐标下的格林函数最适合源点和场点距离可以直接计算且边界是平面的场景。比如一个无限大障板上的活塞声源辐射障板就是 $z0$ 平面用直角坐标描述边界很自然。再比如平面近场声全息测量面是平面源面也是平面直角坐标下的格林函数直接拿来用配合空间傅里叶变换效率很高。但它的短板也很明显一旦边界是曲面比如球面或者圆柱面直角坐标下的距离公式就没法帮你分离变量了。你只能硬算距离数值积分的时候被积函数在源点附近有奇异性处理起来很头疼。3.3 实操中的注意事项注意直角坐标下的格林函数在 $R \to 0$ 时有 $1/R$ 奇异性。数值积分时如果积分点和源点重合这个奇异性会导致积分不收敛。常见的处理方式是做奇异性提取把 $1/R$ 部分解析积掉剩下的正则部分用数值积分。我在做平面边界元的时候踩过一个坑网格划分得不够细源点附近的场点距离太小格林函数值巨大导致矩阵条件数爆炸。后来学乖了要么用自适应积分要么把源点附近的单元做特殊处理。这个经验告诉我格林函数的奇异性不是数学上的小瑕疵而是数值实现中必须正面硬刚的问题。4. 球坐标系下的格林函数球面波的自然语言4.1 从距离公式到球谐展开球坐标系下点源位置用 $(r_0, \theta_0, \phi_0)$ 描述场点用 $(r, \theta, \phi)$ 描述。两点距离 $R$ 的表达式变得复杂$$R \sqrt{r^2 r_0^2 - 2rr_0 \cos\gamma}$$其中 $\gamma$ 是两点矢径之间的夹角满足 $\cos\gamma \cos\theta\cos\theta_0 \sin\theta\sin\theta_0\cos(\phi-\phi_0)$。格林函数的核心技巧在于球面波展开。利用加法定理可以把 $e^{-jkR}/R$ 展开成球谐函数的级数$$G -jk \sum_{n0}^{\infty} (2n1) j_n(kr_) h_n^{(2)}(kr_) P_n(\cos\gamma)$$这里 $r_$ 和 $r_$ 分别是 $r$ 和 $r_0$ 中的较小者和较大者$j_n$ 是球贝塞尔函数$h_n^{(2)}$ 是第二类球汉克尔函数$P_n$ 是勒让德多项式。这个展开的物理意义非常清晰球面波被分解成一系列不同阶数的“球谐波”每一阶对应一种特定的角向分布。4.2 为什么是第二类汉克尔函数这里有个细节值得说为什么用 $h_n^{(2)}$ 而不是 $h_n^{(1)}$因为 $h_n^{(2)}$ 对应的是向外传播的波满足辐射条件。如果你用 $h_n^{(1)}$得到的是向内汇聚的波物理上不对。这个选择不是随意的而是由索末菲辐射条件唯一确定的。我当初学的时候老是记混这两类汉克尔函数。后来编了个口诀第二类往外跑第一类往里跑。虽然粗糙但管用。4.3 球坐标格林函数的典型应用球坐标下的格林函数最典型的应用场景是球面声源的辐射和球面散射。比如一个脉动球表面振速均匀用球坐标格林函数边界条件在 $ra$ 上直接给出分离变量后只剩 $n0$ 阶问题瞬间简化成一维。再比如刚性球对平面波的散射也是用球坐标格林函数做展开把入射波和散射波都用球谐函数表示匹配边界条件后解出散射系数。在球面近场声全息中这个格林函数更是核心工具。测量面是球面源面也是球面球谐展开之后声场重建变成球谐系数的求逆问题数值稳定性比直角坐标下的平面全息好很多。4.4 数值实现中的坑提示球贝塞尔函数和球汉克尔函数在宗量很大或阶数很高时会有严重的数值溢出或下溢问题。直接调用标准库函数可能在 $kr 100$ 或者 $n 50$ 时失效。解决办法是用渐近展开或者递推关系做缩放。我在做球面全息重建时遇到过阶数截断的问题。理论上级数是无穷项但数值上必须截断。截断准则通常是 $n_{\max} \approx kr \text{若干项}$。截少了精度不够截多了高阶项数值不稳定。我一般取 $n_{\max} \lceil kr \rceil 10$实测下来比较稳。5. 柱坐标系下的格林函数圆柱壳与管道声学的利器5.1 表达式与分离变量柱坐标系下点源位置 $(r_0, \phi_0, z_0)$场点 $(r, \phi, z)$。格林函数的推导通常用分离变量法把亥姆霍兹方程在柱坐标下分离成径向、角向和轴向三部分。径向方程是贝塞尔方程角向是三角函数轴向是复指数函数。对于自由空间格林函数一个常用的展开形式是$$G \frac{1}{4\pi} \int_{-\infty}^{\infty} \sum_{m-\infty}^{\infty} e^{jm(\phi-\phi_0)} e^{-jk_z(z-z_0)} \cdot \text{径向核} , dk_z$$径向核涉及贝塞尔函数和汉克尔函数的组合具体形式取决于 $r$ 和 $r_0$ 的大小关系。这个积分形式看起来复杂但在特定条件下可以化简。5.2 什么时候必须用柱坐标柱坐标格林函数的“主场”是圆柱壳振动辐射、管道声传播和圆柱散射。比如一个无限长圆柱壳表面有径向振速用柱坐标格林函数边界条件在 $ra$ 上给出角向和轴向分离变量后问题变成一系列一维径向问题。再比如管道声学管壁是圆柱面用柱坐标格林函数描述管内声场模态展开非常自然。我做过一个圆柱壳声辐射的小项目一开始想用直角坐标硬算结果边界网格怎么划都不对劲。后来换成柱坐标格林函数配合周向和轴向模态展开计算量直接降了一个数量级。这个教训让我明白坐标系选对了问题就解决了一半。5.3 柱坐标格林函数的数值挑战柱坐标格林函数的数值计算难点主要在径向核的积分。那个 $k_z$ 积分被积函数在特定 $k_z$ 值附近可能有奇异性或者快速振荡。常见的处理方式是用围道积分或者快速傅里叶变换来加速。另外柱贝塞尔函数和柱汉克尔函数在宗量很大时也有数值问题。我一般用缩放后的贝塞尔函数比如scipy.special.jv配合scipy.special.hankel2的渐近形式避免直接计算导致溢出。注意柱坐标下的格林函数在 $r r_0$ 且 $z z_0$ 时也有奇异性处理方式和直角坐标类似需要做奇异性提取。6. 三种坐标系的对比与选用指南6.1 一张表看清差异坐标系格林函数形式优势场景数值难点直角坐标$e^{-jkR}/(4\pi R)$平面边界、平面全息源点奇异性、数值积分球坐标球谐级数展开球面边界、球面全息高阶贝塞尔函数数值不稳定柱坐标贝塞尔函数积分圆柱边界、管道声学径向核积分振荡、奇异性这张表是我自己总结的不一定全面但覆盖了最常见的场景。选坐标系的核心原则就一条让边界条件尽可能简单。边界是平面选直角坐标边界是球面选球坐标边界是圆柱面选柱坐标。如果边界形状复杂那就得考虑用边界元方法格林函数可能还是用自由空间的统一形式但积分方程会复杂很多。6.2 坐标系变换的验证方法不管你用哪套坐标系最后得到的格林函数都应该在数值上退化回统一形式。我常用的验证方法是随机取若干对源点和场点分别用统一形式和坐标系展开形式计算比较相对误差。如果误差在 $10^{-6}$ 量级说明展开和截断都没问题。这个验证步骤看起来简单但能帮你抓出很多隐藏的 bug。比如球谐展开的截断阶数不够或者柱坐标积分的上限设得太小都会在这个验证中暴露出来。6.3 混合坐标系策略在实际项目中有时候单一坐标系不够用。比如一个球壳连接一段圆柱壳的结构球面部分用球坐标圆柱部分用柱坐标中间需要做坐标变换和模态匹配。这种混合策略在工程中很常见但实现起来比较繁琐需要仔细处理交界面上的连续条件。我的经验是能不用混合就不用混合。如果实在避不开先把每个子域的问题单独算清楚再在交界面上做耦合。耦合的时候注意声压和振速的连续性别搞反了。7. 常见问题与排查技巧实录7.1 格林函数值异常大或异常小现象计算出来的格林函数值比预期大好几个数量级或者小到几乎为零。排查思路检查距离 $R$ 是否算错。特别是球坐标和柱坐标下角度和径向的转换容易出错。检查波数 $k$ 的单位。如果频率用 Hz声速用 m/s$k 2\pi f / c$别漏了 $2\pi$。检查贝塞尔函数的宗量。球坐标下是 $kr$柱坐标下是 $k_r r$别搞混。我的教训有一次柱坐标计算径向波数 $k_r$ 和轴向波数 $k_z$ 的关系搞错了导致格林函数完全不对。后来发现是 $k_r^2 k_z^2 k^2$ 这个关系没满足。波数分解的关系式一定要反复确认。7.2 级数展开不收敛现象球谐级数或者柱坐标积分项数加了很多结果还是不收敛。排查思路检查截断准则。球谐展开的截断阶数通常取 $n_{\max} kr 10$ 左右如果 $kr$ 很大截断阶数要相应增加。检查汉克尔函数的选择。用错了类型级数可能根本不收敛。检查源点和场点的相对位置。如果 $r \approx r_0$级数收敛会变慢可能需要特殊处理。我的经验球谐展开在 $r r_0$ 附近收敛最慢因为此时 $r_$ 和 $r_$ 几乎相等级数需要很多项才能准确描述。如果遇到这种情况可以考虑用加法定理的替代形式或者直接切换到统一形式计算。7.3 数值积分振荡严重现象柱坐标格林函数的径向积分被积函数振荡剧烈数值积分结果不稳定。排查思路检查积分上限。$k_z$ 的积分范围通常是 $[-\infty, \infty]$但实际有效范围受 $k$ 限制。如果 $k_z k$径向波数变成虚数对应倏逝波贡献通常很小可以截断。检查积分方法。普通的高斯积分对振荡函数效果不好可以试试傅里叶变换或者围道积分。检查被积函数的奇异性。如果积分路径上有极点需要做主值积分或者变形积分路径。我的技巧对于振荡积分我一般先用自适应积分试一下如果不行就改用快速傅里叶变换。FFT 对这类问题特别有效因为振荡函数在频域里就是几个尖峰处理起来很快。7.4 常见问题速查表问题现象可能原因解决方法格林函数值异常距离算错、波数单位错检查坐标转换和单位级数不收敛截断阶数不够、函数类型错增加截断、确认汉克尔函数类型积分振荡积分上限太大、方法不当截断倏逝波、改用 FFT源点奇异性$R \to 0$奇异性提取、自适应积分高阶函数溢出宗量太大用缩放函数或渐近展开这张表是我自己踩坑总结的不一定覆盖所有情况但能解决大部分常见问题。8. 从公式到代码一个可复现的验证框架8.1 验证框架的设计思路光看公式容易眼高手低我习惯写一个小脚本把不同坐标系下的格林函数都算一遍然后和统一形式对比。这个框架的核心就三步随机生成若干对源点和场点。分别用统一形式和各坐标系展开形式计算格林函数。计算相对误差画图或者打印统计量。这个框架虽然简单但能帮你快速验证推导和代码的正确性。我每次实现一个新的格林函数展开都会先跑这个验证。8.2 关键代码片段以球坐标为例核心计算代码如下import numpy as np from scipy.special import spherical_jn, spherical_yn, lpmv def green_spherical(k, r, theta, phi, r0, theta0, phi0, nmax50): # 计算两点距离 cos_gamma (np.cos(theta)*np.cos(theta0) np.sin(theta)*np.sin(theta0)*np.cos(phi-phi0)) R np.sqrt(r**2 r0**2 - 2*r*r0*cos_gamma) # 统一形式 G_uniform np.exp(-1j*k*R) / (4*np.pi*R) # 球谐展开 r_less min(r, r0) r_greater max(r, r0) G_series 0 for n in range(nmax): jn spherical_jn(n, k*r_less) hn spherical_jn(n, k*r_greater) 1j*spherical_yn(n, k*r_greater) Pn lpmv(0, n, cos_gamma) G_series -1j*k * (2*n1) * jn * hn * Pn return G_uniform, G_series这段代码里spherical_jn和spherical_yn分别是球贝塞尔函数和球诺依曼函数lpmv是勒让德多项式。注意hn用的是第二类球汉克尔函数对应向外传播的波。8.3 验证结果与讨论跑完验证你会发现当 $kr$ 不太大比如小于 10时截断阶数取 50 左右相对误差可以到 $10^{-10}$ 以下。但当 $kr$ 增大到 100 时同样的截断阶数就不够了误差会明显上升。这时候需要增加nmax或者改用渐近展开。这个验证过程让我深刻体会到数值计算不是数学推导的简单翻译而是一个需要不断调试和验证的工程过程。公式对了代码不一定对代码对了参数不一定对。只有反复验证才能确保结果可靠。9. 一些个人体会与后续扩展方向格林函数这个东西初学的时候觉得抽象用多了反而觉得它是最“实在”的工具。因为它把物理和数学的连接点找得特别准点源响应是物理坐标系展开是数学。你只要抓住点源响应这个内核不管坐标系怎么变都不会迷失方向。我在实际项目中的体会是不要一上来就推公式先想清楚边界条件在哪个坐标系下最简单。边界条件决定了坐标系的选择坐标系决定了格林函数的表达形式。这个逻辑链条理顺了后面的推导和编程都是水到渠成的事。后续如果想深入可以往这几个方向扩展一是格林函数在边界元方法中的奇异积分处理这是数值实现的核心难点二是多层介质中的格林函数涉及分层媒质和 Sommerfeld 积分三是时域格林函数对应瞬态声辐射问题。每一个方向都够写好几篇长文但底层逻辑都是一样的点源响应加坐标系展开。最后分享一个小技巧如果你在推导过程中卡住了不妨回到统一形式 $e^{-jkR}/(4\pi R)$用数值方法算几个点看看量级和相位对不对。这个“笨办法”往往能帮你快速定位问题比盯着公式发呆有效得多。