How I can plot this equation?
1 view (last 30 days)
Show older comments
Franco Belletti
on 16 Dec 2018
Commented: Franco Belletti
on 16 Dec 2018
TB =3.8470608473110835049785363403119;
TI =15.448004476144693589582536852599;
n20av=1.1;
%confinement time [s]
tauE=1.2;
%internal radius [m]
a=2.0;
%plasma ellipticity
k=1.8;
%external radius [m]
R0=6.2;
%inverse aspect rdius
epsy=a/R0;
%plasma volume [m3]
Vp=2*pi*R0*pi*k*a^2;
%plasma power [W/m3]
%magnetic current [MA]
Imax=7.5;
%homic power [W/m3]
%bremmstraulung power [W/m3]
Sb1=6.14*(10^3)*(n20av^2)*TB.^(1/2);
%conduction power
Sk2=5.34*(10^4)*n20av*TI*1/tauE;
T1=0.1;
sigmav1=(10^(-6))*exp(-21.38*(T1.^(-0.2935))-25.20-7.101*(10^(-2))*T1+1.938*(10^(-4))*T1.^2+4.925*(10^(-6))*T1.^(3)-3.984*(10^(-8))*T1.^4);
sigmavn1=sigmav1/10^(-22);
Sa1=(2.31*10^(5))*(n20av^(2))*sigmavn1;
Sohm1=(1/Vp)*(5.6*10^(-2)*((1-1.31*epsy^(1/2)+0.46*epsy)).^-1)*(R0*Imax^(2))*((a^(2)*k*((T1).^(3/2))).^-1);
dwdt = Sa1 + Sohm1-Sb1-Sk2;
syms Tk
W_t=(5.34*10^(4))*n20av*Tk;
bil = W_t +dwdt ==0;
w=vpasolve(bil,Tk);
T=linspace(0,0.5,200);
sigmav1=(10^(-6))*exp(-21.38*(T.^(-0.2935))-25.20-7.101*(10^(-2))*T+1.938*(10^(-4))*T.^2+4.925*(10^(-6))*T.^(3)-3.984*(10^(-8))*T.^4);
sigmavn1=sigmav1/10^(-22);
Sa1=(2.31*10^(5))*(n20av^(2))*sigmavn1;
Sohm1=(1/Vp)*(5.6*10^(-2)*((1-1.31*epsy^(1/2)+0.46*epsy)).^-1)*(R0*Imax^(2))*((a^(2)*k*((T).^(3/2))).^-1);
dwdt = Sa1 + Sohm1-Sb1-Sk2;
bill = w+dwdt;
for T=linspace(1,16,400)
if (TB< T) && T <= TI
Saux = (1-((T-TB)/(TI-TB)));
bill1= bill+Saux;
elseif T>TI
Saux=0;
end
end
plot(T,bill1)
0 Comments
Accepted Answer
Mark Sherstan
on 16 Dec 2018
Edited: Mark Sherstan
on 16 Dec 2018
plot(bill1)
You are currently getting an error as T is only 1x1 and bill1 is 1x200
More Answers (0)
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!