% Thruster \ Simulation
% Matlab code
% Stefan Deucker 10-10-96
% Thomas Hoover 10-11-2006 typed in from hard copy
% Thomas Hoover 11-03-2006 modified to allow for voltage and torque range,
% output to mesh graphs

clear all;
clear figure;
close all;
! del diary
diary on


%*************************************************************************
% ////////////////////////// INPUT  //////////////////////////////////////
%*************************************************************************

%*************************************************************************
% MOTOR AND OPERATING CHARACTERISTICS
%*************************************************************************

% E: applied voltage [V]
E = 0.5: .5 : 24;

% r: stator resistance [ohms]
R = 3.03;

%RT : overall resistanc [ohms]
RT = R + 1;

%KT; torque constan [N*m/A]
KT = 0.2200;

%KE: back emf constant [V/(rad/s)]
KE = 2408.6;

%TF: friction torque [Nm}
TF = 0.033;

% Tmax: maximum winding temperature [degree C]
Tmax = 130;

% TPR: thermal impedance [deg/W]
TPR = 3.13;

%TL: external torque load range [N cm] {n m / 100 = n cm}
TL = 0.1 : 0.1 : 25;


%T: ambient Temperature [degree C]
T = 3;

% D1: rotor diameter [mm]
D1 = 39.37;

% D2: stator diameter [mm]
D2 = 40.37;

% L: length of rotor/stator [mm]
L = 31.75;

% G: gearbox speed reduction ratio
G = 20;

%S: simulate shaft seal? T/F 1/0
S = 1;

%*************************************************************************
% OIL CHARACTERISTICS
%*************************************************************************

% is the motor oil filled? T/F 1/0
oil = 1;

%T1: first temperature of known viscosity 
T1 = 40;

% visc1: viscosity at T1 [Ns/m^2]
visc1 = 0.001569;

%T2: second temperature of known viscosity 
T2 = 20;

% visc2: viscosity at T1 [Ns/m^2]
visc2 = 0.002484;

% viscosity of air, no oil filling [Ns/m^2]
viscair = 0.00001798;

%*************************************************************************
% ////////////////////////// CALCULATE    ////////////////////////////////
%*************************************************************************

%*************************************************************************
% calc. torqu loss due to gearbox + seal for each load torque
%*************************************************************************

TLs = TL*G; % gearbox shaft torque
effgear = G^(-0.0385); % estimated gearbox eff
TMgear = ((TLs/(effgear*G))-(TLs/G))/100; %torque loss due to gearbox
if S==1
    TMgear = TMgear + 0.012; %torque loss due to shaft seal
end

%*************************************************************************
% calc. characteristic motor curves
%************************************************************************

for m=1:size(E,2)
    for n=1:size(TL,2)    
        Tactual(m,n) = T ;
    end
end

%figure

% if no heating up  of the thruster unit is to be simulated, set the loop to
% 1 and the first plot will show the thruster unit performance at room
% temperature

 for a = 1:4
    for m=1:size(E,2)
     for n=1:size(TL,2) % defines max n to be number of elements in TL

        %calculate fill flued viscousity, oil or air
        viscosity(m,n) = visc(viscair, oil, visc2, visc1, T2, T1, Tactual (m,n));
            if viscosity(m,n) < viscair
                viscosity(m,n) = viscair;
            end
        % viscous damping constant, need for speed calculation
        DF(m,n) = (-pi * (D1 * 0.001)^3 * (L*0.001)* viscosity(m,n))/((D1*0.001)-(D2*0.001));

        %calculate speed
        omega(m,n) = -( -E(m)*KT + RT*(TL(n)/100) + RT*(TF+TMgear(n)))/ (DF(m,n)*RT + KT*KE);
        rpm(m,n) = omega(m,n)*30/pi;
        if rpm(m,n) < 0
               rpm(m,n) = 0;
            end
        screw(m,n) = rpm(m,n) / G;
        if screw(m,n) < 0
                screw(m,n) = 0;
            end
        
        % calucuate torqu loss, current, power out, power in and eff.
        %frictional and viscous losses
        TM(m,n) = DF(n) * omega(m,n) + (TF+TMgear(n));
        if TM(m,n) < 0
                TM(m,n) = 0;
        end
        %current
        I(m,n) = (1/KT) * (TL(n)/100) + TM(m,n)/KT;
        if I(m,n) < 0
                I(m,n) = 0;
        end
        % power out
        Pout(m,n) = 1/30 * pi * rpm(m,n) * TL(n)/ 100;
        if Pout(m,n) < 0
                Pout(m,n) = 0;
        end
        %power in
        Pin(m,n) = E(m) * I(m,n);
        if Pin(m,n) < 0
                Pin(m,n) = 0;
        end
        eff(m,n) = Pout(m,n) / Pin(m,n);
        if eff(m,n) < 0
                eff(m,n) = 0;
        end
      end % end n loop
    end % end m loop
    
    %*************************************************************************
    % calc. overall ultimate temp. rise due to power losses
    %************************************************************************
    % electric power losses
    IR = I .* I .* R;
    % viscous and frictional power losses
    loss = TM .* omega;
    % temperature rise due to electrical losses
    deltaTIR = IR .* TPR;
    % temparture rise due to visc and fric losses
    deltaTloss = loss .* TPR;
    % ultimate temp rise
    deltaT = deltaTIR + deltaTloss;
    %actual winding temperature
    Tactual = T + deltaT;
    
    %*************************************************************************
    %  convergence check plot
    %************************************************************************
    
    hold on
    
    %subplot(1,2,1), 
    figure(1)
    meshc(TL, E, viscosity)
    %plot (TL,viscosity,'.y')
    %xlabel ('load torque [N cm]')
    %ylabel ('viscosity [Ns/m2]')
    %title('Convergence Check ')
    
    %subplot(1,2,2), 
    figure(2) %2
    meshc(TL, E, Tactual)
    %plot (TL, Tactual, '-y')
    %xlabel ('load torque [N cm]')
    %ylabel ('Tactual [C]')
    %ylabel ('Temperature [C]')
    
   
    

 end
% determine for which load torque the max winding temp is exceeded
over = 0;
for m=1:size(E,2)
    for n=1:size(TL,2)  
        if Tactual(m,n) < Tmax
            countm = m;
            countn = n;
        end;  
        
    end
end

%*************************************************************************
%  calculate altered motor constants and resulting performance 
%************************************************************************
     
% constant for rare earth magnets, ELCOM series
C = -0.00025;
clear m
clear n
for m=1:size(E,2)
     for n=1:size(TL,2) % defines max n to be number of elements in TL
         %viscosity
    viscosity(m,n) = visc(viscair, oil, visc2, visc1, T2, T1, Tactual(m,n));
    if viscosity(m,n) < viscair
        viscosity(m,n) = viscair;
    end
    %viscous damping constant 
    DF(m,n) = (-pi * (D1*0.001)^3 * (L* 0.001) * viscosity(m,n)) / ...
        ((D1 * 0.001) - (D2 *0.001));
    % motor constants at each operating point
    RTsteady(m,n) = (RT-R) + (R * (234.5 + (23 + deltaT(m,n))) / (234.5 + 23));
    KTsteady(m,n) = KT * (1 + C * (deltaT(m,n)));
    KEsteady(m,n) = KE * (1 + C * (deltaT(m,n)));
    
    % calculate speed
    omegasteady(m,n) = -(-E(m) * KTsteady(m,n) + RTsteady(m,n)*(TL(n)/100) + ... 
        RTsteady(m,n) * (TF + TMgear(n))) / (DF(m,n) * RTsteady(m,n) + ... 
        KTsteady(m,n) * KEsteady(m,n));

    rpmsteady(m,n) = omegasteady(m,n) * 30 / pi;
    if rpmsteady(m,n) < 0
                rpmsteady(m,n) = 0;
    end
    screwsteady(m,n) = rpmsteady(m,n) / G;
    if screwsteady(m,n) < 0
                screwsteady(m,n) = 0;
        end
    %calculate torque loss, current power out, power in and eff.
    TMsteady(m,n) = DF(m,n) * omegasteady(m,n) + (TF + TMgear(n));
    if TMsteady(m,n) < 0
                TMsteady(m,n) = 0;
        end
    Isteady(m,n) = (1/KTsteady(m,n)) * ( TL(n) / 100) + TMsteady(m,n)/KTsteady(m,n);
    if Isteady(m,n) < 0
                Isteady(m,n) = 0;
        end
    Poutsteady(m,n) = 1/30 * pi .* rpmsteady(m,n) * TL(n)/100;
    if Poutsteady(m,n) < 0
                Poutsteady(m,n) = 0;
        end
    Pinsteady(m,n) = E(m) * Isteady(m,n);
    if Pinsteady(m,n) < 0
                Pinsteady(m,n) = 0;
        end
    effsteady(m,n) = Poutsteady(m,n) / Pinsteady(m,n);
    if effsteady(m,n) < 0
                effsteady(m,n) = 0;
        end
        
end
end
%*************************************************************************
%  ////////////////////// WRITE /////////////////////////////////////
%************************************************************************

para1 = ['Voltage:' num2str(E) ' Amb.temp.:' num2str(T) ' Therm.impedance:' num2str(TPR)];
para2 = ['Oil:' num2str(oil) ' Gearbox:' num2str(G) '  Seal:' num2str(S)  '  '];
para3 = ['The max. efficiency without temp. effects is:' num2str(max(eff)) '  '];
para4 = ['The max. efficiency with    temp. effects is:' num2str(max(effsteady)) '  '];

%t = ['The maximum operating voltage (due to the max. winding temp) is ' num2str(E(countm)) ' volts']; 
%u = ['The maximum operating load torque (due to the max. winding temp) is ' num2str(TL(countn)) ' Ncm'];    
%v = ['    TL       - Poutst   - effst   - Ist     - rpmst   - T      - Pout   - eff'];
%w = [ TL' Poutsteady' effsteady' Isteady' rpmsteady' Tactual' Pout' eff'];
disp(para1)
disp(para2)
disp(para3)
disp(para4)
%disp(t)
%disp(u)
%disp(v)
%disp(w)

diary off

%*************************************************************************
%  ////////////////////// PLOT /////////////////////////////////////
%************************************************************************
figure(1)
% plot final value in convergence plot
hold on
%subplot(1,2,1), 
meshc(TL, E, viscosity)
%plot(TL, viscosity, '-r')
xlabel('Load Torque [N cm]')
ylabel('Applied voltage [V]')
zlabel('Viscosity [Ns/m2]')
title('Oil Viscosity Convergence Check')
colorbar
%axis([ 0 TL(countn) 0 0.006])
%legend('Convg Visc 1','Convg Visc 2','Convg Visc 3','Convg Visc 4','Convg Visc 5','Convg Visc 6','Oil Viscosity')


figure(2) %2
hold on
%subplot(1,2,2), 
meshc(TL, E, Tactual)
hold on
title('Actual Winding Temperature [C] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Actual Winding Temperature [C]')
colorbar
%plot (TL,Tactual,'.r')
%xlabel('Load Torque [N cm]')
%ylabel('Actual Winding Temperature [C]')
%title('Actual Winding Temperature vs Load Torque')
%axis([ 0 TL(countn) 0 Tactual(countm, countn) ])

%subplot(2,2,3), 
figure(2) %5
meshc(TL, E, deltaTIR)
hold on
title('Temp.Rise due to I2R and Friction Losses [C] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Temp.Rise due to I2R and Friction Losses [C]')
colorbar
%plot(TL, deltaTIR,'-g')
%hold on
%plot(TL,deltaTloss, '-.g')
%xlabel('Load Torque [N cm]')
%ylabel('Winding Temperature Rise [C]')
%title('Temp.Rise due to I2R and Friction Losses vs. Load Torque [Ncm]')
%axis([ TL(1) TL(countn) 0 deltaTIR(countm, countn)])
% plot power losses and temparture rise

%subplot(2,2,4), 
figure(2) %6
meshc(TL, E, Tactual)
title('Ultim. Temp. [C] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Ultimate Temp [C]')
colorbar
%axis([ TL(1) TL(countn) 0 Tactual(countm, countn)])
%legend('Convg Temp 1','Convg Temp 2','Convg Temp 3','Convg Temp 4','Convg Temp 5','Convg Temp 6','Temp Actual Winding','Temp I2R+','Temp Frict ','UltTemp')

figure(3)
hold on
%subplot(2,2,1), 
meshc(TL, E, loss)
%title('Power [W] vs. Load Torque [Ncm] & Applied Voltage [V]')
%xlabel('Load Torque [N cm]')
%ylabel('Applied Voltage [V]')
%zlabel('Power [W]')
%plot(TL, loss, '-.r')
%hold on
%xlabel('Load Torque [N cm]')
%ylabel('Power [W]')
%title('Visc. + Fric. P - Los [W] vs. Load Torque [Ncm]')
%axis([ TL(1) TL(countn) 0 max(loss)])

%subplot(2,2,2), 
figure(3) %4
meshc(TL, E, IR)
%title('Power [W] vs. Load Torque [Ncm] & Applied Voltage [V]')
%xlabel('Load Torque [N cm]')
%ylabel('Applied Voltage [V]')
%zlabel('Power [W]')

%plot(TL, IR,':r')
%hold on
%xlabel('Load Torque [N cm]')
%ylabel('Power [W]')
%title('I2R Losses [W] vs. Load Torque [Ncm]')
%axis([ TL(1) TL(count) 0 IR(count)])

%subplot(2,2,3), 
figure(3) %9
meshc(TL, E, Pout)
meshc(TL, E, Pin)
meshc(TL, E, Poutsteady)
meshc(TL, E, Pinsteady)
title('Power [W] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Power [W]')
colorbar

%plot (TL, Pout,'-g')
%hold on
%plot(TL,Pin,'--b')
%hold on
%plot(TL,Poutsteady,'-.g')
%hold on
%plot(TL,Pinsteady,':b')
%axis([TL(1) TL(count) 0 Pinsteady(count) ])
%title('Power [W] vs. Load Torque [Ncm]')
%legend('Visc. + Fric. Loss','I2R Losses','Pout' ,'Pin','PoutSteady','PinSteady')

% plot comparison between curves, altered results
%figure
%subplot(2,2,1), 
figure(7)
hold on
meshc(TL, E, I)
meshc(TL, E, Isteady)
title('Current [A] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Current [A]')
colorbar
%plot (TL,I,'-g')
%hold on
%plot(TL,Isteady,'--g')
%axis([TL(1) TL(count) 0 Isteady(count) ])
%title('Current [A] vs. Load Torque [Ncm]')
%legend('Isteady','Isteady')

%subplot(2,2,2), 
figure(8) %8
hold on
meshc(TL, E, rpm)
meshc(TL, E, rpmsteady)
meshc(TL, E, screw)
meshc(TL, E, screwsteady)
title('Shaft& Screw [RPM] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Motor Shaft and Screw [rpm]')
colorbar
%plot (TL,rpm,'-r')
%hold on
%plot(TL,rpmsteady,'-.r')
%hold on
%plot(TL,screw,'-b')
%hold on
%plot(TL,screwsteady,'-.b')
%axis([TL(1) TL(count) 0 rpm(1) ])
%title('Speed [rpm] vs. Load Torque [Ncm]')
%legend('MotorRPM','MotorRPMsteady','ScrewRPM','ScrewRPMsteady')
%legend('ScrewRPM','ScrewRPMsteady')


%subplot(2,2,4), 
figure(9) %10
hold on
meshc(TL, E, eff)
meshc(TL, E, effsteady)
%meshc(TL, E, screw)
%meshc(TL, E, screwsteady)
title('Efficiency [percent] vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Effciency [percent]')
colorbar

figure(10)
hold on
meshc(TL,E,effsteady*max(max(screwsteady)))
surfc(TL,E,screwsteady)
title('Screw Speed & Efficiency vs. Load Torque [Ncm] & Applied Voltage [V]')
xlabel('Load Torque [N cm]')
ylabel('Applied Voltage [V]')
zlabel('Screw Speed [rpm] and Efficiency [percent of max RPM]')
colorbar

%plot (TL, eff,'-m')
%hold on
%plot(TL,effsteady,'-.m')
%axis([TL(1) TL(count) 0 1 ])
%title('Efficiency vs. Load Torque [Ncm]')
%legend('Eff','Effsteady')