Simulation des circuits linéaires

1.a. Transmittance de Laplace

On considère un circuit linéaire comportant N noeuds et on note V_i le potentiel électrique du noeud i. Le circuit comporte une source de tension idéale (l'entrée du filtre) dont le potentiel sera pris par convention égal à 1. Par définition, la masse est un noeud du circuit dont le potentiel est nul.

Pour un noeud k auquel sont attachés des dipôles linéaires passifs, la loi des noeuds s'écrit :

$$\left(\sum _i{Y}_i\right)V_k-\sum _i{Y}_i{V}_i=0$$

où la somme porte sur les noeuds reliés au noeud k par un dipôle d'admittance Y_i. Les admittances dont définies comme des fonctions de la variable s=jω.

Le système d'équations linéaires pour les potentiels se met sous la forme AV+B=0 où A est une matrice $N\times N$ de fractions rationnelles, V et B deux vecteurs colonnes. La ligne k de la matrice A correspond au noeud k.

L'ajout d'une source de tension continue (constante) entre les noeuds na et nb introduit pour ces deux noeuds l'équation :

$$V_{na}-V_{nb}=U$$

Une source de tension commandée en tension (STCT) entre les noeuds na et nb, commandée par la la d.d.p. entre les noeuds nx et ny, introduit l'équation :

$$V_{na}-V_{nb}=g(V_{nx}-V_{ny})$$

pour les noeuds na et nb.

Si la sortie du filtre est la tension au noeud ns, la transmittance de Laplace du filtre est le potentiel de ce noeud :

$$h\left(s\right)=V_{ns}\left(s\right)$$

1.b. Réponse fréquentielle

La réponse en régime sinuoïdal est donnée par la fonction de transfert :

$$H\left(\omega \right)=h\left(j\omega \right)$$

1.c. Réponse temporelle

On suppose qu'à l'instant initial t=0 les condensateurs sont déchargés et les inductances n'ont pas de courant. Dans ce cas, la transformée de Laplace du signal de sortie est :

$$L\left(s\right)=h\left(s\right)F\left(s\right)$$

où $F\left(s\right)$ est la transformée de Laplace du signal d'entrée. Pour un échelon de tension appliqué en entrée, on a :

$$F\left(s\right)=\frac{1}{s}$$

La réponse impulsionnelle s'obtient avec :

$$F\left(s\right)=1$$

2.a. Définition du système

Les fonction suivantes transforment soit la matrice A seulement, soit les matrices A et B. Les sources de tension doivent être définies après les dipôles passifs car elles désactivent la loi des noeuds.

function [P]=ajouter_dipole(A,na,nb,Y)
  P=A
  P(na,na)=P(na,na)+Y
  P(na,nb)=P(na,nb)-Y
  P(nb,nb)=P(nb,nb)+Y
  P(nb,na)=P(nb,na)-Y
endfunction                
                

Ajout de dipôles passifs particuliers :

function [P]=ajouter_resistance(A,na,nb,R)
  P=ajouter_dipole(A,na,nb,1/R)
endfunction

function [P]=ajouter_capacite(A,na,nb,C)
  P=ajouter_dipole(A,na,nb,poly(0,'s')*C)
endfunction

function [P]=ajouter_inductance(A,na,nb,L)
  P=ajouter_dipole(A,na,nb,1/(poly(0,'s')*L))
endfunction
                

Ajout d'une source de tension continue :

function [P,Q]=ajouter_source_tension_continue(A,B,na,nb,U)
  P=A
  Q=B
  N=size(P,'c')
  P(na,:)=zeros(1,N)
  P(nb,:)=zeros(1,N)
  P(na,na)=1
  P(na,nb)=-1
  Q(na)=-U
  P(nb,nb)=-1
  P(nb,na)=1
  Q(nb)=-U
endfunction
                

Ajout d'une source de tension commandée en tension :

function [P,Q]=ajouter_source_STCT(A,B,na,nb,nca,ncb,g)
  P=A
  Q=B
  N=size(P,'c')
  P(na,:)=zeros(1,N)
  P(nb,:)=zeros(1,N)
  P(na,na)=1
  P(na,nb)=-1
  P(na,nca)=-g
  P(na,ncb)=g
  Q(na)=0
  P(nb,nb)=-1
  P(nb,na)=1
  P(nb,nca)=-g
  P(nb,ncb)=g
  Q(nb)=0
endfunction
                

Ajouter une masse, c.a.d. annuler un potentiel :

function [P,Q]=ajouter_masse(A,B,n0)
  P=A
  Q=B
  P(n0,:)=zeros(1,size(P,'c'))
  P(n0,n0)=1
  Q(n0)=0
endfunction
                

Si ne désigne le noeud de l'entrée du filtre, on pose V_ne=1. De cette manière, le potentiel V_k donne directement la fonction de transfert obtenue en plaçant la sortie sur le noeud k :

function [P,Q]=definir_entree(A,B,ne)
  P=A
  Q=B
  N=size(P,'c')
  P(ne,:)=zeros(1,N)
  P(ne,ne)=1
  Q(ne)=-1
endfunction
                

2.b. Réponse fréquentielle

Pour chaque pulsation ω, la matrice de fractions rationnelles A est évaluée avec la fonction horner, puis le système est résolu.

function [w,db,phi]=reponse(A,B,noeud,lgmin,lgmax,np)
  w=logspace(lgmin,lgmax,np)
  H=w
  for i=1:np,
    M=horner(A,%i*w(i))
    [x0,ker]=linsolve(M,B)
    H(i)=x0(noeud)
  end;
  [db,phi]=dbphi(H)
endfunction
                    

2.c. Exemple

Soit le circuit RLC suivant :

Figure pleine page

La fonction suivante définie le circuit :

getf('simulineaire.sci');
function [A,B]=circuitRLC(R,L,C)
  N=4
  A=zeros(N,N)
  B=zeros(N,1)
  A=ajouter_resistance(A,1,2,R)
  A=ajouter_inductance(A,2,3,L)
  A=ajouter_capacite(A,3,4,C)
  [A,B]=ajouter_masse(A,B,4)
  [A,B]=definir_entree(A,B,1)
endfunction
            

Définition du circuit et tracé du diagramme de Bode pour le noeud 3

[A,B]=circuitRLC(10,1e-3,1e-6);
[w,db,phi]=reponse(A,B,3,2,6,200);
            

plot1=scf();
plot2d(w,db,style=5,logflag='ln');
xtitle('Circuit RLC','log f','GdB');
            

plot1.pdf

plot2=scf();
plot2d(w,phi,style=5,logflag='ln');
xtitle('Circuit RLC','log f','phi');
            

plot2.pdf

3.a. Définition du système

ajouterDipole[A_,na_,nb_,Y_]:=Module[{P},
    P=A;
    P[[na,na]]=P[[na,na]]+Y;
    P[[na,nb]]=P[[na,nb]]-Y;
    P[[nb,nb]]=P[[nb,nb]]+Y;
    P[[nb,na]]=P[[nb,na]]-Y;
    Return[P];
]

ajouterResistance[A_,na_,nb_,R_]:=ajouterDipole[A,na,nb,1/R];
ajouterCapacite[A_,na_,nb_,C_]:=ajouterDipole[A,na,nb,C*s];
ajouterInductance[A_,na_,nb_,L_]:=ajouterDipole[A,na,nb,1/(L*s)];

ajouterSourceTensionContinue[A_,B_,na_,nb_,U_]:=Module[{P,Q,d,n},
    P=A;
    Q=B;
    d = Dimensions[A];
    n = d[[2]];
    P[[na]]=Table[0,{n}];
    P[[nb]]=Table[0,{n}];
    P[[na,na]]=1;
    P[[na,nb]]=-1;
    P[[nb,nb]]=-1;
    P[[nb,na]]=1;
    Q[[nb]]=-U;
    Return[{P,Q}];
]

ajouterSourceTensionSTCT[A_,B_,na_,nb_,nx_,ny_,g_]:=Module[{P,Q,d,n},
    P=A;
    Q=B;
    d = Dimensions[A];
    n = d[[2]];
    P[[na]]=Table[0,{n}];
    P[[nb]]=Table[0,{n}];
    P[[na,na]]=1;
    P[[na,nb]]=-1;
    P[[nb,nb]]=-1;
    P[[nb,na]]=1;
    P[[na,nx]]=-g;
    P[[na,ny]]=g;
    P[[nb,nx]]=-g;
    P[[nb,ny]]=g;
    Q[[nb]]=0;
    Return[{P,Q}];
]

ajouterMasse[A_,B_,n0_]:=Module[{P,Q,d,n},
    P=A;
    Q=B;
    d = Dimensions[A];
    n = d[[2]];
    P[[n0]]=Table[0,{n}];
    P[[n0,n0]]=1;
    Q[[n0]]=0;
    Return[{P,Q}];
]

definirEntree[A_,B_,ne_]:=Module[{P,Q,d,n},
    P=A;
    Q=B;
    d = Dimensions[A];
    n = d[[2]];
    P[[ne]]=Table[0,{n}];
    P[[ne,ne]]=1;
    Q[[ne]]=-1;
    Return[{P,Q}];
]
                

3.b. Transmittance

Transmittance de Laplace pour une sortie sur le noeud ns :

transfert[A_,B_,ns_]:=Module[{v},
    v=LinearSolve[A,-B];
    Return[v[[ns]]];
]
                

3.c. Réponse fréquentielle

bodeGain[h_,lgmin_,lgmax_,Gmin_,Gmax_]:=Module[{},
    Return[Plot[20*Log[10,Abs[h/.{s->I*2*Pi*10^p}]],{p,lgmin,lgmax},PlotRange->{{lgmin,lgmax},{Gmin,Gmax}},AxesLabel->{"log f","GdB"}]];
]
bodePhase[h_,lgmin_,lgmax_]:=Module[{},
    Return[Plot[180/Pi*Arg[h/.{s->I*2*Pi*10^p}],{p,lgmin,lgmax},AxesLabel->{"log f","phi"}]];
]
                

3.d. Réponse à un échelon

La fonction suivante prend en argument la transmittance de Laplace et le temps maximal, et renvoit la courbe de réponse temporelle :

reponseEchelon[h_,tmax_,vmin_,vmax_]:=Module[{r},
    r=InverseLaplaceTransform[h/s,s,t];
    Return[Plot[r/.{t->temps},{temps,0,tmax},PlotRange->{{0,tmax},{vmin,vmax}},AxesLabel->{"t","V"}]];
]
                

On définit de même la réponse imulsionnelle :

reponseImpulsion[h_,tmax_,vmin_,vmax_]:=Module[{r},
    r=InverseLaplaceTransform[h,s,t];
    Return[Plot[r/.{t->temps},{temps,0,tmax},PlotRange->{{0,tmax},{vmin,vmax}},AxesLabel->{"t","V"}]];
]
                

3.e. Exemple

On reprend l'exemple traité plus haut :

Get['simulineaire.m']; 
circuitRLC[R_,L_,C_]:=Module[{A,B,n},
    n=4;
    A=Table[0,{n},{n}];
    B=Table[0,{n}];
    A=ajouterResistance[A,1,2,R];
    A=ajouterInductance[A,2,3,L];
    A=ajouterCapacite[A,3,4,C];
    {A,B}=ajouterMasse[A,B,4];
    {A,B}=definirEntree[A,B,1];
    Return[{A,B}];
]
{A,B}=circuitRLC[10,10^-3,10^-6];
                

Fonction de transfert pour le noeud 3 :

h=transfert[A,B,3]

1000000000/(1000000000 + 10000*s + s^2)

Diagramme de Bode :

bodeGain[h,2,6,-80,20]

plot3.pdf

bodePhase[h,2,6]

plot4.pdf

Réponse à un échelon :

reponseEchelon[h,2*10^-3,0,2]

plot5.pdf

Le calcul peut être fait de manière formelle :

{A,B}=circuitRLC[R,L,C];
h=transfert[A,B,3]/.s->I*omega 
                

$$\frac{1}{-CL{\omega }^2+iC\omega R+1}$$