如何在Python中拟合载荷-位移曲线线性段以计算屈服强度与杨氏模量

Hey there! Great job getting your load-displacement curve plotted already. Let's walk through how to calculate Young's modulus and yield strength from this data—these are key mechanical properties, and Python makes it straightforward once you know the steps.

First, a quick recap to align on the basics:

• Young's Modulus (E): Derived from the linear elastic region of your curve, where load and displacement are directly proportional. It’s stress divided by strain.
• Yield Strength: The stress where the material starts plastic deformation. The most widely used method is the 0.2% offset method—we draw a line parallel to the linear fit but shifted by 0.2% strain, then find its intersection with the original curve.

First, install scipy if you haven’t already (run pip install scipy in your terminal) for linear regression. Then we’ll modify your existing code step by step.

1. Load Data and Initial Plot (Tweaked for Clarity) #

Let’s start with your code, plus a grid to help identify the linear region:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress  # Add this for linear regression

# Load data
data = np.genfromtxt('test.txt', delimiter=None, skip_header=2)
displacement_mm = data[:, 0]  # Rename variables for readability
load_kN = data[:, 1]

# Plot full curve with grid
plt.figure(figsize=(8, 6))
plt.plot(displacement_mm, load_kN, 'b-', label='Raw Data')
plt.xlabel('Displacement [mm]')
plt.ylabel('Load [kN]')
plt.ylim(0, 1)
plt.legend(loc='best')
plt.title('Load vs Actuator Displacement')
plt.grid(True)  # Grid helps spot the linear elastic region
plt.show()

Run this, then note the displacement range where the curve is perfectly straight (e.g., 0 to 0.1 mm)—that’s your linear elastic region.

2. Select the Linear Elastic Region #

Filter your data to only include this straight segment:

# Adjust these bounds to match your curve's linear region!
linear_mask = (displacement_mm >= 0) & (displacement_mm <= 0.1)
x_linear = displacement_mm[linear_mask]
y_linear = load_kN[linear_mask]

3. Perform Linear Fit #

Use linregress to get the slope and intercept of the linear line, plus metrics to check fit quality:

# Linear regression on the elastic region
slope, intercept, r_value, p_value, std_err = linregress(x_linear, y_linear)

# Print fit quality (R-squared close to 1 means a good fit)
print(f"Linear Fit R-squared: {r_value**2:.4f}")
print(f"Slope (kN/mm): {slope:.4f}")

4. Calculate Young's Modulus #

To compute E, you need two sample-specific values (update these for your test piece):

• gauge_length_m: Gauge length in meters (e.g., 0.05 m for 50 mm)
• cross_section_area_m2: Cross-sectional area in m² (e.g., 0.0001 m² for a 10x10 mm square)

Here’s the code for the calculation:

# Sample properties (UPDATE THESE FOR YOUR SAMPLE!)
gauge_length_m = 0.05
cross_section_area_m2 = 0.0001

# Convert units to SI (Newtons, meters)
load_N = load_kN * 1000
displacement_m = displacement_mm / 1000

# Young's Modulus = (Slope of load-displacement curve * gauge length) / cross-sectional area
# Slope is kN/mm = 1e6 N/m, so convert accordingly
e_modulus_pa = slope * 1e6 * gauge_length_m / cross_section_area_m2
e_modulus_gpa = e_modulus_pa / 1e9  # Convert to GPa for readability

print(f"Young's Modulus: {e_modulus_gpa:.2f} GPa")

5. Calculate Yield Strength (0.2% Offset Method) #

This method finds where a line parallel to the linear fit (shifted by 0.2% strain) intersects your raw data:

# 0.2% offset strain (0.002 in decimal)
offset_strain = 0.002

# Convert offset strain to displacement (mm)
offset_displacement_mm = offset_strain * gauge_length_m * 1000

# Equation of the offset line
offset_load = slope * (displacement_mm - offset_displacement_mm) + intercept

# Find intersection between raw data and offset line
diff = load_kN - offset_load
intersection_idx = np.where(np.diff(np.sign(diff)))[0][0]

# Interpolate to get exact yield point values
yield_displacement_mm = np.interp(0, diff[intersection_idx:intersection_idx+2], displacement_mm[intersection_idx:intersection_idx+2])
yield_load_kN = np.interp(yield_displacement_mm, displacement_mm, load_kN)

# Convert yield load to yield strength (GPa)
yield_strength_pa = (yield_load_kN * 1000) / cross_section_area_m2
yield_strength_gpa = yield_strength_pa / 1e9

print(f"Yield Strength (0.2% offset): {yield_strength_gpa:.2f} GPa")

6. Visualize Fit and Yield Point #

Add the fit lines and yield point to your plot for clarity:

plt.figure(figsize=(8,6))
plt.plot(displacement_mm, load_kN, 'b-', label='Raw Data')
plt.plot(x_linear, slope*x_linear + intercept, 'r--', label=f'Linear Fit (R²={r_value**2:.4f})')
plt.plot(displacement_mm, offset_load, 'g--', label='0.2% Offset Line')
plt.scatter(yield_displacement_mm, yield_load_kN, color='red', s=50, label=f'Yield Point ({yield_load_kN:.2f} kN)')

plt.xlabel('Displacement [mm]')
plt.ylabel('Load [kN]')
plt.ylim(0, 1)
plt.legend(loc='best')
plt.title('Load vs Displacement with Fit and Yield Point')
plt.grid(True)
plt.savefig('load_displacement_with_fit.png', dpi=300)
plt.show()

• Sample Properties: Double-check your gauge length and cross-sectional area—incorrect values will lead to wrong results.
• Linear Region: If manual selection feels tedious, you could automate it by finding the segment with the highest R-squared value, but manual selection is simpler for most cases.
• Units: Always verify unit conversions—mixing mm and meters is a common pitfall!

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