Tecniche Numeriche per lo Studio delle Protezioni Termiche ......Tecniche Numeriche per lo Studio...

Tecniche Numeriche per lo Studio delle ProtezioniTermiche Ablative

Alessandro Turchi

Roma21 Marzo 2011

Tecniche Numeriche per lo Studio delle Protezioni Termiche Ablative - A. Turchi

Outline

1 IntroduzioneSistemi di protezione termica ablativaMateriali TPSIntroduzione al fenomeno dell’erosione negli ugelli

2 Metodologie di calcolo dell’erosione

3 Trasporto del calore nei materiali ablativiConduzione in caso di parete mobile

Soluzione stazionaria

Il codice CMAMateriali pirolizzantiProcedura risolutiva di CMATabelle termochimicheEsempio applicativo di CMA

4 Approccio CFDModello teorico/numericoCasi test

Il motore BatesZefiro 9

() 21 Marzo 2011, Roma 2 / 30


Caratteristiche delle TPS Ablative

Consumo della TPSI sistemi di protezione termica ablativa (TPS) consentono di evitare il surriscaldamento delle strutture sottostantigrazie al loro consumo⇒ sono largamente utilizzati nei veicoli di rientro e negli ugelli degli endoreattori asolido

Dimensionamento delle protezioni TPSA causa del suo consumo, il materiale della protezione termica deve essere accuratamente dimensionato.

Due principali categorie di materialiDistinguiamo due principali categorie di materiali ablativi: i pirolizzanti e i non-pirolizzanti

() 21 Marzo 2011, Roma 3 / 30


Materiali TPS

non-Pirolizzanti

Il consumo di materiale è superficiale (erosione termochimica e/o meccanica).

Pirolizzanti

Sono composti da fibre di rinforzo in una matrice di resina fenolica;

Quando riscaldata, la resina subisce una serie di reazioni chimiche che liberanoprodotti gassosi (pirolisi) trasformando il materiale restante in un residuo poroso(char);

La pressione all’interno del materiale forza questi gas di pirolisi a fluire attraverso ilchar residuo fino alla superficie del materiale;

Il restante materiale pirolizzato, viene consumato anch’esso da azioni chimichee/o meccaniche agenti sulla superficie.

() 21 Marzo 2011, Roma 4 / 30


Protezioni Termiche negli Ugelli (1/2)

Negli ugelli degli endoreattori a solido le temperature raggiunte sono molto elevate(3000− 3500 K);

A causa dell’assenza di un sistema di raffreddamento attivo, gli ugelli devono essere costruiticon materiali resistenti alle alte temperature;

I materiali utilizzati sono selezionati in base alla loro capacità di assorbimento termico(elevato calore specifico e elevata energia richiesta per la decomposizione) e buonecapacità di isolamento (bassa conducibilità termica);

Il Carbon-Carbon è il materiale più utilizzato per proteggere la zona di gola;

Una certa quantità di materiale TPS viene consumato per evitare che la struttura metallicasuperi le temperature consentite.

() 21 Marzo 2011, Roma 5 / 30


Protezioni Termiche negli Ugelli (2/2)

La conoscenza della variazione del profilo interno dell’ugello èfondamentale per la corretta predizione delle prestazioni del motore

⇓

A causa del bilancio tra la portata generata dal grano equella smaltita dalla gola, la crescita dell’area di gola porta alla

diminuzione della pressione in camera

Perdite nei boosters dello Space Shuttle ∗

Impulso Specifico Teorico (nel vuoto) 278.1 sImpulso Specifico Effettivo (nel vuoto) 268.2 sPerdite (calcolate): (9.9 s totali)

Flusso bidimensionale e bifase (incluse perdite di divergenza) 7.4 sErosione della gola (riduzione del rapporto d’area) 0.9 sBStrato limite (perdite per attrito) 0.7 sSommergenza 0.7 sEtc. . . .

∗ G. P. Sutton & O. Biblarz, Rocket Propulsion Elements, 7th ed.

() 21 Marzo 2011, Roma 6 / 30


Introduzione al Fenomeno dell’Erosione negli Ugelli

Le reazioni chimiche eterogenee (superficiali) sono le principali responsabili delconsumo del materiale TPS

⇓Numerosi studi sperimentali hanno riscontrato un legame tra le specie ossidanti

presenti nei prodotti di combustione (H2O, CO2 and OH),e la velocitàcon cui viene consumato il materiale TPS

⇓

La velocità con cui viene consumato il materiale di parete è influenzata dalla cineticachimica delle suddette reazioni e dalla velocità di diffusione dei reagenti

attraverso lo strato limite fino alla superficie ablativa⇓

Nel caso di elevata velocità di reazione (alte temperature), la velocità di erosione èlegata prevalentemente ai meccanismi di diffusione delle specie ossidanti

() 21 Marzo 2011, Roma 7 / 30




⇓

Numerosi studi sperimentali hanno riscontrato un legame tra le specie ossidantipresenti nei prodotti di combustione (H2O, CO2 and OH),e la velocità

con cui viene consumato il materiale TPS

⇓La velocità con cui viene consumato il materiale di parete è influenzata dalla cinetica

chimica delle suddette reazioni e dalla velocità di diffusione dei reagentiattraverso lo strato limite fino alla superficie ablativa

⇓


() 21 Marzo 2011, Roma 7 / 30




⇓


con cui viene consumato il materiale TPS⇓


attraverso lo strato limite fino alla superficie ablativa

⇓Nel caso di elevata velocità di reazione (alte temperature), la velocità di erosione è

legata prevalentemente ai meccanismi di diffusione delle specie ossidanti

() 21 Marzo 2011, Roma 7 / 30




⇓


con cui viene consumato il materiale TPS⇓


attraverso lo strato limite fino alla superficie ablativa⇓


() 21 Marzo 2011, Roma 7 / 30


Metodologie di Calcolo dell’Erosione

Approccio a Coefficienti

A parete, il flusso di massa è legato al flusso termico convettivo attraverso unarelazione semiempirica

Utilizzo delle tabelle termochimiche (per il bilancio energetico di parete).

Approccio CFD

Risoluzione completa del flusso (2D assialsimmetrico) all’interno dell’ugello

Condizione al contorno ablativa

Modello di equilibrio o cinetico per le reazioni di parete

Ipotesi di “ablazione stazionaria” per il flusso conduttivo all’interno della parete

() 21 Marzo 2011, Roma 8 / 30


Conduzione del Calore nei Materiali Ablativi

Equazione di Fourier “classica”:

∂

∂t(ρhA)y = ∂

∂y

(kA∂T∂y

)t

Equazione di Fourier in coordinate mobili:

∂

∂t(ρhA)x = ∂

∂x

(kA∂T∂x

)t+ s

∂

∂x(ρh)t

dove la coordinata x è solidale con la parete mobile (indica una posizione ad una certa profondità nel

materiale) e la coordinata y indica la posizione nel riferimento fisso.() 21 Marzo 2011, Roma 9 / 30


Soluzione Stazionaria dell’Equazione del Calore (1/2)

distanza dalla superficie [mm]

Tem

pera

tura

[K]

0 2 4 6 8 10 12

1000

2000

3000

4000

0.5 s1 s5 s10 s50 s70 s

Parete fissa

t


Tem

pera

tura

[K]

0 20 40 60 80 100 120

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s90 s100 s

0.125 mm/s, k=30 W/m-K, cp= 2000 J/kg-K

t

Nel caso di parete mobile la soluzione stazionaria è tale per cuiil profilo di temperatura nel materiale rimane costante nel tempo

T − T0

Ts − T0= e−

sxα

() 21 Marzo 2011, Roma 10 / 30


Soluzione Stazionaria dell’Equazione del Calore (2/2)

tempo [s]

Tem

pera

tura

[K]

0 20 40 60 80 100

1000

2000

3000

4000

0 mm1 mm2 mm5 mm10 mm20 mm40 mm

Parete fissa

tempo [s]

Tem

pera

tura

[K]

-20 0 20 40 60 80 100 1200

500

1000

1500

2000

2500

3000

3500

4000

4500

0 mm1 mm2 mm5 mm10 mm20 mm40 mm

0.125 mm/s

() 21 Marzo 2011, Roma 11 / 30


Influenza della Velocità di Regressione (1/2)


Tem

per

atur

a[K

]

0 1 2 3 4

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s90 s

t

0.125 mm/s


Tem

per

atur

a[K

]0 1 2 3 4

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s90 s

t

0.25 mm/s

() 21 Marzo 2011, Roma 12 / 30


Influenza della Velocità di Regressione (2/2)

La velocità con cui si muove la superficie influenzail raggiungimento dello stato stazionario


Tem

per

atur

a[K

]

0 0.5 1 1.5 2

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s

t

0.5 mm/s


Tem

per

atur

a[K

]

0 0.05 0.1

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s90 s

t

1 mm/s

() 21 Marzo 2011, Roma 13 / 30


Influenza della Diffusività Termica

E’ la diffusività termica (k/(ρcp)) ad influenzareil raggiungimento dello stato stazionario


Tem

pera

tura

[K]

0 20 40 60 80 100 120

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s90 s100 s

0.125 mm/s, k=30 W/m-K, cp= 2000 J/kg-K

t


Tem

pera

tura

[K]

0 20 40 60 80 100 120

500

1000

1500

2000

2500

3000

3500

4000

0.5 s1 s2 s5 s40 s90 s100 s

0.125 mm/s, k=60 W/m-K, cp= 4000 J/kg-K

t

() 21 Marzo 2011, Roma 14 / 30


Trasporto del Calore nei Materiali Pirolizzanti

ρcp∂T∂t

∣∣∣x

= 1A∂

∂x

(kA∂T∂x

)t+ (hg − h)∂ρ

∂t

∣∣∣y

+ sρcp∂T∂x

∣∣∣t+

mg

A∂hg

∂x

∣∣∣t

I termini aggiuntivi sono legati ai contributi energetici derivanti dalla decomposizione delmateriale (variazione di densità e generazione di un gas di pirolisi).

() 21 Marzo 2011, Roma 15 / 30


Charring Material Thermal Response and Ablation Program (CMA) (1/2)

La procedura seguita da CMA può essere sintetizzata come segue:

1) Calcolo decomposizione (se presente)

∂ρi∂t = −Bie

( EaiRT

)ρ0i

(ρi−ρriρ0i

)ψi

2) Calcolo conduzione nel sistema mobile (solidale con la parete)ρcp

∂T∂t

∣∣∣x

= 1A∂∂x

(kA ∂T

∂x

)t+ (hg − h) ∂ρ

∂t

∣∣∣y

+ sρcp∂T∂x

∣∣∣t+ mg

A∂hg∂x

∣∣∣t

3) Verifica del bilancio energetico superficiale

ρeueCh(hr − hw)e︸︷︷︸qsen

+ ρeueCm

[ N∑i=1

(Z∗ie − Z∗iw )hTwi − B′hw

]+ mchc + mghg︸︷︷︸

qchem

= qradin− qradout + qcond

() 21 Marzo 2011, Roma 16 / 30



La procedura seguita da CMA può essere sintetizzata come segue:

1) Calcolo decomposizione (se presente)

∂ρi∂t = −Bie

( EaiRT

)ρ0i

(ρi−ρriρ0i

)ψi

2) Calcolo conduzione nel sistema mobile (solidale con la parete)

ρcp∂T∂t

∣∣∣x

= 1A∂∂x

(kA ∂T

∂x

)t+ (hg − h) ∂ρ

∂t

∣∣∣y

+ sρcp∂T∂x

∣∣∣t+ mg

A∂hg∂x

∣∣∣t

3) Verifica del bilancio energetico superficiale

ρeueCh(hr − hw)e︸︷︷︸qsen

+ ρeueCm

[ N∑i=1

(Z∗ie − Z∗iw )hTwi − B′hw

]+ mchc + mghg︸︷︷︸

qchem

= qradin− qradout + qcond

() 21 Marzo 2011, Roma 16 / 30



Il CMA permette di simulare la recessione della parete solida: velocità istantanea dierosione, flussi termici sulla parete e temperatura all’interno del materiale solido. Ilconto può essere effettuato con condizioni al contorno (gas–parete) variabili.

⇓Gli input richiesti all’utente sono i seguenti:

Coefficiente di Scambio Termico: ρeueCh

Entalpia di Recupero: hr = h + r(u2e/2)

Flusso Radiativo: qradin

Tabelle Termochimiche

Cm/Ch

Proprietà dei Materiali: conducibilità termica, densità, calore specifico, etc.

() 21 Marzo 2011, Roma 17 / 30



Le tabelle forniscono gli input necessari per risolvere il bilancio energetico di parete

Variabili indipendenti

Parametro di Blowing (char): B′c = mc/(ρeueCm);

Parametro di Blowing (gas di pirolisi):B′g = mg/(ρeueCm);

Variabili dipendenti

Entalpie richieste nel bilancio energeticosuperficiale

Temperatura di parete

() 21 Marzo 2011, Roma 18 / 30



Le tabelle forniscono gli input necessari per risolvere il bilancio energetico di parete

La portata di materiale eroso dipendefortemente da Cm/Ch:

Cm

Ch= Le

23ww�

m =B′(ρeueCm)

() 21 Marzo 2011, Roma 18 / 30


Esempio: Motore a Propellente Solido (1/2))

1 Motore di lunga durata (>100s);

2 Propellente alluminizzato;

3 Temperatura in camera dell’ordine dei 3500 K;

4 Pressione variabile nel tempo;

5 Parete in Carbon–Carbon.

() 21 Marzo 2011, Roma 19 / 30


Esempio Applicativo: Motore a Propellente Solido (2/2)

∗ A. Turchi, Numerical Techniques for Nozzles Throat Erosion,Final Thesis of Master Course in “Space Transportation Systems”

() 21 Marzo 2011, Roma 20 / 30


Esempio Applicativo: Pressione Media vs Pressione Variabile

La dipendenza dell’erosione dalla pressione risulta pressoché lineare

() 21 Marzo 2011, Roma 21 / 30


Esempio Applicativo: Influenza del Coefficiente di Scambio Termico

L’influenza del Ch sulla soluzione finale è pressoché lineare

Variando il valore di tale coefficiente siottiene infatti:

Tempo [s] Eroso (Ch -20%) Eroso (Ch +20%)

10.00 79.7% 120.3%60.00 79.9% 120.2%120.00 79.8% 120.2%

() 21 Marzo 2011, Roma 22 / 30


Esempio Applicativo: Influenza del Coefficiente di Scambio di Massa

L’influenza del Cm/Ch sulla soluzione finale è pressoché lineare

Variando il valore di tale rapporto si ottieneinfatti:

Tempo [s] Eroso(

CmCh

-20%)

Eroso(

CmCh

+20%)

10.00 80.4% 119.5%60.00 80.3% 119.6%

120.00 80.3% 119.5%

() 21 Marzo 2011, Roma 23 / 30


Esempio Applicativo: Influenza di Ulteriori Parametri

Parametro Influenza su Eroso Influenza su Twall Influenza su Tdepth

Tabelle Termochimiche Alta Alta AltaGas di Pirolisi Bassa Bassa BassaB.C. Esterna Bassa Bassa Bassa

Back–up Pirolizzante Bassa Bassa MediaFlusso Radiativo Bassa Bassa Bassaww�

Il parametro più sensibile è senza dubbio il Cm/Ch

() 21 Marzo 2011, Roma 24 / 30


CFD Approach to Throat Erosion Prediction

Why:

To advance the fundamental understanding of erosion processes

To support the experimental tests during the development

How:

Erosion is handled using an accurate “ablative” boundary condition

A simplified approach is used to manage the conduction inside the solid material

The following assumptions are made in the analysis:

Flow is axisymmetric and steady;

Radiation heat transfer is negligible;

Negligible effect of gas-phase reactions on erosion (due to small concentration of oxygen).

Code Specifications*:

Navier-Stokes 2-D axisymmetric solver

Finite difference method (λ scheme)

One-equation Spalart-Allmaras turbulence model

Le = cost diffusion model

Thermodynamic and transport properties described by curve fits (Gordon & McBride, 1994)∗ D. Bianchi, F. Nasuti, E. Martelli, Coupled Analysis of Flow and Surface Ablation in Carbon-Carbon Rocket Nozzles, Journal of Spacecraft and Rockets, 2009

() 21 Marzo 2011, Roma 25 / 30


Surface Balance Equations

Surface mass balance (SMB):

ρDi∂Yi

∂y

∣∣∣∣w

+ mYs,i = (ρv)Yw,i i = 1,Nc

where Ys,i = mi/m is the mass of species i produced or consumed in the ablation process per unit mass

of ablated material

Surface energy balance (SEB):

k∂T∂y

∣∣∣∣∣w︸︷︷︸

convective

+Nc∑i

hiρDi∂Yi

∂y

∣∣∣∣∣w

− m(hw − hs)︸︷︷︸chemical

+ ©©qrad︸︷︷︸radiative

= ©©©εσT4︸︷︷︸re−radiative

+ m(hs − hs0 )︸︷︷︸solid conduction

where the conduction term is represented by a closed expression available at steady state() 21 Marzo 2011, Roma 26 / 30


Bates Motor Tests

The nozzle geometry reproduces the BATES∗ experimental motor

Five different propellant compositions have been considered for wide variations of aluminum content:

CO CO2 HCl H2 H2O N2 Al2O3 pc (bar) Tc (K) Al

0.175 0.04 0.24 0.02 0.145 0.10 0.28 69 3580 15%0.18 0.025 0.23 0.02 0.105 0.10 0.34 69 3655 18%0.20 0.015 0.195 0.02 0.07 0.10 0.40 69 3715 21%0.20 0.005 0.190 0.02 0.045 0.10 0.44 69 3750 24%0.20 0.005 0.190 0.02 0.025 0.10 0.46 69 3745 27%

∗ R. L. Geisler The History of the BATES Motors at the Air Force Rocket Propulsion Laboratories, AIAA Paper 1998

() 21 Marzo 2011, Roma 27 / 30


Bates Motor Tests

Calculated erosion rates∗ are compared with the experimental data by Geisler

The finite-rate model improves the predictions when the kinetic-limited regime is approached⇒ this alsoconfirms a good prediction of surface temperature

∗ D. Bianchi et al., Coupled Analysis of Flow and Surface Ablation in Carbon-Carbon Rocket Nozzles, Journal of Spacecraft and Rockets 2009 vol.46 no.3

∗ D. Bianchi et al., Thermochemical Erosion Analysis for Carbon-Carbon Rocket Nozzles, accepted for publication on Journal of Propulsion and Power() 21 Marzo 2011, Roma 28 / 30


Zefiro 9 motor

Zefiro 9 SRM

Zefiro 9 is the third stage of the small European launcher VEGA,realized by AVIO Group S.p.a., whose characteristics area:

Grain Mass [tons] ≈ 10Grain Shape Finocyl/Aft StarDiameter [m] ≈ 3Length [m] ≈ 3.8

Burn Time [s] ≈ 120Vacuum Specific Impulse [s] ≈ 295

Max. Vacuum Thrust [kN] ≈ 330MEOP [bar] ≈ 85

Nozzle Expansion Ratio ≈ 61

aThe European Small Launcher-BR257, April 2007, Vega Programme Dept., ESADirectorate of Launchers

() 21 Marzo 2011, Roma 29 / 30


Zefiro 9 motor Erosion Results

Z9 Motor: erosion analysis results (Equilibrium Model)

Final Shape (Error %)Section CFD Equilibrium Expt. data∗

1 +2.5% ±1.6%2 +2.9% ±1.8%3 +2.6% ±1.7%4 +3.7% ±1.4%5 +4.2% ±1.4%6 +4.2% ±1.2%7 +3.8% ±1.1%8 +3.8% ±1.1%9 +0.8% ±0.8%10 −0.3% ±0.7%

The equilibrium ablation model slightly over-estimate the final eroded profile

∗ courtesy of Avio Group

() 21 Marzo 2011, Roma 30 / 30


Zefiro 9 motor Erosion Results

Z9 Motor: erosion analysis results (Finite-Rate Model)

Final Shape (Error %)Section CFD Finite-Rate Expt. data∗

1 −0.5% ±1.6%2 −0.8% ±1.8%3 −1.4% ±1.7%4 −1.0% ±1.4%5 −0.8% ±1.4%6 −0.9% ±1.2%7 −1.2% ±1.1%8 −1.2% ±1.1%9 −3.1% ±0.8%10 −3.1% ±0.7%

The finite-rate ablation model improves the agreement with the experimental data

∗ courtesy of Avio Group

() 21 Marzo 2011, Roma 30 / 30

Tecniche Numeriche per lo Studio delle Protezioni Termiche ......Tecniche Numeriche per lo Studio...

Documents

Transcript of Tecniche Numeriche per lo Studio delle Protezioni Termiche ......Tecniche Numeriche per lo Studio...