【零基础学智能仿真-40】随机有限元入门——材料不确定性如何影响位移预测

课程摘要

上一节讨论了网格划分带来的数值误差，本节把目光转向材料本身的不确定性。我们以两单元拉杆为例，让两个单元的弹性模量随机变化，重复进行有限元求解，观察端部位移的分布，并估计位移超过指定限值的概率。课程同时介绍空间相关性、对数正态分布、Monte Carlo 统计误差，以及随机场向复杂有限元模型扩展的思路。

一、从“一个答案”到“一组可能的答案”

此前，我们给拉杆指定弹性模量 \\(E=210000\\ \\mathrm{MPa}\\)，计算出唯一的端部位移：

\\[ u_3=\\frac{FL}{EA}=0.476190\\ \\mathrm{mm} \\]

但真实材料不可能处处完全相同。即使是同一根杆，不同位置的弹性模量也可能略有差异。

于是，问题变成：

如果材料参数存在波动，端部位移通常是多少？它有多大可能超过允许值？

这正是随机有限元要回答的问题之一。它不是让有限元求解器“随机地求解”，而是为不确定的输入建立概率模型，再研究输出如何变化。

上一节的网格误差与本节的材料不确定性也不要混为一谈：前者来自计算近似，后者来自我们对物理对象的描述。

二、建立一个看得清楚的随机有限元模型

继续使用课程中的两单元一维拉杆：

• 杆长：\\(L=1000\\ \\mathrm{mm}\\)
• 截面积：\\(A=100\\ \\mathrm{mm^2}\\)
• 右端拉力：\\(F=10000\\ \\mathrm N\\)
• 左端固定
• 两个单元等长：\\(L_e=500\\ \\mathrm{mm}\\)

这次不要求两个单元的弹性模量完全相同，而是分别记为 \\(E_1\\) 和 \\(E_2\\)：

固定端 单元 1 单元 2 拉力
u₁=0 ───────── E₁ ─────────●──────── E₂ ─────────→ F
节点 2 节点 3

对于某一组确定的 \\(E_1,E_2\\)，单元刚度仍按熟悉的公式计算：

\\[ k_1=\\frac{E_1A}{L_e},\\qquad k_2=\\frac{E_2A}{L_e} \\]

固定节点 1 后，待求位移为 \\(u_2,u_3\\)，对应的有限元方程是：

\\[ \\begin{bmatrix} k_1+k_2 & -k_2\\\\ -k_2 & k_2 \\end{bmatrix} \\begin{bmatrix} u_2\\\\u_3 \\end{bmatrix} = \\begin{bmatrix} 0\\\\F \\end{bmatrix} \\]

这个小模型还有一个特别有用的解析结果：

\\[ u_2=\\frac{FL_e}{AE_1},\\qquad u_3=\\frac{FL_e}{A} \\left(\\frac{1}{E_1}+\\frac{1}{E_2}\\right) \\]

因此，我们既能练习“随机输入—有限元求解—统计输出”，又能用解析式逐次检查程序有没有算错。

例如，若 \\(E_1=210000\\ \\mathrm{MPa}\\)、\\(E_2=189000\\ \\mathrm{MPa}\\)，则：

\\[ u_3= \\frac{10000\\times500}{100} \\left(\\frac{1}{210000}+\\frac{1}{189000}\\right) \\approx0.502646\\ \\mathrm{mm} \\]

第二个单元变软后，端部位移比名义值 \\(0.476190\\ \\mathrm{mm}\\) 增大了。

三、怎样描述材料的随机性？

本节作三个明确的教学假设：

弹性模量必须为正，因此示例采用对数正态分布。需要留意：对数正态分布参数对应的是 \\(\\ln E\\) 的正态分布，而不是直接把 \\(E\\) 当作正态变量。SciPy 对数正态分布文档

设弹性模量的目标平均值为 \\(\\bar E\\)、变异系数为 \\(c_v\\)，代码中使用：

\\[ \\sigma_{\\ln E}=\\sqrt{\\ln(1+c_v^2)},\\qquad \\mu_{\\ln E}=\\ln\\bar E-\\frac{\\sigma_{\\ln E}^2}{2} \\]

然后由标准正态变量 \\(Z\\) 生成：

\\[ E=\\exp\\!\\left(\\mu_{\\ln E}+\\sigma_{\\ln E}Z\\right) \\]

代码中的 rho = 0.60 用来控制两个对数弹性模量所对应的正态变量的相关性；它不严格等于变换后 \\(E_1,E_2\\) 的相关系数。对入门实验，先分清这一点即可。

四、完整可运行代码

运行环境只需 NumPy 和 Matplotlib。将下面代码保存为 lesson40_random_bar.py，在选定的 Python 环境中运行：

"""第四十节：两单元随机杆的非侵入式 Monte Carlo 有限元。"""

from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

# N、mm、MPa 单位制：1 MPa = 1 N/mm²
length, area, force = 1000.0, 100.0, 10000.0
e_mean, cv, rho = 210000.0, 0.10, 0.60
n_samples, limit = 10000, 0.55
element_length = length / 2

# 先生成相关标准正态变量，再变换为始终为正的弹性模量。
rng = np.random.default_rng(2026)
z = rng.standard_normal((n_samples, 2))
z[:, 1] = rho * z[:, 0] + np.sqrt(1 – rho**2) * z[:, 1]
sigma_log = np.sqrt(np.log1p(cv**2))
mu_log = np.log(e_mean) – 0.5 * sigma_log**2
moduli = np.exp(mu_log + sigma_log * z)

# 节点 1 固定，只求解节点 2、3 的位移。
tip_displacements = np.empty(n_samples)
for i, (e1, e2) in enumerate(moduli):
k1, k2 = e1 * area / element_length, e2 * area / element_length
k_free = np.array([[k1 + k2, -k2], [-k2, k2]])
u_free = np.linalg.solve(k_free, np.array([0.0, force]))
tip_displacements[i] = u_free[1]

# 每次有限元求解都应与串联杆的解析表达式一致。
exact = force * element_length / area * (
1 / moduli[:, 0] + 1 / moduli[:, 1]
)
np.testing.assert_allclose(
tip_displacements, exact, rtol=1e-12, atol=1e-12
)

nominal = force * length / (area * e_mean)
sample_mean = tip_displacements.mean()
sample_sd = tip_displacements.std(ddof=1)
q025, q975 = np.quantile(tip_displacements, [0.025, 0.975])
p_exceed = np.mean(tip_displacements > limit)
mc_se = np.sqrt(p_exceed * (1 – p_exceed) / n_samples)
theoretical_mean = nominal * (1 + cv**2)

print(f"样本数: {n_samples}")
print(f"确定性名义位移: {nominal:.6f} mm")
print(f"Monte Carlo 平均位移: {sample_mean:.6f} mm")
print(f"理论总体平均位移: {theoretical_mean:.6f} mm")
print(f"位移样本标准差: {sample_sd:.6f} mm")
print(f"位移 2.5%～97.5% 分位区间: [{q025:.6f}, {q975:.6f}] mm")
print(f"P(位移 > {limit:.2f} mm) 估计: {p_exceed:.4%}")
print(f"上述概率估计的 Monte Carlo 标准误: {mc_se:.4%}")
print("检查通过：每个样本的有限元位移与解析式一致。")

fig, ax = plt.subplots(figsize=(9, 5.2), dpi=160)
ax.hist(tip_displacements, bins=55, color="#257d9c", alpha=0.86)
ax.axvline(
nominal, color="#183b56", linestyle="–", linewidth=2,
label=f"Nominal: {nominal:.3f} mm",
)
ax.axvline(
sample_mean, color="#36a386", linewidth=2,
label=f"Sample mean: {sample_mean:.3f} mm",
)
ax.axvline(
limit, color="#d56b32", linewidth=2,
label=f"Illustrative limit: {limit:.2f} mm",
)
ax.set(
xlabel="Tip displacement (mm)",
ylabel="Number of samples",
title="Monte Carlo FEM: uncertain modulus, uncertain displacement",
)
ax.grid(axis="y", alpha=0.18)
ax.legend(frameon=False)
fig.tight_layout()

output = Path(__file__).with_name("lesson40_displacement_distribution.png")
fig.savefig(output)
plt.close(fig)
print(f"图片已保存: {output}")

预期输出

以下是上述代码实际运行得到的结果。使用不同的 NumPy 版本时，随机样本的末几位统计值可能略有差别：

样本数: 10000
确定性名义位移: 0.476190 mm
Monte Carlo 平均位移: 0.480983 mm
理论总体平均位移: 0.480952 mm
位移样本标准差: 0.042945 mm
位移 2.5%～97.5% 分位区间: [0.403046, 0.569974] mm
P(位移 > 0.55 mm) 估计: 6.1900%
上述概率估计的 Monte Carlo 标准误: 0.2410%
检查通过：每个样本的有限元位移与解析式一致。
图片已保存: …/lesson40_displacement_distribution.png

程序生成的实际数据图如下；橙线右侧是超过示例位移限值的样本：

五、读懂结果，而不只是读出数字

第一，为什么平均位移比“平均弹性模量算出的位移”大？

名义位移是 \\(0.476190\\ \\mathrm{mm}\\)，随机样本平均位移约为 \\(0.480983\\ \\mathrm{mm}\\)。原因藏在 \\(u_3\\) 与 \\(1/E\\) 的关系里：材料变软时增加的位移，不能简单地被同幅度“变硬”造成的位移减少抵消。

对于本例的对数正态分布，还可推得理论总体平均值：

\\[ \\mathbb E[u_3] = \\frac{FL}{A\\bar E}(1+c_v^2) = 0.480952\\ \\mathrm{mm} \\]

它与样本平均值接近，是对随机采样和统计代码的第二重检查。

第二，区间和“标准误”回答的是不同问题。

\\([0.403046,\\ 0.569974]\\ \\mathrm{mm}\\) 是模拟位移的 2.5%～97.5% 分位区间，描述在既定材料概率模型下，不同材料样本可能产生怎样的位移。它不是“平均位移的95%置信区间”。

0.2410% 则是超限概率估计的 Monte Carlo 标准误，描述只抽取 10000 次时，所报 6.19% 有多大的抽样波动。这里的 0.2410% 应读作约 0.241 个百分点。

第三，不能把示例限值当作工程判据。

\\(0.55\\ \\mathrm{mm}\\) 是为了练习概率计算而设定的演示限值。实际工程的允许位移必须由具体结构、功能要求和适用规范确定。6.19% 也只在本节假定的材料分布、变异系数及相关性成立时才有意义。

还有一个容易忽略的现象：在这个等截面、右端给定拉力的一维杆中，两单元的轴向应力均为 \\(F/A=100\\ \\mathrm{MPa}\\)。随机弹性模量改变的是应变和位移，不会使本例的轴向应力也随机变化。换成位移控制、变截面或更复杂的结构后，情况可能不同。

六、从两个随机数走向随机场

本节实际上把材料划成两段：每段使用一个随机弹性模量。这是一个很粗的、分段常数的空间随机模型。真实材料可能沿位置 \\(x\\) 连续变化，此时要描述 \\(E(x,\\omega)\\)，还要回答“相距多远的两处材料性质仍然相似”——这便涉及空间相关长度。

一种常见的随机场表示方法是 Karhunen–Loève（KL）展开。例如，为保持弹性模量为正，可以先表示其对数场：

\\[ \\ln E(x,\\omega) \\approx m(x)+ \\sum_{j=1}^{M} \\sqrt{\\lambda_j}\\,\\phi_j(x)\\,\\xi_j(\\omega) \\]

其中，\\(\\phi_j(x)\\) 描述空间变化模式，\\(\\lambda_j\\) 表示相应模式的权重，\\(\\xi_j\\) 是随机系数。保留有限项后，就能用有限个随机变量近似连续随机场。本节代码尚未实现 KL 展开；两个单元各取一个模量，是帮助我们先掌握整个计算链条的起点。随机介质有限元与 KL 展开的原始研究

七、本节练习

本节完成了“随机材料参数 → 重复有限元求解 → 位移分布 → 超限概率”的最小闭环。接下来可以把这种方法移植到二维或三维网格：先为各区域生成具有空间相关性的材料场，再对每个样本调用有限元求解器。