Cinemática de la pala de un rotor:
Vamos a considerar ese orden para los movimientos aunque no tiene porqué ser el caso siempre.
from sympy import init_session
init_session(use_latex=True, quiet=True)
IPython console for SymPy 0.7.5 (Python 3.3.4-64-bit) (ground types: python)
Algunas variables:
M_b = Symbol('M_b') # Masa de la pala
g = Symbol('g') # Gravedad
e = Symbol('e') # Excentricidad de batimiento
from sympy.physics.mechanics import ReferenceFrame, Point, dynamicsymbols
# Definimos nuestra propia clase para que los versores sean IJK
class IJKReferenceFrame(ReferenceFrame):
def __init__(self, name):
super().__init__(name, latexs=['\mathbf{%s}_{%s}' % (idx, name) for idx in ("i", "j", "k")])
Sistema de referencia árbol:
A = IJKReferenceFrame("A")
# Definimos la velocidad de giro del rotor principal $\Omega$:
Omega = Symbol('Omega')
Omega
Vamos introduciendo sistemas de referencia y sus orígenes:
A1 = IJKReferenceFrame("A1")
psi = dynamicsymbols('psi')
A1.orient(A, 'Axis', [psi, A.z])
A1.dcm(A) # T_{A1A}
T_a, T_k, T_h = dynamicsymbols('T_a, T_k, T_h')
T_a
t = symbols('t')
omega1, omega2, omega3 = dynamicsymbols('omega1, omega2, omega3')
theta1, theta2, theta3 = dynamicsymbols('theta1, theta2, theta3')
eqs = [omega1 - theta1.diff(t), omega2 - theta2.diff(t), omega3 - theta3.diff(t)]
eqs
psi.diff() * A.z
A1.ang_vel_in(A)
A2 = IJKReferenceFrame("A2")
beta = dynamicsymbols('beta')
A2.orient(A1, 'Axis', [beta, -A1.y])
A2.dcm(A1) # T_{A2A1}
A2.ang_vel_in(A1)
A2.dcm(A) # T_{A2A}
A2.ang_vel_in(A)
A3 = IJKReferenceFrame("A3")
zeta = dynamicsymbols('zeta')
A3.orient(A2, 'Axis', [zeta, A2.z])
A3.dcm(A2) # T_{A3A2}
A3.dcm(A1) # T_{A3A1}
B = IJKReferenceFrame("B")
theta = dynamicsymbols('theta')
B.orient(A3, 'Axis', [theta, A3.x])
B.dcm(A3) # T_{BA3}
B.dcm(A1) # T_{BA1}
B.dcm(A) # T_{BA}
Con SymPy podemos obtener matrices de giro complicadas de forma casi instantánea. Definamos ahora el punto de acoplamiento de la pala $E$:
O = Point("O")
O.set_vel(A, 0)
E = O.locatenew('E', e * A1.x)
E.pos_from(O)
Y su centro de masas $GB$:
x_GB = Symbol('x_GB')
GB = E.locatenew('GB', x_GB * B.x)
GB.pos_from(E)
La velocidad angular de la pala $\omega_B$ es:
omega_B = B.ang_vel_in(A).subs(psi.diff(), Omega)
omega_B
Expresada en ejes pala por componentes:
omega_B
omega_B.to_matrix(B)
El tensor de inercia en ejes pala, asumiendo que estos son ejes principales de inercia, será:
from sympy.physics.mechanics import inertia
I_theta, I_beta, I_zeta = symbols('I_theta, I_beta, I_zeta')
I_B = inertia(B, I_theta, I_beta, I_zeta)
I_B
I_B.to_matrix(B)
El momento cinético $h_B$ será, por tanto:
h_B = I_B.dot(omega_B).subs(psi.diff(), Omega)
h_B.to_matrix(B)
Para la derivada del momento cinético tendremos simplemente:
dh_B_dt = h_B.dt(A)
dh_B_dt
Aparece una expresión bastante compleja que no merece la pena analizar. Para la aceleración del acoplamiento:
E.set_vel(A1, 0)
v_E = E.v1pt_theory(O, A, A1)
v_E
a_E = v_E.dt(A).subs({psi.diff(): Omega, psi.diff().diff(): 0})
a_E
a_E.to_matrix(B).simplify()
EGB = GB.pos_from(E)
EGB.cross(M_b * a_E).to_matrix(B).simplify()
Y ahora el momento del peso, teniendo en cuenta que la gravedad será $\mathbf{g} = -g \mathbf{k}_{A1}$:
EGB.cross(M_b * (-g * A1.z)).to_matrix(B).simplify()
Vamos a estudiar cada uno de los grados de libertad desacoplados entre sí en vacío, es decir, sin acciones aerodinámicas: $\mathbf{M}_b^{a, E} = 0$. Además, se cumplirá la relación $I_{\zeta} = I_{\theta} + I_{\beta}$.
En este caso tenemos $\beta(t) \neq 0$, $\zeta(t) = \theta(t) = 0$. Las magnitudes analizadas anteriormente se simplifican:
omega_B.express(B).subs({theta: 0, zeta: 0})
_.doit()
h_B.express(B).subs({theta: 0, zeta: 0}).doit()
dh_B_dt
(dh_B_dt.express(B)
.subs({
theta: 0,
zeta: 0,
psi.diff(): Omega,
psi.diff().diff(): 0,
I_zeta: I_theta + I_beta})
.doit()
.to_matrix(B)
.simplify()
.expand(trig=True))
EGB.cross(M_b * a_E).express(B).subs({theta: 0, zeta: 0}).doit()
EGB.cross(M_b * (-g * A1.z)).express(B).subs({theta: 0, zeta: 0}).doit()
.doit()
es para que se simplifiquen expresiones del tipo Derivative(0, t)
.La ecuación de momentos para la pala quedaría, por tanto:
MrxB, MrzB = symbols('M^r_xB, M^r_zB')
k_beta = Symbol('k_beta')
Mr = MrxB * B.x + k_beta * beta * B.y + MrzB * B.z
Mr
lhs = Mr + EGB.cross(M_b * (-g * A1.z)).express(B).subs({theta: 0, zeta: 0}).doit()
rhs = (dh_B_dt.express(B).subs({theta: 0, zeta: 0, psi.diff(): Omega, psi.diff().diff(): 0, I_zeta: I_theta + I_beta}).doit()
+ EGB.cross(M_b * a_E).express(B).subs({theta: 0, zeta: 0}).doit())
lhs.simplify()
rhs.simplify()
eqxB = Eq((lhs & B.x) - (rhs & B.x))
eqxB
eqyB = Eq(expand_trig(collect((lhs & B.y) - (rhs & B.y), I_beta)))
eqyB
eqzB = Eq((lhs & B.z) - (rhs & B.z))
eqzB
Si tomamos la ecuación en $y_B$ y linealizamos teniendo en cuenta $|\beta(t)| \ll 1$:
beta_ = Dummy('beta')
eqyB_lin = Eq((eqyB.lhs
.subs(beta, beta_)
.series(beta_, n=2)
.removeO()
.subs(beta_, beta)
/ I_beta)
.expand()
.collect(beta))
eqyB_lin
En general, el término gravitatorio es despreciable frente al término centrífugo:
assert eqyB_lin.lhs.is_Add
# Filtramos término gravitatorio
args = [add for add in eqyB_lin.lhs.args if not add.has(g)]
# Reconstruimos la ecuación
eqyB_lin = Eq(eqyB_lin.lhs.func(*args))
eqyB_lin
Vamos a manipular convenientemente el lado izquierdo de la ecuación para sacar factor común $\Omega^2$:
eqyB_lin = Eq((eqyB_lin.lhs / Omega ** 2)
.expand()
.collect(beta)
* Omega ** 2)
assert eqyB_lin.lhs.is_Mul
for factor in eqyB_lin.lhs.args:
if factor.is_Add:
adds = factor.args
else:
mul = factor
eqyB_lin = Eq(mul * adds[0] + mul * adds[1])
eqyB_lin
Y ahora sustituyo:
Omega_beta = Symbol('Omega_beta')
eqyB_lin = Eq(eqyB_lin.lhs.subs(
Omega * sqrt(1 + M_b * e * x_GB / I_beta + k_beta / (I_beta * Omega ** 2)),
Omega_beta))
eqyB_lin
Esta es la ecuación diferencial de batimiento en vacío, donde $\Omega_{\beta} = \Omega \sqrt{1 + \frac{M_{b} e}{I_{\beta}} x_{GB} + \frac{k_{\beta}}{I_{\beta} \Omega^{2}}}$.
Nota: No se pueden definir ecuaciones vectoriales.
try:
Eq(Omega * A.z, 3 * A.z)
except Exception as e:
print(e)
SympifyError: Omega*A.z
Issue: ¿Funcionamiento extraño? Yo esperaría que Eq(a, 1) + 1 == Eq(a + 1, 2)
. Issue #2646 https://github.com/sympy/sympy/issues/2646
from sympy.abc import a
Eq(a, 1)
Eq(a, 1) + 1
Si utilizo un objeto RigidBody
para definir la pala y calcular el momento cinético, los resultados no son tan útiles.
from sympy.physics.mechanics import RigidBody
Definimos la pala como un RigidBody
de centro $GB$, sistema de referencia fijo $B$, masa $M_b$ e inercia $I_B$ alrededor de $E$.
Pala = RigidBody('Pala', GB, B, M_b, [I_B, E])
SymPy calcula el tensor central de inercia automáticamente:
Pala.central_inertia
E.set_vel(A, E.v2pt_theory(O, A, A1))
Pala.angular_momentum(E, A) # ¿Qué estoy intentando hallar?
Pero al calcular el momento cinético, SymPy no aprovecha el hecho de que el punto $E$ es fijo en el sistema $B$. Además, he tenido que definir explícitamente la velocidad del punto $E$ en el sistema de referencia $A$.
Otros:
Origin
de un ReferenceFrame
.Point
están un poco desligados de los ReferenceFrame
en sí, y esto choca con la intuición al principio. ¿No sería más apropiado definir la posición de un punto en un determinado sistema de referencia, en lugar de la velocidad?angular_momentum
) hay expresiones más inteligentes si el sólido tiene un punto fijo. ¿Esto se puede detectar, o al menos dar una forma de pedirlo explícitamente?