Metodi teorico-numerici per l’identificazione di un impatto su una struttura aerospaziale. Correlazione con dati sperimentali.

Introduzione
Durante il primo volo del lanciatore europeo VEGA uno shock è stato rilevato dagli accelerometri posti sull’ugello del motore per il primo stadio a circa 20 secondi dopo l’accensione. Tra le ipotesi fatte per giustificare il fenomeno vi era la possibilità di un impatto sulla parete dell’ugello da parte di un oggetto interno, come un pezzo di protezione termica separatosi dalla sede. Questo lavoro, su richiesta di AVIO S.p.A., capocommessa del lanciatore VEGA, aveva lo scopo di scartare tale ipotesi oppure dimostrare la sua fattibilità sviluppando una procedura che determinasse le caratteristiche di impatto: livello di eccitazione e posizione.

Analisi preliminari
Per far ciò si è partiti da analisi preliminari sui dati sperimentali disponibili in modo da fornire un’idea di massima sui risultati da ottenere e poter giudicare la bontà della procedura successivamente sviluppata sfruttando diverse potenzialità del software FEM MSC.Nastran.
Dapprima si è fatto diretto riferimento agli 11 canali accelerometrici di cui fossero disponibili i dati riferiti all’anomalia (specificati nelle rispettive posizioni in Figura 1): confrontando ad esempio le storie temporali delle accelerazioni assiali sui due martinetti elettromeccanici (Figura 2) alla loro giunzione con l’ugello, si riconosce che il martinetto EMA1B (sulla sinistra dell’ugello in Figura 1) sperimenta in pieno lo shock legato all’anomalia prima di quanto accada sull’altro, di modo tale da poter supporre che l’impatto sia difatti avvenuto nella regione sinistra del divergente dell’ugello.

 

In queste come nelle altre misure si riconosce che l’anomalia inizia a circa 20.6 s dopo l’accensione e dura per un minimo di 60 ms fino a un massimo di 100 ms in alcune uscite.
Si osserva esplicitamente che qui e nel seguito tutti i dati riportati sono stati normalizzati per motivi legati al segreto industriale.

Procedura per identificazione impatto: analisi di sensitività
Per identificare le regioni di impatto più plausibili data la risposta meccanica rilevata si è ricorsi ad una analisi di sensitività delle energie dei segnali di accelerazione rispetto al carico di impatto; dopodiché tale carico è stato ottimizzato per massimizzare la correlazione numerico-sperimentale delle energie. Entrambi gli approcci richiedono l’applicazione iterativa di analisi FEM transitorie modali.
Al fine di rendere la procedura iterativa più rapida possibile sono stati impiegati tre accorgimenti:
limitare la banda di estrazione degli autovalori, pre-filtrando i dati sperimentali nella banda 10-400 Hz;
fornire un tipo di carico, nell’ambito di una analisi lineare, assimilabile all’effettivo fenomeno di impatto;
ridurre il modello del grano propellente, particolarmente dispendioso in termini di costo computazionale per analisi modali.
Per il secondo punto, è chiaro che la soluzione migliore è quella di ricorrere ad un impulso. Questo è stato modellizzato con una massa concentrata (scheda CONM2 in MSC.Nastran) inizialmente unitaria, nel punto supposto per l’impatto, ed agente su di essa un carico di accelerazione (ACCEL1) con una storia temporale ottenuta, in termini di forma e durata, da una preliminare simulazione non lineare esplicita dell’impatto effettuata con il solutore MSC.Nastran SOL 700, assumendo che l’oggetto impattante fosse un pezzo di protezione termica del motore, ossia un elastomero. Il ruolo della CONM2 è di scalare il valore del carico applicato soprattutto in ottica di ottimizzazione della correlazione numerico-sperimentale, affinché il fenomeno sia riprodotto al meglio dati i risultati sperimentali.
Per il terzo punto invece si sono testati in un’analisi transitoria modale di prova tre diversi modelli: l’originale, un modello di grano condensato dinamicamente con l’approccio di Craig-Bampton, ed uno senza il grano ma solo la sua massa mediante l’impiego di masse non strutturali. Laddove i risultati in termini di energie nei singoli canali sono vicini (Figura 3), la notevole differenza nei tempi computazionali ottenuti (Tabella 1) ha spinto per l’adozione dell’ultimo modello.

 

Modello	Tempo S
Completo	1594
Grano condensato	517
Masse non strutturali	74

Per limitare le possibili posizioni di impatto è stata eseguita l’analisi di sensitività delle energie dei segnali di accelerazione in output rispetto al carico applicato. Dapprima dodici posizioni sono state considerate, divise in tre fasce assiali e quattro angolari uniformemente distribuite sul cono del divergente dell’ugello, e per ognuna di esse 3 distinte componenti di carico: carico tangenziale alla parete dell’ugello, normale ad essa, e circonferenziale. Per quanto detto sul ruolo della massa concentrata unitaria CONM2, l’analisi di sensitività rispetto al carico è stata condotta in realtà in modo lecito rispetto a tale massa.

Si è proceduto dapprima da un punto di vista teorico: partendo dalla formulazione matematica del problema agli elementi finiti è possibile mostrare che esiste la relazione (1) tra sensitività,

energia del segnale di accelerazione  e massa concentrata  per ciascun canale. Sfruttando questa e le analisi transitorie modali, sostenute col solutore MSC.Nastran SOL 112, si ottengono i risultati di sensitività per ognuna delle 36 condizioni di carico.
Per validare l’analisi teorica sviluppata si è studiata la sensitività anche in maniera numerica, dapprima con uno schema alle differenze finite, Equazione (2), implementato in Matlab® sfruttando esternamente la SOL 112 per i calcoli FEM nelle configurazioni originale e incrementata.

                                                                          (2)

 

 

 

 

I risultati che si ottengono sono uguali a quelli teorici entro un errore medio del 13%, che risulta accettabile. I 36 coefficienti ottenuti per uno dei canali sono riportati in Figura 4.

 

 

In generale si può riconoscere come la sensitività sia sempre molto più alta rispetto alla componente normale dell’impatto (N) che non per le altre due (T tangenziale e C circonferenziale) per ogni posizione ed ogni canale di output. Pertanto nel seguito è stato considerato un impatto normale, immaginando di trascurare gli effetti delle altre due componenti ai fini dell’anomalia rilevata.
Un modo più veloce ed elegante di svolgere l’analisi di sensitività è stato quello di impiegare un solutore dedicato, MSC.Nastran SOL 200, pensato per calcoli di ottimizzazione strutturale ed appunto analisi di sensitività. Per la scelta effettuata di separare il carico in una massa concentrata che funge da variabile di progetto ed un carico di accelerazione costante il software non è in grado di calcolare la variazione del carico imposto legata alla variazione della massa concentrata, ma questa limitazione è stata agevolmente risolta sfruttando il linguaggio di programmazione interno al Nastran, DMAP. In conseguenza di questa modifica i risultati ottenuti dalla SOL 200 sono uguali a quelli dello schema alle differenze finite entro un errore medio del 9.2%, anche qui accettabile. Le diverse procedure impiegate si validano quindi l’un l’altra.
Nella regione di impatto più plausibile si devono ritrovare rapporti tra le risposte (qui, energie dei segnali di accelerazione) più vicini possibile a quelli sperimentali. A tal scopo è stato introdotto un “fattore di merito” definito dalle (3), dove con  si è indicata la sensitività, il pedice  indica il canale e gli altri due la regione assiale ed angolare. Il confronto eseguito tra rapporti di energie e sensitività è reso lecito dalla loro relazione di proporzionalità identificata dalla (1).

          (3)

 

 

 

 

 

Resta definita una matrice di fattori di merito in cui ciascun elemento corrisponde ad una posizione di impatto tra le 12 scelte, e laddove esso è minimo corrispondentemente si ha il punto di contatto più probabile, essendo massimizzata la correlazione dei rapporti tra le risposte. Sono state identificate così le 3 migliori posizioni, una per ogni fascia assiale.

Procedura per identificazione impatto: ottimizzazione
Le tre posizioni rimaste sono state poi investigate con l’ottimizzazione strutturale, volta a massimizzare la correlazione tra le energie dei segnali sperimentali e quelle numeriche attraverso la minimizzazione della funzione obiettivo definita dalla (4),


     (4)

 

 

 

rispettando dei vincoli relativi ai massimi delle accelerazioni più alte, forzati a restare sufficientemente vicini ai corrispondenti valori sperimentali. Detti vincoli sono risultati essere molto sensibili alla regione assiale (motivo per il quale tre posizioni sono state processate a valle della sensitività) e necessari per una soddisfacente correlazione dei segnali nel tempo. Anche l’ottimizzazione è stata condotta in due modi, sfruttando un ottimizzatore Matlab® che impiegasse esternamente l’MSC.Nastran SOL 112 per i calcoli FEM, e di nuovo il solutore dedicato MSC.Nastran SOL 200. I risultati per le tre posizioni sono riportati in Tabella 2. Ovviamente per soluzioni “feasible” essi sono gli stessi, e dunque restano validati. La posizione di impatto è stata così individuata approssimativamente da quella indicata come BASSA 270°, insieme con la massa concentrata ottimizzata che si traduce in un carico applicato banalmente rimoltiplicando per il carico di accelerazione fornito con la scheda ACCEL1. La soluzione può essere rifinita andando a sfruttare l’algoritmo di ottimizzazione in altre posizioni in prossimità di quella ottima attuale: si trova allora che la migliore è per una posizione intermedia indicata come MEDIO-BASSA 270°.

 

 

Dall’integrazione della seconda legge di Newton il carico ottimizzato ottenuto può anche tradursi in una quantità di moto trasferita alla parete dell’ugello durante l’impatto (Eq. (5)). Si può fare inoltre l’ipotesi di urto perfettamente elastico, poiché si sta giudicando la possibilità che l’oggetto impattante sia un frammento di gomma, contro la parete molto più rigida che forma l’ugello. In questo modo, indicando con  la velocità iniziale del componente e con  la sua massa reale, la conservazione della quantità di moto si scrive come nella (6).
  

(5)

(6)

 

 

 

A ciò si può aggiungere una stima della massa reale dell’oggetto, fornita dal committente Avio S.p.A. Questa è stata ottenuta ricorrendo alle misure di pressione in camera di combustione durante l’anomalia, in cui si è registrato un picco supposto dovuto al passaggio dell’oggetto nella gola dell’ugello. La sua entità è stata messa in relazione alla sezione dell’oggetto, mentre la sua lunghezza alla durata del picco in camera. Nota questa massa (0.278 kg), dalla (5) e dalla (6) si ha una stima della velocità iniziale dell’oggetto impattante pari a 10.8 m/s: si tratta in realtà della sua sola componente normale, poiché come detto si stanno trascurando gli effetti nell’impatto delle altre due. Tale stima resta vincolata all’assunzione di urto perfettamente elastico, dunque alla supposta natura elastomerica dell’oggetto. Si hanno così tutti gli elementi per mettere a punto una simulazione del fenomeno di impatto con lo scopo di validare quanto ottenuto.

Simulazione esplicita non lineare
La validazione è stata condotta con il solutore esplicito non lineare MSC.Nastran SOL 700, avendo a disposizione tutte le informazioni sull’oggetto impattante (massa, velocità e caratteristiche).
Per far ciò il modello originale è stato adattato alle esigenze del nuovo solutore e ciò ha comportato, tra gli altri, tre problemi principali da curare per una simulazione efficace.
Il primo riguarda il calcolo dell’incremento temporale nel solutore esplicito. Per rendere lo schema esplicito stabile esso deve suddividere il più piccolo periodo naturale nella mesh del modello. Per semplicità, l’MSC.Nastran SOL 700 impiega il criterio di Courant, basato sul tempo minimo necessario ad un’onda di stress per attraversare un elemento e sintetizzabile nella relazione (7), in cui  è la lunghezza dell’elemento,  la velocità del suono all’interno di esso e  un fattore minore di 1 che assicura la stabilità.
     

 

(7)
 

 

 

 

L’applicazione del criterio è eseguita dal solutore in maniere differenti a seconda del tipo di elemento, tuttavia un problema è stato riscontrato con il modello fornito dal committente per il giunto flessibile dell’ugello: questo è composto di un elemento di rigidezza e smorzamento concentrati a sei gradi di libertà (CBUSH), che connette le parti mobile e fissa dell’ugello. I valori di rigidezza forniti erano troppo elevati per questo scopo, fornendo un incremento temporale dell’ordine di 10−20 s, con un tempo computazionale stimato di quasi 280000 ore. Per superare il problema, il modello del giunto è stato sostituito con uno rigido: in tal modo tuttavia si è ovviamente persa la possibilità di riprodurne lo smorzamento.
Il secondo problema è dovuto al fatto che lo schema numerico impiega formule di quadratura gaussiane ad un punto per l’integrazione numerica. Ciò lo rende particolarmente veloce, ma introduce dei modi numerici degli elementi, chiamati modi di hourglass. La loro influenza viene misurata dall’energia di hourglass, che deve quindi risultare trascurabile rispetto a quelle cinetica ed interna del sistema affinché le soluzioni siano fisicamente valide. Da un lato il controllo viene fatto selezionando opportuni schemi numerici di smorzamento dei modi di hourglass già implementati in SOL 700, dall’altro anche una mesh appropriata e simile tra i corpi impattanti fornisce un contributo positivo.
Una terza interessante questione riguarda la definizione dello smorzamento. Nella procedura di soluzione del problema inverso, che sfrutta analisi modali, era stato inserito uno smorzamento percentuale dell’1% su tutti i modi, ma questo non è supportato dal solutore esplicito non lineare. Si è quindi ricorsi allo smorzamento proporzionale di Rayleigh, definito dalla (8) e che può essere correlato con il modale in una banda di interesse sfruttando la decomposizione riportata nella (9).

(8)
   
(9)

 

 

 

La banda adeguata è stata selezionata dopo una analisi FFT delle 11 risposte sperimentali, dopodiché la definizione dello smorzamento per i vari elementi è stata fornita calcolando i due coefficienti  e  da un approccio ai minimi quadrati basato sull’Eq. (9) nella banda di frequenza selezionata.

 

I risultati della validazione in termini di energia cinetica, interna e di hourglass sono riportati in Figura 5. Gli andamenti temporali confermano che l’energia di hourglass per il modello definitivo impiegato è trascurabile, validando il calcolo effettuato; si può inoltre verificare l’effetto dello smorzamento di Rayleigh essendo l’energia totale decrescente, seppure la riduzione sia contenuta poiché l’oggetto rimbalza con velocità prossima a quella iniziale, conservando gran parte dell’energia cinetica. L’andamento di tale velocità per i primi 10 ms della simulazione è riportato in Figura 6 e conferma l’importante ipotesi di urto perfettamente elastico. Questa non è ovviamente esatta (ciò significherebbe assenza di trasferimento di energia all’ugello), ma la velocità finale è solo il 22.5% inferiore a quella iniziale, introducendo un errore contenuto entro l’11.2% sulla stima di quest’ultima fatta mediante la (6).

 

 

 

 

La correlazione per il canale più rilevante (#2) è riportata in Figura 7. È evidente che i massimi sono molto prossimi, il contenuto in frequenza è leggermente spostato verso frequenze più alte per la soluzione numerica, mentre si ha un errore nel tempo di smorzamento per la presenza di due inattesi battimenti acustici. La causa sembra essere il modello rigido del giunto: infatti la bassa frequenza contenuta nei due battimenti è di circa 30 Hz, prossima a quella del primo modo di pendolo dell’ugello. Essendo assente lo smorzamento del giunto, tale modo sembra non essere quindi smorzato.
Passando alle energie, ottimizzate dalla procedura di soluzione del problema inverso, il confronto è riportato per la prima metà della simulazione (Figura 8), per filtrare il problema di smorzamento del giunto. Si riconosce una soddisfacente correlazione, in cui solo i canali 7 e 8 non sono ben in linea con i valori sperimentali, a causa del modello semplicistico di trave utilizzato per i martinetti a fronte di una struttura reale molto complessa che richiederebbe analisi dedicate.

 

 

 

Conclusioni
La correlazione numerico-sperimentale ottenuta con la procedura sviluppata risulta soddisfacente in termini delle quantità globali impiegate (energie dei segnali, massimi e tempi di smorzamento), premesse le incertezze dovute ad una struttura complessa a fronte di un modello agli elementi finiti semplificato. In risposta alla richiesta di valutazione iniziale, la procedura conferma che un problema di impatto è una causa plausibile a spiegare le accelerazioni ottenute, con una buona correlazione accoppiati a valori realistici dell’eccitazione fornita nell’impatto. La posizione di impatto ottenuta si trova sul lato sinistro dell’ugello, come anticipato dalle considerazioni preliminari sulle accelerazioni, a riprova della bontà dei risultati. Il carico di impatto è stato determinato come forza e quantità di moto trasferita all’ugello; inoltre, sulla base di una stima della massa impattante, si è tradotta l’informazione anche in termini di velocità dell’oggetto.
Molto importante per gli scopi del lavoro era confermare la natura del corpo impattante, necessaria per confermare l’anomalia come non critica. L’ipotesi che l’oggetto fosse un frammento di protezione termica, ossia una gomma, è stata impiegata sia nell’assunzione di urto perfettamente elastico sia nella modellizzazione FEM del corpo. L’urto perfettamente elastico è stato dimostrato con un soddisfacente grado di approssimazione ed anche la correlazione, giudicata soddisfacente, ha confermato l’ipotesi sul frammento di protezione termica.
Infine, il modello per analisi esplicita non lineare potrà essere eventualmente impiegato in future simulazioni di impatto a scopi di validazione strutturale da parte di Avio S.p.A., proprietaria del design, sebbene uno sviluppo del modello del giunto flessibile sia consigliabile.

« Indice del n. 67