Masse Feder System, Theoretische Modellbildung

Wir wollen eine dynamische Modellbildung mit Hilfe der Lagrange II Methode durchführen. Dazu wird die Python Bibliothek sympy verwendet. Als Resultat erhalten wir die die Bewegungsgleichung welche in ein Zustandsraumdarstellung überführt werden kann.

Das Masse Feder Dämpfer (2FG) System

Dieses Modell ist ein sehr wichtiges Model der Mechanik. Je nach Parametrierung dient es als Modell für Antriebstränge im Maschinenbau, für die Modalanalyse von Maschinen oder als Schwingungsmodell für Hochhäuser.

Abbildung 1: MFG

Code: Start

from sympy import init_printing
init_printing()

from sympy import * 

from sympy import Matrix, BlockMatrix, eye, zeros
from sympy import symbols, lambdify
from sympy import simplify

Parameter

Anlegen der physikalischen Parameter.

# parameters
m1, m2, d1, d2, d3, c1, c2, c3 = symbols('m1 m2 d1 d2 d3 c1 c2 c3')

Koordinaten und Kräfte

# forces and moments
F1, F2 = symbols('F1, F2')
display(F1,F2)

# coodinates
q1, q2 = symbols('q_1 q_2')
q1d, q2d = symbols('qdot_1 qdot_2')
q1dd, q2dd = symbols('qddot_1 qddot_2')
display(q1, q1d, q1dd)
display(q2, q2d, q2dd)

\(\displaystyle F_{1}\)

\(\displaystyle F_{2}\)

\(\displaystyle q_{1}\)

\(\displaystyle \dot{q}_{1}\)

\(\displaystyle \ddot{q}_{1}\)

\(\displaystyle q_{2}\)

\(\displaystyle \dot{q}_{2}\)

\(\displaystyle \ddot{q}_{2}\)

# states
x1, x2, x3, x4 = symbols('x1 x2 x3 x4')

# coordinates subs lists
list_q2x = [(q1, x1), (q2, x2), (q1d, x3), (q2d, x4)]
list_x2q = [(x1, q1), (x2, q2), (x3, q1d), (x4, q2d)]

# generalized coordinates (lagrangian coordinates)
q_vec =  Matrix([[q1],[q2]])
q_vec

\(\displaystyle \left[\begin{matrix}q_{1}\\q_{2}\end{matrix}\right]\)

# generalized velocities 
qd_vec =  Matrix([[q1d],[q2d]])
qd_vec

\(\displaystyle \left[\begin{matrix}\dot{q}_{1}\\\dot{q}_{2}\end{matrix}\right]\)

# generalized acceleration 
qdd_vec =  Matrix([[q1dd],[q2dd]])
qdd_vec

\(\displaystyle \left[\begin{matrix}\ddot{q}_{1}\\\ddot{q}_{2}\end{matrix}\right]\)

# state vector
x_vec =  Matrix([[x1],[x2],[x3],[x4]])
x_vec

\(\displaystyle \left[\begin{matrix}x_{1}\\x_{2}\\x_{3}\\x_{4}\end{matrix}\right]\)

# state vector in q
q_and_qd = Matrix([[q_vec],[qd_vec]])
q_and_qd

\(\displaystyle \left[\begin{matrix}q_{1}\\q_{2}\\\dot{q}_{1}\\\dot{q}_{2}\end{matrix}\right]\)

# test subs list
q_and_qd.subs(list_q2x)

\(\displaystyle \left[\begin{matrix}x_{1}\\x_{2}\\x_{3}\\x_{4}\end{matrix}\right]\)

# test sub list
x_vec.subs(list_x2q)

\(\displaystyle \left[\begin{matrix}q_{1}\\q_{2}\\\dot{q}_{1}\\\dot{q}_{2}\end{matrix}\right]\)

Bewegungsgleichung - Position und Geschwindigkeits Vektoren

# postion vector
rs1 = Matrix([[q1],[0],[0]])
rs1

\(\displaystyle \left[\begin{matrix}q_{1}\\0\\0\end{matrix}\right]\)

# velocity vector
vs1 = rs1.jacobian(q_vec)@qd_vec
vs1

\(\displaystyle \left[\begin{matrix}\dot{q}_{1}\\0\\0\end{matrix}\right]\)

# postion vector
rs2 = Matrix([[q2],[0],[0]])
rs2

\(\displaystyle \left[\begin{matrix}q_{2}\\0\\0\end{matrix}\right]\)

# velocity vector
vs2 = rs2.jacobian(q_vec)@qd_vec
vs2

\(\displaystyle \left[\begin{matrix}\dot{q}_{2}\\0\\0\end{matrix}\right]\)

Bewegungsgleichung - Energie Terme und Rayleighfunktion

T = nsimplify(1/2*Matrix([vs1.T*m1*vs1+vs2.T*m2*vs2]))
T

\(\displaystyle \left[\begin{matrix}\frac{m_{1} \dot{q}_{1}^{2}}{2} + \frac{m_{2} \dot{q}_{2}^{2}}{2}\end{matrix}\right]\)

V = nsimplify(Matrix([1/2*c1*(q1)**2 + 1/2*c2*(q2-q1)**2 + 1/2*c3*(-q2)**2]))
V

\(\displaystyle \left[\begin{matrix}\frac{c_{1} q_{1}^{2}}{2} + \frac{c_{2} \left(- q_{1} + q_{2}\right)^{2}}{2} + \frac{c_{3} q_{2}^{2}}{2}\end{matrix}\right]\)

R = nsimplify(Matrix([1/2*d1*(q1d)**2 + 1/2*d2*(q2d-q1d)**2 + 1/2*d3*(-q2d)**2]))
R

\(\displaystyle \left[\begin{matrix}\frac{d_{1} \dot{q}_{1}^{2}}{2} + \frac{d_{2} \left(- \dot{q}_{1} + \dot{q}_{2}\right)^{2}}{2} + \frac{d_{3} \dot{q}_{2}^{2}}{2}\end{matrix}\right]\)

Bewegungsgleichung - Ableitungen

T_qd = nsimplify(T.jacobian(qd_vec).T)
T_qd

\(\displaystyle \left[\begin{matrix}m_{1} \dot{q}_{1}\\m_{2} \dot{q}_{2}\end{matrix}\right]\)

T_qd_dt = nsimplify(T_qd.jacobian(qd_vec)@qdd_vec)
T_qd_dt

\(\displaystyle \left[\begin{matrix}m_{1} \ddot{q}_{1}\\m_{2} \ddot{q}_{2}\end{matrix}\right]\)

M_Matrix = T_qd_dt.jacobian(qdd_vec)
M_Matrix

\(\displaystyle \left[\begin{matrix}m_{1} & 0\\0 & m_{2}\end{matrix}\right]\)

T_q = T.jacobian(q_vec).T
T_q

\(\displaystyle \left[\begin{matrix}0\\0\end{matrix}\right]\)

R_qd = nsimplify(R.jacobian(qd_vec).T)
R_qd

\(\displaystyle \left[\begin{matrix}d_{1} \dot{q}_{1} + \frac{d_{2} \left(2 \dot{q}_{1} - 2 \dot{q}_{2}\right)}{2}\\\frac{d_{2} \left(- 2 \dot{q}_{1} + 2 \dot{q}_{2}\right)}{2} + d_{3} \dot{q}_{2}\end{matrix}\right]\)

G_Matrix = R_qd.jacobian(qd_vec)
G_Matrix

\(\displaystyle \left[\begin{matrix}d_{1} + d_{2} & - d_{2}\\- d_{2} & d_{2} + d_{3}\end{matrix}\right]\)

V_q = nsimplify(V.jacobian(q_vec).T)
V_q

\(\displaystyle \left[\begin{matrix}c_{1} q_{1} + \frac{c_{2} \left(2 q_{1} - 2 q_{2}\right)}{2}\\\frac{c_{2} \left(- 2 q_{1} + 2 q_{2}\right)}{2} + c_{3} q_{2}\end{matrix}\right]\)

K_Matrix = V_q.jacobian(q_vec)
K_Matrix

\(\displaystyle \left[\begin{matrix}c_{1} + c_{2} & - c_{2}\\- c_{2} & c_{2} + c_{3}\end{matrix}\right]\)

Qq = Matrix([[F1],[F2]])

EoM = T_qd_dt - T_q + R_qd + V_q - Qq
EoM

\(\displaystyle \left[\begin{matrix}- F_{1} + c_{1} q_{1} + \frac{c_{2} \left(2 q_{1} - 2 q_{2}\right)}{2} + d_{1} \dot{q}_{1} + \frac{d_{2} \left(2 \dot{q}_{1} - 2 \dot{q}_{2}\right)}{2} + m_{1} \ddot{q}_{1}\\- F_{2} + \frac{c_{2} \left(- 2 q_{1} + 2 q_{2}\right)}{2} + c_{3} q_{2} + \frac{d_{2} \left(- 2 \dot{q}_{1} + 2 \dot{q}_{2}\right)}{2} + d_{3} \dot{q}_{2} + m_{2} \ddot{q}_{2}\end{matrix}\right]\)

Mq = EoM.jacobian(qdd_vec)
Mq

\(\displaystyle \left[\begin{matrix}m_{1} & 0\\0 & m_{2}\end{matrix}\right]\)

gqqd = simplify(T_qd_dt-T_q+R_qd+V_q-Mq@qdd_vec)
gqqd

\(\displaystyle \left[\begin{matrix}c_{1} q_{1} + c_{2} \left(q_{1} - q_{2}\right) + d_{1} \dot{q}_{1} + d_{2} \left(\dot{q}_{1} - \dot{q}_{2}\right)\\- c_{2} \left(q_{1} - q_{2}\right) + c_{3} q_{2} - d_{2} \left(\dot{q}_{1} - \dot{q}_{2}\right) + d_{3} \dot{q}_{2}\end{matrix}\right]\)

Bewegungsgleichung - Resultat

Mq

\(\displaystyle \left[\begin{matrix}m_{1} & 0\\0 & m_{2}\end{matrix}\right]\)

gqqd

\(\displaystyle \left[\begin{matrix}c_{1} q_{1} + c_{2} \left(q_{1} - q_{2}\right) + d_{1} \dot{q}_{1} + d_{2} \left(\dot{q}_{1} - \dot{q}_{2}\right)\\- c_{2} \left(q_{1} - q_{2}\right) + c_{3} q_{2} - d_{2} \left(\dot{q}_{1} - \dot{q}_{2}\right) + d_{3} \dot{q}_{2}\end{matrix}\right]\)

Qq

\(\displaystyle \left[\begin{matrix}F_{1}\\F_{2}\end{matrix}\right]\)

Zustandsgleichung

fxu = nsimplify(Matrix([[qd_vec], [Mq.inv()@(-gqqd+Qq)]]).subs(list_q2x))
fxu

\(\displaystyle \left[\begin{matrix}x_{3}\\x_{4}\\\frac{F_{1} - c_{1} x_{1} - c_{2} \left(x_{1} - x_{2}\right) - d_{1} x_{3} - d_{2} \left(x_{3} - x_{4}\right)}{m_{1}}\\\frac{F_{2} + c_{2} \left(x_{1} - x_{2}\right) - c_{3} x_{2} + d_{2} \left(x_{3} - x_{4}\right) - d_{3} x_{4}}{m_{2}}\end{matrix}\right]\)

gx = fxu.jacobian(Matrix([[F1],[F2]]))
gx

\(\displaystyle \left[\begin{matrix}0 & 0\\0 & 0\\\frac{1}{m_{1}} & 0\\0 & \frac{1}{m_{2}}\end{matrix}\right]\)

fx = simplify(fxu - gx@Matrix([[F1],[F2]]))
fx

\(\displaystyle \left[\begin{matrix}x_{3}\\x_{4}\\- \frac{c_{1} x_{1} + c_{2} \left(x_{1} - x_{2}\right) + d_{1} x_{3} + d_{2} \left(x_{3} - x_{4}\right)}{m_{1}}\\\frac{c_{2} \left(x_{1} - x_{2}\right) - c_{3} x_{2} + d_{2} \left(x_{3} - x_{4}\right) - d_{3} x_{4}}{m_{2}}\end{matrix}\right]\)

A = nsimplify(fx.jacobian(x_vec))
A

\(\displaystyle \left[\begin{matrix}0 & 0 & 1 & 0\\0 & 0 & 0 & 1\\- \frac{c_{1} + c_{2}}{m_{1}} & \frac{c_{2}}{m_{1}} & - \frac{d_{1} + d_{2}}{m_{1}} & \frac{d_{2}}{m_{1}}\\\frac{c_{2}}{m_{2}} & \frac{- c_{2} - c_{3}}{m_{2}} & \frac{d_{2}}{m_{2}} & \frac{- d_{2} - d_{3}}{m_{2}}\end{matrix}\right]\)

B = gx
B

\(\displaystyle \left[\begin{matrix}0 & 0\\0 & 0\\\frac{1}{m_{1}} & 0\\0 & \frac{1}{m_{2}}\end{matrix}\right]\)

Parametergleichung für die Identifikation

Das hier gezeigte Vorgehen funktioniert nur dann, wenn die gesuchten Parameter linear in die Bewegungsgleichung eingehen. Im Allgemeinen ist dies aber nicht der Fall.

para_vec = Matrix([[m1, m2, d1, d2, d3, c1, c2, c3]]).T
para_vec

\(\displaystyle \left[\begin{matrix}m_{1}\\m_{2}\\d_{1}\\d_{2}\\d_{3}\\c_{1}\\c_{2}\\c_{3}\end{matrix}\right]\)

EoM.jacobian(para_vec)

\(\displaystyle \left[\begin{matrix}\ddot{q}_{1} & 0 & \dot{q}_{1} & \dot{q}_{1} - \dot{q}_{2} & 0 & q_{1} & q_{1} - q_{2} & 0\\0 & \ddot{q}_{2} & 0 & - \dot{q}_{1} + \dot{q}_{2} & \dot{q}_{2} & 0 & - q_{1} + q_{2} & q_{2}\end{matrix}\right]\)

simplify(EoM - EoM.jacobian(para_vec)@para_vec+Qq)

\(\displaystyle \left[\begin{matrix}0\\0\end{matrix}\right]\)

Fazit

Wir haben nun ein dynamisches Modell in Form einer Bewegungsgleichung und eines Zustandsraummodells. Diese Gleichungen eignen sich für Simulationen und Animationen.

Die Parametergleichung ist hilfreich das Modell mit der Realität abzugleichen. Das Anpassen der physikalischen Parameter mittels Messdaten nennen wir identifizieren.

Vorher wollen wir aber erste Simulationen und Animation durchführen.