UNIVERSITA DEGLI STUDI DI TRIESTE

Facolta di Scienze Matematiche, Fisiche e Naturali

Corso di Laurea Triennale in Fisica

Tesi di Laurea

Simulazione numerica dello stress termomeccanico in unellissometro

Relatore: Laureando:

Prof. Edoardo Milotti Alessandro Candolini

Anno Accademico 2007-2008

Ai miei genitori.

Indice

Prefazione 1

1 Modellizzazione dello stress termomeccanico in uno slab di

vetro ottico 3

1.1 Birifrangenza indotta da stress meccanico in vetri ottici . . . 31.2 Equazione di diffusione per la conduzione termica . . . . . . . 51.3 Modellizzazione dell’assorbimento della radiazione laser . . . 71.4 Integratori numerici per equazioni differenziali stocastiche . . 11

2 Cenni di analisi agli elementi finiti 13

2.1 Metodo degli elementi finiti per la conduzione termica stazio-naria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2 Problemi dipendenti dal tempo al prim’ordine. . . . . . . . . 18

3 Calcolo termico e risultati preliminari 21

3.1 Implementazione . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Risultati preliminari . . . . . . . . . . . . . . . . . . . . . . . 24

Bibliografia 29

v

Prefazione

In ellissometri di precisione come quello dell’esperimento PVLAS e impor-tante considerare tutti gli effetti che possono portare a errori sistematicie fluttuazioni casuali nell’ellitticita misurata. In questa tesi mi occupo inparticolare di un effetto che potrebbe eventualmente spiegare sia alcuni fe-nomeni di deriva osservati nelle misure di ellitticita in PVLAS, sia il livellorelativamente elevato di rumore che finora ha impedito di effettuare ancoraun test di precisione della QED: l’assorbimento della luce laser negli elementiottici dell’ellissometro, e in particolare negli specchi, produce dei gradientitermici e conseguentemente degli stress meccanici che in presenza di rumoredi puntamento del laser possono produrre effetti spuri.

In questa tesi, si espone una modellizzazione del campo di tempera-tura ottenuta risolvendo l’equazione differenziale alle derivate parziali chedescrive il riscaldamento di uno specchio ideale tramite il metodo di solu-zione ad elementi finiti. Nel primo capitolo e descritto piu nel dettaglio ilproblema fisico, nel secondo si espongono i fondamenti teorici e computazio-nali del metodo a elementi finiti e nel terzo e ultimo capitolo sono descrittil’implementazione in C++ e alcuni risultati preliminari.

Ringraziamenti

Un sincero ringraziamento a Carlo, Davide e Roberta, per i molti e prezio-si suggerimenti e il continuo incoraggiamento durante la stesura di questolavoro, e soprattutto perche sono degli amici stupendi!

1

Capitolo 1

Modellizzazione dello stress

termomeccanico in uno slab

di vetro ottico

1.1 Birifrangenza indotta da stress meccanico in

vetri ottici

I vetri ottici possono esibire sia birifrangenza statica che dinamica, que-st’ultima indotta dalla presenza di stress meccanici, ad esempio dovuti aisupporti ottici. Poiche i vetri presentano un debole assorbimento della luce,fasci laser di alta intensita che attraversano elementi ottici possono portare ariscaldamenti locali anche di diversi gradi rispetto alla temperatura ambien-te. La dilatazione locale del vetro produce allora degli stress meccanici nonuniformi e conseguentemente birifrangenza uniassiale con assi che cambianoda punto a punto nell’elemento ottico.

Un esempio di apparato sperimentale in cui e bene tenere presente questieffetti e quello dell’esperimento PVLAS, che si propone di indagare le pro-prieta di polarizzazione del vuoto previste dall’elettrodinamica quantistica.Lo schema dell’apparato e una breve descrizione dell’esperimento sono ri-portati in Figura 1.1. Una descrizione completa dell’apparato si trova in[26]. Lo scopo principale di questa tesi e andare avanti nella comprensionedegli effetti associati alla propagazione della luce negli specchi interferenzialiche costituiscono la cavita di Fabry-Perot. Nel seguito, per la simulazionesi fara riferimento ai valori dei parametri fisici dell’apparato di PVLAS perle caratteristiche dell’ellissometro e del laser.

Nel caso di uno stress uniassiale, la birifrangenza indotta puo essereformalizzata assegnando due diversi valori dell’indice di rifrazione n‖ e n⊥

rispettivamente per direzioni di polarizzazioni della luce incidente parallelae perpendicolare all’asse di stress. La differenza n‖ − n⊥ e legata in ogni

3

4 Cap. 1: Modellizzazione dello stress termomeccanico in uno slab di vetro ottico

Figura 1.1. Schema dell’apparato sperimentale di PVLAS. Il fascio laser entrada sinistra e viene polarizzato linearmente da un polarizzatore di precisione, quin-di il fascio entra nella cavita di Fabry-Perot definita dai due specchi. All’internodella cavita la luce attraversa anche un intenso campo magnetico bipolare. Si puopensare che la cavita agisca come moltiplicatore del percorso dei fotoni che rim-balzano tra i due specchi e in questo modo viene amplificato il minutissimo effettodi interazione tra luce e campo magnetico previsto dall’elettrodinamica quantisti-ca. In presenza di questo campo magnetico il vuoto si comporta come se fosse unmezzo birifrangente uniassiale e quindi dopo l’analizzatore, che e incrociato con ilpolarizzatore, si osserva ancora un segnale luminoso dovuto al cambiamento dellostato di polarizzazione causato dalla birifrangenza magnetica. Questa birifrangen-za viene modulata a frequenza ΩMod ruotando il campo magnetico. Poiche questarotazione avvviene a una frequenza bassa a causa della grande massa del magneteche viene fatto ruotare, il cammino ottico include anche un modulatore di ellitticitache sposta piu in alto la frequenza portando quindi il segnale al di fuori del piccodi rumore 1/f a basse frequenze.

punto all’entita σ dello stress in quel punto tramite la relazione empirica

n‖ − n⊥ = Kσ , (1.1)

dove K e la “costante fotoelastica” del materiale. Questa costante in generedipende sia dalla temperatura che dalla lunghezza d’onda: nel nostro caso ladipendenza dalla lunghezza d’onda e irrilevante visto che in PVLAS si uti-lizza radiazione monocromatica. La dipendenza dalla temperatura potrebbeessere importante, pero per molti vetri le misure di K a diverse tempera-ture non sono accessibili e d’altra parte misure effettuate in quarzo amorfoper fibre ottiche danno variazioni relative dell’ordine di 10−4 K−1 e quindise assumiamo che le variazioni relative per i vetri ottici di interesse sianodello stesso ordine di grandezza si vede che in prima approssimazione possia-mo trascurarle nell’intervallo di temperature effettivamente realizzato neglielementi ottici in PVLAS. Valori tipici di K per diversi vetri ottici sono ri-portati nella Figura 1.2. Dalla figura si vede che i valori di K sono distribuitiperlopiu nell’intervallo tra 1×10−6 e 4×10−6 mm2/N e il coefficiente stessoe quasi indipendente dal valore dell’indice di rifrazione del vetro.

La Eq. (1.1) e valida qualunque sia l’origine dello stress meccanico σ. Inparticolare, lo stress meccanico puo essere indotto da forti gradienti termiciall’interno del vetro. Lo stress meccanico indotto termicamente da un gra-

1.2 Equazione di diffusione per la conduzione termica 5

Figura 1.2. Valori della costante fotoelastica per diversi vetri ottici. I valori sonoriferiti alla temperatura di 21 C e alla lunghezza d’onda di 589.3 nm. Fonte: [23].

diente di temperatura e diretto nella direzione del gradiente e la sua intensitae approssimata dalla legge empirica

σ =αE

1 − ν∆T , (1.2)

dove α e il coefficiente di espansione termica lineare E e il modulo di Young,ν e il rapporto di Poisson e ∆T = T − T0 e la differenza tra la temperaturaT nel punto e una temperatura di riferimento T0 alla quale non si manifestaalcuno stress termico indotto. Per maggiori dettagli tecnici, si puo fare riferi-mento ai data-sheet della SCHOTT all’indirizzo http://www.schott.com/

advanced_optics/english/tools_downloads/download/index.html.

1.2 Equazione di diffusione per la conduzione ter-

mica

Il formalismo matematico per descrivere la conduzione di calore nei materialie quello della diffusione.

All’interno del materiale che si vuole studiare si definisce una densita dicorrente di calore j(t, x) in modo che j(t, x) · ndS sia la quantita netta dicalore trasferita per unita di tempo (i.e., la corrente termica) nell’istante

6 Cap. 1: Modellizzazione dello stress termomeccanico in uno slab di vetro ottico

t attraverso un elemento di superficie piana dS centrato in x; n indica ilversore normale a dS. j ha le dimensioni di una potenza su una superficie,e nel sistema internazionale e misurato in W/m2.

Dal primo principio della termodinamica, il bilancio energetico su unvolume di controllo V chiuso rispetto a scambi di materia e in assenza dicontributi di lavoro impone che

∂

∂t

∫

V

u(t, x) dV =

∫

V

w(t, x) dV −∫

S

j(t, x) · n dS , (1.3)

dove u(t, x) e la densita di energia termica interna, w(t, x) rappresenta ladensita di potenza termica assorbita/generata nel volume di interesse e l’in-tegrale di superficie al membro di destra rappresenta il flusso di calore uscen-te dal volume. Il segno negativo compare perche a una perdita di caloreattraverso S e associata una diminuzione dell’energia interna, mentre unagenerazione (positiva) di calore entro V determina un aumento dell’energiainterna.

Facendo uso del teorema della divergenza, la Eq. (1.3), valendo per ognivolume arbitrario V , puo essere posta in forma differenziale:

∂u(t, x)

∂t= w(t, x) − ∇ · j(t, x) . (1.4)

La Eq. (1.4) e l’equazione di continuita per la conservazione locale dell’e-nergia interna termica in un sistema chiuso in presenza di pozzi/sorgentiinterne di calore.

Almeno in prima approssimazione, e lecito aspettarsi una risposta linea-re della corrente di conduzione termica j(t, x) dal gradiente del campo ditemperatura T (t, x) secondo la relazione

j(t, x) = −κ∇T (t, x) , (1.5)

dove κ e un parametro, eventualmente dipendente dalla posizione, noto come“conducibilita termica”. Le unita di misura di κ sono W · m−1 · K−1. Inmezzi anisotropi, κ diventa un tensore, ma noi qui restringiamo la nostraattenzione ai vetri amorfi e quindi la descrizione con una quantita scalare esufficiente.

Introducendo il calore specifico a volume costante cV , si ottiene infinel’equazione di diffusione del calore in forma differenziale

cV ρ∂T (t, x)

∂t= w(t, x) + ∇ · (κ∇T (t, x)) , (1.6)

dove ρ e la densita di massa. Per mezzi omogenei, cV , ρ e κ hanno lo stessovalore in tutti i punti e la Eq. (1.6) si riduce a

∂T (t, x)

∂t=

1

cV ρw(t, x) + D∇2T (t, x) , (1.7)

1.3 Modellizzazione dell’assorbimento della radiazione laser 7

dove D = k/(cV ρ) e la costante di diffusione termica, espressa in m2 · s−1.Ulteriori effetti non-lineari possono insorgere se la capacita termica e laconducibilita termica variano con la temperatura.

La Eq. (1.6) si puo risolvere considerando condizioni iniziali e al con-trono assegnate. Condizioni al contorno tipiche consistono nell’assegnare ilcampo di temperatura sulla superficie in ogni istante di tempo (condizionedi Dirichlet) oppure il flusso di calore attraverso la superficie (condizione diNeumann) o condizioni miste (di Robin).

Nel seguito, ci interessa risolvere questa equazione per un parallelepipedodi dimensioni assegnate, inizialmente in equilibrio termico con l’ambienteesterno, cV , ρ e κ sono assegnate e si assume che la condizione al contorno siadel tipo di Dirichlet, con temperatura al bordo costante e uniforme su tuttala superficie, e pari alla temperatura ambiente. Condizioni al contorno cheincorporano anche il trasferimento radiativo all’interfaccia sono accennatein sezione 3.2. Il termine di assorbimento della radiazione w e descritto nellasezione successiva.

1.3 Modellizzazione dell’assorbimento della radia-

zione laser

Per la modellizzazione della sorgente di calore considero un fascio laser chesi propaga lungo z nel modo fondamentale gaussiano, e per questo motivosi puo descrivere l’irradianza nel vetro per mezzo della formula seguente:

I(ρ) = I0 exp(

−2ρ2/w20

)

, (1.8)

dove ρ e la distanza radiale dall’asse del fascio, I0 e l’irradianza di piccosull’asse del fascio e w0 e lo spot size. In questa formulazione e nel seguitosi trascura la divergenza del fascio rispetto alla distanza di propagazione.Inoltre, non si prende in considerazione la diminuzione dell’irradianza dipicco con la profondita. Se l’assorbimento di potenza dal fascio e lineare, ladensita di potenza di calore generata internamente e w(t, x) ∝ I(t, x).

In PVLAS, il fascio laser ha una potenza di circa 1 W, tuttavia in corri-spondenza dello strato interferenziale degli specchi l’assorbimento comples-sivo e di quattro/cinque ordini di grandezza maggiore perche occorre tenereconto delle riflessioni multiple che avvengono all’interno del Fabry-Perot.Questo significa che l’assorbimento nello specchio interessa perlopiu lo stra-to interferenziale alla superficie, e siccome lo spessore caratteristico di questostrato e circa 10 µm mentre lo spessore complessivo del vetro e dell’ordinedel centimetro, in prima approssimazione si puo trattare il problema comebidimensionale. Lo spot-size del fascio misura circa 2 mm.

Supponiamo che l’irradianza del fascio non cambi nel tempo e che l’as-sorbimento nel vetro sia costante. In assenza di variazioni nella posizionedell’asse del fascio laser, il campo termico evolve allora fino a raggiungere

8 Cap. 1: Modellizzazione dello stress termomeccanico in uno slab di vetro ottico

x

y

z

Figura 1.3. Illustrazione del fascio laser che attraversa il vetro dell’ellissometro.

una configurazione finale stazionaria, caratterizzata da un flusso di calorecostante dalle regioni piu calde a quelle piu fredde compensato dal calo-re assorbito dal fascio. In realta, la posizione dello spot del laser cambianel tempo, sia con una componente deterministica dovuta alla meccani-ca dell’apparato, sia con una componente stocastica dovuta all’inevitabilepresenza di fluttuazioni, delle quali e importante tenere conto nella nostrasimulazione.

Nel caso di PVLAS la presenza di questi movimenti del fascio si osservachiaramente dagli spettri del segnale riprodotti in Figura 1.4. Le diver-se componenti deterministiche che partecipano al moto dello spot del laserappaiono nello spettro come picchi di risonanze meccaniche, se ne posso-no individuare parecchie, mentre la componente stocastica si riflette in unelevato rumore di fondo.

In questa tesi modellizzo la posizione del fascio per mezzo di un’equa-zione di Langevin: questo significa che la posizione del centroide del fascior(t) = (x(t), y(t)) evolve secondo l’equazione di un oscillatore armonicobidimensionale, con coefficiente di smorzamento Γ. Il contributo forzanteadditivo e rappresentato da un vettore bidimensionale ξ(t) di processi sto-castici stazionari, indipendenti, bianchi, Gaussiani, a media nulla e densitaspettrale unitaria:

d2r(t)

dt2= −ω2

0r(t) − Γdr(t)

dt+ Λξ(t) , (1.9)

dove Λ e una matrice 2× 2 che pesa i contributi delle diverse componenti dirumore. Qui e nel seguito si considera per semplicita oscillatori indipendentisui due assi, e quindi Λ e diagonale; chiamiamo Λx e Λy i due autovalori.

1.3 Modellizzazione dell’assorbimento della radiazione laser 9

0 0.5 1 1.5 2

-100

-80

-60

-40

(a)

0 0.2 0.4 0.6 0.8 1

-100

-80

-60

-40

(b)

506 507 508 509 510 511

-160

-140

-120

-100

-80

(c)

Figura 1.4. (a) Spettro del segnale di PVLAS tra 0 e 2 kHz. Il picco alla frequenzadi 1.012 kHz, ben visibile nello spettro, e dovuta al modulatore (seconda armonica).La figura (b) mostra un primo ingrandimento tra 0 e 1 kHz, sono chiaramente visibilile side-bands a ±70 Hz dal picco di prima armonica del modulatore (a 506 Hz), einoltre si puo notare la presenza di molte risonanze meccaniche spurie non ancoraidentificate. La figura (c) infine mostra lo spettro con un livello di dettaglio ancoramaggiore, relativamente alla finestra spettrale tra 506 e 511 Hz circa; e visibile inquesto caso il picco di prima armonica a (506 + 0.3)Hz, ma non si ha riscontrodel picco previsto a circa (506 ± 0.6)Hz (qui 0.3 Hz e la frequenza di rotazione delmagnete), se c’e e mascherato dall’elevato livello di rumore; il picco largo intorno a510 Hz e dovuto a risonanza meccanica.

10 Cap. 1: Modellizzazione dello stress termomeccanico in uno slab di vetro ottico

Rimane da assegnare un valore agli elementi di Λ in modo da riprodurrela fluttuazione quadratica media dello spostamento sperimentalmente osser-vata. A tal fine, utilizzo la relazione di fluttuazione-dissipazione per questaequazione di Langevin. Per semplicita, considero il solo spostamento lungox, il caso in y e analogo.

Prendendo la trasformata di Fourier della Eq. (1.9) si trova

−ω2X(ω) = −ω20X(ω) − iωΓX(ω) + ΛxN(ω) , (1.10)

da cui

X(ω) =Λx

(ω20 − ω2) + iωΓ

N(ω) . (1.11)

Utilizzando quindi il teorema di Wiener-Kintchin, e la definizione di densitaspettrale, si trova che la fluttuazione quadratica media e

〈x2〉 =1

2πΛ2

x

∫ +∞

−∞

1

(ω20 − ω2)2 + ω2Γ2

dω (1.12)

Quest’ultimo integrale puo essere risolto in maniera esatta, ad esempio conla teoria dei residui. Una stima si puo ottenere piu velocemente come se-gue. La funzione integranda ha due massimi locali in corrispondenza diω = ±

√

ω20 − (Γ2/2), purche Γ2 < 2ω2

0;in questo caso, essa puo riscriversiapprossimativamente come somma di due Lorentziane centrate in ±ω0:

1

(ω20 − ω2)2 + Γ2ω2

=1

(ω0 − ω)2(ω0 + ω)2 + Γ2ω2

≈ 1

4ω20(ω0 − ω)2 + Γ2ω2

0

+1

4ω20(ω0 + ω)2 + Γ2ω2

0

.

Si ha quindi

⟨

x2⟩

≈ 1

2π

2Λ2x

ω20

∫ +∞

0

1

4(ω0 − ω)2 + Γ2dω ≈ 1

2Γ

Λ2x

ω20

,

da cui segue la relazione di fluttuazione-dissipazione

Λx = ω0

√2Γ

√

〈x2〉 = ω0

√2ΓxRMS . (1.13)

In PVLAS ad esempio, si osservano due frequenze di risonanza dominantiin prossimita della regione di interesse fisico (oltre a diverse risonanze mecca-niche minori), i valori numerici corrispondenti sono riportati in Tabella 1.1.L’oscillazione a 2π × 70 ≈ 440 radHz ad esempio, si ritrova (spostata a fre-quenza piu alta per effetto della modulazione) nello spettro in Figura 1.4(b),nelle bande laterali intorno alla portante di modulazione a 506 Hz.

L’equazione di Langevin (1.9) e un esempio di equazione differenzialestocastica, e per l’integrazione numerica di questo tipo di equazioni sonodisponibili diversi algoritmi numerici derivati piu o meno direttamente daanaloghi algoritmi per le equazioni differenziali ordinarie: alcuni di questialgoritmi sono presentati nella sezione successiva.

1.4 Integratori numerici per equazioni differenziali stocastiche 11

Tabella 1.1. Valori numerici per la modellizzazione stocastica del fascio diPVLAS.

(a) Componente a corto range

ω0 ≈ 440 rad · Hzγ ≈ 30 s−1

xRMS ≈ 10−9 myRMS ≈ 10−9 m

(b) Componente a lungo range

ω0 ≈ 0.006 rad · Hzγ ≈ 3 s−1

xRMS ≈ 10−4 myRMS ≈ 10−4 m

1.4 Integratori numerici per equazioni differenzia-

li stocastiche

In questa sezione vengono presentati brevemente due algoritmi per l’inte-grazione numerica di equazioni differenziali stocastiche e come si fa pratica-mente sempre in fisica viene utilizzata la rappresentazione di Ito. Si trattadi schemi espliciti alle differenze finite con passo di integrazione singolo adampiezza fissata. Le loro proprieta sono solo riassunte, per ulteriori dettaglisi vedano le referenze [7, 8, 11].

L’equazione di Langevin (1.9) fa parte di una classe di equazioni chehanno la forma seguente:

dX(t) = f(X(τ)) dτ + g(X(τ)) dW (τ) , (1.14)

con condizioni iniziali assegnate X(t0) = X. Queste equazioni diffenziali sipossono integrare formalmente nel modo conseguente:

X(t + ∆t) = X(t) +

∫ t+∆t

t

f(X(τ)) dτ +

∫ t+∆t

t

g(X(τ)) dW (τ) , (1.15)

nell’intervallo [t, t + ∆t], e quindi se si assume la quasi-stazionarieta dellafunzione f nell’intervallo e si scrive esplicitamente l’incremento del processostocastico nel secondo integrale si ottiene la formula Euler-Maruyama

X(t + ∆t) = X(t) + f(X(t))∆t + g(X(t)) [W (t + ∆t) − W (t)] . (1.16)

Questa formula puo essere applicata ricorsivamente a partire dalle condizioniiniziali ottenendo

Xn+1 = Xn + f(Xn) [tn+1 − tn] + g(Xn) [W (tn+1) − W (tn)] , (1.17)

con tn+1 = tn + ∆t. Per le proprieta di un processo scalare standard diWiener si ha ∆Wn ∼

√∆tN(0, 1), i.e., gli incrementi ∆Wn sono variabili

casuali statisticamente indipendenti tra loro e con distribuzione normale amedia zero e varianza ∆t. Se g = 0, si ritrova l’algoritmo di Eulero perequazioni differenziali ordinarie.

12 Cap. 1: Modellizzazione dello stress termomeccanico in uno slab di vetro ottico

Nel caso dell’algoritmo di Eulero per equazioni differenziali ordinarie,non e difficile dimostrare che l’errore di troncamento globale e O(∆t). Piudifficile e valutare una stima (o almeno una qualche limitazione superiore)dell’errore commesso nel caso stocastico, perche non c’e un’unica soluzioneesatta X(t) con cui confrontare anche solo teoricamente il risultato, bensıX(t) e un processo stocastico di cui si considerano diverse realizzazioni. Echiaro che i criteri per stabilire la convergenza dell’algoritmo devono com-portare un qualche processo di media. Senza entrare troppo nel dettaglio,qui si ricorda solo il cosiddetto criterio di convergenza “forte”: si dice chel’algoritmo di integrazione per SDEs ha ordine di convergenza forte p seesiste una costante positiva K tale che

⟨

|Xn − X(tn)|⟩

≤ K∆tp , (1.18)

per ogni fissato tn ∈ [t0, t0 + T ], le parentesi ad angolo denotano la mediad’ensamble. Un altro criterio di convergenza diffuso e quello “debole”, chequi non e trattato. Si puo dimostrare che il metodo di Euler-Maruyamaha ordine di convergenza forte 1/2, se si assume che le funzioni f e g sianoregolari.

L’algoritmo di Milstein consiste nel modificare la ricorsione (1.17) con

Xn+1 = Xn + f(Xn)∆t + g(Xn) [W (tn+1) − W (tn)]

+1

2g(Xn)g′(Xn)

[W (tn+1) − W (tn)]2 − ∆t

, (1.19)

E importante menzionare l’algoritmo di Milstein perche si puo dimostrareche ha ordine di convergenza forte 1, a patto che le funzioni f e g soddisfinoopportune condizioni di regolarita. Nel caso della Eq. (1.9) che ci interessa,g′ = 0 e quindi in questo caso i due algoritmi forniscono gli stessi risultati

Altri algoritmi standard per la risoluzione numerica di equazioni dif-ferenziali ordinarie non si estendono in modo altrettanto semplice al casostocastico, al riguardo si veda ad esempio [11].

Capitolo 2

Cenni di analisi agli elementi

finiti

Per il calcolo della temperatura all’interno dello specchio modellizzato comespiegato nel capitolo 1, ho deciso di utilizzare l’analisi agli elementi finiti cheincorpora una varieta di approcci numerici per la risoluzione approssimatadi equazioni alle derivate parziali con condizioni al contorno assegnate.

Questa sezione offre una rassegna introduttiva delle basi dei metodi FEAcon l’intento di presentare una panoramica delle idee piu importanti. Unriferimento tipico sui metodi FEA e [28]. Gli aspetti computazionali sonotrattati, ad esempio, in [22]; si veda anche [13]. Le basi matematiche sonosolo accennate in questo capitolo, si trovano sviluppate estesamente in [5].

2.1 Metodo degli elementi finiti per la conduzione

termica stazionaria

In questa sezione, mi occupo di esporre un metodo di soluzione a elementifiniti per un problema di conduzione termica in regime stazionario, limitan-domi in particolare al metodo di Galerkin bidimensionale. L’estensione aproblemi di conduzione dipendenti dal tempo e svolta nella prossima sezione.

L’equazione stazionaria del calore e

−∇ · [κ(x)∇T (x)] = f(x) , (2.1)

e per semplicita considero condizioni al contorno di Dirichlet. Ovviamente,si puo generalizzare il discorso anche ad altri tipi di condizioni al contorno.

Come noto, l’analisi alle differenze finite consiste sostanzialmente neldiscretizzare lo spazio con un reticolo di punti e sostituire le derivate checompaiono nelle Eqs. (2.1) con rapporti incrementali finiti sul dominio di-scretizzato. Procedendo ricorsivamente a partire dalle condizioni di bordo

13

14 Cap. 2: Cenni di analisi agli elementi finiti

assegnate, e possibile pervenire a un’approssimazione del valore della solu-zione esatta del problema (2.1) in corrispondenza dei punti discreti. Per unarassegna dei metodi alle differenze finite, si veda Press e altri [19, cap. 20].

Per contro, l’idea centrale dell’analisi agli elementi finiti e quella di ap-prossimare la soluzione cercata con uno sviluppo finito in termini di uninsieme di funzioni a supporto locale e linearmente indipendenti tra loro,cioe

T (x) ≈ T (x) =M∑

j=1

Tjφj(x) , (2.2)

dove T (x) denota la soluzione esatta. Soprattutto in ingegneria, le funzioniφj(x) sono note con il nome di “funzioni di forma”. L’obiettivo consistenel determinare il valore dei coefficienti incogniti Tj in modo tale che T (x)approssimi “quanto meglio” T (x).

La scelta delle funzioni di forma e un passo importante dell’analisi aglielementi finiti, e tra poco presento un esempio concreto di come si possanoscegliere queste funzioni. Prima pero occorre precisare che cosa significa cheT approssima “quanto meglio” T . L’idea e la seguente. E lecito aspettarsiche, in genere, T non coincida con T ovunque. In altre parole, se sostitui-sco T nell’equazione differenziale Eq. (2.1), quest’ultima in genere non sarasoddisfatta esattamente. Lo scarto e detto “residuo” R, e altro non e che

R = ∇ ·

[

κ(x)∇T (x)]

+ f(x) .

Cambiando i coefficienti Tj , in genere cambia T e quindi anche R. Se fosseT = T , allora R = 0 ovunque; piu in generale, ci si accontenta di determinarei coefficienti Tj nella Eq. (2.2) in modo che il residuo sia “piccolo”.

Che cosa significa che R e “piccolo”? Qui mi limito a considerare ilcosiddetto “metodo dei residui pesati”, che consiste nell’imporre che le mediepesate di R con M funzioni-peso Wi(x) linearmente indipendenti tra loro,siano nulle:

∫

RWi dS = 0 , i = 1, . . . , M . (2.3)

Qui e nel seguito gli integrali si intendono estesi a tutta la regione di in-teresse, in tre dimensioni gli integrali di superficie sono sostituiti da inte-grali di volume. Matematicamente, la Eq. (2.3) esprime l’ortogonalita delcampo residuo con le funzioni di test. A parole, essa equivale sostanzial-mente a richiedere che la Eq. (2.1) sia soddisfatta in media sul dominio diintegrazione.

Si parla di “metodo di Galerkin” quando si scelgono Wi = φi per ognii. Se Wi 6= φi, si parla di “metodo di Petrov-Galerkin”. Per semplicita, nelseguito considero solo il metodo di Galerkin.

Come si possono incorporare le condizioni di bordo? La condizione diDirichlet puo essere soddisfatta automaticamente se lo sviluppo e sostituito

2.1 Metodo degli elementi finiti per la conduzione termica stazionaria 15

con

T (x) = Ψ(x) +M∑

j=1

Tjφj(x) , (2.4)

dove Ψ(x) e una qualunque funzione che soddisfi le condizioni al contorno, erichiedendo che tutte le φj(x) si annullino identicamente al bordo. Esistonostrategie per incorporare anche altri tipi di condizioni al contorno, ma persemplicita qui non sono esposte.

Inserendo lo sviluppo (2.4) nelle Eqs. (2.3) si ottiene

M∑

j=1

∫

κ(x)∇φi(x) · ∇φj(x) dS

Tj =

∫

f(x)φi(x) dS

−∫

κ(x)∇φi(x) · ∇Ψ(x) dS ,

(2.5)

per ogni i = 1, . . . , M . L’insieme delle Eqs. (2.5) per ogni i costituisce unsistema di M equazioni algebriche accoppiate nelle incognite Tj , che si puoriscrivere in forma piu compatta introducendo la notazione matriciale

KT = c , (2.6)

dove K e la “matrice di rigidezza globale” di dimensioni M ×M ed elementi

kij =

∫

κ(x)∇φi(x) · ∇φj(x) dS , (2.7)

T e il vettore colonna M -dimensionale dei coefficienti incogniti, con compo-nenti T1, . . . , TM , infine c e il vettore colonna M -dimensionale dei terminiinomogenei, le cui componenti sono

ci =

∫

f(x)φi(x) dS −∫

κ(x)∇φi(x) · ∇Ψ(x) dS . (2.8)

Si osservi che le condizioni di bordo sono state trasferite nel vettore c comeparte dell’inomogeneita. La determinazione dei coefficienti Tj risulta dallarisoluzione algebrica del sistema lineare: T = K−1c.

Una questione e rimasta aperta: come si scelgono le φj(x)? E que-sto un punto qualificante dell’analisi agli elementi finiti. Le tecniche FEAsi distinguono rispetto ad altri metodi spettrali per la scelta di funzionidi base a supporto locale: il dominio di integrazione viene partizionato insotto-dominii non sovrapposti (gli “elementi finiti”) tipicamente di formageometrica semplice (ad esempio, poligoni in due dimensioni e poliedri intre dimensioni); per ogni elemento, si definiscono dei “punti nodali”, i puntinodali possono coincidere con i vertici degli elementi ma ce ne possono es-sere anche altri disposti lungo i lati degli elementi o addirittuta all’internodi essi; infine, per ogni nodo di ogni elemento si introduce una funzione di

16 Cap. 2: Cenni di analisi agli elementi finiti

forma (generalmente, polinomiale a tratti) a supporto solo in una piccolaregione attorno al nodo corrispondente, e tale da interpolare la variabile dicampo sui singoli elementi contenenti quel nodo.

Senza insistere sulla generalita, viene illustrato il caso di elementi pianitriangolari per un problema di conduzione in due dimensioni. Naturalmente,il discorso puo essere generalizzato a casi piu complicati. I nodi sono ivertici del triangolo e come funzioni di forma si considerano polinomi diinterpolazione lineare a tratti e con supporto locale.

Si consideri un singolo elemento triangolare, i cui tre vertici (nodi) sianoindividuati rispettivamente dai vettori di coordinate piane (globali) xj =(xj , yj), dove j = 0, 1, 2 e l’indice locale di nodo. Per ciascuno dei tre

vertici, si introduce una funzione di forma locale φ(e)j (x) lineare e a supporto

sul singolo elemento. Le φ(e)j sono convenientemente scelte in modo da

soddisfare la condizione ai nodi

φ(e)j (xi) = δi,j , (i, j = 0, 1, 2) (2.9)

dove δi,j e l’usuale simbolo di Kronecher, avendo per ora l’accortezza diignorare la richiesta di annullamento delle funzioni di forma sulla frontieraessenziale. Le condizioni al contorno essenziali potranno essere ripristinatein un secondo momento come spiegato piu avanti. Lo sviluppo (2.4) e oraesteso a tutte le funzioni di forma definite per ogni nodo di ogni elementodel dominio. La proprieta (2.9) non e obbligatoria, e una scelta. Rispettoa questa scelta, i coefficienti Tj nello sviluppo (2.2)-(2.4) si prestano a unasemplice interpretazione: essi coincidono in questo caso con i valori puntualiassunti dalla variabile di campo in corrispondenza dei nodi.

La formula esplicita di φ(e)j puo essere ricavata per via analitica in termini

delle coordinate (globali) dei vertici: ad esempio, sia

φ(e)0 (x) = a0 + b0x + c0y ,

allora la richiesta (2.9) si traduce nel sistema

a0 + b0x0 + c0y0 = 1 ,

a0 + b0x1 + c0y1 = 0 ,

a0 + b0x2 + c0y2 = 0 ,

la risoluzione conduce a

a0 =1

∆det

(

x1 y1

x2 y2

)

, b0 = − 1

∆det

(

1 y1

1 y2

)

, c0 =1

∆det

(

1 x1

1 x2

)

,

dove

∆ = det

1 x0 y0

1 x1 y1

1 x2 y2

,

2.1 Metodo degli elementi finiti per la conduzione termica stazionaria 17

0 0.5 1 1.5 2 2.5 3 3.5 4x

0

0.5

1

1.5

2

2.5

3

3.5

4

y

(a)

0 0.5 1 1.5 2 2.5 3 3.5 4

x

0 0.5 1 1.5 2 2.5 3 3.5 4

y

0

0.2

0.4

0.6

0.8

1

(b)

Figura 2.1. (a) Un esempio di suddivisione di un dominio quadrato in elementidi forma triangolare. Gli elementi triangolari sono stati ottenuti a partire da unagriglia cartesiana equispaziata. I nodi coincidono con i vertici dei triangoli. Sonoevidenziati in rosso gli elementi che hanno il punto di coordinate (2, 2) come nodo.(b) Funzione di forma globale per il nodo di coordinate (2, 2), la funzione e linearesu ciascun elemento che contiene (2, 2) come nodo, ed e identicamente nulla su tuttigli altri elementi.

il cui valore assoluto corrisponde a due volte l’area del triangolo. Il calcoloe analogo per le funzioni di forma relative agli altri due vertici.

Per maggiore concretezza, consideriamo un esempio particolare (si ve-da la Figura 2.1). Il dominio di interesse e il quadrato [0, 4] × [0, 4]. Lasuddivisione del dominio in elementi triangolari e stata eseguita costruen-do dapprima una griglia cartesiana equispaziata in modo da ottenere 4 × 4quadrati, poi ciascun quadrato e stata a sua volta suddiviso lungo una dia-gonale in modo da ottenere due elementi triangolari, come raffigurato inFigura 2.1(a). In rosso sono indicati tutti gli elementi triangolari che hannoil punto di coordinate (2, 2) come nodo. Il grafico della funzione di formaglobale corrispondente al nodo (2, 2) e riportato in Figura 2.1(b). Si vedeche essa e lineare a tratti sui sei elementi che contengono (2, 2) come nodo,vale 1 nel punto (2, 2) ed e identicamente nulla su tutti gli altri elementi.

Gli integrali nelle Eqs. (2.7)-(2.8) possono essere scritti come somme dicontributi sui singoli elementi, ad esempio per la matrice di rigidezza si ha

ki,j =

Ne∑

e=1

k(e)i,j , k

(e)i,j =

∫

e

κ(x)∇φi(x) · ∇φj(x) dS ,

dove la sommatoria e estesa a tutti gli Ne elementi del dominio mentre l’in-

tegrale e esteso al solo elemento di indice e. Chiaramente, i coefficienti k(e)i,j

18 Cap. 2: Cenni di analisi agli elementi finiti

sono diversi da zero solo se i e j sono nodi appartenenti a un medesimoelemento (altrimenti φi o φj risultano identicamente nulle entro l’elemento).Questo ha due implicazioni: il sistema lineare globale puo essere costruito co-me assemblaggio di contributi sui singoli elementi, e la matrici di rigidezzaglobale cosı assemblata e sparsa (i.e., con tanti zeri), che computazional-mente e molto vantaggioso perche esistono algoritmi molto efficienti per larisoluzione di sistemi lineari tramite l’inversione di matrici sparse.

Le formule per i contributi locali per la matrici di rigidezza si possonocalcolare analiticamente in questo caso, dato la forma semplice delle funzionidi forma e della geometria degli elementi. Ad esempio, per la matrice dirigidezza si trova

k(e)i,j =

∫

e

κ(x)∇φ(e)i (x) · ∇φ

(e)j (x) dS

≈ κ(x(e)c ) (bibj + cicj) Se ,

dove Se e l’area del triangolo e x(e)c indica la posizione del baricentro dell’e-

lemento. Si procede in modo analogo per il vettore dei termini noti.Finora le condizioni al contorno sono state ignorate. Una soluzione di

facile implementazione per ripristinare le condizioni al contorno essenziali ela seguente: per ogni nodo di Dirichlet, si azzera la corrispondente riga nellamatrice di rigidezza globale, tranne sulla diagonale dove si impone il valoreunitario, e si forza la corrispondente componente del vettore dei termini ino-mogenei ad assumere proprio il valore prescritto dalla condizione di Dirichletin quel nodo. Questa strategia e generalmente semplice da implementare,anche se non e computazionalmente efficiente perche introduce nel sistemalineare delle incognite aggiuntive che non sono vere incognite e potrebberoessere rimosse gia in partenza perche il loro valore e noto.

2.2 Problemi dipendenti dal tempo al prim’ordi-

ne.

Questa sezione si occupa di presentare una strategia (il “θ-metodo”) pergestire problemi di valori al contorno dipendenti dal tempo al prim’ordine,nell’ambito di un’analisi agli elementi finiti con il metodo di Galerkin. Unadiscussione del θ-metodo e reperibile ad esempio in [21].

In principio, non c’e nessuna ragione per cui il metodo dei residui pe-sati visto alla sezione precedente non possa essere automaticamente estesoper incorporare oltre alle dimensioni spaziali anche la dimensione tempo-rale. Tuttavia, l’approccio piu usuale (anche se certamente non il solo)con cui vengono affrontati problemi dipendenti dal tempo consiste nel trat-tare la dipendenza temporale separatamente, mediante il ricorso a schemidi discretizzazione alle differenze finite, mentre l’integrazione spaziale vienesvolta tramite il metodo degli elementi finiti.

2.2 Problemi dipendenti dal tempo al prim’ordine. 19

Consideriamo l’equazione di diffusione del calore:

∇ · [κ(x)∇T (t, x)] + f(t, x) = ρ(x)cV (x)∂T (t, x)

∂t.

Si intendono assegnate la condizione iniziale T (t0, x) = T0(x) all’istanteiniziale t0 e condizioni al contorno di Dirichlet a ogni istante. Il discorsopuo essere esteso anche ad altri tipi di condizioni al contorno, ma qui persemplicita si considerano solo le condizioni di Dirichlet.

Si procede dapprima alla discretizzazione del problema spaziale. A ogniistante di tempo, la variabile di campo puo essere approssimata espandendosulla base delle funzioni di forma, cioe

T (t, x) ≈ T (t, x) = Ψ(t, x) +M∑

j=1

Tj(t)φj(x) , (2.10)

dove le φj(x) sono sempre le stesse a ogni t, mentre i coefficienti dello svilup-po cambiano nel tempo perche il campo di temperatura in genere sara diver-so a istanti di tempo diversi. Ψ(t, x) e una qualunque funzione che riproducaa ogni istante di tempo le corrette condizioni al contorno essenziali.

Usando il metodo di Galerkin, si ottiene

∫

∇φi(x) ·

[

κ(x)∇T (t, x)]

dS =

∫

f(t, x)φi(x) dS

−∫

κ(x)∇φi(x) · ∇Ψ(t, x) dS −∫

ρ(x)cV (x)∂T (t, x)

∂tφi(x) dS .

(2.11)

Inserendo lo sviluppo (2.10), e siccome

∂T (t, x)

∂t=

∂Ψ(t, x)

∂t+

M∑

j=1

(

dTj(t)

dt

)

φj(x) ,

si ottiene infine il sistema di M equazioni algebriche accoppiate

M∑

j=1

∫

Ω

(

κ(x)∇φi(x) · ∇φj(x) dS

Tj(t)

+M∑

j=1

∫

Ωρ(x)cV (x)φi(x)φj(x) dS

dTj(t)

dt=

∫

f(t, x)φi(x) dS −∫

κ(x)∇φi(x) · ∇Ψ(t, x) dS

−∫

ρ(x)cV (x)∂Ψ(t, x)

∂tφi(x) dS ,

(2.12)

20 Cap. 2: Cenni di analisi agli elementi finiti

per ogni i = 1, . . . , M , che in forma matriciale si puo scrivere come

MT (t) + KT (t) = c(t) , (2.13)

dove T (t) e il vettore dei coefficienti di espansione (relativi alla distribuzio-ne spaziale al tempo t), T = dT (t) / dt, M e la matrice delle masse, conelementi

mi,j =

∫

ρ(x)cV (x)φi(x)φj(x) dS , (2.14)

K e la matrice di rigidezza gia definita nella Eq. (2.7) e c(t) e il vetto-re dei termini inomogenei, le cui componenti si ricavano immediatamenteconfrontando con la Eq. (2.12).

Gli integrali nella Eq. (2.12) sono stati indicati come estesi all’interodominio di interesse, ma in pratica le funzioni di forma sono a supportolocale e pertanto gli integrali si riducono a somme di contributi sui singolielementi. Inoltre, se gli elementi sono sufficientemente piccoli allora κ, ρ, cV

e f si possono considerare costanti all’interno di ogni singolo elemento.Affinche il sistema (2.13) possa essere risolto rispetto al vettore dei coef-

ficienti T (t) ai diversi istanti di tempo, occorre ancora poter esprimere T (t)in funzione di T (t). A tal fine, si ricorre a uno schema iterativo alle diffe-renze finite nel dominio del tempo. In particolare, il θ-metodo consiste nelsostituire il sistema (2.13) con la relazione ricorsiva a tempi discretizzati

MT (tn+1) − T (tn)

∆t+K [θT (tn+1) + (1 − θ)T (tn)] = θc(tn+1)+(1−θ)c(tn) ,

ovvero(

M

∆tn+ θK

)

T (tn+1) =

(

M

∆t+ (θ − 1)K

)

T (tn)+θc(tn+1)+(1−θ)c(tn) ,

(2.15)su una griglia temporale tn+1 = tn + ∆tn. Partendo dalle condizioni inizialiassegnate, occorre risolvere a ogni passo la Eq. (2.15) rispetto al vettore dicoefficienti incognito T (tn+1), noto T (tn) al passo precedente. Nel seguito,si considera per semplicita una griglia equispaziata ad ampiezza costante∆t. In questo caso la matrice A = M

∆t+ θK e la stessa a ogni passo, e se

per risolvere il sistema lineare si fa uso di schemi diretti che comportanola fattorizzazione di A allora e conveniente fattorizzare la matrice una solavolta all’inizio del ciclo temporale.

Per θ = 0 si ottiene lo schema di Euler in avanti (esplicito), per θ = 1 siottiene lo schema di Euler all’indietro (implicito), per θ = 1/2 si ha lo schemadi Crank-Nicolson. Si dimostra che per θ 6= 1/2 l’errore di discretizzazionenel tempo e dell’ordine di ∆t, mentre il metodo di Crank-Nicolson e unoschema del second’ordine. Inoltre, per ogni θ ≥ 1/2 e possile mostrare cheil θ-metodo e incondizionatamente stabile (cioe e stabile per ogni scelta di∆t) mentre per θ < 1/2 il metodo e stabile condizionatamente.

Capitolo 3

Calcolo termico e risultati

preliminari

3.1 Implementazione

In questa sezione, espongo brevemente l’organizzazione del codice in C++che ho messo a punto per la risoluzione numerica del campo termico.

Per la discretizzazione spaziale, ho usato il metodo di Galerkin bidimen-sionale con elementi triangolari quadratici. Gli elementi triangolari sonostati costruiti a partire da una griglia rettangolare equispaziata, secondolo schema illustrato in Figura 3.1. Per l’evoluzione temporale, ho usato loschema di Eulero all’indietro (θ = 1).

La classe principale del programma e FEMSolver, che contiene tutti iparametri e le funzioni necessarie a integrare il problema di campo. Leseguenti funzioni sono private, e vengono chiamate dal costruttore:

• setNodePosition(): definisce le coordinate globali di tutti i nodi;

• setTriangularGrid(): definisce la numerazione globale degli elemen-ti e l’array di connettivita (che restituisce l’indice di nodo globalecorrispondente a un indice di nodo locale di un dato elemento);

• setBoundaryNodes(): stabilisce quali siano i nodi sul bordo.

La matrice di rigidezza globale e la matrice della masse sono assemblate dallafunzione assemble(). L’evoluzione temporale e implementata all’internodella funzione pubblica run(), che esegue una chiamata ad assemble() eil ciclo temporale. A ogni passo temporale, l’array dei termini noti vieneriassemblato, e si risolve numericamente il sistema lineare (2.12).

La risoluzione del sistema lineare a ogni passo temporale e stata esegui-ta con il classico algoritmo di decomposizione LU, implementato in una suavariante per matrici a bande [19, cap. 2.4.2]. Il calcolo della semi-larghezza

21

22 Cap. 3: Calcolo termico e risultati preliminari

0 1 2 3 4 5 6 7 8 9 10

11 12 13 14 15 16 17 18 19 20 21

22 23 24 25 26 27 28 29 30 31 32

33 34 35 36 37 38 39 40 41 42 43

44 45 46 47 48 49 50 51 52 53 54

55 56 57 58 59 60 61 62 63 64 65

66 67 68 69 70 71 72 73 74 75 76

77 78 79 80 81 82 83 84 85 86 87

88 89 90 91 92 93 94 95 96 97 98

99 100 101 102 103 104 105 106 107 108 109

110 111 112 113 114 115 116 117 118 119 120

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

Figura 3.1. Illustrazione di una semplice suddivisione di un dominio rettangolarein elementi di forma triangolare a partire da una griglia cartesiana con spaziaturauniforme. Ciascun elemento contiene sei nodi (elementi quadratici). E riportato unesempio di possibile numerazione globale degli indici di nodo e di elemento. Nellasimulazione si e adottatato un numero di punti nodali maggiore di quello illustratoqua. (Precisamente, la griglia utilizzara conteneva 69 × 69 punti.) Si puo pensareeventualmente di addensare il numero di elementi nelle regioni centrali, dove ilcampo termico varia piu rapidamente per effetto del laser.

3.1 Implementazione 23

di banda della matrice assemblata puo essere effettuato all’inizio, una vol-ta che siano state definite la numerazione globale degli indici e la matricedi connettivita. Dal momento che la matrice dei coefficienti non dipendedal tempo, la fattorizzazione LU puo essere eseguita una volta per tutteprima del ciclo temporale, con un grosso guadagno in termini computazio-nali. Gli algoritmi di inversione sono stati inglobati all’interno della classeBandedLUSolver, e si e provveduto a predisporre una classe BandedMatrix

che faccia da interfaccia per una generica matrice a bande, di dimensionenota e con larghezza di banda superiore e inferiore assegnate.

La sorgente di calore e stata implementata nel modo seguente. E sta-ta realizzata una classe base astratta HeatSource. La sua funzione inter-na getHeat(double aX, double aY, double aTime) e puramente virtua-le ed e pensata per restituire la potenza termica generata nel punto di coor-dinate aX e aY e al tempo aTime. Essendo puramente virtuale, la funzionegetHeat non e implementata nella classe HeatSource, si trova implementa-ta nelle varie sotto-classi derivate da HeatSource, ciascuna delle quali vuolerapppresentare una modalita differente di generazione di calore. La classeFEMSolver accede tramite puntatore a una classe derivata da HeatSource

che implementi una specifica modalita di generazione del calore. La pro-grammazione a oggetti incoraggia enormemente l’uso di funzioni virtuali eclassi base astratte, tuttavia si deve tenere presente che l’uso di funzioni vir-tuali puo essere poco efficiente, in primo luogo perche le funzioni virtuali nonpossono essere rese in-line, e in secondo luogo perche spesso il compilatorenon e in grado di ottimizzare al meglio sezioni del codice che contengonochiamate a funzioni virtuali.

La classe GaussianBeamHeatSource derivata da HeatSource implemen-ta nello specifico la generazione di calore dovuta ad assorbimento di po-tenza da un fascio gaussiano il cui asse e soggetto a fluttuazioni di posi-zione, come spiegato al Capitolo 1. Essa contiene (come puntatori) dueoscillatori di Langevin, descritti dalla classe Oscillator2D. La proceduraadottata per l’integrazione numerica dell’equazione stocastica di un oscilla-tore e il semplicissimo schema di Euler-Maruyama, ed e stato implementatodirettamente all’interno della classe oscillatore, precisamente nella funzio-ne integrate(). Piu in generale, si sarebbero potute realizzare due classibase astratte SDEProblem e SDESolver, la prima delle quali pensata comeinterfaccia per descrivere un qualunque problema a valori iniziali stocasti-co, mentre la seconda per implementare un singolo metodo di integrazionestocastica.

L’integrazione di equazioni stocastiche necessita di una sorgente di nu-meri casuali. La scelta adottata e quella, piuttosto tipica, di utilizzare deglialgoritmi di generazione di numeri pseudo-casuali. Variabili a distribuzionenon-uniforme sono ottenute mediante manipolazione di un campione pseudo-casuale estratto dalla distribuzione uniforme su un intervallo di riferimento,che di solito e l’intervallo [0, 1[ o ]0, 1[.

24 Cap. 3: Calcolo termico e risultati preliminari

Il generatore di numeri psuedo-casuali uniformi utilizzato e il generatorea ricorsione lineare multipla combinato MRG32k3a proposto da [17]. Perla generazione di variabili gaussiane, si e fatto uso del classico algoritmo diBox-Muller, per una descrizione del quale si consulti, e.g., [4] o Press e altri

[19, cap. 7]; si tratta di un algoritmo che oggi riveste importanza perlopiudidattica, un metodo piu competitivo e, e.g., quello suggerito in [19, p. 369].Si e provveduto comunque a predisporre una libreria di generatori uniformie non uniformi (le dichiarazioni sono poste in due header files diversi) inmaniera tale da agevolare la scelta di altre tecniche di generazione. Diver-samente da quanto suggerito in Press e altri [19, cap. 7], si e preferito nonimplementare i generatori non uniformi come sotto-classi derivate da quelliuniformi, in parte per praticita ma anche perche dal punto di vista logi-co un generatore non uniforme non e un esempio particolare di generatoreuniforme! Casomai, si sarebbe potuto implementare i generatori uniformicome sotto-classe di quelli non uniformi. Per alcune considerazioni sullaprogrammazione orientata a oggetti nell’ambito dei generatori di numeripseudo-casuali si consulti [9, 10].

3.2 Risultati preliminari

L’ampiezza di un singolo passo temporale e stata fissata a 0.0001 s, visti ivalori di frequenza dell’oscillatore a corto raggio (si veda Tabella 1.1). Ilcalcolo termico viene rieseguito a intervalli di 50 passi temporali. A ognipasso, il programma registra le posizioni degli oscillatori e del centroidedel fascio su un file di testo, e a ogni calcolo termico i valori nodali dellasoluzione sono pure salvati su file. Le posizioni dei nodi sono salvate soloall’istante iniziale, dal momento che sono sempre le stesse anche agli istantisuccessivi.

La posizione della componente a corto range dell’oscillazione del fascio eriprodotta in Figura 3.2 dopo 4 s di simulazione. Il campo termico dopo 4 sdi simulazione e riprodotto in Figura (3.3). Si e posta convenzionalmentea zero la temperatura ambiente. La grafica nella figura non e purtrop-po ottimale, sembra suggerire un comportamento oscillatorio del campo ditemperatura avvicinandosi ai bordi; non ci sarebbe da stupirsi, oscillazio-ni spurie di natura puramente algoritmica non sono rare in questi schemidi integrazione. D’altro canto, se si prende una sezione a y costante, adesempio y = 0.01 m come in Figura 3.4, si vede che non e presenta alcuncomportamento oscillatorio.

L’analisi svolta finora non e completa. Alcuni aspetti da migliorare sonoi seguenti. Per prima cosa, l’implementazione e da estendersi al caso diun vetro tridimensionale. Dal punto di vista delle tecniche numeriche, unproblema tridimensionale non comporta grosse complicazioni rispetto a unproblema in due dimensioni, i concetti rimangono infatti gli stessi, tuttavia

3.2 Risultati preliminari 25

-2e-09

-1.5e-09

-1e-09

-5e-10

0

5e-10

1e-09

1.5e-09

2e-09

-2e-09 -1.5e-09 -1e-09 -5e-10 0 5e-10 1e-09 1.5e-09 2e-09

y (m

)

x (m)

Position of the short-range oscillator

time = 4 sec.

Figura 3.2. Simulazione della componente a corto range del moto oscillatorio delfascio di PVLAS, limitatamente ai primi 4 s di simulazione, ottenuto con un passotemporale di 0.0001 s.

0 0.005

0.01 0.015

0.02 0

0.005

0.01

0.015

0.02

0

0.002

0.004

0.006

0.008

0.01

0.012

Temperature field

time = 4 sec.

x (m)

y (m)

Figura 3.3. Profilo del campo termico nel vetro in approssimazione 2D e per 4 sdi simulazione.

26 Cap. 3: Calcolo termico e risultati preliminari

0

0.001

0.002

0.003

0.004

0.005

0.006

0.007

0.008

0.009

0.01

0 0.005 0.01 0.015 0.02

x (m)

time = 4 sec.

Figura 3.4. Sezione a y = 0.01 m (meta della profondita in y) del campo termico,le diverse curve si riferiscono ai diversi passi temporali fino alla curva a temperaturamaggiore che corrisponde a 4 s di simulazione.

3.2 Risultati preliminari 27

alcuni dettagli di implementazione cambiano parecchio, e il debugging delcodice diventa piu delicato.

Si puo pensare di includere nella simulazione anche uno stress puramentemeccanico dovuto alla pressione di radiazione del fascio alla superficie, ostress meccanici espliciti dovuti al supporto. Questa richiesta impone nonsolo alcune modifiche del codice, ma anche la modellizzazione preliminaredi questi stress.

Le condizioni al contorno possono essere modificate in modo da inclu-dere piu realisticamente anche l’irraggiamento all’interfaccia. L’espressioneper la densita superficiale di potenza netta scambiata radiativamente conl’ambiente esterno, in tutte le direzioni e a tutte le frequenze, in regime diquasi-equilibrio, si trova dalla ben nota legge di Stefan-Boltzmann [e.g., 27]ed e

wS(t, x) = σ(

eT 4S(t, x) − ebaT 4

0

)

, (3.1)

dove σ = 5.670 × 10−8 W · m−2 · K−4 e la costante di Stefan-Boltzmann, e,a ed eb sono rispettivamente l’emissivita e il coefficiente di assorbimento delcorpo e l’emissivita dell’ambiente esterno, TS e la temperatura superficialedel corpo nel punto x al tempo t, e T0 e la temperatura ambiente, qui assuntaper semplicita uniforme e costante nel tempo.

Tenere conto della Eq. (3.1) nel programma di simulazione non e ba-nale, la Eq. (3.1) introduce infatti una dipendenza non-lineare dal campoincognito di temperatura. Lo stesso problema si incontra cercando di inclu-dere una dipendenza esplicita dei parametri termo-fisici dalla temperatura.Esistono schemi di risoluzione come il classico algoritmo di Newton, mal’implementazione di questi schemi comporta un’ulteriore complicazione nelcodice.

Bibliografia

[1] Alda J. (2003). Laser and gaussian beam propagation and trasformationIn Encyclopedia of Optical Engineering. A cura di Driggers R., pp. 999–1013. Marcel Dekker, Inc., New York.

[2] Barton J. J.; Nackman L. R. (1994). Scientific and Engineering C++. An

Introduction with Advanced Techniques and Examples. Addison-Wesley,Reading, MA.

[3] de Grooth B. G. (1999). A simple model for Brownian motion leadingto the Langevin equation. Am. J. Phys., 67(12).

[4] Devroye L. (2006). Nonuniform random variate generation In Simula-

tion. A cura di Henderson S. G., Nelson B. L., Handbooks in Opera-tions Research and Management Science, capitolo 4, pp. 83–121. Elsevier,Amsterdam, The Netherlands.

[5] Ern A.; Guermond J.-L. (2004). Theory and Practice of Finite Elements.Springer-Verlag, Berlin.

[6] Feynman R. P. (2001). Elettromagnetismo e Materia, volume 2 di La

fisica di Feynman. Zanichelli, Bologna.

[7] Greiner A.; Strittmatter W.; Honerkamp J. (1988). Numerical inte-gration of stochastic differential equations. J. of Stat. Phys., 51(1–2),95–108.

[8] Higham D. J. (2001). An algorithmic introduction to numerical si-mulation of stochastic differential equations. SIAM Review, 43(3),525–546.

[9] Joines J. A.; Roberts S. D. (1997). An introduction to object-orientedsimulation in C++ In Proceedings of the 1997 Winter Simulation Con-

ference. A cura di Andradottir S., Healy K. J., Withers D. H., NelsonB. L., pp. 78–85. IEEE Press.

[10] Joines J. A.; Roberts S. D. (1998). Fundamentals of object-orientedsimulation In Proceedings of the 1998 Winter Simulation Conference. A

29

30 BIBLIOGRAFIA

cura di Medeiros D. J., Watson E. F., Carson J. S., Manivannan M. S.,pp. 97–104. IEEE Press.

[11] Kloeden P. E.; Platen E. (1992). Numerical Solution of Stochastic

Differential Equations. Springer-Verlag, Berlin.

[12] Landau R. H.; Jose Pa M.; Bordeianu C. C. (2008). A Survey of

Computational Physics. Introductory Computational Science. PrincetonUniversity Press, Princeton.

[13] Langtangen H. P. (2003). Computational Partial Differential Equations.

Numerical Methods and Diffpack Programming. Texts in ComputationalScience and Engineering. Springer-Verlag, Berlin, 2 edizione.

[14] L’Ecuyer P. (1994). Uniform random number generation. Annals of

Operations Research, 53, 77–120.

[15] L’Ecuyer P. (1997). Uniform random number generators: A reviewIn Proceedings of the 1997 Winter Simulation Conference. A cura diAndradottir S., Healy K. J., Withers D. H., Nelson B. L., pp. 127–134.IEEE Press.

[16] L’Ecuyer P. (1998). Uniform random number generators In Proceedings

of the 1998 Winter Simulation Conference. A cura di Medeiros D. J.,Watson E. F., Carson J. S., Manivannan M. S., pp. 97–104. IEEE Press.

[17] L’Ecuyer P. (1999). Good parameters and implementations for combi-ned multiple recursive random number generators. Operations Research,47(1), 159–164.

[18] Paschotta R.; Aus der Au J.; Keller U. (2000). Thermal effects in high-power end-pumped lasers with elliptical-mode geometry. IEEE Journal

of selected topics in quantum electronics, 6(4).

[19] Press W. H.; Teukolsky S. A.; Vetterling W. T.; Flannery B. P. (2007).Numerical Recipies. The Art Of Scientific Computing. Cambridge Univ.Press, Cambridge, 3 edizione.

[20] Priesley M. B. (1981). Spectral Analysis and Time Series. Series onProbability and Mathematical Statistics. Academic Press, London.

[21] Quarteroni A. (2000). Modellistica Numerica per Problemi Differenziali.Springer, Milano.

[22] Smith I. M.; Griffiths D. V. (1998). Programming the Finite Element

Method. John Wiley & Sons, Chichester, England, 3 edizione.

[23] Sun L.; Edlou S. (2006). Low birefringence lens design for polarizationsensitive optical systems. Proc. SPIE, 6288.

BIBLIOGRAFIA 31

[24] Van Kampen N. G. (2007). Stochastic Processes in Physics and

Chemistry. Elsevier, Amsterdam, 3 edizione.

[25] Yevick D. (2005). A First Course in Computational Physics and Object-

Oriented Programming with C++. Cambridge Univ. Press, Cambridge.

[26] Zavattini E.; Zavattini G.; Ruoso G.; Raiteri G.; Milotti E.; Lozza V.;Karuza M.; Gastaldi U.; Di Domenico G.; Della Valle F.; Cimino R.;Carusotto S.; Cantatore G.; Bregant M. (2008). New pvlas results andlimits on magnetically induced optical rotation and ellipticity in vacuum.Physical Review D, 77.

[27] Zemansky M. W.; Dittman R. M. (1997). Heat and Thermodynamics.Mechanical Engineering Series. McGraw-Hill, New York, 7 edizione.

[28] Zienkiewicz O. C.; Taylor R. L. (2000). The Finite Element Method,volume 1: the basis. Butterworth-Heinemann, Oxford, UK, 5 edizione.