% 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 = 1: 2 : 60;

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

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

%KT; torque constant [N*m/A]
KT = 0.5800;

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

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

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

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

%TL: external torque load range [N cm] {n m / 100 = n cm}
TL = .1 : .1 : 75;


%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 = 1;

%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;

%*************************************************************************
% Screw and Vehicle Characteristics
%*************************************************************************

% Screw values based on current Dorado class vehicle
% Vehicle values based on JGB research

ScrewAdvance = 0.25; % meters of advance per revolution
RPMatFullSpeed = 60 / ScrewAdvance; % RPM needed for speed of 1 Meter/sec
VehicleSpeed = 0.1:.1:1.0; %speed range for vehicle in meters/sec
VehicleMultiplier = 1.26; % JGB calc
VolumeBase = 0.33; % JGB calc
CDV = 0.08; % JGB calc
PropEff = 0.5; %JGB calc

for t = 1:size(VehicleSpeed,2)
    PropulsionPower(t)=((1/PropEff)* 0.5 * CDV * (VolumeBase*VehicleMultiplier)^2 * ... 
    1027 * VehicleSpeed(t)^3);
    RadPerSec(t) = (RPMatFullSpeed * 2 * 3.1415926 * VehicleSpeed(t)/60);
    TorqueSpeed(t) = PropulsionPower(t)/RadPerSec(t);
end
    



%*************************************************************************
% ////////////////////////// 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);
        if omega(m,n) < 0
               omega(m,n) = 0;
            end
        rpm(m,n) = omega(m,n)*30/pi;
        screw(m,n) = rpm(m,n) / G;
        
        % 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));
        
        %current
        I(m,n) = (1/KT) * (TL(n)/100) + TM(m,n)/KT;
        
        % power out
        Pout(m,n) = 1/30 * pi * rpm(m,n) * TL(n)/ 100;
        
        %power in
        Pin(m,n) = E(m) * I(m,n);
        
        eff(m,n) = Pout(m,n) / Pin(m,n);
        
      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
clear t
clear ESpeed
zero = 0
ESpeed = zeros(size(VehicleSpeed))
for n=1:size(TL,2) % defines max n to be number of elements in TL
    for m=1:size(E,2)
            %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));
    
        if omegasteady(m,n) < 0
                omegasteady(m,n) = 0;
        end

        rpmsteady(m,n) = omegasteady(m,n) * 30 / pi;
    
        screwsteady(m,n) = rpmsteady(m,n) / G;
    
        for t = 1:size(VehicleSpeed,2)
            %if n < size(TL,2)
                if TL(n) < TorqueSpeed(t)
                    if screwsteady(m,n) < (RPMatFullSpeed * VehicleSpeed(t))
                        ESpeed(t) = E(m); 
                    end
                 end
            %end
        end
    
    
        %calculate torque loss, current power out, power in and eff.
        TMsteady(m,n) = DF(m,n) * omegasteady(m,n) + (TF + TMgear(n));
    
        Isteady(m,n) = (1/KTsteady(m,n)) * ( TL(n) / 100) + TMsteady(m,n)/KTsteady(m,n);
    
        Poutsteady(m,n) = 1/30 * pi .* rpmsteady(m,n) * TL(n)/100;
    
        Pinsteady(m,n) = E(m) * Isteady(m,n);
    
        effsteady(m,n) = Poutsteady(m,n) / Pinsteady(m,n);
    
        speed(m,n) = screwsteady(m,n) .* ScrewAdvance ./ 100;    
    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(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), 
contour3(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
contour3(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
contour3(TL, E, Pout)
contour3(TL, E, Pin)
contour3(TL, E, Poutsteady)
contour3(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
contour3(TL, E, I)
contour3(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)
contour3(TL, E, screw)
contour3(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
contour3(TL, E, eff)
contour3(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
contour3(TL,E,effsteady * max(max(screwsteady)))
contour3(TL,E,screwsteady)
contour(TorqueSpeed,ESpeed)
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')