课程摘要上一节生成了沿杆长变化的随机弹性模量场本节将这一方法扩展到二维。我们以受拉矩形薄板为例在三角形单元上生成空间相关的弹性模量反复组装刚度矩阵并求解统计右边界位移。课程还通过均匀材料解析解和整体力平衡检查程序说明随机有限元中“生成样本、映射网格、验证求解、解释概率”缺一不可。一、从一维求和进入真正的二维求解第四十一节的拉杆有一个便利之处每个单元承受同样的轴力。因此端部位移可以写成各单元伸长量之和。二维薄板则不同。某块区域变软后周围区域可能分担更多载荷位移和应力必须通过整体刚度矩阵共同求解。我们要完成的计算链条是\[ \boxed{ \text{生成二维随机材料场} \longrightarrow \text{赋值给三角形单元} \longrightarrow K^{(s)}u^{(s)}f \longrightarrow \text{统计结构响应} } \]上标 \(s\) 表示第 \(s\) 次材料抽样。载荷和几何保持不变每次改变的是单元弹性模量因而刚度矩阵也随之改变。本节沿用第二十三节介绍过的 NumPy 常应变三角形单元CST思路让新增的“随机场—单元映射”过程容易核查。**本节代码是二维 NumPy 教学求解器并非声称已经在 FEniCSx 中运行。**文末再说明如何迁移到 FEniCSx。二、案例随机材料矩形薄板建立长 \(100\ \mathrm{mm}\)、高 \(20\ \mathrm{mm}\)、厚 \(1\ \mathrm{mm}\) 的薄板采用平面应力假设左边界约束水平位移左下角额外约束竖向位移以消除整体刚体运动。右边界施加 \(100\ \mathrm{MPa}\) 的均布水平拉应力对应合力 \(2000\ \mathrm N\)。泊松比固定为 \(\nu0.3\)。弹性模量平均值设为 \(210000\ \mathrm{MPa}\)变异系数设为 \(10\%\)。材料对数场的相关长度设为 \(25\ \mathrm{mm}\)。教学网格由6×3个矩形组成每个矩形剖分成两个三角形共36个三角形单元。每个三角形取一个弹性模量值。我们关心的输出不是某个任意节点的位移而是右边界各节点水平位移的平均值\[ Q^{(s)} \frac{1}{n_R} \sum_{a\in\Gamma_R}u_x^{(s)}(a) \]这样可以得到一个明确的、可重复统计的响应量。请注意这并不表示右边界所有节点的位移都相同。三、先建立确定性“标准答案”在均匀材料、单轴均布拉伸的理想条件下平面应力解析位移场为\[ u_x(x,y)\frac{p}{E}x,\qquad u_y(x,y)-\nu\frac{p}{E}y \]其中 \(p100\ \mathrm{MPa}\)。所以右端水平位移应为\[ u_x(L,y) \frac{pL}{E} \frac{100\times100}{210000} 0.047619\ \mathrm{mm} \]这一步很重要如果程序连均匀材料都算不对随后运行500次随机求解也只会得到500个值得怀疑的数字。本节代码会先用这个解析场检查单元刚度、载荷与边界约束。四、把空间随机场映射到二维单元对第 \(e\) 个三角形取其中心坐标 \(\boldsymbol c_e\)。与上一节类似我们为对数弹性模量设置距离相关的协方差\[ C_{ef} \sigma_G^2 \exp\left( -\frac{\|\boldsymbol c_e-\boldsymbol c_f\|}{\ell} \right) \]这一次距离是二维欧氏距离而不再只是沿杆长的距离。对协方差矩阵做特征值分解再生成各单元的随机弹性模量 \(E_e^{(s)}0\)。本课保留全部离散模态避免在讲“映射与求解”的同时混入 KL 截断误差。平面应力单元刚度为\[ K_e^{(s)} tA_e B_e^{\mathsf T}D(E_e^{(s)},\nu)B_e \]其中 \(t\) 为厚度\(A_e\) 为三角形面积\(B_e\) 是应变—位移矩阵。把每个 \(K_e^{(s)}\) 放到对应的整体自由度位置便得到本次样本的整体刚度矩阵。五、完整可运行代码第四十二节二维平面应力 CST 有限元与随机弹性模量场。 from pathlib import Path import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np # 单位N、mm、MPa右侧均布拉应力为 100 MPa。 L, H, THICKNESS, PRESSURE 100.0, 20.0, 1.0, 100.0 E_MEAN, NU, CV, CORRELATION_LENGTH 210000.0, 0.3, 0.10, 25.0 NX, NY, N_SAMPLES, LIMIT 6, 3, 500, 0.050 x, y np.meshgrid( np.linspace(0, L, NX 1), np.linspace(0, H, NY 1) ) nodes np.column_stack((x.ravel(), y.ravel())) triangles [] for j in range(NY): for i in range(NX): a j * (NX 1) i b, c, d a 1, a NX 1, a NX 2 triangles.extend(((a, b, d), (a, d, c))) triangles np.asarray(triangles, dtypeint) centers nodes[triangles].mean(axis1) ndof 2 * len(nodes) # 平面应力的单位弹性模量本构矩阵。 d_unit np.array([ [1, NU, 0], [NU, 1, 0], [0, 0, (1 - NU) / 2] ]) / (1 - NU**2) element_dofs, element_k_unit [], [] for tri in triangles: xy nodes[tri] xx, yy xy[:, 0], xy[:, 1] area np.linalg.det( np.column_stack((np.ones(3), xx, yy)) ) / 2 assert area 0 b np.roll(yy, -1) - np.roll(yy, 1) c np.roll(xx, 1) - np.roll(xx, -1) bmat np.zeros((3, 6)) bmat[0, 0::2], bmat[1, 1::2] b, c bmat[2, 0::2], bmat[2, 1::2] c, b bmat / 2 * area element_dofs.append(np.array([ 2 * node k for node in tri for k in (0, 1) ])) element_k_unit.append( THICKNESS * area * bmat.T d_unit bmat ) # 右侧边均匀施加拉应力。 load np.zeros(ndof) for j in range(NY): lower j * (NX 1) NX upper (j 1) * (NX 1) NX edge_force PRESSURE * THICKNESS * (H / NY) / 2 load[2 * lower] edge_force load[2 * upper] edge_force left np.flatnonzero(np.isclose(nodes[:, 0], 0.0)) right np.flatnonzero(np.isclose(nodes[:, 0], L)) fixed np.r_[2 * left, 1] free np.setdiff1d(np.arange(ndof), fixed) def solve(modulus_by_element): stiffness np.zeros((ndof, ndof)) for modulus, dofs, k_unit in zip( modulus_by_element, element_dofs, element_k_unit ): stiffness[np.ix_(dofs, dofs)] modulus * k_unit displacement np.zeros(ndof) displacement[free] np.linalg.solve( stiffness[np.ix_(free, free)], load[free] ) return displacement, stiffness # 均匀材料先用解析位移场校核。 uniform_u, _ solve(np.full(len(triangles), E_MEAN)) exact_u np.column_stack(( PRESSURE * nodes[:, 0] / E_MEAN, -NU * PRESSURE * nodes[:, 1] / E_MEAN )).ravel() np.testing.assert_allclose( uniform_u, exact_u, rtol1e-8, atol1e-10 ) nominal_tip uniform_u[2 * right].mean() # 在三角形中心构建二维随机材料场。 distance np.linalg.norm( centers[:, None, :] - centers[None, :, :], axis2 ) sigma_log np.sqrt(np.log1p(CV**2)) covariance sigma_log**2 * np.exp( -distance / CORRELATION_LENGTH ) eigenvalues, eigenvectors np.linalg.eigh(covariance) eigenvalues np.maximum(eigenvalues, 0.0) rng np.random.default_rng(2042) normal_samples rng.standard_normal( (N_SAMPLES, len(triangles)) ) log_modulus ( np.log(E_MEAN) - 0.5 * sigma_log**2 normal_samples (eigenvectors * np.sqrt(eigenvalues)).T ) modulus_samples np.exp(log_modulus) assert np.all(modulus_samples 0) tips np.empty(N_SAMPLES) for sample_no, modulus in enumerate(modulus_samples): u, stiffness solve(modulus) tips[sample_no] u[2 * right].mean() if sample_no 0: reaction stiffness u - load np.testing.assert_allclose( reaction[2 * left].sum(), -PRESSURE * H * THICKNESS, atol1e-7 ) first_right_range np.ptp(u[2 * right]) p_exceed np.mean(tips LIMIT) probability_se np.sqrt( p_exceed * (1 - p_exceed) / N_SAMPLES ) print( f节点/三角形/自由度: f{len(nodes)}/{len(triangles)}/{ndof} ) print(f右边界合力: {load[0::2].sum():.3f} N) print(f均匀材料右边平均位移: {nominal_tip:.6f} mm) print( f首个随机样本右边位移极差: f{first_right_range:.6f} mm ) print(fMonte Carlo 样本数: {N_SAMPLES}) print( f随机材料右边平均位移的样本均值: f{tips.mean():.6f} mm ) print( f随机材料右边平均位移的样本标准差: f{tips.std(ddof1):.6f} mm ) print( fP(右边平均位移 {LIMIT:.3f} mm): f{p_exceed:.2%} ) print( f上述概率的 Monte Carlo 标准误: f{probability_se:.2%} ) print( 检查通过均匀材料解析解及首样本整体力平衡。 ) out Path(__file__).parent mesh mtri.Triangulation( nodes[:, 0], nodes[:, 1], triangles ) fig, ax plt.subplots(figsize(9, 3.3), dpi160) field ax.tripcolor( mesh, facecolorsmodulus_samples[0] / 1000, edgecolors#64748b, linewidth0.6, cmapviridis ) fig.colorbar( field, axax, labelElastic modulus E (GPa) ) ax.set( xlabelx (mm), ylabely (mm), titleOne realization of the element-wise random material field ) ax.set_aspect(equal) fig.tight_layout() fig.savefig(out / lesson42_material_field.png) plt.close(fig) fig, ax plt.subplots(figsize(9, 5), dpi160) ax.hist(tips, bins35, color#287d9d, alpha0.85) ax.axvline( nominal_tip, color#183b56, linestyle--, linewidth2, labelfUniform-material result: {nominal_tip:.4f} mm ) ax.axvline( LIMIT, color#c2410c, linewidth2, labelfIllustrative limit: {LIMIT:.3f} mm ) ax.set( xlabelMean right-edge displacement (mm), ylabelNumber of samples, titleMonte Carlo response of the 2D tensile plate ) ax.grid(axisy, alpha0.18) ax.legend(frameonFalse) fig.tight_layout() fig.savefig(out / lesson42_displacement_histogram.png) plt.close(fig) print( 图片已保存lesson42_material_field.png、 lesson42_displacement_histogram.png )实际运行输出节点/三角形/自由度: 28/36/56 右边界合力: 2000.000 N 均匀材料右边平均位移: 0.047619 mm 首个随机样本右边位移极差: 0.002392 mm Monte Carlo 样本数: 500 随机材料右边平均位移的样本均值: 0.047999 mm 随机材料右边平均位移的样本标准差: 0.002855 mm P(右边平均位移 0.050 mm): 25.00% 上述概率的 Monte Carlo 标准误: 1.94% 检查通过均匀材料解析解及首样本整体力平衡。 图片已保存lesson42_material_field.png、lesson42_displacement_histogram.png不同 NumPy 版本可能使随机统计值的末几位略有差别。六、先看材料场再看结构响应下图是一次抽样得到的单元弹性模量。颜色只表示各三角形单元在这次抽样中的材料参数不是应力图也不能凭颜色直接断定哪里已经失效。把500次求解得到的右边界平均位移放在一起就形成下图的分布。虚线是均匀材料结果橙线是本课为了练习概率统计而设的位移限值。这次样本中有 \(25\%\) 的位移超过 \(0.050\ \mathrm{mm}\)即500次中约125次。但估计值的 Monte Carlo 标准误约为1.94个百分点如果工程上关心的是很小的超限概率500次抽样远远不够。更重要的是这个概率完全依赖于本课假设的材料分布、变异系数和相关长度不能当成真实构件的失效概率。首个随机样本中右边各节点的水平位移极差为 \(0.002392\ \mathrm{mm}\)。这提醒我们二维模型的右边界不再必然像一维杆端点那样只对应一个位移值。因此在统计前必须先说清楚关注的是平均位移、最大位移还是某个特定位置的位移。七、这段程序检查了什么还没检查什么已完成的检查有两项当所有单元取相同弹性模量时数值解与上述解析位移场一致。随机材料首个样本的左端水平反力与右端 \(2000\ \mathrm N\) 拉力平衡。但这不等于模型已经达到工程使用精度。6×3矩形、36个三角形的网格主要用于教学500次抽样也只够展示方法。进一步使用前至少应分别改变网格密度、单元尺度与相关长度的比例、抽样次数并检查所关心响应量是否稳定。尤其要注意网格细化后不能把“重新抽到一组不同随机材料”造成的差异全部误判为有限元网格误差。八、如何迁移到 FEniCSx本节最值得迁移的不是 NumPy 的矩阵循环而是清晰的数据对应关系\[ \text{网格单元编号} \longleftrightarrow \text{随机场样本中的 }E_e \longleftrightarrow \text{弱形式中的材料系数} \]在 FEniCSx 中逐单元常值材料参数可用DG0零次间断函数空间表示官方教学示例说明了这种空间的单元值与自由度之间的对应方式。FEniCSx 多材料子区域教程 然后对每个随机场样本更新材料系数、求解弹性问题、提取同一个响应量。二维或三维线弹性弱形式及应力定义可参考 FEniCSx 线弹性教程。具体实现时还须核对所用版本的 API并行计算下尤其要处理好单元编号与本地自由度的映射。九、课后练习将关注量从“右边界平均位移”改成“右边界最大位移”。同一个 \(0.050\ \mathrm{mm}\) 限值下超限比例会怎样变化先预测再修改程序验证。将CORRELATION_LENGTH分别设为10 mm和60 mm对比位移标准差。注意每次都要保持其他假设不变。增加N_SAMPLES观察概率标准误再增加NX、NY观察网格变化。尝试分别记录这两类变化不要把它们合并为一次实验。本节完成了从二维随机材料场到二维有限元响应分布的闭环。下一节可以继续研究随机有限元的收敛与可信度多少单元、多少随机场模态、多少抽样次数才足以支撑一个具体结论 
长河君 · 2026/10/7 16:30:32 · 阅读 1.2k 【零基础学智能仿真-42】二维随机有限元实战——把随机材料场赋给网格单元
课程摘要上一节生成了沿杆长变化的随机弹性模量场本节将这一方法扩展到二维。我们以受拉矩形薄板为例在三角形单元上生成空间相关的弹性模量反复组装刚度矩阵并求解统计右边界位移。课程还通过均匀材料解析解和整体力平衡检查程序说明随机有限元中“生成样本、映射网格、验证求解、解释概率”缺一不可。一、从一维求和进入真正的二维求解第四十一节的拉杆有一个便利之处每个单元承受同样的轴力。因此端部位移可以写成各单元伸长量之和。二维薄板则不同。某块区域变软后周围区域可能分担更多载荷位移和应力必须通过整体刚度矩阵共同求解。我们要完成的计算链条是\[ \boxed{ \text{生成二维随机材料场} \longrightarrow \text{赋值给三角形单元} \longrightarrow K^{(s)}u^{(s)}f \longrightarrow \text{统计结构响应} } \]上标 \(s\) 表示第 \(s\) 次材料抽样。载荷和几何保持不变每次改变的是单元弹性模量因而刚度矩阵也随之改变。本节沿用第二十三节介绍过的 NumPy 常应变三角形单元CST思路让新增的“随机场—单元映射”过程容易核查。**本节代码是二维 NumPy 教学求解器并非声称已经在 FEniCSx 中运行。**文末再说明如何迁移到 FEniCSx。二、案例随机材料矩形薄板建立长 \(100\ \mathrm{mm}\)、高 \(20\ \mathrm{mm}\)、厚 \(1\ \mathrm{mm}\) 的薄板采用平面应力假设左边界约束水平位移左下角额外约束竖向位移以消除整体刚体运动。右边界施加 \(100\ \mathrm{MPa}\) 的均布水平拉应力对应合力 \(2000\ \mathrm N\)。泊松比固定为 \(\nu0.3\)。弹性模量平均值设为 \(210000\ \mathrm{MPa}\)变异系数设为 \(10\%\)。材料对数场的相关长度设为 \(25\ \mathrm{mm}\)。教学网格由6×3个矩形组成每个矩形剖分成两个三角形共36个三角形单元。每个三角形取一个弹性模量值。我们关心的输出不是某个任意节点的位移而是右边界各节点水平位移的平均值\[ Q^{(s)} \frac{1}{n_R} \sum_{a\in\Gamma_R}u_x^{(s)}(a) \]这样可以得到一个明确的、可重复统计的响应量。请注意这并不表示右边界所有节点的位移都相同。三、先建立确定性“标准答案”在均匀材料、单轴均布拉伸的理想条件下平面应力解析位移场为\[ u_x(x,y)\frac{p}{E}x,\qquad u_y(x,y)-\nu\frac{p}{E}y \]其中 \(p100\ \mathrm{MPa}\)。所以右端水平位移应为\[ u_x(L,y) \frac{pL}{E} \frac{100\times100}{210000} 0.047619\ \mathrm{mm} \]这一步很重要如果程序连均匀材料都算不对随后运行500次随机求解也只会得到500个值得怀疑的数字。本节代码会先用这个解析场检查单元刚度、载荷与边界约束。四、把空间随机场映射到二维单元对第 \(e\) 个三角形取其中心坐标 \(\boldsymbol c_e\)。与上一节类似我们为对数弹性模量设置距离相关的协方差\[ C_{ef} \sigma_G^2 \exp\left( -\frac{\|\boldsymbol c_e-\boldsymbol c_f\|}{\ell} \right) \]这一次距离是二维欧氏距离而不再只是沿杆长的距离。对协方差矩阵做特征值分解再生成各单元的随机弹性模量 \(E_e^{(s)}0\)。本课保留全部离散模态避免在讲“映射与求解”的同时混入 KL 截断误差。平面应力单元刚度为\[ K_e^{(s)} tA_e B_e^{\mathsf T}D(E_e^{(s)},\nu)B_e \]其中 \(t\) 为厚度\(A_e\) 为三角形面积\(B_e\) 是应变—位移矩阵。把每个 \(K_e^{(s)}\) 放到对应的整体自由度位置便得到本次样本的整体刚度矩阵。五、完整可运行代码第四十二节二维平面应力 CST 有限元与随机弹性模量场。 from pathlib import Path import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np # 单位N、mm、MPa右侧均布拉应力为 100 MPa。 L, H, THICKNESS, PRESSURE 100.0, 20.0, 1.0, 100.0 E_MEAN, NU, CV, CORRELATION_LENGTH 210000.0, 0.3, 0.10, 25.0 NX, NY, N_SAMPLES, LIMIT 6, 3, 500, 0.050 x, y np.meshgrid( np.linspace(0, L, NX 1), np.linspace(0, H, NY 1) ) nodes np.column_stack((x.ravel(), y.ravel())) triangles [] for j in range(NY): for i in range(NX): a j * (NX 1) i b, c, d a 1, a NX 1, a NX 2 triangles.extend(((a, b, d), (a, d, c))) triangles np.asarray(triangles, dtypeint) centers nodes[triangles].mean(axis1) ndof 2 * len(nodes) # 平面应力的单位弹性模量本构矩阵。 d_unit np.array([ [1, NU, 0], [NU, 1, 0], [0, 0, (1 - NU) / 2] ]) / (1 - NU**2) element_dofs, element_k_unit [], [] for tri in triangles: xy nodes[tri] xx, yy xy[:, 0], xy[:, 1] area np.linalg.det( np.column_stack((np.ones(3), xx, yy)) ) / 2 assert area 0 b np.roll(yy, -1) - np.roll(yy, 1) c np.roll(xx, 1) - np.roll(xx, -1) bmat np.zeros((3, 6)) bmat[0, 0::2], bmat[1, 1::2] b, c bmat[2, 0::2], bmat[2, 1::2] c, b bmat / 2 * area element_dofs.append(np.array([ 2 * node k for node in tri for k in (0, 1) ])) element_k_unit.append( THICKNESS * area * bmat.T d_unit bmat ) # 右侧边均匀施加拉应力。 load np.zeros(ndof) for j in range(NY): lower j * (NX 1) NX upper (j 1) * (NX 1) NX edge_force PRESSURE * THICKNESS * (H / NY) / 2 load[2 * lower] edge_force load[2 * upper] edge_force left np.flatnonzero(np.isclose(nodes[:, 0], 0.0)) right np.flatnonzero(np.isclose(nodes[:, 0], L)) fixed np.r_[2 * left, 1] free np.setdiff1d(np.arange(ndof), fixed) def solve(modulus_by_element): stiffness np.zeros((ndof, ndof)) for modulus, dofs, k_unit in zip( modulus_by_element, element_dofs, element_k_unit ): stiffness[np.ix_(dofs, dofs)] modulus * k_unit displacement np.zeros(ndof) displacement[free] np.linalg.solve( stiffness[np.ix_(free, free)], load[free] ) return displacement, stiffness # 均匀材料先用解析位移场校核。 uniform_u, _ solve(np.full(len(triangles), E_MEAN)) exact_u np.column_stack(( PRESSURE * nodes[:, 0] / E_MEAN, -NU * PRESSURE * nodes[:, 1] / E_MEAN )).ravel() np.testing.assert_allclose( uniform_u, exact_u, rtol1e-8, atol1e-10 ) nominal_tip uniform_u[2 * right].mean() # 在三角形中心构建二维随机材料场。 distance np.linalg.norm( centers[:, None, :] - centers[None, :, :], axis2 ) sigma_log np.sqrt(np.log1p(CV**2)) covariance sigma_log**2 * np.exp( -distance / CORRELATION_LENGTH ) eigenvalues, eigenvectors np.linalg.eigh(covariance) eigenvalues np.maximum(eigenvalues, 0.0) rng np.random.default_rng(2042) normal_samples rng.standard_normal( (N_SAMPLES, len(triangles)) ) log_modulus ( np.log(E_MEAN) - 0.5 * sigma_log**2 normal_samples (eigenvectors * np.sqrt(eigenvalues)).T ) modulus_samples np.exp(log_modulus) assert np.all(modulus_samples 0) tips np.empty(N_SAMPLES) for sample_no, modulus in enumerate(modulus_samples): u, stiffness solve(modulus) tips[sample_no] u[2 * right].mean() if sample_no 0: reaction stiffness u - load np.testing.assert_allclose( reaction[2 * left].sum(), -PRESSURE * H * THICKNESS, atol1e-7 ) first_right_range np.ptp(u[2 * right]) p_exceed np.mean(tips LIMIT) probability_se np.sqrt( p_exceed * (1 - p_exceed) / N_SAMPLES ) print( f节点/三角形/自由度: f{len(nodes)}/{len(triangles)}/{ndof} ) print(f右边界合力: {load[0::2].sum():.3f} N) print(f均匀材料右边平均位移: {nominal_tip:.6f} mm) print( f首个随机样本右边位移极差: f{first_right_range:.6f} mm ) print(fMonte Carlo 样本数: {N_SAMPLES}) print( f随机材料右边平均位移的样本均值: f{tips.mean():.6f} mm ) print( f随机材料右边平均位移的样本标准差: f{tips.std(ddof1):.6f} mm ) print( fP(右边平均位移 {LIMIT:.3f} mm): f{p_exceed:.2%} ) print( f上述概率的 Monte Carlo 标准误: f{probability_se:.2%} ) print( 检查通过均匀材料解析解及首样本整体力平衡。 ) out Path(__file__).parent mesh mtri.Triangulation( nodes[:, 0], nodes[:, 1], triangles ) fig, ax plt.subplots(figsize(9, 3.3), dpi160) field ax.tripcolor( mesh, facecolorsmodulus_samples[0] / 1000, edgecolors#64748b, linewidth0.6, cmapviridis ) fig.colorbar( field, axax, labelElastic modulus E (GPa) ) ax.set( xlabelx (mm), ylabely (mm), titleOne realization of the element-wise random material field ) ax.set_aspect(equal) fig.tight_layout() fig.savefig(out / lesson42_material_field.png) plt.close(fig) fig, ax plt.subplots(figsize(9, 5), dpi160) ax.hist(tips, bins35, color#287d9d, alpha0.85) ax.axvline( nominal_tip, color#183b56, linestyle--, linewidth2, labelfUniform-material result: {nominal_tip:.4f} mm ) ax.axvline( LIMIT, color#c2410c, linewidth2, labelfIllustrative limit: {LIMIT:.3f} mm ) ax.set( xlabelMean right-edge displacement (mm), ylabelNumber of samples, titleMonte Carlo response of the 2D tensile plate ) ax.grid(axisy, alpha0.18) ax.legend(frameonFalse) fig.tight_layout() fig.savefig(out / lesson42_displacement_histogram.png) plt.close(fig) print( 图片已保存lesson42_material_field.png、 lesson42_displacement_histogram.png )实际运行输出节点/三角形/自由度: 28/36/56 右边界合力: 2000.000 N 均匀材料右边平均位移: 0.047619 mm 首个随机样本右边位移极差: 0.002392 mm Monte Carlo 样本数: 500 随机材料右边平均位移的样本均值: 0.047999 mm 随机材料右边平均位移的样本标准差: 0.002855 mm P(右边平均位移 0.050 mm): 25.00% 上述概率的 Monte Carlo 标准误: 1.94% 检查通过均匀材料解析解及首样本整体力平衡。 图片已保存lesson42_material_field.png、lesson42_displacement_histogram.png不同 NumPy 版本可能使随机统计值的末几位略有差别。六、先看材料场再看结构响应下图是一次抽样得到的单元弹性模量。颜色只表示各三角形单元在这次抽样中的材料参数不是应力图也不能凭颜色直接断定哪里已经失效。把500次求解得到的右边界平均位移放在一起就形成下图的分布。虚线是均匀材料结果橙线是本课为了练习概率统计而设的位移限值。这次样本中有 \(25\%\) 的位移超过 \(0.050\ \mathrm{mm}\)即500次中约125次。但估计值的 Monte Carlo 标准误约为1.94个百分点如果工程上关心的是很小的超限概率500次抽样远远不够。更重要的是这个概率完全依赖于本课假设的材料分布、变异系数和相关长度不能当成真实构件的失效概率。首个随机样本中右边各节点的水平位移极差为 \(0.002392\ \mathrm{mm}\)。这提醒我们二维模型的右边界不再必然像一维杆端点那样只对应一个位移值。因此在统计前必须先说清楚关注的是平均位移、最大位移还是某个特定位置的位移。七、这段程序检查了什么还没检查什么已完成的检查有两项当所有单元取相同弹性模量时数值解与上述解析位移场一致。随机材料首个样本的左端水平反力与右端 \(2000\ \mathrm N\) 拉力平衡。但这不等于模型已经达到工程使用精度。6×3矩形、36个三角形的网格主要用于教学500次抽样也只够展示方法。进一步使用前至少应分别改变网格密度、单元尺度与相关长度的比例、抽样次数并检查所关心响应量是否稳定。尤其要注意网格细化后不能把“重新抽到一组不同随机材料”造成的差异全部误判为有限元网格误差。八、如何迁移到 FEniCSx本节最值得迁移的不是 NumPy 的矩阵循环而是清晰的数据对应关系\[ \text{网格单元编号} \longleftrightarrow \text{随机场样本中的 }E_e \longleftrightarrow \text{弱形式中的材料系数} \]在 FEniCSx 中逐单元常值材料参数可用DG0零次间断函数空间表示官方教学示例说明了这种空间的单元值与自由度之间的对应方式。FEniCSx 多材料子区域教程 然后对每个随机场样本更新材料系数、求解弹性问题、提取同一个响应量。二维或三维线弹性弱形式及应力定义可参考 FEniCSx 线弹性教程。具体实现时还须核对所用版本的 API并行计算下尤其要处理好单元编号与本地自由度的映射。九、课后练习将关注量从“右边界平均位移”改成“右边界最大位移”。同一个 \(0.050\ \mathrm{mm}\) 限值下超限比例会怎样变化先预测再修改程序验证。将CORRELATION_LENGTH分别设为10 mm和60 mm对比位移标准差。注意每次都要保持其他假设不变。增加N_SAMPLES观察概率标准误再增加NX、NY观察网格变化。尝试分别记录这两类变化不要把它们合并为一次实验。本节完成了从二维随机材料场到二维有限元响应分布的闭环。下一节可以继续研究随机有限元的收敛与可信度多少单元、多少随机场模态、多少抽样次数才足以支撑一个具体结论