(*Zacatek*)

<<VectorAnalysis`
SetCoordinates[Cylindrical[r,ϕ,z]];
r1=0.115;(*poloměr induktoru [m]*)
r2=0.0085;(*poloměr vsázky [m]*)
omega=2*Pi*4000;(*frekvence [(s^(-1))]*)
delkaPrim=0.35;(*delka induktoru [m]*)
delkaSekund=0.2;(*delka vsázky [m]*)
alfaPrim=0.7;(*součinitel tvaru induktoru [-]*)
alfaSekund=0.3;(*součinitel tvaru vsázky [-]*)
k=1.2;(*vzrůst odporu vinutí s více zavity [-]*)
miVak=4×3.14×0.0000001;(*permeabilita vakua [(H/m)]*)
F=0.470;(*součinitel vzájemné indukčnosti [-]*)
gamaPrim=59*10^6;(*elektrická vodivost induktoru[(S/m)] *)
alfaPrestup=1;(*koeficient přestupu tepla se vsázky[(W/((m^2)*K))] *)
miPrim=1.002;(*permeabilita induktoru [-]*)
L1=0.000133;(*indukčnost generátoru [H]*)
C1=0.000012;(*kapacita příčného kondenzátoru [C]*)
C2=0.000015;(*kapacita kondenzátorové baterie [C]*)
U=400;(*napětí generátoru [V]*)
koefOdporuMed=0.0004;(*teplotní koeficient odporu [(K^(-1))]*)
sigma=5.67×(10^(-8));(*stefan-boltzmanova konstanta[(W/((m^2)*(K^(-4))))]*)
zavit=13;(*pocet zavitu ohrivaci civky [-]*)

proud2[{R1_,R2_,L2_},C2_]:=Module[{prvniSmycka,druhaSmycka,I1,I2},

prvniSmycka=U+I*omega*L1*I1+(I1/(I*omega*C1))+((I1-I2)/(I*omega*C2))==0;
druhaSmycka=I*omega*L2*I2+((I2-I1)/(I*omega*C2))+(R1+R2)*I2==0;

I2/.(Solve[{prvniSmycka,druhaSmycka},{I1,I2}][[1]])]

fceMir2[T_,miDvacetStupnu_,Tcurie_]:=(miDvacetStupnu/π*(-ArcTan[0.1*(T-Tcurie)]+Pi/2)+1)(*funkce relativní permeability na teplotě*)

Plot[fceMir2[T,1000,1000],{T,300,1300}]

fceGama[T_,gamaDvacetStupnu_,koefOdporu_]:=(gamaDvacetStupnu/((1+(T-293)*koefOdporu)))(*funkce elektrické vodivosti na teplotě*)

obvParamVsazky[miSekund_,gamaSekund_]:=Module[{hloubkaPrim,hloubkaSekund,elPrumPrim,elPrumSekund,indukcnostJedenZavitPrim,indukcnostSekund,vzajemInduknostJedenZavit,prevodJedenZavit,odporCelkemJedenZavit,odporVsazky,odporJedenZavitPrim,odporSekund,indukcnostCelkemJedenZavit},
hloubkaPrim=(Sqrt[(2/(omega×gamaPrim×miVak×miPrim))]);
hloubkaSekund=(Sqrt[(2/(omega×gamaSekund×miVak×miSekund))]);
elPrumPrim=2*r1+hloubkaPrim;
elPrumSekund=2*r2-hloubkaSekund;
odporJedenZavitPrim=k×((π×elPrumPrim)/(gamaPrim×hloubkaPrim×delkaPrim));
odporSekund=((π×elPrumSekund)/(gamaSekund×hloubkaSekund×delkaSekund));
indukcnostJedenZavitPrim=miVak×miPrim×(((π×(elPrumPrim^2))/4)/delkaPrim)×alfaPrim;
indukcnostSekund=miVak×(((π×(elPrumSekund^2))/4)/delkaSekund)×alfaSekund;
vzajemInduknostJedenZavit=miVak×(((π×(elPrumSekund^2))/4)/delkaSekund)×F;
prevodJedenZavit=(Sqrt[(((omega^2)×(vzajemInduknostJedenZavit^2))/((odporSekund^2)+(omega^2)×(indukcnostSekund^2)))]);
odporVsazky=(((prevodJedenZavit))^2)×odporSekund;
odporCelkemJedenZavit=odporJedenZavitPrim+(((prevodJedenZavit))^2)×odporSekund;
indukcnostCelkemJedenZavit=indukcnostJedenZavitPrim-(((prevodJedenZavit))^2)×indukcnostSekund;
{odporJedenZavitPrim×(zavit^2),(zavit^2)*odporVsazky,(zavit^2)*indukcnostCelkemJedenZavit}]

proudovaHustota2[gama_,miR_,omega_,proud_]=Block[{gama,miR,omega,proud},proud*(Exp[1000*r])/Integrate[(2*π*r*(Exp[1000*r])),{r,0,r2},Assumptions->gama>0&&miR>0&&omega>0](*rozložení proudu ve vsázce*)
]

gamaDvacetStupnu=37.8*10^6;(*elektricka vodivost vsazky pri 20 stupnich[(S/m)]*)
koefOdporu=6.5*10^-3;(*teplotni koeficient odporu vsazky [(K^(-1))]*)
mirDvacetStupnu=500; (*relativni permeabilita pri 20 sttupnich [-]*)
Tcurie=780+273; (*Curieho teplota vsazky [K]*)
lamda=100;(*tepelna vodivost vsazky []*)
ro=7000;(*hustota vsazky [](kg/(m^3))*)
CpCo=400;(* merna tepelna kapacita vsazky [(J/(K*kg))]*)
Tokoli=293;(*teplota okoli [K]*)
Tpocatecni=293;(*pocatecni teplota vsazky [K]*)
emisivita=0.7;(*emisivita [-]*)

(*Plne sdruzene*)

scenar={{1000,0.000015,True},{1000,0.00025,True},{1000,0.00015,False}}(*{{Δt,C2,zapnuto}..}*)

resJedenPrubehT[{casPocatecni_,casKoncovy_,TpocatecniFce_,C2_,zapnuto_}]:=
Module[{rovTepla,okrVpravo,okrVlevo,pocPodm,Qvivin,gradT,resT},
Print[{casPocatecni,casKoncovy,TpocatecniFce[r2/2]}];
gradT=Grad[T[r,t]][[1]]/.r->r2;
okrVpravo=-lamda×gradT==sigma×emisivita×(((T[r2,t])^4)-(Tokoli^4))-alfaPrestup×2×Pi×r2×delkaSekund×(T[r2,t]-Tokoli);(*odvod tepla na rozhraní*)
okrVlevo=Derivative[1,0][T][0.00001,t]==0;
pocPodm=T[r,casPocatecni]==TpocatecniFce[r];
Qvivin[r_,t_]:=((Abs[proudovaHustota2[fceGama[T[r,t],gamaDvacetStupnu,koefOdporu],fceMir2[T[r,t],mirDvacetStupnu,Tcurie],omega,proud2[obvParamVsazky[fceMir2[T[r,t],mirDvacetStupnu,Tcurie],fceGama[T[r,t],gamaDvacetStupnu,koefOdporu]],C2]]])^2)/fceGama[T[r,t],gamaDvacetStupnu,koefOdporu]/;zapnuto==True;(*výkon ve vsázce*)
Qvivin[r_,t_]:=0/;zapnuto==False;
rovTepla=lamda×Div[Grad[T[r,t]]]+Qvivin[r,t]==ro×CpCo×(∂_t)T[r,t];(*tepelná rovnice*)
resT=NDSolve[{rovTepla,okrVpravo,okrVlevo,pocPodm},T,{r,0.00001,r2},{t,casPocatecni,casKoncovy}][[1]];
(*řešení tepelné rovnice*)
(*reseniT=(T/.resT);*)
T/.resT
]

resScenarT[scenar:{{_,_,_}..}]:=Module[{novyScenar,resMinule,cas,cas0},
cas=0;
resMinule=Tpocatecni&;
Table[cas0=cas;
cas=cas+scenar[[ii,1]];resMinule=resJedenPrubehT[{cas0,cas,resMinule[#,cas0]&,scenar[[ii,2]],scenar[[ii,3]]}],{ii,Length[scenar]}]
]

reseniScenarT=resScenarT[scenar]

reseniScenarT[[1]][r2/2,100]

cas=0;
plots=Table[
cas0=cas;
cas=cas+scenar[[ii,1]];
Plot[reseniScenarT[[ii]][r2/2,t],{t,cas0,cas}],
{ii,Length[scenar]}];
Show[plots,PlotRange->All]

cas=0;
plots=Table[
cas0=cas;
cas=cas+scenar[[ii,1]];
Plot3D[reseniScenarT[[ii]][r,t],{r,0,r2},{t,cas0,cas}],
{ii,Length[scenar]}];
Show[plots,PlotRange->All]

cas=0;
Table[{cas,cas=cas+scenar[[ii,1]],scenar[[ii,2]],scenar[[ii,3]]},{ii,Length[scenar]}]