Matlab间断有限元求解声波方程：高阶格式实现与稳定计算

简介这套Matlab程序面向数值计算与物理仿真方向的学习者和科研人员以二维声波方程为例完整演示了间断有限元DG方法的编程实现适用于课程教学、毕业设计或科研预研。程序采用间断伽辽金格式进行空间离散并配合三阶龙格库塔法完成时间积分适合需要从公式推导走向代码实操的读者。压缩包为RAR格式共包含6个文件其中5个为M脚本、1个为自动保存的ASV备份文件整体仅有4KB代码精简、无冗余模块便于逐行理解网格划分、基函数构造、空间离散化、时间迭代与边界条件处理等核心步骤。借助正弦初始波场程序可在单位矩形区域内模拟声波传播过程并通过Matlab绘图直观呈现解随时间的演化。目前已有1564人学习下载对正在钻研偏微分方程数值解或着手波动模拟的入门及中级开发者具有很好的参考价值与复用性。1. 间断有限元求解声波方程波动模拟里值得认真考虑的一种格式做声波正演的人经常会在连续有限元和高阶格式之间纠结。波动问题对色散误差特别敏感高频分量在长距离传播中如果被数值格式吃掉或者拖慢波形就会“散”掉。间断有限元求解声波方程的做法是把计算域切成大量单元每个单元内部用多项式近似单元之间允许解本身不连续通过数值通量把信息跨界面传过去。这套思路既保留了高阶精度又让质量矩阵天然分块便于显式时间推进尤其适合处理非协调网格和复杂分界面。如果你正在做地震波模拟、超声检测、室内声场这类需要长时间、宽频带计算的场景又希望用 matlab 快速验证算法间断有限元是一个投入产出比很高的选项。下面按从格式选择到可运行代码、再到参数与踩坑的路径展开。2. matlab里实现间断有限元求解声波方程先做一阶系统还是直接离散二阶方程2.1 常规连续有限元在波动问题上的两个硬伤连续有限元处理声波方程时必须在单元边界上强制解连续。这个约束对静力学问题很自然对波动问题却并不总是合理。声波方程真正沿特征线传播的信息是 p±Zv 这两个特征量而不是压力 p 或速度 v 本身。强制 p 跨单元连续等于给格式加了一个物理上并不存在的约束结果往往是高频部分色散误差变大波形尾部拖出振荡。另一个实际问题是质量矩阵。连续有限元的全局质量矩阵是耦合的显式时间推进时每一大步都要解一个稀疏线性方程组。虽然质量集中技巧可以缓解但在一维高阶单元上质量集中后的精度损失很直接。间断有限元在每个单元上独立做 Galerkin 离散单元之间只靠数值通量联系质量矩阵是分块对角甚至完全对角的matlab 里用逐单元矩阵运算就能把主循环写得很干净这对想快速看到结果的从业者很友好。2.2 一阶压力-速度系统与迎风数值通量的推导常见做法是先把二阶声波方程拆成一阶系统。一维均匀介质中声波方程可以写成$$ \frac{\partial p}{\partial t} \rho c^2 \frac{\partial v}{\partial x} 0,\qquad \frac{\partial v}{\partial t} \frac{1}{\rho}\frac{\partial p}{\partial x} 0 $$其中 p 是声压v 是质点速度ρ 是密度c 是声速。也可以写成矩阵形式 U_t A U_x 0其中 U[p; v]。矩阵 A 的特征值就是 ±c对应的右特征向量分别是 [ρc; 1] 和 [ρc; -1]所以 pρcv 以速度 c 向右传播p-ρcv 以速度 c 向左传播。这个分解直接决定了单元边界的数值通量形式。对单元 e 做分部积分后单元边界上会出现一个需要单独定义的界面值。取迎风通量时界面处的压力数值通量和速度数值通量分别为$$ \hat p \frac{1}{2}\left(p_L p_R Z(v_L - v_R)\right),\qquad \hat v \frac{1}{2}\left(v_L v_R \frac{1}{Z}(p_L - p_R)\right) $$其中 Zρc 是声阻抗下标 L、R 分别表示界面左侧和右侧单元的内部极限值。这组表达式看起来两边都取了平均但实际上加了阻抗项作用是把信息按特征方向偏置向右传播的波主要从左侧取信息向左传播的波主要从右侧取信息。注意这两个公式里 p̂ 和 v̂ 的符号特别容易写反后面避坑章节会专门说。2.3 LDG 的适用场景什么时候值得多引入一个变量如果坚持离散二阶声波方程本身就要引入辅助变量 q∇p得到 p_tr、r_tc²∇·q、q∇p 这样的三变量一阶系统。这种局部间断 Galerkin 写法的关键点是 q 的数值通量必须与 p 的数值通量按“交替方向”配对否则可能不稳定或者损失收敛阶。LDG 的好处是能直接处理带二阶算子的方程形式某些弹性波、频散波方程里更自然。但就纯声波方程而言我一般不推荐入门直接上 LDG。三套变量意味着自由度和内存开销接近翻倍一维里还需要额外推导交替通量的组合工程收益并不明显。先做 p-v 一阶系统跑通之后如果遇到必须保持二阶方程形式的场景再回头对照 LDG 的辅助通量选择这样踩坑最少。matlab 里这两类格式的骨架几乎一样差别只在残差函数里变量数量和通量表达式的不同。3. 在matlab中搭最小可运行的间断有限元求解声波方程脚本3.1 模态DG的最小数据结构单元、基函数和系数矩阵这里用模态间断有限元实现比节点型更容易理解也更容易写出紧凑代码。把一维区间 [0,L] 均分成 N 个单元每个单元长度 hL/N单元中心记作 xmid(e)。每个单元内部用 P 阶勒让德多项式做基函数参考区间取 [-1,1]。物理坐标与参考坐标的变换是 x xmid(e) h/2·ξ求导时需要多乘一个 2/h。所有单元按列存成一个 (P1)×N 的矩阵 p 和 v每一列是该单元内的模态系数。勒让德基函数正交单元质量矩阵是对角阵 M_i2/(2i1)所以时间导数项的求解只是简单的逐点除法这是 matlab 里速度最快也最好调试的部分。基函数在左右边界上的取值是固定的φ_i(-1)(-1)^iφ_i(1)1可以提前算好。3.2 残差函数与RK3主循环代码下面这段代码是整个求解器的核心。第一步先准备高斯-勒让德节点和基函数矩阵第二步写残差函数第三步用强稳定保留的三阶 Runge-Kutta 推进时间。% gauss_legendre.m 用三对角特征值法生成节点和权重 function [x,w] gauss_legendre(n) beta (1:n-1)./sqrt(4*(1:n-1).^2-1); J diag(beta,1) diag(beta,-1); [V,D] eig(J); x diag(D); w 2*V(1,:).^2; % 第一个特征向量分量生成权重 [x,idx] sort(x); w w(idx); end % legendreBasisDeriv.m 返回 P 阶以内所有勒让德基函数值与导数值 function [L,dL] legendreBasisDeriv(P, xi) p P 1; m length(xi); L zeros(p,m); dL zeros(p,m); L(1,:) 1; if p 1 L(2,:) xi; dL(2,:) 1; end for n 2:p-1 L(n1,:) ((2*n-1)*xi.*L(n,:) - (n-1)*L(n-1,:))/n; dL(n1,:) ((2*n-1)*(L(n,:)xi.*dL(n,:)) - (n-1)*dL(n-1,:))/n; end end % residual_pv.m 一维声波方程模态DG残差迎风数值通量 function [dpdt, dvdt] residual_pv(p, v, phiL, phiR, Minv, S, rho, c, h) P1 size(p,1); N size(p,2); Z rho*c; pL phiL*p; pR phiR*p; vL phiL*v; vR phiR*v; p_hat zeros(N1,1); v_hat zeros(N1,1); % 左边界一阶无反射 p_hat(1) 0.5*(pL(1) Z*vL(1)); v_hat(1) 0.5*(vL(1) pL(1)/Z); % 内部界面迎风通量 for j 2:N p_hat(j) 0.5*(pR(j-1) pL(j) Z*(vR(j-1) - vL(j))); v_hat(j) 0.5*(vR(j-1) vL(j) (pR(j-1) - pL(j))/Z); end % 右边界一阶无反射 p_hat(N1) 0.5*(pR(N) Z*vR(N)); v_hat(N1) 0.5*(vR(N) pR(N)/Z); dpdt zeros(P1,N); dvdt zeros(P1,N); for e 1:N bnd_v v_hat(e1)*phiR - v_hat(e)*phiL; bnd_p p_hat(e1)*phiR - p_hat(e)*phiL; dpdt(:,e) -rho*c^2 * Minv * (bnd_v - (2/h)*S*v(:,e)); dvdt(:,e) -(1/rho) * Minv * (bnd_p - (2/h)*S*p(:,e)); end end这里 Minv 是参考单元质量矩阵的逆因为勒让德基正交Minv diag((2*(0:P)1)/2)。S 是参考单元上的刚度贡献矩阵S(i,j)∫φ_j φ_i dξ在高斯节点上通过数值积分一次性算好。残差函数里每一列 p(:,e)、v(:,e) 是单元 e 的模态系数phiL*p 是把模态系数换算成界面处物理值这一步相当于把单元边界上的“跳变量”显式暴露出来。RK3 主循环采用 SSP 格式比普通 RK3 在强间断附近更稳。每步先对当前状态计算一次残差再对中间状态计算残差最后加权更新。子步代码短但能保证总变差不增这对声波方程里的光滑波包不是必需遇到界面反射系数突变时却能避免明显过冲。% 主循环参数 rho 1.0; c 1.0; L 20; N 200; P 2; h L/N; xmid linspace(h/2, L-h/2, N); % 提前计算算子矩阵 [gq, gw] gauss_legendre(P2); [phiq, dphiq] legendreBasisDeriv(P, gq); phiL legendreBasisDeriv(P, -1); phiR legendreBasisDeriv(P, 1); Minv diag((2*(0:P)1)/2); S phiq * diag(gw) * dphiq; % 数值积分组装刚度贡献 % 初值高斯压力脉冲速度为零 x0 5; sigma 0.4; p zeros(P1,N); v zeros(P1,N); for e 1:N xq xmid(e) h/2*gq; p0q exp(-((xq-x0)/sigma).^2); p(:,e) (h/4)*(2*(0:P)1).*(phiq*(gw(:).*p0q(:))); end % 时间推进SSP-RK3 cfl 0.25; dt cfl*h/(c*(P1)); T 3.0; nstep ceil(T/dt); dt T/nstep; for n 1:nstep [dp1,dv1] residual_pv(p,v,phiL,phiR,Minv,S,rho,c,h); p1 p dt*dp1; v1 v dt*dv1; [dp2,dv2] residual_pv(p1,v1,phiL,phiR,Minv,S,rho,c,h); p2 0.75*p 0.25*(p1 dt*dp2); v2 0.75*v 0.25*(v1 dt*dv2); [dp3,dv3] residual_pv(p2,v2,phiL,phiR,Minv,S,rho,c,h); p (1/3)*p (2/3)*(p2 dt*dp3); v (1/3)*v (2/3)*(v2 dt*dv3); end3.3 初值投影与参数设置高斯脉冲的初值投影必须用足够高阶的高斯积分。这里取了 P2 个高斯点对 P 阶勒让德基而言积分精确度足够不会因为投影不准而在初始时刻就引入毛刺。投影系数公式里那个 (h/4) 因子容易漏物理单元长度是 h参考区间积分要乘 h/2而投影系数本身还有 (2i1)/2二者合起来就是 h/4。如果你改成二阶中心通量残差函数里通量公式要同步换但算子矩阵不用改。dt 由 cflh/(c(P1)) 给出。P 越大允许时间步越小这是高阶 DG 的普遍代价。参数 rho1、c1 是归一化介质实际算空气、水或钢材时直接换这两个标量即可。p 和 v 的物理单位会随之变化通量公式里 Z 也自动跟着变。4. 数值通量、边界条件和时间步让声波方程计算长期稳定4.1 三种常用数值通量的选择迎风通量不是唯一选择。中心通量也很常见表达式是界面两侧简单平均实现最容易但耗散偏小对网格质量和时间步长更敏感容易出现高频毛刺。局部 Lax-Friedrichs 通量在一维常系数声波方程里和迎风通量表达式完全相同区别在于它的阻抗项来自对特征速度最大值的估计因此在变系数介质中更常用。三种通量的实际差异在光滑介质中很小在密度或声速突变界面才能明显体现出来。做介质分界面时迎风通量的阻抗项会更自然地表达“波从低声阻抗一侧向高声阻抗一侧传播时反射系数为正是反相还是同相”这类问题。如果只做均匀介质验证中心通量也够用但建议从一开始就用迎风通量后面扩展到两层介质时不改框架。4.2 边界条件自由表面、刚性壁与吸收边界声学里最常见的两类反射边界一维里实现起来非常简洁。刚性壁要求界面处 v0在残差函数里把对应边界的数值通量直接改成 p̂p_L、v̂0 即可自由表面要求 p0则改成 p̂0、v̂v_L。这两类边界条件在 matlab 里改的是 p_hat(1) 和 v_hat(1) 两行代码不需要动主循环。无反射边界稍微复杂一点。第3章代码里用的是“一阶无反射”本质是假定边界处只有出射波、没有入射波推导出来的界面通量$$ \hat p \frac{1}{2}(pZv),\qquad \hat v \frac{1}{2}(vp/Z) $$对正入射平面波完全透明但对斜入射和大角度波前会有明显反射。如果模拟时间不长、边界离源足够远这个边界够用。长时间计算或需要贴边界激发源时一阶边界反射会污染目标区域的波形那时候就得考虑 PML。PML 在一维的实现是在控制方程里加阻尼项和辅助变量篇幅不小建议先跑通一阶边界再升级。4.3 SSP-RK3与CFL条件稳定步长的经验取值时间积分用标准 RK4 也可以但 SSP-RK3 在 DG 框架里更常用。三种值得注意的写法是时间格式稳定性每步残差调用次数适用场景显式欧拉最差1只做格式验证不实用SSP-RK2中等2低阶基函数快速试探SSP-RK3最好3默认选择尤其 P≥2CFL 条件里 P 的影响是真实的。一维标量波的 DG 稳定性分析给出 dt ≤ C·h/(c(2P1))其中 C 是常数SSP-RK3 大概能到 0.4~0.8。第3章代码取 0.25 是很保守的值目的是让读者第一次运行必然稳定。实际调试时我习惯先按 0.25 跑通确认波形形态正确后再逐步放大到 0.5 左右一旦出现界面振荡就退回。用同样的网格把 P 从 1 提高到 3稳定时间步可能缩到原来的 1/3 甚至更小这是高阶 DG 绕不开的成本。时域点源可以用右端项投影实现。下面这段代码把力源加到速度方程上适用于模拟激振源或换能器辐射% 点源投影在单元 s 内、物理坐标 xsrc 处注入 f(t) dphi_src legendreBasisDeriv(P, (xsrc - xmid(s)) / (h/2)); source_shape (h/2) * Minv * dphi_src; % (P1)x1 % 每个时间步更新时叠加 dvdt(:,s) dvdt(:,s) source_shape * f_t / rho;f_t 是当前时刻的激励幅值。这里 source_shape 相当于把狄拉克源投影到基函数上相当于该点处的等效体力。用这种方法加源比在初值里塞波形更灵活可以方便地换成雷克子波、正弦扫频或任意实测信号。注意源点在单元边界上时投影会遇到左右单元归属问题通常把源整体放进左单元即可。5. 间断有限元求解声波方程的常见翻车点与排查5.1 现象一调大时间步就发散减小后长时间不收敛出现这个现象时第一反应是检查 dt 是否超过 CFL 限制。时间步略超临界值波形通常不是立即爆炸而是运行几百步后界面处出现高频振荡然后整场发散。很多人会把责任推给 RK3实际上问题常在通量表达式里的 Z 用错。Zρc但在某些 matlab 程序里 c 被赋成网格速度 h/dt或者 ρ 被写成 1c² 导致阻抗项整体偏大格式的耗散和耗散方向会变得不正常。解决方法是把通量公式里的 Z 单独抽出来打印验证。在均匀介质里让左右两个单元状态完全相等数值通量应该严格等于单元内的物理通量即 p̂p_L、v̂v_L。如果这个基本一致性不满足时间步再小也救不回来。先做这种半数值半解析的检查再去调 dt至少能排除一半“玄学发散”。5.2 现象单元界面附近出现锯齿振荡锯齿振荡通常集中在界面处单元中心曲线看起来还光滑。最常见的原因是初值投影的高斯节点数太少。如果只用 P 阶基函数对应的 P 个节点做积分高斯求积本身就截断了高频初始模态系数容易带有虚假振荡。把高斯点从 P 提升到 P2 一般就够了。另一个直接原因是模态系数里高频分量在单元边界处过冲。勒让德基底在边界上取值不为零高 P 时吉布斯现象放得更明显。处理办法是把初值先按更高阶解析表达式采样再投影回 P 阶空间而不是直接用低阶插值。很多网上流传的 matlab 有限元编程求解实例用节点值直接当模态系数这在模态 DG 里会造成界面永远有跳变看起来就会像“只有界面处有锯齿”。5.3 现象吸收边界产生明显反射一阶无反射边界在很多模型里完全够用但反射总是存在的。它的反射系数在正入射时接近零在波前掠射时明显增大。如果模型中源贴近边界或者波从多层介质返回时边界附近角度不再规则反射波幅值可能达到入射波的 5%~15%在声压云图上表现为边界处一圈圆弧慢慢向内部扩散。这个现象的根源不是代码 bug是边界数学模型本身。一维情况下可以用 PML 替换一阶边界需要新增两个辅助变量和阻尼剖面。二维里 PML 要处理角点区域参数设置更敏感。我在用一阶无反射边界时会把介质区向外延伸至少两个主波长让反射波回来得足够晚再用时间窗截掉它这是最省事的工程做法。5.4 现象加密网格后误差不下降理论上 P 阶基函数配 SSP-RK3误差应该随网格加密呈现 P1 阶收敛。如果加密后误差几乎不动先检查时间步是否同步缩小了。只加密空间但 dt 不变时间误差占主导收敛阶会被锁死。另一个隐蔽问题是质量矩阵逆写错Minv 写成 diag(2*(0:P)1) 而不是除以 2这会让高阶模态在时间推进时被放大波形整体失真。还有一个常见误用是把跨单元的平均值当成节点值来画图导致看起来误差巨大。模态系数不是节点值绘图前要先用基函数在采样点上重建物理场。一个简单检查是让 P0此时每个单元只有一个常数格式退化为有限体积法如果这个退化版都算不对问题一定出在通量或边界条件上如果退化版对而 P 增大后乱套再检查刚度矩阵 S 的组装方向。6. 用一个简单技巧验证求解器对称性和积分误差同时看写完一个间断有限元求解声波方程程序后最优先做的不是看波形动画而是看两个量总能量衰减和 wave arrival time。一维无吸收均匀介质中高斯脉冲从 x0 出发经过时间 T 后中心应到达 x0cT。先用单元平均画一个粗略的 A-scan如果峰的位置对说明通量方向、边界条件、时间步整体正确如果峰的位置差得远多半是质量矩阵或特征方向的错误。波形峰值高度对可进一步说明耗散误差小能量损失主要来自一阶无反射边界。收敛阶验证可以做一个简单测试表固定 dt 随 h 同步缩小算 L2 误差并记录误差比。P1 时加密一倍误差应缩小约 4 倍P2 时缩小约 8 倍。如果比例明显偏低检查残余函数里 h/2 因子和时间步缩放。这个验证习惯能把“看起来对”变成“量起来也对”。从一维走向二维时matlab 的矩阵化写法仍然适用但有几个点会变。张量积网格上每个矩形单元里的基函数是 x 和 y 方向勒让德基的乘积自由度从 P1 变成 (P1)²。界面积分沿每条边分别计算数值通量公式里要多一个法向速度分量。二维下 PML 参数调起来比一维麻烦得多我一般的做法是先做均匀介质行波验证再做圆形波前验证最后才加 PML 和复杂介质。遇到二维算不动的情况先把一维程序里所有数组预分配和向量化全部做掉因为 matlab 的循环在二维网格上很慢每一层单元循环都会让运行时间成倍增长。最终我的习惯是保留一份带打印开关的最小验证脚本每次改动基函数阶数或边界条件后都先跑一次收敛阶测试再开始算大模型。这才把间断有限元的高精度真正握在手里。希望帮到你。本文还有配套的精品资源点击获取