Soluzione PDE con metodo Crank-Nicolson

per L’argomento equazioni alle derivate parziali è stato trattato qualche mese fa dall’autore su Matematicamente.it (R1). Ad esso si rimanda per una introduzione sulla PDE (Partial Derivative Equation) e per la descrizione/soluzione di un esempio applicativo di interesse tecnico: lo studio del transitorio termico per la valvola di alimentazione di un motore a scoppio. Questo articolo si propone di riesaminare e risolvere lo stesso problema con un metodo numerico diverso e più efficiente.

Qui sotto è riportato il disegno del motore con la valvola di alimentazione (in verde).

Sezione motore a scoppio con valvola di alimentazione

 

 

 

 

 

 

 

 

 

 

Mentre la valvola è disegnata sotto:

Valvola motore

 

 

 

 

 

La valvola di alimentazione riceve alla sua estremità sinistra (sx) un flusso termico dalla camera di combustione, che si è supposto abbia una temperatura media $ T_m = 1300 °C $, viene lambita dal flusso di miscela in alimentazione (freccia in blu nel disegno del motore), che ha temperatura $ T_a = 20 °C $, la quale asporta calore lungo tutta la superficie laterale della valvola. Infine l’estremità destra (dx) della valvola è a contatto con la camma che pesca nell’olio del carter (che assumiamo avere temperatura costante pari a $ T_c = 80 °C $) e quindi avviene uno scambio termico (in entrata o uscita) tra la valvola e l’olio.

La valvola ha la forma di un cilindro sottile e viene schematizzata come un’asta, vale a dire un corpo monodimensionale.

L’equazione che descrive il fenomeno in forma indefinita:

\[ \begin{equation}\frac{\partial T}{\partial t}=D \frac{\partial^2 T}{\partial x^2}-G(T-T_a) \label{eq1}\tag{1}\end{equation}\]

è di tipo lineare parabolico e correla l’andamento della temperatura (T) lungo l’asse della valvola in funzione del tempo (t). Il termine $ G(T-Ta) $nel membro a destra rappresenta la perdita di calore della valvola mentre viene lambita dal flusso della miscela di alimentazione al motore. La derivazione della \(\eqref{eq1}\) è descritta in dettaglio nell’articolo citato.

Risolvere l’equazione vuol dire determinare (nel nostro caso per via numerica) la funzione T (x, t) sulla base della condizione iniziale (t=0) ed al contorno alle estremità (x=0 e x= L).

Nel precedente articolo questo esempio è stato risolto approssimando le derivate parziali con uno schema molto semplice di differenze finite che si presta ad una soluzione diretta, con il metodo di Eulero esplicito. In pratica si è costruito un reticolo spazio-temporale

Schema differenze finite - Metodo di Eulero esplicito

 

 

 

 

 

 

 

 

 

 

 

dove lo spazio è riportato in ascissa ed il tempo in ordinata. Con il metodo esplicito di Eulero si risolve ogni cella del reticolo spostandosi da sinistra a destra fino a completare le riga per passare poi alla riga temporale successiva. Il metodo, descritto in dettaglio nell’articolo citato, è molto semplice e si presta a essere facilmente usato con EXCEL, ma ha un serio problema di convergenza. Questo richiede di usare un passo temporale molto piccolo. In tal modo per ottenere la convergenza, ed anche un ragionevole precisione, è necessario costruire un reticolo di grandi dimensioni. Ciò che comporta un notevole consumo di memoria del computer ed un tempo di elaborazione abbastanza lungo.

METODO DI CRANK-NICOLSON

Per rimediare alle difficoltà di convergenza gli esperti hanno sviluppato nuovi schemi alle differenze finite per approssimare le derivate parziali evitando i problemi del metodo di Eulero. In letteratura (R2) si trovano diversi metodi, uno dei più usati è quello di due matematici statunitensi John Crank e Phyllis Nicolson (in seguito C-N) (R3, R2) che vogliamo applicare in questo articolo. Vediamo come è fatto e come si usa.

Con il metodo C-N le derivate sono approssimate come segue:

\[ \begin{equation} \frac{\partial T}{\partial t} \approx \frac{T^{j+1}_i – T^j_i}{\Delta t} \label{eq2} \tag{2} \end{equation} \]

\[ \begin{equation} \frac{\partial^2 T}{\partial x^2} \approx \frac{1}{2(\Delta x)^2} \Big( T^{j+1}_{i-1} -2 T^{j+1}_{i} + T^{j+1}_{i+1} + T^{j}_{i-1} -2 T^j_i + T^j_{i+1}\Big) \label{eq3} \tag{3} \end{equation} \]

(dove: pedice = spazio, apice = tempo)

Notiamo che la derivata temporale è molto semplice (come nel metodo di Eulero esplicito) mentre la derivata spaziale è di tipo più complesso, in quanto coinvolge sei elementi del reticolo spazio-temporale, suddivisi in due gruppi di tre elementi per ogni valore temporale, che è rappresentato dal seguente schema (stencil):

Stencil di Crank-Nicolson

 

 

 

 

 

 

 

L’approssimazione \(\eqref{eq3}\) è dunque un po’ complicata ma comporta un errore di troncamento quadratico nello spazio e nel tempo quindi una maggior precisione rispetto allo schema di Eulero esplicito. Il metodo è convergente ed incondizionatamente stabile. Ha lo svantaggio di essere implicito, il che vuol dire che la soluzione dello stencil richiede la soluzione di un sistema di equazioni lineari, come vedremo.

Per completare lo schema risolutivo è necessario approssimare alle differenze finite anche il secondo termine al secondo membro della \(\eqref{eq1}\):

\[ \begin{equation} G(T-T_a) \simeq G(T^{j+1}_i – T_a) \label{eq4} \tag{4} \end{equation} \]

Infine dobbiamo fornire la condizione iniziale:

\[ t = 0 \,\,\,\,\,\,\,\,, \,\,\,\,\,\,\,\, T (x, t= 0) = T_0 \]

Vale a dire:

\[ \begin{equation}T^1_i = T_0\,\,\,\,\,\,\,\,\, \mbox{ per } i=1, n \label{eq5}\tag{5} \end{equation} \]

E le condizioni al contorno, che nel nostro caso, come nell’articolo precedente, sono del tipo di Neumann (R4), vale a dire sono espresse da derivate:

Estremità sx ( lato motore) : \( \gamma_m (T_m – T) = -k \frac{\partial T}{\partial x} \)

Estremità dx (lato carter): \( \gamma_c (T_c – T) = k \frac{\partial T}{\partial x} \)

Nota: sx e dx significano rispettivamente sinistra e destra.

(dove $ k $ è la conduttività termica della valvola, mentre \( \gamma_m \) e \(\gamma_c\) e sono i coefficienti di scambio termico convettivo lato motore e lato carter)

Entrambe le condizioni al contorno traducono il principio fisico per cui il flusso termico convettivo, tra l’estremità della valvola ed il fluido che la lambisce, eguaglia il flusso termico conduttivo nel metallo alle estremità della valvola: in definitiva il principio di conservazione dell’energia.

Anche per le condizioni al contorno le derivate sono approssimate dalle loro differenze finite:

\[ \begin{equation} \gamma_m (T_m – T_{1,j+1}) = \frac{k}{\Delta x} (T^{j+1}_1 – T^{j+1}_2) \label{eq6}\tag{6} \end{equation} \]

\[ \begin{equation} \gamma_c (T_c – T_{n,j+1}) = \frac{k}{\Delta x}(T^{j+1}_{n-1}-T^{j+1}_n) \label{eq7}\tag{7} \end{equation} \]

Inoltre è conveniente concentrare i parametri creandone cinque nuovi, in questo modo:

\( \mu = G\Delta t \)

\( \Theta_c = \gamma_c \frac{\Delta x}{k} \)

\( \Theta_m = \gamma_m \frac{\Delta x}{k} \)

\( \lambda = \frac{D}{2} \frac{\Delta t}{(\Delta x)^2} \)

\( \alpha = 1 + 2\lambda + \mu \)

Non è irrilevante notare che tali parametri sono tutti adimensionali.

Ora le \(\eqref{eq2}\) e \(\eqref{eq3}\) possono essere sostituite nella \(\eqref{eq1}\) facendo uso dei parametri concentrati:

\[ \begin{equation} -\lambda T^{j+1}_{i-1} + \alpha T^{j+1}_{i} – \lambda T^{j+1}_{i+1} = \lambda T^{j}_{i-1} + (1 -2\lambda) T^{j}_{i} + \lambda T^{j}_{i+1} + \mu T_a \label{eq8}\tag{8}\end{equation} \]

Osserviamo che la \(\eqref{eq8}\) è stata scritta portando a sx le temperature relative allo step temporale $ j+1 $, a dx quelle per lo step $ j $.

Analogamente si possono riscrivere e semplificare le \(\eqref{eq6}\) e \(\eqref{eq7}\) in forma sintetica:

\[ \begin{equation} (1 + \Theta_m ) T^{j+1}_1 – T^{j+1}_2 = \Theta_m T_m \label{eq9}\tag{9} \end{equation} \]

\[ \begin{equation} -T^{j+1}_{n-1} + (1+\Theta_c) T^{j+1}_n = \Theta_c T_c \label{eq10}\tag{10} \end{equation} \]

In definitiva le equazioni \(\eqref{eq8}\), \(\eqref{eq9}\), \(\eqref{eq10}\) con la condizione iniziale \(\eqref{eq5}\) risolvono il nostro problema. Vediamo come.

Notiamo che la \(\eqref{eq8}\) è in realtà un sistema di $ n – 2 $ equazioni lineari (con $i = 2, n -1$), mentre le \(\eqref{eq9}\) e \(\eqref{eq10}\) sono equazioni singole. Quindi il sistema \(\eqref{eq8}\), \(\eqref{eq9}\), \(\eqref{eq10}\) è costituito da equazioni in altrettante variabili \(T^{j+1}_i \) dove $ i=1, n$ per ogni valore $j+1$. In altri termini abbiamo creato un sistema lineare a matrice quadrata \( n \times n\). Se risolviamo il sistema a partire dalla condizione iniziale \(\eqref{eq5}\) la prima riga temporale è nota: tutti i termini valgono $T_0$. Vediamo quindi che i membri di $dx$ del sistema di equazioni sono noti e costituiscono il termine noto del sistema di equazioni. Non c’è che da risolvere il sistema lineare \( n \times n\) per calcolare il profilo di temperatura della valvola lungo l’asse x per la riga temporale $j+1$.

Tuttavia a questo punto vien da pensare che la faccenda diventi complessa. Infatti è il numero di punti nel reticolo spaziale, che si ottiene dividendo la lunghezza della valvola in $n-1$ intervalli uguali di lunghezza \(\Delta x\). Ora per avere una ragionevole precisione di calcolo è facile immaginare che l’intervallo \(\Delta x\) debba essere abbastanza piccolo. In tal modo è facile rendersi conto che la \(\eqref{eq8}\) sarà costituita da decine e forse centinaia di equazioni lineari. Ma non c’è motivo di preoccuparsi. Esistono algoritmi di soluzione semplici e molto precisi per la soluzione iterativa di tali sistemi( Jacobi e meglio ancora Gauss-Seidel). Nel nostro caso, poi, la situazione è ancora più semplice perché la matrice della \(\eqref{eq8}\) è tridiagonale ed esiste un metodo diretto (vale a dire non iterativo) per la soluzione esatta. Vediamo di visualizzare come mai il sistema risolutivo è tridiagonale. Intanto notiamo che ogni equazione del gruppo \(\eqref{eq8}\) contiene solo tre variabili, quella sulla diagonale centrale e quelle a lati sx e dx. Di seguito vediamo come è fatta la matrice dei coefficienti del sistema:

Matrice dei coefficienti del sistema

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(Le temperature riportate in testa al disegno sono prive dell’apice temporale j+1, che viene sottinteso).

Notiamo che sulla diagonale principale abbiamo tutti coefficienti uguali, pari a α, tranne il primo e l’ultimo, che sono rispettivamente i coefficienti delle equazioni al contorno \(\eqref{eq8}\) e \(\eqref{eq9}\). Le diagonali immediatamente a dx ed a sx sono piene di valori uguali, pari a -λ, tranne gli estremi.

In conclusione la matrice dei coefficienti del sistema è tridiagonale ed i valori delle componenti matriciali sono delle costanti, che rimangono invariate passando da una riga temporale all’altra.

La matrice dei termini noti ovviamente deve essere variabile per ogni riga temporale. È un vettore colonna il cui elemento tipico è dato in forma generale da:

\[ \lambda T^j_{i-1} + (1-2\lambda) T^j_i + \lambda T^j_{i+1} + \mu T_a, \,\,\,\,\,\,\,\, i = 2, n-1 \]

Mentre i componenti $1$ ed $n$ sono rispettivamente:

\[ \Theta_m T_m \,\,\,\,\,\,\,\, \mbox{ e } \,\,\,\,\,\,\,\, \Theta_c T_c \]

In breve il meccanismo di calcolo è facilmente intuibile: si parte dalla prima riga temporale, dove

\( T^1_i = T_0 \,\,\,\,\,\,\,\, \mbox{ per } i = 1, n \)

Ora si calcolano le componenti del vettore dei termini noti usando le \(T^{1}_i\). A questo punto si risolve il sistema \(n \times n\) ottenendo le temperature per la seconda riga temporale $j = 2$. Si ripete il procedimento per tutte le righe temporali successive. Notiamo che mentre il reticolo spaziale è definito dalla lunghezza della valvola, il reticolo temporale è teoricamente illimitato. Tuttavia è evidente che dopo un tempo abbastanza lungo si raggiunge una situazione di equilibrio termico: il profilo di temperatura lungo l’asse x non cambia, quindi la valvola ha raggiunto lo stato stazionario e non c’è motivo di aggiungere altre righe.

Ci manca di precisare come si risolve il sistema tridiagonale: ci ha pensato il fisico e matematico Llewelynn Thomas che ha costruito un semplice algoritmo (R5). Perché esso funzioni è necessario che la diagonale della matrice dei coefficienti sia strettamente dominante per righe, ciò che si verifica nel nostro problema.

Le equazioni del modello sono state incorporate in un semplice codice VBA che contiene la routine dell’algoritmo di Thomas. Il codice è stato usato per simulare il comportamento della valvola a partire dalla condizione iniziale per diverse dimensioni del reticolo. Converge sempre perfettamente.

La simulazione più precisa ottenuta con il modello ha richiesto un reticolo \(i \times j = 150 \times 256 =38400\) celle ed il profilo di temperatura a tempi crescenti è riportato sotto

 

 

 

 

 

 

 

 

 

 

 

Riproponiamo sotto il diagramma da articolo precedente (R1), risolto con metodo di Eulero esplicito:

 

 

 

 

 

 

 

 

 

 

 

È utile confrontare i profili di temperatura all’equilibrio termico, che in pratica è raggiunto quando t=1024s, ottenuti con Eulero e C-N rispetto alla soluzione esatta, di tipo analitico (quest’ultima è descritta in APPENDICE). Questo si vede nel diagramma successivo.

 

 

 

 

 

 

 

 

 

 

 

Si nota che i profili di Eulero e C-N sono abbastanza vicini, tuttavia all’estremità sx C-N è molto più preciso di Eulero, mentre all’estremità dx i due metodi sono quasi coincidenti e sovrastimano di circa 9 °C la soluzione esatta. Complessivamente C-N è più preciso di Eulero, ha richiesto un reticolo enormemente più piccolo (38 mila contro 768 mila punti) e tempo di elaborazione molto inferiore. Si può concludere che è sicuramente preferibile.

Il codice VBA e tutti i dati, parametri, calcoli, diagrammi, sono riportati nel file EXCEL associato all’articolo.

APPENDICE-SOLUZIONE ANALITICA ALL’EQUILIBRIO TERMICO

All’equilibrio termico, che in teoria si raggiungerebbe dopo un tempo infinito, l’equazione \(\eqref{eq1}\) si semplifica dal momento la derivata temporale è nulla: \( \frac{\partial T}{\partial t} = 0\).

Questo vuol dire che scompare la variabile tempo. Il modello si semplifica e la derivata spaziale non è più parziale:

\[ \begin{equation}\frac{d^2 T}{dx^2} = \frac{G}{D} (T-T_a) \label{eq11} \tag{11}\end{equation} \]

Che è del tipo: \(\frac{d^2 T}{dx^2} = a^2 (T-b)\) dove \(a^2=\frac{G}{D}, b = T_a \).

È questa una equazione differenziale lineare del secondo ordine la cui soluzione è:

\[ \begin{equation} T=k_1 \exp(ax)+k_2 \exp(-ax)+b \label{eq12}\tag{12} \end{equation} \]

Si impongono ora le condizioni al contorno di Neumann:

Estremità sx, x=0: \[ {\frac{dT}{dx}}_{x=0} = -\frac{\gamma_m}{k}(T_m-T(0)) \]

Estremità dx, X =L: \[ {\frac{dT}{dx}}_{x=L} = -\frac{\gamma_c}{k}(T_c-T(L)) \]

In questo modo si ottiene un sistema lineare \(2 \times 2\) nelle variabili $k_1$ e $k_2$ che si risolve ad es. con la regola di Kramer. La soluzione esatta nel caso della nostra valvola determina le seguenti temperature alle estremità: \( T(0) = 218°C, T(L) = 27°C\).

Profilo di temperatura all'equilibrio termico - Soluzione analitica

 

 

 

 

 

 

 

 

 

 

Il profilo di temperatura, riportato nel diagramma sopra, sembra di tendenza asintotica, ma si tratta di un inganno visivo. Nella realtà esiste un minimo relativo, situato a circa 0,1 mm dall’estremità dx della valvola. Il significato fisico del fenomeno è semplice da capire: il flusso in controcorrrente della miscela in alimentazione raffredda la valvola a temperature inferiori a 30°C. Tuttavia la valvola all’estremità dx è in contatto con l’olio caldo del carter ($T_c = 80° C$) che provoca un leggero aumento di temperatura in tale zona.

RIFERIMENTI

 

Commenti

commenti