1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
|
L=(10^(-6))*[0.805775 0.538783; 0.538783 1.07757];
C=(10^(-12))*[134.693 -67.3467;-67.3467 97.6102];
Lg=400;
vit=1.20476*10^8;%inv(sqrt(L(2,2)*C(2,2)))
dz=200;
dt=dz/vit;
Rs=[50 0;0 50];
RL=[50 0;0 50];
Nx=40;%Lg/dz;
itmax=1000;
nc=2;
a1=inv(dz/dt*Rs*C+1);
a2=dz/dt*Rs*C-1;
a3=dt*inv(dz*L);
a4=inv(dz/dt*RL*C+1);
a5=dz/dt*RL*C-1;
a6=dt*inv(dz*C);
%INITIALISATION it = 1
V(1:nc,1:Nx+1,1)=0;
I(1:nc,1:Nx,1)=0;
Vs(1:nc,1)=0;
for it=2:itmax;
if it<=20
Vs(1,it)=0.05*(it-1);
else
Vs(1,it)=1;
end
end
for it=2:itmax
V(1:nc,1,it)=a1*(a2*V(1:nc,1,it-1)-2*Rs*I(1:nc,1,it-1)+(Vs(1:nc,it-1)+Vs(1:nc,it)));
V(1:nc,Nx+1,it)=a4*(a5*V(1:nc,Nx+1,it-1)+2*RL*I(1:nc,Nx,it-1));
for k=2:Nx
V(1:nc,k,it)=V(1:nc,k,it-1)-a6*(I(1:nc,k,it-1)-I(1:nc,k-1,it-1));
end
for k=1:Nx
I(1:nc,k,it)=I(1:nc,k,it-1)-a3*(V(1:nc,k+1,it)-V(1:nc,k,it));
end
end
subplot(2,1,1),plot(Vs(1,1:itmax)),title('Vs')
subplot(2,1,2),plot(V(1,1,:)) ,title('V1(1)') % la tension V du 1er conducteur, à l'endroit 1, et à tt les instants |
Partager