SymPy单位自动化简抵消问题求助：弹簧属性计算场景

背景 #

用SymPy进行弹簧属性计算，程序包含10-11个用户可控变量、十余条混合单位的方程，支持求解任意缺失的非材料属性变量（共5种可选）。此前尝试用pint定义变量，虽能正常化简单位，但无法转换为SymPy符号，只能转回SymPy单位系统。

复现代码 #

from sympy.physics.units import *
from sympy import pi, sqrt, N, Eq, symbols, solve

lbf=Quantity('lbf', abbrev='lbf')
lbf.set_global_relative_scale_factor((convert_to(pound*acceleration_due_to_gravity,newton))/newton, newton)

F, Y, L, A, m, ns, xi, d = symbols('F Y L A m ns xi d')
ssy, alpha, beta, C = symbols('ssy alpha beta C')

F= 20*lbf
Y= 2*inch
L= 3.25*inch
d= .08*inch

m= .145
A= 201.00
A*=kilo*psi*inch**m
ns= 1.20
xi= .15

eqSsy=Eq(ssy,.45*A/(d**m))
ssy=solve(eqSsy)[0]
eqAlpha=Eq(alpha,eqSsy.rhs/ns)
alpha=solve(eqAlpha)[0]
eqBeta=Eq(beta,(8*(1+xi)*F)/(N(pi)*(d**2)))
beta=solve(eqBeta)[0]
eqC=Eq(C,((2*eqAlpha.rhs-eqBeta.rhs)/(4*eqBeta.rhs))+sqrt((((2*eqAlpha.rhs-eqBeta.rhs)/(4*eqBeta.rhs))**2)-(3*eqAlpha.rhs)/(4*eqBeta.rhs)))
C=solve(eqC)[0]
print(ssy, '\n', alpha, '\n', beta, '\n', C)

核心问题 #

计算结果C本该是无量纲值（正确结果为10.5334875999498），但实际输出为带复杂单位的表达式：

1.62e-28*(3.66645442100299e+28*inch**2*psi - 1.54320987654321e+27*lbf + 3.666454421003e+28*sqrt(-0.252539870841386*inch**2*lbf*psi + (inch**2*psi - 0.0420899784735643*lbf)**2))/lbf

原因是计算过程中单位未自动抵消：ssy和alpha的单位是psi，而beta的单位是lbf/inch²，两者未自动统一，导致后续运算单位无法消去。该问题与自定义的lbf单位无关，使用原始单位时同样存在。

临时解决方案（需避免） #

手动修改eqBeta的定义，强制将结果转换为psi单位：

eqBeta=Eq(beta,convert_to((8*(1+xi)*F)/(N(pi)*(d**2)),psi))

但程序支持英制/公制单位切换，硬编码psi（或强制C为无量纲）不符合需求。

需求 #

1. 如何让beta自动化简为合适的压力单位，无需手动指定目标单位？
2. 有没有更简洁高效的多变量方程求解方法？不想为每种可能的缺失变量手动编写分步求解流程。

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