求直线与应力-应变曲线交点以确定规定屈服强度的代码问题

我需要找到直线与应力-应变曲线的交点，以此确定规定屈服强度（offset yield strength）。目前代码能成功绘制曲线与直线，但计算出的交点位置错误。strain1和stress1是从Excel提取的大型数组，m和c是对应弹性阶段直线的参数。

原代码如下：

strain_offset = strain1 + 0.002
stress_offset = m[0]*strain1 + c[0]
yield_strength_index = np.argwhere(np.diff(np.sign(stress1 - stress_offset)))[0][0]

def intersection_point(Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2):
    d = (By2-By1)*(Ax2-Ax1)-(Bx2-Bx1)*(Ay2-Ay1)
    if d:
        uA = ((Bx2-Bx1)*(Ay1-By1)-(By2-By1)*(Ax1-Bx1))/d
        uB = ((Ax2-Ax1)*(Ay1-By1)-(Ay2-Ay1)*(Ax1-Bx1))/d
    else:
        return None

    x_intersection = Ax1 + uA * (Ax2 - Ax1)
    y_intersection = Ay1 + uA * (Ay2 - Ay1)
    
    return x_intersection, y_intersection

    
first = yield_strength_index
second = first + 1
# A points from the stress strain curve
Ax1 = strain1[first]
Ay1 = stress1[first]
Ax2 = strain1[second]
Ay2 = stress1[second]
# B points from the offset line
Bx1 = strain_offset[first]
By1 = stress_offset[first]
Bx2 = strain_offset[second]
By2 = stress_offset[second]

# run our function that finds the intersection point
x_intersection, y_intersection = intersection_point(Ax1,Ay1,Ax2,Ay2,Bx1,By1,Bx2,By2)

print(x_intersection,y_intersection)

fig,ax = plt.subplots()


ax.plot(strain1,stress1)
ax.plot(strain_offset,stress_offset)
ax.plot(x_intersection,y_intersection,'go')

ax.set_ylim([0,50000000])
ax.set_xlim([0,0.035])

plt.show()

问题根源 #

1. 偏移直线与曲线的x轴不匹配：原代码中stress_offset对应strain1的应力值，但绘制偏移直线时用了strain_offset = strain1 + 0.002作为x轴，导致stress1（对应x=strain1）和stress_offset（对应x=strain1+0.002）的差值比较无意义，无法正确定位交叉索引。
2. 交点计算的线段端点不对应：计算交点时，偏移直线的端点用了偏移后的x值，和应力应变曲线的端点x值不一致，导致计算出的交点不在两条曲线的实际交叉位置。

修正方案 #

1. 统一偏移直线的x轴 #

规定屈服强度的偏移直线是弹性阶段直线向右平移0.002应变，正确的直线方程为σ = m*(ε - 0.002) + c，对应每个strain1的x值计算应力，确保和曲线的x轴一致：

stress_offset = m[0]*(strain1 - 0.002) + c[0]

2. 修正交点计算的线段端点 #

偏移直线的端点需使用和曲线相同的x值（strain1），通过直线方程计算对应应力：

# 偏移直线的线段端点
Bx1 = strain1[first]
By1 = m[0]*(Bx1 - 0.002) + c[0]
Bx2 = strain1[second]
By2 = m[0]*(Bx2 - 0.002) + c[0]

3. 增加交点有效性判断 #

在交点计算函数中，判断参数uA和uB是否在[0,1]区间内，确保交点确实在选取的线段之间，避免噪声导致的错误。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 假设strain1, stress1, m, c已从Excel加载完成
# 修正偏移直线计算：对应每个应变strain1的应力值
stress_offset = m[0]*(strain1 - 0.002) + c[0]

# 查找曲线与直线交叉的索引（差值符号变化的位置）
cross_indices = np.argwhere(np.diff(np.sign(stress1 - stress_offset)))
if len(cross_indices) == 0:
    print("未找到交点，请检查直线参数或曲线数据")
else:
    yield_strength_index = cross_indices[0][0]

def intersection_point(Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2):
    d = (By2 - By1)*(Ax2 - Ax1) - (Bx2 - Bx1)*(Ay2 - Ay1)
    if d == 0:
        return None  # 线段平行或重合
    
    uA = ((Bx2 - Bx1)*(Ay1 - By1) - (By2 - By1)*(Ax1 - Bx1)) / d
    uB = ((Ax2 - Ax1)*(Ay1 - By1) - (Ay2 - Ay1)*(Ax1 - Bx1)) / d
    
    # 确保交点在两条线段的范围内
    if 0 <= uA <= 1 and 0 <= uB <= 1:
        x_intersection = Ax1 + uA * (Ax2 - Ax1)
        y_intersection = Ay1 + uA * (Ay2 - Ay1)
        return x_intersection, y_intersection
    else:
        return None

first = yield_strength_index
second = first + 1

# 应力应变曲线的线段端点
Ax1 = strain1[first]
Ay1 = stress1[first]
Ax2 = strain1[second]
Ay2 = stress1[second]

# 偏移直线的对应线段端点（同一x值）
Bx1 = strain1[first]
By1 = m[0]*(Bx1 - 0.002) + c[0]
Bx2 = strain1[second]
By2 = m[0]*(Bx2 - 0.002) + c[0]

# 计算交点
intersection = intersection_point(Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2)
if intersection is not None:
    x_intersection, y_intersection = intersection
    print(f"交点坐标：应变={x_intersection:.6f}，屈服强度={y_intersection:.0f}")
else:
    print("所选线段未找到有效交点，请检查索引或数据")

# 绘图
fig, ax = plt.subplots()
ax.plot(strain1, stress1, label='应力-应变曲线')
ax.plot(strain1, stress_offset, label='0.2%偏移直线')
if intersection is not None:
    ax.plot(x_intersection, y_intersection, 'go', markersize=8, label='屈服点')

ax.set_ylim([0, 50000000])
ax.set_xlim([0, 0.035])
ax.set_xlabel('应变')
ax.set_ylabel('应力')
ax.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Clem Entwistle