一维弹性杆数值验证：解析解、误差与能量检查

一根杆来回振动，曲线平滑、计算没有发散，能否说明求解器正确？还不够。应力符号可能反了，振动周期可能偏了，也可能每次往返都损失一部分振幅。只截取一个时刻，甚至会把耗散更强的结果判为“更准确”。

本文用固定两端的一维线弹性杆，把验证拆成三个互相补充的方向：解析场与相位是否一致、细化后误差如何变化、整个过程中能量怎样交换。所有数值表对应弹性杆解析验证实验的实际推进结果，不是把解析曲线加扰动后当作数值解。

需要的背景是一维波动方程、线弹性关系和粒子—网格插值。速度梯度为何控制材料变形，可先读物质点法运动学；本文只保留验证所需的一维更新式。

一、把参考问题约定清楚

取参考坐标$X\in[0,L]$向右，位移$u(X,t)$向右为正，速度$v=u_t$，轴向应力$\sigma$以拉伸为正。均匀杆的截面积为$A$，密度为$\rho$，杨氏模量为$E$。忽略体力、阻尼与横向运动，小应变关系为$\varepsilon=u_X$、$\sigma=E\varepsilon$。轴向平衡给出

$$\rho A u_{tt}=\frac{\partial(A\sigma)}{\partial X} \quad\Longrightarrow\quad \rho u_{tt}=E u_{XX}. \tag{1}\label{bar-wave}$$

这里截面积恒定，因此从运动方程中约去，但计算质量、内力和能量时仍要保留。$X,u$的单位为$\mathrm{m}$，$t$为$\mathrm{s}$，$v$为$\mathrm{m/s}$，$E,\sigma$为$\mathrm{Pa}$，$\rho$为$\mathrm{kg/m^3}$，应变无量纲。

两端固定，初始无位移、无应力，只给定第一模态的速度：

$$\begin{aligned} &u(0,t)=u(L,t)=0,\\ &u(X,0)=0,\\ &v(X,0)=v_0\sin\!\left(\frac{\pi X}{L}\right). \end{aligned}$$

固定端同时意味着$v(0,t)=v(L,t)=0$，并不意味着端部应力为零。约束反力可以存在；把固定端误写成零应力端，就已经换了一个参考问题。

实验采用以下固定 SI 数据，不在界面中调节材料参数：

参数	数值	含义
$L$	$1\ \mathrm{m}$	杆长
$E$	$10^6\ \mathrm{Pa}$	杨氏模量
$\rho$	$1000\ \mathrm{kg/m^3}$	密度
$A$	$0.01\ \mathrm{m^2}$	截面积
$v_0$	$0.001\sqrt{E/\rho}$	初速度幅值，单位$\mathrm{m/s}$

二、从单模态推到位移、速度与应力

令波数$k=\pi/L$，写成$u(X,t)=q(t)\sin(kX)$，两端条件自动满足。代入式$\eqref{bar-wave}$，得到

$$\ddot q+\omega^2q=0, \qquad c=\sqrt{E/\rho},\qquad \omega=kc,$$

其中$c$是波速，单位$\mathrm{m/s}$；$k$的单位为$\mathrm{m^{-1}}$，角频率$\omega$的单位为$\mathrm{s^{-1}}$。由$q(0)=0$、$\dot q(0)=v_0$，有

$$\begin{aligned} u(X,t)&=\frac{v_0}{\omega}\sin(kX)\sin(\omega t),\\ v(X,t)&=v_0\sin(kX)\cos(\omega t),\\ \sigma(X,t)&=\frac{Ev_0}{c}\cos(kX)\sin(\omega t),\\ T&=\frac{2\pi}{\omega}=\frac{2L}{c}. \end{aligned} \tag{2}\label{bar-exact}$$

应力式来自对位移求空间导数，不是凭振型猜测：$\partial_X\sin(kX)=k\cos(kX)$，且$k/\omega=1/c$。因此在第一四分之一周期，杆的左半段受拉、右半段受压，中心应力为零。

本例$c\approx31.6228\ \mathrm{m/s}$，$\omega\approx99.3459\ \mathrm{s^{-1}}$，$T\approx0.0632456\ \mathrm{s}$，$v_0\approx0.0316228\ \mathrm{m/s}$。记位移、应力幅值为$u_{\mathrm{amp}}$、$\sigma_{\mathrm{amp}}$，则

$$u_{\mathrm{amp}}=\frac{v_0}{\omega}\approx318.3099\ \mathrm{\mu m}, \qquad \sigma_{\mathrm{amp}}=\frac{Ev_0}{c}=1\ \mathrm{kPa}.$$

最大应变为$v_0/c=0.001$，与小应变假设一致。下图画的是$t=T/4$的解析场：横轴为参考位置，位移呈正弦形，应力呈余弦形；此时全杆解析速度为零。

检验相位时，不必先看复杂误差曲线，可以先核对四个特征时刻：

时刻	位移	速度	应力
$0$	零	正向正弦，幅值$v_0$	零
$T/4$	正向正弦峰值	零	左拉右压
$T/2$	零	反向正弦，幅值$v_0$	零
$T$	零	回到初始速度场	零

数值解无需逐位等于解析解，但方向、相位顺序和量纲必须先对。若这一层失败，讨论收敛阶没有意义。

三、实验究竟推进了什么

这个求解器是固定参考网格上的线性小应变物质点离散，不是随粒子当前位置重算权重的 updated-Lagrangian MPM。粒子有随时间变化的位移记录，但参考坐标$X_p$不移动；权重、梯度、质量与体积全程固定。不能用本例验证网格穿越误差或有限变形能力。

设单元数为$N$，网格宽度$h=L/N$，每单元粒子数为$n_{pc}$。粒子放在均匀子区间中点，总数$n_p=Nn_{pc}$；其体积与质量为

$$V_p=\frac{Ah}{n_{pc}},\qquad m_p=\rho V_p.$$

$N_i=N_i(X_p)$是节点$i$的线性形函数值，无量纲；$N'_i$是对$X$的导数，单位$\mathrm{m^{-1}}$。初始$u_p,\varepsilon_p,\sigma_p$均为零，$v_p=v_0\sin(kX_p)$。每一步按以下顺序计算。

先映射动量并组装内力。 上标$n$表示当前时间层，节点质量$m_i$只需预计算一次；动量与力则每步重新累加：

$$m_i=\sum_p m_pN_i,\qquad \bar v_i=\frac{\sum_p m_pv_p^nN_i}{m_i},\qquad f_i=-\sum_pV_p\sigma_p^nN'_i,\qquad a_i=\frac{f_i}{m_i}.$$

$f_i$的单位是$\mathrm{N}$，负号对应拉正应力的弱形式。端点处理发生在节点更新前：把两端映射速度$\bar v_i$与加速度$a_i$都置零，再计算

$$v_i^{n+1}=\bar v_i+\Delta t\,a_i.$$

再回传速度并更新材料状态。 两种回传方式只在粒子速度更新处不同：

$$v_p^{n+1}= \begin{cases} v_p^n+\Delta t\sum_iN_i a_i,&\text{FLIP},\\ \sum_iN_i v_i^{n+1},&\text{PIC}. \end{cases}$$

这里 FLIP 只加上已施加约束的加速度增量，不把未约束端点的映射速度差混入回传。两种方法都用更新后的节点速度推进位移与应变：

$$\begin{aligned} u_p^{n+1}&=u_p^n+\Delta t\sum_iN_i v_i^{n+1},\\ \varepsilon_p^{n+1}&=\varepsilon_p^n+\Delta t\sum_iN'_i v_i^{n+1},\\ \sigma_p^{n+1}&=E\varepsilon_p^{n+1}. \end{aligned}$$

因此，FLIP 情况下位移并不是直接用粒子存储速度$v_p^{n+1}$推进。把这两个速度混为一谈，会复现出另一套算法。循环中的应力与内力都由数值状态得到；解析式只用于初始化速度和事后对照，没有逐步“喂入”精确力或速度，也没有位置裁剪、反弹或额外人工阻尼。回传机制的含义另见 PIC 与 FLIP 速度回传。

四、过零时刻怎样定义误差

常见相对误差把误差范数除以当前解析场的范数。但本例在$t=T$时$u=\sigma=0$，在$t=T/4$时$v=0$；这些分母为零或只剩浮点舍入残差。得到一个极大百分数，并不代表数值解真的爆炸。

实验改用固定解析峰值 RMS 尺度。对$q\in\{u,v,\sigma\}$，定义

$$e_q(t)= \frac{\displaystyle\sqrt{\frac{1}{n_p}\sum_p \left[q_p(t)-q_{\mathrm{exact}}(X_p,t)\right]^2}} {q_{\mathrm{amp}}/\sqrt{2}}, \qquad q_{\mathrm{amp}}\in \left\{\frac{v_0}{\omega},v_0,\frac{Ev_0}{c}\right\}. \tag{3}\label{bar-error}$$

分子是同一批参考坐标上的均方根误差；均匀等体积粒子允许这里使用等权平均。分母是正弦或余弦振型达到峰值时的空间 RMS，不随时间改变。界面与表格显示$100e_q\%$，并非以“该时刻真实值”为分母的相对百分比。三种物理量虽有不同单位，经各自尺度归一化后都可比较。

过零问题解决了，单时刻偏好却没有消失。$t=T$的解析位移为零，一个振幅被强烈削弱的数值解也可能接近零；它的位移误差看起来很好，速度却远没有恢复初值。必须把三个误差一起看，并检查周期中间的相位。

五、一组可以复现的联合细化结果

取 FLIP、每单元$4$个粒子、CFL 为$0.20$，终止于$t=T$。界面的 CFL 控制的是时间步上限：若终止时间为$t_f$，则

$$n_t=\left\lceil\frac{t_f}{\mathrm{CFL}\,h/c}\right\rceil, \qquad \Delta t=\frac{t_f}{n_t}.$$

$n_t$为步数，实际$\Delta t$可能略小于上限，以精确到达相同终止相位。本组四行的实际时间步与上限相等。下表来自求解器的终态采样，能量记号$\mathcal E$与杨氏模量$E$区分。

单元数$N$	步数$n_t$	$\Delta t$/ s	位移误差 / %	速度误差 / %	应力误差 / %	$\mathcal E/\mathcal E_0$
8	80	$7.905694\times10^{-4}$	11.8967	1.8325	11.9781	0.978034
16	160	$3.952847\times10^{-4}$	3.0351	0.1872	3.0403	0.997183
32	320	$1.976424\times10^{-4}$	0.7618	0.0206	0.7621	0.999646
64	640	$9.882118\times10^{-5}$	0.1906	0.0024	0.1906	0.999956

图中横、纵轴均为对数尺度，绘制的是表中同一组数据。位移和应力误差随单元数加倍约降为四分之一；速度在此终止相位下降更快。这个现象只能描述为本组参数、当前相位下的联合细化趋势，不能直接宣称“空间二阶、速度更高阶”。固定 CFL 时，$h$与$\Delta t$同时减小，粒子数也增加；特殊相位还可能使某些误差项抵消。

PIC 对照更能说明只看位移的风险。同样取$N=16$，PIC 在$t=T$的位移误差为$1.2771\%$，小于 FLIP 的$3.0351\%$；但其速度误差为$41.0798\%$，总能量只剩初始的$0.347324$。到$N=64$，PIC 位移误差降至$0.1185\%$，速度误差仍有$12.4103\%$，能量比为$0.767198$。这不是“位移指标错误”，而是零位移终态没有充分约束振幅和速度。

六、能量是另一条证据，不是准确性的替身

没有阻尼和外部功输入，固定端速度为零，约束反力不做功。连续参考解的动能$K$、应变能$U$与总能量满足

$$\begin{aligned} K(t)&=\int_0^L\frac12\rho A v^2\,\mathrm dX =\mathcal E_0\cos^2(\omega t),\\ U(t)&=\int_0^L\frac{A\sigma^2}{2E}\,\mathrm dX =\mathcal E_0\sin^2(\omega t),\\ \mathcal E(t)&=K(t)+U(t)=\mathcal E_0, \qquad \mathcal E_0=\frac{\rho ALv_0^2}{4}=0.0025\ \mathrm J. \end{aligned}$$

$t=0$时全为动能，$T/4$时全为应变能；两者周期性交换，总量不变。数值诊断则从粒子状态独立计算

$$K_h=\sum_p\frac12m_pv_p^2,\qquad U_h=\sum_p\frac{V_p\sigma_p^2}{2E},\qquad \mathcal E_h=K_h+U_h.$$

界面以离散初始动能归一化；对本例均匀中点粒子分布，它与上述$0.0025\ \mathrm J$一致。能量历史的横轴为$t/T$，纵轴为各能量除以初始能量，解析总能量是一条数值为$1$的水平线。

至少查看一到两个完整周期：动能与应变能的峰谷是否按顺序交换？总能量是有界波动、持续下降，还是增长？终态回到接近$1$，不代表中途没有偏离；总能量接近常数，也不能排除相位滞后。PIC 的反复投影可造成显著耗散，而本实现中的 FLIP 同样不是严格能量守恒格式。

因此应分别作出判断：稳定关注误差是否失控，守恒关注离散预算，准确关注数值场是否接近参考解。三者有关联，但不能互相替代。

七、在实验室中按同一口径复现

1. 打开弹性杆解析验证，设“单元数”为$16$、“每单元粒子数”为$4$、CFL 为$0.20$、“求解时长 / 周期 T”为$1.00$，回传选择 FLIP，点击“运行验证”。核对“已运行算例”记录的$160$步与实际时间步，而不只看左侧待运行参数。
2. 拖动“已保存快照”，以显示的$t/T$为准检查初态、四分之一周期、半周期与终态。空间横轴始终是参考坐标$X$；位移单位为微米，速度为米每秒，应力为千帕。展开“当前快照粒子数据”可逐点核对数值与解析值。粒子连线只辅助阅读，不代表连续应力重构。
3. 点击“网格收敛对照”，复现上一节四行结果。该表使用当前粒子数设置、CFL、求解时长与回传方式，独立计算四种网格的终态，不随快照滑块改变。
4. 切换 PIC 后重新运行并生成对照，重点比较$t=T$的位移、速度与能量，而不是只挑最小的一列。随后将时长改为$2.00T$并重新运行，查看完整能量历史与中间相位。

修改参数后，已有图表会标为过期，但仍保留上一次算例；快照只是选择保存结果，不会重新求解。网格对照缓存会随参数修改失效。把“当前控件”误当成“当前图表所用参数”，是复现时很容易漏掉的一步。

八、下一步怎样把误差来源分开

上面的实测结果是联合细化，不是完整收敛研究。若继续分析，可分别组织三类实验，并始终保持相同物理终止时间和误差尺度：

• 时间细化：固定网格与每单元粒子数，只改变 CFL，比较多个相位及整个区间的误差。PIC 减小时间步会增加固定时间内的重映射次数，投影耗散可能反而加强；更小的$\Delta t$不保证每个指标都改善。
• 空间细化：为各网格选取相同且足够小的时间步，先确认时间误差不主导，再比较网格。现有界面通过 CFL 控制步长，需要换算并核对实际值；若控件范围不足，就不能把固定 CFL 的按钮结果冒充此实验。
• 粒子积分细化：固定网格和时间步，改变每单元粒子数，区分积分与粒子采样的影响。粒子增加不会提高线性网格基函数的阶次。

区间误差可以采用保存时刻上的最大值或时间 RMS，同时报告采样方式；稀疏快照可能漏掉峰值。本文没有给出上述独立研究的实测阶次，也不把推荐步骤当作已经取得的结果。

最后，验证结论应附上范围：本例只包含均匀杆、线弹性、小应变、线性基函数与固定参考离散，不涉及移动粒子过网格、接触、塑性、不同基函数优劣或 APIC。界面允许的 CFL 区间$[0.05,0.5]$是这个均匀固定网格算例的保守选区，不是通用 MPM 稳定性保证。

提交一份可复核的验证记录前，逐项检查：边界与拉正符号是否一致；单位和振幅是否正确；数值推进是否独立于解析场；误差分母是否避开瞬时零范数；比较是否处于同一相位；是否记录实际步长、粒子数与回传方式；是否同时检查三个场和能量历史。只有这些证据相互支持，“算出了振动”才进一步成为“这套离散在明确条件下得到了可信结果”。

本文含有 AI 生成内容