Python实战：低周疲劳试验数据可视化与滞回环分析

1. 从数据文件到第一张图：快速上手

如果你手头有一份低周疲劳试验的原始数据，比如一个CSV文件，里面密密麻麻记录着时间、应力、应变，你的第一反应可能是：“这数据怎么看？” 别急，用Python把它变成直观的图形，其实就几行代码的事。

我处理过不少这类数据，第一步永远是先“看看它长什么样”。假设你的数据文件叫fatigue_data.csv，通常会有“time(s)”、“stress(MPA)”、“strain”这几列。直接用pandas读进来，然后用matplotlib画个最原始的应力-应变散点图或连线图，这是最直接的诊断。

import pandas as pd import matplotlib.pyplot as plt # 读取数据，注意列名可能有空格，最好先打印看看 data = pd.read_csv('fatigue_data.csv') print(data.head()) # 看看前几行，确认列名 print(data.columns) # 确认列名到底是什么 # 提取数据，列名根据实际情况调整，比如可能是 'strain '（带空格） time = data['time(s)'] stress = data['stress(MPA)'] strain = data['strain '] # 注意这里列名可能有个尾随空格 # 画最原始的应力-应变曲线 plt.figure(figsize=(10, 6)) plt.plot(strain, stress, 'b-', linewidth=0.5, alpha=0.7) plt.xlabel('Strain') plt.ylabel('Stress (MPa)') plt.title('Raw Stress-Strain Data from Fatigue Test') plt.grid(True, linestyle='--', alpha=0.5) plt.show()

运行这段代码，你大概率会得到一张“毛线团”一样的图。所有的循环加载卸载数据叠加在一起，线条密集交错，根本分不清单个循环，更别说观察滞回环的演变趋势了。这张图的价值在于告诉你数据是连续的，没有明显的断点或异常值，但它离我们想要的分析结果还差得远。很多新手到这一步就卡住了，觉得数据太乱无从下手。其实这正是低周疲劳数据的特点——它是在一个样本上反复加载卸载成百上千次记录下来的，我们需要做的，就是把这个“毛线团”按照时间周期，一根一根地拆解开来。

这里有个小坑我踩过：数据采集的频率可能不是绝对均匀的。虽然实验设定可能是每2秒记录一次，但实际设备可能会因为通讯、存储等原因有微小的时间漂移。所以，不要想当然地认为每个循环的数据点数量完全一样。在后续处理中，我们需要用时间列作为分组的依据，而不是简单地按数据点数量等分。另外，绘图时如果发现图形异常，比如应力值突然跳到零或者出现不合理的峰值，一定要回头检查原始数据，可能是传感器瞬时报错或数据记录有误，需要清洗或剔除。

2. 拆解“毛线团”：数据分段与周期识别

拿到“毛线团”图之后，下一步关键操作就是按加载周期把数据拆开。低周疲劳试验通常是应变控制或应力控制，以固定的周期（比如200秒）进行加载和卸载。我们的目标是把连续的时间-应力-应变数据，切割成一个个独立的循环。

最可靠的方法是利用时间列。假设我们知道一个完整的加载-卸载周期是 T=200 秒，那么第 k 个循环的数据，就满足时间 t 在[k*T, (k+1)*T]这个区间内。用numpy或pandas的布尔索引可以轻松实现这个筛选。

import numpy as np T = 200 # 已知的周期，单位秒 total_cycles = int(time.iloc[-1] // T) # 计算总循环次数 plt.figure(figsize=(12, 8)) colors = plt.cm.viridis(np.linspace(0, 1, total_cycles)) # 使用渐变色区分不同循环 for k in range(total_cycles): # 创建布尔掩码，筛选出第k个周期的数据 mask = (time >= k * T) & (time < (k+1) * T) strain_cycle = strain[mask] stress_cycle = stress[mask] # 绘制单个循环 plt.plot(strain_cycle, stress_cycle, color=colors[k], linewidth=1, label=f'Cycle {k+1}' if k < 5 else "") plt.xlabel('Strain') plt.ylabel('Stress (MPa)') plt.title('Individual Cycles (First 5 labeled)') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) plt.show()

运行这段代码，你会看到很多条颜色各异的曲线叠在一起，每一条代表一次完整的加载卸载过程，也就是一个滞回环。这时你应该能初步看出一些规律了：比如环的形状是否对称？环的面积（代表耗散能）是否在变化？最大应力点是否在逐渐下移或上移？不过，如果循环次数很多（比如100次），这张图仍然会显得非常拥挤，我们很难看清从第一次循环到最后一次循环，滞回环是如何一步步演变的。我们真正需要的，是一种能够清晰展示演变趋势的可视化方法，而不是把所有原始循环都堆上去。

这里就引出了核心问题：当循环次数很多时，如何既能减少绘图线条避免重叠，又能真实反映材料性能随循环次数的变化？直接每隔几个循环画一个，可能会错过重要特征；把所有循环画出来，又太乱。这就需要用到“曲线平均化”的策略。我试过几种方法，比如简单地对每N个循环的数据点取平均，但效果很差，因为不同循环的数据点在应变轴上并不对齐。后来发现，更合理的思路是：在时间（或相位）维度上进行对齐和平均。我们可以把每个200秒的周期，映射到一个标准化的“相位时间”上（比如0到200秒），然后对相同相位时间点的应力值进行平均，从而生成一条“代表”某几个循环的典型滞回环曲线。

3. 化繁为简：滞回环的“平均化”绘制技巧

这是本文最核心也最实用的一部分。直接绘制所有原始循环线条会糊成一团，而只画首尾几个循环又可能丢失中间过程信息。“平均化”绘制的目的，是在保留演变趋势的前提下，大幅减少绘制线条的数量，让图像清晰可读。

我常用的策略是“步进平均法”。假设原始有 n=100 个循环，我们只想画出 m=10 个有代表性的滞回环来展示演变过程。除了第一个循环（初始状态很重要，单独画）和最后一个循环（最终状态），中间的8个环，每个都不是来自单一循环，而是由附近多个循环平均合成而来。这样既能平滑掉单次循环的微小波动，又能显著减少线条数量。

关键是如何选择哪些循环参与平均。我设计了一种“滑动窗口”平均的方法，效果比较稳定。思路如下：我们希望画 m 个环，就把 n 个循环分成 m 组。第一组就是第一个循环（不平均），最后一组是最后一个循环（不平均），中间的 m-2 组，每一组都包含连续 k 个原始循环，对这 k 个循环在相同相位时间点的应力、应变值进行平均，生成一个新的代表性循环。k 的大小和分组方式需要计算，确保所有原始循环都被用到，且分组均匀。

def generate_averaging_scheme(total_cycles, desired_plots): """ 生成平均化方案。 total_cycles: 原始总循环数 n desired_plots: 想要绘制的代表性滞回环数量 m (包括首尾) 返回一个列表，列表长度=m，每个元素是一个整数列表，包含用于平均的原始循环编号。 """ n = total_cycles m = desired_plots # 第一个和最后一个循环单独处理 scheme = [[0]] # 第一个循环 # 中间需要平均的组数 groups_in_middle = m - 2 # 可用于分配的循环数 (去掉首尾) cycles_to_distribute = n - 2 # 计算每组大致包含的原始循环数 cycles_per_group = cycles_to_distribute // groups_in_middle remainder = cycles_to_distribute % groups_in_middle start_idx = 1 for i in range(groups_in_middle): # 确定当前组包含的循环数，处理余数 size = cycles_per_group + (1 if i < remainder else 0) end_idx = start_idx + size scheme.append(list(range(start_idx, end_idx))) start_idx = end_idx # 添加最后一个循环 scheme.append([n-1]) return scheme # 示例：100个循环，想画10个代表性环 plot_scheme = generate_averaging_scheme(100, 10) for i, group in enumerate(plot_scheme): print(f"代表性环 {i}: 由原始循环 {group} 平均生成")

有了分组方案，下一步是对每组内的多个循环进行平均。难点在于：不同循环在同一相位时间点（比如周期开始后第50秒）的数据点可能不存在，因为采样时间点可能对不齐。我的解决方法是插值。具体步骤是：1. 为每个原始循环，根据其时间列，生成一条应力-相位时间、应变-相位时间的插值函数。2. 定义一个标准化的相位时间轴，比如从0到T（200秒），均匀取100个点。3. 对于每个标准相位时间点，从该组所有循环的插值函数上获取对应的应力和应变值，然后求算术平均。这样就得到了一组新的、在标准相位时间点上对齐的（应力，应变）数据，可以画出一个光滑的、具有代表性的滞回环。

from scipy.interpolate import interp1d def average_cycles(data, cycle_indices, T=200, interp_points=100): """ 对一组原始循环进行平均，生成一个代表性滞回环。 data: 包含'time(s)', 'stress(MPA)', 'strain '的DataFrame cycle_indices: 需要平均的原始循环编号列表 T: 周期 interp_points: 标准化相位时间轴的点数 返回: (avg_strain, avg_stress) 代表平均后的滞回环数据 """ # 标准化的相位时间轴 phase_time_axis = np.linspace(0, T, interp_points) all_stress_interp = [] all_strain_interp = [] for cycle_idx in cycle_indices: # 1. 提取单个循环的数据 mask = (data['time(s)'] >= cycle_idx * T) & (data['time(s)'] < (cycle_idx + 1) * T) cycle_data = data[mask].copy() if len(cycle_data) < 2: continue # 数据太少的循环跳过 # 2. 计算该循环内的相对相位时间 cycle_data['phase_time'] = cycle_data['time(s)'] - cycle_idx * T # 3. 创建插值函数 (使用线性插值，简单稳定) # 注意：需要确保phase_time是单调递增的，通常实验数据都是 f_stress = interp1d(cycle_data['phase_time'], cycle_data['stress(MPA)'], kind='linear', bounds_error=False, fill_value='extrapolate') f_strain = interp1d(cycle_data['phase_time'], cycle_data['strain '], kind='linear', bounds_error=False, fill_value='extrapolate') # 4. 在标准相位时间轴上获取插值 stress_interp = f_stress(phase_time_axis) strain_interp = f_strain(phase_time_axis) all_stress_interp.append(stress_interp) all_strain_interp.append(strain_interp) # 5. 对所有循环的插值结果进行平均 avg_stress = np.mean(all_stress_interp, axis=0) avg_strain = np.mean(all_strain_interp, axis=0) return avg_strain, avg_stress

通过这种方式，我们可以用清晰的几条曲线，展示出滞回环从试验开始到结束的完整演变过程。比如，如果曲线整体向左下方移动，说明在相同的应变幅值下，应力响应降低了，这可能是循环软化的特征；如果向右上方移动，则可能是循环硬化。环的面积变化则反映了材料耗散能量能力的变化。

4. 让图表说话：专业级可视化与解读

数据处理好了，画出专业的图表才能让分析结果更有说服力。对于低周疲劳，通常需要并排展示两个核心图：滞回环演变图和应力-循环次数曲线。

滞回环演变图就是我们上一节用平均化方法得到的。在绘制时，我习惯用渐变色（如viridis或plasma色彩映射）来表征循环的先后顺序，从冷色（开始）到暖色（结束），一目了然。同时，一定要标注坐标轴和单位，应力单位通常是MPa，应变是无量纲量或百分比。

应力-循环次数曲线则反映了材料强度随循环加载的退化情况。通常我们提取每个循环（或每半个循环）的最大应力幅值或最大拉应力，将其与循环次数作图。这条曲线可以直观判断材料的疲劳寿命趋势，是评估材料抗疲劳性能的关键。

def plot_fatigue_analysis(data, T=200, num_representative_cycles=10): """ 综合绘制滞回环演变图和应力-循环次数曲线 """ fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # --- 图1: 平均化滞回环演变 --- total_cycles = int(data['time(s)'].iloc[-1] // T) scheme = generate_averaging_scheme(total_cycles, num_representative_cycles) colors = plt.cm.plasma(np.linspace(0, 0.8, len(scheme))) # 使用plasma渐变色 for i, cycle_group in enumerate(scheme): avg_strain, avg_stress = average_cycles(data, cycle_group, T) ax1.plot(avg_strain, avg_stress, color=colors[i], linewidth=2, label=f'Cycles {cycle_group[0]+1}-{cycle_group[-1]+1}' if i==0 or i==len(scheme)-1 else "") ax1.set_xlabel('Strain', fontsize=12) ax1.set_ylabel('Stress (MPa)', fontsize=12) ax1.set_title('Hysteresis Loops Evolution (Averaged)', fontsize=14, fontweight='bold') ax1.legend(loc='best', fontsize=9) ax1.grid(True, linestyle=':', alpha=0.6) # 添加箭头指示演变方向 ax1.annotate('', xy=(0.02, 300), xytext=(0.015, 320), arrowprops=dict(arrowstyle='->', lw=1.5, color='gray')) ax1.text(0.025, 310, 'Cycle Progression', fontsize=10, color='gray') # --- 图2: 最大应力-循环次数曲线 --- max_stresses = [] cycle_numbers = [] for k in range(total_cycles): mask = (data['time(s)'] >= k * T) & (data['time(s)'] < (k + 1) * T) stress_cycle = data.loc[mask, 'stress(MPA)'] if not stress_cycle.empty: # 取每个循环中应力的最大值（绝对值或拉应力） max_stress = stress_cycle.max() # 也可以用绝对值 max(abs(stress_cycle)) max_stresses.append(max_stress) cycle_numbers.append(k+1) # 循环编号从1开始 ax2.plot(cycle_numbers, max_stresses, 'o-', color='darkred', linewidth=2, markersize=5) ax2.set_xlabel('Cycle Number (N)', fontsize=12) ax2.set_ylabel('Max Stress in Cycle (MPa)', fontsize=12) ax2.set_title('Stress Amplitude vs. Cycle Number', fontsize=14, fontweight='bold') ax2.grid(True, linestyle=':', alpha=0.6) ax2.set_xlim(0, total_cycles+1) # 尝试添加趋势线（例如指数衰减） if len(cycle_numbers) > 3: z = np.polyfit(cycle_numbers, max_stresses, 1) # 线性拟合，可改为更高阶或指数 p = np.poly1d(z) ax2.plot(cycle_numbers, p(cycle_numbers), 'k--', alpha=0.7, label=f'Trend') ax2.legend() plt.tight_layout() plt.show() # 调用函数绘图 plot_fatigue_analysis(data, T=200, num_representative_cycles=8)

这样，我们就得到了一张信息丰富的分析图。左图可以清晰看到滞回环形状、位置和面积的变化，判断材料的硬化/软化行为以及稳定性。右图则定量展示了材料承载能力随循环次数的衰减，我们可以从中估算疲劳寿命（比如应力幅值下降到初始值的某个百分比时所对应的循环次数）。在实际项目中，将这样的图表放入报告，能极大地提升分析的专业性。

5. 进阶分析：从滞回环中提取材料参数

绘制出漂亮的滞回环曲线并不是终点，我们的目标是从中提取有价值的材料力学参数。对于低周疲劳分析，以下几个参数至关重要，而Python可以帮助我们自动计算：

1. 循环应力幅 (Δσ/2)与循环应变幅 (Δε/2)：对于一个稳定的滞回环，可以计算其应力范围和应变范围的一半。这需要识别每个环的上下包络线或转折点。一种简单的方法是找到每个环的最大、最小应力和对应的应变。
2. 塑性应变幅 (Δε_pl/2)：这是评估低周疲劳寿命的核心参数。它可以通过计算滞回环的宽度（在零应力点对应的应变差值）来近似，更准确的方法是分离总应变幅中的弹性部分。弹性应变幅 ≈ 应力幅 / 杨氏模量 (E)，那么塑性应变幅 ≈ 总应变幅 - 应力幅 / E。
3. 滞回能 (ΔW)：即单个循环加载卸载过程中，应力-应变曲线所围成的面积。这个面积代表了材料在每个循环中耗散的能量（主要以热能形式），与疲劳损伤累积直接相关。计算面积可以使用数值积分方法，比如梯形法则 (numpy.trapz)。

def extract_cycle_parameters(cycle_strain, cycle_stress, E=210000): # E: 杨氏模量，单位MPa，示例为钢 """ 从一个滞回环数据中提取关键参数。 假设输入的是单个循环的应力和应变数组。 """ # 找到最大最小应力及对应应变 max_stress_idx = np.argmax(cycle_stress) min_stress_idx = np.argmin(cycle_stress) stress_amplitude = (cycle_stress[max_stress_idx] - cycle_stress[min_stress_idx]) / 2 strain_amplitude = (cycle_strain[max_stress_idx] - cycle_strain[min_stress_idx]) / 2 # 计算滞回环面积（近似耗散能） # 使用梯形数值积分，注意确保数据构成闭合环（首尾点应力应变接近） hysteresis_energy = np.trapz(cycle_stress, cycle_strain) # 单位: MPa (即 MJ/m^3) # 注意：面积可能为负，取决于积分方向，取绝对值 hysteresis_energy = abs(hysteresis_energy) # 估算塑性应变幅 (简化方法) elastic_strain_amp = stress_amplitude / E plastic_strain_amp = strain_amplitude - elastic_strain_amp # 防止计算误差导致负值 plastic_strain_amp = max(plastic_strain_amp, 0) return { 'stress_amp': stress_amplitude, 'strain_amp': strain_amplitude, 'plastic_strain_amp': plastic_strain_amp, 'hysteresis_energy': hysteresis_energy } # 对每个平均化后的代表性循环提取参数 all_params = [] for i, cycle_group in enumerate(plot_scheme): avg_strain, avg_stress = average_cycles(data, cycle_group, T) params = extract_cycle_parameters(avg_strain, avg_stress, E=210000) params['cycle_group'] = f'{cycle_group[0]+1}-{cycle_group[-1]+1}' all_params.append(params) # 将结果转为DataFrame便于查看和绘图 params_df = pd.DataFrame(all_params) print(params_df[['cycle_group', 'stress_amp', 'plastic_strain_amp', 'hysteresis_energy']].round(4))

有了每个代表性循环的关键参数，我们就可以进一步绘制例如“塑性应变幅-循环次数”、“滞回能-循环次数”等曲线，深入分析材料的疲劳行为。更进一步，可以拟合著名的Coffin-Manson方程（Δε_pl/2 = ε_f' * (2N_f)^c）来估算材料的疲劳延性系数 ε_f' 和疲劳延性指数 c。这些拟合工作可以利用scipy.optimize.curve_fit等工具完成，将数据分析推向更深层次。

6. 避坑指南与实战经验分享

在实战中处理低周疲劳数据，我踩过不少坑，这里分享几个最常见的，希望能帮你节省时间。

第一个坑是数据清洗。实验数据很少是完美的。可能会有明显的噪声尖峰、由于设备启动/停止造成的异常点、或者某次循环数据记录不完整。在分段和平均之前，一定要做基本的检查。我常用的方法是快速绘制每个循环的应力-时间曲线，肉眼观察是否有异常。也可以用简单的统计方法，比如计算每个循环数据点的数量，如果某个循环的点数远少于或远多于平均值，就需要警惕。对于明显的噪声尖峰，可以使用滑动中值滤波 (scipy.signal.medfilt) 进行平滑，但要注意滤波窗口不能太大，以免扭曲滞回环的形状。

第二个坑是相位对齐。在平均化时，假设每个循环都从完全相同的应变（或应力）状态开始。但实际数据中，由于控制精度或材料蠕变/松弛，每个循环的起点可能略有漂移。如果直接按固定时间窗口截取然后平均，会导致平均后的环“变胖”或失真。一个改进方法是：以应变（或应力）的过零点或极值点作为循环的起点进行对齐，而不是死板地用时间。这需要先识别每个循环的转折点，然后以这些特征点为基准重新采样数据，再进行平均。虽然复杂一些，但结果更精确。

第三个坑是材料模量E的选取。在计算塑性应变幅时，需要知道材料的杨氏模量E。这个值最好从你本次试验的初始加载段（通常是第一个循环的弹性部分）通过线性回归拟合得到，而不是简单地用一个教科书上的理论值。因为不同批次的材料、不同的试验机刚度都会影响测量到的E值。用试验得到的E值进行计算，结果会更可靠。

第四个是关于绘图细节。发表论文或撰写报告时，图的质量很重要。记得设置合适的图形尺寸 (figsize)、字体大小、线宽和标记大小。给不同的曲线加上清晰的图例。如果循环很多，用颜色映射 (cmap) 比用离散的颜色更专业。可以尝试将关键参数（如最后几个循环的应力幅值、累积塑性应变）以文本形式标注在图上，让读者一目了然。

最后，整个分析过程最好封装成一个类或一系列函数，这样当你拿到新的数据集时，只需要修改文件路径和少数几个参数（如周期T、材料E值），就能一键生成所有分析图表和参数表格，效率会大大提高。数据处理脚本的复用性和可读性，在长期科研或工程项目中非常重要。