MATLAB-SIMULINK - Università di Cagliaripisano/SimulinkFeb_2011.pdf · Realizzazione di un modello...
Transcript of MATLAB-SIMULINK - Università di Cagliaripisano/SimulinkFeb_2011.pdf · Realizzazione di un modello...
2
IntroduzioneLibrerie e blocchi elementari (1) Realizzazione di un modello Esempio: costruzione e visualizzazione di una sinusoide Scelta del solutore Modifica delle impostazioni predefinite Utilizzo di variabili dal workspaceLibrerie e blocchi elementari (2) Esportazione dati verso il Workspace e su file esterno Esempio: filtro passa basso Integratore Esecuzione automatizzata di test Analisi spettrali (FFT) Filtraggio digitale Esempio: sistema termico ad 1 e 2 gradi di libertàUtilizzo di blocchi Trasnfer functionEsempio: Simulazione di un sistema di regolazione di temperaturaCreazione di sottosistemiMaskEsempio Distribuzione di temperatura nel rotore di una turbinaRealizzazione di sistemi MIMO LTIVariabili popup e checkboxEsempio: sistema di frenatura con ABSEmbedded Matlab FunctionToolbox avanzati. SimMechanics e SimDriveline
379111620272831374043445963687381859099109112123133
Indice
Alessandro Pisano - [email protected]
3
Finestra di avvio (v. 7.8.0)
Cartella corrente
Avvio SIMULINK
Editor M-files
Alessandro Pisano - [email protected]
5
Possibilità di suddividere il modello complessivo in sottosistemi paralleli o embricati (un sottosistema può essere importato direttamente in un modello di simulazione differente)
Possibilità di definire finestre di parametrizzazione (Masks)
Esportazione nel workspace Matlab dei risultati della simulazione
Esecuzione automatizzata di test
Toolbox avanzati: SimDrivelline e SimMechanics
Alessandro Pisano - [email protected]
9
Pagina di lavoro
Realizzazione di un modello Simulink
10
Realizzazione di un modello Simulink
1. Importare nella pagina di lavoro i blocchi elementari Simulinknecessari, trascinandoli con il mouse dalla rispettiva libreria (drag-and-drop)
3. Collegare tra loro i blocchi Simulink tracciando le opportune linee di interconnessione in modo da realizzare le funzionalità desiderate
2. Parametrizzare i blocchi Simulink nelle rispettive finestre di parametrizzazione, alle quali si accede dalla pagina di lavoro facendo doppio click con il mouse sopra il blocco stesso.
3 fasi
Alessandro Pisano - [email protected]
11
Esempio introduttivo: costruzione e visualizzazione di un segnale sinusoidale
Sono sufficienti due blocchi elementari: un blocco che generi il segnale desiderato, ed un blocco che ne permetta la visualizzazione.
Il primo blocco lo troveremo nella libreria “Sources” (blocco Sine Wave),Il secondo blocco (blocco Scope), si trova nella libreria “Sinks”.
I blocchi necessari vannoimportati nella pagina di lavoroUntitled trascinando con ilmouse (drag-and-drop) l’iconadel blocco all’interno dellapagina di lavoro. Il risultato ditale procedura è mostrato inFigura.
Salvare il modello e attribuire unnome al file con estensione .mdl
Alessandro Pisano - [email protected]
12
Si deve ora collegare l’uscita del generatore di funzione “Sine Wave” conl’ingresso del blocco di visualizzazione “Scope”.
Per effettuare un collegamento tra due blocchi vi è una procedura rapida. Sideve selezionare il blocco di origine (cliccandovi sopra), e si devesuccessivamente selezionare il blocco di destinazione con il tasto ctrlpremuto.
Un collegamento correttamente eseguito vieneindicato come in Figura
In alternativa, si può portare la freccia del mouse nel punto di origine del collegamento e quindi “tracciarlo” tenendo premuto il tasto sinistro del mouse, portandosi fino al punto di destinazione.
Alessandro Pisano - [email protected]
13
Devono ora essere impostati i parametri di ampiezza, frequenza e sfasamentoche definiscono la particolare sinusoide che si desidera generare. A tal fine ènecessario fare doppio click sul blocco “Sine Wave”, e come risultato si apreuna finestra di dialogo all’interno della quale vanno impostati i parametri difunzionamento.
Ampiezza
Bias
Frequenza
Sfasamento
Tasto OKAlessandro Pisano - [email protected]
14
Si deve ora impostare la durata (cioè l’intervallo temporale) della simulazione.
Cliccando sul tasto RUN vieneeseguita la simulazione.Dopo che è stata eseguita lasimulazione si può visualizzare ilsegnale generato cliccando sul bloccoScope.
La durata si può impostare direttamente dai menù della pagina di lavoro
Durata(valore di default 10.0)
Tasto RUN
Alessandro Pisano - [email protected]
15
Ora si aumenti la frequenza della sinusoide da 1 rad/s a 2 rad/s
Si ripeta la simulazione. Si ri-aggiorni il grafico della finestra grafica Scope
cliccando sul pulsante nella barra dei menu della finestra Scope
Grafico “spigoloso”
Si deve andare a modificare il “metododi integrazione”, che definisce il passo didiscretizzazione temporale che vieneimpiegato nella simulazione delmodello.
16
Il metodo di integrazione (Solver) si imposta selezionando il menu della pagina dilavoro Simulation->Configuration Parameters
Alessandro Pisano - [email protected]
17
La scelta del solutore con il quale di risolvono numericamente le equazionidifferenziali del modello è ovviamente irrilevante per il semplice esempio inesame che non coinvolge alcun legame differenziale.
Il motivo della spigolosità del grafico sta nel fatto che il metodo proposto di default(ode45 a passo variabile) ha “scelto” dei passi di discretizzazione temporalepiuttosto elevati, e sono stati quindi generati “pochi campioni” del segnale
In modelli di simulazione complessi, la scelta del solutore numerico (Runge-Kutta,Dormand-Prince, Eulero,…) e delle relative caratteristiche (passo fisso/variabile,etc) va fatta con criterio.
Alessandro Pisano - [email protected]
18
Runge-Kutta e Dormand-Prince sono degli ottimi solutori “general-purpose”.
Per simulazioni “multi-domain” con la copresenza di costanti di tempo moltodifferenti tra loro sono consigliati i metodi a passo variabile dedicati ai problemiStiff (es. ode15s/stiff).
Per simulazioni con elementi discontinui (non-smooth dynamics) i metodi a passovariabile talvolta forniscono risposte non veritiere. Il solutore Eulero a passo fisso,con un passo sufficientemente piccolo, è ritenuto affidabile per sistemi non-smooth.
Si scelga il solutore ode1 (Eulero) a passo fisso, e se ne imposti il Fixed-Step size a0.001
Alessandro Pisano - [email protected]
19
Una scelta ottimale per il solutore bilancia, per il problema in esame, la precisione
della soluzione e la mole di calcoli richiesta, che influenza il tempo di simulazione.
Per identificare un solutore adeguato serve esperienza. Il passo di campionamento
deve essere commisurato alla rapidità di variazione dei segnali in gioco. Quando si
sceglie un solutore a passo variabile si può pensare di introdurre un limite massimo
per il passo adattativo.
Quando il modello non contiene stati continui (non vi sono cioe blocchi “dinamici”
come Integratori, blocchi Transfer Fcn, etc.) Simulink usa il solutore “discrete”
anche se viene specificato un solutore differente,.
Considerazioni aggiuntive
Alessandro Pisano - [email protected]
20
Modifica delle impostazione predefinite per i files Simulink all’apertura
Alessandro Pisano - [email protected]
23
Si ripeta la simulazione e si riaggiorni il grafico
La sinusoide viene però mostrata a partiredall’istante t=5.Sono stati “persi” i campioni precedenti.
Il motivo è che, al fine di non saturare rapidamente la memoria del programma,vige in Simulink una impostazione di default in base alla quale nei blocchi di tipoScope vengono visualizzati e mantenuti in memoria solo gli ultimi 5000 campionidel segnale.
Il grafico della sinusoide è ora correttamenterappresentato.In base alla scelta fatta per il passo fisso delsolutore, vengono ora generati, e interpolati dalgrafico, 1000 campioni per ogni secondo dievoluzione del segnale.
Alessandro Pisano - [email protected]
24
Per modificare tale impostazione per uno specificoblocco si deve cliccare sul pulsante Parameters nellafinestra del blocco Scope
La finestra “Scope Parameters” ha due sottomenu: “General” e “Data History”.
Dal sottomenu Data History si deve disselezionare la check-box “Limit data points to last ..”
Alessandro Pisano - [email protected]
25
Ora ripetendo la simulazione e riaggiornando il grafico la sinusoide vienevisualizzata per intero
Per visualizzare un segnale costituito dalla somma di tre sinusoidi importiamonella pagina di lavoro due nuove istanze del blocco elementare Sine Wave, edimportiamo anche un blocco che rappresenti un nodo sommatore (blocco Sum
dalla libreria dei Commonly Used Blocks)
Alessandro Pisano - [email protected]
26
Dopo aver cancellato la linea di collegamentopreesistente tra il primo blocco Sine Wave edil blocco Scope, si realizzi la connessioneriportata in Figura.
Il blocco Sum deve esserepreliminarmente parametrizzatospecificando il numero di segnali iningresso, ed il segno con il qualeconcorrono alla sommatoria, permezzo di una stringa (es. +++ +)
Scegliamo +++L’aspetto del blocco diventa
Alessandro Pisano - [email protected]
27
Ora si possono assegnare i parametri di ampiezza, frequenza, bias e sfasamento, delletre sinusoidi, rieseguire la simulazione e visualizzare il grafico prodotto dal bloccoScope .
I parametri di un modello possono essere espressi utilizzando delle variabili (es.A1,f1,b1, …) alle quali si può assegnare un valore con un file script Matlab daeseguirsi prima della esecuzione della simulazione.
Tutte le variabili definite nel workspace di Matlab sono disponibili e accessibili daparte dei blocchi Simulink.
Il modello può essere in questomodo riparametrizzato con estremafacilità.
A1=1;
f1=1;
b1=5;
phi1=0;
A2=4;
f2=pi;
b2=2;
phi2=pi/2;
A3=1;
f3=4*pi;
b3=2;
phi3=0;
Alessandro Pisano - [email protected]
28
Libreria “Continuous”
Libreria “Discontinuities”
Libreria “Discrete”
29
Libreria “Lookup Tables”
Libreria “Math operations”
31
Il blocco To Workspace riceve in ingresso il segnale(scalare o vettoriale) che salva nel workspace. Ilblocco si interconnette agli altri come in Figura.Per tracciare un collegamento a partire da uncollegamento preesistente si deve portare il mousenel punto di diramazione, premere il tasto destro, epoi allontanarsi e tracciare il collegamento tenendoil tasto destro premuto, fino a giungere al punto didestinazione del collegamento
Vediamo come esportare in Matlab i dati prodotti eseguendo i modelliSimulink.
Serve il blocco To Workspace dalla libreria Sinks
Alessandro Pisano - [email protected]
32
Il blocco To Workspace va parametrizzato specificando ilnome della variabile che verrà creata nel workspace diMatlab (scegliere y) ed il formato di salvataggio (èopportuno modificare il formato di default Structure eselezionare invece Array) .Per ottimizzare l’impiego della memoria del programma sipuò anche impostare un fattore intero di decimazione(es. con decimation = 10 i dati vengono salvati nelworkspace con uno step temporale 10 volte superiore,quindi si avranno meno elementi nel vettore y).Si mantenga il valore unitario di default.
Il vettore dei tempi viene salvato didefault sotto forma di array con ilnome tout.Bisogna però disabilitare unaimpostazione che limita a 1000 ilnumero massimo di elementi pertout. Si deve andare ne Simulation-
>Configuration Parameters, e nel
menu Data Import/Exportdisselezionare la check-box Limit
Data Points to Last
Alessandro Pisano - [email protected]
33
Lanciando la simulazione vengono creati nel workspace di Matlab gli array y e
tout.
Verificarlo digitando il comando whos
Si può visualizzare in Matlab il grafico del segnale con il comando
plot(tout,y),grid
Alessandro Pisano - [email protected]
34
Se si desidera che i dati siano non solo esportatinel workspace di Matlab ma anche salvati inmodo permanente su un file, si può generare unfile dati nel formato binario .mat con il blocco To
File, dalla libreria Sinks
I files con estensione mat sono detti mat-files.
Esportazione su file dei dati prodotti da modelli Simulink.
Il blocco deve essere parametrizzato specificandonella apposita finestra di configurazione:
-Il nome (comprensivo dell’estensione) del mat-fileche verrà creato (es. y_test1.mat).
-Il nome che verrà assegnato alla variabilequando il mat-file sarà successivamente
aperto in Matlab (scegliere y)E’ possibile impostare una decimazione dei dati.Se si sceglie una “Decimation” > 1 è bene generareanche un vettore dei tempi “sincrono” con lavariabile sottocampionata
Alessandro Pisano - [email protected]
35
Lanciare la simulazione, e verificare come nella cartella di lavoro sia ora presente ilfile y_test1.mat
Possono essere esportati segnali vettoriali. Per mezzo del blocco Mux (libreriaCommonly Used Blocks) si possono “aggregare” i tre segnali sinusoidali in un unicosegnale vettoriale con tre componenti.
Alessandro Pisano - [email protected]
36
La struttura interna dei mat-files prevede la memorizzazione dei dati in una strutturarettangolare
...
321
2322212
1312111
1321
Nmmmm
N
N
NN
tytytyty
tytytyty
tytytyty
ttttt
I mat-files possono essere aperti successivamente in Matlab con il comando load.
>> load y_test1
Viene generata nel workspace di Matlab una variabile matriciale avente il nomespecificato nella finestra di configurazione del blocco To File e la struttura rettangolareriportata sopra.
Alessandro Pisano - [email protected]
37
Filtro passabasso
Consideriamo un filtro passa-basso RC
+Vin Vout
tVtVtVRC inoutout Equazione differenziale
Equazione differenziale esplicitata rispettoalla derivata di ordine più elevato tVtV
RCtV outinout
1
Posso realizzare uno schema di simulazione utilizzando un blocco Integrator, un bloccoSum e un blocco Gain, oltre che ovviamente un generatore di segnale per costruire latensione di ingresso ed un blocco Scope per visualizzare la tensione di uscita.
Alessandro Pisano - [email protected]
38
Modello Simulink.Grazie al blocco Mux è possibile visualizzare i segnali Vin e Vout nel medesimo bloccoScope. Inseriamo anche un blocco “To Workspace” (nome variabile vout, tipo array)
R=1e4; % 10k OhmC=1e-5; % 10 pF
Assegniamo un valore ai parametriscrivendo un semplice script
La costante di tempo del filtro vale RC=0.1 s
tVtVRC
tV outinout 1
Alessandro Pisano - [email protected]
39
Finestra di parametrizzazione del Signal Generator.4 tipologie di segnali disponibili (sinusoidale, onda quadra,dente di sega, random)
Il filtro può essere implementato in forma più compattamediante un blocco Transfer Function
Si devono specificare icoefficienti dei polinomi anumeratore edenominatore della FdTutilizzando la notazioneMatlab per larappresentazione deipolinomi
1
1
sRCsV
sVsF
in
out
]1 [ 1 RCsRC
Rappresentazioni equivalenti
tVtVtVRC inoutout
Alessandro Pisano - [email protected]
40
Integratore
Nella finestra di parametrizzazione dell’integratore il parametro piu importante da settare èla condizione iniziale (Initial Condition), che di default viene impostata pari a zero
Possono essere introdotte saturazioni inferiorie/o superiori sulle uscite dell’integratore.
Se l’integratore riceve in ingresso un segnalevettoriale, genera in uscita un vettore di paridimensione che contiene l’integrale delle diversecomponenti del vettore di ingresso
Può essere anche applicato un reset sull’uscitadel’integratore. Il reset riporta l’uscitadell’integratore al valore della condizione iniziale
Alessandro Pisano - [email protected]
41
Integratore con ingresso vettoriale
Condizioni iniziali diverse per le uscite
Alessandro Pisano - [email protected]
42
Integratore con saturazione superiore
43
% VERIFICA CHE IL MODELLO SIA APERTO, E IN CASO CONTRARIO LO APRE
if isempty(find_system('Name','filtropassabasso'))
open_system('filtropassabasso')
end
R=1e4; % 10k Ohm
C=1e-5; % 10 pF
sim('filtropassabasso');
y1=vout;
figure(1)
plot(tout,vout)
R=1e4; % 10k Ohm
C=2e-5; % 20 pF
sim('filtropassabasso');
y2=vout;
figure(2)
plot(tout,vout)
figure(3)
plot(tout,y1,tout,y2)
Esecuzione automatica di test
Il seguente codice lancia in sequenza due simulazioni con valori diversi dei parametri R eC, memorizza la Vout nelle due prove nei vettori y1 ed y2, e traccia dei grafici delle variesoluzioni in tre finestre grafiche distinte
IMP. Utilizzo della funzione sim(‘model’)
Alessandro Pisano - [email protected]
44
Risposte in frequenza e analisi spettrali
Impariamo ora a visualizzare la risposta in frequenza di filtri lineari, e a visualizzare lospettro di frequenza di un segnale campionato.
Con i valori consigliati per R e C, la pulsazione di taglio è pari a 1/RC = 10 rad/sec ft=1.6 Hz
1
1
sRCsV
sVsF
in
out1.0RC
Per il filtro considerato
omega_t=1/(R*C);
disp(['La frequenza di taglio è: ', num2str(omega_t/(2*pi)),‘ rad/sec']);
R=input('Inserire il valore di R [Ohm] (valore consigliato: R=1e4): \n');
C=input('Inserire il valore di C [Farad] (valore consigliato C=1e-5): \n');
Le seguenti istruzioni richiedono all’utente l’inserimento da tastiera dei parametri R e C
45
Posso definire in Matlab un oggetto di tipo “Transfer Function”
e visualizzarne quindi la suarisposta in frequenza (piùprecisamente: i diagrammisemilogaritmici del modulo indB e della fase della Funzionedi Risposta Armonica F(jw) infunzione della pulsazione w)con il comando Bode
bode(F),grid;
Cambiamo il valore dellacostante di tempo RC, eritracciamo i diagrammi
46
Ora riferiamoci allo schema modificato
Il generatore di segnali è sostituito da un blocco Fcn (libreria User-Defined Functions) che ricevein ingresso il segnale prodotto dal blocco Clock (libreria Sources) , cioè il tempo corrente.Il blocco To Workspace scrive la variabile out, di tipo array.
Il blocco Fcn puo implementare una qualunque funzione statica, e si parametrizzacompilando, nella apposita finestra di configurazione del blocco, una casella di testo inlinguaggio Matlab, con la variabile standard “u” che denota la variabile in ingresso al blocco
Il codice A1*sin(omega1*u)+A2*sin(omega2*u) definisce un segnale somma di duesinusoidi con ampiezza e pulsazione parametrizzate dai coefficienti A1, omega1, A2, omega2.
Il blocco Fcn consente di implementare facilmente segnali con una espressione analiticaanche complessa, che potrebbero richiedere un elevato numero di blocchi elementari.
Alessandro Pisano - [email protected]
47
A1=1;
omega1=1.6*(2*pi); %rad/s (pari alla pulsazione di taglio)
A2=0;
omega2=0;
Assegniamo un valore alle costanti.
Vogliamo visualizzare l’ingresso e l’uscita del filtro, ed i relativi spettri di potenza
figure(1)
plot(tout,out(:,1),'k',tout,out(:,2),'k--'),grid,
title('Segnale di ingresso V_{in} e segnale di uscita V_{out}'),
xlabel('Tempo [s]'),
legend('V_{in}','V_{out}')
axis([0 10 -2 2])
ZOOM
48
function spettro(t,x,n)
% calcolo del vettore delle frequenze
f=0:1/t(length(t)):1/t(2);
f=f';
% calcolo della Fast Fourier Transform
Y=fft(x);
% calcolo dello densità spettrale di potenza normalizzato
% che permette di ottenere un'ampiezza unitaria dello spettro
% per una sinusoide di ampiezza unitaria
P=2*abs(Y)/length(Y);
% creazione grafico nella finestra n-esima
figure(n),
plot(f(1:ceil(length(f)/2)),P(1:ceil(length(P)/2)))
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
Si salvi nella cartella di lavoro il seguente codice nel file spettro.m
La funzione spettro riceve come argomenti, nell’ordine: il vettore dei tempi, il vettore del
segnale, ed il numero della finestra nella quale tracciare il diagramma
Alessandro Pisano - [email protected]
49
Ora si può utilizzare la funzione spettro() per produrre i grafici desiderati
>>spettro(tout,out(:,1),1), axis([0 3 0 1.1])
>>spettro(tout,out(:,2),2), axis([0 3 0 1.1])
Vin Vout
Alessandro Pisano - [email protected]
50
Per ottenere uno spettro maggiormente fedele a quello, ideale, a larghezza nulla, si deveaumentare il tempo di simulazione.
Tsim=100
Tsim=10
Vin Vout
Vin Vout
51
function [freq data]=spettro2(t,x)
% calcolo del vettore delle frequenze
f=0:1/t(length(t)):1/t(2);
f=f';
% calcolo della Fast Fourier Transform
Y=fft(x);
% calcolo dello densità spettrale di potenza normalizzato
% che permette di ottenere un'ampiezza unitaria dello spettro
% per una sinusoide di ampiezza unitaria
P=2*abs(Y)/length(Y);
freq=f(1:ceil(length(f)/2));
data=P(1:ceil(length(P)/2));
Si salvi nella cartella di lavoro il seguente codice nel file spettro2.m
Rispetto alla funzione spettro, la funzione spettro2 non produce il grafico, ma restituisce
all’esterno i due vettori che consentono di produrre il grafico successivamente (ad esempio,all’interno di una struttura subpplot)
Alessandro Pisano - [email protected]
52
[F,X]=spettro2(tout,out(:,1));
figure(4),
subplot(2,1,1),
plot(tout,out(:,1)),title('Segnale di ingresso'), xlabel(‘Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(2,1,2)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
[F,X]=spettro2(tout,out(:,2));
figure(5),
subplot(2,1,1),
plot(tout,out(:,2)),title('Segnale di uscita'), xlabel(‘Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(2,1,2)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
Analizziamo il seguente codice
Alessandro Pisano - [email protected]
53
L’output del precedente codice è il seguente
Concludiamo questo esempio puntualizzando come il modello possa essere resocompletamente parametrizzabile e gestibile da script avendo cura di specificare le grandezzecome la durata della simulazione, le condizioni iniziali degli integratori, il passo didiscretizzazione, etc, per mezzo di costanti simboliche.
Definiamo la variabile Tsim
Alessandro Pisano - [email protected]
54
% VERIFICA CHE IL MODELLO SIA APERTO, E IN CASO CONTRARIO LO APRE
if isempty(find_system('Name','filtropassabasso2'))
open_system('filtropassabasso2')
end
A1=1;
omega1=1.6*(2*pi);
A2=0;
omega2=0;
Tsim=20;
sim('filtropassabasso2');
[F,X]=spettro2(tout,out(:,1));
figure(6),
subplot(4,1,1),
plot(tout,out(:,1)),title('Segnale di ingresso'), xlabel('Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(4,1,2)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
[F,X]=spettro2(tout,out(:,2));
subplot(4,1,3),
plot(tout,out(:,2)),title('Segnale di uscita'), xlabel('Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(4,1,4)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
Analizziamo il seguente codice
55
A1=1;
omega1=3*(2*pi);
A2=0;
omega2=0;
Tsim=20;
sim('filtropassabasso2');
[F,X]=spettro2(tout,out(:,1));
figure(7),
subplot(4,1,1),
plot(tout,out(:,1)),title('Segnale di ingresso'), xlabel('Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(4,1,2)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
[F,X]=spettro2(tout,out(:,2));
subplot(4,1,3),
plot(tout,out(:,2)),title('Segnale di uscita'), xlabel('Tempo[s]')
grid,axis([0 10 -1.5 1.5])
subplot(4,1,4)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 5 0 1.1])
Analizziamo il seguente codice (cont.)
Alessandro Pisano - [email protected]
57
% VERIFICA CHE IL MODELLO SIA APERTO, E IN CASO CONTRARIO LO APRE
if isempty(find_system('Name','filtropassabasso2'))
open_system('filtropassabasso2')
end
R=1e4; % 10k Ohm
C=1e-5; % 10 pF
A1=1;
omega1=1.6*(2*pi);
A2=2;
omega2=10*(2*pi);
Tsim=20;
sim('filtropassabasso2');
[F,X]=spettro2(tout,out(:,1));
figure(8),
subplot(4,1,1),
plot(tout,out(:,1)),title('Segnale di ingresso'), xlabel('Tempo[s]')
grid,axis([0 10 -4 4])
subplot(4,1,2)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 15 0 2.1])
[F,X]=spettro2(tout,out(:,2));
subplot(4,1,3),
plot(tout,out(:,2)),title('Segnale di uscita'), xlabel('Tempo[s]')
grid,axis([0 10 -4 4])
subplot(4,1,4)
plot(F,X),grid
xlabel('Frequenza [Hz]')
ylabel('X(j2 \pi f)')
title('Spettro di potenza normalizzato')
axis([0 15 0 2.1])
2 armoniche distinte
59
Filtraggio digitale
Molto importante nei sistemi di acquisizione dati
m
j
jkj
n
j
jkF
jkF ybyay
01
kkF
kF ybyay 011 1n 0mEs.
Un filtro digitale del primo ordine rappresenta una implementazione discreta del filtro passa-basso studiato nell’esempio precedente
/
1sT
ea
/
0 1 sTeb
= costante di tempo del filtro
Ts = passo di campionamento della sequenza di input
60
Schema SIMULINK
Blocco UNIT DELAY(libreria discrete)
A1=1;
omega1=2*(2*pi);
Tc=0.01;
tau=0.1;
a1=exp(-Tc/tau)
b0=1-exp(-Tc/tau)
Per ruotare di 90° un blocco lo si deveselezionare e si devono premeresuccessivamente i tasti CTRL + R
Alessandro Pisano - [email protected]
61
A1=1;
omega1=2*(2*pi);
Tc=0.01;
tau=0.1;
a1=exp(-Tc/tau)
b0=1-exp(-Tc/tau)
62
clear all
close all
clc
Tc=0.01;
tau=0.1;
a1=exp(-Tc/tau);a2=0;
b0=1-exp(-Tc/tau);b1=0;
ordine=2;
t=0:0.01:20;
y=sin(t);
rum=0.1*rand(1,length(t));
y_rum=y+rum;
figure(1)
plot(t,y);
grid,title('segnale senza rumore')
figure(2)
plot(t,rum);
grid,title('rumore')
figure(3)
plot(t,y_rum);
grid,title('segnale rumoroso')
yf=zeros(1,length(t));
for i=(ordine+1):length(t)
yf(i)=a1*yf(i-1)+a2*yf(i-2)+b0*y_rum(i)+b1*y_rum(i-1)+;
end
figure(4)
plot(t,yf);
grid,title('segnale FILTRATO ')
Implementazione del filtraggio digitale in Matlab
Alessandro Pisano - [email protected]
63
Ora simuliamo un sistema termico
Consideriamo un sistema termico rappresentato da un volume V circondato da unaparete e contenente un fluido
Sia Te(t) [K] la temperatura esterna alla parete, Tf(t) [K] la temperatura del fluido interno al volume, e q(t) [J/s] una sorgente di calore interna al volume. Sia Cf [J/K] la capacita termica del fluido, e sia Kie (J/K s) il coefficiente di scambio termico tra interno ed esterno.
tT f tTe tq
tTtTKtqtTC feieff
Alessandro Pisano - [email protected]
64
Per tradurre una equazione differenziale in termini di una combinazione tra blocchi Simulinksi deve esplicitare l’equazione differenziale rispetto alle derivate di ordine piu elevato.
tTKtTKtqC
tTC
KtT
C
Ktq
CtT fieeie
f
i
f
iee
f
ie
f
f 11
Con un blocco integratore (Integrator) e con dei blocchi “Constant” e “Gain” si può realizzareil seguente modello di simulazione che ipotizza dei valori costanti per q e Te
Si noti come all’ingresso dell’integratore, punto cui corrisponde il segnale dTf/dt, venga
“costruita” elemento per elemento la formula tTtTKtqC
tT feie
f
f 1
Alessandro Pisano - [email protected]
65
Per eseguire la simulazione bisogna assegnare un valore ai parametri Cf [J/K] e Kie [J/Ks],
alla condizione iniziale Tf(0), ed al valore costante dei segnali Te e q.
La condizione iniziale dell’integratore si specifica nella relativa finestra di dialogo
Tf_zero=298.16; %[K], pari a 25°C
q=2000; % J/s q=2kW
Te=323.16; %[K], pari a 50°C;
Cf=2e3; % J/K per 1 kg di sostanza
Kie=1e2; % J / K s , per 1 m^2 di superficie di scambio
Script di parametrizzazione
Alessandro Pisano - [email protected]
66
Si desidera visualizzare Tf sia in gradi Kelvin che in gradi centigradi in 2 finestre Scope separate.La conversione da °K a °C può essere realizzata con il blocco Fcn (libreria User Defined Functions
Si desidera anche modificare il profilo di q(t)
100t
q
q(t)
Si può utilizzare il blocco Step, libreria Sources.
Il blocco Fcn consente di realizzare una funzione statica tra un parametro di ingresso (scalare o vettoriale) e un parametro di uscita scalare.Si deve scrivere l’espressione della funzione, denotando con u (parola riservata) la variabile in ingresso al blocco
Alessandro Pisano - [email protected]
67
Si realizzi lo schema seguente
Grafico della temperatura Tf in gradi centigradi
Alessandro Pisano - [email protected]
68
Realizziamo lo stesso modello in maniera più compatta, utilizzando il blocco “Transfer Function” (libreria Continuous).
sTKsQsTKsTsC eiefieff sTKsQsTKsC eiefief
tTKtqKsC
sT eie
ief
f
1
Notazione impropria, ma chiara
ief KsC
1
q(t)
+
tTK eie
+
Schema Simulink compattoSchema a blocchi
tTKtTKtqtTtTKtqtTC iieeieieieff
sTKsQKsC
sT eie
ief
f
1
Trasformiamo secondo Laplace l’equazione differenziale
tT f
69
Parametrizzazione del blocco Transfer Function
Si devono specificare i coefficienti dei polinomi anumeratore e denominatore della FdT utilizzando lanotazione Matlab per la rappresentazione dei polinomi(un vettore che contiene i coefficienti dl polinomio inordine decrescente rispetto alle potenze di s)
ief KsC
sF
1
Il blocco deve rappresentare la Funz. di Trasf.
] [ iefief KCKsC
tutTKtTC fieff
tTKtqtu eie ] [ ief KC
]1[ tu tT f
Si può “aggirare” la trasformazione secondo Laplace ragionando sui coefficienti della equazione differenziali
70
][ 01 bb tu tT
tubtubtTatTatTa 01012
012 aaa
FdT del secondo ordine
72
Far variare nel tempo la temperatura esterna ed osservare le corrispondenti fluttuazionidella temperatura interna
Nell’ipotesi che q(t) possa assumere solo due valori (il valore nullo q=0, ed un valorecostante q=Q* ) realizzare un sistema di controllo ON-OFF per regolare ad un valoredesiderato la temperatura del fluido interno al volume
Confrontare i risultati ottenuti impiegando come controllore un rele’ con e senza isteresi
Obbiettivi aggiuntivi
Alessandro Pisano - [email protected]
73
Schema SIMULINK FILE: termico1dof_feedback.mdl
Tf_zero=298;.16 %[K], pari a 25°C
q=2000; % J/s
Te=323;.16 %[K], pari a 50°C;
Cf=2e3; % J/K
Kie=1e2; % J / K s
75
Ampiezzaisteresi = 4°C
Ampiezzaisteresi = 2°C
76
Complichiamo il modello
Trattiamo in maniera distinta gli accumuli termici nel volume e nella parte di contorno. Definiamo quindi un modello più complesso che mi fornisca anche l’evoluzione temporale della temperatura della parete Tp(t).
Cf [J/K] è la capacita termica del gas interno al volumeCp [J/K] è la capacita termica del materiale che costituisce la pareteKip [J/K s] è il coefficiente di scambio termico tra l’interno del volume e la parete. Kpe [J/K s] è il coefficiente di scambio termico tra la parete e l’esterno.
tT f tTe tq
tTptemperatura della parete
tTtTKtqtTC fpipff
tTtTKtTtTKtTC fpippepepp
77
Esplicitiamo il sistema di equazioni rispetto alle derivate di ordine più elevato
tTtTC
Ktq
CtT fp
f
ip
f
f 1
Importiamo nella pagina di lavoro un ulteriore blocco integratore la cui uscita sarà la temperatura Tp della parete.
Utilizziamo anche un visualizzatore a Display (libreria Sinks)Lo schema può essere realizzato come segue
tTtTC
KtTtT
C
KtT fp
p
ip
pe
p
pe
p
78
Script di configurazione dei parametri
Tf_zero=25+273.15; %[K], pari a 25°C
Tp_zero=10+273.15; %[K], pari a 10°C
Cf=2e3; % J/K
Cp=10e3; % J/K
Kip=1e2; % J/K s
Kpe=1e2; % J/K s
q=2000; % J/s
Te=20+273.15 %[K], pari a 20°C;
Risultati della simulazione
79
Si desidera acquisire il segnale q(t) da un file esterno.
File “dati_problema.dat”
Il segnale q(t) sia disponibile nella forma di un file dati di tipo ASCII
La prima riga riporta i tempi e la seconda riporta il valore del segnale
Bisogna creare un mat-file (lo chiamiamo dati_mat.mat) che contenga tali informazioni.
%clear all
load dati_problema.dat
M=dati_problema;
tempi=M(1,:); %non utilizzato
segnale=M(2,:); %non utilizzato
save dati_mat M %si puo’ fare direttamente >> save dati_mat dati_problema
Tf_zero=298.16; %[K], pari a 25°C
Tp_zero=10+273.16; %[K], pari a 10°C
Cf=2e3; % J/K
Cp=10e3; % J/K
Kip=1e2; % J/K s
Kpe=1e2; % J/K s
q=2000; % J/s
Te=20+273.16; %[K], pari a 20°C;
80
Per importare il segnale si può utilizzare il blocco Simulink “From File”, dalla libreria Sources
Schema Simulink
In accordo con i dati acquisiti, la durata della simulazione deve essere posta pari a 2.
Alessandro Pisano - [email protected]
81
Creazione di sottosistemi
Si vuole rendere più compatta la rappresentazione del modello attraverso la definizione di un macroblocco come in figura
tT f
tTe
tq
tTp
tTp
Simulink consente di definire dei macroblocchi (sottosistemi) che rappresentano una particolare interconnessione tra altri blocchi, alla quale si accede esplorando il contenuto del sottosistema.
Alessandro Pisano - [email protected]
82
Creazione di sottosistemi
Si dispongano i blocchi del modello come in figura, in modo che sia possibile tracciare un rettangolo nel quale entrano i segnali di input del sottosistema e dal quale escono i segnali di output.
83
Con il mouse si deve “tracciare” nella pagina di lavoro un rettangolo come quello nella slide precedente. Poi, dal menu Edit della pagina di lavoro Simulink, si deve selezionare il comando Create Subsystem
Ridisponendo i blocchi dello schema si puo realizzare la seguente configurazione
Cliccando sul blocco Subsystem si accede al suo contenuto.Si notino i blocchi “In1”,
“In2” , “Out1” ed “Out2”
Alessandro Pisano - [email protected]
84
Rinominando i blocchi “In1”, “In2”,
“Out1” ed “Out2” si può fare inmodo che nelle porte di input edi Output del Subsystemcompaia il nome della grandezzaassociata
La parametrizzazione del modello è attualmente effettuata per mezzo del file di scriptprecedentemente illustrato, che assegna una valore alle costanti simboliche usate nel modelloSimulink definendo opportune variabili omologhe nel workspace di Matlab. Tale m-file discript deve sempre accompagnare il file Simulink (che ha estensione .mdl) e deve esserelanciato prima di quest’ultimo
Può essere conveniente disporre di una maschera di parametrizzazione interna al modello Simulink.
In questo modo tutte le informazioni associate al modello sono contenute in un unico file (il file .mdl) e la riparametrizzazione del modello avviene in maniera più semplice, senza dovere ogni volta rilanciare il file script con estensione .m
Alessandro Pisano - [email protected]
85
Creazione di MASK
Bisogna portarsi con il mouse sul Subsystem“Camera di Combustione”, premere il tasto destro del mouse, e selezionare MaskSubsystem dal menu che compare
Alessandro Pisano - [email protected]
86
Dopo avere selezionato Mask Subsystem si deve nuovamente portarsi con il mouse sul Subsystem “Camera di Combustione”, premere il tasto destro del mouse, e selezionare stavolta Edit Mask dal menu che compare.
Finestra di configurazione della mask
Andare nel sottomenù:Parameters
Alessandro Pisano - [email protected]
87
Nel sottomenù Parameters completare le voci “Dialog Parameters” come in Figura, creando una riga per ognuno dei 6 parametri da settare .
Pulsante “Add”
Per ciascun parametro , oltre al nome della relativa variabile (Cf, Cp, …) si può riportare una frase descrittiva che comparirà nella maschera di configurazione, che riportiamo nella slide seguente.
Aggiunge una riga.
Alessandro Pisano - [email protected]
88
Maschera di parametrizzazione del Subsystem
Viene visualizzata facendo doppio click sul blocco Simulink del Sottosistema
E’ possibile inserire manualmente il valore dei parametri nelle opportune caselle di testo. I valori inseriti vengono memorizzati al salvataggio del file, e riproposti alla sua riapertura.
%clear all
load dati_problema.dat
M=dati_problema;
tempi=M(1,:); %non utilizzato
segnale=M(2,:); %non utilizzato
save dati_mat M
q=2000; % J/s
Te=20+273.16; %[K], pari a 20°C;
Script semplificato
Per accedere al contenuto del blocco “mascherato” premere il tasto destro del mouse e poi selezionare “Look Under Mask”
FILE:
termico_2dof_matfile_mask01.mdl
89
Si inserisca una ulteriore mask che consenta l’impostazione manuale di q e Te
FILE:
termico_2dof_matfile_mask02.mdl
%clear all
load dati_problema.dat
M=dati_problema;
tempi=M(1,:); %non utilizzato
segnale=M(2,:); %non utilizzato
save dati_mat M
Script
Alessandro Pisano - [email protected]
90
Esercizio
91.
92
X0=0.1;Y0=0.1;
A=1;B=1;C=5;D=1;
93
Sistemi termici spazialmente distribuiti
Variabile di interesse: la temperatura
Barra metallica filiforme di lunghezza L = 50 cm
x
Alessandro Pisano - [email protected]
0x Lx
txT ,
Ipotesi : barra uniforme, termicamente isolata, in assenza di sorgenti di calore
txTk
x
txTk
t
txTCp ,
,, 2
2
2
Equazione del calore monodimensionale
94
pC
k
è il coefficiente di diffusione [m2/s]
pC è la capacita termica a pressione costante per unita di massa [J/kg K]
è la conduttività termica [J/s m K ; W / mK ]k
è la densità [kg /m3 ]
txTk
x
txTk
t
txTCp ,
,, 2
2
2
2
2 ,,
x
txT
t
txT
Valori per il ferro puro (rif. Schaum Trasmissione del calore, Tabella B-1 pag. 306)
J/kgK 435pc C) 100 (a W/mK 63 kkg/m 7856
C) 100 (a /m 108.1 25 s
PDE
95
Analisi in regime stazionario
0
2
2
x
xT
x
0x Lx
h h h h h 5/Lh
4 “NODI SOLUZIONE”
1xx 2xx 3xx 4xx
La soluzione è calcolabile in forma chiusa, ma scegliamo di risolvere l’equazionenumericamente perché questo modo di ragionare di consentirà di risolvere agevolmenteanche i problemi distribuiti in regime transitorio , che vedremo in seguito.
xTtxT ,
Siano noti e mantenuti costanti i valori di temperatura agli estremi 00 TT LTTL
xL
TTTxT L
0
0Soluzione in forma chiusa
96
2
012
2
2 2
1
h
TTTxT
xxx
Definiamo il vettore che contiene le temperature incognite nei nodi soluzione.
TTTTT 4,321 .,,T
2
123
2
2 2
2
h
TTTxT
xxx
2
234
2
2 2
3
h
TTTxT
xxx
2
34
2
2 2
4
h
TTTxT
x
L
xx
Approssimiamo l’equazione in un intorno dei nodi soluzione, e mettiamo a sistema
02 012 TTT
02 123 TTT
02 234 TTT
02 34 TTTL
97
Isoliamo alla sinistra dell’uguale le temperature incognite, ordinandole
0212 TTT
02 321 TTT
02 432 TTT
LTTT 43 2
bAT
2100
1210
0121
0012
A
bAT1
4
3
2
1
T
T
T
T
Sistema di equazioni lineari
LT
T
0
0
0
b
98
T0=10; %temperatura nel punto x=0TL=50; %temperatura nel punto x=LL=0.5; %Lunghezza della sbarraN=4; %numero di nodi di discretizzazioneh=L/(N+1);
A=[-2 1 0 0;1 -2 1 0; 0 1 -2 1; 0 0 1 -2];
b=[-T0;0;0;-TL];
disp('il vettore incognito è:')T=inv(A)*b
disp('Vettore completo:')T_tot=[T0 T' TL]'
x=0:h:L;T_teorica=T0+(TL-T0)/L*x
plot(x,T_tot,'*',x,T_teorica,'r'),gridxlabel('Ascissa lungo la sbarra [m]')ylabel('Temperatura [°C]')
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.510
15
20
25
30
35
40
45
50
Ascissa lungo la sbarra [m]
Tem
pera
tura
[°C
]
99
Modificare il codice per includere un numeroarbitrario N di nodi soluzione.
100
T0=10; %temperatura nel punto x=0TL=50; %temperatura nel punto x=LL=0.5; %Lunghezza della sbarraN=10; %numero di nodi di discretizzazioneh=L/(N+1);
A=zeros(N,N);for i=1:N,
A(i,i)=-2;endfor j=1:N-1,
A(j,j+1)=1;A(j+1,j)=1;
end
b=zeros(N,1);b(1)=-T0;b(N)=-TL;
disp('il vettore incognito è:')T=inv(A)*b%ISTRUZIONE ALTERNATIVAT=A\b;
disp('Vettore completo:')T_tot=[T0 T' TL]'
x=0:h:L;T_teorica=T0+(TL-T0)/L*x
plot(x,T_tot,'*',x,T_teorica,'r'),gridxlabel('Ascissa lungo la sbarra [m]')ylabel('Temperatura [°C]')
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.510
15
20
25
30
35
40
45
50
Ascissa lungo la sbarra [m]
Tem
pera
tura
[°C
]
101
La soluzione in regime permanente è ancora calcolabile in forma chiusa
E’ noto e prefissato il valore di temperatura T0 all’estremo sinistro, mentre latemperatura all’estremo destro non è imposta.
Ora studiamo un problema differente
Poiché la sbarra è termicamente isolata si avrà che il flusso termico nei bordodestro è pari a zero
0
Lx
xTx
0TxT Soluzione in forma chiusa
2
43
34
2
2
4
4
h
TT
h
h
TT
h
xTx
xTx
xTx
xxLx
xx
0
102
0212 TTT
02 321 TTT
02 432 TTT
043 TT
Modificare il codice Matlab
bAT
1100
1210
0121
0012
A
0
0
0
0T
b
103
021 TT
02 321 TTT
02 432 TTT
043 TT
Ora analizziamo il caso in cui anche la temperatura all’estremo sinistro non èimposta.
Il sistema lineare non puo essere risolto numericamente percheammette infinite soluzioni
.costxTSoluzione in forma chiusa
Seguendo il medesimo ragionamento adottato in precedenza si ottiene il sistema:
bAT
1100
1210
0121
0011
A
0
0
0
0
b
104
Cosa succede se le due temperature agli estremi, ancorche’ impostedall’esterno, siano variabili nel tempo ?.
L’analisi in regime stazionario perde di significato !
2
2 ,,
x
txT
t
txT
L’equazione del calore deve ora essere considerata nella sua interezza.
Il ragionamento svolto in precedenza va opportunamente modificato.
x
0x Lx
h h h h h 5/Lh
1xx 2xx 3xx 4xx
La discretizzazione del dominio è sempre valida.
105
TtTtTtTtTt 4321 ,,,T
tT0 tTL
Definiamo il vettore che contiene le temperature incognite nei nodi soluzione, che adifferenza dal caso stazionario non sono più delle costanti, ma sono delle funzioni del tempo.
Condizioni al contorno NOTE
2
2 ,,
x
txT
t
txT
tTtT 0,0
tTtLT L,Questa tipologia di condizioni al contorno sono dette di “DIRICHLET”.
txTtT ii , 4,3,2,1i5
Lihixi
106
ixx
i
x
txT
t
txT
2
2 ,,
2
012
2
2 2,
1
h
tTtTtTtxT
xxx
2
123
2
2 2,
2
h
tTtTtTtxT
xxx
2
234
2
2 2,
3
h
tTtTtTtxT
xxx
2
34
2
2 2,
4
h
tTtTtTtxT
x
L
xx
ixx
ix
txTtT
2
2 , 4,3,2,1i
Le derivate seconde vengono approssimate mediante differenze finite, in maniera analoga aquanto fatto in precedenza ma coinvolgendo funzioni del tempo anzichè costanti
107
ixx
ix
txTtT
2
2 , 4,3,2,1i
tTtTtTh
tT 01221 2
tTtTtTh
tT 12322 2
tTtTtTh
tT 23423 2
tTtTtTh
tT L 3424 2
Sostituendo le approssimazioni alle differenze finite si ottiene un sistema di equazionidifferenziali ordinarie (ODE), di ordine pari al numero di nodi soluzioni (4 nel caso in esame).
108
tTh
tTtTh
tT 022121 2
tTtTtTh
tT 32122 2
tTtTtTh
tT 43223 2
tTh
tTtTh
tT L24324 2
Risulta conveniente separare alla destra dell’uguale le quantità incognite dalle quantità noteT0(t) e TL(t)
TN tTtTt 10 , u
ttt BuATT
TtTtTtTtTt 4321 ,,,T
2100
1210
0121
0012
2h
A
1 0
0 0
0 0
0 1
2h
B
109
Devono essere note le temperature nei nodi soluzione all’istante iniziale t=0, cioè il vettorecostante
Ora sia il coefficiente di diffusione che lo step h di discretizzazione spaziale influenzanoesplicitamente la soluzione. L’analisi in regime stazionario prevedeva la scomparsa di talicoefficienti dalle relazioni risolutive
TTTTT 0,0,0,00 4321T
Va rimarcato che trascorso un transitorio sufficientemente lungo i profili di temperatura nondipendono più dalle condizioni iniziali, e tendono ad una soluzione di “regime dinamico” lecui caratteristiche sono indipendenti dalle condizioni iniziali
110
clear allclc
L=0.5;N=20;
rho=7856;Cp=435;k=63;alfa=k/(rho*Cp);
h=L/(N+1);
A=zeros(N,N);for i=1:N,
A(i,i)=-2;endfor j=1:N-1,
A(j,j+1)=1;A(j+1,j)=1;
end
B=zeros(N,2);B(1,1)=1;B(N,2)=1;
Ac=(alfa/h^2)*A;Bc=(alfa/h^2)*B;
Tcamp=0.1;Tfinale=100;
x=linspace(h,L-h,N);t=0:Tcamp:Tfinale;numsample=length(t);
Tinit=10*ones(1,N);
T_soluz=zeros(numsample,N);T_dot=zeros(numsample,1);
T_soluz(1,:)=Tinit;
T0=30+5*sin(0.1*t)';TL=30*ones(numsample,1);
for i=1:numsample-1Tdot=Ac*T_soluz(i,:)'+Bc*[T0(i);TL(i)];T_soluz(i+1,:)=T_soluz(i,:)+Tcamp*Tdot';
end
x_ext=linspace(0,L,N+2);T_soluz_ext=[T0 T_soluz TL];
[X,Y] = meshgrid(x_ext,t);
h=mesh(X,Y,T_soluz_ext);title('Distribuzione temperatura lungo la sbarra.');xlabel('Coordinata spaziale x [m] ','FontName','times','FontSize',14);ylabel('Tempo [s]','FontName','times','FontSize',14);zlabel('T(x,t)','FontName','times','FontSize',14);set(gca,'FontSize',14,'FontName','Times');
Eulero esplicito
Cond. al contorno
Grafica
111
Profilo spaziotemporale della soluzione
112
clear allclc
L=0.5;N=5;
rho=7856;Cp=435;k=63;
alfa=k/(rho*Cp);
h=L/(N+1);
A=zeros(N,N);for i=1:N,
A(i,i)=-2;endfor j=1:N-1,
A(j,j+1)=1;A(j+1,j)=1;
end
B=zeros(N,2);B(1,1)=1;B(N,2)=1;
Ac=(alfa/h^2)*A;Bc=(alfa/h^2)*B;
Tcamp=0.01;Tfinale=100;
x=linspace(h,L-h,N);t=0:Tcamp:Tfinale;numsample=length(t);
Tinit=30*ones(1,N);
T_soluz=zeros(numsample,N);T_dot=zeros(numsample,1);
T_soluz(1,:)=Tinit;
T0=30+20*sin(0.1*t)';TL=30*ones(numsample,1);
for i=1:numsample-1Tdot=Ac*T_soluz(i,:)'+Bc*[T0(i);TL(i)];T_soluz(i+1,:)=T_soluz(i,:)+Tcamp*Tdot';
end
x_ext=linspace(0,L,N+2);T_soluz_ext=[T0 T_soluz TL];
[X,Y] = meshgrid(x_ext,t);
h=mesh(X,Y,T_soluz_ext);title('Distribuzione temperatura lungo la sbarra.');xlabel('Coordinata spaziale x [m] ','FontName','times','FontSize',14);ylabel('Tempo [s]','FontName','times','FontSize',14);zlabel('T(x,t)','FontName','times','FontSize',14);set(gca,'FontSize',14,'FontName','Times');
113
114
2
2 ,,
x
txT
t
txT
tTtT 0,0
0,
tL
x
TQuesta tipologia di condizione al contorno viene detta di “NEUMANN”.
Ora applichiamo una condizione al contorno che non “fissi” il valore della temperaturanell’estremo destro, ma imponga l’isolamento termico dell’estremo destro x=L inanalogia con quanto fatto nella analisi a regime.
Si ricava il modello discretizzato
h
tTtTtL
x
T NL
,
115
TtTt 0 , 0u
ttt BuATT
TtTtTtTtTt 4321 ,,,T
1100
1210
0121
0012
2h
A
1 0
0 0
0 0
0 1
2h
B
tTt 0u
ttt BuATT
TtTtTtTtTt 4321 ,,,T
1100
1210
0121
0012
2h
A
0
0
0
1
2h
B
116
clear allclc
L=0.5;N=20;
rho=7856;Cp=435;k=63;alfa=k/(rho*Cp);
h=L/(N+1);
A=zeros(N,N);
for i=1:N,A(i,i)=-2;
endfor j=1:N-1,
A(j,j+1)=1;A(j+1,j)=1;
endA(N,N)=-1;
B=zeros(N,2);B(1,1)=1;B(N,2)=1;
Ac=(alfa/h^2)*A;Bc=(alfa/h^2)*B;
Tcamp=0.1;Tfinale=1000;
x=linspace(h,L-h,N);t=0:Tcamp:Tfinale;numsample=length(t);
Tinit=10*ones(1,N);
T_soluz=zeros(numsample,N);T_dot=zeros(numsample,1);
T_soluz(1,:)=Tinit;
T0=10+5*sin(0.1*t)';TL=0*ones(numsample,1);
for i=1:numsample-1Tdot=Ac*T_soluz(i,:)'+Bc*[T0(i);TL(i)];T_soluz(i+1,:)=T_soluz(i,:)+Tcamp*Tdot';
end
x_ext=linspace(0,L,N+2);T_soluz_ext=[T0 T_soluz T_soluz(:,end)];
[X,Y] = meshgrid(x_ext,t);
h=mesh(X,Y,T_soluz_ext);title('Distribuzione temperatura lungo la sbarra.');xlabel('Coordinata spaziale x [m] ','FontName','times','FontSize',14);ylabel('Tempo [s]','FontName','times','FontSize',14);zlabel('T(x,t)','FontName','times','FontSize',14);set(gca,'FontSize',14,'FontName','Times');
117
k=63
118
k=630
119
Sistema LTI (Linear Time Invariant)MIMO (multi-input-multi-output)
tT tT
tAT
tBu tu
ttt BuATT
Uso di Matrix Gain Integratore saturato
Alessandro Pisano - [email protected]
120
Un processo termico 2D a parametri distribuiti
Cilindro cavo.
Parti grigie in acciaio.
Parte bianca: volume con vapore ad alta temperatura.
Parte nera: piccolo volume interno
Es. Sezione di una turbina a vapore
statore
rotorerCoordinata radiale
maxr
Alessandro Pisano - [email protected]
122
Ipotesi: simmetria angolare della distribuzione di temperatura
trT ,
Si desidera calcolare la distribuzione di temperatura nel rotore (in uno dei suoi raggi)
maxmin rrr
Mediante misure acquisite in una turbina in esercizio, si suppone nota la temperatura nella “parte bianca” (regione del vapore in alta temperatura)
Problema complesso (anche nella formulazione semplificata sotto esame) perche il modello matematico è una equazione alle derivate parziali (sistema a parametri distribuiti, sistema infinito-dimensionale)
In una modellazione più dettagliata, si potrebbe essere interessati a valutare la distribuzione di temperatura nella superficie delle pale di rotore, onde valutare gli stress termici sui materiali.
Alessandro Pisano - [email protected]
123
Eq. di diffusione (Equazione del calore) monodimensionale in coord. cilindriche , con unica variabile spaziale la coordinata cilindrica radiale r
pC
kK
è il coefficiente di diffusione [m2/s]
pC è la capacita termica a pressione costante per unita di massa [J/kg K]
è la conduttività termica [J/K s m]k
è la densità [g /m3 ]
124
minrr
maxrr
hirri min
1..., ,2 ,1 ,0 Ni
1
minmax
N
rrh
trTtT ii ,
r =0
Alessandro Pisano - [email protected]
125
Approssimazione delle derivate spaziali mediante differenze finite
h
tTtTtT
r
iii
1
2
11
2
2 2
h
tTtTtTtT
r
iiii
Sistema di ODE
tTtTr
htT
r
h
h
K
tTh
tThrh
tThrh
K
h
tTtTtT
h
tTtT
rKtT
ii
i
i
i
ii
i
i
i
iiiii
i
i
112
12212
2
111
21
11211
21
Ni ..., ,2 ,1
Alessandro Pisano - [email protected]
126
Sistema di ODE
tTtT
r
htT
r
h
h
KtT 01
1
2
1
21 21
tTtT
r
htT
r
h
h
KtT 12
2
3
2
22 21
tTtT
r
htT
r
h
h
KtT NN
N
N
N
N 11221
…
TN tTtTt 10 , u
TNN tTTTtTtTt ,...,,, 2,321 T
tT0 tTN 1 Boundary conditions ttt BuATT
Alessandro Pisano - [email protected]
127
N
N
r
h
r
h
r
h
r
h
r
h
r
h
r
h
h
K
210000
1000
0100
00210
000121
000012
1
3
22
11
2
A
ttt BuATT
Nr
h
h
K
10
00
00
01
2
B
TNN TTTTT 0,...,,0,00 12,321 T Condizioni iniziali
Alessandro Pisano - [email protected]
128
Modello Simulink
Due diverse modalità di specificare i due segnali di input.
Linee spesse che rappresentano segnali multidimensionali
Contenuto del Subsystem “Modello termico Rotore”
FILES
TurbinaVapore01.mdl
T_centrorotore.mat
T_vapor.mat
Esportazione dati verso il Workspace
Alessandro Pisano - [email protected]
132
Matrici A e B (per n = 6)
n=40;
r_min=0.05;
r_max=0.392;
rho=7900;
C=0.45;
k=73;
T_in=80;
h=(r_max-r_min)/n;
K=k/(rho*C*1e3)/h^2;
A=zeros(n,n);B=zeros(n,2);
for i=1:n,
A(i,i)=-(2+h/(r_min+i*h));
end
for j=1:n-1,
A(j,j+1)=1+h/(r_min+j*h);
A(j+1,j)=1;
end
B(1,1)=1;B(n,2)=1+h/r_max;
A=K*A; B=K*B;
Codice copiabileed seguibile
134
Profili temporali degli elementi del vettore T Uscita del blocco Scope “T”
Analisi dei risultati
Impostiamo preliminarmente due valori costanti per le temperature al contorno
CtTtTtT
N
300
30,0u
136
Con dei profili differenti per le boundary conditions, l’evoluzione del profilo di temperatura èdifferente.
Ora processiamo in Matlab i risultati della simulazione, creando dei grafici 3D.
r_min=0.05;
r_max=0.392;
n=40;
[X,Y] = meshgrid(linspace(r_min,r_max,n),tout);
h=mesh(X,Y,Trot)
title('Distribuzione temperatura rotore.')
xlabel('Coordinata radiale r [m]
','FontName','times','FontSize',14)
ylabel('Tempo [s]','FontName','times','FontSize',14)
zlabel('T(r,t)','FontName','times','FontSize',14)
set(gca,'FontSize',14,'FontName','Times')
Alessandro Pisano - [email protected]
137
Invertiamo la direzione dell’asse dei tempi
138
Completiamo questo esempio mostrando l’impiego, nella Mask, di variabili di configurazionetipo “popup” o “checkbox” e un loro possibile impiego.
Si apportino le seguenti modifiche alla lista dei Parameters
139
Si apportino le seguenti modifiche alle istruzioni di Initialization
140
Si apportino le seguenti modifiche allo schema Simulink
Ora si esegua il modello con diverse scelte per le variabili pop up e checkbox, e si analizzino irisultati.
FILES:
TurbinaVapore02.mdl
T_centrorotore.mat
T_vapor.mat
Alessandro Pisano - [email protected]
141
Anti-lock bracking system (ABS)
Simulazione di un sistema di frenatura con ABS Modello quarto di veicolo
psrf FRTI = velocità ruota
I = inerzia ruota= forza da interazione
pneumatico/stradapsF
= coppia frenantefTpsFvm
4
smg
Fps 4
rRv
s/
1
s = scorrimento
Rr = raggio ruota
smg
RTI rf 4
smg
vm
44
Caratteristica -s nota in forma tabellare
Sistema dinamico NON LINEARE, del secondo ordine
Alessandro Pisano - [email protected]
142
Modello SIMULINK smg
RTI rf 4
smg
vm
44
clear all
g = 9.81;
v0 = 40; % velocita iniziale (m/s)
Rr = 0.6; % raggio
m = 200; % massa
J = 5; % inerzia
% Curva mu-slip
slip = 0:.05:1.0;
mu = [0 .4 .8 .97 1.0 .98 .96 .94 .92 .9 .88 .855 .83 .81 .79 .77 .75 .73 .72 .71 .7];
Blocco Fcn per il calcolo dello scorrimento
Lookup Table
Stop a veicolo fermo
143
Lookup Table (libreria Lookup Tables)
La caratteristica viene disegnata sullamaschera del blocco
Blocco Stop Simulation(libreria Sinks )
Interrompe la simulazione quando il veicolo si ferma
Integratori saturati
Alessandro Pisano - [email protected]
145
Soluzione semplificata per interrompere la simulazione quando
il veicolo si ferma
Integratori saturati
Alessandro Pisano - [email protected]
146
Ora modelliamo l’ABS
La dinamica del sistema di frenatura idraulico è approssimata da un filtro del primo ordine(dinamica cassetti distributori) e da un integratore saturato (pressurizzazione).Il controllo è un relè (controllo bang-bang). L’apposito blocco Relay si trova nella libreriaDiscontinuities
Kf = 3; %guadagno
PBmax = 1500; % saturazione
TB = 0.01; %costante di tempo
Alessandro Pisano - [email protected]
147
Modello complessivo
E’ stato aggiunto un blocco “Manual Switch” (libreria Signal Routing) per poter simulare unafrenatura non controllata (Tf=-1000).Eseguiamo una simulazione, e visualizziamo i risultati nei blocchi Scope
FILES:
ABS_OpenClosedLoop.mdl
ABS_OpenClosedLoop_DATI.m
Alessandro Pisano - [email protected]
148
Frenatura non controllata (ABS disattivato)
Il pneumatico si blocca dopo 1.5 secondi circa, mentre la marcia del veicolocontinua per altri 4 secondi.
Lo scorrimento diventa unitario quando il pneumatico si blocca
149
Frenatura con ABS attivato
Il pneumatico ora si blocca solo nell’ultima parte della frenata.
Il bloccaggio del pneumatico avviene quando ormai la velocita di marciadel veicolo è prossima a zero, la marcia del veicolo si arresta infatti dopopochi decimi di secondo.
Lo scorrimento viene regolato attorno al set point desiderato 0.2.diventando unitario solo a frenata ormai conclusa.
Nella prossima slide si confrontano due diversi test del sistema ABS conuna diversa velocita di marcia iniziale V0.
Alessandro Pisano - [email protected]
150
Con ABS
v0 = 40 m/s 145 km/h
Con ABS
v0 = 70 m/s 250 km/h
151
Embedded Matlab Function block
Consideriamo nuovamente il sistematermico del secondo ordine
tTtTC
Ktq
CtT fp
f
ip
f
f 1
tTtTC
KtTtT
C
KtT fp
p
ip
pe
p
pe
p
tT f tTe tq
tTptemperatura della parete
Alessandro Pisano - [email protected]
152
FILE: termico_2dof_matfile_mask01_EMF.mdl
Realizziamo il modello SIMULINK in maniera differente, generando i segnali q e Te con deiblocchi Signal Builder (libreria Sources) …..
Alessandro Pisano - [email protected]
153
….. e realizzando le equazioni per mezzo di un blocco Embeddbed MATLAB Function (EMF)
Integratore“vettoriale”
(bidimensionale)
y = usciteu = ingressip = parametri
154
Il blocco Camera di Combustione conserva la medesima maschera dellaimplementazione precedente
Alessandro Pisano - [email protected]
155
Codice del blocco EMF
function ydot = fcn(y,u,p)
Cf=p(1);
Cp=p(2);
Kip=p(3);
Kpe=p(4);
q=u(1);
Te=u(2);
Tf=y(1);
Tp=y(2);
Tfdot=q/Cf+(Kip/Cf)*(Tp-Tf);
Tpdot=(Kpe/Cp)*(Te-Tp)-(Kip/Cp)*(Tp-Tf);
ydot=[Tfdot Tpdot];Variabili globali non concesseall’interno di blocchi EMF.
I parametri della mask non sonodirettamente accessibili da partedella EMF
Codice analogo a quellodi un Function file
Alessandro Pisano - [email protected]
157
Zoom on T
Spostiamo verso sinistraquesta linea(dopppio click su di essa)
Alessandro Pisano - [email protected]
158
Spostiamo verso destraquesta linea con modalitàanaloghe
Spostiamo verso l’alto questalinea con modalita analoghe
Alessandro Pisano - [email protected]
159
Profilo completo q(t) Profilo più complicato per q(t)
Proflo piu complicato per T(t)
Alessandro Pisano - [email protected]
160
Non serve piu un file script in abbinamento
Risultati della simulazione
Alessandro Pisano - [email protected]
161
SimMechanics
Modellazione di sistemi meccanici multi-body
Toolbox avanzati
Alessandro Pisano - [email protected]
164
Variare il controllo sul FUEL (Subsystem Control)
Osservare la velocità dell’albero motore in RPM, e correlarne le variazioni ai corrispondenti incrementi e decrementi del fuel rate
Alessandro Pisano - [email protected]
165
Toolbox avanzati
SimDriveline
Componenti e modelliarea automotive
Animazione 3D
Alessandro Pisano - [email protected]