Massa e molla

g.schgor1
Su richiesta di marco83, pongo un problemino che
comporta un’equazione differenziale del secondo ordine.



Una massa M di 1 kg e’ appesa ad una molla, ed e’ mantenuta
nella sua posizione iniziale (y=0, con molla non in tensione)
da una forza antagonista del peso.
La molla ha un coefficiente K1 = 5 ( occorrono 5 Newton per
allungare la molla di 1 cm)

Il movimento di M e’ poi contrastato da un ammortizzatore
con coefficiente K2, tarabile fra 0 e 10 (questo contrasto e’
proporzionale alla velocita’ di spostamento).

Se al tempo t=0 si toglie la forza , la massa tende ad abbassarsi
fino alla nuova posizione di equilibrio, con movimento
oscillatorio, oppure aperiodico, in funzione del valore di K2.

Si chiede di determinare appunto il valore di taratura di K2 in
modo che che corrisponda al piu’ rapido spostamento aperiodico
senza sovraelongazioni (smorzamento ‘critico’).

Ci si aspetta un programma con calcolatore che disegni il
movimento del baricentro di M in funzione di t (vedi grafico),
ovviamente con possibilita’ di prefissare di volta in volta K2.

Questo e’ un esempio di simulazione al calcolatore della ‘messa
a punto’ di un sistema meccanico, senza dover ricorrere ad una
ricerca sperimentale (ovviamente piu’ lunga e costosa).

Risposte
Marco831
Il programma l'avevo giusto preparato nel week end e funziona quasi per ogni equazione del secondo ordine, lineare e non, con alcuni piccoli accorgimenti. Descrive quindi sia il moto in grande sia il moto in piccolo intorno alla posizione di eqilibrio.

Per non ammazzare il discorso, aspetto uno o due giorni a postarlo, altrimenti sembra un discorso a due tra me e g.schgor

GIOVANNI IL CHIMICO
Ciao, io avrei anch'io la solzuzione, avevo infatti preparato già in passato un programma simile, con due molle...ma non so come postare immagini e soluzioni...

Marco831
Con due molle inclinate diversamente? Hai provato a plottare la soluzione su un piano x-y? Visto che belle le figure di lisauox (o qualcosa del genere...) che ottieni con rigidezze il cui rapporto è un numero irrazionale?

GIOVANNI IL CHIMICO
No, con due molle lungo lo stesso, asse, molto simile al nostro caso, per le figure di Lisajeux non mi sono ancora attrezzato, cmq le conosco abbastanza bene...

Marco831
Beh, ottenerle è semplice:
risolvi l'equazione mx''+kx=0 con due valori diversi di K e stesso valore di m(come dicevo, il cui rapporto è un numero irrazionale). Detti x1 e x2 i due vettori soluzione, esegui il comando:
plot(x1,x2)

GIOVANNI IL CHIMICO
Grazie, cmq lo sapevo già, intendevo dire che non ho ancora preparato un m-file all'uopuo..
Visto che sei molto appassionato di fluido dinamica, non è che puoi consigliarimi un testo magari scaricabile come pdf sulle correnti a pelo libero?
Sono spesso rimasto affascinato a guardare le strane cose...vortici, sollevamenti, onde, zone che restono ferme ed altre che accelerano, nei canali e nei fiumi

Marco831
Allora, per quello che riguarda le correnti a pelo libero on line, puoi andare sul sito del diiar del poli al corso di Idraulica 2b:
http://www.diiar.polimi.it/larcan/idr2-NO/
lì c'è un po di roba, anche se in questo non sono molto ferrato; la mia vera passione è la fluidodinamica delle macchine e la fluidodinamica numerica.

Come testo da consigliarti, visto che hai seguito un solo corso di meccanica dei fluidi, andrei per prima cosa su:
Idraulica autori: D.Citrini G.Noseda edizioni CEA
Al poli è un must lo trovi dappertutto, in ogni biblioteca.
Se ti piace qualcosa di un po più matematico puoi guardare
Meccanica dei fluidi autori: Marchi Rubatta edizioni: UTET

Se invece vuoi proprio ammazzarti a suon di calcolo tensoriale:
Theoretical Hydrodynamics autori ed edizioni non me li ricordo.

Visto che ti affascinano i moti a pelo libero, quando hai tempo passa al laboratorio Fantoli in Leonardo. E' nella parte sinistra della sede principale, vicino al laboratorio di Fisica Sperimentale; lì hanno delle strumentazioni che permettono proprio di vedere i sollevamenti e le onde di cui parli.

Comunque, per la fluidodinamica in generale guarda quei libri disponibili sul sito di aerospaziale che ti avevo consigliato...

Marco831
Beh, nessun'altro a parte me e Giovanni ha idea di come risolvere questo problema???

Dimostratemi che oltre a sapere sapete anche fare!!!!!!!!!!!

GIOVANNI IL CHIMICO
ciao, grazie, per adesso sto guardando gli appunti sul sito di aerospaziale, ho scoperto che sul mio testo di fluido ci sono un paio di capitoli sul moto delle correnti a pelo libero...
Cmq io sono a Ingegneria Chimica a Bologna

g.schgor1
Dopo giorni di vana attesa, sono perplesso per l’assoluta
mancanza di risposte.

Non ditemi che il problema e’ difficile
(se no cosa vi insegnano all’universita’?).
Non ditemi che non e’ interessante
(e’ il fondamento dell’analisi strutturale dinamica).
Non ditemi che non sapete usare il calcolatore
(o e’ questo il punto debole?).

Ma sapete almeno impostare il problema?
Tento di dare una spinta.

Perche’ la massa M sia in equilibrio statico (in posizione y = 0),
occorre che la sommatoria delle forze applicate sia nulla.
(primo semplice principio della statica).

La massa e’ sicuramente sottoposta alla forza di gravita’ terrestre
(il peso P = M*g), che quindi deve essere bilanciato da una forza F
antagonista.

La molla non agisce perche’ in questa posizione non e’ ne tesa ne
compressa: solo spostando M in un’altra posizione y diversa da 0,
si ha una forza = K1*y.
Anche l’ammortizzatore non agisce finche’ M sta fermo: infatti
si e’ detto che questo esercita una forza che si oppone alla ‘rapidita
di variazione di posizione’ = K2*y’ (non fatemi dire che la derivata
prima di y e’ la velocita’!)

Quando pero’ al tempo 0 si toglie la forza F, compaiono al suo posto
la forza d’inerzia (M*y’’), quella elastica (K1*y) e quella dell’ammortizzatore
(K2*y’), per un nuovo equilibrio ‘dinamico’.

In definitiva, dividendo per M si ha : y’’ = -g - (K1/M)*y – (K2/M)*y’
Questa e’ dunque l’equazione del movimento che dobbiamo studiare.

A presto (sperando che ora qualcuno ci provi).

Marco831
Io la soluzione l'ho già preparata.. come avevo detto, aspettavo a postarla per non restringere il tutto ad un discorso a due.
Visto che comunque non si è associato nessuno, entro stasera la pubblico (ora non ho il tempo).

Marco831
Prima di tutto riscrivo l’equazione che descrive l’equilibrio del sistema:

mx’’+rx’+kx=mg;

m = massa in kg;
r = coefficiente di smorzamento viscoso in Ns/m;
k = rigidezza della molla in N/m (pertanto il valore di k sarà 500 N/m);

Faccio notare che, g.schgor ha richiesto che r appartenga al campo 0:10 e ha poi posto la richiesta:
determinare lo smorzamento critico.

Dato che ho cambiato la scala di valutazione di k, cambierò anche gli estremi dell’intervallo di variabilità di r, ossia:
r appartenente a 0:100.

(Si noti che non è errato moltiplicare solo per 10, come sarà chiaro in seguito).

Presento ora il programma che è stato utilizzato per risolvere tale problema:

[x,dx,tempo]=modified(m,r,k,f,y0,dy0,T0,T,h)
n = (T-T0);
z(1)=dy0;
x(1)=y0;
for i=1:n
t=T0+h*I;
z(i+1)=z(i)+h/eval(m)*(eval(f)-r*z(i)-eval(k)*y(i);
y(i+1)=y(i)+z(i+1)*h;
end
dx=z;
tempo=[T0:h:T];
plot(tempo,y)
return

A questo punto ho implementato una routine del tipo :
hold on
m=’1’;
k=’500’;
f=’1*9.81’;
y0=0;
dy0=0 ;
T0=0 ;
T=2 ;
h=0.001 ;
for i=1:10:101
r(i)=i-1;
[x,dx,tempo]=modified(m,r,k,f,y0,dy0,T0,T,h);
end
mi è stato visualizzato un grafico con tutti gli andamenti delle soluzioni con i vari parametri r da 0 a 100.
Dal grafico si evince che lo smorzamento critico doveva essere compreso tra r=40 e r=50;
Ho eseguito un ciclo analogo al precedente ma con i=41:1:51
Si nota dal grafico che rc era compreso tra r=44 e r=45;
Lanciando un altro ciclo con i=45:0.1:46 ho ottenuto un’approssimazione del valore di smorzamento critico del tipo rc =(circa) 44.7.

Ricordando la teoria, si sa che lo smorzamento critico di un sistema vibrante può essere calcolato come rc = 2*m*sqrt(k/m) che nel nostro caso restituisce rc = 44.7214.

Pertanto con il metodo utilizzato, in pochi minuti si è ottenuta un’approssimazione corretta a meno del 4 per diecimila.

Risulta ora utile spiegare nel dettaglio i programmi utilizzati.

La soluzione è stata ottenuta con matlab, il cui linguaggio è praticamente analogo al C++, però non è necessario dichiarare le variabili.
Ci sono principalmente 6 cose che possono risultare astruse; andiamo con ordine:
-quando si scrive … quellochevuoi(i).. significa che si sta considerando l’i-esima colonna del vettore quellochevuoi.
-l’utilizzo del comando eval fa si che le espressioni tra parentesi, come eval(f), vengano valutate al valore attuale della variabile indipendente. Nel mio caso ho scelto come variabile indipendente t, pertanto le espressioni di m,k e f vanno espresse in funzione di t. Inoltre, dato che il comando eval funziona solo su variabili carattere, ho dovuto esprimere m,k,e f non come numeri ma come stringhe, infatti ho scritto m=’1’; e non m=1;
-perché per il parametro r non ho usato eval? Perché per poter implementare una routine ho dovuto variarlo ad ogni passo e non è possibile farlo con una stringa di caratteri. (qualche cultore di matlab potrebbe avere da ridire.. lo so anch’io come si fa ma non volevo incasinare eccessivamente la trattazione…).
-perché le mie routine partono da 1 anziché da 0 (in generale, perché le mie routine partono da un valore superiore di uno rispetto al valore necessario? Vedi i=45:0.1:46…) Pechè matlab non accetta vettori con indicazione zero sulla prima colonna. Tutti i vettori vanno da 1 ad n, quindi per poter definire correttamente le cose ho dovuto implementare i for con l’aggiunta di 1, che viene immediatamente sottratto nella definizione di r(i).
-perché il programma si chiama modified? Perché nella versione originale era stato pensato con l’eval anche per r, in modo da poter considerare equazioni del tutto non lineari. Il nome originale era vibrazioni1gdl.
-perché uso il metodo di eulero in avanti? Perché è esplicito e non richiede risoluzione di sistemi non lineari. Attualmente sto preparano dei programmi runge-kutta, predictor-corrector e multistep combinati per problemi di dinamica caotica.

Veniamo all’ultimo punto delicato: perché ho modificato il campo di variabilità di r da 0:10 a 0:100? Perché il valore dello smorzamento critico dipende con esponente ½ dalla rigidezza. Avendo riscalato k moltiplicando per 100, ho riscalato r moltiplicando per 10.

E con questo è tutto. Se avete dei dubbi……CHIEDETE!!!

g.schgor1
Complimenti a marco83 per la soluzione.
Ed ora posso dare la mia con MathCad.

In altre occasioni ho cercato di mostrare come si
‘integra numericamente’ una funzione, scrivendo una
semplice relazione fra la variabile indicizzata (n) e
la sua precedente (n-1).
Ma in questo caso abbiamo a che fare con 3 funzioni:
la derivata seconda (che per esigenze di sintassi chiameremo
y2), la derivata prima (y1) ed infine la funzione cercata (y).
Tutte 3 sono funzioni del tempo e sono legate fra loro da
3 equazioni che formano un sistema.
Non e’ quindi possibile integrare separatamente ciascuna
di queste, ma dobbiamo procedere ad un calcolo ‘simultaneo’
delle 3.
Per far questo in MathCad, occorre passare ad una struttura
matriciale, che permette di procedere al calcolo contemporaneo
in corrispondenza di ogni istante t(n), scandito ogni ms, fino
ad 8 s (quindi 8000 iterazioni!)
Questo richiede di definire gli stati iniziali (al tempo t(0)) di
ciascuna funzione e di iniziare il calcolo al tempo t(1).
La procedura qui riportata dovrebbe ora essere quindi chiara.

La scelta di K2 = 1 e’ casuale, e si vede che comporta per y un
andamento ondulatorio smorzato che tende ad una posizione
stabile (y = cica 2 cm).
Per renderci conto di questo valore, ragioniamo sul fatto che
al termine del transitorio (teoricamente per t->inf., ma in pratica
dopo pochi secondi) sia y2 che y1 rimangono stabilmente a
zero, quindi nella nostra equazione rimangono solo i termini
del peso e della molla (da cui y = -M*g/K1 = 1.962 cm).

A questo punto comincia la parte piu’ interessante della ricerca.
Variando il valore di K2, varia l’andamento di tutte le 3 funzioni:
verso lo 0 si esaltano le oscillazioni (a 0 si ha un’oscillazione
sinusoidale permanente), verso il valore massimo (10) non si
hano piu’ oscillazioni, ma un andamento puramente esponenziale.

Vedremo prossimamente come simulando la variazione di K2,
sia possibile trovarne il valore ‘critico’ richiesto dal problema.

g.schgor1
Il bello della simulazione con calcolatore e’ che si puo’ agire
sui parametri del modello variandoli a piacimento, ottenendo
immediatamente la corrispondente nuova soluzione (proprio
come se agissimo sul sistema fisico reale).

Nel caso in esame, variando K2 si ottiene un nuovo grafico
dell’andamento delle funzioni studiate, in particolare quella
della posizione, istante per istante (y, in rosso).
Poiche’ l’obiettivo e’ quello di trovare il valore di K2 che fa
raggiungere la nuova posizione nel minor tempo (e senza
oltrepassarla), con pochi tentativi potremmo constatare che
il valore di K2 = 4.5 raggiunge lo scopo:

(un piccolo trucco per facilitare la ricerca: con l’istruzione
min(y) si ottiene il valore minimo raggiunto dalla funzione.
Se quindi cambiamo K2, possiamo controllare qual’e’ il
primo valore che fa coincidere questo minimo con la posizione
di equilibrio finale).

Per i piu’ esigenti, potremo vedere la prossima volta un
procedimento piu’ matematico.

g.schgor1
Se qualcuno puo’ sentirsi sminuito con la soluzione
‘per tentativi’ adottata precedentemente, puo’ sempre
ricorrere al calcolatore per arrivare ad una soluzione
‘calcolata’ (con tanto di decimali, anche se questo,
ingegneristicamente parlando, non ha alcuna importanza).

Ho gia’ avuto modo di presentare soluzioni di equazioni
differenziali con l’ausilio di MathCad , adottando il metodo
di Laplace. Quello che segue e’ un ulteriore esempio
della facilita’ con cui in pochi minuti e senza fare alcun
passaggio ‘manuale’, si risolva un problema non proprio
cosi’ banale.

(come si puo’ vedere il risultato coincide con quanto trovato
da marco83, con il suo programma in matlab).

Che dire a conclusione?
Il calcolatore rimane uno strumento (da solo non risolve i
problemi), ma sapendolo usare se ne possono trarre
enormi vantaggi.
Ma attenzione.
Non basta saper programmare.(non basta essere ‘informatici’).
Occorre prima una solida preparazione nei concetti di base
delle varie materie (matematica, fisica, ecc.) e soprattutto
saper impostare correttamente e concretamente
le soluzioni.
Il resto facciamolo fare al calcolatore.
E’ questo che non sembra cosi’ diffuso.

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.