TP DAE : pendule sur une courbe

5.1. TP DAE : pendule sur une courbe#

Marc BUFFAT dpt mécanique, Université Lyon 1

vous devez écrire les fonctions dont le nom est fixé, mais dont vous devez spécifier les arguments en fonction de la question posée.

Chaque fonction validée rapporte des points

Attention: executer toutes les cellules depuis le début en utilisant le bouton run

# initialisation
import os,sys
import numpy as np
import matplotlib.pyplot as plt
from validation.validation import check_function,liste_functions,info_etudiant
from IPython.display import Markdown, display, HTML
def printmd(string):
    display(Markdown(string))
# test si numero étudiant spécifier
try: NUMERO_ETUDIANT
except NameError: NUMERO_ETUDIANT = None 
if type(NUMERO_ETUDIANT) is not int :
    printmd("**ERREUR:** numéro d'étudiant non spécifié!!!")
    NOM,PRENOM,NUMERO_ETUDIANT=info_etudiant()
# parametres spécifiques
_uid_    = NUMERO_ETUDIANT 
_precis_ = 1.0e-5
printmd("**Login étudiant {} {} uid={}**".format(NOM,PRENOM,_uid_))
np.random.seed(_uid_)
_m_=np.round(1.0+2*np.random.rand(),2)  #masse suivant la trajectoire
_l_=np.round(1.0+2*np.random.rand(),2)  #longueur de la barre
_M_=np.round(3.0+4*np.random.rand(),2)  #masse au bout du pendule
_a0_=np.round(2.0/(2+np.random.randint(4)),2)
printmd("**Parametres:** m={}   M={}  l={}  a0={}".format(_m_,_M_,_l_,_a0_))

# bibliotheque scikit
from scikits.odes.dae import dae
from scikits.odes.ode import ode

5.1.1. Pendule sur une courbe#

../../../_images/pendule_courbe.png

5.1.1.1. modélisation#

Les coordonnées du système sont \(Q=[x(t),y(t),\theta(t)]\), (coordonnées de m et angle de rotation du pendule). Le lagrangien du système est \(L(Q,\dot{Q},t)\) , et les équations de Lagrange sous contrainte s’écrivent

\[\frac{d}{dt}(\frac{\partial L}{\partial\dot{Q}_{i}})-\frac{\partial L}{\partial Q_{i}}=H_{i}+\lambda G_{i}\]

où le vecteur \(\vec{H}\) représente la contribution du travail \(\vec{H}.\delta\vec{Q}\) de la force de frottement \(\vec{F}\), \(\lambda\) le multiplicateur de Lagrange (i.e. la force de liaison) associé à la contrainte \(f(x,y)=0\), et \(\vec{G}\) la contribution due travail de cette force de liaison \(\lambda\vec{G}.\delta\vec{Q}\). $\( G_i = \frac{\partial f}{\partial Q_i}\)$

La courbe \(f(x,y)=0\) est le chameau \(y=x^{2}+a_0\cos(\omega x)\)

A l’instant initial le pendule se trouve la courbe en \(x_0, y_0=f(x_0)\) et on le lâche sans vitesse initiale avec un angle \(\theta=0\). On cherche à déterminer la position d’équilibre du système en fonction de \(x_{0}\).

5.1.1.2. Définition des repères et des points#

  • un référenciel R0 et un point O de référence:

  • système de deux points:

    • Q(xp,yp) de masse M glissant sur la courbe auquel est accroché l’autre point P

    • P distant de l de Q

# sympy avec Lagrange
from sympy import init_printing, simplify, symbols, sin, cos, Matrix
from sympy.physics.mechanics import dynamicsymbols, Point, ReferenceFrame, \
                                    Particle, LagrangesMethod
init_printing()

# parametres et ddl
t, M , m, l, g, a0, omega, K = symbols('t M m l g a0 omega K')
x, y, theta    = dynamicsymbols('x y theta')
xp, yp, thetap = dynamicsymbols('x y theta',1)
# définition des points
## BEGIN SOLUTION
## END SOLUTION

# definition des points matériels
PaP = Particle('Pa_P',P,m)
PaQ = Particle('Pa_Q',Q,M)
PaP.potential_energy = m*g*P.pos_from(O).dot(R0.y)
display(PaP.potential_energy, PaP.kinetic_energy(R0))
PaQ.potential_energy = M*g*Q.pos_from(O).dot(R0.y)
display(PaQ.potential_energy, PaQ.kinetic_energy(R0))

5.1.1.3. Calcul du lagrangien#

# Lagrangien
L = PaP.kinetic_energy(R0) + PaQ.kinetic_energy(R0) - PaP.potential_energy - PaQ.potential_energy
display("L=",L)

# contrainte
f = y - (x**2 + a0*cos(omega*x))
display("contrainte=",f)

# force frottement
FL = [(P, -K*xp*R0.x - K*yp*R0.y)]
display("Force de frottement=",FL)

5.1.2. Bilan des équation#

  • Énergie cinétique

\[T=\frac{1}{2}m(\dot{x}^{2}+\dot{y}^{2})+\frac{1}{2}M\left((\dot{x}+l\dot{\theta}\cos\theta)^{2}+(\dot{y}+l\dot{\theta}\sin\theta)^{2}\right)\]

  • Énergie potentielle

\[U=mgy+Mg(y-lcos\theta)\]

  • Lagrangien

\[L=T-U\]

  • Force de frottement généralisée

\[\begin{split}H=\left[\begin{array}{c} -K\dot{x}\\ -K\dot{y}\\ 0 \end{array}\right]\end{split}\]

  • contrainte $\( f(x,y) = y-(x^{2}+a_0\cos\omega x)\)$

  • gradient de la contrainte $\( G = \left[\begin{array}{c} \frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\\ 0 \end{array}\right]\)$

Ecrire les équations du mouvement, et les transformer en un système d’ordre 1 de la forme:

\[ M \dot{Y} = F(Y,t) \]

avec

\[Y=[x,y,\theta,\dot{x},\dot{y},\dot{\theta},\lambda]\]

on dérivera la contrainte autant de fois que nécessaire pour pouvoir résoudre.

LM = LagrangesMethod(L, [x,y,theta], hol_coneqs=[f], forcelist=FL, frame=R0)
eqs = LM.form_lagranges_equations()
display(eqs)

5.1.2.1. Mise sous forme matricielle#

\[ A X = B(Y) \mbox{ avec } Y=[x,y,\theta,\dot{x},\dot{y},\dot{\theta}] \mbox{ et } X=[\dot{Y},\lambda]\]
A=LM.mass_matrix_full
display(A)

B=LM.forcing_full
display(B)

5.1.2.2. linéarisation#

pour des cooordonnées \(q\) et des vitesses généralisées \(u\), dont \(q_i\) et \(u_i\) sont indépendantes et \(r\) un forcage externe

\begin{split}M \begin{bmatrix} \delta \dot{q} \ \delta \dot{u} \ \delta \lambda \end{bmatrix} = A \begin{bmatrix} \delta q_i \ \delta u_i \end{bmatrix} + B \begin{bmatrix} \delta r \end{bmatrix}\end{split}

xe=0
ye=xe**2+a0*cos(omega*xe)
print(xe,ye)
P0 = { x:xe, y:ye, theta:0, xp:0, yp:0, thetap:0, xp.diff(t):0,yp.diff(t):0,thetap.diff(t):0}
l0 = LM.solve_multipliers(op_point=P0)
P0.update(l0)
LIN = LM.to_linearizer(q_ind=[x,theta],qd_ind=[xp,thetap], q_dep=[y], qd_dep=[yp])
ML,AL,BL = LIN.linearize(op_point = P0)
display(ML)
display(AL)
display(BL)

lambda1 = symbols('lambda1')
X = Matrix([xp,yp,thetap,xp.diff(t),yp.diff(t),thetap.diff(t),lambda1])
Y = Matrix([x,theta,xp,thetap])
ML*X - AL*Y

AL,BL, inp_vec = LM.linearize([x,theta],[xp,thetap],[y],[yp],op_point=P0, A_and_B=True)
display(AL)
display(BL)

Yp = Matrix([xp,thetap,xp.diff(t),thetap.diff(t)])
Yp - AL*Y

# A.N
A1=np.array(AL.subs([(M,_M_),(m,_m_),(l,_l_),(g,10),(K,2),(omega,3.3),(a0,_a0_)])).astype(np.float64)
display(A1)
np.linalg.eig(A1)

5.1.2.3. changement de variables Y#

from sympy import MatrixSymbol,ImmutableMatrix
Y = MatrixSymbol('Y',7,1)
# met derniere colonne de A (xcateur lagrange) dans B
B1 = B - Y[6]*A[:,6]
A1 = A[:,:]
A1[:,6] = [0,0,0,0,0,0,0]
# penalisation de la contrainte
beta = symbols('beta')
B1[6] = B1[6] - beta*f
# substitution
A1=A1.subs([(xp,Y[3]),(yp,Y[4]),(thetap,Y[5]),(x,Y[0]),(y,Y[1]),(theta,Y[2])])
A1 = ImmutableMatrix(A1)
display(A1)
B1=B1.subs([(xp,Y[3]),(yp,Y[4]),(thetap,Y[5]),(x,Y[0]),(y,Y[1]),(theta,Y[2])])
B1 = ImmutableMatrix(B1)
display(B1)

# generation automatique des fonctions
from sympy.utilities.autowrap import autowrap
!rm -rf wrapp*
matA = autowrap(A1,args=(Y,M,m,l,g,K,a0,omega),backend='cython',tempdir='./autowrap')
smbB = autowrap(B1,args=(Y,M,m,l,g,K,a0,omega,beta),backend='cython',tempdir='./autowrap')

5.1.3. Programmation#

5.1.3.1. parametres#

  • les valeurs de m,l,M et a0 sont fixés au début.

  • on choisira la valeur des parametres de penalisation \(\beta_1\) \(\beta_2\)

  • on prendra \(K=2.0\) , \(\omega=3.3\), \(g=10.0\)

  • le seul paramêtre variable est donc la position initiale \(x_0\)

Définir la valeurs de ces parametres et ecrire une fonction résidu qui calcul le résidu du système:

\[ res = M.\dot{Y} - F(Y,t) \]
# parametres
## BEGIN SOLUTION
# penalisation
## END SOLUTION

def residu(t,Y,dY,res):
    '''calcul res=M dY -F(Y,t)'''
    global nit
    nit += 1
    YY = np.zeros((7,1))
    YY[:,0] = Y[:]
    A = matA(YY,M,m,l,g,K,a0,omega)
    B = smbB(YY,M,m,l,g,K,a0,omega,beta2)
    res[:] = A.dot(dY) - B[:,0] 
    return

5.1.3.2. Vérification#

vérifier que pour les 2 positions d’équilibres le résidu est bien nul. On tracera aussi la courbe f(x) et la position des 2 points d’équilibre stables.

## BEGIN SOLUTION
# verification
# tracer
## END SOLUTION

5.1.4. Résolution#

ecrire une fonction solution qui calcule la solution en fonction d’une position initiale x0

def solution(x0):
    '''calcul solution avec CI x0. Renvoie T et Y'''

5.1.4.1. Vérification#

en choissisant une valeur de x0 proche de la position d’équilibre, vérifier que le mouvement du pendule est celui attendu.

## BEGIN SOLUTION
#trace de la trajectoire
## END SOLUTION

5.1.5. Analyse pour x0 grand#

pour une valeur de x0 assez grande vérifier que le pendule peut passer d’une position déquilibre à l’autre.

## BEGIN SOLUTION
#trace de la trajectoire
## END SOLUTION

5.1.6. Etude en fonction de x0#

en faisant varier x0 de 0 2.0, déterminer la position finale en fonction de x0

## BEGIN SOLUTION
## END SOLUTION

Tracer la position finale en fonction de x0

## BEGIN SOLUTION
## END SOLUTION

5.1.7. Conclusion#

écrire vos commentaires et conclusion

5.1.7.1. BEGIN SOLUTION#

5.1.7.2. END SOLUTION#
Contenu