求解四分之一圆环变形的二阶非线性微分方程Python代码异常问题

推导得到四分之一圆环变形相关公式：

• 应变公式：$\epsilon=(\theta'−1/R)*t/2$
• 控制微分方程：$\theta(s)''=−F/(4EI)*\cos\theta$
• 边界条件：$\theta(0)=\pi/2$，$\theta(L/4)=0$

其中$L$为平均周长，$R$为平均半径，$t$为厚度，$E$为杨氏模量，$I$为惯性矩，$s$为沿圆弧的空间坐标。

编写Python代码求解后，得到的应变结果全部为负值，但理论上应变应在某一角度处改变符号，以下是原代码：

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# Costanti note
F = 10  # Forza massima applicata in Newton
E = 70000  # Modulo di Young in MPascal
R_e = 28  # Raggio esterno in millimetri
R_i = 27.5  # Raggio interno in millimetri
t = R_e - R_i  # Spessore dell'anello in millimetri
h = 5  # Altezza della sezione in millimetri
I = 1/12 * h * t**3  # Momento d'inerzia della sezione in mm^4
R = (R_e + R_i) / 2  # Raggio medio in millimetri
L = 2 * np.pi * R  # Circonferenza media in millimetri

# Equazioni differenziali
def equations(s, y):
    theta = np.pi / 2 - s / R  # Calcolo dell'angolo theta per ogni punto s
    theta_prime = y[1]
    dtheta_prime_ds = - (F / (4 * E * I)) * np.cos(theta)
    return np.vstack((theta_prime, dtheta_prime_ds))

# Condizioni al contorno
def bc(ya, yb):
    return np.array([ya[0] - np.pi / 2, yb[0]])

# Intervallo di integrazione
s = np.linspace(0, L / 4, 100)

# Condizioni iniziali stimate
y_init = np.zeros((2, s.size))
y_init[0] = np.pi / 2

# Risoluzione del problema ai valori al contorno
solution = solve_bvp(equations, bc, s, y_init)

# Verifica della soluzione
if not solution.success:
    raise RuntimeError("La soluzione dell'equazione differenziale non ha avuto successo.")

# Calcolo della soluzione per questi valori di s
s_values = np.linspace(0, L / 4, 100)
theta_values = solution.sol(s_values)[0]
theta_prime_values = solution.sol(s_values)[1]

# Stampa dei valori di theta e theta'
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)

# Calcolo della deformazione epsilon
epsilon_values = (theta_prime_values - (1 / R)) * t / 2

# Stampa dei risultati
print("s:", s_values)
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)
print("epsilon(s):", epsilon_values)

# Visualizzazione dei risultati
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(s_values, theta_values)
plt.xlabel('s (mm)')
plt.ylabel('theta(s)')
plt.title('Soluzione theta(s)')

plt.subplot(3, 1, 2)
plt.plot(s_values, theta_prime_values)
plt.xlabel('s (mm)')
plt.ylabel('theta\'(s)')
plt.title('Derivata theta\'(s)')

plt.subplot(3, 1, 3)
plt.plot(s_values, epsilon_values)
plt.xlabel('s (mm)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon(s)')

plt.tight_layout()
plt.show()

# Grafico separato di epsilon vs theta (in gradi)
theta_degrees = np.degrees(theta_values)
plt.figure(figsize=(8, 6))
plt.plot(theta_degrees, epsilon_values)
plt.xlabel('theta (degrees)')
plt.ylabel('epsilon(s)')
plt.title('Deformazione epsilon vs theta (degrees)')
plt.grid(True)
plt.show()

问题分析与修正 #

核心错误 #

原代码中，微分方程函数里错误地将$\theta$固定为未变形的几何角度$\theta = \pi/2 - s/R$，而非求解的未知变量$y[0]$。这导致求解器没有计算变形后的$\theta(s)$，而是用初始几何代入方程，得到的$\theta'$完全不符合变形规律，最终应变全为负值。

同时，初始猜测值设置不合理：$y_init[0]$全设为$\pi/2$，没有体现$\theta$从$\pi/2$到0的变化趋势，可能影响求解器收敛。

修正后的代码 #

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# 已知参数
F = 10  # 最大作用力（牛顿）
E = 70000  # 杨氏模量（兆帕）
R_e = 28  # 外半径（毫米）
R_i = 27.5  # 内半径（毫米）
t = R_e - R_i  # 圆环厚度（毫米）
h = 5  # 截面高度（毫米）
I = 1/12 * h * t**3  # 截面惯性矩（mm^4）
R = (R_e + R_i) / 2  # 平均半径（毫米）
L = 2 * np.pi * R  # 平均周长（毫米）

# 微分方程组
def equations(s, y):
    theta = y[0]  # theta是待求解的未知变量，而非固定几何角度
    theta_prime = y[1]
    dtheta_prime_ds = - (F / (4 * E * I)) * np.cos(theta)
    return np.vstack((theta_prime, dtheta_prime_ds))

# 边界条件
def bc(ya, yb):
    return np.array([ya[0] - np.pi / 2, yb[0]])

# 积分区间
s = np.linspace(0, L / 4, 100)

# 优化初始猜测值：theta从pi/2线性过渡到0，theta'初始设为接近理论斜率的值
y_init = np.zeros((2, s.size))
y_init[0] = np.linspace(np.pi/2, 0, s.size)
y_init[1] = np.full(s.size, -1/R)  # 初始斜率猜测

# 求解边值问题
solution = solve_bvp(equations, bc, s, y_init)

# 验证求解成功
if not solution.success:
    raise RuntimeError("微分方程求解失败。")

# 获取求解结果
s_values = np.linspace(0, L / 4, 100)
theta_values = solution.sol(s_values)[0]
theta_prime_values = solution.sol(s_values)[1]

# 计算应变epsilon
epsilon_values = (theta_prime_values - (1 / R)) * t / 2

# 打印结果
print("s:", s_values)
print("theta(s):", theta_values)
print("theta'(s):", theta_prime_values)
print("epsilon(s):", epsilon_values)

# 绘图展示
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(s_values, theta_values)
plt.xlabel('s (mm)')
plt.ylabel('theta(s)')
plt.title('theta(s) 解曲线')

plt.subplot(3, 1, 2)
plt.plot(s_values, theta_prime_values)
plt.xlabel('s (mm)')
plt.ylabel('theta\'(s)')
plt.title('theta\'(s) 导数曲线')
plt.axhline(y=1/R, color='r', linestyle='--', label='1/R')
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(s_values, epsilon_values)
plt.xlabel('s (mm)')
plt.ylabel('epsilon(s)')
plt.title('应变epsilon(s)')
plt.axhline(y=0, color='r', linestyle='--')

plt.tight_layout()
plt.show()

# epsilon与theta（角度）关系图
theta_degrees = np.degrees(theta_values)
plt.figure(figsize=(8, 6))
plt.plot(theta_degrees, epsilon_values)
plt.xlabel('theta (角度)')
plt.ylabel('epsilon(s)')
plt.title('应变epsilon vs theta')
plt.axhline(y=0, color='r', linestyle='--')
plt.grid(True)
plt.show()

修正效果 #

修正后，求解器会正确计算变形后的$\theta(s)$，$\theta'(s)$会在某个位置超过$1/R$，从而让应变$\epsilon$出现正负变化，符合理论预期。

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