Dinamica di un circuito mutuamente accoppiato

Modello analitico in python

Questo articolo è ridondante dal momento che si è parlato di come risolvere circuiti elettrici tramite il calcolo analitico dell variabili di stato, utilizzando però Matlab. In questa sede verrà riproposto il medesimo approccio ma utilizzando uno strumento open source, ovvero Python.

Una volta importate le apposite librerie, si crea il modello della funzione e si chiama la funzione odeint.

Si consideri un circuito con mutui accoppiamenti (quindi una prima approssimazione di un trasformatore) come quello in figura:

Le cui equazioni sono le seguenti:

v_1(t) = L_1 \frac{d i_1}{dt} + M \frac{d i_2}{dt}
v_2(t) = M \frac{d i_1}{dt} + L2 \frac{d i_2}{dt}

Si osserva subito che:

v_1(t) = e(t) - R_1 * i_1(t)
v_2(t) =- R_L * i_1(t)

Il cambio di segno nella seconda equazione è dovuta alla convenzione di segno (fisico) adottata. Mettendo a sistema le suddette equazioni, si riesce a scrivere:

\frac{d i_1}{dt} = \frac{\frac{MR_L}{L_2}}{L_1-\frac{M^2}{L_2}}i_2(t)+ \frac{e(t)}{L_1-\frac{M^2}{L_2}}- \frac{R_1}{L_1-\frac{M^2}{L_2}}i_1(t)
\frac{d i_2}{dt} = (- \frac{R_L}{L_2} - \frac{M}{L_2}\frac{\frac{MR_L}{L_2}}{L_1-\frac{M^2}{L_2}}))i_2(t) - \frac{M}{L_2(L_1-\frac{M^2}{L_2}) }e(t) + \frac{M R_1}{L_2 (L_1-\frac{M^2}{L_2})} i_1(t)

In sostanza si è arrangiato il sistema per ottenere la seguente forma:

\frac{d x}{d t} = f(x,y,t)
\frac{d y}{d t} = f(x,y,t)

Il modello viene così scritto in Python:

def trafo(Z,t,e):
r1 = R1
rl = RL
l2 = L2
l1 = L1
m = M
a = l1-(pow(m,2)/l2)
b = m*rl/l2
#Condizioni iniziali (al tempo [i])
i1 = Z[0]
i2 = Z[1]
#Sistema di equazioni differenziali
di1_dt = (b/a)i2-(r1/a)*i1+e/a
di2_dt = (- (rl/l2)-m*b/(l2*a))*i2+(m*r1/(l2*a))i1-(m/(l2*a))*e
#Vettore soluzioni del sistema
Z =[di1_dt ,di2_dt]
return Z

Le variabili a e b sono variabili di comodo (per semplificare la scrittura delle equazioni). I parametri del circuito sono stati dichiarati come variabili globali (in maiuscolo) e richiamati all’interno del modello come variabili locali (in minuscolo).

Il solver ODE viene fatto lavorare all’interno di un ciclo for: all’iterazione i-esima viene fornito il relativo valore della forzante e(t) (che quindi diventa e(i))) e le nuove condizioni iniziali che saranno i valori di i1 ed i2 all’terazione “i-1” , ovvero il risultato dell’iterazione precedente. Notare che il contatore “i” non parte da 0 ma da 1: questo perchè, all’istante 0 i valori di i1 ed i2 sono forniti dall’utente (condizioni iniziali al tempo 0).

#condizioni iniziali al tempo 0
z0 = [0,0]

la funzione forzante e il vettore “tempo” vengono generati a monte della suddetta iterazione, in funzione nel numero “n” di iterazioni:

t = np.linspace(0,t_max,n)
e = np.zeros(n)
for i in range(1,n):
e[i] = 5 + V*math.sin(Wt[i])

Quindi si procede con l’iterazione:

for i in range(1,n):
#passo temporale
tspan = [t[i-1],t[i]]
#risolvi step
z = odeint(trafo,z0,tspan,args=(e[i],))
i_1[i] = z[0][0]
i_2[i] = z[0][1]
#Aggiorna le "condizioni iniziali"
z0 = z[1]

Notare che in tspan viene definito l’intervallo di tempo nel quale si vuole risolvere il sistema: in altre parole, questo approccio, consente una risoluzione a tratti. Aumentado il numero passi “n”, mantenendo invariato il tempo totale di cacolo (t_max), si otterranno risultati migliori.

Si è dato per assodato, dato che si simula per basse frequenze, che il solver stesso, di default, presenti passo di calcolo adeguato.

Di seguito il codice completo:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
#Parametri
R1 = 1
RL = 3
L1 = 0.01
L2 = 0.1
M = 0.9* pow((L1L2),0.5)
n = 1000 #Numero di campioni
V = 10 #Valore di picco tensione in ingresso
f = 10 #Frequenza
W = 2 * 3.14 * f
t_max = 0.6
#Inizzializza vettrori correnti
i_1 = np.zeros(n)
i_2 = np.zeros(n)
#Modello del circuito con accoppiamenti mutui
def trafo(Z,t,e):
r1 = R1
rl = RL
l2 = L2
l1 = L1
m = M
a = l1-(pow(m,2)/l2)
b = mrl/l2 #Condizioni iniziali (al tempo [i])
i1 = Z[0]
i2 = Z[1]
#Sistema di equazioni differenziali
di1_dt = (b/a)i2-(r1/a)i1+e/a
di2_dt = (- (rl/l2)-mb/(l2a))i2+(mr1/(l2a))i1-(m/(l2a))*e
#Vettore soluzioni del sistema
Z =[di1_dt ,di2_dt]
return Z
#condizioni iniziali al tempo 0
z0 = [0,0]
#time points
t = np.linspace(0,t_max,n)
#input
e = np.zeros(n)
for i in range(1,n):
e[i] = 5 + Vmath.sin(Wt[i])
for i in range(1,n):
#passo temporale
tspan = [t[i-1],t[i]]
#risolvi step
z = odeint(trafo,z0,tspan,args=(e[i],))
i_1[i] = z[0][0]
i_2[i] = z[0][1]
#Aggiorna le "condizioni iniziali"
z0 = z[1]

Una nota sulla definizione dei parametri è necessaria. Dato un circuito la cui morfologia è quella qui presentata, è necessario che:

M^2 < L_1L_2

Questa condizione fisica deve essere rispettata per amor della Fisica: se così non fosse, l’accumulo energetico sarebbe negativo. Notare che non è una disuguaglianza in senso stretto, è ammessa anche l’uguaglianza, ma per una maggiore coerenza del modello matematico si è scritto:

M = 0.9* pow((L1L2),0.5)

Il coefficiente 0.9 è stato scelto di autorità, ma non deve mai superare l’unità (per le suddette ragioni).

Verifica bontà del modello

I parametri sono arbitrari, ma il comportamento del sistema “reale” è noto, pertanto si procede con simulazioni i cui risultati sono ben noti. Se il risultato della simulazione è coerente con la fisica del problema, allora il modello sarà valido.

Transitori impulsivi

Applicando un gradino di tensione ci si aspetta che:

a) La corrente iniziale sia nulla: le induttanze non ammettono gradini di corrente;

b) Il valore finale della corrente i1 (a regime) sia dato dal rapporto tra l’ampiezza del gradino e la resistenza R1;

c) Il valore di corrente al secondario sia nullo a transitorio estinto: non essendovi, a regime, variazioni di tensione sul primario, non può sussistere f.e.m. al secondario.

Quindi, dato il seguente gradino di tensione:

e = np.zeros(n)
for i in range(1,n):
if i > n/10:
e[i] =V

Il risultato è il seguente:

Risposta a gradino

la tensione v1 (in viola) assume il valore di “e” all’istante iniziale, dal momento che la corrente i1 (in verde) è nulla. Negli istanti successivi la corrente i1 inizia a crescere fino a raggiungere asintoticamente il suo valore finale (che è pari a V/R1). Notare che al secondario i parametri tensione e corrente sussistono solamente durante le variazioni dei parametri al primario: la corrente i2 (in rosso) e la tensione v2 (in arancione) si evolvono secondo le dinamiche di un transitorio del primo ordine ma si estinguono una volta che sul primario i parametri diventano costani.

NOTA: la corrente i2 assume valori negativi. Questo è coerente per la convenzione di segno adottata. In caso contrario la resistenza produrrebbe energia rendendo necessario l’intervento delle autorità ecclesiastiche.

Disaccoppiando i circuiti, ponendo M=0, ci si aspetta un transitorio RL al primario (al secondario non accade nulla).

Infatti:

Circuiti disaccoppiati (M=0)

Chiusura di un generatore (ideale) di f.e.m. sinusoidale

Alimentando il tutto da una rete prevalente ci si aspetta la presenza di componente unidirezionale della corrente.

La tensione sarà così formata:

e = np.zeros(n)
for i in range(1,n):
e[i] = 5 + V*math.sin(W*t[i])

Nota: la presenza di offset serve a testare ulteriormente il modello, dal momento che non la si deve ritrovare al secondario (per ovvie ragioni fisiche).

In questo caso ci si aspetta che:

a) la componente continua della tensione e(t) non abbia effetti sul secondario;

b) La componente unidirezionele di corrente deve estinguersi;

c) Le correnti siano in ritardo (carico induttivo):

Lanciando la simulazione, il risultato è il seguente:

A colpo d’occhio si può affermare che il risultato è coerente con la Fisica. Nello specifico:

Grandezze al secondario

La componente unidirezionale si annulla. Al secondario, a regime, non vi è traccia di offset.

Grandezze al primario

Al primario è presente un offset di tensione e quindi di corrente: si osservi che, a regime, il valor medio della tensione sull’indutture (mutuamente accoppiato) è nullo. Se così in fosse, la corrente divergerebbe senza limite.

Si può affermare che il modello matematico risponde bene.

Matteo Gentileschi

Lascia un commento