Simulazione Del Modello Dinamico Di Un Serbatoio.pdf

CAPITOLO 5 INTRODUZIONE Nel seguente capitolo saranno descritti alcuni esempi di modellazione di sistemi che possono fare parte degli impianti navali. Dopo un’analisi del sistema si propongono diverse soluzioni, di tipo lineare e non lineare, di tipo analitico o numerico.

SIMULAZIONE DEL COMPORTAMENTO DINAMICO DI UN SERBATOIO Efflusso libero Come primo esempio di simulazione si propone il sistema idraulico mostrato in fig.5.1, si supponga di studiare il comportamento di un serbatoio con efflusso libero, senza essere limitati da ipotesi semplificative di linearizzazione richiesta dai sistemi LTI. Per studiare il comportamento dinamico di questo sistema si parte dal principio di conservazione della massa, per cui: QE − QU = A⋅dh/ dt (1) dove A è la sezione del serbatoio, supposta costante. Dal principio di Torricelli si ha che la velocità media di efflusso dovuta ad un carico h è:

V = k 2 gh da questa relazione deriva la portata delle bocche a battente:

QU = µ ⋅ AB 2 gh

Fig.5.1

(2)

dove AB è la sezione della bocca e µ è il coefficiente di efflusso che dipende dalla forma e dalla sezione della bocca, dalle caratteristiche del fluido e da h: per bocche circolari e per l’acqua questo coefficiente si scosta di poco da 0.6 anche per ampie variazioni di h. Per la simulazione si veda il caso numerico seguente: h(0)=2 [m], livello iniziale acqua AB =5.07e-4 [m2], corrispondente ad un tubo da 1” A= 2 [m2], sezione serbatoio µ= 0.6 L’espressione (2), ponendo:

k B = µ ⋅ AB 2 g =1.347e-3, diventa: QU = k B ⋅ h (2’) L’equazione (1) diventa:

dh 1 = QE − k B h (3) dt A

(

)

Fig.5.2

E’ un’equazione differenziale non lineare che si risolve facilmente con Simulink, il cui modello è mostrato in fig.5.2a. Come si vede il sistema ha una retroazione naturale propria. La portata iniziale in uscita sarà: QU(0)= 1.905 e-3 [m3/s], supponendo di avere una portata QE costante, se questa è superiore a QU(0) il sistema non è stabile, il livello cresce sempre più. In fig.5.2b si vedono le curve del livello e della portata in uscita nel caso di QE= 1 e-3. Facendo girare il programma per il caso QE= 0, ossia per la risposta libera, si verifica che il serbatoio si svuota in un tempo di circa 4000 secondi.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

1

Modello lineare Per rendere il sistema lineare, si ipotizza che la portata in uscita sia funzione del battente h e della resistenza R provocata dalla strozzatura della valvola: QU= h /R (4) L’equazione differenziale (1) diventa:

dh QE h = − dt A RA

(5)

Il modello Simulink di questa equazione è mostrato in fig.5.3a. Questo modello lineare è valido solo per un limitato intorno di h, questa semplificazione è accettabile per un controllo di livello, per cui si può linearizzare la curva nei pressi del valore di riferimento, per esempio supponendo di voler lavorare intorno al valore iniziale di h(0)=2 [m], si può calcolare R per dare la stessa portata calcolata nella forma non approssimata: QU(0)= 1.905 e-3=h(0)/R ⇒ R=1.05e+3 [s/m2], la costante di tempo è τ=RA≈ 2100[s]. In fig.5.3b è mostrata la risposta con le stesse condizioni iniziali e portata QE del modello precedente, si vede che il sistema risponde in modo molto diverso perché opera in un intervallo troppo elevato di h. Dall’espressione (5) si ricava la funzione di trasferimento del serbatoio:

GV ( s ) =

H ( s) R 1.05e3 = = QE ( s ) R ⋅ A ⋅ s + 1 2.1e3 ⋅ s + 1

Fig.5.3

(6)

Fig.5.4 L’ingresso è la portata QE e l’uscita il livello h del serbatoio. Come si sa, il modello che utilizza la FdT non può considerare che condizioni iniziali nulle. L’uso della FdT rende però agevole la sintesi del regolatore. Il modello Simulink che utilizza la FdT del serbatoio è schematizzato in fig.5.4, dal punto di vista della simulazione non è possibile neanche descrivere variabili interne, come la portata in uscita. La regolazione del livello Il sensore più classico è costituito da un galleggiante che misura il livello del serbatoio ed è collegato all’asta di apertura della valvola mediante un sistema di leve, come si vede in fig.5.5. L’attuatore è realizzato con una valvola proporzionale, supponendo lineare la legge tra spostamento della valvola a e la portata: QE = KV⋅a , dove KV è la costante della valvola. Il sistema di leve serve a posizionare la valvola in funzione dell’errore sul livello, misurato dal galleggiante, per cui: QE = KP⋅(hR – h) (7)

Fig.5.5

Fig.5.6

dove hR è il livello di riferimento, in sostanza esso è il nuovo ingresso del sistema e la costante KP tiene conto di tutte le altre (quelle della valvola e delle leve). L’insieme valvola-galleggiante ha delle inerzie che si possono assimilare ai ritardi di un sistema del primo ordine, la cui funzione di trasferimento complessiva sarà:

Gg ( s ) =

QE ( s ) K g KV a KP = = H ( s ) Tg ⋅ s + 1 0.2 ⋅ s + 1

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

(8)

2

Dove nella costante proporzionale KP si sono conglobate tutte le altre. Come caso numerico si è ipotizzata una costante di tempo Tg=0.2 [s]. Potendo Simulink risolvere anche le equazioni differenziali non lineari, si utilizzerà per la vasca il modello più vicino al vero, ossia quello di fig.5.2a. L’ingresso del sistema è ora il livello di riferimento hR, lo schema a blocchi del regolatore proporzionale di livello è mostrato in fig.5.6, dove si è imposta un’azione proporzionale KP=1. Il sistema , nel caso lineare era di tipo uno, per cui esiste un errore statico che si riduce aumentando il guadagno, ma con il rischio di sovraoscillazioni, questo caso era stato visto e discusso nel Corso di Condotta Automatica. Serbatoio con portata dipendente dall’utenza esterna Questo è un caso più frequente nella pratica rispetto a quello precedente. Esso non è più un sistema SISO (Single input-Single Output), ma MISO (Multiple input-Single Output), ma non ci sono problemi usando Simulink. In questo caso la portata QU non è più una variabile dipendente da h, ma varia con legge casuale a secondo della

Fig.5.8 richiesta dell’utenza, da zero ad un valore massimo imposto. Essa si può quindi assimilare ad un disturbo esterno cui il regolatore deve opporsi. Per questo sistema si ritorna all’espressione (1), da cui, applicando Laplace si ottiene:

H ( s) =

QE − QU s⋅ A

(9)

La portata sulla valvola è espressa dalla (8) e dalla legge del galleggiante (7), per cui si ha: QE = (HR−H)⋅Gg (10) Per simulare una variazione di portata si è usato un blocco di libreria Simulink chiamato “Step”, che crea un segnale a gradino di ampiezza assegnata partendo da un certo tempo tS. Nel caso in esame si è impostata un’ampiezza QU=3 litri/s=3e−3 m3/s e tS =5s. In fig.5.8a è mostrato lo schema a blocchi del sistema che è diventato lineare. La fig.5.8b mostra la serie temporale del livello e della portata in uscita. Come si vede, esiste un piccolo errore statico (0.3 mm), la cui presenza non si sarebbe dedotta per via dell’integratore, questo però non opera sull’errore ma sulla ∆Q, sull’errore opera invece il blocco Gg(s) di tipo zero, caratterizzato da un errore statico non nullo. Controllo ad azione proporzionale ed integrale Per eliminare l’errore statico occorre aggiungere un’azione legata all’integrale dell’errore, lo schema di fig.1.9a

Fig.5.9 mostra un tipico esempio e la fig.5.9b mostra le serie temporali, dove si nota che l’errore tende a zero. Si è aggiunto un polo nell’origine ed uno zero: z =− KI , i valori numerici del regolatore sono KI=50 e KP=100, le oscillazioni sono estremamente piccole ed il transitorio si esaurisce in meno di due secondi. La fig.5.10 mostra una realizzazione pratica di un controllo di livello di un serbatoio con un attuatore idraulico (o pneumatico) che ha un’azione di tipo integrativa mentre il sistema di leve che aziona il distributore è di tipo proporzionale, il tutto è quindi un regolatore meccanico-idraulico PI. © Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

Fig.5.10 3

SIMULAZIONE DI UN PENDOLO SEMPLICE E DEL ROLLIO DI UNA NAVE IL PENDOLO Le oscillazioni di un pendolo, se di modesta ampiezza, sono un fenomeno periodico molto regolare, gli orologi meccanici funzionano su questo principio. Galileo Galilei, intuì la caratteristica essenziale dei moti pendolari nel 1583, probabilmente su osservazioni dell'oscillazione di un lampadario nel duomo di Pisa. Egli poté verificare che la durata di quelle lente oscillazioni rimaneva immutata nonostante la loro ampiezza diminuisse sempre di più. I suoi numerosi esperimenti lo condussero a formulare il principio dell'isocronismo delle piccole oscillazioni del pendolo: per spostamenti di qualche grado di escursione, l'oscillazione è indipendente dall'ampiezza. Galileo estese poi, erroneamente, questo principio, anche ad oscillazioni di qualsiasi ampiezza . Secondo le esperienze di Focault, il pendolo tende a conservare immutato il proprio piano di oscillazione, rispetto a un sistema di riferimento inerziale (stelle fisse). Poiché la Terra ruota, un osservatore ad essa solidale vede ruotare il piano di oscillazione del pendolo rispetto al suolo terrestre, di una quantità (velocità angolare) che dipende dalla latitudine del luogo in cui si compie l'esperimento. Il moto del pendolo Il moto di un pendolo è un classico esempio di moto armonico. Si chiama pendolo semplice un sistema formato da una massa m, sospesa a un filo rigido di lunghezza l e di massa trascurabile fissato a un estremo O, detto centro di sospensione. Una volta allontanato dalla posizione di riposo il pendolo, come si vede in fig.5.11, oscilla tra i due punti della traiettoria circolare di centro O e raggio l che corrispondono all'elongazione massima α (con 0 < α < π/2 ): ovvero, nel riferimento polare con origine O ed asse OX, tra il punto G di coordinate polari (l,−α) ed il punto G’ di coordinate (l,+α)). Siano θ la coordinata angolare di P in un generico istante (−α≤θ≤α) e ω=dθ/dt la sua velocità angolare e d ω/dt= d2θ/dt2, la sua accelerazione angolare. A P è associato un vettore velocità avente in ogni istante direzione tangente alla traiettoria e modulo v = ω⋅l. L'equazione del moto del pendolo semplice si trova applicando il teorema della conservazione della quantità di moto. Imponendo che il momento (rispetto ad O) delle forze applicate sia uguale alla derivata rispetto al tempo del momento (sempre rispetto ad O) della quantità di moto, si ha:

OP ∧ (T + mg ) =

d OP ∧ mv dt

(

)

(1)

Dove T è la tensione del filo, che è diretta lungo OP ed ha dunque momento nullo rispetto ad O; nella configurazione descritta dalla fig.5.11 il momento delle forze agenti è diretto normalmente al piano di figura, ha verso che si avvicina all'osservatore e modulo:

m⋅g⋅l⋅senθ

Fig.5.11

La velocità v del pendolo è sempre perpendicolare ad OP

quindi il momento della quantità di moto ha modulo:

m⋅v⋅l = m⋅l2ω e direzione perpendicolare al piano di figura, mentre il verso della sua derivata è diretto lontano dall'osservatore. Sostituendo nella (1) queste espressioni, si ricava l'equazione del moto del pendolo:

d 2θ g + senθ = 0 dt 2 l

(2)

Questa è un’equazione differenziale non lineare del secondo ordine di non facile soluzione. Per piccole oscillazioni si può assumere che senθ ≈ θ , per cui si ottiene il sistema lineare:

d 2θ g + θ =0 dt 2 l

(3)

si è operata una linearizzazione, questa equazione è nota come equazione differenziale del moto armonico.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

4

Equazione differenziale del moto armonico Un'equazione differenziale del tipo:

d 2x + ω2x = 0 2 dt

(3’)

è nota in analisi come equazione differenziale del moto armonico, la cui soluzione generale è del tipo: x(t)= x(0)⋅cos(ω⋅t) (4) l’ampiezza delle oscillazioni è data dalla condizione iniziale. La (4) corrisponde ovviamente ad un moto periodico di periodo:

T= L’espressione (3) si riduce alla (3’) ponendo:

2π

ω

ω2 =

g l

Nei limiti della sua validità, il moto del pendolo semplice è quindi un moto armonico di periodo:

T0 = 2π

l g

(5)

Come si può notare il periodo di oscillazione è indipendente dalla massa e dall’ampiezza α delle oscillazioni, per cui si può affermare che le piccole oscillazioni del pendolo sono isocrone. Gli orologi meccanici funzionano su questo principio. La trattazione fin qui fatta non tiene conto delle perdite di energia dovute agli attriti ed alla resistenza dell’aria, quindi le oscillazioni non sono smorzate come nella realtà. Gli orologi necessitano dell’energia accumulata in una molla per sopperire alle perdite per attrito. Periodo di oscillazione per grandi escursioni Per grandi escursioni α di oscillazione occorre impiegare l’equazione differenziale non lineare (2), la cui soluzione numerica risulta essere molto elaborata. Si può dimostrare che il periodo di oscillazione è dato dalla seguente serie di Taylor:

T0 = 2π

l g

  1 1 32 α 2α sen sen 4 + ..  (6) 1 + +  22 2 2 2 2 4 2  

Simulazione del pendolo semplice con Simulink Simulink si presta molto bene anche per i sistemi non lineari, per cui si non si adotterà il modello linearizzato ma l’equazione (2), che si riscrive ponendo a sinistra il termine di ordine più elevato:

g l

θ ′′ = − senθ Si supponga di volere un periodo di oscillazione T0=1s, si utilizzi in prima approssimazione l’espressione (5) per ricavare la lunghezza :l= g⋅ (T0/2π)2= 0.2485 m, per cui: g/l= (2π)2= 39.48 [s−2]. Per le condizioni iniziali, si supponga θ’(0)=0 e θ(0)= 3°= 0.05236 rad. In fig.5.12a si vede il semplice modello Simulink e nella fig.1.12b l’oscillazione del pendolo, si osservi che il periodo è pari ad un secondo e che la funzione è proprio quella dell’espressione (4). Simulando con una condizione iniziale di 60°=1.047 rad/s, il sistema è fortemente non lineare, per cui il periodo è sensibilmente aumentato, come si vede nella fig.1.2c, usando i primi tre termini della serie di Taylor si ricava T0≈1.07 s. Per rendere il sistema più realistico occorre introdurre uno smorzamento viscoso, per cui l’equazione diventa:

g l

θ ′′ = − senθ − bθ ′ Fig.5.12

La fig.5.12d mostra il modello corrispondente Simulink.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

5

SIMULAZIONE DEL MOTO DI ROLLIO DI UNA NAVE In prima approssimazione si può assimilare il moto di rollio di una nave in acqua calma ad un pendolo semplice. Indicando con ρ il raggio d’inerzia della nave rispetto all’asse longitudinale e con ∆ il dislocamento, il momento d’inerzia sarà:

J=

∆ 2 ρ g

Indicando con θ l’inclinazione trasversale della nave, ossia l’angolo di rollio, il momento MI dovuto alle forze d’inerzia è:

M I = Jθ ′′ =

∆ 2 ρ θ ′′ (7) g

Indicando con GZ il braccio del momento di stabilità e con GM l’altezza metacentrica, il momento Mr raddrizzante sarà: M r = ∆ ⋅ GZ (θ ) (8) Per piccoli angoli si può ipotizzare che la retta d’azione della spinta, Fig.5.13 passante per B, passi anche per il metacentro trasversale M, come si vede in fig.5.13, per cui: M r ≈ ∆ ⋅ GM ⋅ sen(θ ) (9) Per piccoli angoli si può confondere il seno con l’angolo, a questo punto l’espressione diventa lineare: M r ≈ ∆ ⋅ GM ⋅ θ (9’) Eguagliando a zero la somma dei momenti (7) e (9) si ottiene una espressione del tutto simile a quella del pendolo semplice (2):

∆

θ ′′ +

g ⋅ GM

ρ2

senθ = 0 (10)

Indicando con B lo smorzamento dovuto all’acqua, e considerando la forma più precisa (8) per il momento raddrizzante, si ottiene l’equazione che governa il moto di rollio: J ⋅ θ ′′ + B ⋅ θ ′ + ∆ ⋅ GZ (θ ) = 0 (11) Passando a Simulink, come si è già visto, conviene mettere al primo membro il termine con il grado più elevato, quindi:

B J

θ ′′ = − θ ′ −

∆ ⋅ GZ (θ ) J

(11’) Fig.5.14

Il modello Simulink è visibile in fig.5.14, l’angolo θ della tavola è espresso in gradi, mentre è sempre bene che il modello utilizzi variabili del Sistema Internazionale, quindi occorre un blocco di conversione da radianti a gradi. In fig.5.15 è mostrato il diagramma di stabilità della nave GZ(θ), ricavato dalla tavola. Esso è stato calcolato per una nave di piccole dimensioni (20 m di lunghezza e 4 di larghezza), i dati numerici utilizzati per la simulazione sono: ∆ = P = 4.172825e+5 [N] peso della nave, ρ = 1.494 [m] raggio d’inerzia, B = 3700 [kg⋅m2/s] smorzamento idrodinamico, GM = 0.44 [m] altezza metacentrica, da questi dati si ricava il momento d’inerzia della nave: J = 9.5e+4 [kg⋅m2]. Per piccoli angoli dall’equazione differenziale del moto armonico si ricava la frequenza di oscillazione:

ω=

g ⋅ GM

ρ

2

=

g ⋅ GM

ρ

Fig.5.15

(12)

Passando ai valori numerici, si ottiene: ω=1.39 [rad/s] cui corrisponde una frequenza f= 0.221 [Hz] ed un periodo di rollio pari a 4.52 [s]. Per lanciare la simulazione occorre un piccolo programma che assegni le costanti, in particolare la tabella dei bracci di stabilità, il cui listato è il seguente:

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

6

% programma rollsimr.m di simulazione del rollio con Simulink % con modello rollior.mdl % G.Carrera 7 febbraio 2003 Tsim = 200; % durata simulazione J = 9.5e+4; % momento di inerzia della nave [kg*m^2] P = 4.172825e+5; % peso della nave [N] B = 3.7e+3; % smorzamento idrodinamico [kg*m^2/s] % tavola dei bracci di stabilità teta = linspace(-60,60,25); % vettore teta, 25 valori tra +/-60° teta0 = 0; % angolo iniziale di inclinazione trasversale GZ = [-0.07,-0.115,-0.16,-0.2,-0.225,-0.24,-0.235,-0.205,-0.17,-0.13,-0.085,-0.04,... 0.0,0.04,0.085,0.13,0.17,0.205,0.235,0.24,0.225,0.2,0.16,0.115,0.07]; plot(teta,GZ),grid,title('Diagramma di stabilità'), xlabel('inclinazione trasversale'),ylabel('Braccio del momento di stabilità'); teta0= 3; % angolo [°] d'inclinazione iniziale teta0=teta0*pi/180; % conversione in radianti rollior; % chiama il programma simulink di simulazione

Le costanti di integrazione sono θ’=0 e θ=3°, la fig.5.16 mostra la risposta (al rollio) libera della nave in acqua tranquilla, partendo da uno sbandamento iniziale teta0= 3° e velocità angolare nulla. Dalla figura si contano circa 9

Fig.5.16

Fig.5.17 periodi nei primi 40 secondi, per cui risulta un periodo di rollio di 4.44 secondi, che è abbastanza vicino a quello calcolato con il modello linearizzato (T0=4.52). Per angoli iniziali molto maggiori, ad esempio con teta0=60° (ammesso che la nave non affondi), dalla fig.5.17 si possono contare 13 cicli in circa 60 s, per cui il periodo diventa paria a 4.62 s. Utilizzando questo modello ed incrementando teta0 si arriva alla instabilità intorno ai 68°, naturalmente non si considerano qui le altre cause che provocherebbero molto prima un affondamento della nave (imbarchi d’acqua, etc). Se la nave è soggetta a momenti forzanti, indotti per esempio da un’accostata o da spostamenti di carico, l’instabilità si raggiunge molto prima. Questo modello risulta quindi assai utile per lo studio della stabilità. Bibliografia 1. Maurizio Loreti , “Note sul pendolo”, Dipartimento di Fisica - Università degli Studi di Padova, 1999-2000. 2. David Halliday, Robert Resnick, “Fisica generale”, vol.I, Casa Editrice Ambrosiana Milano, 1967 3. A. Cavallo, R. Setola, F. Vasca, “Guida operativa a MATLAB SIMULINK E CONTROL TOOLBOX”, Liguori Editore, 1994. 4. Antonio Fiorentino, “Fondamenti di Automazione Analogica e Numerica”,Edizioni Scientifiche Italiane, Napoli 1981.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

7

SIMULAZIONE DI UN IMPIANTO DI RISCALDAMENTO Il modello Simulink di simulazione dell'impianto è quello illustrato in fig.5.19, l'escursione giorno/notte della temperatura esterna è stata simulata con un segnale sinusoidale con ampiezza di 5°C attorno alla temperatura media di +15°C. ON L'elemento non lineare del sistema è il termostato, questo, come si vede in fig.5.18, è un dispositivo a relè a due soglie. Nel caso illustrato si sono poste le soglie a ±1°C sul segnale di errore, l'uscita del termostato è: Terr • ON = 1 bruciatore acceso dall'istante in cui Ti = 19°C (Terr= +1°C) fino OFF +1°C −1°C all'istante in cui Ti = 21°C (Terr= −1°C). • OFF = 0 bruciatore spento dall'istante in cui Ti = 21°C (Terr= −1°C) fino Fig.5.18- Termostato all'istante in cui Ti = 19°C (Terr= +1°C). Nella fig.5.21 si vede l'andamento della temperatura interna Ti (variabile controllata) e quella esterna To, avendo supposto che questa vari in modo sinusoidale tra il giorno e la notte. Come si vede nel particolare di fig.5.22, l'accensione del riscaldatore dell'aria porta ad un rapido aumento della temperatura (la costante di tempo del riscaldamento è molto minore di quella di raffreddamento), al raggiungimento della soglia superiore (21°C) viene spento il bruciatore. L'aria viene lentamente a raffreddarsi (secondo la classica curva di scarica esponenziale: e−t/τ), per via del calore che, attraverso le pareti e le finestre della casa, viene ceduto all'esterno. Alla temperatura corrispondente alla soglia inferiore (19°C), il termostato scatta ed accende il bruciatore. Come si vede, la temperatura media è quella desiderata, ma si hanno continue oscillazioni. In questo modello si è trascurata la costante di tempo dello stesso termostato, essendo molto inferiore a quella del sistema da controllare. Riducendo il valore di soglia si riducono le ampiezze delle oscillazioni ma si aumenta la loro frequenza con

Fig.5.19

Fig.5.21

Fig.5.20

Fig.5.22

conseguenti frequenti accensioni e spegnimenti del bruciatore ed aumento dell'usura delle parti meccaniche interessate (elettrovalvole, etc). Quando la calderina è accesa (l'uscita del termostato è =1) ci sarà un flusso termico: QC= portata aria calda [kg/sec] × entalpia aria a 40° [joule/kg]= 0.1 × 3.114 ×105 = 31140 [joule/sec] I valori numerici sono quelli relativi ad un’abitazione di ampie dimensioni, essi sono quelli calcolati dal programma thermdat.m allegato. Il modello matematico della casa si ricava dal bilancio energetico:

QC − QE = Ct dove

QE =

dTi Q − QE , applicando Laplace: QC − QE = Ct ⋅ s ⋅ Ti ⇒ Ti = C dt Ct ⋅ s

Ti − To è il flusso termico ceduto all'esterno, Req è la resistenza termica equivalente delle pareti e delle Req

finestre, Ct è la capacità termica. Come si vede questo è un tipico sistema di ordine uno. Nella fig.5.20 si vede lo schema a blocchi del modello termodinamico della casa che soddisfa le equazioni precedenti. Questo sistema è MISO perché ha due ingressi: QC e To ed una uscita Ti. La capacità termica è il prodotto: © Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

8

Ct= massa aria × calore specifico a pressione costante = M⋅cp = 1.788×106 [joule/°C]

I tempi di riscaldamento sono molto inferiori a quelli di raffreddamento se la casa è sufficientemente isolata termicamente (Req alta). Per averne un’idea, si può calcolare facilmente il tempo che impiega il bruciatore a scaldare l'aria di 2 gradi (dalla soglia inferiore alla superiore) nel caso in cui Ti−To sia piccola (QE << QC), per cui si ha:

trisc= Ct⋅∆T/QC = 1.788×106 [joule/°C] ×2[°C]/ 31140 [joule/sec]= 114.83 [sec]

To

1 Req

1 Ct s

Ti

Quando la calderina è spenta QC=0 per cui la funzione di trasferimento è SISO ed il sistema di fig.5.20 si riduce a quello fig.5.20b, da cui si ricava facilmente la funzione di trasferimento ad anello chiuso:

Fig.5.20b

Ti 1 = To sReq Ct + 1

per cui la costante di tempo per la fase di raffreddamento vale:

τ = Req⋅ Ct = 1.5371×10−3 [°C× sec/ joule] × 1.788×106 [joule/°C]= 2748.3 [sec] ≈ 45’48’’ Di seguito è visibile il programma thermdat.m che viene pre-caricato (pre-load) da Simulink al momento del lancio del modello termo.mdl. Per fare questo si procede nel seguente modo: dal menu file del desktop di Simulink si sceglie Model properties , quindi Callbacks. Nel campo Model pre-load function si inserisce il nome del file script (senza estensione *.m). Questo è un utile espediente per fare partire la simulazione dal file modello anziché dal file di lancio, ma nella gran parte dei casi sono proprio le costanti della simulazione a variare e non il modello, per cui è preferibile il sistema del file di lancio da Matlab. % THERMDAT % This script runs in conjunction with the "thermo" % house thermodynamics demo. % Copyright 1990-2000 The MathWorks, Inc. % $Revision: 1.13 $ % Ned Gulley 4-9-92 % ------------------------------% Problem constant % ------------------------------r2d = 180/pi; % ------------------------------% Define the house geometry % ------------------------------% House length = 30 m lenHouse = 30; % House width = 10 m widHouse = 10; % House height = 4 m htHouse = 4; % Roof pitch = 40 deg pitRoof = 40/r2d; % Number of windows = 6 numWindows = 6; % Height of windows = 1 m htWindows = 1; % Width of windows = 1 m widWindows = 1; windowArea = numWindows*htWindows*widWindows; wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ... 2*(1/cos(pitRoof/2))*widHouse*lenHouse + ... tan(pitRoof)*widHouse - windowArea; % ------------------------------% Define the type of insulation used % ------------------------------% Glass wool in the walls, 0.2 m thick kWall = 0.038; LWall = .2; RWall = LWall/(kWall*wallArea); % Glass windows, 0.01 m thick kWindow = 0.78; LWindow = .01; © Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

9

RWindow = LWindow/(kWindow*windowArea); % ------------------------------% Determine the equivalent thermal resistance for the whole building % ------------------------------Req = RWall*RWindow/(RWall + RWindow); % c = cp of air (273 K) = 1005.4 J/kg-K c = 1005.4; % ------------------------------% Enter the temperature of the heated air % ------------------------------% THeater = 30 deg C = 303 K THeater = 303; % ha = enthalpy of air at 40 deg C = 311400 J/kg ha = THeater*c; % Air flow rate Mdot = 0.1 kg/sec Mdot = 0.1; % ------------------------------% Determine total internal air mass = M % ------------------------------% Density of air at sea level = 1.2250 kg/m^3 densAir = 1.2250; M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse*lenHouse)*densAir; % ------------------------------% Enter the cost of electricity and initial internal temperature % ------------------------------% Assume the cost of electricity is $0.09 per kilowatt/hour % Assume all electric energy is transformed to heat energy % 1 kW-hr = 3.6e6 J % cost = $0.09 per 3.6e6 J cost = 0.09/3.6e6; % TinIC = initial indoor temperature = 20 deg C TinIC = 20; Il programma calcola il volume della casa e la resistenza termica equivalente delle pareti e delle finestre, la massa totale dell’aria e le varie costanti termiche. Nel modello di simulazione è anche conteggiato il consumo ed il costo complessivo, noto il costo a kWh di energia termica. È interessante e molto educativo vedere come sale sensibilmente il consumo al crescere della temperatura di riferimento anche di uno o due gradi. Modello lineare Per passare al modello lineare basta eliminare il termostato e sostituirlo con un regolatore proporzionale, come si vede in fig.5.23. Essendo il sistema di tipo zero, esisterà un errore statico, con un guadagno Kp=10 esso è di circa 0.01 °C, come si vede in fig.5.24, molto inferiore alla precisione del sensore di temperatura che in questo caso occorre utilizzare (il termostato faceva Fig.5.23 le funzioni di sensore e dispositivo di soglia). Nel caso lineare si richiede anche un attuatore proporzionale che dosi il calore generato dalla calderina in modo continuo e lineare. Questa funzione è realizzata, per impianti medio- piccoli, per mezzo di elettro-valvole che variano la portata del combustibile al bruciatore. Per impianti di maggiori dimensioni l’attuatore è idraulico Fig.5.24 o pneumatico. Normalmente per questo tipo di impianti è sufficiente usare il sistema a termostato, che ha elementi più semplici e robusti.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

10

SIMULAZIONE DI UN ATTUATORE IDRAULICO Modello matematico Rispetto agli altri, gli attuatori idraulici sono quelli che possono trasmettere le forze o le coppie più elevate con ingombri limitati, essi hanno quindi una potenza specifica molto elevata. Per queste caratteristiche, unite alla prontezza di risposta ed alla accuratezza, hanno trovato sempre più largo impiego a partire dalla fine della seconda guerra mondiale. Le loro prime applicazioni a bordo delle navi sono state in campo militare per gli asservimenti di puntamento dei cannoni e dei missili. La fig.5.25 mostra un tipico esempio di attuatore idraulico, il pistone del martinetto è del tipo a doppio effetto, pilotato da una Fig.5.25 valvola di comando del tipo proporzionale realizzata con un distributore a cassetti. Un piccolissimo spostamento x dell'asta di comando produce uno spostamento y del pistone con notevole forza impegnata. Si supponga di avere un carico elastico, per esempio una struttura in acciaio da sollecitare. L'uscita sia la posizione y del pistone, ossia la freccia imposta al carico e l'ingresso la posizione x della valvola di comando. Indicando con A1 ed A2 le sezioni a sinistra ed a destra del pistone e P1 e P2 le rispettive pressioni, la forza agente sul pistone è: FP = P1 A1 − P2 A2 Indicando con M la massa del pistone e con b la costante di attrito viscoso, l'equazione della dinamica del pistone sarà:

y′′ =

1 ( P1 A1 − P2 A2 − by& − ky ) M

(1)

In un sistema idraulico, l’equazione di continuità, in termini di portate volumetriche, è espressa come:

∑ Q −∑ Q in

out

=

dV V dp + dt β dt

Dove Qin e Qout sono le portate entranti ed uscenti, V è il volume, p la pressione e β =−Vdp/dV è il modulo di compressibilità volumetrica del fluido. Partendo da questa relazione applicata al circuito di fig.1, supponendo nulle le perdite per trafilamento: Qout=0, sapendo che dV/dt =Ady/dt, si ottengono le seguenti equazioni:

p1′ = p2′ =

β V1

β V2

(Q1 − A1 y& )

(2)

(Q2 − A2 y& )

(3)

Con V1= V01+y⋅A1 (4) V2= V02−y⋅A2 (5) dove V01 e V02 sono i volumi delle due zone a pistone centrato (posizione di riposo). Per i fluidi usati negli impianti idraulici β ≈ 1.5169×109 [N/m2], questo è il valore che sarà usato nella simulazione. Le equazioni che permettono di ricavare la portata in funzione, delle pressioni e della posizione x della valvola derivano dall'equazione delle strozzature per regimi di moto turbolento:

Q = C ⋅ A( x ) 2 ⋅

∆p

ρ

dove: C = coefficiente dell'orifizio, A = area dell'orifizio (dipende da x), ρ = densità dell'olio, ∆p = caduta di pressione attraverso l'orifizio. Considerando i movimenti verso destra (x ≥ 0) o verso sinistra (x < 0) dell'asta della valvola del distributore, si hanno due espressioni per ogni portata:

  x ⋅ C f 1 ( x ) Ps − P1 Q1 =   x ⋅ C f 1 ( x ) P1 − Pr

per x ≥ 0 per x < 0

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

(6)

11

 − x ⋅ C f 2 ( x ) P2 − Pr  Q2 =   − x ⋅ C f 2 ( x ) Ps − P2

per x ≥ 0 per x < 0

(7)

Cf sono coefficienti che inglobano anche le altre costanti, Ps è la pressione della mandata della pompa e Pr quella di ritorno. Queste equazioni non sono lineari ed i coefficienti non sono costanti ma dipendono dalla posizione x della valvola di comando. Occorre quindi fare qualche semplificazione, prima di tutto Pr≈0 perché si scarica a pressione atmosferica, per piccole variazioni si possono considerare indipendenti da x i coefficienti d'orifizio per le simmetrie della valvola, per cui le equazioni diventano:

 x ⋅ C f 1 Ps − P1 per x ≥ 0 Q1 =  per x < 0  x ⋅ C f 1 P1 per x ≥ 0  − x ⋅ C f 2 P2 Q2 =   − x ⋅ C f 2 Ps − P2 per x < 0

(6')

(7')

Modello Simulink L'equazione (1) è facilmente modellabile, essendo già nella forma canonica, il risultato è quello mostrato in fig.5.26. Questo modello può essere inserito in un blocco che ha come ingresso le pressioni P1 e P2 e come uscite la posizione y del pistone attuatore e la sua velocità y'.

Fig.5.26 - Dinamica del pistone Ora si passa ai modelli Simulink delle equazioni (2), (3), (4) e (5) . Anche qui non ci sono stati molti problemi, la fig.5.27 mostra il modello corrispondente alle equazioni (2) e (4) e la fig.5.28 quello relativo alle equazioni (3) e (5). Per passare alle pressioni occorrerà inserire in uscita un blocco integratore..

Fig.5.27

Fig.5.28

Ora non mancano che le portate, esse si ricavano dalle equazioni semplificate (6') e (7'), per valutarle in funzione del segno dello spostamento si usano i blocchi switch, la fig.5.29 mostra il modello relativo alla portata Q1 e la fig.5.30

Fig.5.29 – Modello valv_q1 Fig.5.30 Modello valv_q2 quello per la portata Q2. Il blocco switch ha il compito di inviare in uscita il segnale all'ingresso superiore se l'ingresso di controllo (segnale centrale), è positivo, altrimenti manda in uscita il segnale inferiore. © Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

12

Il modello della pompa è molto semplice: si ipotizza una pompa a portata costante con una perdita di carico tale da realizzare una pressione PS di circa 1800 psi, il modello è quello di fig.5.31. La fig.5.32 mostra lo schema a blocchi Simulink dell'attuatore idraulico completo di retroazione di posizione e di blocchi di monitoraggio. Il Fig.5.31 - Modello della pompa. regolatore è di tipo PID e alla sua uscita si nota un blocco ritardo che simula la dinamica della servovalvola con un sistema di ordine uno, all’uscita vi è un blocco saturatore per impedire che l’ingresso della valvola distributrice sia eccessivo per errori grandi. Il segnale

Ps P1 Q1 x

beta

A1

valv_q1 ydot

P1dot

Gain1

Gain A1

Valvola di comando y

A1

P1

Integrator

Product pressioni

Gain2 V01 Portate

1 s

area A1 A2

vol 1

area A2

Ps

Gain6

P2 Q2 pompa

x

beta

A2

valv_q2

Gain4

Gain3 A2

psi

Gain5

-K-

P2dot

Saturation 1 ritardo

Product1

1 s

P2

1 s

y

Integrator 2 Integrator 3

smorzamento b

ydot molla

Integrator 1

V02 vol 2

1 s

1/M

k

Dinamica del pistone PISTONE

y

1 0.03s+1

PID

Sine Wave

Ps

Spostamenti

PID

1797

errore

Fig.5.32 - Modello di attuatore idraulico retroazionato in spostamento di riferimento è una sinusoide di 0.33 Hz di frequenza e di 15 mm di ampiezza con una costante aggiuntiva (bias) di 20 mm. I sottosistemi "valv_q1" e "valv_q2" sono rispettivamente quelli mostrati in fig.5.29 e fig.5.30. Gli altri blocchi sono stati collegati in modo aperto, in modo da dare una più immediata visione del funzionamento del sistema. La difficoltà maggiore per simulare questi sistemi sta nel trovare i giusti coefficienti della valvola distributrice, che sono difficili da reperire nella bibliografia o dai costruttori. Si potrebbero confrontare i dati simulati con quelli sperimentali, ma anche qui non è facile misurare i piccoli spostamenti del cassetto distributore e le pressioni nelle immediate vicinanze degli orifizi. Si propongono quindi dei valori molto approssimativi, necessari per provare la simulazione, il programma di lancio è il seguente: % % % %

******************************************************** programma idr_att.m di lancio del modello Simulink att_idraul.mdl G. Carrera re. 15 genn 2003 ********************************************************

% pompa: Ps= 1800; % pressione tipica impianto idraulico in [psi] Ps= Ps*6895; % in [Pa] Qp=0.00183 ; % portata (costante) della pompa in [m^3/s] Cp= 1.4771e-10; % perdita di carico della pompa % Regolatore Kp= 10; % costante proporzionale % Costante elastica del carico: k= 3e4; % in [N/m] % Costante smorzamento: b= 0.2; % in [N*s/m] © Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

13

% Costanti cilindro: A1 = 0.0123; % sezione del pistone [m^2] D=0.125 m (lato senza asta) A2 = 0.0123-0.00332; % sezione del pistone [m^2] (lato asta, D=0.065 m ) yt = 0.3; % corsa del pistone [m] V01=A1*yt/2; % volume camera 1 a pistone cetrale V02=A2*yt/2; % volume camera 2 a pistone cetrale beta= 1.5169e9 ; % bulk modulus [Pa] olio M= 10; % massa [kg] complessiva pistone, asta e carico % Costanti valvola di comando: Cd1= 0.61; % costante strozzatura Cd2= Cd1; rho= 800; % densità olio [kg/m^3] circuito idraulico Cf1= Cd1*sqrt(2/rho); Cf2= Cd2*sqrt(2/rho); % lancia la simulazione: att_idraul Le dimensioni del cilindro e la pressione di alimentazione sono quelle dell’impianto MTS del Laboratorio Strutture del DINAV. Bibliografia 1. Umut Genç, 1 Keith Glover, Richard Ford - "Nonlinear Control Of Hydraulic Camshaft Actuators In Variable Cam Timing Engines" - Department of Engineering, University of Cambridge, UK 2. Salvador Esqué, Jose LM Lastra - "Modeling and simulation of a servopneumatic gripper" - 10 Dec 99.

© Simulazione degli impianti navali – G. Carrera – Cap.5– rev.27/04/08

14