paralelne[Z1_,Z2_]:=((Z1*Z2)/(Z1+Z2));
  
Funkce pro tvorbu admitanční matice vycházející ze základní formule:
(Subscript[(OverHat[I]), k])=(Subscript[(OverHat[U]), k])(Subscript[(OverHat[Y]), kk])-(∑_(m=1
m≠k)^n)(Subscript[(OverHat[U]), m])(Subscript[(OverHat[Y]), km])
Tedy:
(({{({{({{(Subscript[(OverHat[I]), k])}, {(Subscript[(OverHat[I]), m])}})}, {..}})}, {(Subscript[(OverHat[I]), n])}}))=(({{({{(Subscript[(OverHat[Y]), kk]), -(Subscript[(OverHat[Y]), km])}, {-(Subscript[(OverHat[Y]), mk]), (Subscript[(OverHat[Y]), mm])}}), ({{..., -(Subscript[(OverHat[Y]), kn])}, {..., -(Subscript[(OverHat[Y]), mn])}})}, {({{..., ...}, {-(Subscript[(OverHat[Y]), nk]), -(Subscript[(OverHat[Y]), nm])}}), ({{..., ...}, {..., (Subscript[(OverHat[Y]), nn])}})}})).(({{({{({{(Subscript[(OverHat[U]), k])}, {(Subscript[(OverHat[U]), m])}})}, {..}})}, {(Subscript[(OverHat[U]), n])}}));
Ymatice[uzel1_,uzel2_,Z_,Y_,T_]:=(
Ypom=(1/Z);
Ykm=(Ypom/T);
Yk0=((1-T)/(T^2))*Ypom;
Ym0=(1-(1/T))*Ypom;
(*Nejprve je nutno vložit na diagnálu uzlové admitance*)
m=DiagonalMatrix[Y];
For[i=1,i<Length[Y],i++,
(*Dále pak admitance všech okolních větví*)
m[[uzel1[[i]],uzel1[[i]]]]=m[[uzel1[[i]],uzel1[[i]]]]+Ykm[[i]]+Yk0[[i]];
m[[uzel2[[i]],uzel2[[i]]]]=m[[uzel2[[i]],uzel2[[i]]]]+Ykm[[i]]+Ym0[[i]];
(*nediagonalni prvky*)
m[[uzel1[[i]],uzel2[[i]]]]=m[[uzel1[[i]],uzel2[[i]]]]-Ykm[[i]];
m[[uzel2[[i]],uzel1[[i]]]]=m[[uzel2[[i]],uzel1[[i]]]]-Ykm[[i]];
];
m
);
proudyVetvemiKM[Z_,U_,uzel1_,uzel2_]:=(
Ikm=Table[0,{Length[Z]}];
For[i=1,i<Length[Z]+1,i++,
Ikm[[i]]=((U[[uzel1[[i]]]]-U[[uzel2[[i]]]])/(Z[[i]]));
];
Ikm
);
GVA=10^9;
MVA=10^6;
kV=kA=10^3;
km=1;
μS=10^-6;
a=(E^(I*((2π)/3)));
Amatice={{1,1,1},
{a\[Conjugate],a,1},
{a,a\[Conjugate],1}};
(*Vstupní hodnoty*)
S1={Ss->1.GVA,Ik1s->4.9kA,Ik3s->5.02kA,Us->115kV};
V1={lv->10.km,Rv->0.125,Xv->0.407,Bv->2.817μS,X0X1->3};
V2={lv->0.2km,Rv->0.125,Xv->0.407,Bv->2.817μS,X0X1->3};
T1={U1->115.kV,U2->10.5,uk->0.1,Sn->10MVA};
G1={Sg->10.MVA,Ug->10.5kV,Xd1->0.198,Xd2->0.208,Xd0->0.059,Rn->100};
U1n=115kV;
U2n=10.5kV;
uzel1={1,2,3};
uzel2={2,3,4};

(*Výpočet parametrů*)
Uv=115kV;
Sv=1;
Iv=(Sv/((Sqrt[3])*Uv));
Zv=((Uv^2)/Sv);
(*Síť*)
Ze znalosti náhradních schémat trojfázového a jednofázového zkratu dopočteme náhradní impedance:
(Subscript[(OverHat[I]), k1])=3*((OverHat[E])/((OverHat[Z])1+(OverHat[Z])2+(OverHat[Z])0))*Iv
(Subscript[(OverHat[I]), k3])=((OverHat[E])/((OverHat[Z])1))*Iv
kde (OverHat[E])=(((OverHat[U])s)/Uv) a Iv=(Sv/((Sqrt[3])*Uv))
pak pokud Zv=((Uv^2)/Sv):
(OverHat[Z])1=((OverHat[E])/(Subscript[(OverHat[I]), k3]))*Iv=((((OverHat[U])s)/Uv)/(Subscript[(OverHat[I]), k3]))*(Sv/((Sqrt[3])*Uv))=(((OverHat[U])s)/(Subscript[(OverHat[I]), k3]))*(Sv/((Sqrt[3])*(Uv^2)))=(((OverHat[U])s)/((Sqrt[3])*(Subscript[(OverHat[I]), k3])))*(1/Zv)
(OverHat[Z])0=3*((OverHat[E])/(Subscript[(OverHat[I]), k1]))*Iv-((OverHat[Z])1+(OverHat[Z])2)=3*(((OverHat[U])s)/((Sqrt[3])*(Subscript[(OverHat[I]), k1])))*(1/Zv)-(OverHat[Z])1-(OverHat[Z])2
Protože se jedná o hladinu 115kV, uvažujeme pouze reaktance.;
Zs1=Zs2=I*(Us/(Ik3s*(Sqrt[3])))*(1/Zv)/.S1;
Zs0=I*3*(Us/(Ik1s*(Sqrt[3])))*(1/Zv)-Zs1-Zs2/.S1;
(*Vedení V1*)
Zv11=Zv12=lv*(Rv+I*Xv)*(1/Zv)/.V1;
Zv10=Zv11*X0X1/.V1;
(*Vedení V2*)
Zv21=Zv22=lv*(Rv+I*Xv)*(1/Zv)*((((U1n/U2n)))^2)/.V2;
Zv20=Zv21*X0X1/.V2;
(*trafo T1*)
Zt1=Zt2=Zt0=I*uk*((U1^2)/Sn)*(1/Zv)/.T1;
(*Generátor G1*)
Zg1=I*Xd1*((Ug^2)/Sg)*(1/Zv)*((((U1n/U2n)))^2)/.G1;
Zg2=I*Xd2*((Ug^2)/Sg)*(1/Zv)*((((U1n/U2n)))^2)/.G1;
Zg0=(3*Rn+I*Xd0*((Ug^2)/Sg))*(1/Zv)*((((U1n/U2n)))^2)/.G1;
(*viděné impedance v uzlu 2*)
Z1=paralelne[(Zs1+Zv11),(Zt1+Zv21+Zg1)];
Z2=paralelne[(Zs2+Zv12),(Zt2+Zv22+Zg2)];
Z0=paralelne[(Zs0+Zv10),(Zt0)];

Tvorba matic:;
Poznámka:
Ve funkci pro tvorbu admitanční matice je zakomponován také převod. Impedanční matice získaná inverzí z takové matice má v každém diagonálním prvku přesně viděnou impedanci daného uzlu, respektující rozdílné napěťové hladiny ostatních impedancí a není tedy třeba ostatní impedance přepočítávat na napěťovou hladinu v místě zkratu. Tato matice pak není symetrická a neplatí tedy Zkm=Zmk.
V našem příkladu počítáme zkratové příspěvky v uzlu 2. Protože převod v našich maticích nerespektujeme a T={1,1,1}, musíme impedance ležící za transformátorem (hladina U2n), přepočítat do hladiny napětí námi zkoumaného uzlu (hladina U1n). Takto získanou impedanční matici pak můžeme použít pouze pro výpočty v uzlech 1 a 2 na hladině U1n. Na druhou stranu se tak vyhneme složitému přepočítávání impedancí při výpočtech proudů ve větvích.;
Ymatice1=Ymatice[
uzel1,
uzel2,
{Zv11,Zt1,Zv21},
{(1/Zs1),0,0,(1/Zg1)},
{1,1,1}
];
Ymatice2=Ymatice[
uzel1,
uzel2,
{Zv12,Zt2,Zv22},
{(1/Zs2),0,0,(1/Zg2)},
{1,1,1}
];
Ymatice0=Ymatice[
uzel1,
uzel2,
{Zv10,Infinity,Zv20},
{(1/Zs0),(1/Zt0),0,(1/Zg0)},
{1,1,1}
];
Zmatice1=Inverse[Ymatice1];
Zmatice2=Inverse[Ymatice2];
Zmatice0=Inverse[Ymatice0];

Print["Ymatice1 = ",MatrixForm[Ymatice1]];
Print["Ymatice2 = ",MatrixForm[Ymatice2]];
Print["Ymatice0 = ",MatrixForm[Ymatice0]];
Print["Zmatice1 = ",MatrixForm[Zmatice1]];
Print["Zmatice2 = ",MatrixForm[Zmatice2]];
Print["Zmatice0 = ",MatrixForm[Zmatice0]];

Trojfázový zkrat:;
;
(*Hodnota na souřadnicích 2,2 Zmatice1 je viděnou impedancí v uzlu 2*)
Ik=((U1n/Uv)/Z1);
Print["Hodnota zkratového proudu: ",Abs[Ik]*Iv];
vektIk={0,Ik,0,0};
vektUk=((U1n/Uv)-Zmatice1.vektIk);
Print["Vektor napětí při poruše: ",MatrixForm[ vektUk*{1,1,(U2n/U1n),(U2n/U1n)}*Uv]//Chop];
Toky proudů:;
Ikm=proudyVetvemiKM[{Zv11,Zt1,Zv21},vektUk,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Print["vektory proudů ve fázi A: ",MatrixForm[Ikm]," jejich abs hodnoty: ",MatrixForm[Abs[Ikm]], " a úhly: ", MatrixForm[Arg[Ikm]*(180/π)]];
Print["Sečtením proudů v první a druhé větvi, dostaneme opět zkratový proud v uzlu 1: ",Abs[Ikm[[1]]-Ikm[[2]]]];
Print["Protože je trojfázový zkrat symetrický, je pro jeho výpočet použita pouze sousledná složka a velikosti proudů v ostatních fázích budou stejné, pouze otočené o příslušný úhel."]

Jednofazový zkrat ve fázi A:;
;
Ik=3*((U1n/Uv)/((Z1+Z2+Z0)));
Print["Hodnota zkratového proudu: ",Abs[Ik]*Iv];
(*Nesymetrický proud (je pouze ve fázi A!) v uzlu 2 musíme rozdělit do složek*)
vektIk120=Inverse[Amatice].{Ik,0,0};
vektIk1={0,vektIk120[[1]],0,0};
vektIk2={0,vektIk120[[2]],0,0};
vektIk0={0,vektIk120[[3]],0,0};
vektUk1=((U1n/Uv)-Zmatice1.vektIk1);
vektUk2=(-Zmatice2.vektIk2);
vektUk0=(-Zmatice0.vektIk0);
vektUkA=vektUk1+vektUk2+vektUk0;
vektUkB=a\[Conjugate]*vektUk1+a*vektUk2+vektUk0;
vektUkC=a*vektUk1+a\[Conjugate]*vektUk2+vektUk0;
Print["Vektor napětí souslené složky ve všech uzlech: ",MatrixForm[ vektUk1*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí zpětné složky ve všech uzlech: ",MatrixForm[ vektUk2*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí netočivé složky ve všech uzlech: ",MatrixForm[ vektUk0*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]//Chop];
Print["Vektor napětí ve fazi A ve všech uzlech: ",MatrixForm[ vektUkA*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]//Chop];
Print["Vektor napětí ve fazi B ve všech uzlech: ",MatrixForm[ vektUkB*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí ve fazi C ve všech uzlech: ",MatrixForm[ vektUkC*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Toky proudů:;
Ikm1=proudyVetvemiKM[{Zv11,Zt1,Zv21},vektUk1,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm2=proudyVetvemiKM[{Zv12,Zt2,Zv22},vektUk2,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm0=proudyVetvemiKM[{Zv10,Infinity,Zv20},vektUk0,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
IkmA=Ikm1+Ikm2+Ikm0;
IkmB=a\[Conjugate]*Ikm1+a*Ikm2+Ikm0;
IkmC=a*Ikm1+a\[Conjugate]*Ikm2+Ikm0;
Print["vektory proudů ve fázi A: ",MatrixForm[IkmA]," jejich abs hodnoty: ",MatrixForm[Abs[IkmA]], " a úhly: ", MatrixForm[Arg[IkmA]*(180/π)]];
Print["vektory proudů ve fázi B: ",MatrixForm[IkmB]," jejich abs hodnoty: ",MatrixForm[Abs[IkmB]], " a úhly: ", MatrixForm[Arg[IkmB]*(180/π)]];
Print["vektory proudů ve fázi C: ",MatrixForm[IkmC]," jejich abs hodnoty: ",MatrixForm[Abs[IkmC]], " a úhly: ", MatrixForm[Arg[IkmC]*(180/π)]];
Print["Protože poměr nulové složky je na každé straně transformátoru různý, sečtením proudů v první a druhé větvi dosteneme ",Abs[IkmA[[1]]-IkmA[[2]]]," a ne ",Abs[Ik]*Iv];

Dvoufázový zkrat ve fázich B a C:;
;
Ik2C=I*(Sqrt[3])*((U1n/Uv)/(Z1+Z2));
Ik2B=-Ik2C;
Print["Hodnota zkrat. proudu ve fázi C je :",Abs[Ik2C]*Iv]
Print["Hodnota zkrat. proudu ve fázi B je :",Abs[Ik2B]*Iv]
Print["Proudy ve fázích B a C jsou stejně veliké a opačné"];
Graphics[{Red,{Arrow[{{0,0},{Re[Ik2B],Im[Ik2B]}}],Green,Arrow[{{0,0},{Re[Ik2C],Im[Ik2C]}}]}}]
vektIk120=Inverse[Amatice].{0,Ik2B,Ik2C};
vektIk1={0,vektIk120[[1]],0,0};
vektIk2={0,vektIk120[[2]],0,0};
vektIk0={0,0,0,0};
vektUk1=((U1n/Uv)-Zmatice1.vektIk1);
vektUk2=(-Zmatice2.vektIk2);
vektUk0=(-Zmatice0.vektIk0);
vektUkA=vektUk1+vektUk2+vektUk0;
vektUkB=a\[Conjugate]*vektUk1+a*vektUk2+vektUk0;
vektUkC=a*vektUk1+a\[Conjugate]*vektUk2+vektUk0;
Print["Vektor napětí souslené složky ve všech uzlech: ",MatrixForm[ vektUk1*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí zpětné složky ve všech uzlech: ",MatrixForm[ vektUk2*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí netočivé složky ve všech uzlech: ",MatrixForm[ vektUk0*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))], "musí být nulový"];
Print["Vektor napětí ve fazi A ve všech uzlech: ",MatrixForm[ vektUkA*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí ve fazi B ve všech uzlech: ",MatrixForm[ vektUkB*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí ve fazi C ve všech uzlech: ",MatrixForm[ vektUkC*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Toky proudů:;
Ikm1=proudyVetvemiKM[{Zv11,Zt1,Zv21},vektUk1,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm2=proudyVetvemiKM[{Zv12,Zt2,Zv22},vektUk2,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm0=proudyVetvemiKM[{Zv10,Infinity,Zv20},vektUk0,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
IkmA=Ikm1+Ikm2+Ikm0;
IkmB=a\[Conjugate]*Ikm1+a*Ikm2+Ikm0;
IkmC=a*Ikm1+a\[Conjugate]*Ikm2+Ikm0;
Print["vektory proudů ve fázi A: ",MatrixForm[IkmA]," jejich abs hodnoty: ",MatrixForm[Abs[IkmA]], " a úhly: ", MatrixForm[Arg[IkmA]*(180/π)]];
Print["vektory proudů ve fázi B: ",MatrixForm[IkmB]," jejich abs hodnoty: ",MatrixForm[Abs[IkmB]], " a úhly: ", MatrixForm[Arg[IkmB]*(180/π)]];
Print["vektory proudů ve fázi C: ",MatrixForm[IkmC]," jejich abs hodnoty: ",MatrixForm[Abs[IkmC]], " a úhly: ", MatrixForm[Arg[IkmC]*(180/π)]];
Print["Sečtením proudu v první a druhé větvi, dostaneme zkratový proud v uzlu 1: ",Abs[IkmB[[1]]-IkmB[[2]]]];

Dvoufázový zemní zkrat ve fázích B a C:;
;
I1=((U1n/Uv)/(Z1+paralelne[Z2,Z0]));
Ik2nB=(((a\[Conjugate]-1)Z2+(a\[Conjugate]-a)Z0)/(Z2+Z0))*I1;
Ik2nC=(((a-1)Z2+(a-a\[Conjugate])Z0)/(Z2+Z0))*I1;
Print["Hodnota zkrat. proudu ve fázi B je :",Abs[Ik2nB]*Iv," a ve fázi C :",Abs[Ik2nC]*Iv];
vektIk120=Inverse[Amatice].{0,Ik2nB,Ik2nC};
vektIk1={0,vektIk120[[1]],0,0};
vektIk2={0,vektIk120[[2]],0,0};
vektIk0={0,vektIk120[[3]],0,0};
vektUk1=((U1n/Uv)-Zmatice1.vektIk1);
vektUk2=(-Zmatice2.vektIk2);
vektUk0=(-Zmatice0.vektIk0);
vektUkA=vektUk1+vektUk2+vektUk0;
vektUkB=a\[Conjugate]*vektUk1+a*vektUk2+vektUk0;
vektUkC=a*vektUk1+a\[Conjugate]*vektUk2+vektUk0;
Print["Vektor napětí souslené složky ve všech uzlech: ",MatrixForm[ vektUk1*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí zpětné složky ve všech uzlech: ",MatrixForm[ vektUk2*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí netočivé složky ve všech uzlech: ",MatrixForm[ vektUk0*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]//Chop];
Print["Vektor napětí ve fazi A ve všech uzlech: ",MatrixForm[ vektUkA*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]];
Print["Vektor napětí ve fazi B ve všech uzlech: ",MatrixForm[ vektUkB*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]//Chop];
Print["Vektor napětí ve fazi C ve všech uzlech: ",MatrixForm[ vektUkC*{1,1,(U2n/U1n),(U2n/U1n)}*(Uv/(Sqrt[3]))]//Chop];
Toky proudů:;
Ikm1=proudyVetvemiKM[{Zv11,Zt1,Zv21},vektUk1,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm2=proudyVetvemiKM[{Zv12,Zt2,Zv22},vektUk2,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
Ikm0=proudyVetvemiKM[{Zv10,Infinity,Zv20},vektUk0,uzel1,uzel2]*{1,1,(U1n/U2n)}*Iv;
IkmA=Ikm1+Ikm2+Ikm0;
IkmB=a\[Conjugate]*Ikm1+a*Ikm2+Ikm0;
IkmC=a*Ikm1+a\[Conjugate]*Ikm2+Ikm0;
Print["vektory proudů ve fázi A: ",MatrixForm[IkmA]," jejich abs hodnoty: ",MatrixForm[Abs[IkmA]], " a úhly: ", MatrixForm[Arg[IkmA]*(180/π)]];
Print["vektory proudů ve fázi B: ",MatrixForm[IkmB]," jejich abs hodnoty: ",MatrixForm[Abs[IkmB]], " a úhly: ", MatrixForm[Arg[IkmB]*(180/π)]];
Print["vektory proudů ve fázi C: ",MatrixForm[IkmC]," jejich abs hodnoty: ",MatrixForm[Abs[IkmC]], " a úhly: ", MatrixForm[Arg[IkmC]*(180/π)]];
Print["Protože poměr nulové složky je na každé straně transformátoru různý, sečtením proudů v první a druhé větvi dosteneme ",Abs[IkmB[[1]]-IkmB[[2]]]," a ne ",Abs[Ik2nB]*Iv];