Doubts about constraint torque in a 1 D.O.F. system using auxiliary speeds

[wlllind]

I encountered an interesting concept in literature and decided to tackle it using Kane's equations of motion. The problem is best described by the following figures which is from 'How A Teetered Rotor With Delta-3 Really Works (2000), D.J. Malcolm'.

In the appendix of that paper the equations of motion are derived using Euler's equations of motion for rigid bodies. Moreover, they also include the reaction torque on the teeter pin. This is the moment along the axis that is perpendicular to the rotor rotation axis and the teeter axis. To achieve this using Kane's equations the procedure requires an auxiliary speed. In the code and sketch below this is called $u_2$. When simulating the equations using some basic input and comparing it with another approach for blade flapping the results match very well.

Note that the method labeled with 'flap equation' relies on the classic helicopter blade flap equation of motion, which simply has the form of $I\ddot{x} + c\dot{x} + k x = M_{ext}$. I could provide that code as well, but I deem it unnecessary for now.

The results give me confidence that the system has been set up correctly. Then, I decided to add the auxiliary speed $u_2$ which is used to find the reaction torque $T_h$. To my surprise this gives unexpected results. For example, when assuming $d_3 = 0\degree$ and $I_{zz}=I_{yy}$ and $I_{xx}=0$ where the latter follows from the assumption that the rotor is a stick. The solution for $T_h$ is printed at the end of the file for this situation.

Now, when finding the solution for $T_h$ it appears that it depends on $I_{yy}$, which to me appears as a very strange result. Moreover, in the aforementioned paper the reaction torque is only a function of $I_{xx}$ (which is referred to as $I_z$ in that paper).

• Are there by any chance known bugs in Sympy that can affect the computed result?
• Or am I simply making some errors either in definition of frames or in the computation of the generalized active and inertia forces? (Note that in the code below F1.xreplace({d3:0}) yields $T_a -c_a u_1 - k_a q_1$ and for F2 $T_h$. This seems locical for the active forces).

import sympy as sm
import sympy.physics.mechanics as me


#%% Setting up the model

#establishing symbols
m, hh, g, Ixx, Iyy, Izz = sm.symbols('m, hh, g, I_{xx}, I_{yy}, I_{zz}')
Ixy, Iyz, Ixz = sm.symbols('I_{xy}, I_{yz}, I_{xz}')

eta, Xi, d3, Omega = sm.symbols('eta, Xi, delta_3, Omega')
ka, ca = sm.symbols('k_a, c_a')
F_aero, T_aero = sm.symbols('F_a, T_a')

t = me.dynamicsymbols._t
q1, q2, u1, u2, T_h = me.dynamicsymbols('q1, q2, u1, u2, T_h')

#Establishing reference frames
N = me.ReferenceFrame('N') # Inertial reference frame
A = me.ReferenceFrame('A') # Tower reference frame
B = me.ReferenceFrame('B') # Rotating reference frame of hinge
C = me.ReferenceFrame('C') # Flapping reference frame
D = me.ReferenceFrame('D') # Flapping reference frame with d3

#Establishing points
O = me.Point('O') #Tower base
D0 = me.Point('D0') #Rotor centre

#Set up velocities of points that are known on forehand
O.set_vel(N, 0) #Tower base can' move

#Set up locations of points on forehand
D0.set_pos(O, hh*A.x) #Fix rotor with respect to tower base

#Set up relation between flap in blade frame of reference to hinge' q1
beta = sm.atan(sm.cos(d3)*sm.tan(q1)) #flap angle, necessary for flap stiffness
beta_d = beta.diff(t).xreplace({q1.diff(t) : u1}) #flap velocity, necessary for flap damping

#Orient reference frames with respect to eachother
A.orient_axis(N, Xi, N.z)
B.orient_axis(A, Omega*t, A.x)
C.orient_axis(B, q1, B.z)
D.orient_axis(C, d3, C.y)

#Give the reference frames the correct velocities
B.set_ang_vel(A, Omega*A.x) #Rotational speed of rotor
C.set_ang_vel(B, u1*B.z + u2*B.y) #Flapping velocity and aux. speed for T_h
# D.set_ang_vel(C, u2*C.y) #Aux. speed added in a different way. Commented out on purpose.

#Set inertias
I_D_D0 = me.inertia(D , Ixx, Iyy, Izz, -0, -0, -0) #Inertia of rotor in rotor frame of reference D. Non-diagonal terms are 0 on purpose.

#Set the replacement dictionary
replc = {q1.diff(t) : u1,
         q1.diff(t,2) : u1.diff(t),
         q2.diff(t) : u2,
         q2.diff(t,2) : u2.diff(t)}

#Calculate resultant forces on all points
R_D0 = -m*g*N.y + F_aero * D.y

#Calculate torques on the frames
T_D = T_aero*D.z - ka*beta*D.z - ca*beta_d*D.z + T_h*B.y #last term is reaction torque on the hinge!

#Calculate velocities of points in reference frame N 
N_v_D0 = D0.vel(N)

#Calculate accelarations of points in reference frame N 
N_a_D0 = D0.acc(N)

#Calculate angular velocity of franes in with respect to frame N
N_w_D = D.ang_vel_in(N)

#Calculate angular accelaration of frames w.r.t. N
N_alph_D = D.ang_acc_in(N)

#Calculate partial velocities in reference frame N
N_v_D0_1 = N_v_D0.diff(u1, N)
N_v_D0_2 = N_v_D0.diff(u2, N)

N_w_D_1 = N_w_D.diff(u1, N)
N_w_D_2 = N_w_D.diff(u2, N)

#Compute the generalized active forces on the particles
F1 = N_v_D0_1.dot(R_D0) + N_w_D_1.dot(T_D)
F2 = N_v_D0_2.dot(R_D0) + N_w_D_2.dot(T_D)

#Compute the generalized force vector
Fr = sm.Matrix([F1.xreplace(replc), F2.xreplace(replc)])

#Compute the generalized inertia forces
TDs = -(N_alph_D.dot(I_D_D0) + me.cross(N_w_D, I_D_D0).dot(N_w_D))

F1s = N_v_D0_1.dot(-m*N_a_D0) + N_w_D_1.dot(TDs)
F2s = N_v_D0_2.dot(-m*N_a_D0) + N_w_D_2.dot(TDs)

#Construct the generalized inertia force vector
Frs = sm.Matrix([F1s.xreplace(replc), F2s.xreplace(replc)])

#construct Kane's equation with auxiliary force
fa = Fr + Frs

#make sure the auxilary speeds and accel. are set to 0 because they are!
fa = fa.xreplace({u2 : 0,
                  u2.diff(t) : 0})

#Create the matrices for the system i.e. a udot, udot=0 matrix and find jacobina
ua = sm.Matrix([u1.diff(t), T_h])
ua_zero = {udr:0 for udr in ua}

#Dynamic diff.
Ma = fa.jacobian(ua)
ga = fa.xreplace(ua_zero)

#Kinnematic diff.
Mk = sm.Matrix([-1])
gk = sm.Matrix([u1])

#Solve the system; is not necessary per say but funcitons as a check that it works.
udr_sol = -Ma.LUsolve(ga)

#Print the free symbols such that 'p' can be expressed
Mk.free_symbols|gk.free_symbols|Ma.free_symbols|ga.free_symbols
p = sm.Matrix([Ixx,Ixy,Ixz,Iyy,Iyz,Izz,Omega,T_aero,ca,d3,ka]) #based on the line above

#Define EOM evaluation
u = sm.Matrix([u1])
q = sm.Matrix([q1])
eval_eom = sm.lambdify((q,u,p),[Ma,ga,Mk,gk])

T_h_sol = (udr_sol[1])

print(T_h_sol.xreplace({d3:0, Izz:Iyy, Ixx:0}))

[Peter230655]

I play around with Kane's equations a lot, but I do it a bit differently: I use KanesMethod(...) to do most of the stuff for me.
Is your program on Jupyter notebook? If so, could you send your complete program to me? my email is: peter.stahlecker@gmail.com.
I would try to rewrite it 'my way', then we could compare.
NB: I am retired, so surely I will not steal any of your ideas.

[wlllind]

I think if you just paste the code line for line that it will work. BR

[Peter230655]

I used your parameters, but then used Kane's method. It gave the same result for the reaction torque T_h
So, while I cannot judge whether your set up is correct for a propeller, given your set up the reaction torque seems correct.
I made the same substitutes {d3:0, Izz: Iyy, Ixx:0} you made, but did not study what they mean.

import sympy as sm
import sympy.physics.mechanics as me

#%% Setting up the model

#establishing symbols
m, hh, g, Ixx, Iyy, Izz = sm.symbols('m, hh, g, I_{xx}, I_{yy}, I_{zz}')
Ixy, Iyz, Ixz = sm.symbols('I_{xy}, I_{yz}, I_{xz}')

eta, Xi, d3, Omega = sm.symbols('eta, Xi, delta_3, Omega')
ka, ca = sm.symbols('k_a, c_a')
F_aero, T_aero = sm.symbols('F_a, T_a')

t = me.dynamicsymbols._t
q1, q2, u1, u2, T_h = me.dynamicsymbols('q1, q2, u1, u2, T_h')

#Establishing reference frames
N = me.ReferenceFrame('N') # Inertial reference frame
A = me.ReferenceFrame('A') # Tower reference frame
B = me.ReferenceFrame('B') # Rotating reference frame of hinge
C = me.ReferenceFrame('C') # Flapping reference frame
D = me.ReferenceFrame('D') # Flapping reference frame with d3

#Establishing points
O = me.Point('O') #Tower base
D0 = me.Point('D0') #Rotor centre

#Set up velocities of points that are known on forehand
O.set_vel(N, 0) #Tower base can't move

#Set up locations of points on forehand
D0.set_pos(O, hh*A.x) #Fix rotor with respect to tower base

#Set up relation between flap in blade frame of reference to hinge' q1
beta = sm.atan(sm.cos(d3)*sm.tan(q1)) #flap angle, necessary for flap stiffness
beta_d = beta.diff(t) # .xreplace({q1.diff(t) : u1}) #flap velocity, necessary for flap damping

#Orient reference frames with respect to eachother
A.orient_axis(N, Xi, N.z)
B.orient_axis(A, Omega*t, A.x)
C.orient_axis(B, q1, B.z)
D.orient_axis(C, d3, C.y)

#Give the reference frames the correct velocities
B.set_ang_vel(A, Omega*A.x) #Rotational speed of rotor
C.set_ang_vel(B, u1*B.z + u2*B.y) #Flapping velocity and aux. speed for T_h
# D.set_ang_vel(C, u2*C.y) #Aux. speed added in a different way. Commented out on purpose.

#Set inertias
I_D_D0 = me.inertia(D , Ixx, Iyy, Izz, -0, -0, -0) #Inertia of rotor in rotor frame of reference D. Non-diagonal terms are 0 on purpose.

#Calculate resultant forces on Do
R_D0 = -m*g*N.y + F_aero * D.y

#Calculate torques on the frame D
T_D = T_aero*D.z - ka*beta*D.z - ca*beta_d*D.z + T_h*B.y #last term is reaction torque on the hinge!

#==============================================================================
# Standard Kane's method machinery
#==============================================================================
# the rotor is a rigid body with center D0, body fixed frame D, and mass m, and inertia I_D_D0 about D0 in D frame
Rotor = me.RigidBody('Rotor', D0, D, m, (I_D_D0, D0))
BODY = [Rotor]

# forces and torques acting on the rotor
FL = [(D0, R_D0), (D, T_D)]

q_ind = [q1]
u_ind = [u1]
u_aux = [u2]
kd = [q1.diff(t) - u1]

KM = me.KanesMethod(N, q_ind=q_ind, u_ind=u_ind, kd_eqs=kd, u_auxiliary=u_aux)
(fr, frstar) = KM.kanes_equations(BODY, FL)
MM = KM.mass_matrix_full
force = KM.forcing_full
reaction_force = KM.auxiliary_eqs
Torque = sm.solve(reaction_force, T_h)
print(f'force dynamic symbols {me.find_dynamicsymbols(force)}')
print(f'force free symbols {force.free_symbols}')
print(f'MM dynamic symbols {me.find_dynamicsymbols(MM)}')
print(f'MM free symbols {MM.free_symbols}', '\n')

print('The reaction torque T_h =', Torque[T_h].subs({d3:0, Izz: Iyy, Ixx:0}).simplify())

[Peter230655]

What do you mean?

[moorepants]

Sorry, I started typing something and realized it didn't make sense, but I guess I accidentally submitted the comment.

[Peter230655]

@wlllind:
This Kane's method machinery is really very powerful. All yo uhave to do is to describe your system correctly, and it will do all the rest:

• mass matrix and force term for the numeric integration
• reaction forces
  I have played around with it a lot over the past few years.

[wlllind]

Allright, so it might be in the description of the physical model. One way to achieve another answer is by
uncommenting

# D.set_ang_vel(C, u2*C.y) #Aux. speed added in a different way. Commented out on purpose.

and changing

C.set_ang_vel(B, u1*B.z + u2*B.y) #Flapping velocity and aux. speed for T_h

to

C.set_ang_vel(B, u1*B.z) #Flapping velocity and aux. speed for T_h

Basically, I was defining the auxiliary speed $u_2$ as the speed of frame C w.r.t. frame B. I still do not know why this would be wrong, nevertheless, I am now defining $u_2$ as speed of frame D w.r.t. C. The answer is different with a $cos(q_1)^2$ which is not surprising as $u_2$ is now in a different reference frame. That being said, I am really not understanding why still the $I_{yy}$ appears.

Especially as the aforementioned paper yields $0$ for the values used.

[wlllind]

Could you give it a try starting from blank? The important idea is that blade $D$ has an aerodynamic torque $T_a$ acting on it and a damper and spring constant $c_a,k_a$ restraining the movement. It is a relatively simple system and even the solution I found for $q_1$ matches the flap-equation approach perfectly.

[Peter230655]

I am no expert in propellers, so I have no way to judge whether your description is right or wrong.
I do not really see why it should not depend on I_yy, do you have a reason, why it should not?

[wlllind]

After taking your suggestion into account that if modeled correctly, the method yields the correct EOMs, I found that I made a mistake in the FBD . See below the correct one. When I now work out Kane's equations then I obtain the expected result..

I am grateful for your time.

[Peter230655]

Pleasure to be of help!
I have played around with Kane's method for a few years now, just modelling (mostly useless) systems.
Some of them are available at Github's pst-notebooks repository, should you want took take a closer look at how Kane's mehod really helps to set up EOMs. Of course, you still have to model your system, but I am convinced if you model it right, it gives the right OEMs and the right reaction forces / torques.