Analisi alle differenze finite di un cilindro con generazione di calore

Categoria: Laboratorio e progetti. Scritto da Valerio Rossi Martedì, 26 Ottobre 2010

Il seguente articolo riguarda lo studio di un cilindro con generazione di calore interno attraverso un approccio analitico ed uno numerico.

L'equazione che regola il fenomeno fisico è l'equazione generale della conduzione di Fourier:

$\nabla^2T+q_g/k=1/\alpha*{ \partial T}/{\partial t}$

che viene poi semplificata attraverso le seguenti ipotesi :

• condizioni stazionarie;
• problema monodimensionale;
• temperatura esterna costante;
• proprietà termofisiche del problema costanti.

L'equazione espressa in coordinate cilindriche e semplificata con le precedenti ipotesi diventa:

${\partial}/{\partial r}(r* {\partial T}/{ \partial r})= -q_g/k *r$

Si considera una generazione interna di calore $q_g=4*10^6 W/m^3$, un coefficiente di scambio convettivo sulla superficie esterna $h=500 W/{m^2*K}$, una conduttività termica$k=40 W/{m*K}$, una temperatura esterna $T_\infty=20 °C$ ed un raggio del cilindro $R= 10 cm$. Risolvendo l'equazione differenziale con le condizioni al contorno sopra esposte si ottiene:

$T(r)=T_\infty +{q_g}/{4*k}*(R^2-r^2)+{q_g}/{2*h}*R$

Per la soluzione numerica invece si procede utilizzando il metodo numerico delle differenze finitee scomponendo l'equazione differenziale con le differenze centrali si ottiene:

$(1+1/{2*i})*T_{i+1} - 2*T_i + (1-1/{2*i})*T_{i-1}=-q_g/k*\Delta r^2$

E' stato quindi ideato un programma in codice Fortran in grado di calcolare la soluzione analitica e la soluzione numerica considerando per le condizioni al contorno un'approssimazione sia del primo che del secondo ordine, permettendo quindi un confronto sulla precisione dei due metodi. Il programma quindi dopo aver creato una griglia di calcolo, suddividendo il dominio fisico in tanti nodi, crea una sistema di equazioni che sarà poi risolto da una subroutine contenente un algoritmo di Gaussper la risoluzione di sistemi lineari. Sotto sono visibili i risultati per una griglia di 5 nodi (a sinistra) e per una griglia di 30 nodi (a destra).

Alla fine vengono anche calcolati gli errori assoluti e relatividella soluzione numerica:

$E_{ass}=T_{ana}-T_{n um}$

$E_{rel}={E_{ass}}/{T_{ana}}*100$

dove $T_{ana}$ rappresenta la temperatura analitica e $T_{n um}$ la temperatura numerica. Di seguito gli errori relativi e assoluti per una griglia di calcolo di 30 nodi.

Il post-processing dei risultati è stato eseguito tramite il software Matlab permettendo tramite dei plot grafici il confronto dei vari risultati e mostrando una indipendenza dalla griglia di calcolo per valori di circa 30 nodi; ciò segnifica che oltre i 30 nodi l'accuratezzadella soluzione numerica risulta indipendente dal numero di nodi.

Qui per vedere la relazione completa del progetto.