Cella di Yee e il metodo degli elementi al contorno

Cella di Yee

Se consideriamo, ora, un problema di natura tridimensionale, la configurazione della griglia spazio-tempo che abbiamo utilizzato nel metodo FDTD risulta non più sufficiente per descrivere e caratterizzare il problema. Si ricorre, allora, al caso a celle tridimensionali elementari.

La cella elementare, che rappresenta l’elemento di volume infinitesimo, è quella di Yee. In questa configurazione il campo elettrico è posto al centro di ciascuno spigolo mentre la componente di campo magnetico al centro di ogni faccia (Figura 3.4).
Seguendo le notazioni di Yee, ogni punto del dominio spaziale, che è stato così discretizzato, viene identificato da un tripla ordinata

i,j,k=(i∆x,j∆y,k∆z)
in cui ∆x, ∆y e ∆z rappresentano gli incrementi spaziali nelle rispettive direzioni coordinate. Tutte le funzioni che variano con lo spazio e col tempo vengono, inoltre, rappresentate nella forma:

Fni,j,k=F(i∆x,j∆y,k∆z)
in cui ∆t altro non è se non l’incremento temporale.

Vista la
Rh(x)=fx+h-f(x)h

che è la definizione di rapporto incrementale e grazie alla disposizione geometrica delle componenti di campo nella cella di Yee, si può utilizzare la formula delle differenze centrate per ottenere una precisione del secondo ordine.

A questo punto, discretizziamo le equazioni di Maxwell sostituendo le espressioni relative alla derivata spaziale e a quella temporale.
∂Fni,j,k∂t=Fni+12,j,k-Fn(i-12,j,k)∆x+(∆x2)

∂Fni,j,k∂t=Fn+12i,j,k-Fn-12(i,j,k)∆t+(∆t2)

Poiché il campo in ciascuna cella ed in un certo istante di tempo è funzione della configurazione relativa alla cella adiacente, si può risolvere il problema per approssimazioni successive (di modo che, anche se solo in senso limite, la soluzione converga verso un valore esatto).

Di seguito è riportato un esempio di equazione alle differenze finite che evidenzia come il metodo FDTD sia efficace e versatile anche nel descrivere situazioni e geometrie complesse [30]:

Eyt+∆t,x,y+∆y2,z=Eyt,x,y+∆y2,z-∆tϵ∙1∆xHzt+∆t2,x+∆x2,y,z-Hzt+∆t2,x-∆x2,y,z-Hxt+∆t2,x,y,z+∆z2-Hzt+∆t2,x,y,z-∆z2

Tuttavia, trattandosi di un metodo che prevede un’implementazione al calcolatore, questo presenta alcune problematiche dovute alla precisione con cui vengono seguiti i contorni geometrici della struttura, al troncamento ed all’assestamento delle soluzioni stabili e stazionarie. Il troncamento rappresenta, o introduce, una discontinuità nella propagazione, che non corrisponde ad un caso reale; sebbene questa sia una considerazione legittima, però, una discretizzazione risulta necessaria perché altrimenti il problema è riferito, e va risolto, nell’ambito di un dominio infinito.

Per convenzione imponiamo due condizioni, una temporale ed una spaziale, al fine di rendere i risultati delle simulazioni il più simili e congruenti possibile con quelli analitici. La condizione di stabilità numerica prescrive che

∆t≤1c1∆x2+1∆y2+1∆z2

ovvero che ogni componente di Fourier della radiazione elettromagnetica presente all’interno del dominio di calcolo deve mantenersi stabile per tutta la durata della simulazione.

All’interno del reticolo FDTD, inoltre, la velocità di fase dei modi d’onda numerici può differire dalla velocità della luce nel vuoto con la direzione di propagazione all’interno del reticolo e con il passo di discretizzazione poiché, di fatto, essa varia con la lunghezza d’onda. Proprio a questo scopo viene imposta la condizione di non dispersione numerica:

min∆x,∆y,∆z≤λ10
che rappresenta il secondo vincolo.

È interessante notare che, poiché viene inglobata direttamente anche la sorgente nella cella elementare, i problemi di mutuo accoppiamento vengono trattati direttamente da questa procedura.

Da notare è, infine, che, proprio perché si tratta di un metodo definito nel tempo, permette di trattare segnali variabili nello stesso dominio e, quindi, risulta piuttosto utile nell’analisi di problematiche complesse, come l’esposizione del corpo umano alla radiazione elettromagnetica [da 32 a 35].

Metodo degli elementi al contorno

Come è noto, quando all’interno di un conduttore è presente un campo elettrico, esso viene attraversato da una certa densità di corrente di campo. Perché si possano studiare le prestazioni di un dispositivo, è necessario determinare la densità di corrente in esame nonché i campi termici, la cui distribuzione è dipendente non solo dalla geometria del problema o dalle condizioni al contorno ma anche da relazioni delle conduttività da variabili di campo.

La presenza di una densità di corrente è strettamente legata all’energia termica esistente, dissipata in un conduttore generico. Risulta, quindi, evidente che modi elettrici e termici devono essere trattati contemporaneamente [36].

L’approccio più comunemente adottato è basato sull’applicazione del metodo degli elementi finiti, che ha dimostrato la sua validità per un gran numero di problemi [da 37 a 39]. Tale metodo si presta molto bene a descrivere corpi anche piuttosto grandi e dalla geometria complessa che, pur tuttavia, possono essere suddivisi in elementini infinitesimi. Nel continuum, ogni singolo elemento finito viene considerato un campo di integrazione numerica di caratteristiche omogenee. Sebbene le applicazioni prevalenti siano afferenti al campo dell’ingegneria meccanica (analisi dei modelli di motori d’aereo o autovetture), negli ultimi anni l’analisi di dispositivi e campi elettromagnetici attraverso il metodo degli elementi al contorno (Boundary Element Method ovvero BEM) si è affermata [40]. Questa tecnica prevede la trasformazione delle equazioni differenziali alle derivate parziali iniziali, le quali descrivono il campo all’interno del contorno e al contorno stesso del dominio,nelle equazioni integrali relative solo ai valori al contorno delle variabili sconosciute.

Per quanto riguarda i problemi non lineari, invece, applicare questo metodo risulta ostico per via del fatto che le funzioni di Green, che descrivono i campi nei materiali non lineari, sono impossibili da trovare [41]. Per questo motivo, in questi casi, l’applicazione di questo metodo non è possibile. Data la necessità di studiare comunque certo tipo di materiali, è stato necessario lo sviluppo di un approccio parziale al problema, grazie al quale le difficoltà di cui sopra sono state superate ed è stata possibile un’efficace applicazione di questo metodo [42].

Lo studio in oggetto concorda con l’analisi di tipo “steady mode” (modo stazionario) di accoppiare la densità di corrente ed i campi termici con l’andamento non lineare della conduttività (sia termica sia dielettrica) mediante il BEM. Perché questo studio possa essere condotto, risulta necessario identificare le cosiddette funzioni di Green.

Si tratta di uno strumento matematico atto alla risoluzione delle equazioni differenziali non omogenee. A partire da un operatore differenziale lineare L, infatti, possiamo agire su un opportuno spazio (topologico o vettoriale) di funzioni, e scrivere Lxu(x) = f(x) che altro non è che un’equazione differenziale ordinaria. Definiamo, quindi, la funzione di Green come: LxG(x,y) = δ(x-y) in cui δ è la funzione impulsiva di Dirac. Restando ferma, tuttavia, la impossibilità a caratterizzare il problema mediante questo strumento, e per eliminare la necessità delle funzioni di Green per ottenere le equazioni differenziali alle derivate parziali non lineari, verrà applicata la trasformazione di Kirchhoff alle variabili di campo.

Dagli studi del fisico tedesco è emersa la possibilità di stabilire che la trasmissione del calore che avviene per irraggiamento è il risultato di un bilancio fra energia emessa ed assorbita. È, allora, possibile affermare che il rapporto tra emittanza monocromatica (e) e coefficiente di assorbimento (a) di uno stesso corpo è una funzione universale, in quanto indipendente dal corpo stesso, ed è funzione della sola lunghezza d’onda e della temperatura: e(λ,T)a(λ,T)=f(λ,T). Se in questa relazione poniamo a=1, essa viene riferita al corpo nero, ovvero ε0 = f(T,λ). Nello spazio libero, invece, la funzione di Green può essere utilizzata senza particolari restrizioni.

Formulazione delle equazioni integrali
Indichiamo ora con Ω un conduttore che abbia conducibilità termica e conducibilità termica λ=λ(T) in cui indichiamo con T la temperatura assoluta. All’interno del conduttore in esame ci sono due campi, il vettore di campo relativo alle densità di corrente ed un campo termico scalare T. La distribuzione di corrente che si ha all’interno del conduttore può essere ricavata mediante la seguente relazione:

∙J=
-σT∇V=0
in cui V è, chiaramente, il potenziale ed indichiamo con Γ il dominio Ω. Per risolvere l’equazione in oggetto è necessario dare una descrizione completa del problema attraverso le condizioni al contorno. Possiamo, allora, scrivere:
V=V su Γ1
Jn=-σT∂V∂n=0 su Γ2
avendo scomposto il contorno in due porzioni, resta comunque verificata la relazione. Per risolvere l’equazione possiamo applicare la funzione di Green nello spazio libero all’interno dell’equazione integrale risultante, in modo che diventi:

2V=-1σ
σ∙
V=-1σ∂σ∂T∇T∙
V
Visto che in moltissimi casi il valore della conducibilità è ottenuto analiticamente, il coefficiente ∂σ∂T può essere valutato attraverso il calcolo differenziale. Il vettore densità di corrente può essere calcolato come:
J=-σT
V
La distribuzione del campo termico all’interno del conduttore (Ω) è descritto dalla relazione:

-λT∇T=γT(∇V)2
dove (∇V)2=∇V∙∇V. All’equazione che descrive la distribuzione di campo vanno associate le relative condizioni al contorno:
T=T su Γ3
qn=-λT∂T∂n=qn(T) su Γ4
in cui T, q, λ(T) e n sono, rispettivamente, temperatura, flusso di calore, conducibilità termica (che dipende sempre dalla temperatura) ed il vettore normale uscente dalla superficie. Risulta evidente che, anche per Γ3∪Γ4=Γ. Il normale flusso di calore qn sulla porzione Γ4 del contorno può essere una funzione non lineare della temperatura T. Il flusso termico all’interno del conduttore viene, quindi, espresso come:
q=-λ(T)∇T
Per alcune delle funzioni ottenute, λ(T) ad esempio, risulta molto difficile, se non impossibile, ricavare la funzione di Green come una soluzione dell’equazione

-λT
G*r,r0=-δr,r0
in cui r è un punto di osservazione ed r0 è un punto di sorgente (integrazione). Allora si sfrutta la trasformazione di Kirchhoff per fare in modo da portare all’equazione

-λT∇T=γT(∇V)2 in una forma più semplice e che consenta l’uso della funzione di Green nello spazio libero [43]

∇U=λ(T)∇T
dU=λ(T)dT
in cui U è un nuovo problema variabile legato alla temperatura. Infatti, integrando entrambi i membri tra T0 e T si ottiene la relazione
U=KT=T0Tλ(T)dT
e T0 risulta essere il valore di riferimento arbitrario. Possiamo, quindi, scrivere:

2U=γT(∇V)2
Quella che abbiamo ottenuto in questo modo è un’equazione di Poisson che è stata scritta, però, in termini della nuova variabile U. Quando il valore di U dovesse essere noto, quello di T potrà essere calcolato attraverso la trasformazione di Kirchhoff inversa
T=K-1(U)
Le condizioni al contorno per la nuova variabile U sono ottenute attraverso le relazioni:
U=U=T0Tλ(T)dT
Qn=∂U∂n=∂U∂T∂T∂n=λT∂T∂n=-qn
Se poniamo, ora,

2V=bV

2U=bU
le equazioni che descrivono i campi elettrici ed i campi termici possono essere scritti nella forma
bV=-1σ(U)∂σ(U)∂n(∇V)∙(∇V)
bU=-σ(U)(∇V)2
in conclusione, le equazioni integrali al contorno associate con le relazioni scritte sopra vengono ricavate da [44]:
c1V+Γ2∂G∂nVdΓ-Γ1G∂V∂ndΓ=Γ2V∂G∂ndΓ-ΩbVGVdΩ
c2U+Γ4∂G∂nUdΓ-Γ3G∂U∂ndΓ=-Γ3U∂G∂ndΓ-Γ3qnGdΓ+ΩbVGVdΩ
in cui G(r,r0) e la funzione di Green nello spazio libero e le costanti c1 e c2 vengono determinate dalle variabili di integrazioni principali di Cauchy e sono relative alle singolarità della funzione di Green. È importante notare che esistono due domini integrali (si tratta di Ω e Γ) per ciascuna equazione sul lato destro; questo rende le equazioni integrali sopracitate non totalmente dipendenti dal contorno. Se da un lato, quindi, siamo riusciti ad utilizzare la funzione di Green nello spazio libero, questa dipendenza è d’altronde imprescindibile.

Calcolo degli integrali di dominio
Le equazioni che costituiscono il risultato appena ottenuto, possono essere risolte con un metodo iterativo. Si ottiene, infatti, che le costanti c1 e c2 sono uguali, per i punti r ϵ Ω ed entrambe unitarie. Nonostante ciò, in questo caso, per entrambi i termini delle relazioni precedenti, l’operatore divergenza può essere applicato e le relazioni risultanti sono:
∇V=-Γ2
∂G∂nVdΓΓ1∇G∂G∂n-∇V∂G∂ndΓ+ΩbVGVdΩ
∇U=-Γ4
∂G∂nU+qn∇GdΓΓ3∇G∂G∂n-∇V∂G∂ndΓ+ΩbVGVdΩ
e risulta evidente che l’operatore divergenza dipenda dalla variabile r, mentre V e ∂V∂n dipendono da r0. Così, sotto il segno di integrale, un operatore divergenza è applicato alla funzione G(r,r0) e ∂G(r,r0)∂n. Entrambe le funzioni bU e bV dipendono anche dalla variabile r0 pur non essendo differenziate.
Se analizziamo l’equazione associata a ∇V, notiamo che la funzione integrale bV dipende dal gradiente di V, come è evidente dalla presenza del termine differenziale. È facile verificare, inoltre, che non è possibile valutare direttamente il gradiente di V ed U dalle relazioni precedenti, ed è per questo che si è resa necessaria la sostituzione dei valori bV e bU, ricavati dal passaggio precedente, negli integrali di dominio, sostituzione dettata dal metodo iterativo che siamo chiamati ad adottare.
Per un problema bidimensionale, la funzione di Green viene definita da:
G(r,r0)=-12πlnr0-r
I gradienti di queste espressioni sono dati dalle seguenti equazioni:

rGr,r0= R(r0)2πR2

r∂Gr,r0∂n(r)= 12πR2- n(r0)+2R(r0)∙n(r0)R2R(r0)
in cui sono verificate le relazioni
R(r0)=r0-r=x0-xi+y0-yj
n(r0)=cosnr0,ii+cosnr0,ij
perché sia possibile la valutazione degli integrali di dominio, è necessario dividere l’intero Ω in più elementi di superficie
ΩbVr0
rG(r,r0)dΩ(r0)=i=1NbVri,0
rG(r,ri,0)∆S(ri,0)
in cui la sommatoria comprende gli elementi di superficie sopracitati e ∆S(ri,0), che altro non è che un valore degli elementi per un punto medio ri,0. Per semplificare i conti e per moltiplicare la discretizzazione del dominio, assumiamo che gli elementi al contorno siano costanti.
Le equazioni possono essere discretizzate mediante la tecnica BEM; il risultato che otterremo saranno due equazioni algebriche mutuamente accoppiate
Visto che all’interno del dominio i termini bV e bU dimostrano un comportamento non lineare, ancora una volta il sistema di equazioni dovrà essere risolto in maniera iterativa. Si procede, quindi, assumendo, per prima cosa, lineare il valore della conduttanza σ. In questo modo spariscono i termini di dominio e l’equazione risulta risolvibile per il potenziale V. Esso risulterà, quindi, valutato per primo. Successivamente, e dalle equazioni già viste, è possibile calcolare, bU, qn; fatto ciò, si può risolvere il campo scalare U ed infine si calcolano i valori di
V e ∇U insieme con ∂σ∂U. Il valore di bV è calcolabile ed è possibile fare una seconda stima della distribuzione di potenziale di V in Ω. Il processo iterativo continuerà fino a quando non sarà ottenuta la convergenza dei risultati. Per ottenere, infine, la temperatura di campo principale T è necessario applicare la trasformazione inversa di Kirchhoff T=K-1(U).

In allegato i risultati numerici.

L’indice completo degli articoli relativi alla tesi di laurea sulla interazioni e sugli effetti delle radiazioni sul corpo umano, è disponibile qui

Leave a Reply