SolidWorks_动画模拟与仿真15_仿真结果可视化

仿真结果可视化：图解、矢量图与动态高亮的工程实践

本文深入探讨工程仿真（FEA/CFD）结果可视化的核心技术，从数据映射到动态交互，完整呈现应力、位移等物理场的图形化表达方案，并附赠可直接运行的Python代码示例。

摘要

仿真计算产生海量节点数据，但唯有通过高效的可视化才能转化为工程洞察。本文将系统讲解应力/位移场的伪彩图（云图）、矢量箭头图、变形放大图及动态高亮交互的实现原理与代码实践，涵盖从有限元数据格式到交互式Web可视化的完整链路，帮助工程师与开发者构建自己的仿真后处理工具。

1. 引言：为什么可视化是仿真的“最后一公里”

任何有限元分析（FEA）或计算流体力学（CFD）求解器输出的都是离散的数值矩阵——节点坐标、单元连接、应力张量、位移矢量。这些数据本身毫无直观性可言。仿真结果可视化的本质，是将高维数值映射为人类视觉可感知的图形元素。

一个典型的工程场景：你完成了一个悬臂梁的受力分析，求解器告诉你最大应力为235MPa，位于固定端下缘。但：

• 应力集中区域的具体形状是什么？
• 位移的渐变趋势是线性还是非线性？
• 危险截面到底在哪个精确位置？

这些问题无法从数字中直接获得答案，必须依赖可视化技术。优秀的可视化不仅能“看见”结果，更能揭示物理规律，甚至发现求解错误（如网格畸变导致的应力奇异）。

2. 基础准备：仿真数据的三层结构

在动手绘图前，必须理解有限元结果的数据组织方式。几乎所有商业软件（ANSYS、Abaqus）和开源框架（FEniCS、OpenFOAM）都遵循以下三层结构：

2.1 节点坐标（Nodal Coordinates）

# 节点表：每行 [节点ID, x, y, z]nodes=np.array([[0,0.0,0.0],# 节点1[1,0.1,0.0],# 节点2[2,0.1,0.05],# 节点3...])

2.2 单元连接（Element Connectivity）

# 单元表：每行 [单元ID, 节点1, 节点2, 节点3, ...]elements=np.array([[0,0,1,2],# 三角形单元（2D）[1,2,3,0],# 另一个三角形...])

2.3 物理场数据（Field Data）

# 节点位移：每行 [节点ID, ux, uy, uz]displacements=np.array([[0,0.001,-0.002,0.0],[1,0.002,-0.003,0.0],...])# 单元应力：每行 [单元ID, sx, sy, sz, sxy, syz, sxz]stresses=np.array([[0,120.5,80.2,0.0,35.1,0.0,0.0],...])

关键区别：位移是节点变量（连续），应力通常是单元变量（分段常数），可视化时需要特殊处理。

3. 核心技法一：应力/位移云图（伪彩图）

云图是仿真可视化的“主力军”，通过颜色映射表达标量场的分布。

3.1 实现原理

1. 插值：将单元应力转换到节点（面积加权平均）
2. 归一化：将物理值映射到[0,1]区间
3. 颜色映射：应用colormap（如’jet’、‘viridis’）

3.2 完整代码示例（使用Matplotlib）

importnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.triasmtrifrommatplotlib.colorsimportNormalizedefplot_stress_contour(nodes,elements,stress_values,title="应力云图"):""" 绘制2D三角形网格的应力云图 :param nodes: (N,2) 节点坐标 :param elements: (M,3) 三角形单元连接 :param stress_values: (M,) 单元应力值（如von Mises） """# 1. 将单元应力插值到节点（简单面积加权平均）node_stress=np.zeros(len(nodes))node_count=np.zeros(len(nodes))foreleminelements:area=0.5*abs(np.cross(nodes[elem[1]]-nodes[elem[0]],nodes[elem[2]]-nodes[elem[0]]))fornidinelem:node_stress[nid]+=stress_values[elem[0]]*area node_count[nid]+=area node_stress/=(node_count+1e-12)# 避免除零# 2. 创建三角剖分对象triang=mtri.Triangulation(nodes[:,0],nodes[:,1],elements)# 3. 绘图fig,ax=plt.subplots(figsize=(10,8))tcf=ax.tripcolor(triang,node_stress,shading='gouraud',# 平滑着色cmap='jet',norm=Normalize(vmin=np.min(node_stress),vmax=np.max(node_stress)))# 4. 添加网格线（可选）ax.triplot(triang,'k-',lw=0.3,alpha=0.5)# 5. 装饰plt.colorbar(tcf,ax=ax,label='Stress (MPa)')ax.set_aspect('equal')ax.set_title(title)ax.set_xlabel('X (m)')ax.set_ylabel('Y (m)')plt.tight_layout()returnfig,ax# 示例：生成一个带孔平板模型if__name__=="__main__":# 生成简单网格（实际工程中从求解器读取）frommesh_generatorimportgenerate_plate_with_hole nodes,elements=generate_plate_with_hole()# 模拟应力结果（真实数据来自求解器）stress=100*np.random.rand(len(elements))plot_stress_contour(nodes,elements,stress,"带孔平板 von Mises 应力分布")plt.show()

3.3 进阶技巧

• 对数色标：当应力跨越多个数量级时（如1~10^6），使用LogNorm
• 等值线叠加：使用ax.contour(triang, node_stress, levels=10)叠加等值线
• 透明处理：对低于阈值的区域设置alpha=0.3，突出危险区

4. 核心技法二：矢量箭头图（位移/流动方向）

应力是张量，位移是矢量。矢量场的可视化需要方向+大小双重信息。

4.1 箭头图实现

defplot_displacement_vectors(nodes,displacements,scale=1000,step=1):""" 绘制位移矢量场 :param displacements: (N,2) 每个节点的位移矢量 :param scale: 箭头放大倍数（位移通常很小） :param step: 每隔几个节点绘制一个箭头（避免过密） """fig,ax=plt.subplots(figsize=(10,8))# 提取节点位置（原始坐标）x=nodes[::step,0]y=nodes[::step,1]# 提取位移分量（放大显示）ux=displacements[::step,0]*scale uy=displacements[::step,1]*scale# 绘制箭头q=ax.quiver(x,y,ux,uy,angles='xy',scale_units='xy',scale=1,color='b',width=0.002,alpha=0.7)# 添加颜色映射（表示位移大小）magnitude=np.sqrt(ux**2+uy**2)q.set_array(magnitude)plt.colorbar(q,label='Displacement Magnitude (×scale)')# 绘制原始网格（浅色背景）ax.triplot(mtri.Triangulation(nodes[:,0],nodes[:,1],elements),'k-',lw=0.2,alpha=0.2)ax.set_aspect('equal')ax.set_title('位移矢量分布')ax.set_xlabel('X (m)')ax.set_ylabel('Y (m)')returnfig,ax

4.2 流线图（CFD专用）

对于流体仿真，箭头图过于杂乱，应使用流线：

# 使用matplotlib的streamplot（需要规则网格）defplot_streamlines(velocity_field,x_grid,y_grid):fig,ax=plt.subplots()strm=ax.streamplot(x_grid,y_grid,velocity_field[:,:,0],velocity_field[:,:,1],density=1.5,color=velocity_field[:,:,2],cmap='coolwarm',linewidth=1.5)plt.colorbar(strm.lines,label='速度大小 (m/s)')returnfig,ax

5. 核心技法三：变形放大图与动画

结构仿真的位移往往远小于几何尺寸（如0.1mm vs 100mm），直接绘制无法观察。变形放大是必备技术。

5.1 静态变形图

defplot_deformed_shape(nodes,elements,displacements,scale_factor=100):""" 绘制变形前后对比图 :param scale_factor: 位移放大倍数 """# 计算变形后坐标deformed_nodes=nodes+displacements*scale_factor fig,(ax1,ax2)=plt.subplots(1,2,figsize=(14,6))# 原始形状ax1.triplot(mtri.Triangulation(nodes[:,0],nodes[:,1],elements),'b-',lw=1,label='原始形状')ax1.set_title('未变形')ax1.set_aspect('equal')# 变形形状ax2.triplot(mtri.Triangulation(deformed_nodes[:,0],deformed_nodes[:,1],elements),'r-',lw=1.5,label=f'变形 (×{scale_factor})')ax2.set_title(f'变形放大{scale_factor}倍')ax2.set_aspect('equal')# 可选：叠加云图# ...（结合第3节代码）returnfig

5.2 动态高亮：使用Matplotlib Animation

importmatplotlib.animationasanimationdefanimate_vibration(nodes,elements,mode_shape,frequency=5,cycles=3):""" 动态显示振型（模态分析结果） :param mode_shape: (N,2) 模态位移 :param frequency: 动画频率（Hz） :param cycles: 显示多少个周期 """fig,ax=plt.subplots(figsize=(10,8))triang=mtri.Triangulation(nodes[:,0],nodes[:,1],elements)# 初始化绘图对象tri_plot=ax.tripcolor(triang,np.zeros(len(nodes)),cmap='RdYlBu_r',vmin=-1,vmax=1)ax.set_aspect('equal')# 时间参数fps=30total_frames=int(fps*cycles/frequency)t=np.linspace(0,cycles/frequency,total_frames)defupdate(frame):# 当前时刻的位移 = 模态位移 × sin(ωt)phase=np.sin(2*np.pi*frequency*t[frame])current_disp=mode_shape*phase# 更新节点位置current_nodes=nodes+current_disp*0.05# 放大系数# 更新三角剖分new_triang=mtri.Triangulation(current_nodes[:,0],current_nodes[:,1],elements)# 更新颜色（按当前位移大小着色）magnitude=np.linalg.norm(current_disp,axis=1)tri_plot.set_array(magnitude)# 更新网格ax.clear()ax.tripcolor(new_triang,magnitude,cmap='RdYlBu_r',vmin=0,vmax=1)ax.set_title(f'时间:{t[frame]:.3f}s')ax.set_aspect('equal')return[tri_plot]anim=animation.FuncAnimation(fig,update,frames=total_frames,interval=1000/fps,blit=False)returnanim# 保存动画# anim.save('vibration.gif', writer='pillow', fps=fps)

6. 核心技法四：动态高亮与交互式查询

静态图片无法满足工程分析的深度需求，我们需要交互式探索。

6.1 使用Plotly实现鼠标悬停查询

importplotly.graph_objectsasgodefinteractive_stress_plot(nodes,elements,stress):"""创建可交互的应力云图"""# 构建三角形网格数据tri_points=[]foreleminelements:fornidinelem:tri_points.append([nodes[nid,0],nodes[nid,1],stress[nid]])fig=go.Figure(data=[go.Mesh3d(x=nodes[:,0],y=nodes[:,1],z=np.zeros(len(nodes)),i=elements[:,0],j=elements[:,1],k=elements[:,2],intensity=stress,colorscale='Jet',showscale=True,hovertemplate='<b>X</b>: %{x:.3f}<br>'+'<b>Y</b>: %{y:.3f}<br>'+'<b>应力</b>: %{intensity:.2f} MPa<extra></extra>')])fig.update_layout(title='交互式应力分布 (悬停查看数值)',scene=dict(aspectmode='data'),template='plotly_white')returnfig

6.2 动态高亮：点击单元显示详细信息

fromipywidgetsimportinteract,widgetsdefhighlight_extreme_elements(nodes,elements,stress,threshold=0.9):""" 高亮超过阈值的危险单元 :param threshold: 应力阈值比例（0~1） """max_stress=np.max(stress)danger_mask=stress>threshold*max_stress fig,ax=plt.subplots(figsize=(12,8))# 绘制全部单元（灰色）triang=mtri.Triangulation(nodes[:,0],nodes[:,1],elements)ax.tripcolor(triang,stress,cmap='gray',alpha=0.3)# 高亮危险单元（红色）danger_elements=elements[danger_mask]iflen(danger_elements)>0:danger_triang=mtri.Triangulation(nodes[:,0],nodes[:,1],danger_elements)ax.tripcolor(danger_triang,stress[danger_mask],cmap='autumn_r',edgecolors='k',linewidth=1)# 标记最大应力点max_idx=np.argmax(stress)max_node=elements[max_idx][0]ax.plot(nodes[max_node,0],nodes[max_node,1],'r*',markersize=15,label=f'最大应力:{max_stress:.1f}MPa')ax.legend()ax.set_aspect('equal')ax.set_title(f'危险区域高亮 (阈值:{threshold*100:.0f}% 最大应力)')returnfig

7. 工程实践：完整案例——悬臂梁受力分析可视化

现在将所有技术整合，完成一个完整的工程案例。

7.1 问题定义

• 悬臂梁尺寸：1m × 0.2m
• 材料：钢（E=210GPa, ν=0.3）
• 载荷：自由端施加垂直向下100kN

7.2 完整流程代码

importnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.triasmtrifromscipy.sparseimportlil_matrixfromscipy.sparse.linalgimportspsolvedefcantilever_beam_analysis():""" 2D悬臂梁有限元分析（简化版，仅演示可视化） 实际工程请使用专业FEA软件 """# ---------- 1. 网格生成 ----------nx,ny=20,5# 网格密度x=np.linspace(0,1.0,nx+1)y=np.linspace(0,0.2,ny+1)nodes=[]forjinrange(ny+1):foriinrange(nx+1):nodes.append([x[i],y[j]])nodes=np.array(nodes)# 生成三角形单元elements=[]forjinrange(ny):foriinrange(nx):n0=j*(nx+1)+i n1=n0+1n2=n0+(nx+1)n3=n2+1elements.append([n0,n1,n2])# 三角形1elements.append([n1