圆孔周围应力分布（极坐标图）代码问题排查求助

尝试使用受拉无限大板圆孔的Kirsch方程绘制孔周围径向和切向应力随角度θ的分布，假设径向距离r等于孔半径a。运行代码后生成的极坐标图与预期不符，求排查代码疏漏。

给出的Kirsch方程：

radial stress(θ)= applied uniaxial tensile stress[1+ (a²/r²)(1.5cos(2θ)+cos(4θ)+(3a⁴/(2r⁴))cos(4θ)]
tangential(θ)=applied uniaxial tensile stress[1+ (a²/r²)(0.5cos(2θ)-cos(4θ)+(3a⁴/(2r⁴))cos(4θ)]

运行的代码：

import numpy as np
import matplotlib.pyplot as plt

applied_stress = 178.53  # Applied uniaxial tensile stress
hole_radius = 0.042   # Radius of the circular hole


theta = np.linspace(0, 2 * np.pi, 1000)


r = hole_radius  # Assume r = a
radial_stress = applied_stress * (1 + (hole_radius**2 / r**2) * (1.5 * np.cos(2 * theta) + np.cos(4 * theta) + (3 * hole_radius**4 / (2 * r**4)) * np.cos(4 * theta)))
tangential_stress = applied_stress * (1 + (hole_radius**2 / r**2) * (0.5 * np.cos(2 * theta) - np.cos(4 * theta) + (3 * hole_radius**4 / (2 * r**4)) * np.cos(4 * theta)))


fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='polar')

ax.plot(theta, radial_stress, label='Radial Stress')
ax.plot(theta, tangential_stress, label='Tangential Stress')


ax.set_rlabel_position(135)
plt.title('Radial and Tangential Stress Distribution')
plt.legend()


plt.show()

核心问题是使用了错误的Kirsch方程，标准单轴拉伸无限大板圆孔的应力解与你给出的公式不符，导致计算结果偏离预期。

正确的Kirsch方程（单轴拉伸） #

当无限大板受远程单轴拉伸应力σ时，距孔中心r、与拉伸轴夹角θ处的径向和切向应力公式为：

σ_r = (σ/2) * [1 - (a²/r²) + (1 - 4a²/r² + 3a⁴/r⁴) * np.cos(2θ)]
σ_θ = (σ/2) * [1 + (a²/r²) - (1 + 3a⁴/r⁴) * np.cos(2θ)]

其中：

• θ=0°对应拉伸载荷方向
• 当r=a（孔边界）时，径向应力σ_r=0（符合自由表面无径向应力的边界条件），切向应力化简为σ_θ = σ*(1 - 2*np.cos(2θ))，这是孔边应力集中的经典结果（θ=90°时应力达到3σ，θ=0°时为-σ）。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

applied_stress = 178.53  # Applied uniaxial tensile stress
hole_radius = 0.042      # Radius of the circular hole

theta = np.linspace(0, 2 * np.pi, 1000)
r = hole_radius  # r = a，孔边界位置

# 使用正确的Kirsch方程计算应力
radial_stress = (applied_stress / 2) * (
    1 - (hole_radius**2 / r**2) + 
    (1 - 4*(hole_radius**2 / r**2) + 3*(hole_radius**4 / r**4)) * np.cos(2 * theta)
)
tangential_stress = (applied_stress / 2) * (
    1 + (hole_radius**2 / r**2) - 
    (1 + 3*(hole_radius**4 / r**4)) * np.cos(2 * theta)
)

# 绘制极坐标图
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='polar')

ax.plot(theta, radial_stress, label='Radial Stress')
ax.plot(theta, tangential_stress, label='Tangential Stress')

ax.set_rlabel_position(135)
plt.title('Radial and Tangential Stress Distribution (Hole Boundary)')
plt.legend()
plt.show()

简化说明（r=a时） #

当r=a时，hole_radius**2 / r**2 = 1，hole_radius**4 / r**4 = 1，公式可进一步简化，代码更直观：

# r=a时的简化计算
radial_stress = (applied_stress / 2) * (1 - 1 + (1 - 4 + 3) * np.cos(2 * theta))  # 结果恒为0
tangential_stress = (applied_stress / 2) * (1 + 1 - (1 + 3) * np.cos(2 * theta))  # 等价于σ*(1-2cos2θ)

此时径向应力曲线会与θ轴重合（值为0），切向应力曲线完全符合孔边应力集中的预期分布。

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