Metodi Numerici v2.2

Indice 1 Generalit` a sui metodi numerici 1.1 Analisi degli errori . . . . . . . . . . . . . . . . . . . . . . 1.2 Condizionamento di un problema numerico . . . . . . . . 1.3 Stabilit` a di un algoritmo . . . . . . . . . . . . . . . . . . . 1.4 Criteri di convergenza e troncamento di processi iterativi .

. . . .

. . . .

3 3 5 9 12

2 Metodi numerici 16 2.1 Metodi di risoluzione di equazioni lineari con sorgente . . . . 16 2.1.1 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . 16 2.1.2 Metodo di Gauss-Seidel . . . . . . . . . . . . . . . . . 23 2.1.3 Metodo di Jacobi vs Metodo di Gauss . . . . . . . . . 24 2.1.4 Metodo SOR . . . . . . . . . . . . . . . . . . . . . . . 25 2.1.5 Shift tra metodo di Jacobi e metodo SOR . . . . . . . 28 2.1.6 Ordinamento dei nodi di flusso . . . . . . . . . . . . . 29 2.2 Risoluzione di sistemi di equazioni non lineari . . . . . . . . . 30 2.2.1 Metodo di Newton per la risoluzione di sistemi non lineari . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Metodi per problemi agli autovalori . . . . . . . . . . . . . . . 31 2.3.1 Scomposizione QR di una matrice . . . . . . . . . . . 31 2.3.2 metodo QR . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.3 Metodo delle potenze . . . . . . . . . . . . . . . . . . 33 2.4 Propriet` a spettrali delle matrici di interazione . . . . . . . . . 38 2.5 Approssimazione di dati sperimentali . . . . . . . . . . . . . . 39 2.5.1 Funzioni spline . . . . . . . . . . . . . . . . . . . . . . 39 2.6 Integrazione numerica . . . . . . . . . . . . . . . . . . . . . . 42 2.7 Equazioni integrali . . . . . . . . . . . . . . . . . . . . . . . . 50 3 Metodi numerici per l’ingegneria nucleare 52 3.1 Introduzione alla tecnologia ed alla fisica nucleare . . . . . . . 52 3.1.1 Il concetto di criticità e le sue implicazioni . . . . . . . 55 3.1.2 Fenomeni non stazionari . . . . . . . . . . . . . . . . . 56 3.1.3 L’importanza neutronica ed i metodi per lo studio di problemi pseudo-stazionari . . . . . . . . . . . . . . . 57

1

3.2 3.3

3.4

3.5 3.6

Il consumo del combustibile nucleare . . . . . . . . . . . . . . Differenze finite fini . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Diffusione monodimensionale monogruppo . . . . . . . 3.3.2 Problemi con sorgente . . . . . . . . . . . . . . . . . . 3.3.3 Non omogeneità dello spazio . . . . . . . . . . . . . . 3.3.4 Diffusione bidimensionale monogruppo . . . . . . . . . 3.3.5 Diffusione multigruppo . . . . . . . . . . . . . . . . . . 3.3.6 Accelerazione della convergenza delle iterazioni esterne Metodi Coarse Mesh . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Metodi Quabox . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Metodi Cubbox . . . . . . . . . . . . . . . . . . . . . . Elementi di cinetica neutronica . . . . . . . . . . . . . . . . . Elementi di Trasporto neutronico . . . . . . . . . . . . . . . . 3.6.1 Formulazione integrale . . . . . . . . . . . . . . . . . .

A Aggiustamento di sezioni d’urto A.1 Introduzione alle basi di dati . . . . . . . . A.2 Elementi di statistica . . . . . . . . . . . . . A.3 Aggiornamento di dati differenziali a partire A.4 Teoria delle perturbazioni . . . . . . . . . .

2

. . . . . . . . . . . . . . . . . . . . da dati integrali . . . . . . . . . .

60 61 61 65 66 69 72 85 88 89 92 93 101 102 106 106 108 112 116

Capitolo 1

Generalit` a sui metodi numerici 1.1

Analisi degli errori

L’analisi numerica di un problema fisico porta necessariamente alla produzione di un certo numero di errori, dei quali è necessario conoscere l’entità al fine di ridurla al minimo Esistono varie tipologie di errore: Errori di idealizzazione : Derivano dalla decisione di trascurare determinati elementi del problema fisico di partenza per semplicità Errori di matematizzazione : Derivano dalla semplificazione delle equazioni allo scopo di ottenere un modello matematico pi` u facile da trattare. Questa perdita di informazioni è di difficile recupero Errori numerici : Risiedono intrinsecamente nella risoluzione numerica di un problema. Sono costituiti ad esempio dalla discretizzazione temporale e spaziale del problema Errori sui dati di input : Derivano come facilmente comprensibile da eventuali incertezze o imprecisioni sui dati di input Errori di calcolo numerico : Derivano da tutte le approssimazioni effettuate all’interno della risoluzione del problema numerico, come il troncamento di processi iterativi Errori di misura : Derivano dal fatto che nel momento in cui si vanno a confrontare i dati numerici con quelli sperimentali questi ultimi saranno inevitabilmente soggetti ad un errore di misura. In questa fase è opportuno ricordare che grande attenzione deve essere posta nell’utilizzo degli strumenti, in particolare evitando che essi vengano utilizzati al di fuori del loro range

3

Errori di arrotondamento : Derivano dal fatto che un calcolatore può memorizzare solo un numero finito di cifre decimali. Viene qui approfondito il discorso legato agli errori di arrotondamento. L’insieme detto dei numeri di macchina è di fatto un sottoinsieme dell’insieme dei numeri reali. All’interno di un calcolatore i numeri sono immagazzinati in virgola mobile, ovvero sottoforma di un modulo e di una potenza in base 10 1 che diventa virgola mobile normalizzata ove il modulo della mantissa è sempre 0.1 < |a| < 1. Definito t il numero di cifre decimali della mantissa, si avrà ovviamente che le dimensioni dell’insieme dei numeri di macchina aumentano con l’aumentare di t. Si ha banalmente che per t → ∞ =⇒ A → R, ovvero se t tende ad infinito non commetto errori di arrotondamento2 . Vediamo meglio come è possibile definire l’errore di arrotondamento. Si definisce u l’operatore somma per il calcolatore, dato da x u y = (x + y)(1 + )

con|| ≤ eps

Si ricordi per` o che un calcolatore opera solo con numeri di macchina. Si definisce dunque che se x non appartiene a tale insieme, esisterà un suo arrotondamento opportuno che lo porterà a diventare un numero di macchina In generale avremo che eps ∼ 10−22 . Apparentemente è una quantità molto piccola, ma se le operazioni sono associate a fattori di amplificazione dell’errore derivanti dal cattivo condizionamento del problema può risultare rilevante la propagazione di tale errore. Si supponga ora y =a+b+c

con

y˜ = f l((a + b) + c)

ove l’operazione “fl” indica operazione reale eseguita in virgola mobile. Chiamando η = f l(a + b) risultato in virgola mobile della somma tra parentesi, esso varr` a η = (a + b)(1 + 1 ) da cui y˜ = f l(η + c) = = (η + c)(1 + 2 ) = = (a + b + c)[1 +

a+b 1 (1 + 2 ) + 2 ] a+b+c

Si definisce ora y tale per cui y˜ = (a + b + c)(1 + y ) 1

ad esempio il numero 439 diventa 0.439 · 102 Per maggiore completezza è necessario precisare che è necessario che anche il numero di cifre per l’esponente di 10 deve tendere ad infinito affinchè l’insieme dei numeri di macchina tenda a quello dei numeri reali 2

4

Si pu` o ragionevolmente pensare di trascurare il termine contenente la moltiplicazione di due infinitesimi, in modo da ottenere y =

a+b 1 + 2 a+b+c

Si vede dunque che i due errori sono pesati con pesi diversi. In particolare, nel caso in esame, l’errore sulla seconda operazione pesa pi` u di quello sulla prima. Due algoritmi diversi accumuleranno errori di arrotondamento in modo diverso. Tuttavia questa considerazione perde di validità ove il valore di t tenda ad infinito. Questa condizione è definita di computer ideale ovvero algoritmo indipendente. In un computer reale invece, anche avendo in input solo numeri di macchina, sono in genere sufficienti pochi passaggi per rendere nuovamente necessaria una approssimazione, dando luogo ad errori di arrotondamento. La quantificazione di tali errori è data da eps = 5 · 10−t , ove tale valore è di fatto il massimo errore di arrotondamento che posso avere per ogni calcolo, ` del computer. ed è detto sensibilita Va ricordato che un problema numerico è disgiunto dalla strategia che viene utilizzata per risolverlo. In sostanza il problema numerico e l’algoritmo che verr` a ad esso applicato per ricavare la soluzione sono slegati l’uno dall’altro.

1.2

Condizionamento di un problema numerico

Elemento importante da studiare è il condizionamento di un problema numerico, che di fatto esprime la sua suscettibilità ad oscillazioni ed errori sui dati di input. Un problema mal condizionato non può essere risolto a cuor leggero, in quanto non si può essere sicuri dell’affidabilità dei risultati ottenuti. Questo problema risulta di particolare interesse nei problemi nucleari soprattutto per via di alcune grandezze sperimentali, come le sezioni d’urto e le lunghezze di diffusione, il cui valore è facilmente soggetto ad incertezze anche importanti. La maggior parte dei problemi di condizionamento è legata alla sensibilit` a del computer. Essa è definita in maniera tale per cui a + b = a se b < , ove è proprio la sensibilità del computer. Si prenda una generica funzione:   ϕ1 (x1 ....xn )   .. ϕ(x) =   . ϕn (x1 ....xn ) Si chiami ora x ˜ il valore approssimato, affetto da incertezza, tale che ∆x = x| |x − x ˜| e definisco xi = |x−˜ la perturbazione o incertezza relativa. x 5

Ovviamente nella realt` a il valore di tale perturbazione relativa non è noto, in quanto x stesso non è noto. In generale tuttavia, quantomeno sui dati in ingresso, si pu` o effettivamente quantificare tanto x quanto x ˜ Il risultato desiderato è una quantificazione del condizionamento del problema, vale a dire una relazione del tipo ∆y = f (∆x) Un esempio per procedere è utilizzare uno sviluppo in serie di Taylor troncato al prim’ordine. ∆yi = |˜ yi − yi | = ϕi (˜ x ) − ϕi x =

n X

(˜ xj − xj )

j=1

∂ϕi + [...] ∂xj

Il poter supporre che le entit` a |x − x ˜| siano piccole permette di approssimare la serie al prim’ordine. Evidentemente questo non è sempre vero: ove le siano dell’ordine del 30% del valore nominale, questa approssimazione non è pi` u valida. Il vettore ∆y dato dalle singole componenti delle perturbazioni in uscita sar` a dunque dato da:  ∂ϕ1  ∂ϕ1   · · · ∂x ∆x1 ∂x1 n  ..   ..  .. ∆y =  ... . .  .  ∂ϕm ∂x1

∂ϕm ∂xn

···

∆xn

Che, in forma matriciale, diventa: ∆y = Dϕ(x) · ∆x Ove Dϕ(x) rappresenta la jacobiana della trasformazione calcolata in x. Si ricorda a questo proposito che la jacobiana, calcolata in un punto preciso, è una matrice di numeri a tutti gli effetti, ed in questo caso costituisce il fattore di amplificazione della perturbazione assoluta in ingresso e dipende solo dalla definizione del problema numerico, e non dall’algoritmo scelto per la risoluzione. In particolare l’elemento {ai,j } rappresenta il fattore di amplificazione con cui la componente i-esima del vettore in uscita viene influenzata dalla perturbazione della j-sima componente del vettore di input Dunque, come gi` a detto in precedenza e qui dimostrato, il condizionamento di un problema è algoritmo indipendente. Effettuando discorsi analoghi sugli errori relativi si ottiene il fattore di amplificazione delle perturbazioni relative in ingresso: yi =

n X

xj (

i=1

∂ϕi xj ) ∂xj ϕi

Il condizionamento del problema tuttavia non è uguale per tutti i valori di input, ma potrebbe essere buono rispetto ad alcuni e meno rispetto ad 6

altri. Poichè il fattore di amplificazione non è un numero ma una matrice, utilizzando un tipo di norma piuttosto che un’altra si potranno evidenziare in modo pi` u o meno marcato eventuali picchi di condizionamento. Si noti che, seppure la libert` a di scelta della norma da utilizzare è ampia, l’operatore “norma” di una matrice deve rispettare determinate proprietà che ne permettono il confronto con vettori ed altre matrici. Un esempio classico di norma qP 2 è kAk = ai,j , ma una possibile alternativa è anche kAk = max(|ai,j |). Si prenda un esempio di problema: ϕ(x, y) = x − y, da cui avremo x−y =

y x x − y x−y x−y

Da questo si vede che se (x − y) tende a zero i fattori di amplificazione aumentano enormemente. Se il problema è lineare, esso è scrivibile come Ax = b Lo studio del condizionamento di un problema lineare sarà esprimibile dalla soluzione del problema A(x + δx) = (b + δb), ove δb dipende dalla sensibilità dello strumento. La matrice A invece è una caratteristica del problema numerico e non varia con il termine sorgente. Ovviamente indipendentemente dal valore della perturbazione, il problema è risolvibile (in questo caso, con problema linere e non omogeneo) solo se det(A) 6= 0 kδxk kδbk ≤ kAkkA−1 k · kxk kbk Infatti si avr` a che: Ax + Aδx = b + δb Che, ricordando la dicitura iniziale del problema (che rimane ovviamente valida anche in seguito alla perturbazione) diventa: Aδx = δb

→

δx = A−1 δb

Passando al confronto tra norme si ottiene: kδxk = kA−1 δbk ≤ kA−1 kkδbk Ricordando che vale anche che: kbk = kAx ≤ kAkkbk Si ottiene che, dividendo ambo i membri della disequazione precedente per la norma di b, si ottiene: kdeltaxk kδbk ≤ kA−1 k kAkkxk kbk 7

Dalla quale si deriva facilmente quanto si voleva dimostrare Il valore kAkkA−1 k è definito fattore di amplificazione, e dipende da A. Pi` u tale fattore è alto, peggiore sarà il condizionamento del problema. Si tenga ben presente il fatto che finora non è mai stato fatto alcun cenno ad algoritmi o a metodi di risoluzione del problema. Il condizionamento di un problema è infatti algoritmo indipendente, ed è invece strettamente legato al problema fisico di partenza. In tutti questi passaggi si noti che gli errori considerati non sono necessariamente quelli sui dati in ingresso. Il vettore x di input può essere contenuto in una qualsiasi fase del problema, e dunque l’errore ad esso associato potrebbe essere, ad esempio, l’errore di arrotondamento derivante dall’approssimazione del risultato del calcolo precedente. Questo insieme di considerazioni prende il nome di backward analysis ed ha la funzione di validare problemi ed algoritmi in funzione del condizionamento del problema. Il punto chiave del problema risiede ove un eventuale errore di calcolo porti ad ottnere comunque risultati plausibili, il che rende difficile accorgersi della presenza di un errore dalla semplice analisi dei risultati. Prendiamo un esempio molto comune: ci si trova davanti ad un polinomio di grado n del quale si vogliono identificare gli zeri. Questo è un problema spesso non lineare, il cui condizionamento raramente è buono, e tanto peggiore quanto pi` u alto è il grado n del polinomio, anche se non tutti gli zeri sono condizionati allo stesso modo. In generale infatti avremo che un generico zero µi di un polinomio è mal condizionato se: ∂Pn →0 ∂µ µi In generale ove ci troviamo di fronte ad un problema mal condizionato non è possibile risolverlo cos`ı com’è tramite algoritmi. Si rivela in questi casi necessario andare a variare ci` o che è stato fatto a monte, cercando ad esempio di ridurre gli errori iniziali (idealizzazione, matematizzazione, ecc.) o a riconsiderare il problema dall’inizio. Tuttavia nessuna di queste opzioni fornisce garanzie sul risultato finale. Spesso accade dunque che un problema del tipo Pn (x) = 0 venga trasformato in un problema Ax = λx, ove i due non hanno niente in comune se non la soluzione.3 3

Si noti come tutti i discorsi fatti fino ad ora per gli errori possano essere applicati in maniera perfettamente analoga a delle perturbazioni “volontarie”, come potrebbe essere ad esempio il voler testare come il punto di funzionamento di un sistema è influenzato dal variare dei parametri operativi

8

1.3

Stabilit` a di un algoritmo

Si prenda ora un passaggio qualunque del processo risolutivo di un problema numerico. Ogni singola operazione di un algoritmo si presenterà come una specie di micro-problema numerico, in cui in ingresso si ha il risultato di tutte le operazioni che stanno a monte che sarà, evidentemente, soggetto ad una incertezza legata agli arrotondamenti precedenti. Sarà dunque necessario anche lo studio del condizionamento di tale micro-problema, definito condizionamento locale. In questa visione dunque un algoritmo non è altro che una serie di micro-problemi numerici, che vanno studiati localmente dal punto di vista del condizionamento. Un esempio banale di questo concetto è l’operazione x−y. Tale problema è fortemente mal condizionato. Supposto che in una qualsiasi operazione {x1 , x2 , ..., xn } sia il vettore dei dati in ingresso ed {y1 , y2 , ..., ym } sia il vettore di uscita, si può definire y = ϕ(x) il problema numerico, che sarà caratterizzato da: ϕ : D → Rm

D ⊆ Rn

y ∈ Rm

x ∈ Rn

yi = ϕi (x1 , ....xn ) i=1:m

Cosa è invece l’algoritmo? Si avrà che ϕ(i) : D(i) → D(i+1)

con

i=0:r

e

D(j) ⊆ Rnj

ϕ = ϕ(r) ϕ(r−1) · · · ϕ(1) ϕ(0) Dunque se il dato di ingresso è x0 , si avrà che il dato in uscita sarà ϕ(0) x0 = x1 7−→ ϕ(1) x1 = x2 ... 7−→ ϕ(r) xr = xr+1 Per risolvere un problema numerico si possono applicare molti algoritmi diversi. Al di l` a delle questioni che vedremo in seguito sul costo e la velocità di convergenza di un algoritmo, è importante a priori assicurarsi che esso sia ben fatto, dall’inglese robust. Per verificare la stabilit` a di un algoritmo si tengono in conto gli arrotondamenti in virgola mobile ed il condizionamento del problema che si vuole andare a risolvere. Si supponga di avere un algoritmo composto da r funzioni ϕ(i) , ciascuna delle quali sia derivabile e continua nel suo dominio. Si definisce inoltre una nuova funzione ψ (i) come la trasformazione residua, data dalla composizione di tutte le ϕ(j) con i ≤ j ≤ r. Ne deriva che ψ (0) ≡ ϕ Sapendo che vale D(f ◦ g)x = D(f )x D(g)x

9

è possibile definire le jacobiane di ϕ e ψ, ed in particolare si può dire che: Dϕ(x) = Dϕ(r) (x(r) )Dϕ(r−1) (x(r−1) ) · · · Dϕ(0) (x(0) ) Dψ (i) (x(i) ) = Dϕ(r) (x(r) )Dϕ(r−1) (x(r−1) ) · · · Dϕ(i) (x(i) ) Si ricorda che, nel momento in cui si parla di un algoritmo, si parla di un computer reale , il che porta a dire che il numero di cifre della mantissa è finito. Si avr` a dunque a che fare con grandezze approssimate del tipo: x ˜(i+1) = f l(ϕ(i) (˜ x(i) )) Quello che ora è necessario calcolare è il valore di ∆xi ad un i-esimo passaggio di un algoritmo, che sar` a dato dal valore approssimato ottenuto meno il valore esatto: ∆x(i+1) = [f l(ϕ(i) x ˜(i) ) − ϕ(i) x(i) ] Sommando e sottraendo la quantità [ϕ(i) x ˜(i) ] si ottiene: ∆x(i+1) = [f l(ϕ(i) (˜ x(i) )) − ϕ(i) (˜ x(i) )] + [ϕ(i) (˜ x(i) ) − ϕ(i) (x(i) )] Il secondo termine tra parentesi quadre rappresenta il condizionamento locale del problema ϕ(i) , e sar` a dunque uguale a [ϕ(i) (˜ x(i) ) − ϕ(i) (x(i) )] = Dϕ(i) (x(i) ) · ∆x(i) Ne consegue che l’errore finale non sarà semplicemente la somma degli errori di arrotondamento, ma conterrà anche il contributo del condizionamento di ogni step dell’algoritmo considerato. Per calcolare invece il contributo dato dal primo termine tra parentesi quadre si deve svolgere come segue. Essendo ϕ(i) una funzione a pi` u variabili, si potr` a scrivere:  (i)  ϕ1 (u)  (i)   ϕ2 (u)  (i)   ϕ = ..  .   (i)

ϕn (u) Prendendo la j-esima di queste funzioni e analizzando il suo operare su u in virgola mobile si ottiene un risultato affetto dal seguente errore: (i)

(i)

con |j | ≤ eps

f l(ϕj (u)) = (1 + j )(ϕj (u))

Mettendo il tutto in una equazione pi` u compatta si potrà dunque scrivere: f l(ϕ(i) (u)) = (I + E(i+1) )(ϕ(i) (u)) Con 

(i+1)

1

0 ..

 E(i+1) =  0 10

. (i+1)

n

  

Questa è dunque l’espressione che interessa, ricordando che al posto di u ci sar` ax ˜(i) . In realt` a si adotterà qui l’approssimazione di supporre l’errore di arrotondamento piccolo, andando cioè a confondere x ˜(i) con x(i) . Dunque in definitiva [f l(ϕ(i) x ˜(i) − ϕ(i) (˜ x(i) )] ' Ei+1 ϕ(i) (x(i) ) Il termine a destra dell’uguale è chiamato anche αi+1 ed è chiamato errore assoluto di arrotondamento che nasce dal calcolo in virgola mobile di ϕ In definitiva, sommando i due fattori, si ottiene: ∆x(i+1) ' Ei+1 ϕ(i) (x(i) ) + Dϕ(i) (x(i) ) · ∆x(i) Quanto varr` a in definitiva ∆y ? Sarà legato sia al microcondizionamento che al macrocondizionamento, nonchè agli errori di arrotondamento. Chiamando αi l’errore di arrotondamento all’i-esima iterazione si avrà che: • per i=0 ∆x(0) = ∆x • per i=1 ∆x(1) = Dϕ(0) (x(0) )∆x + α1 • per i=2 ∆x(2) = Dϕ(1) (x(1) )[Dϕ(0) (x)∆x + α1 ] + α2 Si vede dunque che gli errori non si sommano ma si propagano. Dunque infine ricordando della proprietà della jacobiana di una composizione di funzioni, varr` a che: ∆y = Dϕ(x)∆x +

r X

Dψ (i) (x(i) )αi + Er+1 y

i=0

Si avr` a dunque un primo contributo, dato dal macrocondizionamento del problema ed indipendente dall’algoritmo, pi` u un termine dato dalle operazioni interne dell’algoritmo costituito da una composizione degli errori e della jacobiana della funzione composta ψ. Questa seconda parte è quella legata dunque alla stabilit` a dell’algoritmo, che sarà considerato accettabile ove l’incertezza legata all’algoritmo è confrontabile con quella generata a priori dal problema numerico. Si dirà in questo caso che l’algoritmo è ben fatto o robusto. Si faccia attenzione tuttavia al fatto che il concetto di robustezza e stabilit` a di un algoritmo è indipendente da quello di convergenza del problema, che è invece calcolata in termini ideali. Venendo infine al termine Er+1 y, si vede che esso rappresenta la coda dell’errore, l’arrotondamento sull’ultimo step dell’algoritmo, è evidentemente sempre presente ed il suo valore è maggiorato da Er+1 y ≤ eps|y|. Si avrà inoltre che qualunque sia il problema o l’algoritmo, a meno che in ingresso non si abbiano solo numeri di macchina esisterà anche un arrotondamento sui dati in ingresso. Si definisce in definitiva errore inevitabile la quantità: ∆(0) y = eps(|Dϕ(x)||x| + |y|) 11

Comunque sia fatto l’algoritmo, non è possibile andare al di sotto di questo errore, ed è per questo usato come termine di paragone. Si tenga conto che, in generale, non è possibile analizzare per intero, step per step, la stabilit` a di un algoritmo che abbia qualche migliaio di passaggi. Quel che si fa è, in generale, analizzare soltanto le situazioni considerate potenzialmente critiche.

1.4

Criteri di convergenza e troncamento di processi iterativi

Uno dei problemi principali in analisi numerica è quello di stabilire un criterio per determinare la convergenza di un problema iterativo. Tale procedimento è purtroppo effettuabile a priori solo nella risoluzione di problemi lineari. Per quanto riguarda invece i problemi non lineari questo non è sempre vero. Altra questione non da poco conto è quella di stabilire un criterio per il troncamento di un processo iterativo. Posto infatti che l’algoritmo porti a convergenza, ci si chiede quale sia il momento per arrestare il calcolo, definendo cos`ı la soluzione che, pur non essendo quella esatta, verrà considerata come definitiva. Bisogna fare attenzione alle caratteristiche intrinsecamente diverse di differenti algoritmi, che potranno convergere pi` u o meno velocemente gli uni rispetto agli altri. Risulta dunque in genere inopportuno arrestare un processo iterativo dopo un numero arbitrario di iterazioni, poichè la soluzione potr` a essere troppo o troppo poco precisa a seconda della velocità di convergenza dell’algoritmo utilizzato. Per questo è opportuno stabilire dei criteri di convergenza che diano informazioni su quale sia il momento migliore per il troncamento4 . Riguardo la velocit` a di convergenza è opportuno fare un discorso ulteriore. In numerica ogni cosa si paga, e di conseguenza un algoritmo che converge pi` u velocemente sarà, in genere, pi` u “costoso“, ovvero presenterà all’interno di ogni iterazione calcoli pi` u complessi od in numero maggiore. Il numero totale delle operazioni effettuate da un algoritmo è infatti quantificabile come NxM, ove N è il numero di iterazioni ed M è il numero di operazioni per iterazione. La scelta dell’algoritmo da utilizzare non è dunque banale e richiede accurate considerazioni di rapporto costi-benefici. Una strategia spesso attuata è inoltre quella di utilizzare differenti algoritmi in cascata per la risoluzione dello stesso problema. 4

Esistono poi casi particolari, come si vedr` a in seguito a proposito della teoria della diffusione, ove invece viene imposto un numero fissato di iterazioni. Le ragioni di questa scelta saranno pi` u evidenti quando ne verranno chiarite le implicazioni nel relativo capitolo

12

Esistono inoltre procedure di accelerazione di algoritmi lenti, che permettono di ridurre il numero di iterazioni senza incrementare eccessivamente il numero di operazioni per iterazione. Importante è segnalare che le zone dell’asse delle iterazioni non sono tutte uguali, ed in particolare è possibile definirne due: Zona delle prime iterazioni : qui si ha elevata differenza tra due soluzioni successive. Tale zona non è assestata, avvengono sconvolgimenti notevoli, e spesso non è possibile effettuare considerazioni sulla convergenza di un algoritmo, che è influenzata anche dal termine forzante. Zona asintotica : qui vi sono pochi sconvolgimenti, ma la convergenza è lenta e vi è differenza molto limitata tra due soluzioni successive. Dall’altro lato è per` o la velocità di convergenza è molto spesso stimabile a priori e dipende solo dalla matrice di iterazione dell’algoritmo. Soffermandosi sulla zona asintotica si nota che un algoritmo molto lento potrebbe indurre a credere di trovarsi già prossimi alla convergenza, quando invece se ne è ancora ben lontani Una soluzione a tale problema è quella di utilizzare il residuo, ovvero la differenza dall’eguaglianza che si ottiene andando a sostituire la soluzione approssimata nell’equazione di partenza. L’entità di tale residuo darà dunque, generalmente, informazioni utili sulla qualità della soluzione ottenuta. Si prenda ad esempio un generico problema del tipo Ax = b. Si definsice ~x soluzione esatta del problema, ove invece ~xn ne è la soluzione approssimata all’n-esima iterazione. Ne viene che: A¯ x − ~b = ~0 A~xn − ~b = ~r 6= ~0 Calcolare il valore di k ~r k può come detto aiutare a valutare il livello di convergenza della soluzione. Il problema del residuo tuttavia è che, per come esso è stato sinora definito, si tratta di una quantità assolulta, il che implica che la sua interpretazione non è immediata. Tale criterio verrà in generale utilizzato in AND con altri criteri di convergenza. Per quel che riguarda i criteri sulla soluzione in sé, esistono in generale due approcci: Norma : Non conoscendo la soluzione si utilizza uno pseudo-errore relativo: kΦm+1 − Φm k ≤ N kΦm+1 k Tramite questo metodo si perde il dettaglio sulla singola componente, concentrandosi invece sull’insieme dei valori. Esistono ovviamente numerose norme che possono essere applicate per questo criterio 13

Puntuale : In questo caso non si mescolano tra loro le componenti, ma per ciascuna viene valutata la convergenza: |Φm+1 − Φm i | i ≤ P,i m+1 |Φi |

∀i

Il problema è ora chiaramente spostato sul valore dell’errore massimo da imporre per il troncamento del processo iterativo, la cui scelta sarà ovviamente data dal rapporto costo/benefici. Il criterio puntuale è, ad esempio, pi` u stringente e porter` a in generale ad ottenere una soluzione pi` u precisa ma al costo di un numero maggiore di iterazioni. La scelta verrà comunque effettuata soprattutto in base a considerazioni fisiche. Nel caso, ad esempio, del calcolo delle sezioni d’urto medie sull’intervallo energetico non è richiesta la convergenza puntuale. Nel caso invece di sistemi di sicurezza ove vi sono vari parametri sotto controllo, sarà necessario applicare un criterio di convergenza puntuale per evitare imprecisione elevata su parametri molto delicati. Si noti in ogni caso che è possibile anche utilizzare pi` u criteri di convergenza, ad esempio puntuale su alcune elementi del vettore soluzione e in norma sulla soluzione presa interamente. Quali sono le possibilit` a che si hanno nel caso di fallimento di un algoritmo? • Nel caso in cui l’algoritmo diverga, ovvero non fornisca una soluzione, la causa pu` o avere due nature distinte: – Natura fisica: c’è un errore nella definizione del problema fisico che rende la soluzione non consistente – Natura numerica: c’è un errore nell’algoritmo • L’altro caso è quello in cui l’algoritmo fornisce comunque una soluzione, diversa tuttavia da quella esatta. Questi sono in effetti i casi pi` u pericolosi, perchè se in un contesto di modellizzazione di fenomeni osservati sperimentalmente è facile effettuare una verifica accorgendosi dell’errore5 , dall’altro lato se sono in fase di progettazione non è possibile fare questo controllo diretto fino a che non è realizzato un prototipo6 . 5

Oggigiorno il grande sviluppo dei metodi numerici e la loro sempre maggiore affidabilit` a stanno portando al verificarsi di condizioni opposte: talvolta è, in effetti, il metodo numerico che viene utilizzato per validare una determinata procedura sperimentale e non viceversa. Altro aspetto del potenziamento dei metodi numerici è che essi hanno talvolta permesso di prevedere fenomeni che solo successivamente è stato possibile osservare sperimentalmente 6 Anche in questo caso comunque non sempre le condizioni sono agevoli, anzi. Le situazioni tipiche dei reattori nucleari sono infatti particolarmente negative: non soltanto i prototipi sono particolarmente dispendiosi, ma anche una volta che questi vengono realizzati non è affatto facile effettuare misure di flusso neutronico per validare i calcoli numerici tramite dati sperimentali. Le uniche apparecchiature che permettono di effettuare tale op-

14

erazione sono i TIP (Travelling Incore Probes), sonde che viaggiano in canali dedicati e che permettono l’ottenimento di misure discrete. Dal punto di vista sperimentale spesso vengono realizzati dei sistemi di prova, poco corrispondenti alla versione industriale del progetto, ma la cui natura semplice e ben nota permette di agevolare la validazione dei modelli numerici. Sistemi di questo tipo sono detti benchmarks

15

Capitolo 2

Metodi numerici 2.1 2.1.1

Metodi di risoluzione di equazioni lineari con sorgente Metodo di Jacobi

Il metodo di Jacobi è un algoritmo iterativo che converge alla soluzione di un problema lineare non omogeneo. Esso costituisce l’algoritmo base, di cui si vedranno in seguito le evoluzioni. Viene in genere utilizzato solo nella fase turbolenta delle prime iterazioni per sgrossare la soluzione. Supponiamo di operare con una matrice a diagonale dominante. Si vedrà nelle prossime righe come effettivamente questa condizione è garanzia della convergenza dell’algoritmo. Il metodo di Jacobi prevede di spezzare la matrice A di interazione del problema in due matrici: • Una matrice D, diagonale e facilmente invertibile1 • Una matrice E, tale che A = D - E Dunque il problema diventa: DΦ = EΦ + S Il calcolo di D−1 è molto semplice, dunque potrò scrivere per questo problema: Φ = D−1 EΦ + D−1 S Ricordiamo che la moltiplicazione di una matrice o vettore qualsiasi per una matrice diagonale mi fornisce come risultato il prodotto di ciascun elemento di ogni riga per l’elemento della diagonale corrispondente. 1

si ricorda che una matrice diagonale è invertibile solo se tutti i suoi elementi sulla diagonale sono diversi da 0. Al tempo stesso è anche dimostrabile che, data una matrice invertibile, è sempre possibile effettuare permutazioni di riga e di colonna che portino all’ottenimento di una matrice avente solo elementi non nulli sulla diagonale

16

Chiamer` o ora: • B = D−1 E matrice di iterazione2 • q = D−1 S Da cui risulter` a Φ = BΦ + q Si vede che, per come è costruita, la matrice B è a diagonale nulla ed i suoi elementi diversi da zero sono sempre negativi ed inferiori ad 1 in modulo, a essendo essi definiti come bij = − aijii ∀ i 6= j Da tutto questo segue che kBk∞ < 1. La norma infinito è infatti definita come ilP valore massimo tra le somme degli elementi delle righe di una matrice (max( i |aij |). Nel momento in cui si trasforma quanto detto in un processo iterativo opero un’assegnazione, che mi porta ad una scrittura di questo tipo: Φ(m+1) = BΦ(m) + q Dunque si utilizza ogni volta il flusso dell’iterazione m-esima per calcolare quello dell’iterazione successiva. E’ stato dunque dato il via ad un processo iterativo. Si pu` o dimostrare che tale metodo porterebbe, dopo un numero infinito di iterazioni, alla soluzione esatta. Si prenda l’equazione di iterazione e si sottragga ad entrambi i membri la soluzione esatta. Ne verr` a: Φ(m+1) − Φ = BΦ(m) + q − Φ Ma essendo Φ = BΦ + q, si avrà Φ(m+1) − Φ = BΦ(m) + q − BΦ − q = B(Φ(m) − Φ) Lo scopo è chiaramente quello che sia kΦ(m+1) − Φk < kΦ(m) − Φk Partendo ora con Φ(0) = 0, Ne seguirà che:  (1) Φ =q     Φ(2) = BΦ(1) + q ..  .    (n+1) Φ = BΦn + q 2 da non confondere con la matrice di interazione, la matrice del problema orginiario, che nel nostro caso è A

17

Per come sono definite, le iterazioni possono anche essere scritte come:  (1) Φ =q     Φ(2) = Bq + q ..  .    (n+1) Φ = q + Bq + B 2 q + · · · + B n q Si vede dunque che il flusso neutronico all’n-esima iterazione dipende solo da B e da q. Da qui si vede che si ha di fatto una specie di serie geometrica matriciale, che nella pratica converge solo se kBk < 1. Dunque la soluzione di partenza è proprio il termine forzante, mentre ogni iterazione andr` a ad affinare la soluzione aggiungendo ulteriori termini. Attenzione per` o: il fatto che B sia cos`ı importante non è necessariamente una buona notizia. La matrice B è infatti strettamente dipendente da A, la quale è legata a sua volta alla fisica del problema. Dunque se B non converge, ci sar` a poco da fare: sarà necessario trovare una soluzione a monte dell’impostazione del problema fisico. Si pu` o ora notare un altro dettaglio. Andando infatti a prendere l’ultimo termine di ogni iterazione e li mettiamo a vettore, ho di fatto costruito il metodo delle potenze (vedi paragrafo 2.3.3 a pagina 33), che mi permette di calcolare gli autovalori di una matrice. Dunque l’ultimo termine dell’ultima iterazione tender` a, per numero infinito di iterazioni, a divenire parallelo all’ ϕ0 autosoluzione fondamentale del problema Bϕ = µϕ Ove µ è l’autovalore di modulo massimo e ϕ l’autovettore corrispondente. Questo discorso, per ora apparentemente di scarsa utilità, verrà utilizzato quando si parler` a del metodo SOR che costituisce, di fatto, un upgrade del metodo di Jacobi. Si dimostra ora un teorema di fondamentale importanza per questi problemi: Teorema 1 Condizione necessaria e sufficiente affinché il metodo di Jacobi converga è che ρ(B) < 1 Da qui deriva il fatto della convergenza del problema se la sua matrice di interazione è a diagonale dominante: si ha infatti una serie di consequenzialità a cascata che ci dicono che: • Se la matrice A è a diagonale dominante, allora vale che kBk∞ < 1 • Se kBk∞ < 1 allora ρ(B) < 13 • Se ρ(B) < 1 allora il metodo di Jacobi converge, tanto pi` u velocemente quanto pi` u il valore del raggio spettrale è piccolo 3

Questo vale in realt` a per ogni norma di matrice compatibile con una di vettore

18

Si vuole ora dimostrare il concetto fondamentale che sta alla base degli studi di convergenza del metodo di Jacobi: che la convergenza di tale metodo è associata al raggio spettrale della matrice di iterazione Per farlo è necessario definire ed enunciare le proprietà di una categoria di matrici, dette matrici irriducibili Definizione 1 Una matrice quadrata è irriducibile se per ogni coppia di i,j distinti esiste almeno una successione finita di interi, compresi tra 1 ed n (dove n è la dimensione della matrice) che va da i a j del tipo i = m0 → m1 → m2 → · · · → mq = j tale che: am0 ,m1 , am1 ,m2 , · · · , am(q−1) ,m1 sono tutti diversi da 0 Non è facile dimostrare che una matrice è irriducibile, mentre molto pi` u facile è dimostrare che una matrice non lo è. Di fatto comunque, a livello logico, l’irriducibilit` a di una matrice indica il fatto che tutti i nodi sono in qualche modo collegati tra loro. Ecco dunque perchè irriducibile: se un problema fisico è rappresentato da una matrice riducibile, esso può essere spacchettato (e dunque, ridotto) in pi` u problemi di ordine inferiore. Teoricamente sarebbe necessario dimostrare che per ogni coppia di elementi della discretizzazione esite un percorso fatto di linee verticali od orizzontali che collegano tali elementi passando solo attraverso elementi non nulli. Il problema è che dimostrarlo per tutte le coppie possibili non è affatto facile. Quale potrebbe essere un problema riducibile? Si tratta ad esempio del caso, sempre in ambito nucleare, di due zone di materiali diversi separate da una barriera di un materiale assorbitore perfetto. La matrice legata a tale problema sar` a dunque riducibile in quanto ciò che accade ai neutroni nella zona 1 non influenza in alcun modo ciò che invece avviene nella zona 2. 4 . Teorema 2 Se una matrice è riducibile è sempre possibile, permutando opportunamente righe e colonne, portarla ad una forma di tipo diagonale a blocchi: Aˆ1 ˆ0 ˆ0 Aˆ2 Viene ora mostrata una prima applicazione della proprietà di irriducibilità di una matrice: 4

Un esempio pratico di sistemi caratterizzati da matrici di riducibili sono i cosiddetti rubbiatroni, ove si opera in modo da creare delle specie di “gabbie stagne” per i neutroni per determinate energie

19

Teorema 3 Se una matrice A è irriducibile, allora affinchè valga che kBk∞ < 1 è sufficiente che la dominanza stretta della diagonale della matrice A sia valida per uno dei suoi elementi, mentre gli altri basta che soddisfino l’uguaglianza Dunque la convergenza di un problema avente come matrice di interazione una matrice irriducibile è assicurata anche se non è a diagonale dominante Si procede ora alla dimostrazione vera e propria della convergenza del metodo di Jacobi Noto che: (m) = Φ − Φ(m)

(m+1) = B(m) = · · · = B m (1)

→

Se B è una matrice simmetrica, esiste una n-upla di assi ortogonali rispetto ai quali essa assume forma diagonale. Essa è dunque diagonalizzabile5 , ove gli elementi di tale diagonale risultano essere gli autovalori del problema associato alla matrice stessa mentre i versori degli assi della n-upla ortogonale saranno gli autovettori, che formeranno un sistema completo ortonormale (vedi teorema 4 a pagina 32). Si prenda ora il vettore (1) . Poichè se la matrice è reale e simmetrica gli autovettori del problema agli autovalori Bϕ = λϕ formano una base ortonormale (vedi capitolo 2.4 a pagina 38) sarà sempre possibile scrivere il vettore errore assoluto come combinazione lineare di tali autovettori: (1)

=

n X

(1)

ch ϕh

h=1

Dunque (m+1)

= B

m

n X

(1)

ch ϕh =

h=1

=

n X

(1)

ch B m ϕh =

h=1

=

n X

(1)

ch λm ϕh

h=1

ove la j-esima componente di sarà dunque (m+1)

j

=

n X

(1)

ch λm h ϕh,j

h=1

Ma se ρ(B) < 1, allora qualunque autovalore λ ha modulo inferiore ad 1. Ne deriva evidentemente che, per m → ∞, il valore di (m) tenderà a 0, 5

tali autovalori saranno ripetuti ove essi abbiano molteplicit` a>1

20

e dunque la soluzione tenderà alla soluzione esatta. Attenzione che dalla stessa equazione si evince anche che non tutti gli elementi del vettore degli errori andranno a 0 con la stessa velocità, che dipende infatti dall’autovalore corrispondente (pi` u il modulo è elevato, pi` u lenta è la convergenza). Dimostrato che la condizione di ρ(B) < 1 è necessaria, di dimostra ora che essa è sufficiente. Si prenda un termine forzante q tale per cui l’errore alla prima iterazione sia tale che (1) k ϕ1 Ma se effettivamente questo è vero, allora si può scrivere: (m+1)

j

(1)

= c1 λm 1 ϕ1,j

Che ovvimante tende a 0 per m che tende ad infinito se e solo |λ1 | < 1 6 . Si pu` o fare questo ragionamento per ogni ϕi il che porta a poter dire che tale condizione è necessaria.7 Dopo aver dimostrato che la condizione di ρ(B) < 1 è, in generale, necessaria e sufficiente per la convergenza del metodo di Jacobi, accertiamoci ora che il raggio spettrale della matrice di iterazione mi dà anche una quantificazione della velocità di convergenza del metodo (nella zona asintotica). Si supponga ora che λ1 sia l’autovalore di modulo massimo. Questo implicher` a dunque che |λ1 | = ρ(B) e che |λj | < ρ(B) ∀j 6= 1. Si può dunque scrivere, come fatto in precedenza ma esplicitando il termine relativo 6

Cosa accade se la matrice di partenza A é reale e simmetrica? Questo implica che gli autovalori del problema Aϕ = λϕ sono a coppie in modulo uguale e segno opposto. Avremo cos`ı dunque che, per la convergenza, dovr` o considerare sempre due autovalori invece di uno. Al posto della equazione appena scritta avremo che il calcolo dell’errore sar` a effettuato sulla norma, e potremo scrivere che: k(m+1) k ∼ Cρ(B)m 7

Si noti che questo discorso è stato effettuato per dei termini forzanti specifici, ovvero tali per cui l’errore risulti parallelo ad un autovettore. Questo per potersi porre nel caso pi` u sfortunato, e di conseguenza definire una regola generale. Esisteranno dunque dei valori di q per cui tale condizione non sar` a, invece, necessaria. Basti pensare al caso in cui, per assurdo, λ1 > 1 ma (1) ⊥ ϕ1 . In questo caso il valore dell’errore all’m-esima iterazione sar` a dato da: n X (m+1) (1) j = ch λm h ϕh,j h=2

Ma dunque l’unico autovalore maggiore di 1 non interverr` a nella sommatoria e di conseguenza avr` o convergenza pur essendo la condizione di ρ(B) < 1 non rispettata. Tuttavia nel mondo reale bisogna ricordarsi che non solo che il calcolatore commette errori di arrotondamento e non conosce la quantit` a “0”, e di conseguenza il concetto di ortogonalit` a è mantenuto solo per poche iterazioni

21

all’autovalore di modulo massimo, che : (m+1) j

=

(1) λm 1 [c1 ϕ1,j

+

n X

(1)

ch (

h=2

λh m ) ϕh,j ] λ1

Si vede chiaramente che, per m → ∞, sarà (1)

m+1 ' λm 1 c1 ϕ1,j j m

In quanto il termine λλh1 tende a 0. Maggiore sarà il valore dell’h-esimo autovalore, maggiore sar` a la resistenza al contributo corrispondente a calare con le iterazioni. A questo proposito si definisce rapporto di dominanza il modulo del rapporto tra il secondo autovalore di modulo maggiore e l’au2| tovalore di modulo massimo ( |λ u tale rapporto di dominanza |λ1 | ). Tanto pi` è elevato tanto pi` u è necessario aspettare affinchè l’errore alla m-esima iterazione sia parallelelizzato al contributo dell’autovalore fondamentale. Anche il rapporto di dominanza, come il raggio spettrale, è definito dalla fisica del problema (oltre che dalla scelta dell’algoritmo). Il punto è che, in fisica della diffusione, i due autovalori di modulo maggiore sono molto vicini tra loro. ` di convergenza asintotica il valore8 : Si definisce velocita R∞ = − ln ρ(B) Si suppunga a questo punto di trovarsi in zona asintotica, ad una iterazione m0 . Com’è possibile capire quante altre iterazioni bisogna fare per ridurre l’errore di un fattore 100? Questo sarà come dire che ρ(B)m < 0.01 da cui si otterr` a che dovr` a essere: m>

4, 76 − ln ρ(B)

Il metodo di Jacobi è dunque cos`ı esplicitato: molto semplice e lineare. Si vuole ora per` o analizzarne i punti deboli, in modo tale da poter capire ove si possa lavorare per cercare di migliorarlo. Il suo handicap fondamentale è legato alla sua scarsa ottimizzazione nell’utilizzo delle informazioni disponibili. Si prenda infatti una generica m(m+1) esima iterazione. Arrivati al calcolo di Φi si avrà a disposizione già una serie di informazioni “di qualità”, vale a dire tutti i flussi della m+1-esima (m+1) iterazione fino al Φi−1 . Il metodo di Jacobi è dunque detto “delle iterazioni successive” poichè per ogni iterazione uso solo informazioni “vecchie”. Il punto di novità del metodo di Gauss-Seidel è proprio quello di utilizzare, per il calcolo del flusso aggiornato, anche le informazioni nuove. 8

Si noti che questa non costituisce un calcolo della velocit` a di convergenza ma solo una sua stima, utile per confrontare tra loro differenti algoritmi

22

2.1.2

Metodo di Gauss-Seidel

Come è dunque possibile effettuare un upgrading del metodo di Jacobi, con poco sforzo e grande risultato? In un generico nodo (prendiamo ad esempio l’11◦ , supponendo che i nodi adiacenti nelle direzioni N S W E siano rispettivamente 15, 7, 10, 12) sarà funzione del flusso nei 4 nodi adiacenti, ovvero 15 10 11 12 14 Φ11 = f (Φ7 , Φ10 , Φ12 , Φ15 ) Trasportato nella riga corrispondente all’interno delle iterazioni del metodo di Jacobi avremo: (m+1)

Φ11

(m)

(m)

(m)

(m)

= f (Φ7 , Φ10 , Φ12 , Φ15 )

Tuttavia, poichè all’interno di ogni iterazione i flussi sono calcolati in modo ordinato dal primo all’ultimo, al momento di calcolare il flusso in 11 avremo già noti i valori all’m+1-esima iterazione del flusso in 7 e 10. Si può dunque sostituire a tali valori il valore del flusso m+1-esimo La matrice B del metodo di Jacobi è del tipo:   0 T∗   .. B = T + T∗ =   . T

0

Ove le matrici T e T ∗ sono evidentemente rispettivamente triangolare inferiore e superiore. Con il sistema originario dato da: Φ = BΦ + q Che pu` o essere riscritto come: Φ = (T + T ∗ )Φ + q Sfruttando quanto detto riguardo il metodo di Gauss Seidel, ovvero che per il calcolo della i-esima componente sono già state calcolate le i-1 componenti precedenti, il processo iterativo relativo al metodo di Gauss-Seidel diverrà, a seguito del processo di assegnazione, della forma: (I − T )Φ(m+1) = T ∗ Φ(m) + q Da cui Φ(m+1) = (I − T )−1 T ∗ Φ(m) + (I − T )−1 q 23

Risulter` a dunque, sottraendo membro a membro l’espressione iterativa e l’espressione esatta: (m+1) = (I − T )−1 T ∗ (m) Ove (I − T )−1 T ∗ è la matrice di iterazione del metodo di Gauss-Seidel. Si vedano ora le equazioni per componenti: per il metodo di Jacobi si aveva: n X aij (m) (m+1) Φi =− Φ + qi aii j j=1

Con Gauss-Seidel bisogna invece spezzare le sommatorie, poichè una sarà riferita al “nuovo” ed un’altra al “vecchio”: (m+1)

Φi

=−

i−1 X aij j=1

aii

(m+1)

Φj

−

n X aij (m) Φ + qi aii j

j=i+1

Il metodo di Gauss-Seidel risulta dunque, asintoticamente, pi` u efficace. Si dimostra infatti che la sua velocità di convergenza è pari a ρ(B)2 Si noti che l’utilizzo del metodo di Gauss-Seidel non incrementa lo sforzo computazionale: il numero di operazioni per iterazione è lo stesso del metodo di Jacobi.

2.1.3

Metodo di Jacobi vs Metodo di Gauss

Dunque il metodo di Gauss-Seidel converge pi` u velocemente di quello di Jacobi, e lo fa a costo zero. La domanda che sorge legittima è la seguente: perchè si parla ancora di metodo di Jacobi? Un motivo esiste: il metodo di Gauss-Seidel introduce una asimmetria nel sistema, in quanto il flusso in ogni nodo verr` a calcolato utilizzando informazioni di qualità diversa a seconda della direzione. Questo è un difetto che risulta pienamente accettabile nella zona asintotica, ma non in quella delle prime iterazioni ove si andrà invece in generale ad utilizzare il metodo di Jacobi. In un caso particolarmente felice ove il valore di ρ(B) fosse particolarmente ridotto, ci si potrebbe addirittura ritrovare nella condizione di convergenza talmente veloce da rendere inutile se non dannosa l’applicazione del metodo di Gauss-Seidel. Esistono per` o altre condzioni in cui si utilizza Jacobi piuttosto che GaussSeidel, e sono quelle in cui vi sono due catene di iterazioni una dentro l’altra. In questo caso è possibile che il processo iterativo interno sia tale da richiedere un’unica iterazione per ogni iterazione esterna9 . In questo caso si andr` a ovviamente ad applicare Jacobi, sempre per la stessa ragione: nelle fasi iniziali di applicazione di un algoritmo la convergenza non è garantita ed è necessario evitare di introdurre asimmetrie nella risoluzione. L’applicazione del metodo di Gauss-Seidel in questo frangente potrebbe portare all’ottenimento di soluzioni non fisiche. 9

Si tratta come vedermo del caso della teoria della diffusione

24

Altra situazione caratteristica è quella dell’applicazione dei due metodi a scacchiera. In questo modo, calcolando il flusso tramite Jacobi solo sui nodi “neri” ci si ritrova che tutti i nodi “bianchi” sono circondati solo da informazioni nuove che possono cos`ı essere utilizzate per calcolarvi il flusso. Questo metodo di fatto è migliore di quello di Jacobi anche alle prime iterazioni. Si vede inoltre che la velocità di convergenza del metodo “a scacchiera” nella zona asintotica è esattamente la stessa di quella del metodo di Gauss-Seidel standard.

2.1.4

Metodo SOR

Il metodo SOR lavora tramite estrapolazione, ovvero un processo di creazione di informazioni sulla base di altre di cui si è già in possesso. Questo passaggio, a differenza di quello che portava dal metodo di Jacobi a quello di Gauss-Seidel, include un costo. Ne segue che tale metodo, dotato di velocità di convergenza superiori, è adatto ad un uso nelle fasi finali di affinamento della soluzione, dopo aver sgrossato il problema tramite l’uso dei metodi visti in precedenza. Il metodo SOR è anche detto di sovrarilassamento, e si basa sul forzare l’upgrade della soluzione da un passaggio all’altro. Riscrivendo infatti l’evoluzione dell’errore relativo tramite metodo di Gauss-Seidel: [Φ(m+1) − Φ(m) ]GS = T Φ(m+1) + T ∗ Φ(m) + q − Φ(m) Partendo da questo presupposto, l’upgrade della soluzione tramite metodo SOR diventa: [Φ(m+1) − Φ(m) ]SOR = ω[Φ(m+1) − Φ(m) ]GS Ove ovviamente il caso di ω = 1 porta al riottenimento al metodo di GaussSeidel. ω è detto parametro di forzamento Espandendo quanto scritto si ottiene: Φ(m+1) = ωT Φ(m+1) + [(1 − ω)I + ωT ∗ ]Φ(m) + ωq Da cui, in forma esplicita:

Φ(m+1)

matrice di iterazione z }| { = (I − ωT )−1 [(1 − ω)I + ωT ∗ ] Φ(m) + (I − ωT )−1 ωq

Si pu` o dimostrare che il valore ottimale per ω è dato da: ωopt =

2 p 1 + 1 − ρ(B)2

E visto che si pu` o fare, facciamolo. 25

L’obbiettivo è trovare un’espressione per il raggio spettrale della matrice di iterazione del metodo SOR. Una volta fatto questo, il passaggio successivo sar` a evidentemente quello di minimizzare tale valore, che come sappiamo rappresenta la velocit` a con cui l’algoritmo tende a convergenza. Il problema sar` a del tipo Jω ψ = γψ ove ci` o che ci si prefigge di calcolare è il γ di modulo massimo. Richiamando la definizione della matrice di iterazione: Jω = (I − ωT )−1 [(1 − ω)I + ωT ∗ ] Se ne pu` o ricavare una ulteriore equazione agli autovalori. Infatti vale Jω ψ = γψ (I − ωT )

−1

[(1 − ω)I + ωT ∗ ]ψ = γψ [(1 − ω)I + ωT ∗ ]ψ = (I − ωT )γψ ψ − ωψ + ωT ∗ ψ = γψ − ωT ψ γ−1+ω (γT + T ∗ )ψ = ψ ω

Si pu` o ora riscriverla componente per componente, tramite il doppio indice, ove con i simboli dei punti cardinali si indicheranno gli elementi delle matrici posizionati adiacenti all’elemento di interesse. Dunque γWi,j ψi−1,j + γSi,j ψi,j−1 + Ei,j ψi+1,j + Ni,j ψi,j+1 =

ω−1+γ ψi,j ω

Viene ora sfruttata l’arbitrarietà degli autovettori per scrivere un legame tra ψ e Φ, al fine di ottenere una espressione simile a quella scritta per il metodo di Jacobi e poter dunque confrontare i due algoritmi. Si userà in i+j particolare che ψij = Φij γ − 2 ( per Jacobi Wi,j Φi−1,j + Ni,j Φi,j+1 + Ei,j Φi+1,j + Si,j Φi,j−1 = λΦi,j √ per SOR Wi,j Φi−1,j + Ni,j Φi,j+1 + Ei,j Φi+1,j + Si,j Φi,j−1 = ω−1+γ ω γ Φi,j Confrontando i due ottengo: λ=

ω−1+γ √ ω γ

Risolvendo per γ, diventa √ λω γ = ω − 1 + γ λ2 ω 2 γ = ω 2 + 1 + γ 2 + 2ωγ − 2ω − 2γ γ 2 − (2 + λ2 ω 2 − 2ω) + (ω − 1)2 = 0 26

da cui si avranno due soluzioni: q 2 2 λ2 ω 2 γ1,2 = (1 − ω − 2 ) ± (1 − ω − λ 2ω )2 − (1 − ω)2 = q 2 2 4 4 = (1 − ω − λ 2ω ) ± λ 4ω + λ2 ω 2 − ωλ2 ω 2 = q 2 2 2 2 = (1 − ω − λ 2ω ) ± λω λ 4ω + 1 − ω Ponendosi nel caso tipico di matrice simmetrica gli autovalori di B saranno tutti reali. Essendo inoltre che, per ipotesi, λ2 < 1 (altrimenti la convergenza del metodo di Jacobi non è assicurata) e ω > 1 (altrimenti il metodo SOR rallenterebbe la convergenza invece di accelerarla10 ), è cos`ı definito l’intervallo di variazione di ω Si studino ora i tre casi separatamente, avendo chiamato ω0 il valore di ω che annulla il discriminante dell’equazione: ω0 =

2 √ 1 + 1 − λ2

si ottengono 3 casi possibili: • Per ω = ω0 avremo che: γ = 1 − ω0 +

λ2 ω02 2

che, sostituendo il valore di ω0 , diviene γ = ω0 − 1 • Per ω0 < ω < 2 si ha ∆ < 0, il che porta ad ottenere per γ due radici complesse coniugate. Attraverso gli opportuni passaggi si ottiene che |γ1 | = |γ2 | = ω − 1 che, essendo per ipotesi ω > ω0 , darà luogo ad un valore di γ maggiore che nel caso precedente • Per 1 < ω < ω0 si ottengono due soluzioni reali distinte, delle quali quella di modulo maggiore sarà legata al discriminante preso come positivo, che ci permette di ottenere: r 1 2 2 λ2 ω 2 γ1 = (1 − ω + λ ω ) + λω 1 − ω + 2 4 Si vede che questo genera un andamento il cui minimo risulta essere situato proprio in ω = ω0 10 Attenzione, questo non è sempre vero. A volte infatti gli algoritmi vengono rallentati invece che accelerati, in particolare in quei casi ove ho fenomeni che si svolgono su scale temporali distinte, ed il rischio è che le “qualit` a” dei risultati non risultino coerenti. Un esempio tipico è l’accoppiamento di fenomeni elettro-magnetici, estremamente veloci, ed i fenomeni termo-fluidodinamici, pi` u lenti

27

Dunque è stato dimostrato come il valore ottimale per ω sia dato da: ωopt = ω0 =

2 √ 1 + 1 − λ2

Ecco dunque identificato il costo aggiuntivo del metodo SOR: se fino ad ora si era parlato del raggio spettrale della matrice B come qualcosa che identificava la velocit` a di convergenza dell’algoritmo, senza però mai aver manifestato la necessit` a di calcolarlo, ecco che qui è invece necessario conoscerne il valore al fine di applicare il metodo SOR in maniera ottimale assicurandosi cos`ı la migliore velocit` a di convergenza possibile. Bisogna in questa fase ricordare le ipotesi fatte per ottenere questo risultato: • Griglia strutturata • Matrice simmetrica Si pu` o dimostrare, ma non verrà fatto in questa sede, che il valore della velocit` a di convergenza del metodo SOR dipende dalla numerazione scelta, ed in particolare quella che garantisce i risultati migliori è la numerazione naturale.

2.1.5

Shift tra metodo di Jacobi e metodo SOR

Ci si pone nella condizione tipica in cui si utilizza il metodo di Jacobi nella fase iniziale, per lo sgrossamento della soluzione, per poi passare al metodo SOR per l’affinamento finale. Sorgerà dunque la questione su quando abbandonare il primo metodo per passare al secondo. Si vede che in base al problema da risolvere11 varierà la velocità di convergenza del metodo di Jacobi. Pi` u tale velocità è elevata, minore sarà il numero di iterazioni necessario per il passaggio al metodo SOR. Esiste tuttavia un altro elemento. Ricordiamo infatti quanto detto in modo molto innocente all’inizio della descrizione del metodo di Jacobi, ovvero che gli ultimi elementi di ogni iterazione successiva costituiscono di fatto l’applicazione del metodo delle potenze sulla matrice B che porta al calcolo degli autovalori della stessa. Appare ora evidente che dunque, nell’applicare Jacobi, si opera involontariamente il calcolo iterativo di B n q, il che permette di calcolare un valore approssimato di ρ(B) che si potrà dunque utilizzare per il calcolo di ωopt . Questa riflessione dice dunque che al fine di determinare il momento per lo shift da metodo di Jacobi e metodo SOR non ci si concentrer` a soltanto sulla convergenza del metodo di Jacobi in sè, ma anche su quella del metodo delle potenze per il calcolo di ρ(B), il cui valore è fondamentale per la corretta applicazione del metodo SOR 11

e dunque in base alla matrice A → la matrice B → ρ(B)

28

2.1.6

Ordinamento dei nodi di flusso

Si accenner` a nelle pagine successive alla necessità, una volta che si passa ad una logica multidimensionale, della scelta di una numerazione per i nodi di flusso in modo da passare da una forma del tipo Φi, j ad una Φl , che mi permette l’utilizzo standard di vettori e matrici. A questo proposito sono mostrate di seguito le tre tipologie pi` u utilizzate, a ciascuna delle quali è associata la corrispondente forma della matrice risultante: Numerazione naturale 

3 4 1 2

 0 X X 0  X 0 0 X     X 0 0 X  0 X X 0

→

Numerazione a scacchiera  0 0 X X  0 0 X X     X X 0 0  X X 0 0 

4 2 1 3

→

Numerazione ciclica  0 X 0 X  X 0 X 0     0 X 0 X  X 0 X 0 

4 3 1 2

→

Tutte e tre le matrici avranno gli stessi autovalori (0 con molteplicità 2 e ±2X con molteplicit` a 1. Per l’applicazione del metodo di Jacobi le tre strategie sono equivalenti, nel senso che portano alla stessa identica velocità a convergenza. Per quanto riguarda invece il metodo di Gauss-Seidel avremo che, per ρ(B) = 0.5, sar` a: • JGS = 0.25 con numerazione naturale (come da considerazioni precedenti pari a ρ(B)2 ) • JGS = 0.25 con numerazione a scacchiera (come da considerazioni precedenti pari a ρ(B)2 ) • JGS = 0.28 con numerazione ciclica

29

2.2

Risoluzione di sistemi di equazioni non lineari

I processi non lineari vengono trattati, in analisi numerica, tramite una linearizzazione iterativa. L’idea è dunque quella di, partendo da una f (x) = 0 con f (x) non lineare, andare a risolvere una fˆ(x) = 0, con fˆ(x) ovviamente parente stretta di f (x). Questo procedimento comporta ovviamente una complicazione della risoluzione, che aumenta con l’aumentare del numero delle equazioni.

2.2.1

Metodo di Newton per la risoluzione di sistemi non lineari

Si supponga di avere un generico sistema di equazioni non lineare.  f1 (x1 , x2 , · · · , xn )     f2 (x1 , x2 , · · · , xn ) f (x) = =0 ..  .    fn (x1 , x2 , · · · , xn ) La prima cosa che viene in mente di fronte ad un sistema non lineare è quella di linearizzarlo tramite uno sviluppo in serie di Taylor troncato al prim’ordine. Si chiami ξ la soluzione esatta e x(i) la soluzione approssimata alla iesima iterazione. Si avr` a ovviamente che f (ξ) = 0 ed invece f (x(i) ) 6= 0 . (i) Essendo dunque x una approssimazione di ξ, sarà lecito supporre che i due vettori siano molto “vicini” tra loro, di modo da poter calcolare tramite linearizzazione in serie di Taylor la funzione nell’intorno di x(i) . Vale infatti che, trascurando i termini di grado superiore al primo, si ottiene un sistema del tipo:  (i) ∂f1 (i) ∂f1  f1 (x(i) ) + (ξ1 − x1 )( ∂x ) (i) + · · · + (ξn − xn )( ∂x ) (i) ' 0  n x=x 1 x=x    f2 (x(i) ) + (ξ1 − x(i) )( ∂f2 ) (i) + · · · + (ξn − xn(i) )( ∂f2 ) (i) ' 0 1 ∂x1 x=x ∂xn x=x ..   .    (i) (i) ∂fn n fn (x(i) ) + (ξ1 − x1 )( ∂f ∂x1 )x=x(i) + · · · + (ξn − xn )( ∂xn )x=x(i) ' 0 Se in questo sistema operiamo una assegnazione potremo innescare un procedimento iterativo. In particolare otterremo questo scrivendo:  (i) ∂f1 (i) ∂f1  f1 (x(i) ) + h1 ( ∂x ) (i) + · · · + hn ( ∂x ) (i) ' 0  n x=x 1 x=x   ∂f2  f2 (x(i) ) + h(i) ( ∂f2 ) (i) + · · · + h(i) n ( ∂xn )x=x(i) ' 0 1 ∂x1 x=x .  ..     (i) ∂fn (i) n fn (x(i) ) + h1 ( ∂f ∂x1 )x=x(i) + · · · + hn ( ∂xn )x=x(i) ' 0

30

(i)

Di questo sistema l’unica incognita è proprio h(i) , ove si è definito hj (i+1) (i) (xj − xj ) Definendo il passaggio da x(i+1) = x(i) + h(i) è possibile innescare

=

una iterazione all’altra in modo che un processo iterativo, i cui passaggi

saranno: 1. Calcolo dello Jacobiano nel punto x(i) 2. Calcolo delle f (x(i) ) 3. Calcolo di h tramite soluzione del sistema 4. Aggiornamento di x(i+1) Un problema di tale procedimento è che, ove il numero delle incognite sia elevato, J pu` o diventare molto grande. Se le sue dimensioni sono ridotte sar` a possibile scriverlo in forma analitica (di fatto manualmente), mentre ove le dimensioni aumentino si passerà in genere ad una scrittura in forma numerica: (

fj (x(i) + ek hjk ) − fj (x(i) ) ∂fj (i) )x=x(i) = Bj,k = ∂xk ek hjk

Importante è la valutazione dell’intervallo di discretizzazione. Se questo è troppo grande, si perde ovviamente in precisione, mentre ove questo sia troppo piccolo si corre il rischio di cancellazione numerica. Come sempre la scelta risulter` a dal compromesso giudicato migliore. Di fatto la scrittura dello Jacobiano e la risoluzione del problema associato sono le fasi pi` u critiche del processo iterativo. Per questo una possibilità è quella di non aggiornarlo ad ogni iterazione. La risoluzione del problema avente lo Jacobiano come matrice di iterazione dipende dallo Jacobiano stesso: se si tratta di una matrice si usano metodi diretti, in caso contrario si privilegiano invece metodi iterativi. Tuttavia, qualunque sia il metodo scelto ed il caso specifico, questo metodo di risoluzione converge con estrema difficolt` a.

2.3 2.3.1

Metodi per problemi agli autovalori Scomposizione QR di una matrice

Si definisce trasformazione di similitudine una qualsiasi trasformazione del tipo: B = SAS −1 Caratteristiche peculiari di questa trasformazione sono che: • La matrice A e la matrice B hanno gli stessi autovalori

31

• Dato x autovalore di A, y= Sx è autovettore della nuova matrice B L’utilit` a di tali trasformazioni di similitudine è data anche dai seguenti teoremi: Teorema 4 Data una matrice A reale e simmetrica, esiste una matrice ortogonale V1 tale che D = V1 AV1T ove D è una matrice diagonale Teorema 5 Data una matrice A reale, esiste una matrice ortogonale V2 tale che H = V2 AV2T ove H è una matrice di Hessenberg, che diventa tridiagonale se A è simmetrica Si definisce a questo punto una nuova tipologia di matrice, detta riflettore, nella maniera seguente: √ U = I − 2uuT con uT u = 1 Si dimostra facilmente che tale matrice presenta le seguenti proprietà: • è simmetrica (U T = U ) • è ortogonale (U T = U −1 ) • è involutoria (U 2 = 1) Ulteriore propriet` a di tale matrice, che utilizzeremo in seguito, è enunciata dal seguente teorema: Teorema 6 Dato un generico vettore non nullo x, esiste un riflettore elementare U tale che U x = −σe1 Dove

  u = x + σe1 T e1 = (1, √0, · · · , 0)  σ = ± xt x

Dalle considerazioni precedenti si è visto che è possibile costruire delle matrici U tali che, moltiplicate per un vettore, ne annullano tutte le componenti tranne la prima. Si può dunque pensare di costruire una matrice U1 tale che moltiplicata per una generica matrice A fornisca, in uscita, una matrice con tutti gli elementi della prima colonna nulli fatta eccezione per la prima riga. Moltiplicando poi per una analoga U2 , si possono annullare 32

anche tutti gli elementi della seconda colonna ad eccezione dei primi due. Il metodo di costruzione della generica matrice Uk è il seguente: ! Ik−1 0 (n−k+1) 0 U1 Il risultato finale di tutte le moltiplicazioni sarà dunque una matrice triangolare superiore, che verr` a chiamata R. Moltiplicando invece tra loro tutti gli n-1 riflettori otterremo una matrice QT che risulterà ortogonale, in quanto prodotto di matrici ortogonali. Avrò quindi che: QT A = R

7−→

A = QR

Si tratta dunque di un processo che permette di ottenere un risultato analogo alla decomposizione di Gauss, avendo come unica pecca l’utilizzo di 2n3 /3 operazioni, il doppio esatto di quelle necessarie per la decomposizione di Gauss.

2.3.2

metodo QR

Il metodo QR è di tipo iterativo. Si definisce A1 la matrice A, e calcoliamo delle nuove matrici secondo il processo seguente: 1. Decomposizione della matrice all’i-esima iterazione: Ai = Qi Ri 2. Calcolo della matrice per l’iterazione successiva: Ai+1 = Ri Qi Tutte le matrici cos`ı calcolate risulteranno simili ad A, in quanto risulteranno di fatto dalla moltiplicazione a destra per le Qi ed a sinistra per le loro trasposte: Ai+1 = (Q1 Q2 · · · Qi )T A(Q1 Q2 · · · Qi ) Trattandosi di operazioni di similitudine, si conservano gli autovalori! Si può dimostrare che la matrice An tenderà all’infinito a diventare una matrice triangolare superiore, i cui autovalori saranno gli elementi della diagonale principale. Il problema di questo calcolo è che è lungo e dispendioso. Il procedimento pi` u generalmente utilizzato è quello di trattare preventivamente la matrice A, cercando di farla diventare tridiagonale o di Hessenberg, per poi applicare il metodo QR ad una condizione ove si ha un numero consistente di elementi nulli.

2.3.3

Metodo delle potenze

Il metodo delle potenze è un algoritmo di tipo iterativo che converge ad un vettore parallelo all’autovettore fondamentale della matrice di interesse. Pu` o essere applicato a problemi di tipo: 33

• AΦ = λΦ • AΦ = λBΦ • D(λ)Φ = 0 Ove nel terzo caso non è garantita la linearità del problema rispetto agli autovalori. Il metodo delle potenze è un algoritmo iterativo che, a convergenza, fornisce una stima dell’autovettore fondamentale, utile in particolare per due classi specifiche di problemi: • Problemi agli autovalori di fusione neutronica in mezzo moltiplicante • Utilizzo in metodi per la risoluzione di sistemi di equazioni non omogenei Si supponga che la matrice M abbia un insieme completo di autovettori. A partire da tale ipotesi si può dimostrare che la successione M m converge all’autovalore di modulo massimo qualunque sia q. Si supponga inoltre che la matrice M abbia tutti autovalori distinti: esister` a dunque uno ed un solo autovalore avente modulo superiore a tutti gli altri. Avendo inoltre imposto che la matrice M abbia un set completo di autovettori, è ragionevolmente possibile supporre di poter scrivere q=

n X

ch ϕh

h=1

Ove i ϕh sono evidentemente gli autovettori della matrice M.12 Si avr` a dunque che: M mq =

n X

ch M m ϕh =

h=1

=

n X

ch µm h ϕh =

h=1

=

µm 1 [c1 ϕ1

+

n X µh h=2

µ1

ch ϕh ]

Ma dunque, poichè la sommatoria tenderà a 0 per indice m che tende ad infinito, il vettore M m q tenderà a diventare parallelo a ϕ1 . Come per quanto detto relativamente al metodo di Jacobi, si avrà anche qui un termine (quello legato al secondo autovalore di modulo massimo) che resisterà a questa 12

Tale scrittura equivale a dire che è sempre possibile scrivere un qualunque vettore q come combinazione lineare degli autovettori della matrice M

34

convergenza pi` u di tutti gli altri con una “forza” quantificata dal rapporto di dominanza Attenzione per` o. Nei passaggi precedenti è stato dimostrato che il metodo converge ad un vettore parallelo all’autovettore fondamentale, ma si vede che compare nella formula risolutiva l’elemento µm 1 che moltiplica tutta l’equazione. Risulta dunque evidente che per µ1 > 1 la formula tenderà a divergere per m grandi, mentre per µ1 < 1 essa tenderà a 0. Per ottenere dunque una soluzione è opportuno, ricordandosi che gli autovettori sono definiti a meno di una costante moltiplicativa arbitraria, normalizzare la soluzione ad ogni step iterativo. Per quanto visto prima, è evidente che affinchè la soluzione converga ad un valore reale diverso da 0 è necessario che il termine a destra elevato all’m-esima potenza sia il pi` u possibile vicino ad uno. Verr` a dunque introdotto come elemento normalizzante per ogni iterazione l’autovalore di modulo massimo calcolato all’iterazione precedente. A questo punto per` o sar` a necessario capire come calcolare il raggio spettrale. Si è infatti dimostrato che il vettore M m q tende a diventare parallelo all’autovettore fondamentale, ma non è ancora stato fatto il passo successivo. Infatti se il problema di partenza è lineare sarà possibile sfruttare il fatto che se ϕ associato a µ è soluzione al problema agli autovalori, allora lo stesso varr` a per cϕ, ovvero qualunque vettore è soluzione del problema agli autovalori a patto che sia parallelo ad uno dei suoi autovettori. Sulla base di queste considerazioni si potr` a dunque dire che, supponendo di essere giunti a convergenza all’m-esima iterazione: M (M m q) ' µ1 (M m q) Ove il vettore M m q corrisponde al generale multiplo di autovettore cϕ Dall’equazione precedente si può dunque dedurre che µ1 '

(M m+1 q)1 (M m+1 q)2 (M m+1 q)n ' ' · · · ' (M m q)1 (M m q)2 (M m q)n

Si ha infatti che a convergenza i due vettori M m q ed M M m q devono essere paralleli in quanto uguali a meno di una costante moltiplicativa (l’autovalore). Questo fornisce dunque un efficace criterio di convergenza per il mio metodo iterativo, che si potrà basare sulla similitudine tra i vari rapporti tra le componenti. Quando essi saranno sufficientemente simili tra loro si potr` a finalmente considerare di essere giunto a convergenza e poter dunque troncare il processo iterativo. Occorre tuttavia quantificare questa vicinanza. Si ricorre in generale ad una imposizione sul numero delle cifre della mantissa che devono essere uguali. Si ha in questo modo imposto una sorta di errore relativo, che infatti non dipende dal valore in sè della quantità.13 . Occorre tuttavia fare come 13

Si tratta evidentemente in questo caso di un criterio di convergenza puntuale. Esistono anche criteri di convergenza in norma che saranno in generale meno stringenti

35

sempre grande attenzione agli algoritmi molto lenti, che portano ad avere le soluzioni di due iterazioni successive molto vicine tra loro. In questo caso si rischia che il criterio di convergenza, per quanto stringente, non porti a selezionare una soluzione realmente vicina a quella cercata. In casi come questo è dunque sempre opportuno andare a sostituire i risultati ottenuti nell’equazione di partenza per verificarne l’effettiva veridicità. Rimane tuttavia il dubbio legato alla scelta di quale degli n valori possibili scegliere per l’autovalore, ognuno dei quali è di fatto teoricamente equivalente (se mi trovassi di fronte alla soluzione analitica essi sarebbero infatti tutti uguali). Una possibilità è data dall’operazione seguente: Pn m+1 q) < M m+1 q, I > i i=1 (M P µ1 ' = n m q) < M m q, I > (M i i=1 Tale operazione fornisce infatti una stima dell’autovalore fondamentale pi` u precisa, a parit` a di iterazione, di qualunque stima effettuata dal semplice rapporto delle componenti. Si nota per` o che questo calcolo sta equiparando tra loro tutte le componenti dell’autovettore fondamentale. Tuttavia, se per la matematica le variabili sono tutte uguali, cos`ı non è per la fisica e l’ingegneria. Esisterà dunque un modo per sfruttare ci` o che si conosce riguardo la fisica del problema al fine di rendere ancora pi` u precisa la soluzione. Si introducono dunque i quozienti di Rayleigh, definiti in modo che: µ1 '

< q (m+1) , q (m) > < q (m) , q (m) >

oppure

µ1 '

< q (m+1) , q (m+1) > < q (m) , q (m+1) >

avendo utilizzato quella che viene detta una pesatura alla Galerkin, che prevede l’utilizzo della stessa classe funzionale della soluzione per la pesatua delle componenti. In neutronica questo corrisponde al tenere in considerazione il fatto che i neutroni, in un sistema critico, hanno importanza differente a seconda della posizione del reattore in cui si trovano, ed in particolare risulteranno in generale pi` u importanti nelle zone ove essi sono pi` u numerosi. Sarebbe dunque in teoria pi` u corretto usare l’importanza come funzione peso, ma questo richiederebbe la soluzione del problema aggiunto, aumentando cos`ı il costo computazionale in modo non giustificato dai risultati ottenuti. Solo nel caso monodimensionale, poichè l’importanza coincide con il flusso neutronico, pesando secondo Galerkin si ottiene una pesatura tramite importanza neutronica. Veniamo all’analisi di un caso particolarmente sfortunato. Per dare il via all’algoritmo iterativo sar` a necessario imporre una soluzione di tentativo, vale a dire il termine q. Cosa accade se tale termine è, sfortunatamente, perpendicolare all’autovettore fondamentale? Si avrà in questo caso infatti

36

che: q=

n X

ch ϕh

h=2

che non ha componente lungo l’autovettore fondamentale e dunque non potra mai convergere ad esso. Il metodo delle potenze convergerà in questo caso verso il secondo autovalore di modulo massimo. Questo sar` a vero per` o solo dal punto di vista matematico. Dal punto di vista numerico infatti la perpendicolarità perfetta non esiste (non esiste infatti il concetto di zero, ma al massimo delle quantità molto piccole), ed infatti baster` a una iterazione affinchè gli errori di arrotondamento del calcolatore portino ad ottenere come M q un vettore non perfettamente perpendicolare a ϕ1 . Si avr` a in questo caso tuttavia una convergenza molto lenta, caso non favorevole ma comunque migliore della convergenza ad un valore sbagliato. Attenzione però che il rischio esiste comunque: troncando il processo dopo un numero insufficiente di iterazioni, il vettore soluzione sar` a molto pi` u vicino a ϕ2 che a ϕ1 , con le dovute conseguenze. Si ricorda ora per` o che, all’inizio del paragrafo, si è detto che il metodo delle potenze è utilizzato per risolvere ogni tipo di problema agli autovalori, mentre per quanto visto finora si è visto solo come risolvere problemi agli autovalori classici ma non. Per poter estendere tale metodo anche ai problemi generalizzati è necessario essere in grado di provare l’esistenza della matrice B −1 , ove B è la matrice del problema agli autovalori generalizzato: Ax = λBx Come risolvere a questo punto il problema? Sarà necessario ricondursi al caso precedente, andando a scrivere dunque: B −1 Ax = λx ed applicando il metodo delle potenze nella maniera seguente: B −1 Aq (0) = q (1) Ecco dunque che il problema si è complicato notevolmente. Si avranno due cicli di iterazioni concatenati l’uno dentro l’altro: • Un ciclo di iterazioni esterne: B −1 Aq (0) = q (1)

→

B −1 Aq (1) = q (2)

• Un ciclo di iterazioni interne: Aq (0) = Bq (1)

37

Ritrovandosi un processo iterativo interno, anche il valore di q (1) non è esatto ma approssimato rispetto all’aggiornamento della soluzione. Si potrebbe a questo punto pensare di fare un’operazione di questo tipo: calcolare all’inzio, una volta per tutte, l’inversa di B e utilizzarla per tutte le iterazioni esterne. Ci si libererebbe cos`ı di un processo iterativo tramite un quantitativo di calcoli che, seppure elevato, sarebbe presumibilmente inferiore alla somma di tutte le iterazioni interne. Per capire perchè tale procedimento spesso non è attuato bisogna andare ad analizzare il contenuto fisico del problema. Se dal punto di vista matematico il discorso è teoricamente corretto, considerando le implicazioni fisiche ci si accorge che si andrebbe a risolvere i due problemi su due piani diversi, e dunque a mischiare informazioni di qualità molto differenti. Perdendo il parallelismo fisico, la convergenza non è pi` u garantita. Si vedrà in seguito come questa problematica appare in maniera molto evidente nel caso della risoluzione numerica delle equazioni della diffusione neutronica multigruppo. Questo discorso si attua però solo ove la matrice B sia sparsa, e dunque risolvibile tramite algoritmi di tipo iterativo. Ove invece B sia una matrice densa verr` a utilizzato, come già detto in precedenza, un metodo diretto, che porter` a forzatamente all’ottenimento della soluzione esatta.

2.4

Propriet` a spettrali delle matrici di interazione

Si definisce matrice di interazione una matrice i cui coefficienti hanno un significato fisico. Si enunciano di seguito alcune proprietà e teoremi di matrici che risulteranno utili nel seguito: • Se la matrice M è reale e simmetrica o complessa ma Hermitiana, allora tutti gli autovalori di M sono reali e tutti gli autovettori corrispondenti sono tra loro perpendicolari • Se la matrice M è reale e simmetrica e definita positiva14 , allora tutti i suoi autovalori sono positivi e reali • Se M è reale e simmetrica e semi-definita positiva15 , allora tutti i suoi autovalori saranno non negativi Si definicono invece matrici positive (Da notare la differenza con le matrici definite positive) matrici i cui elementi sono tutti positivi. Riguardo le matrici positive esiste un teorema di rilevante importanza: Teorema 7 (di Perr` on) Data una matrice M positiva, anche non simmetrica, vale che: 14 15

~ 6= ~0, si ha che (< M Φ, ~ Φ >) > 0 una matrice è definita positiva se, ∀ Φ ~ ~ Φ >) ≥ 0 ~ una matrice è semi-definita positiva se, ∀ Φ 6= 0, si ha che (< M Φ,

38

1. l’autovalore fondamentale è positivo ed ha sempre molteplicit` a 1 (ovvero ∃λi con m(λi ) = 1 tale che λ1 > |λj | ∀j 6= 1) 2. All’autovalore fondamentale λ1 corrisponde, a meno di una costante moltiplicativa arbitraria, un unico autovettore ϕ1 avente componenti tutte positive16 . L’autovettore fondamentale è l’unico a godere di tale propriet` a. 3. Se un qualsiasi elemento della matrice M cresce, allora λ∗1 autovalore fondamentale della matrice M ∗ avr` a modulo maggiore di λ1 Possiamo allo stesso modo definire matrici ad elementi non negativi o matrici semipositive le matrici aventi solo termini positivi o nulli. Per esse vale invece il seguente teorema: Teorema 8 (di Frobenius) Data una matrice semipositiva valgono le propriet` a enunciate dal teorema 7 a patto che la matrice sia anche irriducibile. Si faccia per` o attenzione che la condizione espressa dal punto 1 del teorema di Perron include una maggiorazione semplice e non stretta: λ1 ≥ |λj |.

2.5

Approssimazione di dati sperimentali

Esistono numerosi modi per approssimare dei dati sperimentali. Uno dei pi` u comuni è quello dell’interpolazione, ove essendo in possesso di n punti sperimentali potremo approssimarli con un polinomio di grado n-1. Purtroppo la determinazione dei coefficienti di tali polinomi costituisce un sistema di equazioni lineare e non omogeneo Altra via per l’approssimazione dei dati è il metodo dei minimi quadrati, dal quale non si otterr` a una funzione interpolante. Tale sistema è generalmente utilizzato ove sia presente una quantità molto elevata di dati sperimentali, caso in cui non è pensabile il ricorso ad una interpolazione. Si tenga conto comunque che nel caso della determinazione dei coefficienti della funzione approssimata tramite il metodo dei minimi quadrati il problema è tanto pi` u mal condizionato quanto pi` u grande è il numero n di punti sperimentali a disposizione. Useremo in questi casi i polinomi di Hermite e le funzioni spline.

2.5.1

Funzioni spline

Le funzioni spline sono dei metodi utilizzati per l’interpolazione di punti. Abbiamo infatti visto come una semplice interpolazione polinomiale spesso non è adeguata ed è dunque opportuno ricorrere ad altri metodi. 16

ovviamente a meno di una costante moltiplicativa arbitraria. Questo corrisponde a dire che hanno tutte lo stesso segno

39

Supponiamo dunque di avere N+1 punti, ciascuno identificato da una coppia (xi , yi ), compresi in un intervallo [a,b] tale che a ≡ x0 < x1 < · · · < xn ≡ b In linea teorica è sempre possibile, come detto, trovare un polinomio di grado N in grado di interpolare esattamente gli N+1 punti sperimentali, i cui coefficienti sono determinabili risolvendo il sistema di equazioni avente come matrice di coefficienti una matrice di Van der Monde. In effetti però, oltre al fatto che tale problema è spesso e volentieri mal condizionato, tale tipo di approssimazione non è necessariamente la pi` u adatta. Un alternativa pu` o essere quella di usare una serie di funzioni lineari, il che mi fornisce di fatto una linea spezzata. La funzione interpolante sarà in questo modo del tipo: fn (x) =

(xi+1 − x)yi + (x − xi )yi+1 xi+1 − xi

per xi < x
Una funzione di questo tipo è semplice da implementare, poco costosa, ma ha un grosso problema: la sua derivata non è continua. Prendiamo ora un altro esempio, quello in cui la soluzione interpolante è un arco di parabola: prenderò cos`ı i nodi a tre a tre invece che a due a due, ottenendo per` o in generale una rappresentazione poco realistica e continuando ad avere la discontinuità nella derivata, anche se questa volta a nodi alterni. Definiamo ora le funzioni interpolanti di tipo spline. Diremo che una funzione Sd (x), con d ≥ 1 è una funzione spline di ordine d associata ai punti {xi , yi }su di un intervallo [a,b] se: • Sd (x)è un polinomio di grado d in ogni intervallo [xi , xi+1 ] • Ogni derivata k-esima per k ≤ d − 1 è continua Le funzioni che vedremo in quanto in genere maggiormente usate sono le spline di ordine 3, per le quali avremo discontinuità a partire dalla derivata terza. Nelle spline cubiche avremo che  per i = 0 : n  S3 (xi ) = yi S3 (x) = ai + bi x + ci x2 + di x3 ∀x ∈ [xi−1 , xi ] per i = 1 : n  (k)+ (k)+ S3 (xi ) = S3 (xi ) per i = 1 : n − 1, k = 0, 1, 2 Ove di fatto le prime condizioni corrispondono al passaggio della funzione per i nodi (n+1 condizioni), le seconde all’imporre che la funzione interpolante sia un polinomio di terzo grado (4n incognite) e la terza alla continuità della funzione e delle sue derivate prima e seconda (3n-3 condizioni). Avremo dunque in totale 4n incognite e 4n-2 condizioni. Rimangono dunque due

40

condizioni da imporre a piacere, che saranno generalmente ad esempio l’imposizione dei valori della derivata della funzione agli estremi del dominio di integrazione. Si tratta ora dunque apparentemente di risolvere un sistema abbastanza impegnativo di 4n equazioni ed incognite. In realtà tramite opportune trasformazioni è possibile ridurre tali dimensioni sensibilmente Definiamo una nuova incognita: Mi = S300 (xi ). Essendo la mia spline un polinomio di terzo grado, il polinomio risultante da una derivata seconda sar` a lineare. Dunque avremo che: S300 (x) =

(xi − x)Mi−1 + (x − xi−1 )Mi xi − xi−1

Potr` o dunque scrivere la mia funzione integrando due volte ed ottenendo cos`ı S3 (x) =

(xi − x)3 Mi−1 + (x − xi−1 )3 Mi + Ci (x − xi−1 ) + Di 6(xi − xi−1 )

Imponendo la condizione di passaggio per i punti sperimentali ottengo ( i−1 i−1 − hi Mi −M Ci = yi −y hi 6 Di = yi−1 −

h2i 6 Mi−1

Da cui l’equazione diventa:

S3 (x) =

(xi −x)3 Mi−1 +(x−xi−1 )3 Mi + 6hi yi −yi−1 Mi −Mi−1 ](x +[ hi − hi 6

− xi−1 ) + yi−1 −

h2i 6 Mi−1

Non resta altro che definire i valori degli Mi , che si ottengono imponendo la continuit` a della derivata prima. Si vede che apparentemente due condizioni non sono state utilizzate. In realt` a la continuit` a della derivata 0-esima (ovvero della funzione interpolante stessa) è assicurata dall’imposizione del passaggio per gli n+1 punti, mentre la continuit` a della derivata seconda negli stessi nodi è assicurata dal fatto che essa stessa è stata scelta come incognita. Imponiamo dunque la continuità della derivata prima: (xi −x)2 Mi−1 +(x−xi−1 )2 Mi i−1 i−1 + [ yi −y − hi Mi −M ]= 2hi hi 6 (xi+1 −x)2 Mi +(x−xi )2 Mi+1 yi+1 −yi Mi+1 −Mi = + hi+1 − hi 2hi+1 6

Semplificando si ottiene: hi Mi−1 + 2(hi + hi+1 )Mi + hi+1 Mi+1 = 41

6 hi+1

(yi+1 − yi ) −

6 (yi − yi−1 ) hi

Il sistema che ne deriva è di tipo tridiagonale simmetrico e dunque sostanzialmente facile da risolvere tramite il metodo della doppia passata. Si noti inoltre come la matrice del sistema sia a diagonale dominante, il che assicura la convergenza dell’algoritmo. Vediamo ora quali sono le due condizioni aggiuntive che possiamo imporre. In generale ci saranno 3 possibili strategia adottabili: 1. Imporre il valore della derivata prima agli estremi 2. Imporre il valore della derivata seconda agli estremi. In questo caso, viste le incognite che sono state scelte, questo è il sistema che permette la risoluzione pi` u semplice in quanto consiste nell’imporre i valori di due delle incognite a priori. Si parla in questo caso di spline naturali. Tale condizione è ad esempio impostata di default nel software “Matlab”, ma si faccia attenzione che non sempre è fisicamente rappresentativa 3. Nel caso di andamento periodico si impone l’uguaglianza agli estremi delle derivate prime e seconde. In questo caso tuttavia si perde la tridiagonalit` a della matrice

2.6

Integrazione numerica

Come si calcola un integrale in modo numerico? Avremo Z b I= f (x)dx a

Il modo pi` u semplice per tale calcolo è dato dalle formule di NewtonCotes, ove scelti dei nodi xi vengono calcolati dei pesi ωi tali per cui il mio integrale è approssimato da una sommatoria del tipo: I=

n X

f (xi )ωi

i=1

L’entit` a dell’errore generato dall’approssimazione dipende dal numero n di nodi. Queste sono formule di quadratura di tipo interpolatorio polinomiale, di conseguenza un numero n di nodi corrisponderà ad un polinomio interpolante di grado n-1 che passi esattamente attraverso tutti questi punti. Se i nodi sono equidistanti, come è nel caso delle formule di Newton-Cotes, tale grado di precisione è pari a (n-1). Tale informazione corrisponde di fatto a dire che tale metodo mi consente di integrare con esattezza equazioni di grado pari o inferiore ad (n-1). L’obbiettivo di studi avanzati nel calcolo numerico è quello di aumentare tale grado di precisione. Le formule Newton-Cotes sono infatti ormai obsolete. L’aumento di precisione tuttavia, come sempre, non è a costo nullo. 42

Per applicare un metodo numerico è necessario di fatto sostituire la funzione integranda con un’altra funzione pi` u facile da integrare. Per fare questo, è necessario scegliere tale funzione ed applicarvi le opportune condizioni al contorno per il calcolo dei coefficienti. In particolare nel caso delle formule Newton-Cotes tale operazione era effettuata imponendo n condizioni nel calcolo dei pesi. Tali condizioni sono del tipo: Z

b k

x dx = a

n X

ωi xki

i=1

Questo corrisponde a pretendere che i miei pesi mi permettano di interpolare esattamente tutte le funzioni di base dello spazio vettoriale dei polinomi di grado n-1. Ho in questo modo scritto una matrice di tipo Van der monde, ovvero del tipo:      b−a 1 1 ··· 1 ω1   x1   b2 −a2  x2 · · · xn    ω2     2  ..  ..   ..  =  .. ..  .  . .  .   . ··· xn−1 xn−1 2 1

xn−1 n

ωn

bn −an n

Si vede dunque che il problema si riduce alla risoluzione di una matrice di Van der Monde, che sappiamo essere fortemente mal condizionata per n elevati. In particolare tale sistema è generalmente utilizzato per n ≤ 4. Tale formula di quadratura ha dunque potenzialità limitate. Veniamo ora ad ulteriori considerazioni. Supponiamo di dover calcolare l’integrale Z b I= ω(x)f (x)dx a

Questa volta ho due funzioni. In realtà sono funzioni create in maniera tale da avere: • una parte (ω(x)) ove ho posto tutte le singolarità della funzione integranda, ma anche tutto il contenuto fisico del problema. • una parte (f(x)) regolare Stiamo dunque di fatto trattando del caso visto inizialmente, ove a p corrisponde K ed a f corrisponde ϕ. L’idea è per` o quella di vedere come effettuare un upgrade delle formule Newton-Cotes, ovvero di ottenere un grado di precisione nell’approssimazione dell’integrale superiore a (n-1). Una possibilit` a è quella di abbandonare l’imposizione a priori dei nodi. Questi dunque non soltanto saranno incognite di equazioni, ma non saranno pi` u evidentemente equidistanti. Si tratta dunque a questo punto di imporre non pi` u soltanto n condizioni per il calcolo dei pesi, ma anche ulteriori n 43

condizioni per il posizionamento dei nodi. Otterrò dunque un sistema del tipo: Z b n X k ω(x)x dx = λi (µi )k per k = 0...2n − 1 a

i=1

Ove i λi sono i pesi ed i µi sono i nodi. Il problema, che da un lato fornirà generalmente una soluzione pi` u precisa, è però diventato pi` u complesso: non soltanto le condizioni sono diventate 2n, e dunque è raddoppiata la dimensione del sistema da risolvere, ma è diventato non lineare. In questo diventa dunque fondamentale l’utilizzo di polinomi ortogonali, ai quali sono associate formule di quadratura dette gaussiane L’obbiettivo che ci si pone, aumentando il numero dei gradi di libertà, è quello di ottenere un gradi di precisione maggiore, in particolare passando da n-1 a 2n-1, ovvero essere in grado di integrare in maniera esatta polinomi di grado pari o minore Dunque andrò ad imporre tale condizione R b a 2n-1. k che, chiamati mk = a ω(x)x dx, permette di ottenere il seguente sistema di equazioni:  λ1 + λ2 + · · · + λn = m0     λ1 µ1 + λ2 µ2 + · · · + λn µn = m1 ..  .    + λ2 µ22n−1 + · · · + λn µn2n−1 = m2n−1 λ1 µ2n−1 1 Ho cos`ı imposto 2n condizioni. A supporto della funzionalità di questi sistemi esistono alcuni teoremi che enunceremo qui di seguito: Teorema 9 Nell’ipotesi in cui ω(x) ≥ 0 ∀x ∈ [a, b] , ovvero sempre non negativo, e in cui tutti gli mk esistano per ogni k tra 0 e 2n-1, il sistema di equazioni non lineari che ne deriva ha una soluzione unica ed univoca, ovvero esiste uno ed un solo set di {λi } e {µi } che forniscano un polinomio interpolante di grado di precisione 2n-1 Si noti che procedendo in questa maniera i nodi non sono pi` u equidistanti tra loro. Questo diventa di importanza molto rilevante ove essi vadano a corrispondere con la discretizzazione di un dominio. Tali formule di quadratura dovrebbero per` o, quantomeno in linea teorica, attuare un processo di discretizzazione “intelligente”, ovvero tale da addensare la discretizzazione ove la funzione ha variazioni pi` u sensibili e a renderla meno fitta ove invece la funzione sia pi` u regolare. Il problema è che ora la matrice da risolvere è ancora pi` u complessa e, comunque, mal condizionata. A questo proposito entra in gioco un teorema fondamentale per l’integrazione numerica: Teorema 10 Affinchè la formula di quadratura sia gaussiana, ovvero abbia grado di precisione pari a 2n-1, deve essere che L’insieme dei µi coincida 44

con l’insieme degli zeri del polinomio Pn (x) ortogonale in [a, b] rispetto alla funzione peso ω(x) Definiamo ora meglio i polinomi ortogonali Prendiamo un intervallo [a,b], non necessariamente finito, sul quale definiamo una funzione peso ω(x) che sia a segno costante su tutto l’intervallo e della quale esistano i momenti, definiti come: Z Mk =

b

ω(x)xk dx < ∞

per k = 1, 2...

a

Individuo quindi una classe di polinomi {P0 (x), P1 (x)....P( x)} che sarà detta ortogonale su [a,b] rispetto alla funzione peso ω(x) se Z

b

ω(x)Pn (x)Pm (x)dx =

a

0 n 6= m cost n = m

Tale costante sar` a uguale per tutti gli n, ed un set di polinomi è definito ortonormale ove tale costante sia uguale ad 1. In ogni caso tale costante sar` a sempre positiva in quanto abbiamo supposto la non negatività della funzione peso p(x) su tutto l’intervallo [a,b]. La scelta dell’intervallo e della funzione peso individua inequivocabilmente una classe di polinomi ortogonali a meno di coefficienti costanti e non nulli. In generale i polinomi ortogonali sono definiti in base ad una formula di ricorrenza a tre termini, del tipo:   P−1 (x) = 0 P0 (x) = 1  Pn+1 (x) = (x − an )Pn (x) − bn Pn−1 (x) Si vede dunque che una volta definiti l’intervallo [a,b], la funzione peso e i valori dei coefficienti an e bn tutti i polinomi ortogonali sono determinabili in maniera univoca. Si vede inoltre come dentro la definizione di un generico Pn+1 siano coinvolti tutti i polinomi di ordine inferiore. L’interesse che risiede nei polinomi ortogonali è dato dalle proprietà dei loro zeri. Si ha infatti che: • Per ogni n ≥ 1, Pn (x) possiede n zeri reali e distinti tutti contenuti nell’intervallo [a,b] • Tra due zeri consecutivi di Pn (x) ne è incluso uno solo di Pn−1 (x) Prendiamo ora un generico polinomio Qm (x). Esso sarà rappresentabile in modo univoco tramite una combinazione lineare di polinomi ortogonali: Qm (x) =

m X i=0

45

ck Pk (x)

Per il calcolo dei coefficienti posso semplicemente sfruttare l’ortogonalità dei polinomi: per il calcolo di un generico ck è sufficiente moltiplicare ambo i membri per il corrispondente polinomio Pk (x) e per la funzione peso ed integrare da entrambe le parti su [a,b]. Se ne ottiene: Rb ω(x)Qm (x)Pk (x)dx ck = a R b 2 a ω(x)Pk (x)dx Dalla propriet` a di ortogonalità dei polinomi si può anche ricavare che: Z b ω(x)Pn (x)qm (x)dx = 0 ∀m < n a

Questo è dovuto al fatto che io posso sempre scrivere Qm (x) =

m X

ck Pk (x)

i=0

Ma essendo ogni volta i diverso da n, per l’ortogonalità dei polinomi l’integrale è sempre nullo. Vediamo ora alcune classi di polinomi ortogonali: I Polinomi di Legendre hanno • ω(x) = 1 • [a, b] = [−1, 1] Il fatto che la funzione peso abbia valore unitario mi dice che tutti gli intervalli in cui posso suddividere la mia funzione originaria hanno la stessa importanza. I polinomi di Legendre sono definiti dalla seguente formula di ricorrenza a tre termini:   P0 (x) = 1 P1 (x) = x  (n + 1)Pn+1 (x) = (2n + 1)xPn (x) − nPn−1 (x) Sono dei polinomi particolarmente importanti perchè sono molto adatti nella discretizzazione di particelle libere in geometria sferica, soprattutto ove si voglia discretizzare la variabile angolare. Un esempio è l’equazione del ~ trasporto, ove la velocit` a di un neutrone sarà data in generale da ~v = v Ω. Gli zeri dei polinomi di Legendre corrisponderanno alle direzioni di volo delle particelle. I Polinomi di Jacobi sono un’altra categoria di polinomi ortogonali, aventi: • ω(x) = (1 − x)α (1 + x)β • [a, b] = [−1, 1] 46

Si vede evidentemente che per α = β = 0 la scrittura degenera nei polinomi di Legendre, che costituiscono dunque un caso particolare di polinomi di Jacobi. Per polinomi di Jacobi con α = β = − 12 ho i Polinomi ti Tchebychev di prima specie, che vedremo in seguito utilizzati come supporto per l’accelerazione della convergenza di algoritmi lenti. Per invece α = β = 12 si hanno i polinomi di Tchebychev di seconda specie che servono ad approssimare funzioni speciali, come la E1 . I Polinomi di Laguerre hanno invece: • ω(x) = e−x • [a, b] = [0, +∞[ e sono usati in problemi di rilevazione satellitare e di fenomeni radiativi (si noti a questo proposito la presenza del termine esponenziale) I Polinomi di Hermite, hanno infine: • ω = e−x

2

• [a, b] =] − ∞, +∞[ Dunque ci` o che ci interessa è trovare gli zeri di tali polinomi. Per farlo conosciamo, di base, almeno quattro modi diversi: • Bisezione • Regula falsi • Secanti • Newton In realt` a non si andranno mai ad usare tali metodi, ma ne esistono di appositi per i polinomi ortogonali. Dunque dal problema Z b X ω(x)f (x)dx = λi f (µi ) + err a

i

sono passato prima a:  λ1 + λ2 + · · · + λn = m0     λ1 µ1 + λ2 µ2 + · · · + λn µn = m1 ..  .    λ1 µ2n−1 + λ2 µ22n−1 + · · · + λn µn2n−1 = m2n−1 1 ed infine a: Pn (µi ) = 0 47

Ove Pn è il polinomio ortogonale scelto all’inizio del problema in funzione dell’intervallo di integrazione e della funzione peso. La risoluzione del problema intermedio è dunque inutile, visto che le soluzioni sono date dagli zeri del polinomio. Dunque di fatto posso ottenere i valori dei nodi, grazie ai quali andr` o a calcolare i pesi. Andiamo per esempio a vedere cosa succede utilizzando i polinomi di Legendre: otterremo una equazione nella forma: Z 1 X 1 · f (x)dx ' λi f (µi ) −1

i

Si vede che in questo caso la funzione peso è pari ad uno, e dunque di fatto non d` a alcun peso differente ai vari intervalli. Si parla in questo caso di formule di Gauss-Legendre Stesso discorso pu` o essere effettuato per altri sistemi, ottenendo ad esempio le formule di Gauss-Laguerre: Z ∞ X e−x f (x)dx ' λi f (µi ) 0

i

ove i nodi sono in questo caso dati dagli zeri del polinomio di Laguerre. I risultati portano a situazioni differenti. Nel caso delle formule di GaussLegendre ottengo un addensamento dei nodi ai bordi del dominio, con tendenza tanto pi` u marcata quanto pi` u n è grande (nonostante la funzione peso sia costante). Attenzione dunque. La scelta della formula di quadratura non è casuale, ma effettuata in base al legame fisico ed alle esigenze computazionali17 Ho dunque trasformato il mio problema in un problema di ricerca di zeri di polinomi omogenei. Attenzione però: tale problema è molto mal condizionato, soprattutto ove la derivata prima del polinomio sia piccola e dove il numero di zeri n sia elevato. Maledizione! Dunque neanche questo terzo problema numerico sarà quello che andremo, in definitiva, a risolvere. Dovremo di fatto ripensare il 17

Ritornando a quanto detto sulla propagazione degli errori, le formule di ricorrenza a tre termini sono esattamente uno dei punti potenzialmente problematici di un algoritmo dei quali andranno dunque studiati condizionamento e stabilita. Dunque in effetti quando andremo ad inserire una formula di ricorrenza a tre termini, non lo faremo secondo la definizione vista in precedenza, ma tramite l’algoritmo di Clanshaw, che trasforma quello precedente in uno pi` u stabile:   ym+2 = ym+1 = 0 yk = (Ak x ¯ + Bk )yk+1 − Ck yk+2 + ck  y0 k0 = Qm (¯ x) Si tratta di un algoritmo ricorsivo all’indietro, ove devono essere noti i coefficienti ricorsivi del polinomio ortogonale. Il risultato che si ottiene utilizzato questo tipo di algoritmi è affidabile

48

problema numerico dall’inizio, in modo tale da ottenere le stesse soluzioni del problema che ci siamo posti, rendendolo però pi` u facile da risolvere. La prima operazione da effettuare è una normalizzazione dei polinomi, in modo tale che Z b

2 ω(x)Pm (x)dx = 1

a

Ottenendo cos`ı un insieme di polinomi ortonormali. In questo modo la formula di ricorrenza assume la forma: Pm+1 (x) = (Am x + Bm )Pm (x) −

Am Pm−1 Am−1

Che si pu` o riscrivere nella forma xPm (x) =

1 Bm 1 Pm+1 − Pm (x) + Pm−1 Am Am Am−1

che, rinominando opportunamente i coefficienti come: αm =

1 Am

e

βm = −

Bm Am

diventa xPm (x) = αm Pm+1 + βm Pm (x) + αm−1 Pm−1 Posso dunque riscrivere il tutto in forma matriciale, del tipo: xP~ (x) = TˆP~ (x) + αn−1 Pn (x)en

con

en αn−1 = {0, 0, · · · , 1}

Se calcoliamo tale espressione nei nodi µi , sapendo che questi sono tali da essere zeri di Pn , otteniamo la forma finale: µi P~ (µi ) = TˆP~ (µi ) Questo non è altro che un problema agli autovalori di tipo classico ove gli zeri del polinomio Pn (x) coincidono con gli autovalori del problema stesso. Questo è dunque, in definitiva, il vero problema che andremo a risolvere. La matrice T è in effetti tridiagonale e simmetrica, il che come detto rende possibile ottenere la risoluzione tramite l’applicazione del metodo della doppia passata. I pesi sono infine determinati, una volta noti gli autovalori del problema e di conseguenza i nodi µi , dall’equazione: 1 2 j=0 [Pj (µi ) ]

λi = Pn−1

49

2.7

Equazioni integrali

Ci si riferisce in particolare alle cosiddette equazioni di Frehdolm di seconda specie. Non è ovviamente l’unico tipo esistente, ma ve ne sono altre come le equazioni integrali di Volterra. Caratteristica peculiare delle equazioni di Frehdolm è quella di basarsi su integrali definiti. Un tipico esempio di equazione di Frehdolm di seconda specie è il seguente: Z ϕ(x) = h(x) + K(x, y)ϕ(y)dy Con riferimento alla simbologia dell’equazione precedente, abbiamo che: ϕ(x) è la funzione incognita h(x) è il termine forzante K(x, y) è il nucleo o kernel della funzione di integrazione, ove in pratica è situata la fisica del problema. 18 Ove il termine forzante sia assente, l’equazione di partenza viene ricondotta (come visto nel paragrafo precedente) ad una equazione integrale agli autovalori, del tipo: Ax = λx

oppure

Ax = λBx

tali casi si riconducono entrambi alla risoluzione di un problema del tipo D(λ)x = ~0, ove D(λ) è una matrice dipendente da parametro. Il caso pi` u agevole, ma anche in geneale quello meno interessante, è quello ove tale dipendenza sia lineare. Quelli enunciati fino ad ora sono problemi matematici. Per passare al problema numerico sar` a necessario introdurre degli errori di idealizzazione. Partiamo dal kernel, ove come detto è nascosta la fisica del problema. Sulla base di questo dovremo procedere all’integrazione scegliendo opportunamente la formula di quadratura. Una volta effettuata tale operazione la risoluzione dell’equazione integrale diventerà un problema algebrico del tipo: n X ϕ(x) ' h(x) + ωi k(x, yi )ϕ(yi ) i=1

Tuttavia è necessario abbandonare l’idea di calcolare direttamente la soluzione esatta ϕ(x). Ci` o che riuscirò ad ottenere sarà sempre e comunque una soluzione approssimata. Andrò dunque a sostituire alla mia ϕ(x) una ϕn (x) dove l’indice n mi d` a informazioni sulla precisione del polinomio interpolante. 18 Tale nucleo pu` o essere singolare nel dominio di integrazione; in questo caso la risoluzione diventa maggiormente complessa

50

La scelta di n dipende dal problema fisico e, come sempre, dal rapporto costi-benefici. Una volta effettuato tale passaggio otterrò che la mia equazione sarà divenuta della forma: ϕn (x) ' h(x) +

n X

ωi k(x, yi )ϕn (yi )

i=1

Occorre a questo punto ricordare che, nel calcolo numerico, non si opera quasi mai nel continuo, ma si opera un processo di discretizzazione per arrivare ad un numero finito di nodi in cui sarà calcolata la funzione. Si tratter` a dunque di andare a collocare i nodi, che sono i punti in cui imporrò l’esattezza della mia equazione. Otterrò cos`ı una forma del tipo: ϕn (xi ) = h(xi ) +

n X

ωi k(xi , yi )ϕn (yi )

i=1

In questo caso il “circa uguale” è stato sostituito con un uguale in quanto stiamo imponendo proprio che la nostra approssimazione sia esatta nei nodi scelti. Sono dunque di fronte ad un problema lineare non omogeneo del tipo X [δij + ωi k(xi , yi )]ϕn (xi ) = h(xi ) i

che assume di fatto una forma del tipo Ax=b. Trattandosi di un problema non omogeneo esso sar` a risolvibile solo se det(A) 6= 0, e come qualsiasi problema numerico è necessario accertarsi del buon condizionamento della matrice A. Una volta effettuata tale verifica, esisteranno come sempre le due strade: • se la matrice è sparsa, si opererà iterativamente • se la matrice è densa, si opererà tramite metodi diretti In generale sar` a pi` u probabile, nel caso di equazioni integrali, trovarsi di fronte a matrici dense. Un esempio tipico di questo caso è la trasmissione dei segnali, ove nonostante si abbia una attenuazione esponenziale del segnale è teoricamente necessario tenere conto di ogni possibile sorgente, il che porta ad avere matrici quasi interamente non nulle.

51

Capitolo 3

Metodi numerici per l’ingegneria nucleare 3.1

Introduzione alla tecnologia ed alla fisica nucleare

Un reattore nucleare è identificato da un punto di funzionamento. Possiamo supporre che esistono varie zone: • Sicure • di Allarme • Incidentali Il produttore deve tendenzialmente adattarsi alle richieste di carico della rete. Questo dipende tuttavia anche da numerosi altri fattori: in Italia ci si aspetta che eventuali centrali nucleari di nuova costruzione opererebbero in modo tale da supportare il carico di base, lavorando quindi in condizioni teoricamente ottimali. In un contesto pi` u generale tuttavia anche i gestori nucleari dovranno adattarsi alle richieste di rete, il che li porterà a porsi in punti di funzionamento differenti dall’ottimale. Il punto chiave è che è necessario assicurarsi che tali variazioni del punto di funzionamento non portino a spostarsi al di fuori della zona sicura. Si noti che pi` u i miei codici di controllo sono efficaci e precisi, pi` u la zona definita sicura potrà essere estesa. Dunque, come è accaduto di recente in Svezia, a volte per aumentare la potenza prodotta dalle centrli è sufficiente migliorarne il controllo in modo da potersi spingere in sicurezza in zone di funzionamento prima inaccessibili. Nelle pagine che seguono ci si concentrerà sullo studio di piccoli transitori e piccole oscillazioni attorno ad un valore costante. Si noti che, nonostante tali fenomeni verranno visti in riferimento al nucleare, non sono affatto di esso esclusivi, ma anzi esistono altri casi di controllo simile1 . 1

Si vedano ad esempio i problemi legati al ripple di corrente

52

Ammettiamo inizialmente di avere un sistema molto semplice, di cui supporremo di conoscere con precisione la geometria (considerata invariante nel tempo, anche se questo è vero solo fino ad un certo punto) e la composizione materiale (che invece varia nel tempo in modo marcato). Supponiamo inoltre che: • Ogni neutrone uscito è perduto; non vi è dunque rientro di neutroni • Non ci sono neutroni entranti nel sistema dall’esterno: nè, come detto, riflessi nè generati da una sorgente esterno Queste ipotesi ci permettono di considerare un problema omogeneo Utilizzeremo come ultieriore astrazione quella di supporre che i neutroni siano gi` a presenti nel sistema, ovvero ci dimenticheremo dei transitori di avviamento. All’interno di questo sistema i neutroni interagiranno con il materiale nucleare (e non) presente e potranno avere determinate iterazioni piuttosto di altre a seconda della loro probabilità, quantificata dalla sezione d’urto. Sappiamo dalla fisica nucleare che l’interazione di un neutrone con un nucleo pu` o dare luogo alla scomparsa del neutrone stesso ed alla produzione di una trasmutazione del nucleo di partenza. Delle reazioni di questo tipo quella di maggior interesse sarà quella di fissione. In questo caso il neutrone, catturato da un nucleo fissile, porta alla formazione di un nucleo composto ed instabile che, in tempi trascurabili, si scinderà in pi` u nuclei di dimensioni inferiori provocando la produzione di radiazioni ed altri neutroni. Questo tipo di fenomeni ci porta a comprendere come, in effetti, il supporre composizione materiale costante sarebbe una approssimazione troppo forte. Non tutti i tipi di interazione tuttavia danno luogo a questo tipo di fenomeni, ed anzi quelli in generale pi` u frequenti saranno quelli in cui il neutrone rimbalza contro i nuclei del materiale in cui si muove, dando cos`ı luogo a scattering. Lo scattering è, per natura, un fenomeno isotropizzante; Tutte le sezioni d’urto che vedremo, invece, sono da considerarsi in materiale isotropo.2 Al fine del controllo del reattore, e lo ripeteremo pi` u volte, sono fondamentali i neutroni ritardati. I fenomeni di emissione rapida sono infatti troppo veloci per essere controllati, mentre l’incidenza dei fenomeni di emissione ritardata è, seppure piccola, sufficientemente elevata da consentire il controllo delle reazioni. Il nostro obbiettivo sarebbe quello di mantenere la potenza prodotta costante. Tuttavia, per quanto detto, un sistema nucleare per fare ciò avrà 2 La differenza tra isotropo e isotropizzante è sottile ma sostanziale. Lo scattering è un fenomeno isotropizzante nel senso che la direzione di volo del neutrone in uscita ha probabilit` a indipendente dall’angolo solido. Un materiale è invece detto isotropo rispetto ad un fenomeno ove tale fenomeno ha una probabilit` a di avvenire equa indipendentemente dall’angolo solido di volo

53

bisogno di un quantitativo elevato di sistemi di controllo e regolazione, che comunque saranno in grado di assicurare tale costanza di produzione solo in prima approssimazione. Esistono fondamentalmente tre tipi di sitemi di controllo: Barre di controllo : Sono sistemi a movimentazione elettromeccanica, di cui esistono a loro volta diverse tipologie: • Per lo spegnimento rapido (scram) • Per la macro regolazione (abbassamento / innalzamento richieste di rete) • Per la micro regolazione (o regolazione fine, sempre in moto per il controllo del k effettivo) Controllo chimico : Si parla in questo caso del cosiddetto avvelenamento del refrigerante per via chimica, in generale tramite dissoluzione di sali di boro Controllo fisico-chimico Alcune pastiglie inserite all’interno degli elementi di combustibile sono realizzate con dei materiali assorbitori (ad esempio ossido di gadolinio) Dei tre sistemi descritti l’ultimo è l’unico non online, ovvero che non può essere modificato a piacimento una volta che gli elementi di combustibile sono stati montati. Tale sistema è in generale usato per regolare la presenza di neutroni in determinate zone del reattore ed a determinate energie. I veleni di cui si è parlato finora sono bruciabili. Ciò implica che essi si consumano dopo l’uso. Dunque, in fase progettuale, la concentrazione degli ossidi di gadolinio nelle pastiglie deve essere calcolata tenendo in conto della progressiva diminuzione del loro quantitativo nel tempo. Il fenomeno per cui i sistemi di controllo non online sono pi` u efficaci all’inizio della vita del `. reattore che alla fine è detto rilascio di reattivita Si tenga conto comunque che i metodi numerici associati alla risoluzione dell’equazione della diffusione multigruppo non permettono nè il controllo di centrale nè la progettazione. Essi sono infatti applicati a problemi meno critici, ed in particolare alla programmazione della sostituzione delle barre di combustibile3 . 3 Ci si potrebbe in realt` a chiedere per quale ragione il consumo non pu` o essere programmato in anticipo sfruttando i metodi pi` u precisi che vengono utilizzati in fase di progettazione. Questo non accade in quanto durante la vita di un reattore, esso non opera sempre in condizioni nominali, ma seguendo una serie di fasi a potenze differenti ed andr` a incontro ad un largo numero di transitori, pi` u o meno traumatici, che renderanno impossibile la previsione a priori del consumo di combustibile. Questa dunque dovr` a essere programmata tenendo in considerazione, seppure in maniera approssimativa, la vita del reattore che precede la sostituzione. Si ricorda inoltre che transitori veloci sono molto pericolosi, in quanto possono provocare la formazione di cricche all’interno degli elemen-

54

3.1.1

Il concetto di criticit` a e le sue implicazioni

Come pu` o un sistema rimanere, da un istante all’altro, stabile? Quando il numero di neutroni che nascono da fissione per unità di volume e di tempo è sempre uguale a quelli che muoiono (sono assorbiti od escono dal sistema), per unit` a di volume e di tempo diremo che il sistema è critico.4 Ci stiamo dunque preoccupando dello studio di criticità di un sistema moltiplicante tramite le equazioni della diffusione multigruppo in 3D. Cercheremo in questa sede di analizzare lo sforzo computazionale che questo implica. Torniamo ora al concetto di criticità. Stiamo parlando di un sistema moltiplicante ove dobbiamo essere in grado di mantenere le reazioni a catena senza che queste diminuiscano od aumentino di intensità oltre certi livelli. Assumendo, come detto, composizione variabile, se il sistema è critico all’istante t non lo sar` a mai anche all’istante t + ∆t, a meno di non cambiare qualcosa nel sistema. Cosa succede quando porto questi ragionamenti in una logica computazionale? Otterr` o come noto dei risultati rappresentativi della realtà fisica, ma non esatti, e di questo dovr` o ovviamente tenerne in conto poichè si parla di sistemi molto sensibili agli errori, in particolare riguardo lo scostamento dalla criticit` a. Ogni volta che tale spostamento supera una soglia anche apparentemente molto piccola i problemi non possono pi` u essere considerati stazionari, e ci si sposter` a in un’altra classe di problemi, detti di cinetica neutronica. In teoria infatti ogni volta che ci troveremo in condizioni di kef f 6= 1 non potremo pi` u parlare di criticità e ci troveremo in condizioni cinetiche e, dunque, non stazionarie. All’interno degli studi di criticità è importante ritornare alla realtà anche dal punto di vista geometrico. Il nostro dominio non sarà pi` u dunque una sfera perfetta ma un sistema complesso che abbia il pi` u possibile le sembianze del reattore reale che vogliamo andare a studiare. Questo avrà anche implicazioni di tipo energetico legate, ad esempio, alla termodinamica della sottrazione del calore o al variare dell’effetto moderante in funzione della fase del moderatore. Come detto l’obbiettivo sarà quello di portarsi in condizioni di potenza termica costante. Nella realt` a delle cose non è tuttavia possibile raggiungere ti di combustibile che portano alla fuoriuscita dei prodotti di fissione e la conseguente contaminazione del circuito primario. Questo evento pu` o determinare la necessit` a della sostituzione anticipata di un elemento di combustibile. La criticit` a della guaina che circonda le barre d’uranio al fine di contenere i prodotti di fissione è tale che le saldature presenti in detta guaina sono controllate ai raggi X una per una in fase di produzione 4 La criticit` a di un reattore dipende esclusivamente dalla sua geometria e dalla sua composizione materiale in un determinato istante, e non dal numero di neutroni che vi sono contenuto. Quindi, almeno a livello concettuale, un reattore pu` o essere critico anche in totale assenza di neutroni, la cui presenza serve solo (e non è poco) alla produzione di energia. La presenza di neutroni, e le conseguenti iterazioni che essi hanno con il materiale, andranno istantaneamente a modificare la composizione materiale del sistema, provocando cos`ı la perdit` a di criticit` a dello stesso.

55

tale obiettivo, e il sistema tenderà di conseguenza ad oscillare continuamente tra condizioni di lieve (si spera) sotto e sopracriticità. La criticit` a è dunque raggiungibile solo agendo sui sistemi di controllo, che permettono di variare la composizione materiale del sistema a piacimento. Esistono invece sistemi di sicurezza intrinseca e passiva il cui compito sar` a quello di far raggiungere al sistema livelli fortemente sottocritici nel minor tempo possibile. In generale a sistema sovracritico corrisponde una crescita di potenza, mentre ad un sistema sottocritico corrisponde un calo della stessa. Questa corrispondenza, seppur frequente, non è tuttavia univoca. Una delle ragioni dell’incidente di Chernobyl è stata anche una condizione di apparente sottocriticit` a associata tuttavia ad un aumento di popolazione neutronica. ` sono situazioni in cui la criticità di un Gli incidenti di criticita sistema è raggiunta involontariamente, e sono evidentemente molto pericolose. Sono situazioni tipiche ad esempio dei processi di ritrattamento del combustibile, ove si pu` o per errore configurare un sistema (ad esempio di stoccaggio intermedio) in modo tale che questo risulti avere casualmente kef f > 1. 5 6

3.1.2

Fenomeni non stazionari

La semplificazione di pseudo-stazionarietà è accettabile solo in determinati casi, mentre in altri non è realistica. Alcuni esempi sono: Transitorie Si tratta di fasi non-stazionarie non pericolose, generate dall’attivit` a volontaria di regolazione derivante dalla gestione del sistema. Si parla in particolare delle fasi di • Avviamento • Spegnimento • Regolazione del carico Incidentali Si tratta in questo caso delle condizioni in cui è necessario spegnere rapidamente il reattore. Tale operazione è definita scram, ed è sensibilmente differente dallle operazioni di spegnimento “standard”. L’incidente pi` u grave che possa verificarsi in una centrale nucleare è il LOCA (Loss Of Cooling Accident), ove si ha una perdita immediata di praticamente tutto il liquido refrigerante. 5 L’ultimo incidente di questo tipo è avvenuto in Giappone nel 1999. Questi incidenti sono potenzialmente molto gravi e portano all’emissione di forti quantit` a di radiazioni γ e di neutroni, con possibilit` a anche di innesco di reazioni esplosive 6 Un caso molto particolare è quello di una zona mineraria del Gabon, ove nel tempo si sono naturalmente ripetute pi` u volte per pi` u di un milione di anni delle situazioni di criticit` a naturale. Si noti che l’uranio 235 è in generale presente in natura in quantitativi molto bassi, in media con arricchimenti inferiori allo 0,7%

56

Approfondiamo, per completezza, qualche cenno legato ai fenomeni transitori, ai quali sono infatti legate delle condizioni in cui il gestore di una centrale vorr` a portarsi deliberatamente in condizioni sovracritiche. Cosa succede ad esempio all’accensione? Supponiamo di aver appena spento il reattore per la ricarica e di volerne garantire l’attività continuativa per un quantitativo compreso tra 18 e 36 mesi. All’inizio di tale periodo il sistema, privo di barre di controllo, sarà dunque fortemente sovracritico per due ragioni: • Per poter garantire per il numero di mesi richiesti il mantenimento della criticit` a, anche in presenza di sistemi di controllo • Per poter accendere il sistema Supponiamo dunque di avere, in t=0, flusso neutronico nullo. Per accendere il sistema dovr` o, in qualche modo, buttarci dentro dei neutroni. Una volta completata questa fase dovrò poi utilizzare i sistemi di controllo per aumentare la potenza passo a passo fino al valore richiesto.7 Veniamo ora invece al caso della regolazione, ovvero ad esempio del caso in cui la rete richieda un aumento di potenza. Tale procedimento risulta essere analogo alle procedure di accensione e spegnimento, ove dunque il livello di potenza richiesto non è raggiunto tramite un unico passaggio ma attraverso diversi piccoli salti.

3.1.3

L’importanza neutronica ed i metodi per lo studio di problemi pseudo-stazionari

Immaginiamo di avere un sistema estremamente semplice, sul quale faremo le seguenti supposizioni: • sferico • critico • geometria e compisizione materiale ben note e costanti nel tempo • trascuriamo i problemi prettamente tecnici, come la strategia scelta per la sottrazione del calore e via dicendo • poniamo tale sistema sotto al Gran Sasso, ove dunque potremo pensare che non vi siano neutroni 7 teoricamente si dovrebbero avere problemi di questo tipo anche in fase di spegnimento o calo di potenza, ma in pratica tali operazioni sono molto meno delicate in quanto un eccesso di scostamento dalla criticit` a provocherebbe solo uno spegnimento troppo brusco e non una esplosione. Nonostante ci` o si tenga comunque in conto che per l’integrit` a del sistema pi` u i transitori sono lenti, meglio è

57

Supponiamo ora di poter generare all’interno di tale sistema delle sorgenti di neutroni puntiformi in maniera totalmente arbitraria. Chiamerò tali sorgenti S1 ed S2 , e le posizionerò in due sistemi altrimenti identici, in due posizioni diverse (immaginiamo S1 esattamente al centro del sistema ed S2 in prossimit` a del confine). La domanda che ci poniamo è la seguente: i due sistemi si comporteranno in maniera differente? Ricordando la supposizione iniziale di criticità e geometria e composizione materiale costanti, avrò che i due reattori, dopo il transitorio iniziale, si porteranno a due condizioni di equilibrio differenti tra loro, con potenze P1 e P2 prodotte diverse. In effetti il punto chiave sta nel fatto che i neutroni prodotti dalla seconda sorgente, essendo essa posizionata pi` u vicina alla frontiera, tenderanno con maggiore facilità ad uscire dal sistema rispetto a quelli prodotti dalla prima sorgente. Si vede dunque che, a livello neutronico, non tutte le zone del rettore sono uguali. Definiamo dunque importanza di un neutrone come il contributo che esso pu` o dare all’ottenimento della potenza asintotica del sistema. Essa è una funzione dell’energia del neutrone, della sua direzione di volo ma anche della sua posizione, come appena detto. A parit` a dunque delle prime due variabili, un neutrone pi` u vicino al cuore del sistema sarà generalmente pi` u importante di uno posto in prossimit` a della periferia8 . Conoscere l’importanza che i neutroni hanno in una determinata zona del reattore è fondamentale ai fini della progettazione delle barre di controllo ed in generale dei sistemi di sicurezza e regolazione. Il loro intervento sarà infatti tanto pi` u efficace tanto maggiore è l’importanza dei neutroni che essi andranno ad assorbire. Avremo dunque, ad esempio, che le barre di controllo per lo scram sono posizionate al centro del reattore. Torniamo al sistema sferico critico ideale. Cosa accadrebbe se andassimo a modificare la geometria o la composizione materiale? Non è difficile intuire che: • Se andiamo ad aumentare l’arricchimento, si avrà uno spostamento verso la sovracriticit` a del sistema • Se andiamo ad aumentare il raggio della sfera andremo a diminuire il rapporto superficie / volume e dunque a diminuire l’incidenza delle fughe sul bilancio neutronico, il che porta anche in questo caso ad uno spostamento verso la sovracriticità Se passiamo dunque al problema numerico corrispondente, avremo in generale due tipi di autovalori: • Materiale 8 In diffusione monodimensionale importanza e flusso neutronico coincidono. Questo mostra come, almeno come linea generale, che le zone ad importanza maggiore sono quelle ove vi sono pi` u neutroni

58

• Geometrico Questo mi dice che di fatto che per modificare la criticità del sistema potrò andare ad operare sia sulla geometria che sulla composizione geometrica dello stesso. Questa operazione non è tuttavia agevole. Torniamo ora alla buona vecchia analisi numerica. Il fine è quello di gestire i programmi di ricarica del combustibile di un reattore tramite la modellizzazione numerica dei fenomeni che avvengono al suo interno. Dunque vorr` o studiare quale sar` a, istante per istante, il contenuto in U 235 , in P u239 e 241 P u , il quantitativo di veleni ancora presenti, all’interno di una panoramica 3D del reattore, che potremo al pi` u semplificare sfruttandone la geometria circolare ed andando dunque a studiarne solo un quarto della sezione, estrapolando i risultati al resto del sistema. Per fare questo useremo dunque dei codici di calcolo in diffusione 3D multigruppo. Ne esistono fondamentalmente due macro-famiglie: Differenze finite fini ove viene attuata una discretizzazione spinta del dominio spaziale. In questo caso il passo tra due nodi adiacenti è dell’ordine di grandezza del libero cammino medio dei neutroni. Tale sistema è molto oneroso dal punto di vista computazionale, ma effettuando importanti approssimazioni sul flusso neutronico (in particolare considerandolo lineare tra due nodi adiacenti) tale complessità è ridotta a livelli accettabili. Coarse mesh s ove si applica una logica opposta: la discretizzazione spaziale è molto grossolana, molto pi` u grande rispetto al libero cammino medio 9 dei neutroni . Per mantenere un livello di precisione accettabile si sceglie tuttavia di aumentare il contenuto della definizione locale del flusso neutronico, che sarà generalmente una funzione polinomiale di grado 2 o 3. Il problema è che questo porta a problemi di tipo non lineare, per i quali non esiste ad oggi una dimostrazione di esistenza ed unicit` a della soluzione. Nonostante questo tali metodi, introdotti nel 1972, sono oggi molto usati nell’ambito della programmazione delle ricariche del combustibile. Parleremo sempre di condizioni di pseudo-criticità, il che ci porterà a risolvere problemi agli autovalori generati dall’omogeneità del problema di partenza. L’ipotesi di pseudo-criticità ci permetterà quindi di dimenticarci di fatto della dipendenza dal tempo. 9

si ricorda che il libero cammino medio di un neutrone è funzione della sua energia, teoricamente in modo continuo, in pratica in modo discreto nel momento in cui il dominio energetico è suddiviso in un numero finito di gruppi

59

3.2

Il consumo del combustibile nucleare

Veniamo ad alcune considerazioni sul combustibile nucleare. Sappiamo bene che, purtroppo, nel momento in cui è necessario il suo ricambio, esso è ben lontano dall’essere esaurito. La sostituzione degli elementi di combustibile infatti, a differenza della maggior parte degli altri sistemi energetici, non è legata all’esaurimento del combustibile ma all’incapacità di garantire la criticit` a del reattore. Nei momenti immediatamente successivi alla ricarica il reattore è fortemente sovracritico, e la reattività, inizialmente controllata dai sistemi di controllo, è rilasciata con il consumarsi del combustibile.10 La condizione ideale che si vorrebbe raggiungere all’interno di un reattore nucleare è quella di reattività il pi` u possibile uniforme sulle varie direzioni. Tuttavia questa condizione è difficile da raggiungere, sia nei PWR che, soprattutto, nei BWR11 . Prendiamo un generico reattore, arricchito al 2-4%. Dovrò tenere in conto, con il passare del tempo, non solo del consumo di U 235 ma anche della produzione di P u239 e P u241 , a loro volta fissili. Dunque il depauperamento della quantit` a di materiale fissile all’interno del reattore è frutto di un bilancio tra consumo di uranio e produzione di plutonio12 . Tale bilancio non è costante durante la vita del reattore, ed anche se in generale il consumo di uranio è preponderante, esiste una fase della vita del combustibile, la cui intensit` a e durata variano da reattore a reattore, in cui la formazione di plutonio prevale sul consumo di uranio. Tale fase è detta plutonium build-up13 Quando si genera plutonio all’interno di un sistema, esso vi immette reattivit` a. Questo dunque corrisponde ad una sorta di lieve “estrazione” 10

Si vuole qui ricordare inoltre che, in generale, alle operazioni di sostituzione del combustibile sono associate anche delle fasi di movimentazione delle barre non ancora esaurite in modo da ottimizzarne il consumo 11 nei BWR l’ebollizione dell’acqua all’interno del reattore provoca una forte diminuzione di densit` a di moderatore nella parte alta del reattore, alla quale corrisponde un’altrettanto marcata diminuzione di reattivit` a 12 I calcoli legati alla diffusione neutronica sono sempre basati sull’ipotesi che i prodotti di fissione non diffondano all’interno del materiale ma restino l`ı dove sono. Questo, in realt` a, non è propriamente vero, ma il loro livello di diffusione è molto ridotto, e considerato trascurabile rispetto invece a quello dei neutroni 13 Nei reattori per la produzione di energia tale plutonio viene consumato direttamente all’interno del reattore. Esistono tuttavia dei reattori, detti plutonigeni, che sono invece pensati per massimizzare la produzione di plutonio, che viene successivamente estratto ed utilizzato per fini militari. A questo proposito esistono anche reattori pensati per il processo inverso, ovvero per essere alimentati tramite il materiale nucleare di orgine militare. 7 dei 58 reattori francesi sono adatti anche a questo scopo, ma tale possibilit` a richiede ulteriori accortezze in ambito di sicurezza e controllo, in quanto la reazione di fissione del plutonio è molto meno controllabile di quella dell’uranio (β(U 235 ) ' 0.65, β(P u239 ) ' 0.22). Questo implica che in ogni caso i combustibili a base di plutonio (o comunque contenenti parti non trascurabili di plutonio) non sono mai utilizzati da soli all’interno del reattore, ma opportunamente mischiati a combustibile tradizionale

60

delle barre di controllo, il tutto legato anche alla maggiore difficoltà del controllo delle reazioni di fissione dell’uranio. Vediamo anche altri tipi di reattori. Ad esempio sappiamo che i CANDU, per come sono pensati, sono particolarmente facili da controllare. Tuttavia, dall’altro lato, la gestione del combustibile è pessima, ovvero il combustibile è sostituito quando il suo contenuto energetico è ancora molto elevato in relazione al suo contenuto iniziale. Questo gli è tuttavia permesso dal fatto che tali reattori sono alimentati ad uranio naturale e permettono il ricambio del combustibile senza dover ricorrere allo spegnimento del reattore.

3.3

Differenze finite fini

Per quanto visto avremo dunque due livelli di discretizzazione: Spaziale : nel passaggio ad un approccio numerico perdiamo come noto la continuit` a dello spazio, il che ci porterà a calcolare le grandezze fisiche solo in alcuni punti Energetica : ho come detto un numero finito G di intervalli monoenergetici Dunque l’ordine di discretizzazione del mio problema sarà Nx Ny Nz G, il cui valore pu` o facilmente essere molto elevato.

3.3.1

Diffusione monodimensionale monogruppo

Veniamo ora alla scrittura del problema, preoccupandoci quantomeno inizialmente del caso pi` u basilare, monodimensionale e monogruppo. Avrò un problema del tipo:  2 νΣ  D ddxΦ2 − Σa Φ + k f Φ = 0 

Φ(0) = Φ(a) = 0

Ove i differenti termini assumono il seguente significato: 2

D ddxΦ2 rappresenta il depapuperamento del bilancio dato dalle fughe. Tale termine non è intrinsecamente negativo, ma lo sarà sempre (o quasi) nei bilanci che andremo a considerare Σa Φ rappresenta il depauperamento dato dagli assorbimenti. Questo termine è sempre negativo, in quanto risulta dal prodotto di due termini definiti positivi cambiato di segno νΣf Φ rappresenta il termine di fissione, ed è sempre positivo in quanto prodotto di tre quantit` a definite positive.

61

Supponiamo ora di eliminare il k da tale equazione. Tale problema sar` a ovviamente risolto dalla soluzione banale di flusso neutronico nullo. Una soluzione non banale può esistere soltanto se la composizione materiale e la geometria del sistema sono tali da renderlo critico. Da qui nasce appunto la natura del k: esso corrisponde all’autovalore che, fisicamente, non funge ad altro che a riequilibratore del fenomeno fisico. Dunque io prendo un problema che, inizialmente, non rappresenterà una situazione di criticit` a, e lo modificherò tramite un parametro in modo tale da rendere il sistema critico. Una volta fatto questo andrò a determinare il valore di tale parametro che, per come è stato inserito, mi dirà che il sistema è sottocritico se inferiore ad uno e sovracritico se superiore ad uno. Ho cos`ı risolto magistralmente il problema della non stazionarietà, forzando il problema a diventare stazionario: il bilancio neutronico diffusivo è dunque sempre nullo, anche se nella realtà fisica delle cose non sarà cos`ı, ma riesco comunque a quantificare lo scostamento dalla criticità. Chiaramente la veridicit` a di tale approssimazione è garantita solo per k molto vicino ad uno. Uno scostamento eccessivo provocherebbe il verificarsi di fenomeni fortemente non-stazionari, il che richiederebbe lo studio del fenomeno sotto un approccio differente. La geometria del sistema viene in questo caso considerata all’interno delle condizioni al contorno alle pareti, ove abbiamo imposto flusso nullo.14 Dal punto di vista fisico, in questo caso l’autovalore è attribuito al termine di fissione νΣf Φ. Questa non è l’unica soluzione possibile: potrei, ad esempio, inserire il k nelle condizioni al contorno, imponendo che il flusso si annulli in (0 + δx) e (a − δx) ove l’autovalore diventa, in questo caso, δx e viene chiamato autovalore geometrico. Lo scopo è sempre lo stesso: forzare la criticit` a del sistema, andando però questa volta a modificarne la geometria piuttosto che la composizione materiale.15 Il problema della scelta dell’autovalore geometrico è che il sistema che ne risulta è ancora pi` u difficile da studiare dal punto di vista numerico in quanto implica una modifica delle dimensioni del dominio e, di conseguenza, della sua discretizzazione. La dipendenza della matrice del problema da k non è pi` u lineare, e la risoluzione ne risente. Tale modello è utilizzato in generale per altri scopi, come ad esempio per comprendere in che modo un sistema critico reagisce ad una variazione della geometria del sistema, ad esempio in caso di incidente. 14

L’approssimazione di flusso nullo alle pareti è detta di affacciamento al vuoto, ed è realistica solo ove non vi sia un riflettore. In generale, pur non considerando la presenza di un riflettore, andremo a sostituire tale condizione con l’annullamento del flusso su di un confine estrapolato 15 inserire il k all’interno del termine di fissione corrisponde infatti a diminuire il ν o il Σf , il che non corrisponde ad altro che ad una sostituzione ideale del materiale fissile con un altro avente propriet` a differenti

62

Ulteriore posizione che l’autovalore potrebbe assumere è l’andare a modificare la densit` a del materiale, in modo da variare i valori di Σa e Σf . Il problema in questo caso risiede nel fatto che ogni modifica della densità materiale andrebbe ad agire in due direzioni opposte sul sistema, il che complica nettamente la comprensione di quali siano le reali conseguenze. Inoltre si avr` a, anche in questo caso, l’introduzione di una non linearità del problema. Anche in questo caso dunque si ha una maggiore complessità della risoluzione, nonostante il suo significato fisico sia, teoricamente, pi` u realistico. Se andassi a risolvere il problema analiticamente, otterrei π Φ(x) = C sin x a

con

k=

νΣf Σa + D( πa )2

Volendo invece procedere dal punto di vista numerico, andrò ad introdurre una discretizzazione spaziale, identificando cos`ı un certo numero di nodi ove andr` o a calcolare il flusso neutronico Avrò dunque che esisteranno degli xi tali che Φ(xi ) = Φi . In questo modo le derivate diventeranno del tipo: dΦ ( dΦ dx )xi+ 1 − ( dx )xi− 1 d2 Φ 2 2 ( 2 ) xi = dx h

Si parla in questo caso di differenze finite centrate. Ho però in questa forma la presenza di definizioni del flusso al di fuori dei nodi preimpostati. Tale apparente incongruenza viene eliminata tramite la definizione delle derivate prime: (

dΦ Φi+1 − Φi )xi+ 1 = dx h 2

e

(

dΦ Φi − Φi−1 )xi− 1 = dx h 2

da cui in definitiva la derivata seconda rispetto ad x calcolata in xi vale: (

d2 Φ Φi+1 − 2Φi + Φi−1 )xi = 2 dx h2

Considerando D come omogeneo, possiamo definire dunque un sistema di equazioni dato da:  Φ −2Φ +Φ νΣ  D i+1 h2i i−1 − Σa Φi + k f Φi = 0 

Φ0 = Φn+1 = 0

Tuttavia stiamo qui come detto imponendo una condizione poco realistica di flusso nullo al contorno. In questo punto verrà in realtà effettuato un upgrade del problema diffusivo tramite sontenuto fisico proveniente dalla teoria del trasporto. In pratica andrò a definire un contorno estrapolato, diverso dal contorno fisico, sul quale il flusso è nullo, basato sulla 63

estrapolazione lineare dei risultati ottenuti dalla teoria del trasporto. Le ipotesi sulle quali la teoria della diffusione è basata, infatti, non sono verificate alla frontiera del dominio, ove è dunque necessario ricorrere ad altri strumenti. Una volta determinata la frontiera estrapolata, l’equazione della diffusione verr` a risolta non sul dominio fisico ma sul dominio estrapolato, ove potremo finalmente in modo realistico imporre l’annullamento del flusso neutronico. Avremo che, al contorno, le condizioni diverranno del tipo: γ0 Φ(r0 ) + δ0 ( dΦ dr )r0 = 0 γN Φ(rN ) + δN ( dΦ dr )rN = 0 Ove δ e γ dipendono dalla geometria scelta. Scegliendo γ0 = 1 e δ0 = −d otteniamo: Φ(r0 ) − d(

dΦ )r = 0 dr 0

dove chiameremo d distanza estrapolata, ed in essa risiederà il costo aggiuntivo del mio upgrading.16 Ho dunque ottenuto una equazione agli autovalori, ovvero un problema del tipo D(λ)Φλ = 0 La diffusione rappresenta un problema fisico di tipo locale: ciò che accade nel nodo i è direttamente collegato a ciò che accade nei nodi (i+1) ed (i1). Questa considerazione, apparentemente banale, provoca in realtà delle conseguenze importanti sul tipo di problema numerico che dovremo andare a risolvere. Problemi di trasporto sono infatti invece non di tipo locale: tutti i punti sono direttamente collegati agli altri punti. In fisica diffusiva andremo dunque a risolvere equazioni con formulazioni: • a tre punti (1D) → matrici tridiagonali • a cinque punti (2D) → matrici pentadiagonali • a sette punti (3D) → matrici eptadiagonali17 Vediamo dunque che avremo, nel caso diffusivo, matrici sparse, che ci porteranno ad approcci di tipo iterativo per la loro risoluzione. Questa è, in generale, una buona notizia, soprattutto ove la matrice sia sparsa ma ordinata, come nel nostro caso. 16

Purtroppo nemmeno la teoria del trasporto è esatta, e di conseguenza anche il valore di d calcolato in questo modo costituir` a un’approssimazione, e il suo valore cambia a seconda di come ho applicato tale teoria e di come ne ho risolto le equazioni 17 L’ottenimento di matrici multidiagonali dipende dalla numerazione scelta. Nel caso monodimensionale la numerazione è ovvia, mentre nei casi bi e tridimensionale non è cos`ı, e sar` a opportuno scegliere una numerazione adeguata per ottenere effettivamente matrici multidiagonali

64

Dalle equazioni notiamo inoltre che il termine di assorbimento e quello di fissione sono posizionati solo sulla diagonale, mentre il termine legato alla diffusivit` a D è posizionato sia sulla diagonale che fuori. Si vede in effetti che per questo tipo di matrici la condizione di dominanza diagonale è automaticamente assicurata. Andando infatti a trascurare il termine di fissione otteniamo che la condizione di diagonale dominante è soddisfatta se: D D D | − 2 2 − Σa | > | 2 | + | 2 | h h h Si vede chiaramente che questo è sempre vero, il che ci garantisce, come conseguenza, la convergenza dei metodi iterativi utilizzati. Questo vale per problemi: • Monoenergetici • Monodimensionali • Omogenei ove per omogenei si intende il fatto che le proprietà del sistema non variano nello spazio. Nella realt` a sarebbe dunque possibile, in alcuni punti, avere Σa = 0 e dunque ottenere una matrice a semi-diagonale dominante.

3.3.2

Problemi con sorgente

Facciamo ora un excursus sui problemi non omogenei, ovvero in presenza di un termine forzante. Questo ci servirà, in seguito, a trattare i problemi multigruppo, in quanto vedremo che per ogni gruppo energetico i neutroni rallentati nei gruppi precedenti risulteranno, matematicamente, come una sorgente esterna. Avremo dunque di fronte un problema del tipo (considerando problema monodimensionale, monoenergetico e con dominio omogeneo, senza fissione ma con sorgente): −D∇2 Φ + Σa Φ = S Che diventa, in seguito alla discretizzazione:  Φ −2Φ +Φ  −D i+1 h2i i−1 + Σa Φi = Si 

Φ0 = Φn+1 = 0

Dando luogo come detto ad una matrice tridiagonale. In casi come questi, ovvero di fronte ad un problema non omogeneo e con matrice tridiagonale utilizzeremo un metodo diretto: il metodo della doppia passata. Tale soluzione verrà tuttavia applicata soltanto in logica monodimensionale. Quando passeremo a 2 o a 3 dimensioni, dovendo 65

operare con matrici pi` u ricche di elementi non nulli, si ricorrerà in genere ad algoritmi di tipo iterativo. Si rivela utile in questa fase evidenziare come deve essere effettuato l’inserimento delle condizioni al contorno. Si vede infatti che le equazioni discretizzate si risolvono in formule a tre punti, ma nei due nodi agli estremi del dominio si avranno a disposizione solo due valori: il flusso nel punto stesso e quello nel punto pi` u interno del dominio. Si andr` a dunque in questo caso ad imporre la condizione di annullamento sul contorno estrapolato. Per farlo conviene integrare l’equazione della diffusione sul volume centrato sull’ultimo nodo ai confini del dominio (si supponga in questo caso r0 ). Z Z Z d2 Φ nD 2 dx − Sdx = 0 Σa Φdx + dx V0 V0 V0 Da cui si ottiene che il primo integrale diventa: Z d2 Φ dΦ + + dΦ nD 2 dx = D0+ A+ A − D 0 0 0 dx dr r 1 dr r0 V0 2

=

Φ1 D0+ A+ 0

− Φ0 γ0 + D0+ A0 Φ0 h0 δ0

Ove si è utilizzata la condizione al contorno che recita: γ0 Φ(r0 ) + δ0 (

dΦ )r = 0 dr 0

Ne viene che l’equazione discretizzata nell’ultimo nodo risulta nella forma: [D0+ A0

3.3.3

γ0 D + A+ + + Φ1 − 0 0 + Σ+ + S0 a,0 ]Φ0 + D0 A0 δ0 V0 h0 V0 h0 V0

Non omogeneit` a dello spazio

Il generico reattore non sar` a composto interamente dallo stesso materiale. Avremo dunque che i valori di Σa , Σf e D varieranno da punto a punto, ed anche in maniera marcata: basti pensare alle differenze che tali valori presenteranno tra il materiale nucleare e, ad esempio, il moderatore. Faremo in generale la supposizione che tali parametri varieranno come delle funzioni a gradino, in modo tale che esse risultino costanti all’interno di ogni macro-intervallo. Come fare per introdurre tale complicazione? Esistono fondamentalmente due strategie differenti Una possibilit` a è scegliere di posizionare alcuni nodi di flusso in corrispondenza delle discontinuit` a materiali. Esisteranno ovviamente molti altri nodi di discretizzazione, ma la caratteristica di questo metodo 66

è appunto quella di imporre che alcuni di questi nodi siano posizionati ove abbiamo il cambiamento dei valori delle costanti nucleari. Prendiamo un nodo ri centrato su di una discontinuità, e chiamiamo Vi − il volume di materiale attorno a ri e racchiuso dalle superfici A+ i ed Ai ri +ri+1 ri +ri−1 e ri− 1 = . posizionate rispettivamente sulle ascisse ri+ 1 = 2 2 2 2 Prendo ora l’equazione della diffusione e la integro sul tale volumetto Vi , ottenendo un bilancio integrale diffusivo. Z Z Z ¯ ∇Φ)dV ¯ − ∇(D + Σa ΦdV = SdV Vi

Vi

Vi

Il primo integrale diventa, sfruttando il teorema della divergenza Z Z Z dΦ dΦ ¯ ¯ ∇(D∇Φ)dV = D dA + D dA = du du Vi A+ A− i i = Di+ A+ i (

dΦ dΦ )ri+ 1 + Di− A− )r i ( dr dr i− 21 2

Se a questo punto introduciamo la linearità del flusso tra due nodi adiacenti, il che sta alla base dell’approccio diffusivo tramite differenze finite fini, possiamo ottenere delle equazioni che possono essere espresse in forma matriciale come abbiamo fatto sino ad ora, con la sola differenza che non ci sar` a un unico valore per le costanti materiali in tutta la matrice, ma valori differenti. La prima e l’ultima equazione avranno, come ovvio, soltanto due elementi invece che tre. Questo è legato alla scelta di condizioni di annullamento del flusso al contorno: lo stesso accade in 2D, ove avremo in generale 5 elementi, tranne per i nodi di bordo ove saranno 4 e i nodi di spigolo ove infine avremo 3 elementi non nulli. Per la risoluzione del problema è a questo punto necessaria la scelta della formula di quadratura da utilizzare. Esse si dividono in due macro categorie: Chiuse , ove gli estremi del dominio sono utilizzati come nodi Aperte , ove è invece vero il contrario Dunque nel nostro caso non è infatti scontata la scelta: se nell’approccio puramente diffusivo era stato imposto l’annullamento del flusso neutronico al contorno, a seguito dell’upgrading tramite teoria del trasporto questo non è pi` u vero, e di conseguenza si potrà scegliere anche di applicare formule di quadratura aperte. L’efficacia di una discretizzazione di questo tipo dipenderà evidentemente da • regolarit` a della funzione • numero di intervalli 67

Da qui la scelta del libero cammino medio neutronico come ordine di grandezza della discretizzazione18 : per poter supporre flusso lineare tra un nodo e l’altro è necessario che la maglia sia molto fine. Nel secondo caso le discontinuità non sono posizionate sui nodi di flusso, ma nei nodi “fittizi” del tipo ri+ 1 e ri− 1 , in modo che il volume di 2 2 materiale Vi centrato in ri si ritrovi ad avere costanti materiali uniformi al suo interno. La discontinuit` a dunque in questo caso non è concentrata in mezzo al volume, ma ai suoi bordi, ove avremo dunque discontinuità della derivata. Infatti: Z Z − dΦ + dΦ ΦdV + SdV = 0 Di Ai ( )r− − Di Ai ( )r+ − Σa dr i+ 12 dr i− 12 Vi Vi Si vede come in questo caso compaia soltanto un valore delle costanti nucleari D e Σa , in quanto esse sono costanti all’interno del volumetto, ma come al contempo sia necessario il lato in cui viene calcolata la derivata, che è infatti discontinua sulle due superfici di confine. In questo modo compaiono però delle derivate in nodi che non esistono: esse andranno calcolate come (

Φri+ 1 − Φi dΦ 2 )r− ' hi dr i+ 12 2

e

(

Φi − Φri− 1 dΦ 2 )r + ' hi dr i− 12 2

Questa definizione è gi` a migliore della precedente, ma continuo ad avere il calcolo di delle grandezze (questa volta il flusso neutronico) in nodi fittizi. Per riportarmi ad una formula a tre punti canonica utilizzerò l’imposizione della saldatura della corrente neutronica. ~ J~ = D∇Φ

7−→

+ − ~ ~ Di (∇Φ i+ 1 ) = Di+1 (∇Φi+ 1 ) 2

2

e, di conseguenza, essendo

Di 6= Di+1

7−→

~ + 1 6= ∇Φ ~ −1 ∇Φ i+ i+ 2

2

In effetti un qualsiasi problema matematico deve, per essere risolto, avere un numero di condizioni al contorno imposte pari al numero di gradi di libert` a19 . Avendo imposto, per ogni intervallo, la linearità del flusso, dovrò imporre al suo contorno due condizioni, che mi serviranno a determinare le costanti A e B date da: Φ = Ar + B 18

Si noti che a questo punto sar` a necessario il calcolo di λ, che dipender` a dalla sezione d’urto di assorbimento e, di conseguenza, anche dall’energia dei neutroni 19 Da qui naturalmente segue che, andando ad aumentare la precisione della modellizzazione e dunque il grado del polinomio interpolante, dovr` o andare ad imporre una nuova condizione al contorno per ogni volume elementare

68

La prima condizione è gi` a stata imposta tramite l’equazione di bilancio del flusso neutronico, la prima che abbiamo scritto. Ora imporremo la saldatura sulle interfacce del volume. Quello che potrebbe in prima battuta apparire come un ipervincolare il problema (una condizione di continuità del flusso neutronico e due di saldatura della corrente, una per ogni superficie di confine), si rivela ben presto essere in realtà solo un’apparenza. Ogni condizione vale infatti per due volumi, e di conseguenza il numero finale delle condizioni che andremo ad applicare è esattamente quello richiesto. Se ne ottiene il seguente sistema:  Φr 1 −Φi Φi+1 −Φr 1 i+ 2 i+ 2   = Di+1  Di hi hi+1 2

2

Φr 1 −Φi−1 Φi −Φr 1  i− 2 i− 2   Di−1 = Di hi−1 hi 2

2

Tramite sostituzione è poi possibile esplicitare il valore del flusso nei nodi intermedi in funzione di quello nei nodi reali. Si tratta in effetti di un sistema a due equazioni in due incognite.

3.3.4

Diffusione bidimensionale monogruppo

Vediamo ora cosa accade andando a muoverci in un dominio a due dimensioni. Manterremo la maggiore complessità legata alla non uniformità delle costanti nucleari sul dominio, in particolare riferendoci alla prima strategia scelta per tale modellizzazione, ovvero imponendo dei nodi di discretizzazione in corrispondenza delle discontinuità delle costanti nucleari. L’equazione in forma analitica sarà del tipo (con sorgente ma senza termine di fissione): −D(

∂2Φ ∂2Φ + 2 ) + Σa Φ = S ∂x2 ∂ y

Da cui si ottiene, a seguito della discretizzazione del dominio spaziale: −D[

Φi+1,j − 2Φi,j + Φi−1,j Φi,j+1 − 2Φi,j + Φi,j−1 + ] + Σa Φi,j = Si,j 2 hx h2y

Avr` o dunque, per un generico punto posizionato all’interno del dominio, cinque nodi di interesse: il nodo stesso pi` u i quattro ad esso adiacenti nelle due direzioni geometriche. Il punto ora è trovare il modo per passare dalla matrice Φi,j al vettore Φl , scegliendo una numerazione opportuna tale per cui la matrice finale dei coefficienti risulti il pi` u semplice possibile. Cerchiamo inoltre di inserire anche in questo contesto i discorsi effettuati a proposito della non uniformità delle costanti nucleari.

69

Anche in questo caso andrò a definire dei nodi intermedi, attorno ad ri,j = {xi , yj }. Tali nodi saranno gli spigoli del quadrato rappresentante il volumetto elementare costruito attorno ad ri,j . Supponiamo che tale nodo sia sede di una discontinuità delle costanti nucleari, tale per cui il volume elementare costruito attorno ad esso abbia il massimo possibile di zone differenti, vale a dire quattro, che chiameremo I, II, III, IV secondo la stessa logica con cui sono numerati i quadrati in un diagramma cartesiano. II I ri III IV Integriamo ora l’equazione della diffusione sul volume elementare Vi,j . Prendiamo ad esempio la derivata seconda lungo x: Z Z y 1 Z x 1 j+ 2 i+ 2 ∂ ∂Φ ∂ ∂Φ D dxdy = )= dy dx( D ∂x ∂x ∂x Vi,j ∂x yj− 1 xi− 1 2 Z y 21 j+ 2 ∂Φ xi+ 12 = dy(D ' )x ∂x i− 12 yj− 1 Z y 21 j+ 2 Φi+1,j − Φi,j Φi,j − Φi−1,j ' − Di− 1 ,j = dy[Di+ 1 ,j 2 2 hxi hxi−1 y 1 j− 2

Φi+1,j − Φi,j 1 1 = [ hyj−1 Di,j−1 + hyj Di,j ] + 2 2 hxi Φi,j − Φi−1,j 1 1 − [ hyj−1 Di−1,j−1 + hyj Di−1,j ] 2 2 hxi−1 Ove è stata introdotta la non uniformità del coefficiente D dividendo il dominio di integrazione nei quattro settori del volume elementare Vi,j in modo tale che in ogni sotto-dominio il valore di D risulti uniforme20 . Il risultato è assolutamente analogo per quanto riguarda la derivata seconda lungo y. Il termine contenente la sezione d’urto di assorbimento sarà invece molto pi` u agevole nel calcolo, e darà una serie di quattro contributi, differenti per ogni sotto-dominio del volume elementare. Il primo di tali contributi (mostrato a titolo di esempio, essendo gli altri perfettamente anologhi) è 20

Si chiariscono di seguito i significati degli indici e delle zone:  Zona Indicedizona spigoloesterno     (i, j) (xi+ 1 , yj+ 1   I 2 2 II (i − 1, j) (xi− 1 , yj+ 1 2 2   III (i − 1, j − 1) (xi− 1 , yj− 1   2 2   IV (i, j − 1) (xi+ 1 , yj− 1 2

70

2

dato da:

1 h, xi h, yj Σai,j Φi,j 4 Se il volume preso in considerazione è un volume di bordo, dovrò imporvi le condizioni al contorno, come sempre a scelta tra: • Annullamento sul contorno • Annullamento sul contorno estrapolato Prendiamo un caso test: supponiamo in questa occasione di usare un algoritmo di risoluzione chiuso, andando però ad utilizzare la condizione al contorno di annullamento sul contorno estrapolato, ed utilizzando una griglia semplificata dotata di solo 9 nodi disposti in un quadrato 3x3. Dalla scrittura delle equazioni e dalla scelta di un metodo di numerazione naturale ( Φ11 = Φ1 ,Φ12 = Φ2 ,Φ21 = Φ4 ), otteniamo un sistema del tipo seguente (ove alle X corrispondono gli elementi di matrice diversi da 0)   X X 0 X 0 0 0 0 0  X X X 0 X 0 0 0 0  Φ   S  1 1    0 X X 0 0 X 0 0 0  Φ   S  2    2    Φ      3    S3       X 0 0 X X 0 X 0 0   Φ4   S4     Φ5  S5  =  0 X 0 X X X 0 X 0         S6  Φ6   0 0 X 0 X X 0 0 X         S7     Φ7        0 0 0 X 0 0 X X 0   Φ8   S8     0 0 0 0 X 0 X X 0  Φ9 S9 0 0 0 0 0 X 0 X X Ne risulta in questo modo una matrice pentadiagonale, ordinata e simmetrica a blocchi, ove ogni blocco 3x3 è a sua volta simmetrico e tridiagonale. Si tratta inoltre di una matrice centrosimmetrica. L’insieme di queste caratteristiche la rendono, in generale, di agevole risoluzione.21 Andiamo ora a risolvere il problema. Supporremo per A alcune caratteristiche particolari: • Elementi della diagonale non nulli: aii > 0

∀ i

• Matrice simmetrica 21

A tale scopo si veda il paragrafo 2.1 a pagina 16 per una trattazione approfondita dei seguenti algoritmi: • Jacobi • Gauss • SOR

71

• Elementi fuori dalla diagonale negativi: ai,j < 0 ∀ i 6= j P • Matrice a diagonale dominante: |aii | > j |ai,j | ∀ i Opereremo tramite metodi iterativi: ne consegue che, tra i dati in ingresso, dovr` o avere anche un flusso di tentativo Φ(0) . Dovrò inoltre imporre uno o pi` u criteri di convergenza per decidere quando interrompere il processo iterativo

3.3.5

Diffusione multigruppo

Vediamo ora di passare alla fase successiva. Partendo dal caso pi` u standard, con discretizzazione spaziale monodimensionale ed a costanti nucleari uniformi su tutto il dominio, siamo passati al considerare una discretizzazione a due dimensioni e alla introduzione di costanti nucleari variabili all’interno del dominio. Introduciamo ora un ulteriore elemento di complicazione. Sappiamo bene infatti che nei reattori termici i neutroni, emessi con una energia media di circa 2 MeV, subiscono una serie di urti che li portano ad arrivare alle energie termiche (pari a circa 0.025 eV), alle quali essi raggiungono una sorta di equilibrio fino a quando non terminano la loro vita a seguito di fughe, eventi di fissione o assorbimenti. Finora non ci siamo preoccupati di definire l’energia dei neutroni: essa non è infatti mai comparsa esplicitamente nelle equazioni scritte, e non è mai stato necessario fino ad ora tirarla in ballo. Questa situazione è però solo apparente. Ogni volta che scriviamo l’equazione della diffusione, infatti, introduciamo delle grandezze, dette sezioni d’urto, il cui valore non varia soltanto a seconda del materiale considerato, ma anche a seconda dell’energia del neutrone incidente. Di fronte dunque ad un intervallo energetico che va dai meV ai MeV, come dobbiamo comportarci? La soluzione risiede in una discretizzazione energetica. Essa corrisponde a dividere lo spettro energetico in un numero finito di gruppi, a ciascuno dei quali verr` a assegnato un valore medio dell’energia. Da quel momento in poi, tutti i neutroni la cui energia è compresa all’interno dell’intervallo relativo ad un determinato gruppo saranno considerati come se avessero energia esattamente pari a quella media del gruppo stesso. Il passaggio di discretizzazione energetica non è semplice, ed il numero di gruppi dipender` a dal tipo di reattore.22 22

Preso un generico problema, la cui soluzione sar` a Φ(~r, E), esso verr` a discretizzato nel dominio spaziale lungo x,y e z tramite rispettivamente N,M ed L nodi e nello spettro energetico in un numero G di gruppi. L’ordine del problema da risolvere sar` a cos`ı diventato NxMxLxG. Questo evidenzia la necessit` a di ridurre al minimo il numero di gruppi che utilizziamo per la discretizzazione dello spettro energetico, nel rispetto della precisione richiesta sui risultati ottenuti. Tale valore è in genere assegnato pari a 3 - 4 per i reattori termici, 7 - 10 per i reattori a gas e circa 20 per i reattori veloci

72

Indicheremo in generale con 1 il gruppo avente energia pi` u elevata e con G quello ad energia pi` u bassa23 . Il gruppo G sarà dunque il gruppo termico, ovvero nella condizione in cui i neutroni sono in equilibrio termodinamico con il mezzo in cui diffondono. Questo implica che, in linea teorica, essi potrebbero trovarsi ad incrementare la propria energia in seguito ad un urto. In tutte le considerazioni che seguono in realtà ignoreremo questa possibilit` a (che ci provocherebbe, come vedremo in seguito, l’impossibilità di risolvere le equazioni relative ai vari gruppi in cascata), e considereremo soltanto fenomeni di down-scattering, ovvero di urti di scattering che portano ad una diminuzione dell’energia del neutrone.24 Supponiamo ora di porci nella condizione di G = 3: avrò dunque tre gruppi, rappresentativi di tre “zone energetiche”: • Zona di energia da fissione • Zona di rallentamento • Zona di termalizzazione25 Dunque la sorgente dei neutroni è contenuta nel gruppo 1, e per arrivare alla zona di fissione dovranno rallentare fino alle energie termiche evitando le fughe e gli assorbimenti Definisco spettro di fissione χ(E) la distribuzione che mi dice con quale probabilit` a i neutroni verranno emessi da fissione con determinate energie. Di questi neutroni ν(E) indica il numero di quelli emessi istantaneamente, detti neutroni pronti. Entrambi questi valori saranno influenzati non soltanto dall’energia, ma anche (soprattutto per il numero di neutroni) dal tipo di nucleo fissile che li ha generati. Prendiamo un generico gruppo intermedio. Esso si vedrà arrivare dalle energie superiori i neutroni che hanno subito rallentamento a seguito di urti di scattering.26 Il bilancio neutronico relativo a tale gruppo dovrà tenere in conto di: 23 Questa convenzione deriva dalla logica di seguire i neutroni dalla loro emissione fino al loro assorbimento: essi sono infatti emessi alle alte energie per poi essere assorbiti, nella pi` u parte dei casi, ad energie inferiori 24 Tale approssimazione è ben supportata in tutte le zone energetiche tranne in quella termica. Qui i neutroni hanno mediamente la stessa energia della matrice materiale in cui si muovono, e i gli urti di up-scattering sono dunque tendenzialmente pi` u frequenti che nelle altre zone 25 La zona di termalizzazione è quella dove le sezioni d’urto di fissione sono maggiori: stiamo infatti parlando di un reattore termico, progettato in modo che i neutroni vengano rallentati dal moderatore al fine di raggiungere l’energia termica ove la sezione d’urto di fissione dell’U 235 è maggiore. Esiste anche la possibilt` a di avere fissioni ad energie superiori, ma l’incidenza percentuale di tali eventi è di molto inferiore 26 Si noti qui che non tutti gli urti di scattering provocano l’uscita del neutrone dal gruppo energetico. Nel caso di invarianza del gruppo energetico a seguito dell’urto si parla di fenomeni di in-scattering. Pi` u la discretizzazione energetica è fine, migliore sar` a l’approssimazione di urti a salto energetico basso

73

Perdite per leakage : Termine di pozzo. Dovrò tenere in conto di un coefficiente di diffusione Dg diverso per ogni gruppo poichè in generale la capacit` a di diffusione dei neutroni in un mezzo dipenderà dalla loro energia Perdite per assorbimento :Termine di pozzo legato alla sezione d’urto di assorbimento e che include sia la cattura che la fissione Comparsa di neutroni : Termine di sorgente legato sia alle reazioni di fissione che hanno origine nel gruppo g (si faccia riferimento alla funzione di distribuzione χ(E)) sia a tutti i neutroni proveniente da gruppi con indice di gruppo inferiore Perdite per rallentamento : Termine di pozzo che porta i neutroni dal gruppo g ai gruppi di indice superiore Messo in equazioni questo diverrà:                     

1 ∂Φg vg ∂t

leakage z }| { assorbimento z }| { ~ g ∇Φ ~ g ) −Σa,g Φg + = −(−∇D comparsa z }| { rallentamento G G z }| { X X + Σs,g0 →g Φg0 + χg (νΣf )g0 Φg0 + −Σs,g Φg +Sg g 0 =1

g 0 =1

g = 1···G

Ove abbiamo in questo caso introdotto la notazione del tipo Φ(~r, Eg , t) = Φg (~r, t). Si noti come per ora il termine di comparsa neutronica tenga in conto, quantomeno in linea di principio, dell’up-scattering, in quanto la sommatoria è estesa su tutti i gruppi energetici. Considerando il problema stazionario e privo di sorgenti esterne si ottiene  ~ g ∇Φ ~ g ) − Σa,g Φg − Σs,g Φg + PG0 Σs,g0 →g Φg0 +  0 = −(− ∇D  g =1  P  0 0 + χg G (νΣ ) Φ f g g g 0 =1     g = 1···G sistema simil-diffusivo ove tutte le equazioni sono accoppiate tra loro.27 Evidentemente tale sistema ha come soluzione la soluzione banale Φ = 0. Dovr` o dunque risolvere un problema agli autovalori, ed esisterà una soluzione 27

Ovviamente il sistema sar` a chiuso solo a seguito dell’imposizione delle condizioni al contorno. Si evidenzia tuttavia qui il fatto che, ove esse siano imposte sul contorno estrapolato, si ha una ulteriore dipendenza dal gruppo energetico in quanto la distanza di estrapolazione è funzione della sezione d’urto di assorbimento la quale è a sua volta, come detto, funzione dell’energia

74

geometrica e materiale che mi assicurerà la criticità e dunque la presenza di una soluzione diversa dal vettore nullo. Andremo dunque a risolvere il problema: AΦ = kBΦ Che diventer` a un problema del tipo D(k)Φ = 0 lineare, ove dunque la dipendenza della matrice D dal parametro k è lineare. Andiamo ora a vedere come ottenere i valori da assegnare alle varie costanti che dovremo utilizzare per ogni gruppo energetico. avremo che: Z Eg−1 χ(E)dE χg = Eg

Σs,g0 →g =

1 Φg0

Z

Z

Eg0 −1

dE Eg

Z

(νΣf )g =

Eg−1

dE 0 Σs,E 0 →E Φ(~r, E 0 )

Eg0

Eg−1

1 Φg Eg Z Eg−1

ν(E)Σf (E)Φ(~r, E)dE

Φg =

Φ(~r, E)dE

Dg =

Eg R Eg−1 ~ r, E)dE D(E)∇Φ(~ Eg R Eg−1 ~ r, E)dE ∇Φ(~ Eg Z Eg−1

Σt,g =

1 Φg

Σt (E)Φ(~r, E)dE Eg

Definisco inoltre sezione d’urto di rimozione la quantità: Σr,g = Σt,g − Σs,g→g Che di fatto corrisponde alla sezione d’urto totale meno quella relativa all’inscattering. Questo termine riunisce dentro di sè tutti i termini di pozzo all’interno del gruppo g, fatta eccezione per le fughe. Da queste definizioni risultano una serie di problemi: • Le medie che devo calcolare dipendono dalle soluzioni del problema stesso • All’interno del calcolo delle medie mantengo una dipendenza dalla posizione non discretizzata, che rischio mi vada a complicare ulteriormente le cose La soluzione a questi problemi deriva del seguente teorema: 75

Teorema 11 (di separabilit` a spazio-energetica) Φ(~r, E) = ψ(~r), ϕ(E) Questo è vero solo nelle parti pi` u interne del reattore. Nelle zone di bordo, tale approssimazione non è accettabile. Tramite l’applicazione di questo teorema posso cos`ı eliminare la dipendenza spaziale dalle medie. Infatti: Z Z 1 1 ( )g = ( )Φ(~r, E) = ( )ψϕdE = Φg ψϕ E Z 1 ( )ϕdE = ϕ E Liberatomi dunque del problema della dipendenza spaziale, non posso dire lo stesso riguardo la funzione ϕ, tuttora incognita. Per trovarla effettuer` o, tramite teoria del trasporto, delle simulazioni iterative che mi porteranno ad ottenere un valore accettabile per ϕ 28 . Definizione processo iterativo esterno Date le definizioni per le differenti variabili e costanti che andremo ad utilizzare in funzione del gruppo energetico, si può ripartire da una equazione di questo tipo:  ~ g (~r)∇Φ ~ g (~r)) − Σr,g (~r)Φg (~r) + P 0 Σs,g0 →g (~r)Φg0 (~r) + χg ψ(~r) = 0  −(−∇D  g
si utilizzano in questo caso criteri di convergenza in norma, in quanto non interessa il comportamento puntuale ma quello medio

76

• l’operatore A definisce le perdite per fughe e rimozioni  ~ 1∇ ~ − Σr,1 (~r))· ∇(D  ..  ˆ0 .   .. A= .   ..  . 0ˆ ~ − Σr,G (~r))· ~ G∇ ∇(D

        

• l’operatore R rappresenta il rallentamento. Sarà di fatto una matrice triangolare inferiore a diagonale nulla in quanto abbiamo escluso i fenomeni di in-scattering ed up-scattering   0   ..  Σs,1→2 (~r)  ˆ0 .     . .  . R=  Σs,1→3 (~r) Σs,2→3 (~r)    .. . .   . .   .. . 0 • L’operatore M è rappresenta la moltiplicazione ed è legato alle sezioni d’urto di fissione   χ1 (νΣf )1 (~r) χ1 (νΣf )2 (~r) · · · · · · χ1 (νΣf )g (~r)   ..   χ2 (νΣf )1 (~r) .     .. .. M =  . .     .. ..   . . χG (νΣf )1 (~r) χG (νΣf )G (~r) • Chiamo infine L l’operatore A + R operatore di diffusione in ambito multigruppo Il problema è dunque riscrivibile in definitiva in forma matriciale come: 1 (A + R)Φ + M Φ = 0 k Dove, ricordando che (A + R) = L e definendo infine H = −L−1 posso giungere alla scrittura “finale” (HM )Φ = kΦ Giungendo dunque nuovamente ad un problema agli autovalori classico. Ovviamente tale scrittura apparentemente semplice ha un prezzo, contenuto nel calcolo di L−1 . 77

Si prosegue qui con il ragionamento fatto a proposito del metodo delle potenze. Potrei investire un numero pur rilevante di operazioni per il calcolo di H, sapendo per` o cos`ı di aver eliminato uno dei due processi iterativi concatenati ed essendomi cos`ı risparmiato, in definitiva, un quantitativo considerevole di operazioni. Il problema dell’applicazione di questo metodo, già descritto nell’ambito del metodo delle potenze, è qui ancora pi` u evidente: applicando un metodo diretto si perderebbe il parallelo fisico-numerico, ove l’indice delle iterazioni esterne corrisponde alle generazioni neutroniche. Come risolvo dunque tale problema? Utilizzando il metodo delle potenze dovr` o, per quanto detto, normalizzare la soluzione all’autovalore di modulo massimo, per evitare che essa diverga. Ma poichè tale autovalore è proprio il k che stiamo cercando, potrò scrivere (avendo definito C = H M): Φ(n+1) =

C (n) Φ k (n)

Il problema si sposter` a dunque sull’ottenimento di una legge di riaggiornamento per k su base fisica. Dimentichiamoci per un po’ che k corrisponde al kef f del mio problema, e supponiamo ora che esso sia solo il parametro di normalizzazione. Supporremo di definirlo in modo tale che, per evitare che il problema esploda, nel passare da un’iterazione all’altra il totale dei neutroni prodotti da fissione rimanga costante. In quest’ottica la formula di aggiornamento sarà data da: ~ (n+1) , νΣ ~ fe ~ (n) , νΣ ~ fe dΦ dΦ = k (n+1) k (n) da cui: k (n+1) = k (n)

~ (n+1) , νΣ ~ fe dΦ ~ (n) , νΣ ~ fe dΦ

Ove ho definito con d, e il prodotto scalare sia sulla discretizzazione spaziale che su quella energetica: ~ (n+1) , νΣ ~ fe = dΦ

G X N X

Φg(n+1) (~ ri )(νΣf )g (~ ri )

g=1 i=1

Da cui risulta che: ~ (n+1) , νΣ ~ f e =< ψ (n+1) , 1 > dΦ Da cui:29 k (n+1) = k (n)

< ψ (n+1) , 1 > < ψ (n) , 1 >

29

Si noti che qui non stiamo dando alcun peso differente alle varie zone del reattore. Una possibilit` a ulteriore sarebbe quella, come accennato parlando del metodo delle potenze, di introdurre una pesatura, ad esempio secondo Galerkin

78

Che si pu` o scrivere anche come: <

ψ (n) ψ (n+1) , 1 >=< ,1 > k (n+1) k (n)

Che corrisponde a dire che la sorgente effettiva di emissione normalizzata ` e costante iterazione dopo iterazione. In questo modo k diventa il rapporto tra i neutroni generati tra due iterazioni - generazioni successiva. Questo è ancora pi` u visibile se andiamo ad eliminare k (n) ottenendo: k (n+1) = k (n)

1 ~ fe ~ fe d k(n) C Φ~(n) , νΣ dC Φ~(n) , νΣ = ~ fe ~ fe dΦ~(n) , νΣ dΦ~(n) , νΣ

Potremmo anche vederla, per dare una ulteriore dimostrazione, nella maniera seguente: < ψ (n+1) , 1 > k (n+1) = 1 ψ (n) , 1 > < k(n) Che ci dice che, di fatto, il k di una generica iterazione è dato dal rapporto tra il numero di neutroni da fissione prodotti in quella iterazione (non normalizzati) ed il numero di neutroni provenienti dall’iterazione precedente (ovvero il numero di neutroni prodotti all’iterazione precedente, normalizzati). Sar` a ora opportuno dimostrare che il k che stiamo cos`ı calcolando è veramente il kef f del mio sistema nucleare e non un semplice elemento di normalizzazione. Dunque: Φ(2) =

1 k (1)

CΦ(1)

Da cui, all’n-esima iterazione: Φ(n) =

1 k (1) k (2) · · · k (n−1)

C n−1 Φ(1)

Che, sostituendo il valore di Φ(1) vale Φ(n) =

X 1 n−1 C [c ϕ + ch ϕh ] 1 1 k (1) k (2) · · · k (n−1) h≥2

Da cui otteniamo Φ(n) =

X λh k n−1 [c ϕ + ( )n−1 ch ϕh ] 1 1 k k (1) k (2) · · · k (n−1) h≥2

Che, moltiplicando per la matrice C ed applicando nuovamente lo stesso procedimento, ottengo: CΦ(n) =

X λh kn [c ϕ + ( )n ch ϕh ] 1 1 k k (1) k (2) · · · k (n−1) h≥2 79

Dunque il rapporto che stiamo cercando e che definisce l’aggiornamento del mio ciclo di iterazioni diventa: ~ fe dC Φ~(n) , νΣ = k (n+1) = ~ fe dΦ~(n) , νΣ P λ k n c1 dϕ1 , νΣf e + h≥2 dϕh , νΣf e( kh )n = P k n−1 c1 dϕ1 , νΣf e + h≥2 dϕh , νΣf e( λkh )n−1 Tale valore tender` a inevitabilmente a k per infinite iterazioni in quanto i termini legati alle sommatorie, essendo moltiplicati per la potenza di un numero inferiore ad 1 in modulo, tenderanno a zero. Vediamo infine di dimostrare a che valore convergerà Φ(n) Ripartiamo dalla seguente scrittura: Φ(n) =

X λh k n−1 ( )n−1 ch ϕh ] [c ϕ + 1 1 (1) (2) (n−1) k k k ···k h≥2

~f Se moltiplico scalarmente da entrambe le parti per il vettore νΣ dΦ(n) , νΣf e =

X λh k n−1 [c dϕ , νΣ e + ( )n−1 ch dϕh , νΣf e] 1 1 f k k (1) k (2) · · · k (n−1) h≥2

Ricorando che dΦ(n) , νΣf e =< ψ (n) , I > e moltiplicando ambo i membri per k si ottiene: k(n) k k (n)

< ψ (n) , I >=

X λh kn [c < ψ , I > + ( )n−1 ch < ψh , I >] 1 1 k k (1) k (2) · · · k (n) h≥2

Ricordando che, date le condizioni di normalizzazione scelte, la densita di emissione effettiva da fissione normalizzata è costante da una iterazione all’altra, e ricordando che il termine di sommatoria tende a zero per tendenza del numero di iterazioni all’infinito, ottengo che: lim

n→∞

k k (n)

< ψ (n) , I >=

kn c1 < ψ1 , I > k (1) k (2) · · · k (n)

da cui kn c1 = n→∞ k (1) · · · k (n) lim

=

< ψ (n) , I > k = < ψ1 , I > k (n) < ψ (1) , I > k (1) < ψ1 , I > k

Se supponiamo di essere giunti sufficientemente vicini alla convergenza, sappiamo che il nostro Φ(n) sarà sufficientemente parallelo all’autovettore fondamentale da poter scrivere che: CΦ(n) = kΦ(n) 80

che, sostutuito nell’equazione precedente porta alla forma C (n) Φ = n→∞ k (n) k = lim (n) Φ(n) n→∞ k kn = lim (1) c1 φ1 = n→∞ k · · · k (n) < ψ (1) , I > k = ϕ1 (1) < ψ1 , I > k

lim Φ(n+1) =

n→∞

lim

Ma sempre considerando la costanza della densità di emissione effettiva e ricordando che dopo infinite iterazioni i valori n-esimi tendono al valore esatto, ottengo: lim Φ(n) = ϕ1 n→∞

autovettore fondamentale. Dunque la mia soluzione, tramite la condizione di normalizzazione imposta, non tende ad un vettore parallelo all’autovettore fondamentale, ma esattamente ad esso. Definizione del processo iterativo interno Ci siamo fino ad ora preoccupati soltanto delle iterazioni esterne, e ci siamo imbattuti in un problema di importanza rilevante: la soluzione delle equazioni per il problema agli autovalori richiede, per il calcolo della matrice di iterazione per il metodo delle potenze, l’inversione di una matrice di dimensioni raguardevoli. Si è inoltre detto come l’inversione diretta di tale matrice non porti ad ottenere i risultati sperati, perchè impedisce di mantenere il parallello fisico-matematico tra iterazioni e generazioni neutroniche. La soluzione è dunque quella di risolvere un problema iterativo interno che dovrebbe teoricamente portare all’ottenimento della matrice di iterazione cercata. In questo caso tuttavia non si tratta pi` u di un problema agli autovalori. Il termine legato alla fissione verrà infatti considerato, per ogni equazione, come un valore noto e non dipendente dal flusso neutronico. Il problema diventer` a dunque non omogeneo. Vediamo nel dettaglio come procede il calcolo. Supponiamo di conoscere (0) (0) la soluzione di tentativo, dunque tanto ΦG quanto ψG ed anche il k (0) . Questo mi implica la possibilità di calcolare l’equazione per il gruppo 1: ( (0) 1 ~ 1 ∇Φ ~ 1 − Σr,1 Φ1 + (0) ∇D Si = 0 k (0)

S1 = χ1 ψ (0) (~r) (0)

Ove si è indicato con Si il termine ψ(~r)(0) χ0 Si tratta dunque di un prob(1) lema diffusivo non omogeneo, risolto il quale avrò un Φ1 , valore aggiornato 81

per il flusso che verr` a anche utilizzato per il calcolo del termine sorgente per il gruppo successivo. Vediamo dunque che per ogni gruppo energetico potrò scrivere un sistema del tipo AΦ = S Procedo successivamente come per quanto visto per il metodo di Jacobi: divido A in due matrici delle quali una diagonale, che porto dall’altro lato invertendola ed innescando cos`ı un processo iterativo del tipo Φ(n+1) = BΦ(n) + Q Ove la matrice B è la matrice di iterazione. Sorge ora un dubbio importante: quante iterazioni fare per i due cicli?30 • In primo luogo, è inutile aumentare eccessivamente il numero delle iterazioni interne: andrei a migliorare sempre di pi` u la precisione di un problema intrinsecamente sbagliato. Si manifesta quindi qui l’esigenza di ridurre al minimo le iterazioni interne • D’altro canto, trattandosi di due problemi concatenati, non ci basta avere la certezza della convergenza dei due metodi separatamente, ma ci serve essere sicuri che entrambi convergano simultaneamente. Di particolare importanza sarà dunque l’influenza delle iterazioni interne sulla convergenza dell’algoritmo esterno Ci chiediamo dunque se esiste una soglia per il numero delle iterazioni interne per ogni iterazione esterna che ci garantisca la convergenza dell’algoritmo esterno. Si pu` o dimostrare che la condzione meravigliosa che caratterizza il problema della diffusione multigruppo è che tale soglia esiste, ma è fissata ad una sola iterazione. Vediamo perché. Il problema che ci prefiggiamo di risolvere è sempre lo stesso: −LΦ(n+1) =

1 k (n)

M Φ(n)

Che diventa, di fatto: −LΦ(n+1) = Q(n) Ricordiamo ancora una volta che, al momento della risoluzione del problema interno, Q(n) ` e considerato noto pari al termine sorgente. Gli indici di iterazione che compaiono nelle equazioni precedenti sono relativi alle iterazioni esterne 30 Esiste una possibilit` a che potrebbe provocare il fallimento di ogni calcolo: quella per cui l’errore sul termine sorgente sia tale da provocare la non convergenza del problema interno. Questa situazione, che dipende dal condizionamento del problema, non si verifica nel caso diffusivo in quanto il problema numerico legato alla risoluzione multigruppo dell’equazione della diffusione è robusto

82

Ai fini della risoluzione iterativa il problema verrà scomposto nella maniera seguente31 Φ = JΦ + Y Q Da questa forma posso ricavare tramite alcuni passaggi una espressione per Y. Infatti so che: Φ = (I − J)−1 Y Q → −L−1 = (I −J)−1 Y → Y = −L−1 (I −J) Φ = −L−1 Q Dunque posso riscrivere il sistema, tornando all’algoritmo iterativo, nella maniera seguente: Φ[m+1] = JΦ[m] − (I − J)L−1 Q All’iterazione successiva otterrò Φ[m+2] = J 2 Φ[m] + (I + J)Y Q Da cui Φ[m+n] = J n Φ[m] + (I + J + J 2 + · · · + J n−1 )Y Q con (I + J + J 2 + · · · + J n−1 ) =

∞ X i=0

Ji −

∞ X

Ji

i=n

Ma essendo J una matrice avente proprietà particolari, ed in particolare raggio spettrale inferiore ad uno (altrimenti la convergenza non è garantita), posso supporre che la serie scritta nell’equazione precedente converga a ∞ X

J i = (I − J)−1

e

∞ X

J i = (I − J)−1 J n

i=n

i=0

Da cui Φ[m+n] = J n Φ[m] + (I − J)−1 (I − J n )Y Q Che, sostituendo l’espressione per Y, vale: Φ[m+n] = J n Φ[m] − (I − J n )L−1 Q Quindi si vede che in realt` a non faccio n iterazioni: faccio una unica iterazione che conta come n, a patto di saper calcolare la potenza n-esima della matrice J. Questo modo di operare è detto a blocchi di iterazioni. Riprendiamo il nostro algoritmo esterno, ove vale la corrispondenza: Q=

1 k (n)

31

M Φ(n)

non si è qui specificato alcun particolare algoritmo risolutivo: la matrice J potrebbe essere quella caratteristica di uno qualunque tra i metodi visti sino ad ora

83

Supponiamo dunque di fermarci dopo m iterazioni interne, e di procedere esternamente. Dunque M ]Φ(n) k (n) Per m → ∞, poichè abbiamo imposto che la matrice J converge, avremo che Φ(n+1) = [J m − (1 − J m )L−1

M (n) Φ → C = L−1 M k (n) Nel caso invece di m iterazioni avrò che l’operatore matriciale “approssimato” sar` a dato da:32 M Cm = [J m − (1 − J m )L−1 (n) ] k Si noti dunque come l’aver svolto m iterazioni interne corrisponde, nella pratica delle cose, ad aver svolto un calcolo iterativo della matrice C (il cui valore è incognito) della quale avremo dunque utilizzato un valore aggiornato all’m-esima iterazione. Tutto il problema ritorna dunque al punto: si era partiti dall’idea di non saper calcolare la matrice C e si è visto che, tramite uno ciclo iterativo interno, si è riusciti a procedere ad una sua stima, che sarà logicamente tanto pi` u precisa tanto pi` u alto è il numero delle iterazioni33 . Perchè tale espressione non provochi la perdita di fisicità del problema, deve ancora essere valido che applicando tale operatore all’autovettore fondamentale si ottiene nuovamente l’autovettore fondamentale. Questo equivale a verificare che il problema approssimato all’m-esima iterazione sia ancora un problema agli autovalori avente come soluzione gli stessi autovettori del problema esatto. In effetti, ricordando che per l’ipotesi di base del problema ϕ = −L−1 M k∗ ϕ, si ottiene Φ(n+1) = −L−1

ϕ = [J m − (1 − J m )L−1

M ]ϕ k∗

M M ϕ + J m L−1 ∗ ϕ ∗ k k ϕ = J mϕ + ϕ − J mϕ ϕ = J m ϕ − L−1

Dunque anche l’operatore matriciale “finito” ha come soluzione lo stesso autovettore dell’operatore ottenuto dopo infinite iterazioni. Dunque qualunque numero non nullo di iterazioni porta all’ottenimento di un opertore fisicamente accettabile, anche se ovviamente non correto, il che mi permette, se necessario, di ridurre ad un minimo di 1 il numero di iterazioni interne. In definitiva ci` o che si opera è la strategia seguente: Si noti come in tutte le considerazioni fatte fino ad ora compaia il termine L−1 anche nelle espressioni approssimate della risoluzione del problema. Questa è ovviamente solo una scrittura formale, in quanto se fossimo in grado di calcolre il valore di tale matrice inversa non si andrebbe ad innescare alcun processo iterativo 33 Questo procedimento è l’illustrazione nel dettaglio di quanto accennato nel capitolo sul metodo delle potenze riguardo la risoluzione di equazioni agli autovalori generalizzate 32

84

• Nella zona delle prime iterazioni si fa una sola iterazione interna per ogni iterazione esterna, in modo da non sprecare tempo di calcolo per affinare inutilmente una soluzione che so già essere sbagliata • Nella zona di congergenza asintotica dell’algoritmo posso aumentare il numero di iterazioni interne: da un lato perchè è evidente che questo porti ad una maggiore precisione della soluzione, e dall’altro perchè saranno in questa fase presenti delle approssimazioni di k tali da potermi permettere di utilizzare metodi di convergenza accelearti per le iterazioni interne Il modo di procedere in questi casi prevederà dunque di svolgere una sola iterazione interna per ogni iterazione esterna nelle fasi iniziali, per non sprecare tempo per ottenere soluzioni di base sbagliate. Giunto in una fase pi` u stabile della convergenza andrò invece ad aumentare il numero delle iterazioni interne per migliorare la precisione della soluzione. Dunque ripercorriamo le fasi risolutive: 1. Impongo una soluzione di tentativo: k (0) , ψ (0) (1)

2. Utilizzo tale soluzione di tentativo per il calcolo di Φ1 , la cui equazione risulta di fatto disgiunta dalle altre (1)

3. Utilizzo il flusso cos`ı calcolato per il calcolo di Φ2 , tenendo in conto del solo down-scattering 4. Procedendo di questo passo calcolo tutte le componenti del vettore Φ(1) 5. Aggiorno tutti i valori che mi serviranno per ripartire con le iterazioni esterne, vale a dire k e ψ

3.3.6

Accelerazione della convergenza delle iterazioni esterne

Se i sistemi sono particolarmente grandi ci troveremo in una soluzione ove il rapporto di dominanza è molto prossimo ad 1, il che provocherebbe, se non prendessimo alcun provvedimento, una convergenza troppo lenta. Oltre che prestare molta attenzione ai criteri di convergenza, che potrebbero essere fuorviati da una riduzione precoce dello pseudo-errore relativo, Vediamo dunque alcuni sistemi per accelerare la convergenza Metodo di Wielaudt s Si tratta di un metodo concettualmente semplice ma di difficile applicazione. Abbiamo visto fino ad ora che è opportuno normalizzare la soluzione ad ogni iterazione, per evitare che, dopo un numero cospicuo di iterazioni, 85

essa tenda a divergere all’infinito o a convergere al vettore nullo. Supponiamo per` o di interrompere questa pratica dopo un certo numero di iterazioni, e comunque in un momento in cui ci troviamo già nella zona asintotica del metodo iterativo. Il metodo di Wielaudt parte dalla conoscenza di un valore k* approssimazione per eccesso del valore vero di k. Tale valore può provenire sia da precedenti iterazioni che da calcoli indipendenti. Scriviamo dunque un nuovo problema: k ∗ Φ − CΦ = k ∗ Φ − kΦ

→

(k ∗ I − C)Φ = (k ∗ − k)Φ

Sar` a dunque un nuovo problema agli autovalori, ove la matrice (k ∗ I − C) è sempre invertibile in quanto abbiamo supposto k* maggiore dell’autovalore di modulo massimo e dunque la matrice risultante sarà sempre a diagonale dominante. Il problema agli autovalori potrà anche essere riscritto come: (k ∗ I − C)−1 Φ =

1 Φ k∗ − k

Ove la matrice (k ∗ I − C)−1 ha gli stessi autovettori della matrice C, ed i 1 suoi autovalori sono dati da k∗ −k h Facciamo un esempio. Supponiamo: k ∗ = 1.01 k = 1.00

k2 = 0.99

tale per cui il rapporto di dominanza del problema originario vale 0.99. La stessa quantit` a, applicando il metodo di Wielaudt, diventa pari a: 1 1.01−0.99 1 1.01−1.00

= 0.5

Applicando dunque il metodo delle potenze al problema modificato la velocit` a di convergenza è enormemente pi` u alta. Tuttavia c’è per` o un punto, pi` u o meno nascosto, in cui è contenuto il maggior costo dell’algoritmo. Si tratta in questo caso della matrice (k ∗ I − C)−1 per il cui calcolo è necessario utilizzare ancora una volta il metodo delle potenze, cosa che non sempre risulta rapido e semplice. Avrò dunque, in questo caso, tre algoritmi iterativi incantenati l’uno dentro l’altro. Il punto sta ovviamente nell’effettuare le opportune valutazioni per verificare se l’applicazione di tale metodo sia conveniente o meno. Metodo di Tchebychev Il metodo di Tchebychev non si basa sull’alterazione della matrice dei coefficienti. Supponiamo di essere arrivati ad una N-esima iterazione essendo dunque in possesso di un valore k*. 86

In base a quanto detto, avremo che: 1 CΦ(N ) k∗

Φ(N +1) =

Da cui, continuando ad applicare il metodo delle potenze interrompendo per` o il processo di aggiornamento del termine normalizzante, si ottiene: Φ(N +n) = [

1 n (N ) C] Φ k∗

Supponiamo ora di sostituire al posto della potenza n-esima del rapporto [ k1∗ C] un generico polinonmio Pn di n+1 elementi, ciascuno composto da un coefficiente e da una potenza di [ k1∗ C]. Il calcolo dei coefficienti sarà effettuato imponendo come condizione il fatto che Pn (1) = 1, come deve essere perchè sia mantenuto il significato fisico.34 Supponiamo inoltre di poter scrivere, come spesso abbiamo fatto in precedenza, la soluzione alla N-esima iterazione come composizione lineare degli autovettori della matrice C. Sfruttiamo ora una proprietà notevole dei problemi agli autovalori: partendo dal generico problema Aϕ = λϕ posso scrivere che: f (A)ϕ = f (λ)ϕ A patto che f sia una funzione sufficientemente regolare. Nel nostro caso posso allora scrivere che

35

X C C )ϕ + ch Pn ( ∗ )ϕh ' 1 ∗ k k h X k kh ' c1 Pn ( ∗ )ϕ1 + ch Pn ( ∗ )ϕh k k

Φ(N +n) = c1 Pn (

h

Ma poichè sappiamo che il rapporto

k k∗

Φ(N +n) = c1 ϕ1 +

' 1 possiamo scrivere:

X h

ch Pn (

kh )ϕh k∗

Le mie condizioni saranno date da P kh h≥2 ch Pn ( k∗ )ϕh = 0 Pn (1) = 1 Dunque una volta trovati gli n+1 coefficienti del polinomio pn potrò calcolarmi la soluzione esatta, dopo sole N iterazioni. Ma questo è teoricamente impossibile, dove sta l’inganno? 34

Il passaggio è qui molto logico: stiamo approssimando una potenza, che dunque qualunque sia l’esponente varr` a 1 se la base è 1 35 Abbiamo gi` a applicato questa regola pi` u volte, sempre con f (A) = Am

87

Evidentemente nell’utilizzo, come incognite, dei kh che ovviamente non conosciamo. Non posso dunque imporre condizioni tali per cui il contributo degli autovettori di indice superiore ad 1 sia nullo. Però posso fare in modo che esso sia il pi` u piccolo possibile. Si tratta dunque di cercare i coefficienti per il polinomio tali da minimizzare il massimo dei valori del polinomio, ovvero: mina0 ,··· ,an+1 (max(0 ≤ x ≤ σ)|Pn (x)|) Si pu` o dimostrare che la soluzione di questo problema di ricerca di minimo (problemi minmax) porta ai polinomi di Tchebychev ortogonali:   c0 (t) = 1 c1 (t) = t  cn+1 (t) = 2tcn (t) − cn−1 (t) La velocit` a di convergenza del metodo di Tchebychev è oscillante in funzione del rapporto di dominanza, ma si vede che, ad esempio, dopo 5 iterazioni di un metodo avente σ = 0.95, avremo che il metodo non accelerato vedrà il termine legato al secondo autovalore di modulo massimo ridotto a σ 5 = 0.774, mentre il polinomio di Tchebychev corrispondentemente applicato dar` a luogo ad un fattore di attenuazione pari a 0.204 . Uno dei problemi di tale metodo è la presenza di una formula ricorsiva a tre termini, ove in particolare le sottrazioni tra elementi possono generare rischi di esplosione dell’algoritmo.

3.4

Metodi Coarse Mesh

Fino ad ora è stato mostrato come risolvere problemi di fisica neutronica tramite la teoria della diffusione utilizzando metodi alle differenze finite fini. Questo corrispondeva, in buona sostanza, ad approssimare il flusso in ogni maglia nel modo pi` u semplice possibile (lineare), riguadagnando in termini di precisione diminuendo al minimo le dimensioni delle maglie stesse. Il principio che sta alla base dei metodi coarse mesh è esattamente l’opposto. Le maglie sono di dimensioni molto maggiori rispetto al libero cammino medio, è ci` o su cui si punta maggiormente è una qualità migliore della funzione interpolante Questo genera di per sè un problema. Nel momento in cui le maglie aumentano di dimensioni, esse tenderanno inevitabilmente a contenere al loro interno materiali differenti tra loro, perdendo cos`ı l’omogeneità del dominio. Per evitare di risolvere problemi in presenza di costanti nucleari non uniformi viene in generale forzata l’uniformità del sistema: i volumi elementari del metodo coarse mesh vengono “omogeneizzati”,36 ottenendo cos`ı costanti nu36

Una volta ottenuti i valori delle costanti sui volumi omogenei l’approccio utilizzato è quello di tipo b, ovvero in cui non prendo i nodi in corrispondenza delle discontinuit` a materiali ma al centro dei volumi considerati omogenei

88

cleari mediate, per il calcolo delle quali si ricorre ancora una volta alla teoria del trasporto. Si andr` a poi a risolvere il tutto tramite la teoria diffusione multigruppo. Ipotizzeremo dunque che il flusso si comporti in modo pi` u complesso, come un polinomio di 2◦ o di terzo grado, invece che linearmente.37 Da un lato la riduzione del numero di maglie mi riduce enormemente il numero di elementi delle matrici che andrò a risolvere. Dall’altro la non-linearità del flusso provoca una non-linearità anche nei problemi ad esso associati, con conseguente maggiore complicazione della risoluzione e, soprattutto, con difficolt` a molto maggiore di garantire a priori la convergenza del metodo risolutivo. Partiamo ancora una volta dal caso monocinetico, in quanto l’espansione al caso multigruppo presenta le stesse caratteristiche e problematiche viste per le differenze finite fini. Andremo dunque a scrivere, nel caso 3D, un numero di equazioni pari al numero di volumi di controllo presi in considerazione, ciascuna delle quali vedrà comparire al suo interno 7 incognite; si andr` a in seguito ad imporre le condizioni di saldatura e le condizioni al contorno. Esistono fondamentalmente due approcci distinti per i metodi coarse mesh, in base all’ordine dei polinomi scelti come funzioni interpolanti: Quabox se si sceglie di utilizzare polinomi di 2◦ grado Cubox se la scelta ricade invece su polinomi di 3◦ grado

3.4.1

Metodi Quabox

Per ragioni che saranno pi` u chiare in seguito andrò a scrivere il flusso in funzione delle tre variabili adimensionali ξ, η, ζ definite nella maniera seguente: x − xi y − yi z − zi ξ= η= ζ= h xi hyj hzl Ove ri = xi , yi , zi è la coordinata del nodo i-esimo. Di tutte le forme quadratiche che potrei scrivere tramite tali variabili quella che otterrò sarà del tipo: Φ(ξ, η, ζ) = Φijl [1 + axijl ξ + bxijl ξ 2 + ayijl η + byijl η 2 + azijl ζ + bzijl ζ 2 ] Dovr` o dunque andare ad imporre 7 condizioni per i 7 parametri utilizzati: • Una condizione di soddisfacimento del bilancio integrale sul volume di controllo 37

Si noti che una prima conseguenza di questo fatto sar` a la necessit` a di andare ad imporre un numero maggiore di condizioni al contorno

89

• Una condizione di saldatura del flusso neutronico e della corrente neutronica per ogni interfaccia Avrei dunque in totale, per ogni volume di controllo, 1 + (2 · 6) condizioni. Ma ricordando che ogni condizione all’interfaccia è imposta due volte, una su ogni lato della faccia stessa, otteniamo esattamente le 7 condizioni necessarie per la corretta definizione del problema. Cominciamo ad esplicitare il calcolo dei coefficienti scrivendo il valore del flusso all’interfaccia lungo x. Si avrà che:   Φ 1 = Φijl [1 + axijl + bxijl i+ 2 ,j,l 2 4  Φ 1 = Φijl [1 − axijl + bxijl ] i− 2 ,j,l

2

4

Sommando e sottraendo membro a membro si ottengono i valori dei due coefficienti a e b lungo x in funzione del flusso all’interfaccia. axijl =

Φi+ 1 ,j,l − Φi− 1 ,j,l 2

2

Φi,j,l

bxijl =

Φi+ 1 ,j,l − 2Φi,j,l + Φi− 1 ,j,l 2

2

Φi,j,l

Vedremo in seguito che poichè quest’ultimo non è noto e non è presente tra le incognite, si utilizzano le condizioni di saldatura sulla corrente neutronica per riscrivere le costanti in funzione del flusso ai nodi. Riprendiamo l’equazione della diffusione. La prima delle condizioni al contorno impone il bilancio integrale, del quale ci apprestiamo a calcolare le varie componenti. • Partiamo dalla componente delle fughe. Avremo che, applicando il teorema della divergenza: Z Z ~ ∇ΦdV ~ ~ ndS ∇D = (D∇Φ)~ Vijl

Sijl

= Dijl {

∂Φ ∂x

x

1 i+ 2

hyj hzl + [y] + [z]}

xi− 1

2

Ove i termini tra parentesi quadre sono gli omologhi di quanto scritto lungo x nelle direzioni y e z. Potrò a questo punto inserire la dipendenza dalla variabile ξ, per cui si avrà che: ∂Φ ∂x

∂Φ ∂ξ = ∂ξ ∂x 1 ∂Φ = = hxi ∂ξ ax 2bx ξ = Φijl [ + ]= hxi hxi ax 2bx (x − xi ) = Φijl [ + ] hxi h2xi =

90

Sostituendo questa espressione nel calcolo del bilancio integrale delle fughe si avr` a che tutti i coefficienti a scompaiono in quanto non dipendenti da x. Dunque otterremo, in definitiva, che: Z by bx bz ~ ∇ΦdV ~ ∇D = 2Dijl Vijl Φijl [ 2 + 2 + 2 ] hxi hyj hzl Vijl • Passando al termine successivo si ottiene che il bilancio integrale fornisce: Z Z νΣf νΣf ( − Σa )ΦdV = ( − Σa )ijl ΦdV = k k Vijl Vijl Z νΣf = ( − Σa )ijl Φijl Vijl (1 + bx ξ 2 + by η 2 + bz ζ 2 )dξdηdζ k Vijl Ove ancora una volta i termini dispari sono stati eliminati poichè la loro integrazione porta ad un contributo nullo. Si ottiene dunque che: Z νΣf νΣf bx + by + bz − Σa )ΦdV = ( − Σa )ijl Φijl (1 + )Vijl ( k k 12 Vijl Abbiamo cos`ı ottenuto la prima delle 7 condizioni, quella legata all’annullamento del bilancio integrale, che coinvolge solo i coefficienti dei termini quadratici. Se andiamo ora a scrivere il bilancio integrale completo otteniamo, avendo sostituito le espressioni calcolate in precedenza per i coefficienti quadratici ( ) Φi+ 1 ,j,l − 2Φi,j,l + Φi− 1 ,j,l h2xi νΣf νΣf 2 2 [Dijl + ( − Σa )ijl ] +{y}+{z}+( −Σa )ijl Φijl 24 k Φi,j,l k Si nota che, come detto, rimane la dipendenza dalle condizioni all’interfaccia, che dovr` a essere eliminata imponendo la continuità della corrente 38 neutronica Dovremo dunque imporre che: −Dijl (

∂Φ ∂Φ )x− = −Di+1,j,l ( ) + 1 ∂x i+ 2 ∂x xi+1− 12

Andando a scrivere la derivata all’interfaccia nella maniera seguente: ∂Φ ax 2bx (x − xi ) = Φijl + = ∂x x− hxi h2xi x− i+ 1 2

i+ 1 2

=

1 (3Φi+ 1 ,j,l − 4Φi,j,l + Φi− 1 ,j,l ) 2 2 hxi

38 Si noti come la scelta degli indici provoca l’automatica verifica della condizione di continuit` a del flusso neutronico:

Φ(i)+ 1 ,j,l = Φ(i+1)− 1 ,j,l 2

2

91

Che, inserito opportunamente nella condizione di saldatura, fornisce l’espressione (implicita) per il flusso neutronico all’interfaccia:     Φi+ 3 ,j,l Φi− 1 ,j,l 2 2  4(1 − 4Φi+1,j,l )   4(1 − 4Φijl )  Φ + Φ Φi+ 1 ,j,l =    ijl D h hxi+1 Di,j,l  ijl 2 ) 3(1 + hxxi i+1,j,l ) 3(1 + Dijl hx Di+1,j,l i+1

i

Ho dunque in definitiva un legame decisamente complesso del tipo Φi+ 1 ,j,l = 2 f (Φijl , Φi+1,j,l ). Si nota abbastanza velocemente che, se due mesh adiacenti hanno stessa ampiezza e stessi valori per le costanti nucleari (in questo caso in effetti è sufficiente che sia uguale il coefficiente di diffusione, ma questo di fatto a stesso materiale) si ottiene che il flusso all’interfaccia non è altro che la media aritmetica dei flussi calcolati nei due nodi adiacenti. Il problema principale della forma ottenuta è che i coefficienti stessi dipendono, a loro volta, da flussi alle interfacce. Si ha cos`ı un problema non lineare, ove dunque la matrice dei coefficienti D sarà, nel problema agli autovalori che andremo a risolvere, dipendente dagli autovalori stessi in maniera non lineare, ed anche dalla stessa soluzione che dobbiamo ottenere. Risolvere problemi di questo tipo è molto complesso. La sola possibilità è quella di procedere iterativamente, ove per eliminare i flussi alle interfacce per l’n-esima iterazione utilizzo i valori calcolati alla (n-1)-esima iterazione. Dunque la matrice d’iterazione varierà ad ogni ciclo dell’algoritmo, condizione tipica degli algoritmi per la risoluzione di problemi non lineari. Si tratta di calcoli pesanti e complessi, la cui convergenza non è garantita ma che, nonostante tutto questo, funzionano e vengono correntemente utilizzati. In generale i problemi di questo tipo, in logica multigruppo, si risolveranno con lo stesso principio applicato per le differenze finite fini: vi saranno due processi iterativi uno dentro l’altro. Il problema non è pi` u complesso in linea di principio, si operer` a con matrici di dimensioni inferiori ma con maggiore complessit` a di risoluzione. L’unico elemento effettivamente distintivo sta nel fatto che, in questo caso, non è possibile risolvere le equazioni per i vari gruppi energetici “a cascata”. Essse sono infatti accoppiate tra loro ed andranno dunque risolte in modo combinato.

3.4.2

Metodi Cubbox

Cosa accade se l’approssimazione del flusso diviene di terzo grado? Si dimostra che, scegliendo opportunamente la forma del termine di terzo grado come ξ(ξ 2 − 41 ) i coefficienti lineari e quadratici rimangono gli stessi che nel caso Quabox, e che il passo aggiuntivo consiste esclusivamente nel calcolo dei coefficienti di ordine cubico. A nostro vantaggio va il fatto che, utilizzando questa forma particolare, il termine di terz’ordine non aggiunge contributi nel bilancio integrale. Nonostante questa apparente inoffensività tali termini sono tuttavia necessari, poichè si dimostra che spesso una semplice approsimazione al secondo 92

ordine non è sufficiente. Sappiamo in ogni caso che se vogliamo aggiungere un grado di libert` a sar` a necessario imporre una condizione al contorno addizionale. In questo caso utilizzeremo delle funzioni peso: dovrà dunque essere vero che: Z νΣf ~ ~ ∇D∇Φ + ( − Σa ) ΦωdV = 0 k Vijl ove si ritorner` a al semplice bilancio neutronico per ω = 1. Utilizzando invece pesi opportuni potr` o calcolare i tre termini aggiuntivi; si dimostra che tali pesi sono del tipo ξ(ξ 2 − 41 ), che applicati nelle tre direzioni forniscono le tre condizioni al contorno aggiuntive necessarie per la chiusura del problema.39 Otterr` o cos`ı una forma del tipo: cx = (fxijl )

Φi+ 1 ,j,l − Φi− 1 ,j,l 2

2

Φijl

Ove, come era lecito attendersi, compaiono ancora una volta le condizioni alle interfacce, che andranno a complicare ulteriormente il problema.

3.5

Elementi di cinetica neutronica

Si è parlato finora di fenomeni pseudo-stazionari, ovvero dove le operazioni permettessero di distaccarsi solo in maniera marginale dalla criticità. Quando questa ipotesi non è verificata non è pi` u possibile utilizzare l’approccio visto fino ad ora e diventa necessario l’utilizzo di equazioni differenti, dinamiche, bastate sulla cinetica neutronica. Studieremo in particolare la cinetica puntiforme Cominciamo dalla definizione del kef f , il mio “metro” di criticità di un reattore. Esso pu` o essere di fatto definito in due modi equivalenti basati per` o su filosfie diverse: • secondo il life cycle • secondo il neutron balance Il risultato è naturalmente lo stesso, ma l’approccio è diverso. Nel secondo caso mi limito infatti a fare il rapporto tra i neutroni che si generano da fissione ed i neutroni perduti in un certo lasso di tempo, operando cos`ı niente pi` u di un bilancio neutronico. Attenzione per` o. Non bisogna dimenticare che da fissione non escono soltanto i neutroni pronti. Alcuni nuclei figli generati dalla fissione, essendo instabili, possono dare luogo ad ulteriore emissione di neutroni dopo un 39 Si parla ancora una volta di pesatura alla Galerkin. Tale imposizione non ha significato fisico ma solo matematico

93

certo lasso di tempo. Tali neutroni sono detti neutroni ritardati e sono di primaria importanza per il controllo delle reazioni nucleari. I neutroni ritardati vengono alla luce in una situazione in cui il flusso neutronico si è gi` a evoluto, ed è dunque diverso da quello che vi era al momento della generazione dei loro progenitori. Si dice dunque che i neutroni ritardati tendono a ricordare il flusso neutronico che li ha generati. Come ne terremo dunque in conto ai fini del calcolo del kef f ? Considereremo inizialmente due popolazioni: • La popolazione di neutroni pronti, n • La popolazione di progenitori di neutroni ritardati, ci A cosa devo il pedice i? Al fatto che esistono molto tipi diversi di progenitori che danno luogo ad emissione di neutroni, ma che lo fanno con tempi diversi. In generale essi sono divisi in 6 gruppi o famiglie, per le quali il ritardo di emissione va da un massimo di circa 80 s ad un minimo dell’ordine dei ms. Tenendo dunque in conto anche del contributo dei neutroni ritardati avremo che il kef f diventa: P n + i ci kef f = neutroni scomparsi Si vede come la scelta di utilizzare l’approccio di bilancio neutronico per il calcolo del kef f sia dettata anche dalla necessità di tenere in conto dei neutroni ritardati. Applicando una logica del tipo “ciclo di vita” non sarebbe invece facile distinguere da che generazione essi provengono. Indicheremo con l o tmg il tempo medio di generazione, in generale differente dalla vita media dei neutroni. Tali quantità sono invece uguali solo ove il reattore sia critico, ed avremo dunque una nuova definizione del kef f come kef f = ll0 . Tali quantità potranno essere dunque confuse quando si opera in regime di pseudocriticità. In realtà avremo che l0 è una costante, mentre l varia assieme al kef f . Potremmo tuttavia ora chiederci se, dopo tutto quanto detto sino ad ora, le due definizioni di kef f siano effettivamente equivalenti o meno. Supponendo per semplicità di trovarsi all’interno di un reattore termico, suddividiamo la vita di un neutrone in tappe ben definite. Avendo poi definito l0 come la vita media di un neutrone, sarà possibile dire anche percorso che l0 = , ma essendo la lunghezza del cammino percorso da un velocit` a neutrone inversamente proporzionale alla sezione d’urto di assorbimento40 potremo in definitiva dire che l0 ' Σ1a v . Si può dunque vedere che, per 40

Trascureremo in questa fase il rallentamento del neutrone, poichè i tempi caratteristici della sua termalizzazione sono molto pi` u veloci di quelli tipici dell’assorbimento: la velocit` a di un neutrone termico è molto inferiore di quella di un neutrone ancora in fase di rallentamento

94

come sono strutturati i reattori termici, in questo caso è abbastanza agevole identificare le generazioni di neutroni distinte. Possiamo dunque supporre ragionevolmente che tutti i neutroni emessi nella fase 1 saranno assorbiti nella fase 2, di modo che la definizione sul ciclo di vita risulti uguale a quella sul bilancio neutronico. La moltiplicazione del numero di neutroni che abbiamo in un tempo t è data da: t M (t) = k l Facciamo un esempio numerico. Se prendiamo l’acqua come moderatore, avremo, per t = 1s, l = 10−4 s41 Dunque, supponendo inoltre k=1.001, otterremo (utilizzando le opportune approssimazioni): M (t) = 1.00110000 ' (e0.001 )10000 ' e10 Si vede dunque che anche per un valore di kef f molto prossimo ad uno, dopo appena un secondo la cinetica del reattore ha portato alla moltiplicazione della popolazione neutronica di un fattore e10 . Questo calcolo è stato effettuato in assenza di neutroni ritardati, in presenza dei quali il valore di l diventa circa 0.08, il che riduce di molto il fattore di moltiplicazione visto prima. Definiamo ora una serie di quantità (ove per neutroni prodotti si intende la totalit` a dei neutroni prodotti ad una determinata generazione): kef f =

n prodotti n persi

kex = kef f − 1 = ρ=

kef f −1 kef f

n prodotti - n persi (reattività) n prodotti

precursori di n ritardati n prodotti

β =

ρ−β = l0 = l =

=

n prodotti - n persi n persi

n pronti - n persi n prodotti

n esistenti n persi(per unit` a di tempo)

n esistenti n prodotti(per unit` a di tempo)

In base a queste definizioni avremo che: X dn ρ−β = n(t) + λi ci (t) dt l i

41

questo dipende dal fatto che l’acqua ha una sezione d’urto di assorbimento elevata. Utilizzando la grafite avremmo avuto che l = 10−3 s

95

Ovvero la variazione del numero di neutroni per unità di tempo è pari alla produzione di neutroni pronti pi` u la somma delle produzioni dei neutroni ritardati provenienti da tutte le famiglie meno i neutroni persi.42 Il termine legato ai neutroni ritardati contiene la popolazione all’istante t moltiplicata per l’opportuna costante di decadimento, e deriva da: dci βi = −λi ci + n(t) dt l Supponiamo ora di avere un gradino di reattività nell’istante t=0, tale per cui ρ = 0 per t ≤ 0 e ρ = p per t > 0 Fintanto che tale gradino è inferiore a β, sono ancora in grado di controllare il reattore. A risposta di tale gradino ho infatti un incremento rapido per un tempo βl , dopo il quale la popolazione neutronica continua ad aumentare ma in modo molto regolare e in entità inferiore. Il salto iniziale è definito β prompt jump, ed è tale per cui n(t) = β−ρ n0 . La chiave di tale processo è che se fosse invece ρ > β, la risposta della popolazione neutronica non prevederebbe tale assestamento, ma si manterrebbe di rapido incremento. Dunque è qui che risiede il fenomeno fondamentale che permette il controllo delle reazioni a catena solo grazie alla presenza dei neutroni ritardati. Abbiamo detto in precedenza che il tempo di risposta l passa da circa −4 10 a 0.08 nel momento in cui si considerano i neutroni ritardati. Ma perchè? Questo è dovuto al fatto che il valore reale di l è dato da X βi ˆl = l + λi i

ove la sommatoria a secondo membro diventa di fatto preponderante. Ma come si arriva a tale scrittura? Ripartiamo dalla X dn ρ−β = n(t) + λi ci (t) dt l i

supponendo per ora β, ρ e l indipendenti dal tempo, anche se per la reattività si tratta di una approssimazione particolarmente forte. Supponiamo inoltre dn i di avere all’istante iniziale dc dt = 0 e dt = 0 Otterremo di conseguenza (tramite l’applicazione della trasformata di Laplace) X sβi ρ = ls − s + λi i

42

Si tenga conto che nel momento in cui vado a studiare fenomeni dipendenti dal tempo sorgono problemi detti di stiffness, rigidit` a, legati all’accoppiamento di fenomeni aventi tempi caratteristici differenti tra loro. In questo caso appare chiaro che i fenomeni legati ai neutroni pronti sono molto pi` u veloci di quelli legati ai neutroni ritardati, provocando il verificarsi della situazione sopra enunciata, con conseguente grande attenzione che deve essere posta al momento della discretizzazione del dominio temporale

96

Dunque una volta fissato ρ ho tante soluzioni quanti sono i valori di λi pi` u una, dunque di fatto 7. Dallo studio della funzione notiamo inoltre che tali soluzioni sono tutte reali e danno luogo ad un andamento monotono a tratti. Se procediamo con l’antitrasformazione otteniamo che: n(t) =

7 X

Aj e s j t

j=1

Il problema si concentra dunque sulla soluzione sj avente valore pi` u grande. Se supponiamo che tale valore sia, in modulo, trascurabile rispetto a λi otteniamo43 ρ s1 = P βi l + i λi Si ha dunque che: ˆl = l +

X βi i

=

λi

= l(1 − β) +

X

= (1 − β)l +

X

βi l +

X βi

i

i

βi (l +

i

λi

=

1 ) λi

Dunque l è la media pesata sulle popolazioni neutroniche dei tempi medi di vita dei neutroni pronti e dei neutroni ritardati. Si vede da quanto scritto che il reattore ha una cinetica sua, e dunque non è possibile variarne le condizioni operative oltre ad una certa velocità limite. Questo è molto importante soprattutto in condizioni di pericolo, ove è necessario provocare lo spegnimento repentino del reattore. Essa sarà dunque limitata dal valore della costante di decadimento dei neutroni ritardati di modulo inferiore44 43

Cosa accadrebbe se ρ > β? Ci troveremo all’interno della zona asintotica di s1 , il cui valore sar` a dunque molto elevato. Otterremo di conseguenza che nella equazione ρ = β + ls −

X λi βi s + λi i

il grande valore di s porta il contributo della sommatoria a divenire trascurabile, in modo che otteniamo per s1 : ρ−β s1 ' l 44

Questo aspetto deve essere tenuto in conto in particolare quando si sceglie di ridurre il numero delle famiglie di progenitori per la semplificazione dei calcoli, mediando le relative costanti di decadimento. In generale è necessario mantenere isolata la famiglia avente la costante di modulo inferiore in quanto essa è quella avente ruolo chiave nei transitori

97

Vediamo ad famiglie:  r,t) 1 ∂Φ(~  =  v ∂t     ∂c1 = ∂t       ∂c2 = ∂t

esempio cosa accade scegliendo di operare con due macro-

D∇2 Φ(~r, t) − Σa Φ(~r, t) + pνΣf ε(1 − β)Φ(~r, t) + p[λ1 c1 (~r, t) + λ2 c2 (~r, t)] −λ1 c1 (~r, t) + β1 pνΣf εΦ(~r, t) −λ2 c2 (~r, t) + β2 pνΣf εΦ(~r, t)

Con: p fattore di sopravvivenza alle catture di risonanza ε fattore di moltiplicazione veloce45 Se ne ottiene che: −λ1 t

c1 = c10 (~r)e

Z + β1 νΣf ε

t

0

Φ(~r, t0 )e−λ1 (t−t ) dt0

0

Si vede dunque che i neutroni ritardati risentono del flusso calcolato all’istante t-t’, come mostrato dal termine di convoluzione, che rappresenta il “ricordo“ dei neutroni ritardati. Se andiamo ad ipotizzare di poter scrivere:  Φ(~r, t) = F (~r)ϕ(t)    c (~r, t) = F (~r)C (t) i i 2 F (~ 2 F (~ ∇ r ) + B r) = 0  G   Ci Gi = νΣ fε 45

tale fattore è tenuto in conto solo per i neutroni pronti, in quanto i neutroni ritardati sono emessi con energie inferiori e la probabilit` a che essi diano luogo a fissioni veloce è trascurabile

98

Otteniamo che il sistema precedente può essere scritto nella forma46  ke f f dϕ 1   dt = lt [(1 − β)kef f − 1]ϕ(t) + lt [λ1 G1 (t) + λ2 G2 (t)]     dG1 dt = −λ1 G1 (t) + β1 ϕ(t)       dG2 = −λ G (t) + β ϕ(t) 1 2 2 dt La soluzione del sistema sar` a data da    1 ϕ(t) lt [(1 − β)kef f − 1] d    G1 (t) = β1 dt G2 (t) β2

1 lt kef f λ1

1 lt kef f λ1

−λ1 0

0 −λ2

 

 ϕ(t)  ·  G1 (t)  G2 (t)

Che in forma matriciale è esprimibile come: dx(t) ˆ = Ax(t) dt Tale sistema è risolto evidentemente da una forma di tipo esponenziale, ovvero: 3 X x(t) = µTi x(0)eαi t µi i=1

Tale problema è pi` u facilmente risolvibile applicando alle funzioni una trasformazione secondo Laplace. In questo modo, chiamando α la nuova variabile 46

sono qui di seguito esplicitati i passaggi che permettono di arrivare a questa forma, che potranno sembrare semplici ai pi` u ma che vengono mostrati per completezza: 1 ∂Φ(~r, t) v ∂t 1 ∂ϕ(~r, t) v ∂t ∂ϕ t

=

D∇2 F (~r)ϕ(t) − Σa F (~r)ϕ(t) + pνΣf ε(1 − β)F (~r)ϕ(t) + p[λ1 F (~r)Ci (t) + λ2 F (~r)Ci (t) =

=

2 −DBG F (~r)ϕ(t) − Σa F (~r)ϕ(t) + pνΣf ε(1 − β)F (~r)ϕ(t) + p[λ1 F (~r)Ci (t) + λ2 F (~r)Ci (t) =

=

2 −DBG ϕ(t) − Σa F (~r)ϕ(t) + pνΣf ε(1 − β)ϕ(t) + p[λ1 Ci (t) + λ2 Ci (t) =

=

vΣa [−(L2 B 2 + 1) + (1 − β)

= =

νΣf εp]ϕ(t) + pvνΣf ε[λ1 G1 + λ2 G2 ] = Σa (1 − β) νΣf L2 B 2 + 1 εp − 1]ϕ(t) + pvνΣf ε[λ1 G1 + λ2 G2 ] 2 2 = vΣa (L2 B 2 + 1)[ 2 2 (l B + 1) Σa L B +1 1 kef f [(1 − β)kef f − 1]ϕ(t) + [λ1 G1 (t) + λ2 G2 (t)] lt lt

Ove si è usato che: L2

=

D Σa

1 lt

=

vΣa (1 + L2 B 2 )

kef f

=

pενΣf Σa (1 + L2 B 2 )

per la cui dimostrazione si rimanda a testi e corsi di fisica del reattore

99

delle funzioni trasformate, otterrò un sistema del tipo seguente: ˆx αi x ˜ = A˜ Questo risulta essere nè pi` u nè meno che un problema agli autovalori, ove gli αi sono gli autovalori e i µi le autofunzioni. Per la seconda e la terza riga del problema (quelle relative ai progenitori) si ottiene che: β1 µ(1) − (λ1 + α2 )µ(2) = 0 β2 µ(1) − (λ2 + α3 )µ(3) = 0 Da cui si ottiene che, se il primo elemento dell’autovettore è normalizzato ad 1, vale: β2 β1 , } µT = {1, λ1 + α2 λ2 + α3 Per determinare infine gli αk devo risolvere l’equazione determinantale che assume la forma: ( P αβi lt α F (α) = 1+l + 1+l1t α 2i=1 α+λ tα i F (α) =

kef f −1 kef f

Ove ho di fatto eguagliato due espressioni per la reattività del sistema. Questo ci permette (vedi grafico su appunti) di determinare gli autovalori della matrice A. Dall’analisi dei risultati si possono trarre le seguenti conclusioni47 : • Se la reattivit` a è nulla, anche il primo degli autovalori avrà modulo nullo. Il sistema evolver` a quindi in termini molto meno incisivi (ma nonostante ci` o si muove: dunque nemmeno se il sistema fosse perfettamente critico avremmo completa stazionarietà, se non asintoticamente). • Se la reattivit` a è maggiore di 0, il primo degli autovalori è positivo, ed il sistema tende (asintoticamente) ad esplodere • se infine la reattivit` a è negativa anche il primo degli autovalori è negativo, il che porta allo spegnimento del reattore. Si vede che la velocità di tale spegnimento è proporzionale al termine pi` u resistente (ovvero α1 ), il cui valore limite (raggiungibile solo per reattività infinitamente negativa) è dato da λ1 . Ecco dunque spiegato il limite alla velocità di spegnimento del reattore. 47

Si noti come solo il primo dei 3 autovettori ha elementi sempre positivi, in quanto sia λ1 che λ2 sono maggiori dell’autovalore, il cui valore limite negativo è proprio −λ1 . Tale funzione è dunque l’unica ad avere un significato fisico

100

3.6

Elementi di Trasporto neutronico

Fino a questo momento abbiamo parlato di modellizzazione della fisica neutronica del reattore tramite l’equazione della diffusione, facendo riferimento alla teoria del trasporto come ad un sistema che, a costo di maggiore complessit` a, permette di ottenere calcoli di precisione superiore. Vediamo ora di effettuare un discorso introduttivo sul trasporto neutronico che permetta di comprenderne i principi di base. L’equazione del trasporto pi` u generica, applicabile a diversi campi della fisica, è detta equazione di Boltzmann e non è lineare. Tale non linearità nasce in generale dal considerare le interazioni tra le particelle che fanno parte della popolazione che stiamo considerando. Caratterstica del trasporto neutronico è quella di ignorare tali collisioni, di modo che l’equazione del trasporto neutronico risulter` a essere invece lineare. Esistono due tipi di formulazione per l’equazione del trasporto neutronico: Integrale : tale formulazione permette di ottenere i valori mediati delle sezioni d’urto da utilizzare in teoria della diffusione. Integro-Differenziale : tale formulazione è la pi` u usata per la progettazione dei sistemi nucleari48 Vediamo ora alcune differenze tipiche tra l’approccio diffusivo e quello trasportistico: Direzione di volo dei neutroni : Il primo elemento di maggiore precisione legato alla teoria del trasporto rispetto all’approccio diffusivo è il tenere in conto anche della direzione di volo dei neutroni. Avre~ al mo dunque a che fare con un flusso neutronico angolare Φ(~r, E.Ω) posto del consueto flusso neutronico Φ(~r, E). Questo mi permette di tenere in conto del fatto che, dopo ogni evento, un neutrone avrà una direzione di volo diversa da quella di partenza.49 Parallelismo fisico-matematico : Anche l’equazione del trasporto è generalmente risolta tramite tecniche di tipo iterativo50 . Caratteristica tipica della forma integro-differenziale è quella di mantenere un parallelismo tra iterazioni e collisioni neutroniche, ove invece in forma integrale è molto pi` u difficile trovare una correlazione di questo tipo. 48

Il codice attualmente in uso presso il CEA per il trasporto neutronico è l’Apollo II L’equazione del trasporto in forma integrale prevede una integrazione sulla direzione di volo. Il risultato dunque di tale operazione sar` a un flusso indipendente da tale variabile, il che crea un forte parallelo con la teoria della diffusione e permette di utilizzare i risultati ottenuti in teoria del trasporto per migliorarla 50 Come soluzione di tentativo viene spesso utilizzato il risultato ottenuto con modelli diffusivi 49

101

Nel caso diffusivo si era invece visto come tale parallelismo era sempre garantito Localit` a del fenomeno Si è detto pi` u volte in precedenza che l’approccio diffusivo è di tipo locale: tutti i punti del dominio sono collegati tra loro indirettamente, ma solo i nodi tra loro adiacenti presentano una dipendenza diretta. Questo si ripercuote, a livello numerico, nell’avere a che fare con matrici sparse. Nel caso della teoria del trasporto, invece, vengono tenuti in considerazione anche i fenomeni a distanza. Dunque un neutrone posto in un generico punto del reattore potrebbe, seppure con probabilit` a estremamente bassa, raggiungere un qualsiasi altro punto del reattore. Conseguenza di tale approccio è l’avere a che fare con matrici dense.

3.6.1

Formulazione integrale

Sotto le ipotesi di: • Problema stazionario • Mezzo isotropo • Fenomeni di scattering di natura isotropizzante • Sorgente isotropa Si pu` o scrivere l’equazione del trasporto neutronico in forma integrale nel modo seguente: lop (r~0 →~r,E) Z Z ∞ e Φ(~r, E) = [ Σ(r~0 , E 0 )Φ(r~0 , E 0 )f (r~0 , E 0 → E) + S(r~0 , E 0 )]dE 0 dV 0 0 2 ~ 4π|~r − r | V 0 Ove si è dunque supposto di poter scrivere (viste le ipotesi) ( ~ 0 → Ω) ~ = 1 f (r~0 , E 0 → E) f (r~0 , E 0 → E, Ω 4π ~ = 1 S(~r, E) S(~r, E, Ω) 4π L’equazione che compare nell’espressione del flusso tra parentesi quadre è una equazione di Frehdolm di seconda specie. La funzione f rappresenta la probabilità che un neutrone avente energia E’ e direzione di volo Ω0 esca da un evento di interazione con energia E e direzione di volo Ω.51 Risulterà dunque che: f (~r, E 0 → E) è dato dalla somma dei seguenti contributi: ΣS (~r, E 0 ) 0 • Scattering: fS (~r, E → E) Σ(~r, E) 51

In realt` a non si parla mai di avere esattamente energia E o E’, o direzione di volo lungo Ω o Ω0 , ma sempre di valori compresi in un intorno dE o dΩ del punto considerato

102

• Cattura: (Σc )52 Σf (~r, E 0 ) 0 0 • Fissione: ff (~r, E → E)ν(~r, E ) Σ(~r, E) Σn,2n (~r, E 0 ) 0 • Reazioni n.2n: 2 fn,2n (~r, E → E) Σ(~r, E) Dovr` a inoltre valere chiaramente che: Z ∞ f (~r, E 0 )dE 0 = 1 0

Il termine esponenziale è invece legato all’attenuazione del segnale, descrive la diminuzione di probabilità che un neutrone arrivi a distanze sempre maggiori senza subire eventi di qualunque tipo. All’interno di tale termine è contenuto il parametro “lop” lunghezza ottica del sistema, che rappresenta come neutroni aventi energie diverse vedono il sistema in modo diverso. Risulter` a di fatto che: lop (r~0 → ~r, E) = Σ(E)|~r − |r~0 | Se introduciamo l’ipotesi di mezzo omogeneo, possiamo far scomparire la dipendenza delle sezioni d’urto dalla posizione. L’equazione si presenterà dunque nella forma: Z Z Φ(~r, E) = [ V

0

∞

Σ(E)|~ r−|r~0 |

e Σ(E )Φ(r~0 , E 0 )f (E 0 → E)dE 0 + S(r~0 , E 0 )] dV 0 0 2 ~ 4π|~r − r | 0

Il problema sar` a a questo punto andare a discretizzare tale equazione a livello spaziale ed energetico. Tale operazione è tanto pi` u difficile quanto pi` u problematica è la natura del nucleo (o kernel) dell’equazione integrale che è dato, in questo caso, dal termine: ~0

eΣ(E)|~r−|r | 4π|~r − r~0 |2 Effettuando la discretizzazione della variabile energetica dovremo definire delle sezioni d’urto mediate. Si definisce: R R 0 0 r, E 0 )] g 0 dE g dE[Σ(E → E)Φ(~ R Σg0 →g = r, E 0 )dE 0 g 0 Φ(~ 52 Questo termine non compare nell’espressione precedente in quanto, evidentemente, la cattura neutronica non provoca riemissione di neutroni e di conseguenza non dar` a luogo ad alcun contributo

103

Se supponiamo di implementare nell’equazione la discretizzazione energetica, andando ad aggiungere a tutte le semplificazioni fatte sino ad ora il trascurare le reazioni n,2n, si ottiene: ( ~0 R P PG ~0 ~0 eΣ(E)|~r−|r | dV 0 ~0 Φg (~r) = V [ G g 0 =1 (νΣf )g 0 Φg 0 (r ) + Sg (r )] g 0 =1 Σg 0 →g Φg 0 (r ) + ~0 2 4π|~ r−r |

g=1:G Il passo successivo, ultima operazione prima della risoluzione numerica, sar` a quello di discretizzare nello spazio. Per farlo si utilizzano delle formule di quadratura tramite polinomi (Newton Cotes) o tramite polinomi ortogonali (Formule di Gauss). Facciamo un esempio in geometria monodimensionale:   Z a X G G X 1  Φg (x) = Σg0 →g Φg0 (x0 ) + (νΣf )g0 Φg0 (x0 ) E1 (Σg |x − x0 |dV 0 2 −a 0 0 g =1

g =1

Ove E1 è la funzione esponenziale integrale di ordine 1, un tipo di polinomio ortogonale appartenente alle funzioni speciali. Tali funzioni hanno propriet` a interessanti ma sono singolari nel dominio. Tali funzioni sono definite nel modo seguente: Z ∞ −zt e dt En (z) = tn 1 Nel nostro caso vale: E1 (z) = −γ − ln z −

∞ X

(−1)n

i=0

zn n · n!

con γ costante di Eulero Mascherani che vale circa 0.57721 Tale espressione costituisce il nucleo dell’equazione. Esso è singolare (a causa del logaritmo) in x’. Tale singolarità è detta debole in quanto è moderatamente facile da risolvere. Ci` o che si va a fare spesso è tuttavia l’utilizzare una serie di polinomi ortogonali al posto della funzione esponenziale integrale. L’approssimazione è infatti lecita e ben posta, in quanto sono entrambi polinomi ortogonali, ma si ha in pi` u una migliore conoscenza della funzione. Andremo dunque a scrivere che, utilizzando i polinomi di Tchebychev (modificati in modo tale da adattarne il dominio alla funzione di interesse) E1 (z) = −γ − ln z −

∞ X n=0

con:

z an Tn∗ ( ) 8

 ∗  T0 = 1 T ∗ = 2t − 1  1∗ ∗ (t) Tn+1 = (4t + 2)Tn∗ (t) − Tn−1 104

Tale approssimazione diventa effettivamente tale solo nel momento in cui la sommatoria sui polinomi viene portata ad un numero finito di termini. Non considerando pi` u la presenza di una sorgente avremo ottenuto un problema agli autovalori del tipo Hk ΦK = 0 Approfondiamo anche in questo caso la scelta dell’autovalore. Si era accennato, a proposito della diffusione, della possibilità di posizionare l’autovalore in altri punti dell’equazione, attribuendogli cos`ı un differente significato fisico. Questo pu` o essere fatto anche in teoria del trasporto: • Scegliendo come autovalore λ = k si va ad influenzare la molteplic` da fissione ita `. • Scegliendo come autovalore λ = γ si va ad agire sulla collisionalita • Scegliendo come autovlaore λ = δ si va infine ad agire sulla densità del materiale. Quest’ultimo caso è particolarmente interessante per modellizzare situazioni incidentali (ove questa apparente incoerenza diventa effettivamente realistica dal punto di vista fisico). Il problema è che le equazioni diventano, in questo caso, non lineari.

105

Appendice A

Aggiustamento di sezioni d’urto A.1

Introduzione alle basi di dati

La presenza di studi legati all’aggiustamento delle sezioni d’urto deriva principalmente dalla grande incertezza esistente su tali valori, ai quali nei database è sempre associato uno scostamento, il cui valore quantifica di quanto il valore reale pu` o allontanarsi da quello nominale. Ogni valore contenuto all’interno di tale intervallo è considerato in accordo con i dati sperimentali. Questo sistema ci ricorda innanzitutto dell’importanza del condizionamento di un problema: la presenza nota di incertezze nei dati in ingresso mi porta alla necessità di assicurarmi che tale imprecisione non mini la mia possibilit` a di risolvere il mio problema numerico. Negli anni ’70 è stato effettuato un enorme lavoro di aggiustamento statistico delle sezioni d’urto basato sui dati provenienti da decine di laboratori diversi, in cui è stata tenuta in conto l’affidabilità di ogni dato tramite schemi precisi, tenenti in conto ad esempio la qualità dell’attrezzatura e l’affidabilit` a dei dati di input utilizzati e, soprattutto, delle influenze che si generano tra i differenti laboratori.1 Per le sperimentazioni sulle sezioni d’urto vengono spesso utilizzati reattori particolari, di forma sferica, pensati per avere una ottima conoscenza della geometria. Lo scopo è quello di rendere il sistema il pi` u semplice possibile per eliminare ogni tipo di incertezze, in modo da portersi concentrare sulla misurazione delle sezioni d’urto. In generale il risultato di tali esperimenti (detti esperimenti integrali per via del tipo di grandezze che si 1

Si faccia l’esempio di un laboratorio (Lab 1) che effettui degli esperimenti a partire dai dati ottenuti da altri due (Lab 2 e Lab 3). Al fine di quantificare l’affidabilit` a dei risultati ottenuti dal Lab 1 sar` a dunque necessario anche tenere in conto di quale sia l’affidabilit` a dei dati dei Lab 2 e Lab 3 e la loro influenza sul processo utilizzato dal Lab 1 al fine dell’ottenimento dei suoi risultati

106

vogliono misurare) vengono confrontati con simulazioni numeriche, al fine di validare i modelli utilizzati per ottenere queste ultime. La libreria pi` u nota in fisica neutronica è la ENDF, mantenuta dai laboratori del sito di Los Alamos. Vi sono contenuti, per ogni nucleo di interesse nucleare e per ogni tipo di sezione d’uro, dati per circa 50’000 punti energetici. La libreria ENDF si divide in pi` u sezioni: ENDF-A che contiene i dati grezzi provenienti dai vari laboratori. Si tratta di dati sparpagliati e disordinati ENDF-B contiene i dati contenuti in ENDF-A opportunamente ordinati, organizzati, interpolati ed estrapolati al fine di poter essere materialmente utilizzabili. La fase di lavoro che porta dal database di tipo A a quello di tipo B è detta evaluation ENDF-C contiene infine i dati sulla correttezza e l’affidabilità di ogni dato contenuto in ENDF-B Tali dati sono purtroppo ben lontani dall’essere corretti e completi. Se già incertezze sono presenti su quei nuclidi ben noti e facili da studiare, si hanno i maggiori problemi per le specie che hanno vita molto breve. Queste infatti sono difficili da studiare, ma il loro impatto sui bilanci può essere comunque rilevante ed è dunque opportuno tenerne in conto. Attualmente i metodi numerici ed i calcolatori per attuarli hanno migliorato enormemente le loro prestazioni, e sarebbe dunque un peccato sprecare un tale potenziale a causa di dati di partenza imprecisi. Se ad esempio imponessi con successo ad un algoritmo un margine di errore dello 0.001% rischierei di scoprire che tale errore è in realtà dell’ordine del punto percentuale a causa di imprecisione sui dati di input. Il punto principale è che i parametri differenziali, quali sono ad esempio le sezioni d’urto, sono molto difficili da misurare. Molto pi` u agevole (quantomeno in relativo) è invece la misurazione dei parametri integrali, come potrebbe essere la criticit` a di un sistema. In effetti uno degli esperimenti tipici in questo ambito è quello di “costruirsi” dei benchmark in modo tale che questi risultino critici. Essendo nota la loro composizione sia dal punto di vista materiale che geometrico, i dati che li riguardano sono poi utilizzati per “testare” le librerie di sezioni d’urto, tramite le quali è possibile calcolare un valore per il kef f di tali sistemi che, nell’ipotesi che i calcoli siano corretti, dovrebbe essere identicamente pari ad 1. Purtroppo questo non accade quasi mai. Ed allora, posto che in un sistema cos`ı noto e semplificato le interferenze sulla soluzione provocate da incertezze di vario tipo sono pressochè escluse, si può pensare di utilizzare il dato misurato per ottenere un miglioramento delle librerie. In generale tuttavia ogni esperimento contiene, al suo interno, uno spazio ad N-1 dimensioni di possibili nuove librerie che portano al perfetto soddis107

facimento del calcolo ad esso associato (ove N è il numero di parametri contenuti nella libreria). Come fare per scegliere quale di tutti queste infinite possibilit` a è quella da utilizzare per aggiornare la libreria? Nel prossimo paragrafo sono introdotte (o ricordate, a seconda delle conoscenze di chi legge) alcune nozioni basilari di statistica e probabilità, utili per la comprensione di quanto segue. Nel paragrafo successivo sarà invece affrontato in modo diretto il problema dell’aggiornamento della libreria.

A.2

Elementi di statistica

Supponiamo di avere due misure della stessa quantità x effettuate da due laboratori distinti: x=x ¯ 1 ± σ1 x=x ¯ 2 ± σ2 Supponiamo che entrambi i laboratori siano affidabili. Come faccio ora a scegliere che valore usare? Dovr` o utilizzare una media di qualche tipo, sia per determinare il valore medio “reale” che per la deviazione standard. Dunque, nella forma: x ˜ = c1 x ¯1 + (1 − c1 )¯ x2 sto cercando il valore opportuno per c1 . Il mio obiettivo sarà quello di minimizzare la deviazione standard. Essendo che: σx2 = E[(x − x ˜)2 ] si avr` a σx2 = E[x − c1 x ¯1 − (1 − c1 )¯ x2 ]2 = = E[c1 x − c1 x + x − c1 x ¯1 − x ¯2 + c1 x ¯ 2 ]2 = = E[c1 (x − x ¯1 ) + (1 − c1 )(x − x ¯2 )]2 = = c21 E(x − x ¯1 )2 + (1 − c1 )2 E(x − x ¯2 )2 + 2c1 (1 − c1 )E(x − x ¯1 )E(x − x ¯2 ) = c21 σ12 + (1 − c1 )2 σ22 + 2c1 (1 − c1 )E(x − x ¯1 )E(x − x ¯2 )

Se per` o i due laboratori sono indipendenti l’ultimo termine è nullo, e la varianza della media pesata è pari alla media delle varianze pesata sui quadrati dei pesi della media.2 Andiamo dunque ad analizzare come ottenere il peso migliore in questo caso semplificato. Per farlo dovremo annullare la derivata rispetto ai pesi 2

Questo ovviamente non è detto: anzi, nel caso pi` u generale i laboratori saranno collegati tra loro, ed i risultati dell’uno influenzeranno in qualche modo quelli dell’altro e viceversa

108

dell’espressione ottenuta, in modo tale da ottenere il valore del peso c1 che minimizzi la varianza: dσx2 = 2c1 σ12 − 2(1 − c1 )σ22 = 0 dc1 Da cui ottengo: c1 =

σ12

σ22 + σ22

Se invece non elmino la covarianza ottengo: c1 =

σ22 − σ12 σ12 + σ22 − 2σ12

Dunque se i due laboratori non sono indipendenti tra loro potrei perfino avere un peso negativo. Questo corrisponderebbe al caso in cui siamo a conoscenza del fatto che uno dei due laboratori ha influenzato pesantemente l’altro in una direzione specifica, e ci costringerà a muoverci in quella opposta. Generalizzando per n laboratori ottengo la matrice di correlazione:  2  σ1 σ12 · · · σ1n  σ21 σ 2  2    ..  ..  .  . σn2

σn1

Prendiamo ora un caso distinto: supponiamo di avere due misure di due grandezze diverse: x1 = x ¯1 ± σ1

x2 = x ¯2 ± σ2

Supponiamo ora che le grandezze di interesse siano invece: y1 = x1 + x2

e

y2 = x 1 − x 2

Supponiamo ora che le due misure di partenza siano state fatte in modo tale da dare luogo a distribuzioni di probabilità di tipo gaussiano (il teorema del limite centrale ci dice che questo è vero sotto determinate condizioni). Dunque potr` o scrivere: x1 ∼ e

−

(x1 −x¯1 )2 2 2σ1

x2 ∼ e

−

(x2 −x¯2 )2 2 2σ2

La funzione di distribuzione associata ad entrambe le variabili sarà del tipo f (x1 , x2 ) ∼ e

−

(x1 −x¯1 )2 (x −x¯ )2 − 2 22 2 2σ1 2σ2

109

Questo per` o è valido soltanto se le due variabili sono tra loro indipendenti, come abbiamo supposto vero per le due grandezze misurate sperimentalmente. Se per` o vogliamo scrivere una espressione analoga per le variabili dipendenti, vale a dire y1 e y2 x1 =

y1 + y2 2

x2 =

y1 − y2 2

Sostituendo tale espressione in quanto scritto prima ottengo una forma di questo tipo: f (y1 , y2 ) = exp{ − +

1 1 1 1 [( 2 + 2 )(y1 − y¯1 )2 + ( 2 + 4 4σ1 4σ2 4σ1 1 1 1 )(y2 − y¯2 )2 + ( 2 − 2 )(y1 − y¯1 )(y2 − y¯2 )} 2 4σ2 2σ1 2σ1

Ancora una volta esiste un termine legato alla correlazione tra le due misure. Si veda come in questo caso si abbia comunque una correlazione tra le due variabili “dipendenti”, nonostante i dati sperimentali siano indipendenti. Questo è forzatamente generato dal fatto che tanto y1 quanto y2 dipendano dalle stesse variabili. Abbiamo dunque all’esponente una forma quadratica, in cui il termine generalmente noto come “doppio prodotto” è legato alla covarianza. La forma quadratica scritta in precedenza è scrivibile anche in forma matriciale (y − y¯)T M (y − y¯) . Vogliamo ora calcolare la deviazione standard di σy1 2 σy1 = E[(y1 − y¯1 )2 ] =

= E{[(x1 + x2 ) − (¯ x1 + x ¯2 )]2 } = = σ12 + σ22 + 2σ12 Ove, come al solito, se i due eventi non sono correlati, l’ultimo termine scompare. Da calcoli analoghi si pu` o arrivare a scrivere che: 2 2 σy1 = σy2

e σy1,y2 = E[(y1 − y¯1 )(y2 − y¯2 )] = σ12 − σ22 Da cui posso costruire la matrice: 2 2 σy1 σy1,y2 σ1 + σ22 σ12 − σ22 Cy = = 2 σy1,y2 σy2 σ12 − σ22 σ12 + σ22 110

Tale matrice non è una matrice qualunque. Notiamo intanto che essa è definita positiva e dunque sempre invertibile, e tale che Cy−1 = M Ho cos`ı ottenuto l’espressione della funzione di distribuzione di y. Potrò inoltre normalizzarla ad 1, ovvero renderla tale che Z f (y)dx = 1 Si dimostra che, a questo fine, devo moltiplicare per un opportuno fattore, tale per cui ottengo: p Z det(M ) 1 f (y) = ¯)T M (x − x ¯)] exp[− (x − x N 2 2 (2ϕ) Si noti che tale equazione è valida per qualunque vettore y funzione di un generico vettore x. Ci` o che varierà con la dipendenza tra i due insiemi di variabili sar` a la matrice M che genera la forma quadratica nell’esponenziale. Prendiamo ora una generica funzione r = r(p) Voglio ora trovare un modo di esprimere quanto r sia sensibile alle incertezze su p. Scriverò dunque che: r(¯ p + δp ) = r¯ + δr Da cui: δr ' Sδp La matrice S, basata sul troncamento al prim’ordine della serie di Taylor ` del problema, ed è di associata a r(p), è detta matrice di sensitivita fatto costituita dalle derivate parziali delle varie componenti del vettore r secondo le variabili contenute nel vettore p.3 Quello che mi interessa, nota la matrice S, è valutare l’entità di δr . Essa sar` a quantificata dal parametro Cr : Cr = < δr δrT >= = < (Sδp )(Sδp )T >= = < Sδp δpT S T >= = S < δp δpT > S T = = SCp S T Tale procedimento è detto sandwich rule, e definice la propagazione delle incertezze 3

Il calcolo di tale matrice, avente dimensioni I x N, sembra a prima vista un compito arduo. In realt` a posso effettuare tale calcolo in modo molto rapido e semplice

111

A.3

Aggiornamento di dati differenziali a partire da dati integrali

Supponiamo ora di avere una libreria p composta da N elementi, con N molto grande. Tale libreria è utilizzata per ricavare dei parametri integrali r in numero I molto minore di N. Supponiamo per ora che p ed r non siano correlati. Supponiamo di essere in grado di calcolare r con un modello esatto, in modo tale dunque che l’errore sulla risposta risulterà dipendere soltanto da p. Il nostro scopo è quello di aggiornare la libreria. Come accennato in precedenza, il discorso ora si ribalta: misuro r sperimentalmente e mi pongo il problema di utilizzare il confrontro tra il valore numerico e quello sperimentale per ottenere un aggiornamento per i dati di partenza. In particolare il modo utilizzato per scegliere quale sia il valore opportuno per l’aggiornamento della libreria sar` a quello di minima distanza dall’orginale: supponendo dunque che i dati di partenza fossero scorretti ma comunque sensati, l’idea è quella di fare in modo di soddisfare le condizioni sperimentali andando tuttavia a modificare la libreria di partenza il meno possibile. Questo modo di procedere è detto principio della massima verosimiglianza. Definisco dunque: p la libreria originaria p0 la libreria soddisfacente la prova sperimentale r(p) la risposta, ovvero il parametro integrale r calcolato in funzione dell’insieme dei valori contenuti nella libreria p Sn =

∂~ r ∂pn

l’elemento del vettore di sensibilità, che mi dice come la risposta r è influenzata dall’n-esimo parametro della libreria p

d = r(p) − r˜ la differenza tra il valore della risposta ottenuta numericamente utilizzando i dati della libreria non aggiornata e lo stesso valore calcolato sperimentalmente (o tramite libreria aggiornata: i due valori in questo primo caso coincidono) x = p0 − p la differenza tra la libreria “vecchia” e quella aggiornata. Queso è il valore che ci imponiamo di minimizzare. Supponendo che, comunque, la libreria di partenza non sia troppo imprecisa, e che dunque il risultato proveniente da essa non sia troppo lontano da quello reale, posso linearizzare la differenza tra i due valori, scrivendo che: ∂r 0 r(p0 ) ' r(p) + (p − p) ∂p 112

Che, ricorrendo alla notazione introdotta in precedenza, può essere anche scritta come; sx + d = 0 Il mio obbiettivo è, come detto, quello di ridurre al minimo il valore di x. Il mio problema è dunque quello di trovare un x che soddisfi l’equazione sx + d = 0 ed al tempo stesso minimizzi il valore di x2 . Per farlo il metodo principe è quello dei moltiplicatori di Lagrange, che mi permette di passare istantaneamente da un problema di minimo condizionato ad uno di minimo incondizionato. L’espressione da minimizzare diverrà dunque: R(x) = xT x + 2λ(sx + d) Ove il coefficiente 2 è posto per comodità e non influenza la correttezza del metodo. Si ottiene cos`ı che, derivando: ∂R = 2(xT + λs) = 0 ∂x Tale espressione costituisce , insieme all’equazione di partenza, N+1 equazioni per la determinazione delle N+1 incognite: le x e il moltiplicatore λ. Si ottiene dunque che: x = −λsT Che, sostituito nella sx + d = 0 fornisce: −λssT + d = 0 da cui: λ=

d ssT

ed, infine: d T s ssT Questo risultato pu` o e deve essere esteso al caso in cui non si abbia una sola risposta r, ma un vettore di I componenti integrali (in numero inferiore alle N componenti della libreria). Ciò che ne risulta è assolutamente identico al caso precedente, ove al posto del vettore s e del suo trasposto compariranno la matrice S e trasposta. Bisogna tuttavia notare che le considerazioni effettuate sino ad ora sono relative a dati privi di incertezze. Questo non è ovviamente il caso di una situazione reale, ed è opportuno dunque estendere il calcolo al caso pi` u realistico in cui i dati presi in considerazione sono affetti da incertezze. Partiamo da alcune considerazioni teoriche. Ci aspetteremo che, tra i vari paramtri contenuti nella libreria che vogliamo aggiornare, ve ne saranno di pi` u e meno precisi, ovvero alcuni valori avranno dati di incertezza pi` u ampi di altri. In quest’ottica appare chiaro che, al fine di un aggiustamento, x=−

113

cercheremo di concentrare la nostra opera sui dati pi` u incerti, che sappiamo essere da un lato pi` u “bisognosi” e dall’altro pi` u facilmente responsabili dell’errore. Per questo ci` o che andremo ad utilizzare non sarà il dato in sè, ma il suo rapporto con la sua stessa deviazione standard. In questo modo abbiamo da un lato adimensionalizzato i parametri in ingresso, e dall’altro siamo sicuri di aver tenuto in conto anche della incertezza su di essi. Possiamo dunque definire subito: x1 x2 ξ1 = ξ2 = ··· σ1 σ2 Da cui risulta che, in forma matriciale, x = Σξ r(p0 ) = r(p + x) = r(p + Σξ) Risolvendo il problema, ed accorgendosi che la matrice Σ2 non è altro che ` C, otteniamo: la matrice di sensibilita x = −CS T (SCS T )−1 d Finora si è lavorato imponendo l’esatta uguaglianza tra il risultato sperimentale e la risposta ottenuta tramite la libreria aggiornata. Questa imposizione è in verit` a poco realistica, in quanto le incertezze di calcolo mi impediranno generalmente di ottenere una eguaglianza esatta. Sarà dunque necessario supporre tali valori come differenti e cercare successivamente di minimizzare la loro differenza. r(p0 ) =6= r˜ Ovvero che il risultato numerico ottenuto tramite l’utilizzo della libreria aggiornata (e supposto dunque “esatto), definito con r’, è diverso dal risultato sperimentale, definito con r˜. Se a questo punto riprendiamo l’equazione di partenza: r(p0 ) = r(p) + S(p0 − p) e sottraiamo membro a membro il valore del risultato sperimentale: r(p0 ) − r˜ = r(p) − r˜ + S(p − p0 ) Otteniamo, definiendo y = r(p0 ) − r˜ si ottiene: y = Sx + d Andiamo ora a fare come nel caso precedente. Essendovi incertezza anche su r, dovr` o definire un nuovo vettore che tenga conto del fatto che le variabili integrali di r soggette ad incertezza maggiore dovranno di fatto contare di meno. Ecco dunque che scrivo y = Σr η, da cui, ricorrendo ancora una volta ai moltiplicatori di Lagrange: R = ξ T ξ + η T η + 2λT (SΣp ξ − Σr η) 114

Ove ho quindi imposto la condizione di minimizzare non solo il valore di (p − p0 ) ma anche quello di (r(p0 ) − r˜. Si arriva tramite i passaggi visti in precedenza a scrivere che4 : x = −Cp S T (Cr˜ + SCp S T )−1 d

y = Cr˜(Cr˜ + SCp S T )−1 d

Vediamo per` o ora di vedere la cosa da un punto di vista pi` u formale. Definite delle distribuzioni gaussiane per p ed ~r (entrambi sono infatti soggetti ad incertezza, a differenza di p’ che è la libreria ”vera“ ed r’ che è il corrispondente risposta ”esatta“), potremo scrivere la densità di probabilità a due variabili per p ed r, supponendo che entrambe siano delle gaussiane centrate sui valori ”esatti” p’ ed r’ : 1 f (p, ~r) ∼ exp{− [(p − p0 )T Cp−1 (p − p0 )] + [(˜ r − r0 )T Cr−1 (˜ r − r0 )]} 2 Il problema è legato alla determinazione di p’, libreria teoricamente esatta e praticamente incognita. Per procedere definisco x = p − p0 e y = r˜ − r0 , da cui l’esponente dell’esponenziale diventa: Q = xT Cp−1 x + y t Cr˜−1 y Tale valore, per la risoluzione del problema, dovrà essere minimizzato. Utilizzeremo ancora una volta i moltiplicatori di Lagrange, definendo cos`ı un nuovo problema Q’ Q0 = xT Cp−1 x + y t Cr˜−1 y + 2λT (Sx − y) Cercher` o il minimo per x ed y, derivando lungo entrambe le direzioni, il che mi fornisce i valori: ( ∂Q0 −1 T T ∂x0 = 2Cp x + 2S λ → x = −Cp S λ ∂Q −1 → y = −Cr˜λ ∂y = 2Cr˜ y − λ Essendo d = y − Sx ottengo che: d = Cr˜λ + SCp S T λ In questa espressione l’unico elemento che non conosco è λ. Esplicitando l’equazione per λ e sostituendo nel valore di ottimo trovato per x ottengo x = −Cp S T (Cr˜ + SCp S T )−1 d 4

Si noti la differenza tra Cr e Cr˜. Il primo rappresenta la matrice di covarianza della risposta numerica, ottenibile di conseguenza come funzione della matrice di covarianza della libreria Cp ; il secondo rappresenta invece la matrice di covarianza della misurazione sperimentale, che nulla ha a che vedere, a priori, con le librerie numeriche

115

da cui, in definitiva p0 = p − Cp S T (Cr˜ + SCp S T )−1 d Ecco dunque trovata un’espressione per la variazione che dovrò dare alle variabili p in funzione del confronto tra i dati sperimentali e i calcoli numerici sui valori integrali. Questo metodo richiede l’inversione di una matrice, ma si vede come le dimensioni di tal matrice siano di ordine I, quindi molto ridotto rispetto alle dimensioni N del vettore p. La base di tale teoria è la linearità delle perturbazioni, o meglio la possibilit` a di linearizzarle senza troppi sconvolgimenti.

A.4

Teoria delle perturbazioni

Ravetto parte da un problema agli autovalori, tanto per fare le cose semplici. Prendiamo ci` o che dicono le dispense, e partiamo da un banale problema del tipo: Hψ = S ove supporremo, per dare un senso pratico ai conti, che H sia l’operatore dell’equazione del trasporto neutronico, ψ sia il flusso neutronico stesso e S sia una generica sorgente. Definiamo inoltre la variabile sperimentale r (response) che dipender` a dal flusso neutronico e che è di fatto utilizzata per misurarlo. Sar` a dunque: r = Dψ Supponiamo ora che vi siano delle variazioni nell’operatore H che porteranno inevitabilmente ad una diversità nella soluzione. Dunque avremo che: (H0 + δH)(ψ0 + δψ) = S che diventa, essendo S = H0 ψ0 . H0 δψ = −δHψ0 Ove il pedice 0 indica l’operatore non perturbato. La perturbazione sul flusso generer` a anche una perturbazione nella risposta, data da: δr = r − r0 = = (Dψ) − D(ψ0 ) = = (Dδψ) Il punto risiede ora nella valutazione del termine δψ. Per farlo ricorriamo al problema aggiunto, definito come tale per cui (g, Af ) = (f, AT g) 116

ove g ed f sono due generiche funzioni. Dunque se Hψ = 0 è l’equazione per il trasporto in forma omogenea, allora H T ψ ∗ = 0 è l’equazione aggiunta, ove la sua soluzione sarà il flusso aggiunto. A questo punto posso supporre di risolvere anche il problema aggiunto, per il quale devo ancora imporre una sorgente, essendo infatti stato definito fino ad ora in termini omogenei. Supporrò dunque che la mia ψ ∗ sia soluzione del problema aggiunto H T ψ ∗ = D. Potrò cos`ı andarlo a sostituire ottenendo: δr = δψD = δψH T ψ ∗ che, sfruttando la definizione di operatore aggiunto, diventa: δr = ψ ∗ Hδψ Ricordando quanto scritto in precedenza posso infine scrivere: δr = −ψ ∗ δHψ0 ' −ψ0∗ δHψ0 Da cui ineffetti posso trarre il valore cercato di δr senza passare per la risoluzione di problemi perturbati ma soltanto utilizzando le soluzioni non perturbate del problema diretto e del problema aggiunto. Lo stesso pu` o essere applicato ad un generico problema agli autovalori A x = k B x. In questo caso andrò a risolvere il problema: (A + δA )(x + δx ) = (k + δk )(B + δB )(x + δx ) Lo scopo della mia ricerca sarebbe quello di ottenere il valore per δk senza passare per (x + δx ), ovvero senza essere costretto a risolvere il problema perturbato. Se facciamo il prodotto di tutti i termini, tralasciando gli infinitesimi di ordine superiore al primo, otterremo: Ax + Aδx + δA x = kBx + kBδx + kδB x + δk Bx Ricordando che il problema originario recita A x = k B x, posso scrivere: Aδx + δA x = kBδx + kδB x + δk Bx Definisco ora ξ operatore aggiunto, definito in modo tale che, date f e g due generiche funzioni, sarà: < g, ξf >=< ξ T g, f > Vado dunque a risolvere il problema aggiunto: AT x∗ = kB T x∗ 117

ove k rimane tale in quanto risulta proprietà degli operatori aggiunti quella di mantenere gli autovalori degli operatori originari. A questo punto moltiplichiamo tutta l’equazione scalarmente per la soluzione del problema agginto x∗ , ottenendo cos`ı: hx∗ , Aδx i + hx∗ , δA xi = khx∗ , Bδx i + khx∗ , δB xi + δk hx∗ .Bxi Sfruttando le propriet` a degli operatori aggiunti ottengo che il primo termine a destra dell’uguale ed il primo a sinistra diventano: hx∗ , Aδx i = hAT x∗ , δx i

e khx∗ , Bδx i = khB T x∗ , δx i

Che sommati tra loro danno contributo nullo in quanto rappresentano l’eqauzione dell’operatore aggiunto: hAT x∗ − kx∗ B T x∗ , δx i = 0 Ci siamo cos`ı liberati della perturbazione della soluzione δx , e possiamo scrivere che: hx∗ , δA xi − khψ ∗ , δB xi δk = hx∗ , Bxi ove ancora una volta è necessario utilizzare solo le perturbazioni note e le soluzioni del problema di partenza e del suo aggiunto.

118

Metodi Numerici v2.2

Recommend Documents