Python读取comsol仿真导出数据并绘图

当我们使用comsol，ansys等仿真工具进行仿真后，难免需要对仿真结果进行导出并进一步处理分析。

今天小姜以comsol的一个简单磁场仿真为例，详细介绍如何对comsol仿真数据进行导出并在python中处理，绘制云图：

comsol数据导出

以一个简单的矢量磁位预测为例，在导出模块右键新建“数据”项

选择相应导出内容和导出文件位置后导出数据为txt文本格式。博主选择导出数据为区域的矢量磁位Az，位置为D盘

导出后可以看到txt文件格式如下，前几行为描述文件信息的注释行。

python读取文件

得到如上格式的txt文件后，首先我们在读取时需要先跳过前面的注释行

    with open(path, 'r') as f:
        lines = f.readlines()
        # 过滤掉以%开头的行，并去掉空行
        data_lines = [line for line in lines if not line.startswith('%')  and line.strip()]
        # 转换为NumPy数组
        data = np.array([[float(num) for num in line.split()] for line in data_lines])

将坐标和对应Az读取到numpy格式的data变量中后还需要整理为x坐标y坐标，Az值三个变量，以便云图绘制

# 提取各列数据
    x = data[:, 0]
    y = data[:, 1]
    mf_az = data[:, 2]

由于导出数据可能存在计算的奇异点等错误数据，所以还需要排除掉值为NaN的点

# 去掉值为NaN的点
    valid_indices = ~np.isnan(mf_az)
    x = x[valid_indices]
    y = y[valid_indices]
    mf_az = mf_az[valid_indices]

当然这样的x,y和az并不能直接用matplotlib进行云图绘制，因为这些坐标点并不是规则的网格而是仿真软件的网格剖分得到的不规则点。我们使用scipy对这些网格进行插值，转化为规则的网格点后才能绘制。

# 创建规则网格
    xi = np.linspace(min(x), max(x), 300)
    yi = np.linspace(min(y), max(y), 300)
    Xi, Yi = np.meshgrid(xi, yi)

    # 插值到规则网格(Cubic表示三次多项式插值)
    zi = griddata((x, y), mf_az, (Xi, Yi), method='cubic')

python绘制云图

获得了规则网格和对应值，使用imshow即可画热力图

# 绘制热力图
    plt.imshow(zi, extent=[min(x), max(x), min(y), max(y)], origin='lower', cmap='jet')
    plt.colorbar(label='mf.AZ (Wb/m)')
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.title('Heatmap of mf.AZ')
    plt.show()

python进一步分析数据

我们导出的数据为矢量磁位，根据磁有关的基础知识，A的旋度为磁场，二维可以表示为

$$H_x=\frac{\partial A_z}{\partial y}$$$$H_y=-\frac{\partial A_z}{\partial x}$$

使用numpy自带的求导功能求导并绘制矢量图：

    skip = 10 #避免箭头过密，设置绘制间隔
    #计算导数
    Zy, Zx = np.gradient(zi, yi, xi)
    # 绘制磁场矢量图
    plt.figure(figsize=(8, 6))
    quiver = plt.quiver(Xi[::skip, ::skip], Yi[::skip, ::skip], Zy[::skip, ::skip], -Zx[::skip, ::skip],
                         cmap='rainbow')
    plt.colorbar(quiver, label='Magnetic Field Magnitude')
    plt.title('Magnetic Field Vector Plot')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()

完整代码

def read_comsol(path):
    """

    :param path: 仿真导出数据文件的路径
    :return: 无输出
    """
    with open(path, 'r') as f:
        lines = f.readlines()
        # 过滤掉以%开头的行，并去掉空行
        data_lines = [line for line in lines if not line.startswith('%')  and line.strip()]
        # 转换为NumPy数组
        data = np.array([[float(num) for num in line.split()] for line in data_lines])

    # 提取各列数据
    x = data[:, 0]
    y = data[:, 1]
    mf_az = data[:, 2]

    # 去掉值为NaN的点
    valid_indices = ~np.isnan(mf_az)
    x = x[valid_indices]
    y = y[valid_indices]
    mf_az = mf_az[valid_indices]

    # 创建规则网格
    xi = np.linspace(min(x), max(x), 300)
    yi = np.linspace(min(y), max(y), 300)
    Xi, Yi = np.meshgrid(xi, yi)

    # 插值到规则网格
    zi = griddata((x, y), mf_az, (Xi, Yi), method='cubic')

    # 绘制热力图
    plt.imshow(zi, extent=[min(x), max(x), min(y), max(y)], origin='lower', cmap='jet')
    plt.colorbar(label='mf.AZ (Wb/m)')
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.title('Heatmap of mf.AZ')
    plt.show()

    skip = 10
    #计算导数
    Zy, Zx = np.gradient(zi, yi, xi)
    # 绘制磁场矢量图
    plt.figure(figsize=(8, 6))
    quiver = plt.quiver(Xi[::skip, ::skip], Yi[::skip, ::skip], Zy[::skip, ::skip], -Zx[::skip, ::skip])
    plt.colorbar(quiver, label='Magnetic Field Magnitude')
    plt.title('Magnetic Field Vector Plot')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.grid()
    plt.show()