Modellazione e implementazione numerica di un motore a turbina per elicottero (2012-2019)

Indice

Inquadramento e finalità

1.1

Propulsione negli elicotteri

1.1.1

Componenti tipici del motore per elicotteri

1.1.2

Ciclo termodinamico

1.1.3

Componenti connessi al propulsore

1.2

Finalità del lavoro di tesi

Stato dell’arte

2.1

Modellazione di motori a turbina

2.2

Software esistenti

Fisica del problema e modellazione

3.1

Modellazione delle turbomacchine

3.2

Proprietà dei fluidi

3.3

Componenti del modello

3.3.1

Presa d’aria

3.3.2

Compressore

3.3.3

Plenum

3.3.4

Combustore

3.3.5

Turbina generatrice di gas

3.3.6

Turbina di potenza

3.3.7

Modello di carico semplificato

3.3.8

Dinamica degli alberi di trasmissione

3.4

Dinamica del sistema e metodo di soluzione

3.5

Parametrizzazione

3.5.1

Esempio di adattamento delle mappe per un generico

motore

Implementazione

4.1

Requisiti del codice

4.2

Genesi e struttura del codice Matlab

INDICE

4.2.1

Modello γ - Matlab

4.2.2

Modello ISO - Matlab

4.2.3

Funzioni di supporto

4.3

Struttura del codice C

4.3.1

Funzione Principale: EngineModel . .

4.3.2

Funzioni di supporto

4.4

Metodo di soluzione e di integrazione

4.5

Interfaccia grafica

4.6

DLL - Dinamic Link Library

4.7

GSL - Gnu Scientific Libraries

Validazione

5.1

Confronto con dati sperimentali a regime

5.2

Validazione dinamica: risposta ad un gradino di portata di

combustibile

Controllore automatico del propulsore

6.1

Architettura e logica di controllo

6.1.1

Limiti di impiego del motore

6.2

Implementazione

6.2.1

Implementazione in linguaggio Matlab

6.3

Simulazioni del motore in anello chiuso

6.3.1

Simulazione 1 - controllore proporzionale

6.3.2

Simulazione 1 - controllore PI

6.3.3

Simulazione 1 - controllore PID

6.3.4

Simulazione 2 - controllore proporzionale .

101

6.3.5

Simulazione 2 - controllore PI

103

6.3.6

Simulazione 2 - controllore PID

104

Simulazione di una manovra

107

7.1

Illustrazione della manovra e software impiegati .

107

7.2

Confronto utilizzando alcuni parametri del motore GE t700

109

7.3

Confronto utilizzando parametri identificati

112

Conclusioni e sviluppi futuri

117

8.1

Conclusioni

117

8.2

Sviluppi futuri

118

Elenco delle figure

1.1

Schema semplificato di motore turbo-albero connesso a tra-

smissione e rotore

1.2

Schema di ciclo joule Brayton ideale nel piano T-S

1.3

Ciclo Joule-Brayton reale nel piano T-S

2.1

Il motore General Electric t700

3.1

Variazione della costante dei gas R in funzione della tempera-

tura per i gas combusti di kerosene

3.2

Calore specifico dei gas ottenuti dalla combustione di kerosene

in funzione della temperatura per diversi rapporti aria/combustibile. 31

3.3

Rapporto fra calori specifici del kerosene in funzione di tem-

peratura e rapporto aria/combustibile

3.4

Rapporto fra calori specifici k dell’aria e dei gas combusti di

kerosene con rapporto aria/combustibile di

ottenuto con il

modello di Guha

3.5

Calore specifico c_p dell’aria e dei gas combusti con rapporto

aria/combustibile di

, ottenuto con il modello di Guha per

un ampio intervallo di temperature

3.6

Calore specifico di aria e gas combusti ottenuto con modello

di Walsh, con FAR=1/35

3.7

Schema del motore e flusso di informazioni fra i vari componenti.

3.8

Pala di rotore semplificata

3.9

Potenza richiesta dal carico modellato in funzione della veloci-

ta di rotazione n₂ confrontata con i valori di potenza e velocità

di rotazione nominali per il motore GE t700

3.10

Diagramma di flusso per il processo di generazione e adatta-

mento delle mappe delle turbomacchine

3.11

Rapporti di compressione/espansione normalizzati

3.12

Mappa normalizzata della portata d’aria elaborata dal com-

pressore (g2) in funzione del regime di rotazione n₁

ELENCO DELLE FIGURE

3.13

Mappa normalizzata del rapporto di espansione della turbi-

na di potenza in funzione della temperatura in uscita dalla

turbina HPT

3.14

Mappa dei β_compressore e β_HPT adattata ad un generico motore.

3.15

Mappa (adattata ad un ipotetico motore) della portata d’a-

ria elaborata dal compressore (g₂) in funzione del regime di

rotazione n₁

3.16

Mappa (adattata ad un ipotetico motore) del rapporto di espan-

sione della turbina di potenza in funzione della temperatura

in uscita dalla turbina HPT

4.1

Schema del codice γ implementato in Matlab

4.2

Mappa della temperatura in uscita dal compressore in funzione

del regime di rotazione n₁ nel modello γ per il motore GE t700.

4.3

Schema delle funzioni implementate nel codice ISO

4.4

Funzione EngineModel.c : ingressi e uscite

5.1

Andamento dell’ingresso di portata di combustibile a gradino.

5.2

Risposta ad un ingresso a gradino di portata di combustibile:

confronto del regime di rotazione n₁ con i dati del report Ballin

- 1988

5.3

Risposta ad un ingresso a gradino di portata di combustibile:

confronto della pressione all’uscita dal compressore p₂₁ con i

dati del report Ballin - 1988

5.4

Risposta ad un ingresso a gradino di portata di combustibi-

le: confronto dellla temperatura T₃ in uscita dalla camera di

combustione con i dati del report Ballin - 1988

5.5

Risposta ad un ingresso a gradino di portata di combustibile:

confronto dellla temperatura T₄ in uscita dalla turbina di alta

pressione con i dati del report Ballin - 1988

6.1

Architettura di un sistema di controllo di un motore a turbina

per elicotteri

6.2

Logica di controllo del regolatore PID implementato all’inter-

no del FADEC

6.3

Sistema motore-controllore in anello chiuso

6.4

Struttura del controllore FADEC nell’implementazione in lin-

guaggio Matlab

6.5

Simulazione 1 con controllore P. Velocità di rotazione dei due

alberi e variabile di controllo g_b

ELENCO DELLE FIGURE

6.6

Simulazione 1 con controllore Proporzionale.Velocità angolare

dell’albero della turbina di potenza. Confronto con il riferi-

mento da inseguire

6.7

Inseguimento del riferimento di velocità di rotazione costante

con controllore PI

6.8

Simulazione 1 con controllore PI. Andamento dei due sta-

ti meccanici e della variabile di controllo durante il primo

secondo di simulazione

6.9

Velocità di rotazione n₁, n₂ e portata di combustibile g_b duran-

te l’inseguimento del riferimento costante n₂

con controllore

ref

di tipo PI

6.10

Simulazione 1 con controllore PID. inseguimento del riferimen-

to sulla velocità di rotazione n₂

6.11

Simulazione 1 con controllore PID. Velocità di rotazione e

portata di combustibile durante la simulazione

100

6.12

Simulazione 1 con controllore PID. Particolare del primo se-

condo di simulazione

100

6.13

Inseguimento del riferimento da parte del sistema motore -

controllore in anello chiuso. Simulazione 2. Controllore solo

proporzionale

102

6.14

Velocità di rotazione dei due alberi e flusso di carburante in

camera di combustione durante la simulazione 2 di insegui-

mento di un riferimento variabile di n₂ con motore controllato

da un regolatore proporzionale

102

6.15

Simulazione 2 con controllore PI. Inseguimento del riferimento

variabile n₂

103

REF

6.16

Simulazione 2 con controllore PI. Andamento dei due stati

meccanici n₁, n₂ e della variabile di controllo g_b

104

6.17

Simulazione 2 con controllore PID. Inseguimento del riferimen-

to variabile di velocità di rotazione n₂

105

6.18

Simulazione 2 con controllore PID. Storia temporale dei due

stati meccanici n₁ ed n₂ e della variabile di contollo g_b durante

la simulazione

106

7.1

Storia temporale della portata di combustibile g_b durante la si-

mulazione effettuata in Flighlab e poi riprodotta con il modello

ISO

109

7.2

Confronto fra la velocità di rotazione n₁ calcolata dal modello

ISO e quella ottenuta dalla simulazione in Flightlab. Si riporta

anche il valore nominale per il motore considerato. Parametri

incogniti derivati dal motore GE t700

110

ELENCO DELLE FIGURE

7.3

Confronto della temperatura T₄ dei gas fra le due turbine cal-

colata con il modello ISO e il software Flightlab.Parametri

incogniti derivati dal motore GE t700

111

7.4

Potenza espressa dal modello ISO confrontata con quella espres-

sa dal modello FlightLab. Per i parametri incogniti si è fatto

ricorso a quelli del motore GE t700

112

7.5

Confronto della temperatura T₄. Modello ISO con parametri

modificati e simulazione FlightLab

114

7.6

Confronto sulla velocità n₁. Simulazione effettuata con soft-

ware FlightLab e con modello ISO con parametri modificati. . 114

7.7

Confronto sulla potenza erogata all’interfaccia con il carico.

Simulazione condotta con software FlighLab e con modello

ISO con parametri modificati

115

Elenco delle tabelle

3.1

Coefficienti dei polinomi per il calcolo di k e c_p dei gas combusti

nel modello di Walsh [1]

3.2

Dati ricavati da test NASA su motore GE t700 (prima parte).

3.3

Dati ricavati da test NASA su motore GE t700 (seconda parte).

5.1

Errori percentuali dei due modelli rispetto al caso 2 di trim

sperimentale

5.2

Errori percentuali dei due modelli rispetto al caso 3 di trim

sperimentale

5.3

Errori percentuali dei due modelli rispetto al caso 4 di trim

sperimentale

5.4

Errori percentuali dei due modelli rispetto al caso 5 di trim

sperimentale

5.5

Errori percentuali dei due modelli rispetto al caso 6 di trim

sperimentale

5.6

Condizioni iniziali di temperatura nelle varie sezioni del mo-

tore prima dell’ingresso a scalino di portata di combustibile.

5.7

Condizioni iniziali delle portate prima dell’ingresso a scalino

di portata di combustibile

5.8

Condizioni iniziali delle pressioni prima dell’ingresso a scalino

di portata di combustibile

6.1

Taratura dei guadagni con il metodo di Ziegler e Nichols

6.2

Condizioni iniziali della simulazione 1 di motore controllato in

anello chiuso da controllore proporzionale

7.1

Dati del motore P&W PT6C-67C

108

Sommario

In questo lavoro viene proposto un modello numerico di motore per elicotte-

ro semplificato, flessibile e adattabile a più applicazioni attraverso il minore

numero di parametri possibile. Tale modello è stato implementato in diverse

versioni prima in linguaggio Matlab e poi in linguaggio C attraverso lo svi-

luppo di una DLL, allo scopo di essere poi incluso in codici di simulazione e

di meccanica del volo. La modellazione termodinamica è stata validata con-

frontando i risultati ottenuti da alcune simulazioni con dei dati sperimentali

disponibili in letteratura.

É stato poi affiancato al modello di motore anche

un controllore che ne regoli il funzionamento, seguendo una logica simile a

quella di un FADEC. Infine si è validato il codice simulando una manovra

tipica e confrontando i risultati con quelli di riferimento.

The object of this thesis work is the development of a simplified, flexible,

customizable helicopter engine numerical model, adaptable to several existing

engines by tuning a small set of parameters. The model has been implemen-

ted both in Matlab and C computing languages, resulting in a DLL that can

be easily linked with Flight Mechanics simulation codes. The thermodynamic

model has been validated comparing the simulation results with available ex-

perimental data. An automatic controller for the engine has been developed

and integrated to reproduce typical FADEC logics. Finally, the model has

been tested using real inputs from a typical manoeuvre and comparing the

results with reference data.

Capitolo 1

Inquadramento e finalità

1.1

Propulsione negli elicotteri

La propulsione tramite motori a turbina a gas è quella preponderante negli

elicotteri di dimensioni medie e grandi. Il motore a turbina per elicotteri è

di tipo turboalbero (turboshaft): tutta la potenza che si riesce a estrarre dal

fluido di lavoro viene trasferita attraverso ad un’albero di trasmissione che

muoverà i carichi ad esso connessi; i carichi sono costituiti principalmente

dal rotore principale, dal rotore di coda e da alcune utenze secondarie degli

impianti del velivolo. I motori turboalbero sono simili ai propulsori turboe-

lica utilizzati in velivoli ad ala fissa, ma a differenza di questi ultimi non è

necessario progettare il motore in modo che il flusso di gas di scarico conservi

all’uscita dal motore un’energia cinetica importante perchè non se ne sfrut-

ta il contributo di reazione nella spinta, bens`ı tutta la potenza che fornisce

il propulsore viene trasferita tramite l’albero di potenza alla trasmissione e

quindi alle utenze.

1.1.1

Componenti tipici del motore per elicotteri

I componenti di un tipico motore a turbina per elicotteri sono riportati in

figura 1.1, essi sono:

• la presa d’aria,

• il compressore,

• la camera di combustione,

• la turbina generatrice di gas,

• la turbina di potenza.

CAPITOLO 1. INQUADRAMENTO E FINALITÀ

Figura

1.1: Schema semplificato di motore turbo-albero connesso a

trasmissione e rotore.

La presa d’aria è il componente che cattura l’aria esterna per convogliarla

all’interno del motore. Essa può essere di due tipi:

• statica: se l’aria fluisce all’interno del motore solo per l’azione di aspi-

razione effettuata dal propulsore.

• dinamica: se l’aria che fluisce nel propulsore lo fa anche sotto l’effetto

della velocità dovuta al moto dell’elicottero.

Gli obiettivi di una buona presa d’aria dinamica sono di:

• aumentare la pressione, per evitare parte del lavoro di compressione al

compressore;

• realizzare un flusso uniforme a monte del compressore;

• ottenere minime perdite di pressione totale, e quindi minimizzare gli

attriti;

• provocare la minima resistenza aerodinamica esterna;

• avere minimo peso e lunghezza.

Dalla presa d’aria, l’aria è convogliata verso il compressore, dove l’a-

ria subisce una compressione che la porta ad una pressione e temperatura

più elevate. Il compressore è una macchina operatrice pneumofora, ovve-

ro una macchina che innalza la pressione di un gas mediante l’impiego di

energia meccanica. I compressori utilizzati nei motori per elicotteri sono

principalmente di due tipi:

1.1. PROPULSIONE NEGLI ELICOTTERI

• assiale, dove l’aria fluisce parallelamente rispetto all’asse longitudinale

del motore. I compressori assiali sono composti di più stadi, ognuno

formato da un elemento solidale al supporto del motore, lo statore, e

da un elemento rotante, il rotore.

• centrifugo, in cui l’aria, dopo essere stata raccolta dalle palette di una

girante, è compressa per azione della forza centrifuga associata al moto

di rotazione delle palette del compressore.

In molti motori elicotteristici si realizzano soluzioni ibride in cui alcuni sta-

di del compressore sono realizzati tramite moduli assiali, mentre altri con

moduli centrifughi.

L’aria in uscita dal compressore viene poi diretta verso la camera di com-

bustione attraverso un diffusore atto a rallentare il flusso e creare le condi-

zioni migliori per l’accensione e la stabilizzazione della fiamma. In camera di

combustione l’aria è miscelata con il combustibile (generalmente kerosene),

e la miscela risultante viene bruciata liberando l’energia necessaria all’ali-

mentazione del propulsore. La conformazione delle camere di combustione

per elicotteri è generalmente di tipo anulare con flusso inverso, ovvero essa

corre per tutta la circonferenza del propulsore e la direzione del flusso d’aria

è ruotata di 180^◦, consentendo un minore ingombro del propulsore. L’aria in

arrivo dal diffusore del compressore viene divisa in due flussi:

• un flusso primario dove l’aria, miscelata con il combustibile all’interno

del tubo di fiamma, raggiunge temperature dell’ordine dei 2500^◦C;

• un flusso secondario, in cui l’aria è più fredda e che consente di non

danneggiare i materiali con cui è realizzata la camera di combustione.

Il combustibile viene iniettato in camera di combustione attraverso de-

gli iniettori che si dividono in primari, i quali consentono la combustione

nelle normali situazioni operative, e secondari, che provvedono attraver-

so degli igniter anche all’accensione del motore. La combustione avviene

con grande eccesso di aria, generalmente in un propulsore per elicotteri il

rapporto aria/combustibile è dell’ordine di 35:1. All’uscita dalla camera di

combustione i gas combusti prodotti hanno una temperatura dell’ordine di

1500 K.

Successivamente, i gas combusti provenienti dalla camera di combustione

vengono inviati nella prima turbina, la turbina di alta pressione (o turbi-

na generatrice di gas [2]), che ha il compito di estrarre dal fluido l’energia

necessaria ad alimentare il compressore. La turbina generatrice di gas e il

compressore sono collegati rigidamente da un albero, che generalmente viene

indicato come albero della turbina generatrice di gas. La turbina è composta

CAPITOLO 1. INQUADRAMENTO E FINALITÀ

da più stadi a loro volta formati da uno statore e un rotore. Lo statore ha

i compiti di deflettere il flusso cos`ı che l’aria investa ogni paletta del rotore

con l’incidenza adeguata, e di convertire parte dell’energia del flusso ad alta

pressione ed alta temperatura in energia cinetica. Ogni rotore è composto

da un albero e di alcune ruote sulle quali sono fissate le palette, le quali

sfruttano l’energia cinetica del flusso per mantenere in rotazione l’albero.

A valle della turbina generatrice di gas si trova la turbina di potenza (o

turbina libera), anch’essa a uno o più stadi, che sfrutta l’energia residua dei

gas di scarico per fornire la potenza necessaria a mantenere in rotazione l’al-

bero di potenza. Essendo la turbina di potenza collegata alla trasmissione

e quindi al rotore, che opera in regimi di rotazione pressoché costanti, essa

deve mantenere una velocità di rotazione assegnata, nonostante la potenza

del carico ad essa collegata sia variabile in funzione delle condizioni di volo

e dei comandi. Il mantenimento della velocità di rotazione della turbina di

potenza è possibile grazie alla variazione del numero di giri della turbina ge-

neratrice di gas. Il collegamento fra le due turbine è solo di tipo pneumatico,

in quanto non vi è nessun organo meccanico che le colleghi. Le turbine la-

vorano a regimi di rotazione compresi fra i 20000 e i 60000 giri/minuto. La

velocità di rotazione della turbina generatrice di gas è indicata generalmente

con N1 (o Ng), mentre quella della turbina di potenza con N2 (o N_f o N_PT ).

A valle delle turbine i gas combusti sono scaricati esternamente tramite un

ugello.

Il propulsore aeronautico è l’unica fonte di potenza nella maggior parte dei

velivoli, in special modo negli elicotteri, dove il requisito di leggerezza è molto

stringente: per questo motivo oltre a sottrarre potenza all’albero tramite ge-

neratori elettrici/alternatori che alimentino l’impianto elettrico, è necessario

anche spillare una quantità della portata d’aria elaborata dal compressore per

alimentare l’impianto pneumatico e il sistema di controllo ambientale (ECS)

per la pressurizzazione, se è prevista, e il condizionamento. Generalmente

tale spillamento viene effettuato nei primi stadi del compressore, quelli a più

bassa pressione.

La maggior parte dei motori per elicotteri è dotato di due alberi: il primo

(albero della turbina generatrice di gas) collega il compressore alla turbina di

alta pressione (HPT: High Pressure Turbine o turbina generatrice di gas), il

secondo (albero della turbina di potenza) è quello di potenza che collega mec-

canicamente la seconda turbina (PT: Power turbine o turbina di potenza o

turbina libera) al carico tramite la trasmissione. Il collegamento fra la prima

e la seconda turbina, come già detto, è solo di tipo pneumatico: ciò comporta

la possibilità per il motore che i due alberi ruotino a velocità diverse. Il se-

condo albero infatti, essendo collegato al rotore tramite la trasmissione, deve

girare a una velocità pressoché costante (variabile al massimo di pochi punti

1.1. PROPULSIONE NEGLI ELICOTTERI

Combustione

Espansione

Compressione

Figura 1.2: Schema di ciclo joule Brayton ideale nel piano T-S.

percentuali in alcune particolari situazioni di volo), mentre l’albero del com-

pressore e della turbina generatrice di gas, oltre a girare a velocità più elevate

(nell’ordine dei 60 000 RPM), deve essere in grado di variare la sua velocità

per mantenere costante quella del secondo albero a fronte di variazioni della

potenza richiesta dal carico.

1.1.2

Ciclo termodinamico

Analizzando nel dettaglio il ciclo termodinamico del motore a turbina per

applicazioni aeronautiche, che è il ciclo Joule-Brayton, esso si compone di

una compressione, una combustione, un’espansione, che avvengono rispetti-

vamente nel compressore, in camera di combustione e in turbina. Essendo

il ciclo Joule-Brayton un ciclo aperto, poichè il fluido di lavoro è l’aria del-

l’atmosfera, esso si chiude idealmente al di fuori del motore attraverso il

ripristino delle condizioni dell’aria da parte dell’atmosfera. Nel caso ideale

(figura 1.2) le compressioni/ espansioni sono considerate isoentropiche, ov-

vero senza aumento di entropia durante la trasformazione, mentre si assume

che la combustione avvenga senza perdite energetiche. Nella realtà, e dun-

que nell’implementazione di un modello fedele ad essa, si considerano delle

trasformazioni reali nelle turbomacchine, che non sono più isoentropiche, e

un processo di combustione affetto da perdite di pressione (figura 1.3).

CAPITOLO 1. INQUADRAMENTO E FINALITÀ

2is

4is

Figura 1.3: Ciclo Joule-Brayton reale nel piano T-S.

1.1.3

Componenti connessi al propulsore

L’importanza di garantire la potenza necessaria al rotore e la necessità di

mantenere la velocità del rotore costante in ogni situazione impone ai pro-

pulsori per elicotteri l’utilizzo di un sistema di controllo che ne governi il

funzionamento; tale sistema di controllo, dapprima realizzato con tecnolo-

gie analogiche o miste analogico-digitali, oggi è costituito dal FADEC (Full

Authority Digital Engine Control), sistema di controllo digitale che governa

tutte le funzionalità del propulsore e che ha autorità assoluta sul suo fun-

zionamento. L’implementazione di un FADEC reale è un’attività complessa

e costosa che viene effettuata molto spesso dalle stesse ditte produttrici dei

propulsori, che per ragioni di segreto industriale non divulgano le logiche che

sono alla base dell’implementazione del controllore.

Nelle applicazioni elicotteristiche di elicotteri medio/grandi i motori in-

stallati sono usualmente più di uno, ed è quindi fondamentale l’accoppiamen-

to fra più propulsori e la loro gestione. Alcuni degli studi più tipici di dina-

mica del volo di elicotteri sono le condizioni OEI (One Engine Inoperative),

o la cosiddetta Categoria A (si veda a riguardo il riferimento [3]).

Numerosi sono i sistemi e gli impianti che del motore regolano il funzio-

namento (gestione delle temperature, afflusso del carburante, lubrificazione,

anti-ghiaccio, sistemi di accensione), per un’analisi accurata di questi aspetti

si è fatto riferimento all’ottimo [4], che tratta con rigore i sistemi connessi ai

moderni motori turbofan commerciali, ma nell’implementazione del modello

non è stato ritenuto necessario entrare nel dettaglio di questi componenti.

1.2. FINALITÀ DEL LAVORO DI TESI

1.2

Finalità del lavoro di tesi

La simulazione della dinamica del volo di un aeromobile, specialmente se esso

è ad ala rotante, è profondamente influenzata dal sistema propulsivo e dal suo

funzionamento, in particolar modo durante i transitori e nel volo manovrato.

La disponibilità di un modello di propulsore che sappia riprodurre fedelmente

la risposta di un motore reale all’interno di un simulatore per lo studio della

meccanica del volo è un aspetto fondamentale.

Questa tesi si pone l’obiettivo di mettere a punto un modello di motore

semplificato, generale e adattabile con pochi aggiustamenti a vari propulsori

esistenti. Il codice che ne risulta sarà impiegato in simulatori di meccani-

ca del volo, sia all’interno del Dipartimento di Ingegneria Aerospaziale del

Politecnico di Milano che presso la divisione Flight Mechanics della ditta

AgustaWestland. Dovrà essere possibile adattare il codice al maggior nume-

ro di propulsori a turbina per elicottero, cos`ı da rendere agevole il confronto

tra l’istallazione di diversi motori su di una stessa macchina. Il codice dovrà

anche essere eseguito all’interno di simulatori di meccanica del volo e quindi

sarà necessario che soddisfi un requisito real-time (non stretto) e sarà neces-

sario mettere a punto un’interfaccia grafica per rendere agevole la prova del

codice in modalità stand-alone.

Riassumendo gli obiettivi del codice saranno:

• adattabilità a molti motori tramite pochi parametri fondamentali, fa-

cilmente reperibili;

• requisito real-time (non stretto);

• scrittura in linguaggio C per integrazione in sistemi più estesi (Simula-

tori di MdV (Meccanica del volo), ecc..);

• interfaccia grafica per prove di motore stand-alone.

Capitolo 2

Stato dell’arte

2.1

Modellazione di motori a turbina

La disponibilitá di grande potenza di calcolo ha reso possibile in tempi re-

centi il proliferare di studi sulla simulazione del funzionamento di motori

e generatori a turbina a gas, sia per il progetto del propulsore/generatore,

sia per l’integrazione dello stesso all’interno di un sistema piú complesso. I

due bisogni, quello dello studio del propulsore in sé e quello dell’integrazione

all’interno di un sistema piú vasto (ad esempio un velivolo o un impianto

industriale), hanno per forza di cose una diversa esigenza di dettaglio nel-

la simulazione: nel primo caso si deve essere in grado di riprodurre tutti

i fenomeni inerenti il motore con grande precisione cos´ı da poter vagliare

le differenti scelte progettuali, nel secondo invece la richiesta é quella di ri-

produrre il funzionamento del componente propulsore/generatore agli effetti

esterni, ovvero nell’interazione che esso ha con il sistema principale, sia esso

aeroplano, elicottero o impianto industriale.

Il nostro lavoro si riconduce alla seconda tipologia di simulazione, in par-

ticolare la modellazione del motore a turbina che abbiamo messo a punto ha

lo scopo di fornire ad un qualunque programma di simulazione di meccanica

del volo i parametri del motore di interesse per la simulazione della dinamica

del velivolo.

La grande maggioranza dei lavori riguardanti simulazioni di motori a tur-

bina per elicotteri, citando il lavoro di Ballin [5] e quello di Duyar e Gu [6], ba-

sano i loro modelli sulle prove effettuate dai laboratori NASA Lewis negli anni

’80 sul motore General Electric t700 (figura 2.1), largamente diffuso e uti-

lizzato, tra gli altri sull’elicottero Sikorsky UH-60 (due unità per elicottero),

sull’AgustaWestland AW149 (due unità per elicottero) e sull’AgustaWestland

AW101 (tre unità per elicottero).

CAPITOLO 2. STATO DELL’ARTE

Figura 2.1: Il motore General Electric t700.

Un lavoro di riferimento per ogni simulazione di motore per elicotteri è il

memorandum tecnico [5], dove sulla base di esperimenti compiuti dalla NASA

su un motore General Electric t700, si mette a punto un modello semplificato

real-time e lo si confronta con i dati sperimentali rilevati al banco. Il modello

messo a punto è di tipo non lineare: tramite la conoscenza dei valori degli

stati termodinamici del motore per vari casi di trim sperimentali, vengono

dapprima calcolati off-line dei punti di funzionamento che poi il modello

utilizzerà durante la simulazione attraverso una logica di adattamento non

lineare. A questo memorandum si farà riferimento in fase di validazione

perché esso contiene la gran parte dei dati sperimentali a cui si riferisce la

maggioranza dei lavori successivi: tali dati sperimentali saranno quelli con

cui ci si confronterà nelle simulazioni di validazione statica e dinamica della

parte termodinamica del modello. In [5] è anche presente un implementazione

del controllore del motore, cos`ı da poterne simulare il comportamento anche

in anello chiuso insieme ad un modello dinamico semplificato dell’elicottero

UH-60.

Prosecuzione del lavoro di Ballin ([5]) è il report [6], dove, a partire dal mo-

dello di Ballin, viene messo a punto un simulatore semplificato, linearizzato e

in anello aperto del solito motore GE t700. Collegando in una formulazione

agli stati diversi punti di funzionamento del modello trattato in [5], viene

quindi generato un sistema lineare che rende più agevole l’implementazione

e la sua risoluzione in real-time.

Nell’articolo di Camporeale e Fortunato [7] si trova un approccio più re-

2.2. SOFTWARE ESISTENTI

cente al problema: il lavoro si pone lo scopo di mettere a punto in ambiente

Simulink ®un codice modulare che consenta la simulazione di impianti in-

dustriali di turbine a gas. Il modello che viene messo a punto dagli autori

è del tipo “aerotermico”, similmente a quanto fatto successivamente anche

nella tesi di Kocer [8] e nel recente articolo di Kocer [9], con una rappresen-

tazione concentrata delle turbomacchine. La modellazione del componente

combustore da noi implementata è concorde con quella utilizzata in tale lavo-

ro, ovvero quella di un puro accumulatore di energia, delegando il bilancio di

massa ad un componente “plenum” che non trova riscontro in un componente

fisico.

In [9] e nella tesi [8] si ritrova un approccio moderno molto simile a quello

di Camporeale e al.

2.2

Software esistenti

I software disponibili per la simulazione di motori a turbina a gas si possono

raggruppare in due tipologie:

• software di simulazione del motore real-time;

• banchi prova virtuali non real-time (engine deck);

I software non real-time consentono una maggiore accuratezza della

modellazione e nella discretizzazione degli stati termodinamici del motore.

Questi software vengono spesso realizzati dal produttore del propulsore sulla

base di prove sperimentali ottenute al banco, e la loro utilità è quella di avere

a disposizione un modello molto accurato che può sostituire il vero motore

in molte applicazioni dove l’utilizzo del propulsore reale sarebbe impossibile

per questioni di sicurezza o di costo:

• prove con propulsore “in the loop”;

• prove pericolose che potrebbero danneggiare il propulsore;

• verifica dei sistemi di controllo e di diagnostica;

Fra i software esistenti in questo campo, uno dei più conosciuti è GSP:

Gas turbine Simulation Program ([10]), messo a punto dall’NLR (Nationaal

Lucht- en Ruimtevaartlaboratorium: centro di ricerca aeronautico olandese)

e scaricabile liberamente nella sua versione dimostrativa. Esso consente la

simulazione di una vasta gamma di motori/generatori a turbina a gas, es-

sendo un codice modulare con un’interfaccia grafica abbastanza intuitiva. Il

codice sorgente è scritto in Delphi, linguaggio orientato agli oggetti, e non

CAPITOLO 2. STATO DELL’ARTE

è reso disponibile nella versione liberamente scaricabile. GSP è un codice

completo che consente di simulare anche l’usura dei componenti e le perdite

di installazione del motore. Durante l’uso di questo programma abbiamo

registrato una difficoltà a convergere in situazioni lontane dal funzionamen-

to nominale, ed è quindi stato impossibile utilizzare tale programma come

confronto nelle simulazioni eseguite, se non in quelle in condizioni nominali,

dove erano disponibili anche dei dati sperimentali, che sono stati preferiti a

quelli generati da GSP. Va riconosciuto a questo programma una intuitiva

interfaccia grafica e una buona documentazione disponibile on line. Esso è

stato infatti impiegato in numerosi progetti europei di ricerca.

I software real-time sono utili invece per un numero di situazioni dove

è richiesto che il funzionamento avvenga nei tempi caratteristici del fenome-

no reale, rinunciando magari a parte dell’accuratezza garantita dagli engine

deck:

• osservatore del propulsore per i sistemi di bordo del velivolo;

• simulatore di motore da inserire in simulatori di volo;

Un software di questa tipologia è TERTS: Turbine Engine Real-Time

Simulator ([11]), messo a punto dallo stesso team dell’NLR di Amsterdam

che ha ideato anche GSP. TERTS è in effetti un software derivato da GSP,

che rispetta però un requisito real-time. Esso è implementato in ambiente

Matlab-Simulink.

Gran parte dei software real-time per la simulazione di motori implemen-

tano un modello lineare ottenuto dall’identificazione di un sistema. Molto

spesso vengono implementati modelli ”lineari a tratti” dove una serie di re-

lazioni lineari viene utilizzata per coprire lo spazio di stati altamente non

lineare [6]. Modelli lineari separati vengono determinati per differenti con-

dizioni operative (per esempio diverse velocità di rotazione). Questo tipo di

approccio è largamente diffuso per l’utilizzo in simulatori di volo o per il pro-

getto di sistemi di controllo. Per analisi più accurate di ogni nuovo effetto,

è comunque preferibile un approccio termodinamico che si attenga maggior-

mente alla fisica del problema, cos`ı da avere un modello di base che possa

essere impiegato in una maggiore varietà di prove. Il set di equazioni non

lineari che traduce generalmente i principi fisici di conservazione della massa,

bilancio di energia e conservazione del momento angolare, viene poi risolto in

diversi modi, fra i metodi di soluzione più utilizzati in letteratura ricordiamo:

il metodo di Newton-Raphson, la generazione di Jacobiani in prossimità di

alcuni punti di funzionamento noti, e il metodo di Broyden per l’aggiorna-

mento dello Jacobiano. Il calcolo dello Jacobiano, se da una parte consente

una rapida convergenza e un’efficace implementazione del codice real-time, in

2.2. SOFTWARE ESISTENTI

quanto può essere effettuato off-line prima della simulazione, dall’altra parte

può portare a un insieme di problemi quali dei punti di stazionarietà, mal

condizionamento numerico della matrice Jacobiana o problemi di oscillazio-

ne intorno a un punto di equilibrio. Il requisito real-time di avere un tempo

massimo di esecuzione del codice porta nella maggior parte di questi software

a un numero massimo di iterazioni concesse (troncamento delle iterazioni):

dopo un limitato numero di passi temporali (massimo tempo di esecuzione

all’interno di un passo temporale), l’iterazione viene fermata e l’accuratez-

za del calcolo viene accettata per quella che è. All’iterazione successiva si

assume che il modello sia giunto a convergenza Questa assunzione è ragio-

nevole se la simulazione del motore coinvolge transitori con grandi gradienti

che durano brevi periodi. Negli intervalli fra uno di questi transitori veloci

e l’altro si assume che il propulsore lavori quasi stazionariamente, cosicché

ogni errore residuo nelle equazioni venga eliminato.

Capitolo 3

Fisica del problema e

modellazione

La modellazione del problema è avvenuta tenendo sempre presente l’obiet-

tivo di un modello semplice, che riproduca i fenomeni fisici di interesse in

modo completo, senza però entrare troppo nel dettaglio del comportamento

dinamico sia dal punto di vista della gasdinamica che della chimica, cosa che

avrebbe richiesto un impegno eccessivo nella modellazione e nell’esecuzione

del codice, che nelle intenzioni del lavoro deve poter avvenire in real-time.

La determinazione del modello è stata pensata considerando come base di

partenza il motore General Electric t700 al quale si riferisce gran parte della

letteratura in materia (es.

[5, 6]), cos`ı da potersi confrontare su casi test il

più possibile vicini a quelli disponibili. Si è posta attenzione all’estendibilità

del modello a propulsori diversi enfatizzando il requisito di portabilità del

modello, rendendolo dipendente da alcuni parametri fisici caratterizzanti il

motore possibilmente reperibili in letteratura per nuovi e differenti motori.

3.1

Modellazione delle turbomacchine

In accordo con la letteratura più recente [7, 8, 9], si è scelto di modellare

le turbomacchine (compressore e turbine) come dei componenti istantanei,

considerando la loro dinamica interna talmente veloce da non influire sulle

dinamiche più lente del sistema complessivo. Le trasformazioni che avven-

gono all’interno delle turbomacchine sono state modellate considerando la

corrispettiva espansione o compressione isoentropica e correggendola poi con

un rendimento per determinare il punto finale della trasformazione reale. La

determinazione del rapporto di compressione o di espansione per ogni tur-

bomacchina è stato valutato a partire da mappe che parametrizzano queste

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

grandezze in funzione di un parametro caratteristico che identifica il regime

di funzionamento del motore nella condizione considerata. Tale parametro,

per quanto riguarda il compressore e la turbina generatrice di gas, è stato

identificato nel regime di rotazione dell’albero della turbina generatrice di

gas (n₁), e invece nella temperatura di uscita dalla turbina generatrice di gas

(T₄) per la turbina di potenza.

La determinazione delle portata richiamata ed elaborata dal compressore

viene fatta in modo analogo al calcolo del rapporto di compressione, para-

metrizzando la portata in funzione del regime di rotazione dell’albero della

turbina generatrice di gas (n₁).

Le scelte delle grandezze sulle quali parametrizzare le mappe sono state

dettate dal ragionamento sulla fisica del problema: la scelta del parametro

n₁ per la determinazione dei rapporti di compressione e della portata ela-

borata dalle turbomacchine calettate sull’albero della turbina generatrice di

gas (compressore e turbina di alta pressione) sono in accordo con il fatto

che una rotazione a maggiore velocità del compressore comporta un richiamo

maggiore di aria, e quindi una maggiore portata, e che la stessa rotazione

a maggiore velocità comporta un maggiore lavoro di compressione e quindi

un salto maggiore della pressione fra l’imbocco e l’uscita del componente.

Lo stesso ragionamento si può trasferire alla turbina. Si pensi ad esempio

a un compressore che lavori al doppio della velocità di un componente ana-

logo in due identici propulsori: il componente alla velocità doppia, a meno

di regolazioni dell’ingresso dell’aria de parte dell’IGV (Inlet Guide Vane) e

di limiti dovuti al regime di moto del fluido nella turbomacchine, elaborerà

una portata maggiore del suo omologo che gira a più bassa velocità e farà

eseguire al fluido una maggiore variazione di pressione fra l’ingresso e l’uscita

del componente.

Per quanto riguarda la scelta del parametro su cui valutare il regime di

funzionamento della turbina di potenza, la scelta è caduta sulla temperatura

di uscita dalla prima turbina, in quanto questo è uno dei parametri variabili

con i vari regimi del motore ed era disponibile un confronto con dati speri-

mentali in letteratura ([5]). Sarebbe stato possibile parametrizzare le mappe

su altre grandezze (ad esempio la pressione p₄ a valle della turbina HPT), ma

non si avrebbe avuto modo di confrontarsi con quello che è il funzionamen-

to di un motore vero al banco. Infatti in letteratura non è disponibile una

campagna di prove che fornisca maggiori variabili di interesse per un motore

elicotteristico in prove al banco.

3.2. PROPRIETÀ DEI FLUIDI

3.2

Proprietà dei fluidi

Determinazione delle proprietà dell’aria. A causa delle grandi varia-

zioni di temperatura e per mantenere un’elevata accuratezza della simula-

zione, le proprietà dei fluidi di lavoro, aria e gas combusti, devono essere

calcolate per ogni condizione in ogni stadio del propulsore. L’aria è consi-

derata una miscela di gas semi-ideali, similmente a quanto fatto in [7], e le

proprietà termodinamiche della miscela (calori specifici c_p, c_v e loro rapporto

k ) vengono determinate in base alla sola temperatura tramite una tabella.

La pressione della miscela aria è la somma delle pressioni parziali dei gas che

la compongono. Sotto queste ipotesi le proprietà termodinamiche del fluido

sono funzioni della temperatura e della composizione della miscela. Il fluido

aria è quello elaborato all’interno dei componenti presa d’aria, compressore

e plenum.

Determinazione delle proprietà dei gas combusti. Sotto l’ipotesi, ge-

neralmente verificata nei moderni motori, che la combustione avvenga con

grande eccesso d’aria e che il processo di combustione possa essere consi-

derato completo all’uscita della camera di combustione, il fluido di lavoro

considerato durante tutte le trasformazioni all’interno delle turbine e della

camera di combustione è una miscela di gas combusti, con uno specifico rap-

porto aria/combustibile (FAR: Fuel to Air Ratio), noto dall’apporto di com-

bustibile e di aria nel motore all’istante considerato. Per la definizione delle

proprietà termodinamiche dei fluidi di lavoro si è fatto riferimento al capitolo

3 del testo [1]. La costante dei gas R (dimensionalmente [J/KgK]) per la mi-

scela di gas combusti è calcolata in funzione del rapporto aria/combustibile

(FAR) come:

R_kerosene = 287.05 - (0.00990) (FAR) + (10^-7) (FAR)²

dove FAR è definito come:

m_fuel

F AR =

m_air

Si noti ( figura 3.1) che per i gas ottenuti dalla combustione di kerosene

con aria, la costante dei gas R non cambia in modo considerevole per valori

di FAR quali quelli utilizzati nei motori aeronautici (0.01:0.05).

Il calore specifico c_p dei gas combusti provenienti dalla combustione del

kerosene è funzione della temperatura e della composizione della miscela

aria/combustibile. Il kerosene impiegato nella propulsione aeronautica ha

una composizione chimica abbastanza omogenea, e il suo calore specifico per

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

KEROSENE COMBUSTION PRODUCTS GAS CONSTANT

295

294

293

292

291

290

289

288

287

286

285

0.01

0.02

0.03

0.04

0.05

FAR - Fuel to Air Ratio

Figura 3.1: Variazione della costante dei gas R in funzione della temperatura

per i gas combusti di kerosene.

temperature comprese fra 200 K e 2000 K è ben approssimatao dalla equa-

zione 3.1, formula empirica che tiene conto della temperatura e del rapporto

aria/combustibile:

c_p = A₀ + A₁T_Z + A₂T^2Z + A₃T^3Z + A₄T^4Z + A₅T^5Z + A₆T^6Z+

F AR

A₇T^7Z + A₈T^8Z +

(B₀ + B₁T_Z + B₂T^2Z + B₃T^3Z+

(3.1)

1 + FAR

B₄T^4Z + B₅T^5Z + B₆T^6Z + B₇T^7Z)

dove T_Z =^T

[K] e A₀ . . . A₈ e B₀ . . . B₇ sono i coefficienti del polinomio

1000

riportati in tabella 3.2.

Per l’andamento del calore specifico dei gas combusti del kerosene si faccia

riferimento alla figura 3.2

Per il calcolo del rapporto fra calori specifici dei gas combusti del kerosene,

denominato con k si è fatto riferimento alla relazione:

c_p

c_v

c_p - R

L’andamento del rapporto dei calori specifici è riportato nella figura 3.3.

3.2. PROPRIETÀ DEI FLUIDI

A₀

0.992313

B₀

-0.718874

A₁

0.236688

B₁

8.74781

A₂

-1.852148

B₂

-15.863157

A₃

6.083152

B₃

17254096

A₄

-8.893933

B₄

-10.233795

A₅

7.097112

B₅

3.081778

A₆

-3.234725

B₆

-0.361112

A₇

0.794571

B₇

-0.003919

A₈

-0.081873

B₈

0.0555930

Tabella 3.1: Coefficienti dei polinomi per il calcolo di k e c_p dei gas combusti

nel modello di Walsh [1].

KEROSENE COMBUSTION PRODUCTS SPECIFIC HEAT - Cp

1.4

1.35

FAR=0

FAR=0.01

FAR=0.02

1.3

FAR=0.03

FAR=0.04

1.25

FAR=0.05

1.2

1.15

1.1

1.05

200

400

600

800

1000

1200

1400

1600

1800

2000

Temperature [K]

Figura 3.2: Calore specifico dei gas ottenuti dalla combustione di kerosene

in funzione della temperatura per diversi rapporti aria/combustibile.

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

KEROSENE COMBUSTION PRODUCTS SPECIFIC HEAT RATIO - k

1.42

FAR=0

FAR=0.01

1.4

FAR=0.02

FAR=0.03

1.38

FAR=0.04

FAR=0.05

1.36

1.34

1.32

1.3

1.28

1.26

200

400

600

800

1000

1200

1400

1600

1800

2000

Temperature [K]

Figura 3.3: Rapporto fra calori specifici del kerosene in funzione di

temperatura e rapporto aria/combustibile.

Le relazioni presentate sono valide per temperature comprese fra 200 e

2000 K, ma dal momento che in camera di combustione si possono raggiun-

gere per brevi momenti anche temperature dell’ordine dei 2500 K, si rende

necessario un modello più complesso, quale quello messo a punto in [12]: in

esso vengono implementate le seguenti relazioni, valide fra 300 e 3500 K e

anche oltre con errori modesti:

c_{pcombustibile}

1 + FAR(I0 + I1 T + ^I2)

T²

c_paria

1 + FAR

dove i coefficienti I0, I1 e I2 si determinano come:

18.0566 n_C + 8.3485 n_H + 15.1616 n_O

I0=

(3.2)

M_combustibile

0.00223 n_H

I1=

(3.3)

M_combustibile

1077768.4 n_C

I2=

(3.4)

M_combustibile

in cui M_combustibile è la massa molecolare di una molecola di combustibile

(nel nostro caso kerosene, M_kerosene = 170), mentre n_C , n_O, n_H sono il numero

3.3. COMPONENTI DEL MODELLO

di moli di carbonio, ossigeno e idrogeno necessarie alla generica reazione di

combustione:

n_H

n_O

n_H

C_n

H_n

O_n

+ (n_C +

)O₂ = n_cCO₂ + (

)H₂O

Nel caso di combustione di kerosene in aria secca in rapporto stechiometrico

avremo:

C₁₂ H₂₆ +

O₂ = 12 CO₂ + 13 H₂O

e quindi i numeri di moli assumono i valori:

n_C = 12

(3.5)

n_H = 26

(3.6)

n_O = 0

(3.7)

L’utilizzo di questo schema di soluzione permette un’ottima rappresen-

tazione delle proprietà dei gas combusti di kerosene in un intervallo di tem-

perature molto vasto (200-5000 K), con errori quantificati da [12] in qualche

decimo di punto percentuale.

Nelle figure 3.5 e 3.4 si è confrontato il comportamento dei gas di scarico

del modello di Guha e quello dell’aria per temperature fino ai 3000 K. Il

modello polinomiale di [1], invece, per temperature superiori ai 2200 K non

è più utilizzabile perchè fornisce valori di c_p e k palesemente impossibili: si

veda a riguardo la figura 3.6.

3.3

Componenti del modello

In figura 3.7 sono visibili tutti i componenti del modello di motore. Si noti

come non ci sia un’effettiva corrispondenza fra i componenti fisici e quelli del

modello numerico.

Il modello numerico è composto dei seguenti componenti:

• presa d’aria;

• compressore;

• plenum;

• combustore;

• turbina di alta pressione;

• turbina di potenza;

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

RAPPORTO FRA CALORI SPECIFICI k - FAR = 1/35

1.55

ARIA

1.5

GAScombusti

1.45

1.4

1.35

1.3

1.25

500

1000

1500

2000

2500

3000

Temp [K]

Figura 3.4: Rapporto fra calori specifici k dell’aria e dei gas combusti di

kerosene con rapporto aria/combustibile di

ottenuto con il modello di

Guha.

CALORE SPECIFICO A PRESSIONE COSTANTE cp - FAR = 1/35

1700

ARIA

1600

GAScombusti

1500

1400

1300

1200

1100

1000

900

500

1000

1500

2000

2500

3000

Temp [K]

Figura 3.5: Calore specifico c_p dell’aria e dei gas combusti con rappor-

to aria/combustibile di

, ottenuto con il modello di Guha per un ampio

intervallo di temperature.

3.3. COMPONENTI DEL MODELLO

CONFRONTO cp ARIA E GAS COMBUSTI CON MODELLO DI WALSH

2000

1000

KEROSENE

ARIA

-1000

-2000

-3000

-4000

-5000

500

1000

1500

2000

2500

3000

Temp [K]

Figura 3.6: Calore specifico di aria e gas combusti ottenuto con modello di

Walsh, con FAR=1/35.

• carico.

Si noti come i componenti del modello non corrispondano esattamente a

quelli di un motore fisico. E’ evidente che il plenum non è presente in un

motore fisico; questo componente, come verrà spiegato più estesamente nel

seguito, serve a implementare un principio di conservazione della massa non

stazionario all’interno di un volume. Allo stesso modo nel modello sono pre-

senti un generico compressore e due generiche turbine non suddivise in stadi,

e non è specificato se questi componenti siano di tipo assiale o centrifugo.

Questa generalità dei componenti del modello è tesa a riprodurre il mag-

gior numero di motori possibile, senza entrare nel dettaglio della specifica

configurazione. Nell’implementazione del modello si è scelto di riprodurre i

componenti agli effetti esterni e non entrare nel dettaglio delle loro peculia-

rità proprio per mantenere maggiore generalità e portabilità rispetto a diversi

sistemi.

3.3.1

Presa d’aria

La presa d’aria è del tipo dinamico, ideale perché si ipotizza che non vi

siano perdite per attrito. In presenza di una componente di velocità v non

nulla in direzione perpendicolare al piano dell’ingresso dell’aria nel motore, la

3.3. COMPONENTI DEL MODELLO

pressione totale e la temperatura totale all’ingresso del compressore vengono

incrementate di una quantità fornita dalle relazioni:

Δp =

ρv²

(3.8)

ΔT =

(3.9)

2c_p

La 3.8 traduce il teorema di Bernoulli per flussi incomprimibili, la 3.9 deriva

dal primo principio della termodinamica. in questo modo la pressione totale

si aggiorna come:

p_tot = p_s + Δp

(3.10)

e la temperatura totale:

T_tot = T_s + ΔT

(3.11)

dove p_s e T_s sono la pressione statica e la temperatura statica, assunte

omogenee su tutta la sezione all’imbocco della presa d’aria.

3.3.2

Compressore

Il compressore è un componente adiabatico, istantaneo e privo di dimensione

(0-D). Durante i transitori, in accordo con [7], la risposta del compressore

è assunta essere quella di un componente quasi stazionario, ovvero la sua

dinamica è considerata molto più veloce di quella dei componenti più lenti,

e sono quindi utilizzate mappe ricavate dal funzionamento del componente

a regime. Le mappe forniscono la portata elaborata dal compressore g₂₁ e

il rapporto di compressione β_C in funzione della velocità di rotazione n₁,

passando attraverso una normalizzazione e una successiva scalatura come

verrà illustrato in dettaglio nel paragrafo relativo alle mappe.

β_c = f(n₁)

(3.12)

g₂₁ = f(n₁)

(3.13)

La temperatura al termine della compressione viene calcolata attraverso una

trasformazione isoentropica corretta con un rendimento. Chiamando T₁ e p₁

la pressione e la temperatura all’imbocco del compressore e quindi all’uscita

della presa d’aria avremo:

p₂₁

T_21is = T₁(

)^kk1

(3.14)

p₁

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

La 3.14 è la trasformazione che l’aria subirebbe passando dalla pressione P₁

alla pressione P₂₁ idealmente senza perdite e senza aumento di entropia (tra-

sformazione isoentropica). Naturalmente, noto il rapporto di compressione

β_c definito come

p₂₁

β_c =

(3.15)

p₁

è possibile riscrivere la 3.14 come:

T_21is = T₁(β_c)^kk1

(3.16)

Nota la T_21is, che è la temperatura di uscita dalla compressione isoentropica,

è possibile calcolare il lavoro compiuto dal compressore correggendo il lavoro

isoentropico con un rendimento:

L_compressore = c_p(T₁ - T_21is)η_c

(3.17)

Il rendimento η_c, definito rendimento isoentropico del compressore, è assunto

essere costante in tutto il campo di funzionamento della turbomacchina, ma

potrebbe in futuro essere parametrizzato in funzione di una grandezza quale

il regime di rotazione dell’albero. E’ infine possibile calcolare la temperatura

reale al termine della compressione attraverso la 3.18:

L_compressore

T₂₁ = T₁ -

(3.18)

c_p

La temperatura T₂₁ all’uscita dal compressore è quindi dipendente dalle con-

dizioni ambientali di temperatura e pressione esterne, dalla velocità di volo

del velivolo e dal rapporto di compressione che il compressore esprime al

momento.

3.3.3

Plenum

Come accennato in precedenza, utilizzando un plenum si riesce a riprodurre

un effetto di comprimibilità dei gas al passaggio nella camera di combustio-

ne e quindi anche un ritardo nell’adattamento della pressione. L’utilizzo del

plenum, comune a molta letteratura [7, 9], permette di attribuire alla camera

di combustione o combustore la sola funzione di accumulatore di energia. Il

plenum si rende necessario per tenere conto del volume dei vari componenti

che nell’implementazione del modello sono trattati come elementi 0-D ovve-

ro privi di volume. Nella nostra implementazione l’unico plenum impiegato

è posto fra l’uscita del compressore e l’imbocco del combustore, e serve a

considerare il bilancio di massa dovuto all’aria proveniente dal compressore,

3.3. COMPONENTI DEL MODELLO

l’aria spillata per le varie utenze, l’aria diretta al raffreddamento della tur-

bina, l’aria elaborata dal combustore. Il bilancio di massa non stazionario è

dato dall’equazione:

V_plenum dp_out

=g_in -g_out

(3.19)

m R T_out dt

dove V_plenum è il volume attribuito al componente plenum, T_out e p_out la

temperatura e la pressione all’uscita dal plenum, g_in e g_out le portate massiche

in ingresso e in uscita dal componente plenum. La portata d’aria in ingresso

al plenum è quella elaborata dal compressore:

g_in = g_compressore = g₂₁

La portata g_out che abbandona il plenum è divisa in due parti:

g_out = g₂ = g_cool + g₃

• g_cool: portata d’aria spillata e inviata alle utenze e al raffreddamento

della turbina HPT.

• g₃: portata d’aria che viene inviata nel combustore per partecipare alla

combustione.

Possiamo quindi determinare la portata massica in uscita dal plenum g₂

in funzione della variazione della pressione in uscita dal plenum∂p2

e delle

∂t

portate massiche di aria del compressore e spillate:

V_plenum

dp₂

g₂ = g₂₁ - g_cool -

τ_plenum

k₂RT₂

dove la costante di tempo τ_plenum è definita come:

M_plenum

τ_plenum =

k₂g₂

Poichè la pressione p₂ è, per l’ipotesi di nessuna perdita nel plenum, la

stessa pressione all’uscita dal compressore, determinata tramite una mappa

in funzione della velocità di rotazione n₁ e della temperatura p₁, ovvero:

p₂ = p₂₁ = β_C(n₁)p₁ = f(n₁,p₁)

Allora è possibile esprimere la derivata∂p2

in funzione dell’accelerazione rota-

∂t

tiva del primo albero n˙₁, che è uno dei residui del sistema, secondo il seguente

ragionamento:

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

∂p₂

=β˙_Cp₁ +β_Cp˙₁

∂t

Poichè si può assumere che p₁ vari molto lentamente rispetto a β_C , ovvero:

1 ≈0

allora possiamo scrivere:

∂p₂

=β˙_Cp₁

∂t

Dove la derivata temporale

C è in realtà funzione dell’accelerazione di

rotazione, e si può cos`ı esplicitare:

∂β_C

∂β_C ∂n₁

∂β_C

C =

∂t

∂n1 ∂t

∂n₁

Per valutare la derivata del rapporto di compressione rispetto alla velocità

di rotazione,∂βC si può alternativamente:

∂n₁

• valutare la derivata del polinomio con cui si è costituita la mappa di

β_C.

• approssimare tale derivata nel campo di funzionamento.

Nel caso∂βC

si possa valutare dal polinomio che esprime la mappa di β_C

∂n₁

β_C = a₀ + a₁n₁ + a₂n²¹ + a₃n³¹ + ··· + a_nnⁿ¹

allora

∂β_C

=a₁ +a₂n₁ +a₃n²¹ +···+a_nn^n-11

∂n₁

Nel caso invece si voglia approssimare la derivata, si può ad esempio

ipotizzare un andamento lineare del rapporto di compressione con la velocità

di rotazione dell’albero e quindi:

β_Cnom(n₁ = n_1nom) - β_C(n₁ = 0)

β_C = β_C(n₁ = 0) +

n₁

n_1nom

Poichè in condizioni di compressore fermo si può assumere che la pressione a

valle del componente sia la stessa a valle dello stesso (nessuna compressione):

β_c(n₁ = 0) = 1

allora la derivata parziale assume il valore costante:

∂β_C

β_Cnom - 1

∂n₁

n_1nom

3.3. COMPONENTI DEL MODELLO

Otteniamo quindi la portata in uscita dal plenum ad ogni istante come:

p₂V_plenum ∂β_C

g₂ = m_out = g₂₁ - g_cool -

n₁

(3.20)

k₂RT₂

∂n₁

3.3.4

Combustore

La camera di combustione è, nel nostro modello, un puro accumulatore di

energia. Essa viene implementata con il nome di combustore (burner), per

risaltare la funzione di generazione di energia attraverso la reazione di com-

bustione. Alla camera di combustione è associata la sola dinamica termica

del sistema: quella della temperatura in camera di combustione T₃. Essendo

la camera di combustione l’elemento subito a valle del plenum, si assume

che in essa le condizioni all’imbocco e allo sbocco siano omogenee sulla se-

zione ed eguali a quelle in uscita dal plenum, che sono le stesse in uscita dal

compressore:

p_incombustore = p₂ = p₂₁

T_incombustore = T₂ = T₂₁

All’interno del componente combustore avviene la miscelazione dell’aria

con il combustibile. La portata di combustibile g_b è l’ingresso principale del

sistema, mentre la portata d’aria in ingresso è quella proveniente dal plenum,

come già visto nel paragrafo relativo:

g₃ = g₂ - g_cool

Il componente combustore non è privo di perdite come il plenum, ma fra

l’imbocco e l’uscita è prevista una perdita di carico per il fluido quantificabile

con un rendimento y_CC :

p₃ = p_{outcombustore} = p₂ y_CC

Il rendimento y_CC è assunto essere costante.

L’equazione che descrive la dinamica del combustore è ottenuta dal bi-

lancio di energia instazionario:

d(M_CC u_CC )

g_inh_in + g_b(h_b + η_bLHV ) - g_outh_out

(3.21)

g_outc_pout

dove M_CC e u_CC sono rispettivamente la massa di gas presenti all’interno

del componente combustore e l’energia interna specifica dei gas; LHV è il

potere calorifico inferiore del combustibile (Lower Heating Value).

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

Trascurando le variazioni della massa e del calore specifico l’equazione

3.21 può essere riscritta come

g_inh_in + g_b(h_b + η_bLHV ) - g_outh_out

out =

(3.22)

g_outc_poutτ_CC

La costante di tempo τ_CC può essere valutata come:

M_CC

τ_CC =

k g_out

Usando la nomenclatura del modello e della figura 3.7, la dinamica ter-

mica, relativa alla temperatura in camera di combustione può essere espressa

dalla:

g₂ h₂ + g_b(h_b + η_b LHV ) - g₃ h₃

3 =

(3.23)

g₃ c_p3 τ_cc

dove le entalpie specifiche h₃, h₂, h_b sono calcolate come:

h₂ = c_p2 T₂

h₃ = c_p3 T₃

h_b = c_vkerosene T_kerosene

Per le proprietà del combustibile si è fatto riferimento a [13].

Riscriviamo l’equazione 3.23 può essere riscritta esplicitando le relazioni

delle entalpie specifiche:

R T₃ k₃

3 =

[g₂ cp₂ T₂ +g_b(c_vb h_b +η_b LHV )-(g₂ +g_b)c_p3 T₃] (3.24)

p₂ y_CC V_CC c_p3

Nell’equazione 3.24, le grandezze dipendenti da T₃, e quindi dalla stessa

derivata

3 sono:

k₃,c_p3,T₃

e quindi l’equazione risulta implicita.

I calori specifici e i loro rapporti sono calcolati attraverso il modello di

Guha ([12]) per i gas combusti di kerosene, come spiegato nella sezione rela-

tiva ai fluidi e alle loro proprietà. I metodi di integrazione adottati dovranno

tenere in conto l’implicitezza del sistema per la risoluzione delle equazioni

differenziali.

3.3. COMPONENTI DEL MODELLO

3.3.5

Turbina generatrice di gas

La turbina generatrice di gas (HPT: High Pressure Turbine) è modellata

come un componente istantaneo e privo di volume, a singolo stadio. L’ipo-

tesi che sta alla base di una tale semplificazione, similmente a quanto fatto

per il compressore, è che la dinamica del componente turbina sia molto più

veloce delle due dinamiche meccaniche e della dinamica termica del siste-

ma, cosicchè il componente possa essere approssimato in modo quasi-statico.

Similmente al compressore, il rapporto di espansione della turbina β_HPT è

parametrizzato, attraverso una mappa, su un parametro caratteristico che in-

dividui il regime di funzionamento del componente: la velocità di rotazione

n₁, normalizzata rispetto al suo valore nominale.

β_HPT = f(n₁)

(3.25)

Il rapporto di espansione cos`ı ottenuto è un β globale, che tiene conto di tutti

gli eventuali stadi di una turbina reale.

La pressione a valle della turbina di alta pressione è quindi determinata

a partire dal regime di rotazione normalizzato e dalla pressione p₃ in camera

di combustione:

p₄ = p_outHPT = β_HPT p₃ = f(p₃,n₁)

(3.26)

Alla turbina HPT giunge parte dell’aria spillata dal compressore, per il

raffreddamento delle palette, di quest’afflusso di aria si tiene conto con il

bilancio di portata:

g₄ = g_outHPT = g₃ + g_cool

(3.27)

Per il calcolo della temperatura in uscita dalla turbina di alta pressione,

necessaria poi al calcolo della potenza generata, si è fatto ricorso a una espan-

sione isoentropica corretta con un rendimento: si trova prima la temperatura

relativa alla trasformazione isoentropica che porti il fluido alla stessa condi-

zione di pressione prevista, si calcola il lavoro che produce tale espansione

e infine questo valore di lavoro viene corretto con un rendimento isoentropi-

co. Il valore del lavoro calcolato con il rendimento isoentropico viene quindi

diviso per il calore specifico, determinato con la temperatura della trasfor-

mazione isoentropica corrispondente, per ottenere la temperatura finale al

termine dell’espansione reale.

T_4is = T₃(

p4 )kk

=T3(βHPT)kk

(3.28)

p₃

c_p4is = f(T_4is,FAR)

(3.29)

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

L_HPT = cp_4is(T₃ - T_4is)η_HPT

(3.30)

c_p4 = f(T₄)

(3.31)

L_HPT

T₄ = T₃ -

(3.32)

c_p4

La potenza generata dalla turbina è immediatamente calcolabile, noto il

salto entalpico fra l’imbocco e lo sbocco del componente:

P_HPT = g₃ h₃ + g_cool h_cool - g₄ h₄

(3.33)

dove le entalpie sono calcolate nel modo seguente:

h₃ = c_p3 T₃

(3.34)

h_cool = c_p2 T₂ = c_p21 T₂₁

(3.35)

h₄ = c_p4 T₄

(3.36)

L’entalpia viene calcolata utilizzando il modello di Guha per le proprietà

dei gas combusti di kerosene. Si noti che la portata g_cool è composta di

sola aria alla temperatura di uscita dal plenum, in quanto è la portata d’aria

utilizzata per raffreddare le palette della turbina. Questo comporta un nuovo

calcolo del rapporto aria/combustibile FAR, che modifica quindi le proprietà

dei gas di scarico.

3.3.6

Turbina di potenza

La turbina di potenza (PT: power turbine) è stata modellata similmente a

alla turbina generatrice di gas eccetto che per due differenze:

1. la turbina PT a differenza della turbina HPT non riceve una portata

d’aria più fresca per il raffreddamento delle palette;

2. la pressione a valle della turbina PT non può assumere valori inferiori

a quelli della pressione esterna (condizione di flusso adattato).

La portata di gas elaborata dalla turbina PT, non essendoci ulteriori

immissioni di fluido, è la stessa che proviene dall’uscita della turbina HPT:

g₅ = g₄

(3.37)

3.3. COMPONENTI DEL MODELLO

La determinazione della pressione all’uscita della turbina avviene sempre

attraverso una mappa, ma poichè in questo caso il componente deve lavo-

rare in un regime di giri pressochè costante in quanto dalla sua velocità di

rotazione dipende la velocità di rotazione del carico (rotore), non è possibile

parametrizzare il β di espansione sul regime di rotazione n₂. Il parametro

scelto per la parametrizzazione del rapporto di espansione β_PT è la tempera-

tura in uscita dalla turbina generatrice di gas, normalizzata con il suo valore

nominale.

β_PT = f(T₄)

(3.38)

p₅ = β_PT p₄ = f(p₄,T₄)

(3.39)

in cui la pressione p₅ non può mai essere inferiore a quella atmosferica.

Noto il valore di p₅ è possibile determinare il valore di temperatura dei gas

all’uscita della turbina tramite un’espansione isoentropica corretta con un

rendimento, esattamente come fatto nel componente HPT.

T_5is = T₄(

p5 )kk

(3.40)

p₄

cp_5is = f(T_5is, F AR5)

(3.41)

L_PT = c_p5is(T₄ - T_5is) η_PT

(3.42)

c_p5 = f(T₅,FAR)

(3.43)

L_PT

T₅ = T₄ -

(3.44)

c_p5

La potenza, espressa in Watt, della turbina di potenza è proporzionale al

salto entalpico a cavallo del componente ed alla portata d’aria elaborata:

P_PT = g₅ (c_p4 T₄ - c_p5 T₅)

(3.45)

3.3.7

Modello di carico semplificato

Dovendo provare il funzionamento del codice senza avere a disposizione fin

da subito un modello sofisticato di rotore da connettere al propulsore, è stato

messo a punto un modello semplificato di carico, che riproducesse in maniera

approssimata il rotore dell’elicottero Sikorsky UH-60.

I dati del rotore considerato sono:

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

Figura 3.8: Pala di rotore semplificata.

• diametro rotore: D = 16.36 m;

• numero pale: 4;

• superficie del disco rotore: S = 210 m²;

• corda media della pala del rotore principale : c = 0.53 m;

• profilo utilizzato per le pale del rotore principale: SC1095;

• velocità angolare nominale del rotore ω_ROT = 240 RP M;

Sotto le ipotesi di:

• volo in hovering (nessuna componente di velocità longitudinale);

• coppia resistente dovuta solo all alla resistenza delle pale del rotore;

• pala rettangolare non svergolata e rigida;

• movimenti di flappeggio e di comando trascurati;

• rapporto di trasmissione costante in ogni condizione di funzionamento;

• presenza di un solo tipo di profilo in tutta l’apertura;

• coefficiente di resistenza C_D costante e pari a quello di massima effi-

cienza del profilo;

• densità dell’aria costante e pari al valore standard ISA;

si è calcolata la resistenza per unità di lunghezza del profilo:

ρv² cC_D

(3.46)

Integrando la resistenza per unità di lunghezza sul Raggio del rotore R =

otteniamo la resistenza di una pala:

3.3. COMPONENTI DEL MODELLO

∫ R

D_pala =

ρ v² c C_D dy

(3.47)

Si noti che la velocità della sezione infinitesima del rotore dipende in

generale dalla coordinata y e dalla velocità di volo dell’elicottero.

v = f(y,v_volo)

. Ad esempio per la pala che avanza parallela alla velocità di volo:

v(y, v_volo) = v_volo + ω_ROT y

(3.48)

ma considerando come da ipotesi l’elicottero in hover (v_volo = 0):

v(y) = ω_ROT y

(3.49)

Moltiplicando la resistenza di un tratto di pala infinitesimo per la sua

distanza dal centro (coordinata y in figura 3.8), integrando il contributo infi-

nitesimo sul raggio del rotore e moltiplicando per il numero di pale, otteniamo

la coppia esercitata dal rotore :

∫

∫ R

C =n_pale

D_palady = 4

ρ v² c C_D y dy

(3.50)

sostituendo alla velocità v(y) la sua espressione 3.49, otteniamo:

∫ R

C =4

ρ (ω_ROT y)² c C_D y dy

(3.51)

dove, per le ipotesi prima descritte, è possibile estrarre dall’integrale tutti i

termini non dipendenti da y e ottenere:

C =2ρcC_D ω^2ROT [

]^R0 =

ρcC_D R⁴ω^2ROT

(3.52)

Sappiamo che la velocità angolare del rotore è legata alla velocità dell’albero

di potenza n₂ attaverso il rapporto di riduzione della trasmissione:

ω_ROT = n₂ R_trasmissione

(3.53)

Essendo i regimi di rotazione nominali noti n2nom e ω_ROT , allora è noto anche

il rapporto di trasmissione:

n₂

nominale

R_trasmissione =

(3.54)

ω_ROT

nominale

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

x 10⁶

CARICO

2.5

NOMINALE

punto di funzionamento nominale

1.5

0.5

1.5

2.5

Angular Speed [RPM]

x 10⁴

Figura 3.9: Potenza richiesta dal carico modellato in funzione della velocita

di rotazione n₂ confrontata con i valori di potenza e velocità di rotazione

nominali per il motore GE t700.

Sostituendo la 3.53 nella 3.52 è quindi possibile scrivere il valore di coppia

richiesta dal carico:

C_LOAD(n₂) =

ρ c C_D R⁴ (n₂ R_trasmissione)² = K_R n²²

(3.55)

Nella 3.55 è possibile notare che la coppia richiesta dal carico è una funzione

quadratica della velocità di rotazione n₂. La costante K_R è pari a:

K_R =

ρcC_D R⁴ R^{2trasmissione}

(3.56)

La potenza dissipata dal carico è quindi ottenibile moltiplicando la coppia

C_LOAD per la velocità di rotazione dell’albero:

P_LOAD = C_LOAD n₂ = K_R n³²

(3.57)

Si riporta in figura 3.9 l’andamento della potenza richiesta dal carico cosı

modellato in funzione della velocità di rotazione n₂.

3.3.8

Dinamica degli alberi di trasmissione

Note le potenze espresse dalle turbine P_PT e P_HPT , la potenza richiesta dal

compressore P_C , la potenza richiesta dal carico P_LOAD e le inerzie dei due

3.4. DINAMICA DEL SISTEMA E METODO DI SOLUZIONE

alberi di rotazione J_PT e J_HPT (intese come inerzie complessive dell’albero e

degli elementi a lui connessi: compressori, turbine, carichi, ecc..), è possibile

scrivere le dinamica meccaniche del motore, che traduce il principio di con-

servazione del momento angolare. Nella loro modellazione si è tenuto conto

di alcune perdite per attrito, modellate come frazioni delle potenze delle tur-

bine: P_f

eP_f

sono le potenze perdute rispettivamente sull’albero della

HPT

turbina generatrice di gas e sull’albero della turbina di potenza:

P_f

=K_f1 P_HPT

(3.58)

HPT

P_f

=K_f2 P_PT

(3.59)

dove K_f1 e K_f2 sono moltiplicatori che esprimono la frazione di potenza persa

in attriti (attriti meccanici, resistenza sul disco dei rotori, e potenza perduta

per muovere eventuali utenze ausiliarie). Nel modello si è implmentato un

valore di qualche centesimo, in linea con [7].

L’equilibrio dei momenti è esprimibile come:

dω

(3.60)

Jω

dove J è il momento di inerzia del corpo rotante, P la potenza delle forze

agenti sul sistema, ω la velocità angolare.

Le espressioni delle dinamiche meccaniche sui due alberi si possono scri-

vere esplicitando le accelerazioni angolari come:

P_HPT - P_C - P_f

HPT

1 =

(3.61)

J_HPT n₁

P_PT - P_load - P_f

2 =

(3.62)

J_PT n₂

3.4

Dinamica del sistema e metodo di solu-

zione

La dinamica del sistema è governata da equazioni algebriche ed equazioni

differenziali. Le grandezze che sono determinate da equazioni differenziali,

tutte del primo ordine, sono tre:

• la dinamica meccanica determina le accelerazioni angolari dei due alberi

1 en˙2;

• la dinamica termica determina la variazione temporale di temperatura

in camera di combustione

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

Il sistema di 3 equazioni differenziali 3.63 accoppiate e equazioni algebri-

che fornisce quindi i tre residui

3, n1 e n2:



P_HPT - P_C - P_f



HPT



1 =



J_HPT n₁



P_PT - P_load - P

f_PT

2 =



J_PT n₂



R T₃ k₃

T˙₃ =

[g₂ c_p2 T₂ + g_b(c_vb h_b + η_b LHV ) - (g₂ + g_b)c_p3 T₃]

p₂ y_CC V_CC cp₃

(3.63)

Il sistema di 3 equazioni differenziali 3.63 accoppiate e equazioni algebriche

fornisce quindi i tre residui

3, n1 e n2. L’integrazione dei tre residui fornisce

il vettore di stato del sistema.

I residui sono ordinabili in un vettore y^′:

)

y^′ =

n₂

Definiamo y il vettore degli stati:

(n₁

)

n₂

T₃

Alcune variabili del motore (rapporti di compressione ed espansione, por-

tate), come già visto, dipendono attraverso delle mappe dagli stati del si-

stema, poiché alcune di esse sono parametrizzate su n₁, mentre altre su T₄,

che è uno stretto parente di T₃. La valutazione di queste grandezze richiede

quindi la soluzione del sistema, ma a meno di utilizzare gli stati calcolati ad

un’istante precedente, essa non è ancora disponibile all’inizio dell’iterazione:

il sistema è quindi implicito.

Generalizzando quindi il problema si presenta nella forma:

y = f(y,y^′,t)

(3.64)

Per risolvere e integrare ad ogni passo un’equazione tipo la 3.64 è neces-

sario un solutore per equazioni differenziali implicite in quanto il sistema è

implicito, poiché lo stato dipende dal valore di se stesso allo stesso istante

(ad esempio la portata in uscita dal plenum è funzione della accelerazione di

rotazione n₁) ed è quindi necessario procedere numericamente risolvendo il

sistema fino a convergenza.

Per quanto riguarda i solutori impiegati, si rimanda al capitolo sull’im-

plementazione.

3.5. PARAMETRIZZAZIONE

DATI NASA

GENERAZIONE MAPPE

NORMALIZZAZIONE CON

SCALATURA

MAPPE ADATTATE

GE T700

VALORI NOMINALI A REGIME

AL NUOVO MOTORE

PARAMETRI

PARAMETRI MOTORE

NUOVO MOTORE

GE T700

Figura

3.10: Diagramma di flusso per il processo di generazione e

adattamento delle mappe delle turbomacchine.

Test Case FUEL FLOW [lb/h] P1 [PSI] T1 [^◦R] P5 [PSI]

140.1

14.37

516.7

14.37

297.2

14.17

515.6

14.43

372.0

14.16

508.3

14.46

458.4

14.09

508.0

14.60

560.4

14.02

507.2

14.63

694.4

13.92

507.2

14.72

Tabella 3.2: Dati ricavati da test NASA su motore GE t700 (prima parte).

3.5

Parametrizzazione

Per la determinazione di alcune grandezze di interesse per la simulazione delle

turbomacchine sono state impiegate delle mappe, messe a punto partendo

da alcuni dati sperimentali ( riassunti nella tabella 3.2 e nella tabella 3.3)

disponibili dal report NASA [5] e relativi a prove effettuate al banco sul

propulsore General Electric t700 nella configurazione utilizzata sull’elicottero

Sikorsky UH-60.

Ogni mappa è ottenuta interpolando i punti sperimentali attraverso una

funzione lineare a tratti (usando la funzione interp1 in linguaggio Matlab o gsl

interp linear in linguaggio C usando le librerie GSL). Si è scelto di utilizzare

Test Case n₁ [%] n₂ [%] g₂₁ [lb/s] p₂₁ [PSI]

T₂₁[PSI] T₄ [R]

65.9

52.6

3.2

58.0

832.0

1413

84.7

95.7

5.16

113.1

1026.0

1577

87.7

95.7

6.16

139.0

1081

1626

90.4

95.7

6.92

161.1

1127

1731

92.6

95.7

7.66

184.8

1173

1838

95.9

95.7

8.5

211.9

1228

1974

Tabella 3.3: Dati ricavati da test NASA su motore GE t700 (seconda parte).

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

una spezzata interpolante passante per i punti considerati per semplicità, ma

l’utilizzo di una funzione interpolante con derivate continue è possibile e di

agevole implementazione.

Tali mappe sono state utilizzate per la determinazione di:

• portata elaborata dal compressore: g₂₁,

• rapporto di compressione del compressore: β_C ,

• rapporto di espansione della turbina di alta pressione: β_HPT ,

• rapporto di espansione della turbina di potenza: β_PT .

I parametri sui quali sono parametrizzati la portata e i rapporti di espan-

sione appena elencati sono stati scelti in base ai dati sperimentali disponibili

e in base al legame che intercorre fra i componenti del motore.

Per le mappe che decidono il funzionamento della prima parte del motore,

ovvero per i componenti calettati sull’albero della turbina generatrice di gas,

si è scelto di far variare le mappe con la velocità di rotazione n₁. Tale scelta

è pressoché obbligata in quanto il parametro che caratterizza il punto di

funzionamento dei componenti compressore e turbina HPT è proprio il loro

regime di rotazione.Si ha quindi:

β_HPT = f(n₁)

β_C = f(n₁)

g₂₁ = f(n₁)

Per la turbina di potenza la scelta della variabile di parametrizzazione è

stata più complessa, in quanto il regime di rotazione dell’albero della turbina

di potenza si deve mantenere costante in ogni condizione di funzionamento

del motore. Tale peculiarità, e la disponibilità di pochi parametri speri-

mentali (riportati tutti nelle tabelle 3.2 e 3.3), hanno portato a scegliere la

temperatura T₄ all’uscita dalla turbina generatrice di gas come parametro

caratterizzante il regime di funzionamento della turbina di potenza.

β_PT = f(T₄)

Le mappe ricavate dai dati NASA sono state adimensionalizzate rispetto

ai valori nominali a regime, cos`ı da renderle estendibili a ogni possibile mo-

tore, noti i suoi parametri caratteristici (si veda a riguardo il diagramma in

figura 3.10). Le dipendenze dei vari parametri sono quindi riferite ai valori

dei parametri normalizzati:

β_HPT = f(n₁)

3.5. PARAMETRIZZAZIONE

β_C = f(n₁)

g₂₁ = f(n₁)

β_PT = f(T₄)

dove i parametri soprasegnati sono normalizzati rispetto al loro valore

nominale:

a_nominal

Per l’estensione delle mappe ad un qualsiasi propulsore devono essere noti dei

parametri che si riferiscano al suo funzionamento nominale, mentre la forma

della mappa viene generata a partire dal funzionamento del motore GE t700

in 6 situazioni di trim riprodotte al banco sperimentalmente.

L’elenco dei parametri da reperire per la scalatura della mappa ad un

particolare propulsore sono:

• n_1nom: valore nominale della velocità di rotazione dell’albero della

turbina generatrice di gas,

• n_2nom: valore nominale della velocità di rotazione dell’albero della

turbina di potenza,

• β_Cnom: valore nominale del rapporto di compressione totale (tutti gli

stadi) del compressore,

• β_HPTnom: valore nominale del rapporto di espansione complessivo della

turbina generatrice di gas,

• β_PTnom: valore nominale del rapporto di espansione complessivo della

turbina di potenza,

• gcomprnom: valore nominale della portata elaborata dal compressore al

regime nominale,

• gb_nom: valore nominale della portata di carburante,

• T4_nom: valore nominale della temperatura in uscita dalla turbina ge-

neratrice di gas.

Nella figura 3.13 è riportata la mappa normalizzata del rapporto di espan-

sione della turbina di potenza: in ascissa è riportata la frazione di tempera-

tura in uscita dalla turbina di potenza.

Nella figura 3.11 sono riportate le mappe del rapporto di compressione

normalizzato del compressore e del rapporto di espansione normalizzato della

CAPITOLO 3. FISICA DEL PROBLEMA E MODELLAZIONE

NORMALIZED MAPS

HPT

0.2

0.4

0.6

0.8

1.2

1.4

Normalized angular velocity n1/n1

NOM

Figura 3.11: Rapporti di compressione/espansione normalizzati.

turbina HPT, entrambi parametrizzati sul regime di rotazione normalizzato

n₁. Si noti come entrambe le mappe passino per il punto (1,1), che è il punto

di funzionamento nominale.

Nella figura 3.12 è riportata la mappa della portata di aria elaborata dal

compressore in funzione di n₁.

3.5.1

Esempio di adattamento delle mappe per un ge-

nerico motore

Nelle figure 3.15 , 3.16, 3.11 si riporta un esempio di adattamento delle mappe

normalizzate a un ipotetico motore, avente i seguenti parametri:

• n_1nom = 75000 RPM,

• n_2nom = 31000 RPM,

• β_cnom = 25.5,

• β_HPTnom = 0.25,

• β_PTnom = 0.58,

• g_{compressorenom} =7kg/s,

Capitolo 4

Implementazione

4.1

Requisiti del codice

L’implementazione del codice è avvenuta tenendo conto dei requisiti previsti

per la sua esecuzione all’interno di un simulatore o comunque di un’interfaccia

di simulazione :

• esecuzione real-time;

• facilità di interfacciarsi con diversi programmi (simulatori, interfacce

grafiche stand-alone);

• possibilità di adattare rapidamente la simulazione all’applicazione mo-

toristica desiderata.

Si è scelto di implementare il primo modello in Matlab, per i notevoli van-

taggi che un linguaggio di alto livello orientato al calcolo numerico comporta

nella prima parte di stesura di un codice di questo genere:

• velocità di scrittura;

• facilità di elaborazione dei dati;

• possibilità di sfruttare funzioni complesse già implementate (integrato-

ri, ecc..);

• documentazione vasta e chiara.

L’implementazione in Matlab presenta diversi limiti:

• difficoltà di interfacciare il codice Matlab con altri software;

• scarsa efficienza computazionale del linguaggio.

CAPITOLO 4. IMPLEMENTAZIONE

Per questi motivi il codice Matlab, già validato, è stato successivamente

tradotto in linguaggio C attraverso la creazione di una DLL (Dinamic Link

Library: libreria a collegamento dinamico). L’implementazione sotto forma

di DLL in C permette di:

• utilizzare uno dei linguaggi di programmazione più efficienti;

• eseguire il codice chiamando la DLL da un qualsiasi altro programma;

• eseguire facilmente il codice su qualsiasi piattaforma;

• caricare selettivamente le varie DLL e da più programmi contempora-

neamente;

• applicare modifiche al modello modificando soltanto la libreria DLL

senza modificare tutti i codici chiamante.

La traduzione del codice dal linguaggio Matlab al linguaggio C sotto for-

ma di DLL ha necessariamente comportato diversi adattamenti, ad esempio

richiedendo la sostituzione delle funzioni native di Matlab (interpolazione at-

traverso INTERP1, POLYFIT, ecc..), con funzioni implementate in C. Dove

possibile si sono utilizzate funzioni facente parte delle GSL (Gnu Scienti-

fic Libraries), librerie di funzioni per il calcolo scientifico ben documentate e

molto vaste, sebbene non intuitive nell’uso come le funzioni native di Matlab.

Dove non è stato possibile utilizzare funzioni o procedure già disponibili, si è

creata la funzione relativa, inclusa poi nella DLL come funzioni di supporto.

4.2

Genesi e struttura del codice Matlab

Il codice implementato in Matlab ha subito diverse modifiche nel tempo,

sia a livello di struttura che a livello di relazioni fisiche implementate. Si

è dapprima implementato un modello fortemente parametrizzato, in cui la

determinazione degli stati termodinamici delle varie sezioni del motore è de-

terminato da alcune mappe ricavate dal funzionamento a regime del solito

motore GE t700. Le varie versioni di questo tipo di modello molto parame-

trizzato sono state identificate con le lettere dell’alfabeto greco: α, β, γ. Il

codice Matlab nella sua versione γ (dipendente da mappe parametriche per la

determinazione di ogni stato termodinamico) è stata validata confrontando

i risultati che si ottenevano dalle simulazioni con dati sperimentali ricavati

dalla letteratura.

Successivamente, validato il codice γ, si è deciso di rendere il modello meno

dipendente da parametri e quindi di modellare le trasformazioni che il fluido

4.2. GENESI E STRUTTURA DEL CODICE MATLAB

subisce all’interno del ciclo termodinamico come espansioni e compressioni

reali. Il codice a cui si è giunti finalmente al termine delle fasi di modellazione,

implementazione, taratura e validazione è stato denominato codice ISO, a

rimarcare il fatto che le trasformazioni che avvengono nelle turbomacchine

sono compressioni o espansioni isoentropiche corrette con un rendimento.

Tale modello è meno vincolato al sistema di partenza da cui sono stati ricavati

i parametri, il motore GE t700, e quindi ci si aspetta che riproduca meno

fedelmente del modello γ il comportamento di tale propulsore, ma sia meglio

adattabile ad altri sistemi.

Nel capitolo riguardante la validazione si confronteranno i risultati otte-

nuti dalle simulazioni effettuate con il codice ISO e con il codice γ e i dati

sperimentali, mettendo quindi in relazione gli errori dei due modelli.

4.2.1

Modello γ - Matlab

Il modello γ presenta una struttura modulare, dove ogni componente del

sistema è rappresentato da una funzione. Le variabili in uscita dalla funzione

rappresentano gli ingressi per un’altra funzione: ad esempio la portata d’aria

in uscita dalla funzione compressore sarà un ingresso per la funzione plenum,

che viene eseguita dopo la funzione compressore. Ordinando in questo modo

il codice è stato possibile isolare i bachi del codice provando separatamente

le varie funzioni. I vari moduli che compongono il modello sono:

• compressore,

• plenum1,

• combustore,

• HPT,

• PT,

• spool 1,

• spool 2.

I moduli che al loro interno implementano una dinamica (Combustore -

T3, Spool 1 - n₁, Spool 2 - n₂), restituiscono in uscita un residuo, ovvero

la derivata dello stato del sistema, e ad ogni iterazione i residui verranno

integrati dall’integratore per fornire i nuovi stati e quindi la nuova soluzione.

Gli stati termodinamici (pressioni, temperature e portate) relativi ai bloc-

chi che implementano le turbomacchine in questo modello vengono generati

CAPITOLO 4. IMPLEMENTAZIONE

DATI:

INGRESSI:

MAPPE

parametri

g_b

condizioni iniziali

g_cool

PRESA D’ARIA

COMPRESSORE

PLENUM 1

COMBUSTORE

HPT

CARICO

SPOOL 2

SPOOL 1

dN1/dt, dN2/dt

TEMPO:

N1,N2

dutata simulazione t_fin

INTEGRATORE

passo temporale dt

Figura 4.1: Schema del codice γ implementato in Matlab.

tramite delle mappe, che forniscano la grandezza desiderata in funzione del

parametro utilizzato nella modellazione. Cos`ı la temperatura, la pressione e

la portata in uscita dal compressore sono funzioni della velocità di rotazione

dell’albero della turbina generatrice di gas attraverso una mappa:

T₂₁ = f(n₁)

p₂₁ = f(n₁)

g₂₁ = f(n₁)

Anche la temperatura e la pressione (o il rapporto di espansione) in uscita

dalla turbina HPT sono funzioni della velocità di rotazione dell’albero della

turbina generatrice di gas:

T₄ = f(n₁)

β_HPT = f(n₁)

Mentre per quanto riguarda la temperatura in uscita dalla turbina PT, e il

rapporto di espansione della turbina PT , essi sono funzioni della temperatura

T4 in uscita dalla turbina HPT:

4.2. GENESI E STRUTTURA DEL CODICE MATLAB

TEMPERATURA USCITA COMPRESSORE

700

650

600

550

500

450

400

Mappe modificate

350

punti LEWIS

retta su dato singolo Lewis

300

spezzata sperimentale

250

0.5

1.5

2.5

3.5

4.5

N1 [giri/min]

x 10⁴

Figura 4.2: Mappa della temperatura in uscita dal compressore in funzione

del regime di rotazione n₁ nel modello γ per il motore GE t700.

T₅ = f(T₄)

β_PT = f(T₄)

Le mappe sono funzioni continue della variabile su cui si è parametrizzato

il dato di interesse. Esse interpolano i dati sperimentali dei 6 casi di trim del

report [5]. Come funzione interpolante è stata scelta una spezzata (funzione

lineare a tratti): si sono confrontati i risultati ottenuti attraverso l’impiego di

altre funzioni (spline cubiche, polinomi di ordine 2,3,4), ma dal momento che

le differenze di funzionamento del codice erano minime si è scelto di utilizzare

semplici spezzate lineari per semplicità della formulazione. Si è anche pen-

sato di utilizzare come funzioni interpolanti dei polinomi che minimizzassero

l’errore dello stato termodinamico in relazione al dato sperimentale (mappe

modificate), opzione che è stata poi scartata perché introduceva nel modello

un grado di aleatorietà che difficilmente avrebbe trovato riscontro in un’ap-

plicazione diversa da quella del motore GE t700 utilizzato per la generazione

delle mappe e poi la valutazione degli errori. Un esempio di mappa per il

modello γ è visibile in figura 4.2, dove è riportata la temperatura in uscita

dal compressore in funzione del regime di rotazione del primo albero.

La modularità del codice gamma è apprezzabile nel listato che ne richiama

in successione i componenti ad ogni iterazione:

CAPITOLO 4. IMPLEMENTAZIONE

% PRESA D’ARIA

p1(i)=p_tot(i);

T1(i)=T_tot(i);

rho(i)=rho_0;

%COMPRESSORE

[p21(i),T21(i),g21(i),P_c(i)]=compressore(p1(i),T1(i),T21(i-1),...

n1(i-1),pg2s,pbetac,pT21,P_pt(i-1),P_f_hpt(i-1));

%PLENUM 1

[p2(i),T2(i),g2(i),g_cool(i),tau_p(i)]=plenum1(p21(i),T21(i),g21(i),...

g_b(i),dt,p2(i-1),V_p1);

%CAMERA DI COMBUSTIONE

[g_gas(i),p3(i),T3(i),dT3(i),iter_comb(i),err_comb(i)]=combustore(p2(i),T2(i)

g_b(i),T3(i-1),dt,LHV,T3max,dT3(i-1));

%TURBINA GENERATRICE DI GAS

[p4(i),T4(i),P_hpt(i),g4(i)] = hpt(p3(i),T3(i),T4(i-1),g_gas(i),...

g_cool(i),T21(i),pT4, n1(i-1),g_b(i),pbetahpt);

%TURBINA DI POTENZA

[p5(i),T5(i),P_pt(i)]=powerturbine(p4(i),T4(i),T5(i-1),g4(i),p_atm,...

n1(i),g_b(i),pbetapt);

%CARICO

[C_ass(i),vel_rotore(i)]=caricoLEWIS(n2(i-1),pcl,g_b(i));

P_load(i)=C_ass(i)*n2(i-1)*0.104719755; %[W]potenza carico

%DINAMICA ALBERI DI TRASMISSIONE

%modellazione perdite per attrito:

P_f_hpt(i)=abs(P_hpt(i)*0.01); %[W]pot perduta in attr.mecc=1/100 Pot hpt

P_f_pt(i)=abs(P_pt(i)*0.01);

%[W]pot perduta in attr.mecc=1/100 Pot pt

%calcolo accelerazioni angolari:

[n1_p(i)]=spool1(P_hpt(i),P_c(i),P_f_hpt(i),J_hpt,n1(i-1));

[n2_p(i)]=spool2(P_pt(i),P_load(i),P_f_pt(i),J_pt,n2(i-1));

%INTEGRAZIONE

%Crank Nicholson

n1(i)=n1(i-1)+0.5*dt*(n1_p(i-1)+n1_p(i));

n2(i)=n1(i-1)+0.5*dt*(n2_p(i-1)+n2_p(i));

Si noti come essendo la formulazione del modello gamma esplicita, sia pos-

sibile integrare le 2 equazioni dinamiche della meccanica del modello espli-

citamente con il metodo dei trapezi (o di Crank-Nicholson). La dinamica

termica infatti, ritenuta molto più veloce delle dinamiche meccaniche, viene

risolta all’interno del componente Bruciatore, con un limite di iterazioni mas-

sime per non sforare oltre il tempo massimo di calcolo per il passo temporale

4.2. GENESI E STRUTTURA DEL CODICE MATLAB

considerato. Nel caso il componente non converga sulla temperatura T₃ nel

numero di iterazioni massimo con la tolleranza assegnata, l’utente ne viene

informato con un segnale di attenzione.

Il modello cos`ı implementato ha fornito buoni risultati se paragonato ai

dati sperimentali disponibili. La convergenza si è sempre rivelata soddisfa-

cente con passi temporali dell’ordine del decimo di secondo, sebbene sia stata

rilevata una certa difficoltà quando venivano forniti ingressi con variazioni

repentine della portata di carburante.

Si raccomanda quindi, qualora interessi implementare il codice del model-

lo gamma, di tarare il FADEC in modo che la derivata temporale della por-

tata di carburante immessa in camera di combustione sia sempre in modulo

minore del valore che impedisce la convergenza del sistema.

4.2.2

Modello ISO - Matlab

Il modello ISO implementa espansioni e compressioni reali per determinare

gli stati termodinamici del fluido dopo il passaggio nelle turbomacchine ed

utilizza delle mappe parametriche solo per la determinazione di portata d’aria

elaborata dal compressore e rapporti di espansione o compressione delle tre

turbomacchine. Il modello numerico è quindi quello del capitolo relativo alla

modellazione, con le dinamiche riassunte nel sistema 3.63.

Le funzioni implementate nel modello ISO sono:

• ENGINEnew2.m: dinamica del motore;

• ENGINEparameters.m: personalizzazione dei parametri caratteristici

del motore;

• NEWMAPS.m: generazione delle mappe normalizzate utilizzate nelle

turbomacchine;

• INITIALconditions.m: assegnamento delle condizioni iniziali per la

simulazione;

• FUELflow.m: assegnamento degli ingressi del motore (portata di com-

bustibile g_b);

• ENGINEload.m: modello di carico che si interfaccia col motore.

L’implementazione della dinamica del modello ISO è avvenuta, a diffe-

renza del modello γ, in un’unica funzione principale denominata ENGINE-

new2.m la cui sintassi è del tipo:

function [ep] = ENGINEnew2(t,y,yp)

CAPITOLO 4. IMPLEMENTAZIONE

Il vettore [ep] raccoglie la differenza fra le derivate degli stati calcolate dalla

funzione (T

3,n˙1,n˙2) e i residui (yp) calcolati al passo precedente e forniti in

ingresso alla funzione insieme al valore degli stati (y). Il vettore [ep] è quindi

composto da:

n₁ = n˙1 - yp(1)

(4.1)

n₂ = n˙2 - yp(2)

(4.2)

e˙T

=T˙₃ - yp(3)

(4.3)

e risulta quindi

[ep] = [

]

n₁; en˙₂; e

Il vettore [ep] è quello che servirà alla funzione integratrice ode15i per iterare

e integrare i residui del sistema.

Le condizioni iniziali di tutte le variabili del motore vengono assegnate

all’interno della funzione INITIALconditions.m; l’assegnamento delle condi-

zioni iniziali deve essere consistente, ovvero deve essere:

[ep] = 0

o in altri termini, facendo riferimento al sistema 3.63:

f (t0, y0, yp0) = 0

Qualora gli stati termodinamici iniziali siano lontani da un punto consistente

per il modello, durante la fase di integrazione si andrà incontro a problemi

di convergenza. Per ovviare a questo problema si può calcolare un punto di

partenza consistente perturbando il modello fino a fare raggiungere al sistema

una condizione iniziale consistente. In Matlab la funzione decic.m permette

di calcolare la condizione iniziale consistente da cui far partire l’integrazione

con le funzioni ode; è possibile specificare alla funzione decic quali stati è

possibile variare e quali invece devono essere mantenuti fissi.

Le mappe relative al funzionamento delle turbomacchine vengono gene-

rate attraverso la funzione NEWMAPS.m: essa genera delle mappe norma-

lizzate a partire dai dati sperimentali disponibili dal report [5]. Tali mappe

andranno poi scalate attraverso i parametri caratteristici del motore consi-

derato nella singola simulazione. Le mappe necessarie al funzionamento del

modello ISO sono:

• mappa per il rapporto di compressione normalizzato del compressore:

β_Compressore = f(n₁)

4.2. GENESI E STRUTTURA DEL CODICE MATLAB

NEWMAPS.m :

mappe normalizzate

INITIALconditions.m :

enginePARAMETERS.m:

condizioni iniziali

parametri caratteristici del

motore considerato

FUELflow.m:

ingressi del sistema

engineLOAD.m:

ENGINEnew2.m:

modello di carico

dinamica motore

TEMPO DI SIMULAZIONE t_fin

ODE15I:

PASSO TEMPORALE dt

integratore

Figura 4.3: Schema delle funzioni implementate nel codice ISO.

• mappa per il rapporto di espansione normalizzato della turbina HPT:

β_HPT = f(n₁)

• mappa per il rapporto di espansione normalizzato della turbina PT:

β_PT = f(T₄)

• mappa per la portata d’aria normalizzata elaborata dal compressore:

g_Compressore = f(n₁)

Per essere utilizzate le mappe devono essere scalate con i parametri ca-

ratteristici del motore considerato. Tali parametri devono essere assegnati

all’interno della funzione ENGINEparameters.m:

• n_1nom: valore nominale della velocità di rotazione dell’albero della

turbina generatrice di gas;

• n_2nom: valore nominale della velocità di rotazione dell’albero della

turbina di potenza;

• J_HPT: momento d’inerzia dell’albero della turbina generatrice di gas,

comprensivo di tutti gli elementi ad esso collegati;

CAPITOLO 4. IMPLEMENTAZIONE

• J_PT: momento d’inerzia dell’albero della turbina di potenza, compren-

sivo di tutti gli elementi ad esso collegati;

• V_CC: volume della camera di combustione;

• β_Cnom: valore nominale del rapporto di compressione totale (tutti gli

stadi) del compressore;

• β_HPTnom: valore nominale del rapporto di espansione complessivo della

turbina generatrice di gas;

• β_PTnom: valore nominale del rapporto di espansione complessivo della

turbina di potenza;

• gcomprnom: valore nominale della portata elaborata dal compressore al

regime nominale;

• gb_nom: valore nominale della portata di carburante;

• T4_nom: valore nominale della temperatura in uscita dalla turbina ge-

neratrice di gas.

La funzione FUELflow.m assegna la storia della portata di carburante per

la durata della simulazione. Nelle prove di validazione della termodinamica

del motore è stata assegnata direttamente una storia temporale alla portata

del carburante, mentre nell’utilizzo del modello all’interno di un codice di

meccanica del volo o di simulazione, tale funzione verrà sostituita da un

regolatore automatico, che in base ad alcuni ingressi regolerà la portata di

carburante con una specifica logica di controllo.

Il carico che si interfaccia al motore è stato implementato nella funzione

ENGINEload.m. Il modello di carico è quello presentato nel capitolo relativo

alla modellazione, dove la potenza richiesta dal carico è proporzionale alla

velocità di rotazione n₂ tramite la 3.57.

4.2.3

Funzioni di supporto

Per il funzionamento del modello ISO sono necessarie alcune funzioni di sup-

porto, che implementano nel modello le proprietà del gas già descritte nel

capitolo relativo alla modellazione.

La funzione aria2.m permette di determinare le proprietà termodinamiche

dell’aria secca. La funzione ha una sintassi del tipo:

[k,cp,cv]=aria2(T)

4.3. STRUTTURA DEL CODICE C

L’help scritto per la funzione ne spiega il funzionamento:

ARIA2 [k,Cp,Cv]=aria2(T) funzione che fornisce le proprietà

termodinamiche dell’aria, assunta gas perfetto,

ricevendo in ingresso la temperatura dell’aria alla

quale si vogliono valutare.

Cp = [J/KgK] calore specifico a pressione costante

Cv = [J/KgK] calore specifico a volume costante

k = [adim] rapporto calori specifici

T = [K] temperatura dell’aria

ARIA2 è ottenuta interpolando dati sperimentali con un polinomio in-

terpolante di ordine 2.

Similmente a quanto fatto per l’aria, si è impiegata una funzione che

calcolasse le proprietà termodinamiche dei gas di combustione. Il modello

impiegato per i gas combusti è quello messo a punto da Guha nel 2000 (ri-

ferimento [12]). La funzione di riferimento per i gas combusti è la funzione

kerosene.m. Essa ha la sintassi:

[k,cp,R] = kerosene(T,FAR)

Gli ingressi della funzione sono la temperatura del gas T e il rapporto

aria/combustibile F AR (Fuel to Air Ratio). Le uscite della funzione sono: k

il rapporto fra calori specifici,c_p il calore specifico a pressione costante (unità

di misura [^J

]), c_v il calore specifico a volume costante, R la costante dei

kgK

gas combusti.

4.3

Struttura del codice C

Il codice C è stato scritto sotto forma di DLL, ovvero tutte le funzioni neces-

sarie al modello sono raccolte all’interno di una libreria a collegamento dina-

mico che andrà poi caricata ed eseguita in fase di compilazione ed esecuzione

della simulazione da un codice chiamante.

Il file a cui si farà riferimento nel seguito è appunto il file sorgente della

libreria DLL EngineModel.c, che al suo interno contiene le funzioni imple-

mentate nel modello e alcuni file header (con estensione .h) che richiamati

nella libreria mettono a disposizione delle funzioni e delle definizioni definite

al loro interno.

Le funzioni contenute nella libreria EngineModel.dll sono:

• EngineModel: funzione principale che implementa la dinamica del

motore;

CAPITOLO 4. IMPLEMENTAZIONE

• kAIR: riceve in ingresso la temperatura e fornisce in uscita il rapporto

fra i calori specifici dell’aria;

• cpAIR: riceve in ingresso la temperatura e fornisce in uscita il calore

specifico dell’aria a pressione costante;

• cvAIR: riceve in ingresso la temperatura e fornisce in uscita il calore

specifico dell’aria a volume costante;

• kGAS:riceve in ingresso la temperatura e il rapporto aria/combustibile

e fornisce in uscita il rapporto fra i calori specifici dei gas combusti di

kerosene e aria con il modello di Guha;

• cpGAS:riceve in ingresso la temperatura e il rapporto aria/combustibile

e fornisce in uscita il calore specifico a pressione costante dei gas com-

busti di kerosene e aria con il modello di Guha;

• MAPgeneration: genera degli oggetti di tipo gsl interp contenenti le

mappe normalizzate per il funzionamento delle turbomacchine;

• GetEngineStates: memorizza il valore degli stati del sistema;

• GetEngineOutput: memorizza in un vettore il valore delle uscite

desiderate;

• EngineModelInit: inizializza i parametri utilizzati dalla funzione

EngineModel;

• EngineModelStep: fa compiere al sistema un passo di integrazione;

4.3.1

Funzione Principale: EngineModel

La funzione principale, che implementa le relazioni fisiche del modello di

motore è EngineModel. La funzione viene dichiarata nel modo seguente:

ENG_DLL_EXPORT int EngineModel(

double dXdt[ENG_STATES_NUM],

double Y[ENG_OUT],

double X[ENG_STATES_NUM],

double U[ENG_CONTROLS],

gsl_interp * pg2n,

gsl_interp * pbetacn,

gsl_interp * pbetahptn,

gsl_interp * pbetaptn,

4.3. STRUTTURA DEL CODICE C

engine_struct *engine_par,

double Pload)

La funzione EngineModel riceve in ingresso:

• il valore dei residui dXdt,

• il valore degli stati X,

• una struttura engine par contenente i parametri caratteristici del mo-

tore,

• il valore degli ingressi U (portata di carburante e portata d’aria spilla-

ta),

• il vettore Y che raccoglie le uscite del sistema,

• Quattro oggetti di tipo gsl interp che contengono le mappe necessarie

alle turbomacchine,

• il valore di potenza Pload richiesta dal carico.

La funzione EngineModel restituisce in uscita una variabile di tipo int

che rappresenta un codice per la buona esecuzione o meno della funzione. In

realtà le uscite della funzione, coerentemente con la sintassi del linguaggio

C, vengono memorizzate nei vettori in ingresso e sono:

• Il vettore dei residui dXdt, contenente le derivate dei tre stati del

sistema,

• Il vettore delle uscite Y, contenente i valori delle variabili interessanti

da esportare a ogni iterazione,

• Il vettore degli stati X.

I residui andranno poi forniti all’integratore insieme al valore degli stati

per determinare i valori degli stati al passo successivo. Le uscite saranno inve-

ce utili nel caso sia necessario valutarne l’andamento temporale, ad esempio

elaborandole in forma grafica.

CAPITOLO 4. IMPLEMENTAZIONE

POTENZA

MAPPE

PARAMETRI

STATI

CARICO

NORMALIZZATE

MOTORE

EngineModel

OUTPUT PER

RESIDUI

ELABORAZIONE

Figura 4.4: Funzione EngineModel.c : ingressi e uscite.

4.3.2

Funzioni di supporto

All’interno della DLL sono implementate diverse funzioni di supporto:

• kAIR

• cpAIR

• cvAIR

• kGAS

• cpGAS

• MAPgeneration

• GetEngineStates

• GetEngineOutput

• EngineModelInit

• EngineModelStep

Tali funzioni consentono di: determinare le proprietà termodinamiche

dei fluidi di lavoro (funzioni kAIR,cpAIR,cvAIR,kGAS,cpGAS ), generare le

mappe necessarie alla simulazione delle turbomacchine (funzione MAPgene-

ration), interfacciarsi con il programma di simulazione o l’eventuale inter-

faccia grafica (funzioni GetEngineStates e GetEngineOutput), avviare la si-

mulazione (funzione EngineModelInit) o farla procedere nel tempo (funzione

EngineModelStep).

4.4. METODO DI SOLUZIONE E DI INTEGRAZIONE

4.4

Metodo di soluzione e di integrazione

Le equazioni che governano il modello sono un sistema di equazioni mi-

ste algebrico differenziali implicite (eq: 3.63), ed è quindi necessario tenere

conto dell’implicità del problema nella scelta del solutore e dell’integratore

temporale.

Linguaggio Matlab Per quanto riguarda l’implementazione in Matlab,

per la scelta dell’integratore si è fatto riferimento alla documentazione del

programma relativa alle equazioni differenziali ordinarie. Matlab implementa

una famiglia di solutori per equazioni differenziali ordinarie: le funzioni ODE

(Ordinary Differential Equations). Fra tutte le funzioni ODE, l’unica adatta

ai problemi impliciti è la funzione ODE15i; essa ha una struttura del tipo :

[TOUT,YOUT] = ODE15I(ODEFUN,TSPAN,Y0,YP0)

Fra le uscite della funzione: TOUT è un vettore temporale che raccoglie gli

istanti a cui sono state calcolate le soluzioni YOUT. Analizzando gli ingres-

si invece: ODEFUN è la funzione che implementa l’equazione differenziale

implicita da integrare, TSPAN l’intervallo temporale in cui integrare, Y0 un

vettore che raccoglie le condizioni iniziali, e YP0 un vettore che raccoglie le

condizioni iniziali delle derivate.

Le condizioni iniziali devono essere consistenti, ovvero

f (T 0, Y 0, Y P 0) = 0

La funzione utilizza un metodo di integrazione di ordine variabile a seconda

della funzione integranda.

É comunque possibile impostare delle opzioni per

la tolleranza o il tipo di integratore con il comando OPTIONS.

Linguaggio C Nell’implementazione del modello in linguaggio C, data la

flessibilità del codice, si è lasciata libertà di scelta dell’integratore da utiliz-

zare. Nella scelta si dovrà tenere conto dell’implicità del sistema e dei tempi

caratteristici di ognuna delle dinamiche considerate. La funzione Engine-

Model fornisce in uscita i tre residui del sistema: i primi due sono relativi

alle dinamiche meccaniche (n₁ ed n₂) ed il terzo è relativo alla dinamica

termica(T₃). A seconda della specifica applicazione considerata i tempi ca-

ratteristici associati ad una o all’altra dinamica varieranno, ad esempio a

seconda delle inerzie associate agli alberi per quanto riguarda le dinamiche

meccaniche, o a seconda del volume della camera di combustione per quanto

riguarda le dinamiche termiche. Si può comunque affermare che, per applica-

zioni su propulsori simili al GE t700, ci si aspetta che la dinamica più veloce

CAPITOLO 4. IMPLEMENTAZIONE

sia quella termica associata alla camera di combustione. La scelta del tempo

di integrazione e del metodo di soluzione andrà quindi effettuata basandosi

sul tempo caratteristico di tale dinamica ed eventualmente sui tempi carat-

teristici delle altre dinamiche che l’integratore dovrà gestire. Ad esempio in

un programma di simulazione della dinamica di un elicottero i tempi carat-

teristici delle dinamiche del motore andranno messi in relazione con i tempi

caratteristici del rotore principale, che a seconda del tipo di modello impie-

gato richiederà un passo di integrazione più o meno breve, che comunque è

in relazione con la frequenza di rotazione del rotore.

Δt_integratore < min(τ_motore, τ_rotore)

La scelta del passo di integrazione è in questo caso una decisione che verrà

presa durante l’integrazione del modello di motore all’interno del programma

di simulazione.

4.5

Interfaccia grafica

Data la flessibilità dell’implementazione del codice in linguaggio C sotto for-

ma di DLL, si è cercato di mettere a punto un’interfaccia grafica che permet-

tesse di simulare i vari motori off-line senza interfacciarsi con nessun codice

di meccanica del volo. Tale interfaccia grafica deve essere abbastanza intui-

tiva da permettere all’utente di ottenere in breve tempo le informazioni sul

funzionamento del motore, avendo inizialmente personalizzato la simulazione

tramite i parametri caratteristici del propulsore che si desidera simulare.

L’interfaccia grafica sarà composta da:

• un pannello principale con i componenti del motore,

• una serie di pannelli che riportino gli andamenti temporali delle gran-

dezze di interesse durante la simulazione,

• un pannello di selezione dei parametri caratteristici del motore da

simulare.

Per l’implementazione dell’interfaccia grafica si è utilizzato il software

LAZARUS. Lazarus è un ambiente ad interfaccia integrata rilasciato e svi-

luppato dal compilatore a codice sorgente aperto Free Pascal. Il linguaggio

con cui viene sviluppata l’interfaccia grafica è quindi il Pascal orientato agli

oggetti.

4.6. DLL - DINAMIC LINK LIBRARY

4.6

DLL - Dinamic Link Library

In informatica, una dynamic-link library (libreria a collegamento dinamico) è

una libreria software che viene caricata dinamicamente in fase di esecuzione,

invece di essere collegata staticamente a un eseguibile in fase di compilazio-

ne. Queste librerie sono note con l’acronimo DLL, che è l’estensione del file

che hanno nel sistema operativo Microsoft Windows, o anche con il termine

librerie condivise (da shared library, usato nella letteratura dei sistemi Unix).

Nei sistemi che usano ELF come formato dei file eseguibili, come ad esempio

Solaris o Linux, sono anche note come .so, abbreviazione di Shared Object.

La separazione del codice in librerie a collegamento dinamico permette

di suddividere il codice eseguibile in parti concettualmente separate, che ver-

ranno caricate solo se effettivamente necessarie. Inoltre, una singola libreria,

caricata in memoria, può essere utilizzata da più programmi, senza la neces-

sità di essere nuovamente caricata, il che permette di risparmiare le risorse

del sistema. Questo metodo di loading on demand consente, inoltre, instal-

lazioni parziali di un sistema software, in cui sono effettivamente presenti

sulla memoria di massa solo le librerie associate alle funzioni che l’utente

desidera utilizzare, come selezionate in fase di installazione: in questo modo

ad esempio viene caricato il sistema operativo Microsoft Windows. Un altro

vantaggio è la possibilità di aggiornare un programma modificando solo le

DLL: inserendo una versione diversa della DLL, che contiene ad esempio dei

bug-fix, tutti i programmi che la usano saranno automaticamente aggiornati

senza bisogno di essere ricompilati. Il principale svantaggio è legato al fatto

che una nuova versione di una DLL potrebbe effettuare dei cosiddetti brea-

king changes, in modo volontario o inconsapevolmente, a causa di bug nella

nuova versione. Un breaking change è una cambiamento critico nel compor-

tamento del codice della funzione che la rende non più compatibile con le

convenzioni in uso (ad esempio, una funzione che prima restituiva NULL in

caso di errore nei parametri e che ora restituisce un valore non nullo). Ancora

più critico il caso in cui un programma di installazione sovrascriva una DLL

con una versione più vecchia.

Le librerie dinamiche vengono caricate dal sistema operativo all’interno

dello spazio di memoria del processo che le ha richieste. In questo modo

l’accesso al codice della DLL avrà prestazioni quasi equivalenti a quelle del

codice dell’applicazione stessa o del codice delle librerie statiche (in seguito si

vedrà perché quasi equivalenti). Per evitare che il codice dell’applicazione e

quello della DLL occupino la stessa posizione in memoria, il linker dovrà pre-

disporre la DLL per la rilocazione. In pratica, il sistema operativo determina

un’area di memoria disponibile e rimappa ogni riferimento alla memoria con-

tenuto nel codice della DLL. Siccome quest’operazione richiede tempo, ogni

CAPITOLO 4. IMPLEMENTAZIONE

DLL dispone di un proprio indirizzo di base ideale: la rilocazione sarà ne-

cessaria solo se a questo indirizzo predeterminato è già stata mappata una

precedente DLL.

4.7

GSL - Gnu Scientific Libraries

Durante l’implementazione in linguaggio C si è fatto ricorso alle librerie open-

source GSL (Gnu Scientific Libraries), che permettono di utilizzare funzioni

di calcolo numerico già implementate e validate. Si è fatto riferimento al

manuale ([14]) per la documentazione relativa ad ognuna delle librerie utiliz-

zate.Le librerie GSL sono librerie software per calcolo numerico molto usate

nella matematica applicata e nelle scienze; sono scritte in linguaggio C, ma

sono disponibili wrappers per altri linguaggi di programmazione. Le GSL

sono parte del progetto GNU Project e sono distribuite sotto la licensa GNU

General Public License.

Le librerie utilizzate in questo lavoro sono:

• Interpolazione: gsl-interp;

• Soluzione di equazioni implicite:

L’inclusione delle librerie considerate viene operata attraverso un file

header nel nostro file DLL, includendo in modo selettivo le singole librerie

necessarie, ad esempio con il comando:

#include <gsl/gsl_interp.h>

Per l’utilizzo delle funzioni relative si rimanda alla documentazione già citata.

Capitolo 5

Validazione

La validazione del modello termodinamico è stata effettuata confrontando i

dati ottenuti dalle simulazioni con i dati sperimentali disponibili dal report

[5], con cui si è confrontata gran parte della letteratura successiva per la

validazione dei propri lavori. Per la validazione della termodinamica del

modello, in anello aperto, si è proceduto in due modi:

• confronto degli stati termodinamici in diverse condizioni di regime;

• confronto della risposta del motore a un ingresso di carburante a gra-

dino;

I dati sperimentali disponibili non coprono tutte le grandezze caratteristiche

del motore e pertanto la validazione è stata possibile solo per le grandezze di

cui si è trovato un riferimento.

5.1

Confronto con dati sperimentali a regime

Per potere confrontarsi con i dati sperimentali del motore GE t700, relativi

ad alcune prove al banco, è necessario simulare le medesime condizioni in cui

si trovava il motore sperimentale. Il confronto è stato possibile solo sui dati

termodinamici disponibili, e quindi il riscontro della buona aderenza o meno

del modello al comportamento del vero propulsore non è completa, ma viene

effettuata solo nel sottospazio degli stati termodinamici disponibili.

Per ottenere una risposta a regime del propulsore si è scelto di simulare

un funzionamento di 100 secondi, con un passo temporale di Δt = 0.1 s, cosı

da essere sicuri che tutti transitori associati alle dinamiche fossero terminati

e che gli stati avessero ormai raggiunto stabilmente il loro valore a regime.

Si è scelto di confrontare i dati sperimentali disponibili e quelli forniti dai

due modelli messi a punto in Matlab: il modello γ e il modello ISO.

CAPITOLO 5. VALIDAZIONE

L’errore percentuale del generico stato termodinamico x ottenuto dalle

simulazioni, è definito come:

x-x

err_% = |

(5.1)

dove x è il valore sperimentale dello stato termodinamico considerato.

Di seguito si riportano le tabelle relative agli errori a regime dei due

modelli per i casi di trim provati nel report [5].

Grandezza Valore sperimentale err_% γ err_% ISO

T₁

286 [K]

0.00

T₂₁

570 [K]

3.91

9.48

T₄

876 [K]

2.28

13.77

g₂₁

2.34 [kg/s]

4.33

23.17

p₁

97702 [Pa]

0.00

p₂₁

779824 [Pa]

3.78

30.14

p₅

99494 [Pa]

4.15

1.80

n₁

37861 [RPM]

2.58

13.54

Tabella 5.1: Errori percentuali dei due modelli rispetto al caso 2 di trim

sperimentale.

Si noti come gli errori relativi agli stati termodinamici T₁ e p₁ siano

sempre nulli, in quanto in condizioni di velocità di volo nulla, coincidono

con le condizioni dell’aria esterna e quindi con le condizioni al contorno che

si sono imposte per il confronto con i dati sperimentali.

L’errore più marcato è quello riscontrato nella pressione di scarico dei

gas combusti, all’uscita della turbina di potenza. Tale errore, in percentuale

anche molto elevato (ad esempio nel caso del test case 6 raggiunge ancheil

50%), è paragonabile in modulo con gli errori commessi nel valutare le altre

pressioni, ma poichè la pressione p₅ è prossima alla pressione atmosferica,

l’errore percentuale risulta molto più marcato rispetto a quello che si ha ad

esempio nella valutazione di p₂₁, che è di un ordine di grandezza più elevato.

É possibile ridurre l’errore su p₅ modificando il valore dei rendimenti isoen-

tropici delle turbomacchine. Si è scelto comunque di mantenere i rendimenti

ai seguenti valori:

• η_C = 0.96

• η_HPT = 0.96

• η_PT =0.86

5.2. VALIDAZIONE DINAMICA: RISPOSTA AD UN GRADINO DI PORTATA DI COMBUSTIBIL

• η_CC = 0.98

perchè questi sono i valori riportati più spesso in letteratura per il motore

GE t700.

Gli errori commessi nelle simulazioni sono riportate nelle tabelle 5.1, 5.1,

5.1, 5.1, 5.1.

Grandezza Valore sperimentale err_% γ err_% ISO

T₁

282 [K]

0.00

T₂₁

600 [K]

1.09

2.42

T₄

903 [K]

1.56

3.59

g₂₁

2.79 [kg/s]

2.94

0.59

p₁

97633 [Pa]

0.00

p₂₁

958405 [Pa]

0.28

0.79

p₅

99494 [Pa]

1.78

27.69

n₁

39201 [RPM]

0.73

0.15

Tabella 5.2: Errori percentuali dei due modelli rispetto al caso 3 di trim

sperimentale.

Grandezza Valore sperimentale err_% γ err_% ISO

T₁

282 [K]

0.00

T₂₁

626 [K]

0.70

1.04

T₄

961 [K]

1.11

1.56

g₂₁

3.13 [kg/s]

2.12

4.29

p₁

97150 [Pa]

0.00

p₂₁

1110785 [Pa]

0.98

5.53

p₅

1100667 [Pa]

3.01

27.07

n₁

40408 [RPM]

0.41

1.17

Tabella 5.3: Errori percentuali dei due modelli rispetto al caso 4 di trim

sperimentale.

5.2

Validazione dinamica: risposta ad un gra-

dino di portata di combustibile

La grandezza che rappresenta l’ingresso principale del sistema, attraverso cui

si può controllare il regime di rotazione del motore, è l’afflusso di carburan-

te g_b: una validazione dinamica può essere effettuata rilevando la risposta

CAPITOLO 5. VALIDAZIONE

Grandezza Valore sperimentale err_% γ err_% ISO

T₁

282 [K]

0.00

T₂₁

652 [K]

1.46

0.13

T₄

1021 [K]

0.28

0.35

g₂₁

3.47 [kg/s]

3.61

7.46

p₁

96667 [Pa]

0.00

p₂₁

1274196 [Pa]

4.53

10.24

p₅

100873 [Pa]

1.38

37.43

n₁

41392 [RPM]

1.15

1.83

Tabella 5.4: Errori percentuali dei due modelli rispetto al caso 5 di trim

sperimentale.

Grandezza Valore sperimentale err_% γ err_% ISO

T₁

282 [K]

0.00

T₂₁

682 [K]

6.38

1.15

T₄

1096 [K]

4.86

1.41

g₂₁

3.86 [kg/s]

13.47

9.84

p₁

95978 [Pa]

0.00

p₂₁

1461051 [Pa]

19.4

13.43

p₅

101494 [Pa]

21.13

50.52

n₁

42867 [RPM]

4.88

3.42

Tabella 5.5: Errori percentuali dei due modelli rispetto al caso 6 di trim

sperimentale.

5.2. VALIDAZIONE DINAMICA: RISPOSTA AD UN GRADINO DI PORTATA DI COMBUSTIBIL

del sistema a una perturbazione di questo ingresso. Nel report NASA [5],

la validazione dinamica del modello di motore in anello aperto viene effet-

tuata confrontando i dati ottenuti dalla risposta ad un gradino di portata

di combustibile del modello semplificato con quella di un simulatore ad alta

fedeltà proprietario del produttore del propulsore. Essendo i dati del report

NASA disponibili solo in formato grafico, si è proceduto dapprima digitaliz-

zando i grafici e poi elaborandoli per ottenere dei vettori che sono poi stati

confrontati con i risultati ottenuti dalla simulazione.

La simulazione è stata effettuata a parità di condizioni iniziali, di tempi e

di passo di integrazione. Si riportano di seguito le condizioni di simulazione:

• durata temporale t = 5s;

• passo di integrazione Δt = 14 ms;

• afflusso di carburante iniziale g_b

= 400 lb /h = 0.050 399 16 kg s^-1;

iniziale

• afflusso di carburante finale gbstep = 775 lb /h = 0.097 648 37 kg s^-1;

• momento di applicazione del gradino di afflusso di carburante t = 0.5s;

Le condizioni iniziali da cui si è fatta partire la simulazione sono state:

Temperatura Valore [k]

T₁

288.15

T₂₁

633

T₃

1171

T₂

633

T₄

900

T₅

730

Tabella 5.6: Condizioni iniziali di temperatura nelle varie sezioni del motore

prima dell’ingresso a scalino di portata di combustibile.

L’andamento dell’ingresso a gradino è visibile in figura 5.1, mentre il

confronto della risposta è riportato nelle figure 5.2, 5.4,5.5, 5.3.

Si noti la buona riproduzione del comportamento della temperatura di

uscita dalla camera di combustione (figura 5.4), che è uno dei tre stati del si-

stema. La buona riproduzione della dinamica del componente indica una cor-

retta modellazione, che riesce a catturare, nonostante la relativa semplicità,

le dinamiche tipiche con un errore modesto.

Nelle altre figure è possibile notare come lo scostamento dei valori sia

più marcato, sebbene le forme delle risposte siano riprodotte abbastanza

5.2. VALIDAZIONE DINAMICA: RISPOSTA AD UN GRADINO DI PORTATA DI COMBUSTIBIL

T3: COMBUSTOR OUTLET TEMPERATURE

1600

NASA

ISO model

1550

1500

1450

1400

1350

1300

1250

1200

1150

Time [s]

Figura 5.4: Risposta ad un ingresso a gradino di portata di combustibile:

confronto dellla temperatura T₃ in uscita dalla camera di combustione con i

dati del report Ballin - 1988.

T4: HIGH PRESSURE TURBINE OUTLET TEMPERATURE

1250

NASA

ISO model

1200

1150

1100

1050

1000

950

900

850

Time [s]

Figura 5.5: Risposta ad un ingresso a gradino di portata di combustibile:

confronto dellla temperatura T₄ in uscita dalla turbina di alta pressione con

i dati del report Ballin - 1988.

Capitolo 6

Controllore automatico del

propulsore

6.1

Architettura e logica di controllo

Il motore a turbina per elicotteri è generalmente regolato da un controllore

automatico, che ne garantisca la funzionalità entro dei limiti di sicurezza. I

sistemi di controllo moderni sono di tipo digitale, senza possibilità di passare

in modalità manuale: essi prendono il nome di FADEC (Full Authority

Digital Engine Control) proprio perchè il controllo del motore è totalmente

automatico. Il FADEC a sua volta si compone di due componenti principali:

ECS e FCU.

Il sistema di controllo del motore ECS (Engine Control System) ha il

compito di regolare la potenza del motore in funzione delle necessità operative

dell’elicottero. Al variare delle condizioni operative della macchina esso deve:

• regolare la potenza erogata;

• gestire la variazione di flusso di combustibile;

• limitare i regimi di rotazione;

• consentire avviamento e spegnimento;

L’ FCU (Fuel Control Unit) gestisce l’afflusso di carburante al motore e

su comando dell’ECS trasforma il segnale elettrico del controllore in flusso

di carburante attraverso una pompa elettrica.

La figura 6.1 mostra l’architettura del sistema di controllo FADEC nella

sua struttura più generale. Il sistema elabora diversi tipi di informazioni:

88CAPITOLO 6. CONTROLLORE AUTOMATICO DEL PROPULSORE

Cyclic pitch control

Collective pitch control

FADEC

Info from other engines

FCU

ECS

Environmental conditions:

p0,T0

fuel flow

Engine mode:

OFF

IDLE

FLIGHT

ENGINE

Engine states:

MAX

N1, N2, T3, T4, P21

Figura 6.1: Architettura di un sistema di controllo di un motore a turbina

per elicotteri.

• grandezze misurate da sensori posti sul motore: velocità di rotazione

n₁ e n₂, pressione p₂₁ a valle del compressore, temperatura T₃ a monte

della turbina HPT e temperatura T₄ tra le due turbine;

• comandi imposti dal pilota: passo ciclico e passo collettivo;

• condizioni ambientali esterne: p₀, T₀.

Per ogni situazione operativa il FADEC calcola la potenza richiesta dal ca-

rico, la velocità di rotazione necessaria all’albero della turbina generatrice

di gas (n_1desiderata) per generare tale potenza al secondo albero e il flusso di

combustibile da inviare in camera di combustione per variare della quantita

Δn₁ la velocità di rotazione del primo albero.

Il controllore impiegato nella maggior parte dei FADEC è di tipo PID,

ovvero con azione proporzionale, integrativa e derivativa.

Il controllore acquisisce in ingresso un valore dal processo, e lo confronta

con un valore di riferimento. La differenza, il cosiddetto segnale di errore

e(t), viene quindi usata per determinare il valore della variabile di uscita del

controllore, che è la variabile manipolabile del processo. Il PID regola l’uscita

in base a:

• il valore del segnale di errore (azione proporzionale);

• ivaloripassatidelsegnaledierrore(azioneintegrale);

6.1. ARCHITETTURA E LOGICA DI CONTROLLO

TIPO

K_P

K_I

K_D

0.5K_U

0.45K_U

0.54KU

T_U

PID

0.6K_U

0.9KU

0.075K_U T_U

T_U

Tabella 6.1: Taratura dei guadagni con il metodo di Ziegler e Nichols.

• quanto velocemente il segnale di errore varia (azione derivativa).

Chiamando u il generico vettore di controllo, esso viene fornito in uscita

dal regolatore PID come:

∫

k_i

du(t)

u(t) =

u(t)dt + k_P u(t) + k_dT_d

(6.1)

T_i

Oppure, chiamando K_I =ki

, K_D = k_dT_d, nella più comune forma:

T_i

∫

du(t)

u(t) = K_I

u(t)dt + K_P u(t) + K_D

(6.2)

dove le costanti K_I , K_D, K_P sono i guadagni dell’azione integrale, deriva-

tiva e proporzionale; l’equazione 6.2 può contenere tutti i tre termini a destra

o solo alcuni di essi a seconda del tipo di controllore che si vuole implementare

(si veda a riguardo il riferimento [15]).

La taratura dei guadagni può avvenire attraverso il secondo metodo di

Ziegler e Nichols ( spiegata nei riferimenti [15, 16]):

1. Si utilizza in primo luogo un controllo solo proporzionale, ovvero con

K_P = 0, K_I = K_D = 0.

2. Si aumenta il guadagno K_P finché nel sistema non si instaurano oscil-

lazioni sostenute, che cioè non spariscono dopo un transitorio: questa

è una misura dell’effetto dei ritardi e della dinamica del processo.

3. Il valore minimo di K_P che provoca le oscillazioni sostenute del sistema

viene denominato guadagno critico K_U .

4. Il periodo delle oscillazioni sostenute del sistema è il periodo critico T_U .

5. Noti il periodo critico e il guadagno critico è possibile determinare il

valore dei guadagni attraverso la tabella 6.1.

Per avere un’idea dell’ordine di grandezza del guadagno dell’azione pro-

porzionale si è messo in relazione il massimo valore dell’errore, che si avrebbe

90CAPITOLO 6. CONTROLLORE AUTOMATICO DEL PROPULSORE

Collective Pitch Control

LOOK UP TABLE

N1des

PID controller

FUEL FLOW

(N2-N2des)

LOOK UP TABLE

deltaN1

Figura 6.2: Logica di controllo del regolatore PID implementato all’interno

del FADEC.

quando il secondo albero fosse fermo mentre la velocità desiderata fosse quella

nominale, e il massimo valore della portata di combustibile. L’errore massimo

sul regime n₂ risulta:

err_MAX = n₂

-n₂ =n2nom -0=n2nom

(6.3)

des

La determinazione del guadagno K_P può avvenire come:

g_b

M AX

K_P

(6.4)

iniziale

n₂

M AX

Nel caso del motore per aeromobili, la variabile manipolabile è la portata

di combustibile diretta in camera di combustione (g_b o portata di combu-

stibile), mentre l’errore viene calcolato sul regime di rotazione dell’albero di

potenza, la 6.2 va quindi valutata assegnando:

u(t) = g_b(t)

(6.5)

e(t) = n₂(t) - n₂

(6.6)

ref

Dove n₂

è il valore del regime di rotazione da mantenere per la determinata

ref

condizione di volo. Ovviamente per raggiungere il valore di n₂

agendo sul

ref

flusso di combustibile il motore reagirà accelerando n₁ e poi trasferendo il

moto al secondo albero, l’azione è quindi indiretta.

É possibile implementare

delle tabelle (look up table) che dato il valore di Δn₂ e di comando del passo

ciclico, forniscano il valore di n₁ da far raggiungere al motore.

La logica di controllo del PID implementato all’interno di un FADEC in

realtà tiene anche conto dei comandi di passo ciclico e passo collettivo che il

pilota applica al velivolo, aggiornando il valore di n₁

da inseguire in base a

ref

questi dati. La logica di controllo qu`ı descritta, implementata in un generico

FADEC ha lo schema riportato in figura 6.2.

6.1.1

Limiti di impiego del motore

Come detto il FADEC deve regolare la portata di combustibile gb cosı da

mantenere il funzionamento del motore all’interno dei suoi limiti operativi.

6.2. IMPLEMENTAZIONE

Le prestazioni che vengono richieste al motore dal pilota o dai sistemi di ge-

stione del volo devono essere compatibili con numerosi vincoli, sia del velivolo

(e delle parti che lo compongono: sistema di trasmissione, rotore, ecc..) che

del motore stesso. Il sistema di controllo deve assicurare che l’utilizzo del

motore avvenga in ogni condizione nel rispetto di questi limiti, il FADEC

deve infatti:

• garantire che il regime di rotazione n₂ si mantenga in prossimità del

valore nominale in tutte le condizioni operative;

• limitare il regime di rotazione n₁ sia inferiormente (per non avere troppo

ritardo nella successiva accelerazione) che superiormente (overspeed:

regimi troppo elevati potrebbero portare i componenti soggetti alla

forza centrifuga a superare i propri limiti meccanici);

• garantire la composizione ottimale della miscela aria/combustibile sia

in accelerazione (rischio di surge: fiammate che sottopongono i com-

ponenti meccanici a temperature troppo elevate), che in decelerazione

(rischio di flame out: spegnimento della fiamma);

• limitare la potenza P_PT erogata dal propulsore, per non sottoporre ad

eccessive sollecitazioni le componenti a valle del motore;

• limitare la coppia del motore per non portare a rottura la trasmissione;

• limitare la temperatura di ingresso in turbina per non sottoporre le

palette della turbina a un eccessivo stress termico.

6.2

Implementazione

Come effettuato per la parte termodinamica, l’implementazione del control-

lore è avvenuta prima in linguaggio Matlab, e successivamente in linguaggio

C, sfruttando la facilità di scrittura ed elaborazione dei dati di Matlab.

Il controllore del propulsore è, dal punto di vista dell’implementazione,

assimilabile ad un sistema che, noti i parametri del sistema-motore (i suoi

stati termodinamici) e noti gli ingressi (ad esempio un target di potenza da

fornire al carico o un comando di passo ciclico del pilota), genera l’ingresso

di carburante al motore (si veda a riguardo la figura 6.3). L’unica differenza

fra l’implementazione del motore senza controllore e l’implementazione del

motore controllato è quindi la presenza di un sistema (implementato attra-

verso una funzione) che varia il flusso di carburante invece di assegnargli una

storia temporale come veniva fatto ad esempio con la funzione FUELflow.m

92CAPITOLO 6. CONTROLLORE AUTOMATICO DEL PROPULSORE

Y_des

Err

CONTROLLORE

MOTORE

CARICO

Figura 6.3: Sistema motore-controllore in anello chiuso.

nel modello ISO. Ovviamente la presenza di un sistema dinamico come il con-

trollore, modificherà la risposta del sistema complessivo, composto da motore

e controllore.

6.2.1

Implementazione in linguaggio Matlab

Non dovendo il codice scritto in linguaggio Matlab interfacciarsi con alcun

simulatore di meccanica del volo, il controllore sviluppato per il modello ISO

non ha a disposizione le grandezze tipiche con cui un FADEC regola tipi-

camente il suo funzionamento: i comandi di passo ciclico e passo collettivo

forniti dal pilota. Si è quindi implementato un controllore di tipo PID, do-

ve l’obiettivo che viene richiesto è quello di mantenere il regime di rotazione

dell’albero di potenza in prossimità del regime nominale. In una implementa-

zione all’interno di un simulatore sarebbe naturalmente necessario elaborare

i comandi di passo ciclico e passo collettivo per poi gestire l’afflusso di car-

burante in modo che la potenza dei propulsori sia tale da mantenere ad ogni

istante il rotore in rotazione al regime nominale.

Il controllore in Matlab è una funzione del tipo

[gb]=FADEC(Parametri)

che regola il flusso di combustibile da inviare al motore in base al valore che

assumono alcune variabili interessanti per il controllo. Queste variabili con-

sentono al controllore di valutare lo stato in cui si trova ad operare il motore,

e permettono di decidere la portata di carburante da inviare in camera di

combustione.

É richiesto quindi all’utente interessato all’implementazione

di un motore in anello chiuso, di fornire non soltanto quei parametri uti-

li alla simulazione della parte termodinamica del motore, ma anche alcuni

limiti meccanici o termici del sistema da simulare, che saranno necessari a

impostare la logica di controllo del FADEC. Ad esempio per l’applicazione

in un simulatore di meccanica del volo di un elicottero, bisognerà specificare

al controllore:

6.2. IMPLEMENTAZIONE

VINCOLI:

gb_max - gb_min

n1_max -n1_min

n2_max - n2_min

T3_max

Torque_max

Power_max

STATI:

n1, n2, T3

FUEL FLOW:

PID

RIFERIMENTO:

n2_ref

PASSO TEMPORALE:

Figura 6.4: Struttura del controllore FADEC nell’implementazione in

linguaggio Matlab.

• la massima potenza (P_PT

) (o coppia (C_MAX )) trasferibile alla tra-

M AX

smissione,

• il minimo (n₂

) e il massimo regime (n₂

) di rotazione sopportato

MIN

M AX

dal carico,

• il minimo (n₁

) e il massimo (n₁

) regime di rotazione ammesso

MIN

M AX

per l’albero della turbina generatrice di gas: per avere un ritardo ac-

cettabile nella successiva accelerazione del motore il primo, e per non

sforare i limiti strutturali con le forze centrifughe il secondo,

• il minimo (g_b

) e il massimo (g_b

) afflusso di combustibile ammesso

MIN

M AX

in camera di combustione,

• la massima temperatura T₃ accettabile in camera di combustione.

Nell’implementazione la taratura dei guadagni del regolatore PID è stata

effettuata attraverso il secondo metodo di Ziegler e Nichols, andando poi

a modificare leggermente i valori forniti dalla tabella 6.1 per privilegiare la

prontezza della risposta o l’annullamento più rapido possibile dell’errore di

n₂.

94CAPITOLO 6. CONTROLLORE AUTOMATICO DEL PROPULSORE

6.3

Simulazioni del motore in anello chiuso

Si riportano gli andamenti degli stati termodinamici di alcune simulazioni

effettuate con il modello ISO di motore controllato da un PID in ambiente

Matlab. Al propulsore è collegato un carico, la cui potenza è descritta dalla

3.57. L’obiettivo del controllore è quello di mantenere il carico in rotazione

a una velocità prossima a quella del riferimento n₂

REF

6.3.1

Simulazione 1 - controllore proporzionale

I parametri caratteristici della simulazione 1 sono:

• tempo di simulazione: 5 secondi;

• passo di integrazione: 0.01 secondi;

• tipo di controllore: Proporzionale;

• guadagno: K_P = 4.65 10^-3;

• target di velocità di rotazione: n₂ = n₂

= 20900[RP M];

nominale

La simulazione è stata effettuata accoppiando il modello di carico che

simula in modo semplificato il rotore dell’elicottero Sikorsky UH-60 ad un

solo motore. Tale tipo di utilizzo del motore è simile a quello che si avrebbe

durante un guasto di uno dei due propulsori, situazione OEI (One Engine

Inoperative), in cui la potenza necessaria alla macchina deve essere fornita

interamente da una sola delle due unità propulsive. Tale condizione è molto

gravosa per il motore e ci si aspetta che il consumo sia superiore a quello che

si ha durante il funzionamento normale.

Le condizioni iniziali della simulazione 1 sono riportati nella tabella 6.3.1.

Osservazioni L’aggiunta del controllore e della sua dinamica ha richiesto

nell’implementazione in Matlab con la funzione integratrice ode15i di ridurre

il passo di integrazione a Δt = 0.01 s per non avere problemi di convergenza.

Tale passo di integrazione è dell’ordine di grandezza di quello richiesto dalla

dinamica del rotore in un simulatore che implementi un modello di rotore

con grandezze mediate sul giro, ma è inferiore di un ordine di grandezza di

quello richiesto dalla sola dinamica termica-meccanica del motore in anello

aperto.

6.3. SIMULAZIONI DEL MOTORE IN ANELLO CHIUSO

Variabile Valore Unità di misura

n₁

44700^π

rad/s

n₂

20000^π

rad/s

p₁

95978

T₁

281.778

p₂₁

95978

T₂₁

281.778

p₂

95978

T₂

281.778

p₃

95978

T₃

281.778

K/s

p₄

98978

T₄

281.778

p₅

95978

T₅

281.778

P_c

P_hpt

P_pt

P_fpt

P_fhpt

g₂₁

Kg/s

g₂

Kg/s

g₃

Kg/s

g₄

Kg/s

τ_CC

n₁

rad s^-2

n₂

rads^-2

K/s

Tabella 6.2: Condizioni iniziali della simulazione 1 di motore controllato in

anello chiuso da controllore proporzionale.

96CAPITOLO 6. CONTROLLORE AUTOMATICO DEL PROPULSORE

x 10⁴

N2: POWER TURBINE SHAFT SPEED

2.1

NOMINAL

1.9

x 10⁴

N1: GAS TURBINE SHAFT SPEED

4.5

NOMINAL

4.4

4.3

FUEL FLOW TIME HISTORY

0.2

0.1

time [s]

Figura 6.5: Simulazione 1 con controllore P. Velocità di rotazione dei due

alberi e variabile di controllo g_b.

N2: POWER TURBINE SHAFT SPEED

x 10⁴

2.1

2.08

2.06

2.04

REFERENCE

2.02

1.98

1.96

time [s]

Figura 6.6: Simulazione 1 con controllore Proporzionale.Velocità angolare

dell’albero della turbina di potenza. Confronto con il riferimento da inseguire.

6.3. SIMULAZIONI DEL MOTORE IN ANELLO CHIUSO

Utilizzo di un controllore proporzionale L’utilizzo di un controllore

solo proporzionale (di tipo P) con un guadagno adeguato fornisce prestazioni

piuttosto soddisfacenti. L’andamento della variabile di controllo nella simu-

lazione 1 mostra una sovraelongazione di lieve entità e si attesta sul valore

di regime dopo circa 0.18 s. La variabile da controllare raggiunge dopo circa

0.15 s il valore del riferimento, mentre la velocità dell’albero della turbina

generatrice di gas (n₁), si attesta su un valore di poco inferiore al valore

nominale. Il non annullarsi dell’errore a regime è coerente con il tipo di con-

trollore scelto: un controllore solo proporzionale non è capace di annullare

l’errore a regime, ma è pronto nell’inseguire il riferimento. Il consumo di car-

burante per mantenere in rotazione il rotore con un solo propulsore è, come

ci si aspetta, è superiore a quello che si avrebbe in condizioni nominali. Il

surplus di portata di carburante richiesto rispetto alle condizioni nominali è

di circa il 33%.

6.3.2

Simulazione 1 - controllore PI

La simulazione 3 riprende esattamente i dati iniziali della simulazione 1, ma

con l’impiego di un controllore di tipo PI: proporzionale e integrativo. Tale

controllore permette l’annullamento dell’errore a regime. Per la taratura dei

guadagni si è partito dai valori ottenuti dalla regola di Ziegler e Nichols, va-

riando poi il guadagno della componente integrativa perchè il valore calcolato

con la regola della tabella 6.1 provocava oscillazioni sostenute.

Si noti (figura 6.7) come il controllore con azione proporzionale e inte-

grativa riesca ad annullare l’errore a regime, a differenza del controllore solo

proporzionale. L’azione integrativa provoca un ritardo nell’inseguimento del

riferimento: la velocità di rotazione n₂ raggiunge il suo valore target dopo

circa 0.4 s (figura 6.8).

In figura 6.9 è riportato l’andamento delle due velocità di rotazione e della

variabile di controllo g_b, che rappresenta l’uscita del controllore.

6.3.3

Simulazione 1 - controllore PID

La simulazione 1 condotta utilizzando un controllore PID, che impieghi tutti

i tre contributi proporzionale, integrativo e derivativo, fornisce come risultati

le storie temporali visibili nelle figure 6.10, 6.11 e 6.12.

É possibile notare (figura 6.11) come già dopo 0.25s la velocità n₂ sia

giunta già in prossimità del riferimento.

In figura 6.12 si apprezza come nei primi istanti della simulazione la varia-

bile di controllo cambi repentinamente, a causa dell’azione derivativa. Tali

Modellazione e implementazione numerica di un motore a turbina per elicottero (2012-2019) - 1