MATLAB

% REFERENCES:
% https://www.mathworks.com/matlabcentral/fileexchange/70244-polar-phasors
clearvars

Plot the phase to neutral line voltages for two fundamental cycles on the same graph

freq = 50;

[vtAArr, vtBArr, vtCArr, tArr] = getThreePhiSinusoid(freq, 40e-3, 120*sqrt(3));

figure()

hold on

plot(tArr*1e3, vtAArr)

plot(tArr*1e3, vtBArr)

plot(tArr*1e3, vtCArr)

legend('v_{t-A}', 'v_{t-B}', 'v_{t-C}')

grid on

grid minor

xlabel('Time (ms)')

ylabel('V_{phase-neutral} (V)')

title("Phase to Neutral Line Voltages vs. Time Graph")

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

saveas(gcf, "E.3.a.grid.emf")

saveas(gcf, "E.3.a.grid.png")

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fit a polynomial to the 𝐼𝑓 vs 𝐸𝑎𝑓 data obtained from the experiment

% Experimental Data

IfArr = [0.06 0.16 0.18 0.22 0.29 0.36 0.42];

EfArr = [75 176 210 250 320 380 420]/sqrt(3); % divide by sqrt(3) to get line-neutral values

% Polynomial fitting

coeffs = polyfit(IfArr, EfArr, 2)

coeffs = 1×3

 -380.7657  747.3807   -2.2177

fittedPoly = @(x) (coeffs(1)*x.^2 + coeffs(2)*x + coeffs(3));

% Creating virtual data using fitted polynomial

IfApprox = linspace(0, 1, 1e5);

EfApprox = fittedPoly(IfApprox);

% plotting

figure()

hold on

plot(IfArr, EfArr, "-*")

plot(IfApprox, EfApprox)

legend('Experimental', 'Approximated')

grid on

grid minor

xlabel('I{f} (A)')

ylabel('E_{af} (V)')

title("I{f} vs. E{af} Characteristics")

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

saveas(gcf, "E.3.b.grid.emf")

saveas(gcf, "E.3.b.grid.png")

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Calculating 𝐸𝑎𝑓 and corresponding 𝐼𝑓 such that S = 150 kVA at unity power factor is injected to the grid at the ABC points.

S_mag = 150e3; % by definition
VtA = 120; % by definition
pf = 1; % by definition
R_armature = 0.7/500; % measured
R_field = 190/500; % measured
xg = j*0.2;
xs = j*0.1;
ztot = R_armature + xg + xs;
ztotPolar = getPolar(ztot);
phi = acos(pf); % radians
iA_mag = (S_mag/3)/VtA;
iA = iA_mag*cos(phi) - j*sin(phi);
eaf = VtA + iA*ztot
eaf = 1.2058e+02 + 1.2500e+02i
eafPolar = getPolar(eaf)
eafPolar = "173.682∠46.030°"
corrIf = immersiveFinder(EfArr, IfArr, abs(eaf))
corrIf = 0.2900
Ploss = 3 *( ((abs(iA))^2)*R_armature + (corrIf^2)*R_field ) % I^2R_arm + If^2R_field
Ploss = 729.2625
P = S_mag*abs(pf);
efficiency = (P/(P + Ploss))*100
efficiency = 99.5162
maximumAllovedVpcc = VtA*1.10 % %10 more of rated value
maximumAllovedVpcc = 132
vpcc = VtA + iA*xg
vpcc = 1.2000e+02 + 8.3333e+01i
vpccPolar = getPolar(vpcc) % yeah, it exceeds the allowed range
vpccPolar = "146.097∠34.778°"

Plotting the 𝐸𝑎𝑓, 𝑉𝑃𝐶𝐶, 𝑉𝑡, 𝐼𝑎 phasors

% I found that beautiful function from:

% https://www.mathworks.com/matlabcentral/fileexchange/70244-polar-phasors

spplot(getPolarVal(eaf), getPolarVal(vpcc), ...

getPolarVal(VtA), getPolarVal(iA) )

Plotting the 𝐸𝑎𝑓, 𝑉𝑃𝐶𝐶, 𝑉𝑡, 𝐼𝑎 phasors in time domain

% plotting

[eafSig, tSig] = getSigFromPhasor(eaf, 60e-3);

[vpccSig, ~] = getSigFromPhasor(vpcc, 60e-3);

[vtASig, ~] = getSigFromPhasor(VtA, 60e-3);

[iASig, ~] = getSigFromPhasor(iA, 60e-3);

figure()

hold on

yyaxis left

plot(tSig*1e3, eafSig)

plot(tSig*1e3, vpccSig,"b-")

plot(tSig*1e3, vtASig, "c-")

grid on

grid minor

ylabel('Voltage (V)')

xlabel('Time (ms)')

title('E_{af}, V_{pcc}, V_t and I_{a} phasors versus time plot')

yyaxis right

plot(tSig*1e3, iASig)

ylabel('Current (A)')

legend('E_{af}', 'V_{pcc}', 'V_t', 'I_a')

hold off

Let us calculate the price for 10 years opearation @ 100% duty cycle assuming the electric energy price as 1,37 TL/kWh

time = 10*365*24; % (10yr)*(365day/yr)*(24h/yr)
energy = (P/1e3)*time; % convert to KW back
format shortEng
price = 1.37 * energy
price = 
    18.0018e+006
format default

The other losses I ignored are:

Core losses
Friction losses
Windage losses
Stray Losses

function [ph1, ph2, ph3, t] = getThreePhiSinusoid(freq, timeEnd, vllrms)
    omg = 2*pi*freq;
    t = linspace(0, timeEnd, 1e5);
    theta = omg*t;
    vll = vllrms*sqrt(2);
    vln = vll/sqrt(3);
    ph1 = vln*sin(theta);
    ph2 = vln*sin(theta - 2*pi/3);
    ph3 = vln*sin(theta + 2*pi/3);
end

function [sig t] = getSigFromPhasor(phasor, timeEnd)
    freq = 50;
    omg = 2*pi*freq;
    t = linspace(0, timeEnd, 1e5);
    theta = omg*t;
    sig = abs(phasor)*sin(theta - angle(phasor));
end

function strToReturn = getPolar(var)
    precision = 3;
    r = abs(var); theta_deg = rad2deg(angle(var));
    if (abs(r) < 0.1 || abs(r) > 1e3) % display in exponential form
        rstr = num2str(r, "%." + num2str(precision, '%d') + "e"); 
    else
        rstr = num2str(r, "%." + num2str(precision, '%d') + "f");
    end
    if (abs(theta_deg) < 0.1) % display in exponential form
        thetaStr = num2str(theta_deg, "%." + num2str(precision, '%d') + "e");
    else
        thetaStr = num2str(theta_deg, "%." + num2str(precision, '%d') + "f");
    end
    strToReturn = rstr + "∠" + thetaStr + "°";
end

function arrToReturn = getPolarVal(var)
    r = abs(var); theta_deg = rad2deg(angle(var));
    arrToReturn = [r theta_deg];
end

function [valToReturn index] = findClosest(array, value)
    [garbage index] = min(abs(array - value));
    valToReturn = array(index);
end

function valToReturn = immersiveFinder(xArray, yArray, xValue)
    xValue = double(xValue);
    [garbageX xIndex] = findClosest(xArray, xValue);
    valToReturn = yArray(xIndex);
end

https://www.mathworks.com/matlabcentral/fileexchange/70244-polar-phasors

% edited from:
% https://www.mathworks.com/matlabcentral/fileexchange/70244-polar-phasors
function [ pas ] = spplot(varargin)
    %Plotter phasors
    figure1=figure('Position', [0, 0, 1024, 1200]);
    subplot(2,1,1)
    msg =  'Fejl';
    maxval = 0;
    minval = 0;
    clrs = {[0 0.4470 0.7410], 	[0.8500 0.3250 0.0980], [0.9290 0.6940 0.1250], [0.4940 0.1840 0.5560],[0.4660 0.6740 0.1880],[0.3010 0.7450 0.9330],  	[0.6350 0.0780 0.1840]};
    first = 0;
    
    if nargin == 1
        sum = varargin{1};
        Vinkel = sum(2);
        Scale = 1;
        Laengde = sum(1);
        vinkel = @(x) exp(1i*x*pi/180);
        dummyX=Laengde*Scale;
        plot=dummyX*vinkel(Vinkel);
        compass(plot);
        max_lim = abs(Laengde)*Scale ;
        x_fake=[0 max_lim 0 -max_lim];
        y_fake=[max_lim 0 -max_lim 0];
        h_fake=compass(x_fake,y_fake);
        hold on
        p1=compass(plot);
        set(h_fake,'Visible','off');
        set(findall(gcf, 'String', '30', '-or','String','60') ,'String', '  ');
        set(findall(gcf, 'String', '0'),'String', ' ' );
        set(findall(gcf, 'String', '330'),'String', ' ');
        set(findall(gcf, 'String', '300'),'String', ' ' );
        set(findall(gcf, 'String', '270'),'String', ' ');
        set(findall(gcf, 'String', '240'),'String', ' ' );
        set(findall(gcf, 'String', '210'),'String', ' ');
        set(findall(gcf, 'String', '180'),'String', ' ');
        set(findall(gcf, 'String', '150'),'String', ' ' );
        set(findall(gcf, 'String', '120'),'String', ' ');
        set(findall(gcf, 'String', '90'),'String', ' ' );
        
        for m=0:0.5:abs(Laengde)*2
            set(findall(gcf, 'String', ['  ' num2str(m)]),'String', '  ');
        end
        p1(1).LineWidth = 2;
        p1(1).Color='k';
        if 150 > Vinkel && Vinkel >= 0
            text(real(plot)+2,imag(plot)+2,[ num2str(abs(Laengde)) ' {\angle} ' num2str(Vinkel) '^{\circ}' ])
        elseif 360 > Vinkel && Vinkel >= 210
            text(real(plot)-2,imag(plot)-2,[ num2str(abs(Laengde)) ' {\angle} ' num2str(Vinkel) '^{\circ}' ])
        elseif 210 > Vinkel && Vinkel >= 150
            text(real(plot)-2,imag(plot)+2,[ num2str(abs(Laengde)) ' {\angle} ' num2str(Vinkel) '^{\circ}' ])
        elseif Vinkel < 0
            text(real(plot)-2,imag(plot)-2,[ num2str(abs(Laengde)) ' {\angle} ' num2str(Vinkel) '^{\circ}' ])
        else
            msg = 'Error occurred.';
            error(msg)
        end
        
    elseif nargin >= 2
        for m=1:1:nargin
            M = varargin{m};
            Vinkel = M(2);
            Laengde = M(1);
            Scale = 1;
            enhed = ' ';
            pre=inputname(m);
            if first == 0
                for ii=1:1:nargin
                    O = varargin{ii};
                    if abs(O(1)) > maxval
                        maxval = abs(O(1));
                    end
                end
                max_lim = abs(maxval) ;
                x_fake=[0 max_lim 0 -max_lim];
                y_fake=[max_lim 0 -max_lim 0];
                h_fake=compass(x_fake,y_fake);
                hold on
                first = 1;
            end
            vinkel = @(x) exp(1i*x*pi/180);
            dummyX=(Laengde*Scale);
            plot=dummyX*vinkel(Vinkel);
            p1=compass(plot);
            set(h_fake,'Visible','off');
            p1(1).LineWidth = 2;
            p1(1).Color=clrs{m};
            text(real(plot)*1.15,imag(plot)*1.15,[pre])
        end
        hold off
    end
    hold off
    
    firstthis = 0;
    subplot(2,1,2)
    for k=1:1:nargin
        str=inputname(k);
        Scale = 1;
        enhed = ' ';
        pre = strcat(str, ': ');  
        BB = varargin{k};
        Vinkel1 = BB(2);
        Laengde1 = BB(1);
        xt = @(t) Laengde1.*sin(2*pi*t-(deg2rad(Vinkel1)));
        p2=fplot(xt,[0 2]);
        string=strcat(num2str(Laengde1,'%100.2f'),' \angle ',num2str(Vinkel1,'%100.2f'),'^{\circ}');
        string1=strcat(pre, 32 ,num2str(abs(Laengde1),'%100.2f'),  32 , enhed ,  32 ,' {\angle} ', num2str(Vinkel1,'%100.2f'), '^{\circ}' );
        grid on
        ax = gca;
        ax.XTick =  0:0.25:2;
        ax.XTickLabel = {'0^\circ','90^\circ','180^\circ','270^\circ','360^\circ','450^\circ','540^\circ','630^\circ','720^\circ'};
        p2(1).LineWidth = 2;
        p2(1).Color=clrs{k};
        p2(1).DisplayName=string1;
        hold on
    end
    hold off
    legend show
end