一个滚动的圆盘，使用Kane的方法和约束力

我们现在将重新审视滚动圆盘的例子，但这次我们将引入非贡献（约束）力。有关更详细的解释，请参见 [Kane1985]。在这里，我们将开启在进行向量运算时自动简化的功能。它使得小问题的输出更美观，但可能会导致较大的向量运算挂起。:

>>> from sympy import symbols, sin, cos, tan
>>> from sympy.physics.mechanics import *
>>> mechanics_printing(pretty_print=False)
>>> q1, q2, q3, u1, u2, u3  = dynamicsymbols('q1 q2 q3 u1 u2 u3')
>>> q1d, q2d, q3d, u1d, u2d, u3d = dynamicsymbols('q1 q2 q3 u1 u2 u3', 1)
>>> r, m, g = symbols('r m g')

这两行引入了寻找约束力所需的额外量。:

>>> u4, u5, u6, f1, f2, f3 = dynamicsymbols('u4 u5 u6 f1 f2 f3')

大部分主要代码与之前相同。:

>>> N = ReferenceFrame('N')
>>> Y = N.orientnew('Y', 'Axis', [q1, N.z])
>>> L = Y.orientnew('L', 'Axis', [q2, Y.x])
>>> R = L.orientnew('R', 'Axis', [q3, L.y])
>>> w_R_N_qd = R.ang_vel_in(N)
>>> R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z)

无滑滚动的定义要求接触点的速度为零；作为引入约束力的一部分，我们必须在这个点引入速度，根据定义，这些速度总是为零。它们垂直于地面，沿着圆盘滚动的路径，并且沿着地面垂直方向。:

>>> C = Point('C')
>>> C.set_vel(N, u4 * L.x + u5 * cross(Y.z, L.x) + u6 * Y.z)
>>> Dmc = C.locatenew('Dmc', r * L.z)
>>> vel = Dmc.v2pt_theory(C, N, R)
>>> I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
>>> kd = [dot(R.ang_vel_in(N) - w_R_N_qd, uv) for uv in L]

正如我们之前在这个过程中引入了三种速度一样，我们也引入了三种力；它们与速度方向相同，并表示这些方向上的约束力。:

>>> ForceList = [(Dmc, - m * g * Y.z), (C, f1 * L.x + f2 * cross(Y.z, L.x) + f3 * Y.z)]
>>> BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
>>> BodyList = [BodyD]

>>> KM = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3], kd_eqs=kd,
...           u_auxiliary=[u4, u5, u6])
>>> (fr, frstar) = KM.kanes_equations(BodyList, ForceList)
>>> MM = KM.mass_matrix
>>> forcing = KM.forcing
>>> rhs = MM.inv() * forcing
>>> kdd = KM.kindiffdict()
>>> rhs = rhs.subs(kdd)
>>> rhs.simplify()
>>> mprint(rhs)
Matrix([
[(4*g*sin(q2) + 6*r*u2*u3 - r*u3**2*tan(q2))/(5*r)],
[                                       -2*u1*u3/3],
[                          (-2*u2 + u3*tan(q2))*u1]])
>>> from sympy import trigsimp, signsimp, collect, factor_terms
>>> def simplify_auxiliary_eqs(w):
...     return signsimp(trigsimp(collect(collect(factor_terms(w), f2), m*r)))
>>> mprint(KM.auxiliary_eqs.applyfunc(simplify_auxiliary_eqs))
Matrix([
[                                      -m*r*(u1*u3 + u2') + f1],
[-m*r*u1**2*sin(q2) - m*r*u2*u3/cos(q2) + m*r*cos(q2)*u1' + f2],
[                -g*m + m*r*(u1**2*cos(q2) + sin(q2)*u1') + f3]])