From df53d90f572ba80a55a253437c14128033bde062 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Sun, 26 Mar 2017 19:27:34 -0400 Subject: [PATCH 01/21] overhaul data import --- Matlab/sample.asv | 19 +++++ Matlab/sample.ini | 17 +++++ Matlab/sample.tdu | 7 ++ Matlab/tdu-air.asv | 10 +++ Matlab/tdu.asv | 152 +++++++++++++++++++++++++++++++++++++++ Matlab/tdu.m | 161 ++++++++++++++++++++---------------------- Matlab/tduOBSOLETE.m | 164 +++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 444 insertions(+), 86 deletions(-) create mode 100644 Matlab/sample.asv create mode 100644 Matlab/sample.ini create mode 100644 Matlab/sample.tdu create mode 100644 Matlab/tdu-air.asv create mode 100644 Matlab/tdu.asv create mode 100644 Matlab/tduOBSOLETE.m diff --git a/Matlab/sample.asv b/Matlab/sample.asv new file mode 100644 index 0000000..bcc7de1 --- /dev/null +++ b/Matlab/sample.asv @@ -0,0 +1,19 @@ +[Prop] +Name = Air +k = 1.4 +MW = .0289645 +MW_unit = kg/mol +[Chamber] +T_0 = 273 +T_unit = K +P_0 = 101325 +P_unit = Pa +[Nozzle] +Inlet_radius= .0075 +Exit_radius = .00708 +Throat_radius = .005 +Radius_unit =m +15 deg + + + diff --git a/Matlab/sample.ini b/Matlab/sample.ini new file mode 100644 index 0000000..709da93 --- /dev/null +++ b/Matlab/sample.ini @@ -0,0 +1,17 @@ +[Prop] +Name = Air +k = 1.4 +MW = .0289645 +MW_unit = kg/mol +[Chamber] +T_0 = 273 +T_unit = K +P_0 = 101325 +P_unit = Pa +[Nozzle] +Inlet_radius = .0075 +Exit_radius = .00708 +Throat_radius = .005 +Radius_unit = m +Alpha = 15 +Alpha_unit = deg \ No newline at end of file diff --git a/Matlab/sample.tdu b/Matlab/sample.tdu new file mode 100644 index 0000000..70f46ad --- /dev/null +++ b/Matlab/sample.tdu @@ -0,0 +1,7 @@ +sample +Air 1.4 .0289645 kg/mol +273 K +101325 Pa +.0075 .00708 .005 m +15 deg + diff --git a/Matlab/tdu-air.asv b/Matlab/tdu-air.asv new file mode 100644 index 0000000..872e05c --- /dev/null +++ b/Matlab/tdu-air.asv @@ -0,0 +1,10 @@ +[Prop] +Name=Air +k=1.4 .0289645 kg/mol +273 K +101325 Pa +.0075 .00708 .005 m +15 deg + + + diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv new file mode 100644 index 0000000..bcf8df9 --- /dev/null +++ b/Matlab/tdu.asv @@ -0,0 +1,152 @@ +function varargout = tdu % main function + clear all; + clc; + format long + % universal constants + R_0 = 8.3144598; % [J/(mol*K)] universal gas constant + g_0 = 9.81; % [m/s^2] standard gravity + unitless='[-]'; + global debug; + debug = true; +% === THRUSTER PARAMETERS === + +% % MANUAL ENTRY +% % % gas properties of propellant +% % prop_name = 'Air'; +% % k = 1.4; % 1.4 for air +% % mw = .0289645; % .0289645 for air +% % mw_units = '[kg/mol]'; +% % +% % % chamber conditions +% % T_0 = 273; % stagnation temperature +% % temperature_units = '[K]'; +% % P_0 = 101325; % stagnation pressure +% % pressure_units = '[Pa]'; +% % +% % % nozzle geometry +% % inlet_radius = .0075; % radius at inlet of converging section +% % exit_radius = .00708; % radius at exit of diverging section +% % throat_radius = .005; % radius at throat +% % % exit_radius = .00708; % radius at nozzle exit +% % % throat_radius = .005; % radius at nozzle throat +% % length_units = '[m]'; +% % conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal +% % angle_units = '[deg]'; +% ===------------=== + +% Math (the fun part) + +% NOZZLE NOMENCLATURE +% ******************************** +% +% /- +% / +% ===\---/ c = chamber +% (c) (t) (e) t = throat +% ===/---\ e = exit +% \ +% \- +% ******************************** + + R = R_0/mw; % gas constant + R_units = '[J/kg K]'; + A_e = pi*exit_radius^2; % exit area + A_t = pi*throat_radius^2; % throat area + area_units = '[m^2]'; + length = (exit_radius-throat_radius)/tan(deg2rad(conical_half_angle)); % [m] + + % Throat conditions + + + % format & display outputs + linedivider='------------'; + result = {'Propellant','',prop_name; + linedivider,'',''; + 'Specific heat ratio', k, unitless; + 'Molar mass', mw, mw_units; + 'Specific gas constant',R,R_units; + linedivider,'',''; + 'Chamber temperature', T_0, temperature_units; + 'Chamber pressure', P_0, pressure_units; + 'Exit radius', exit_radius, length_units; + 'Throat radius', throat_radius, length_units; + 'Exit area',A_e,'[m^2]'; + 'Throat area',A_t,'[m^2]'; + 'Half-angle',conical_half_angle, angle_units; + linedivider,'',''; + 'Length',length,length_units; + 'Exit area',A_e,area_units; + 'Throat area',A_t,area_units; + 'Throat temperature',throat_temperature,temperature_units; + 'Throat pressure',throat_pressure,pressure_units; + 'Mass flow rate',mass_flowrate,mass_flowrate_units; + linedivider,'',''; + 'Exit temperature',exit_temperature,temperature_units; + 'Exit pressure',exit_pressure,pressure_units; + linedivider,'',''; + 'Exhaust velocity',exit_velocity,velocity_units; + 'Thrust',thrust,force_units; + 'Specific impulse',specific_impulse,isp_units; + linedivider,'',''; + 'Exit Mach',exit_mach,unitless; + 'A/At',A_e/A_t,unitless; + 'T/Tc',exit_temperature/T_0,unitless; + 'P/Pc',exit_pressure/P_0,unitless; + 'v/at',exit_velocity/throat_velocity,unitless; + }; + display(result); +end + +function Mach = solve_mach(A,At,k) +% solve Mach number from area ratio by Newton-Raphson Method. (assume +% supersonic) +% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html + P = 2/(k+1); + Q = 1-P; + R = (A/At).^((2*Q)/P); + a = Q.^(1/P); + r = (R-1)/(2*a); + X = 1/((1+r)+sqrt(r*(r+2))); % initial guess + diff = 1; % initalize termination criteria + while abs(diff) > .0001 + F = (P*X+Q).^(1/P)-R*X; + dF = (P*X+Q).^((1/P)-1)-R; + Xnew = X - F/dF; + diff = Xnew - X; + X = Xnew; + end + Mach = 1/sqrt(X); +end +function loadfromfile(filein) +% read data from input file + ini = fopen(filein,'r'); + %%fprintf('reading data from input file %s...\n',filein); + header = fscanf(ini,'%s',[1 1]); % descriptive header (no quotes, no spaces) + prop_params = fscanf(ini,'%s%g',[1 4]); % scan propellant parameters + prop_name = prop_params(1,1); % name of propellant + k = prop_params(1,2); % specific heat ratio + mw = prop_params(1,3); % molecular weight + mw_units + + g = fscanf(ini,'%g',[1 1]); % body force acceleration (gravity) + nodes=zeros(numnod,4); + for i=1:numnod + nodes(i,:) = fscanf(ini,'%g',[1 4]); % node xcoord bcflag bcvalue + end + %%fprintf('\tnodes imported.\n'); + elems=zeros(numel,6); + for i=1:numel + elems(i,:) = fscanf(ini,'%g',[1 6]); % elem node1 node2 dia E dens + end + %%fprintf('\telements imported.\n'); + fclose('all'); %close input file + %fprintf('\tinput file closed.\n'); +end +function display(result) + [n,~]=size(result); + for i = 1:n + fprintf('\n%24s\t%15.8f\t%s',result{i,:}); + end + fprintf('\n') +end +[ k mw mw_units] = \ No newline at end of file diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 45fe8e3..9225c53 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -2,35 +2,39 @@ clear all; clc; format long - %% universal constants - universal_gas_constant= 8.3144598; % [J/(mol*K)] - standard_gravity = 9.81; % [m/s^2] + % universal constants + R_0 = 8.3144598; % [J/(mol*K)] universal gas constant + g_0 = 9.81; % [m/s^2] standard gravity unitless='[-]'; - -%% === EDIT THESE === - % gas properties of propellant - propellant_name = 'Air'; - specific_heat_ratio = 1.4; % 1.4 for air - molecular_weight = .0289645; % .0289645 for air - molar_mass_units = '[kg/mol]'; - - % chamber conditions - chamber_temperature = 273; % 273 K = 0C - temperature_units = '[K]'; - chamber_pressure = 101325; % 1 atm = 101325 Pa - pressure_units = '[Pa]'; - - % nozzle geometry - exit_radius = .00708; % radius at nozzle exit - throat_radius = .005; % radius at nozzle throat -% exit_radius = .00708; % radius at nozzle exit -% throat_radius = .005; % radius at nozzle throat - length_units = '[m]'; - conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal - angle_units = '[deg]'; + global debug; + debug = true; +% === THRUSTER PARAMETERS === + +% % MANUAL ENTRY +% % % gas properties of propellant +% % prop_name = 'Air'; +% % k = 1.4; % 1.4 for air +% % mw = .0289645; % .0289645 for air +% % mw_units = '[kg/mol]'; +% % +% % % chamber conditions +% % T_0 = 273; % stagnation temperature +% % temperature_units = '[K]'; +% % P_0 = 101325; % stagnation pressure +% % pressure_units = '[Pa]'; +% % +% % % nozzle geometry +% % inlet_radius = .0075; % radius at inlet of converging section +% % exit_radius = .00708; % radius at exit of diverging section +% % throat_radius = .005; % radius at throat +% % % exit_radius = .00708; % radius at nozzle exit +% % % throat_radius = .005; % radius at nozzle throat +% % length_units = '[m]'; +% % conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal +% % angle_units = '[deg]'; % ===------------=== -%% Math (the fun part) +% Math (the fun part) % NOZZLE NOMENCLATURE % ******************************** @@ -44,74 +48,35 @@ % \- % ******************************** - specific_gas_constant = universal_gas_constant/molecular_weight; - specific_gas_constant_units = '[J/kg K]'; - exit_area = pi*exit_radius^2; - throat_area = pi*throat_radius^2; + R = R_0/mw; % gas constant + R_units = '[J/kg K]'; + A_e = pi*exit_radius^2; % exit area + A_t = pi*throat_radius^2; % throat area area_units = '[m^2]'; length = (exit_radius-throat_radius)/tan(deg2rad(conical_half_angle)); % [m] - %% Throat conditions - % We are designing a nozzle that "chokes" the fluid flow at the throat, - % so we assume that at the throat, Mach number is isentropically - % brought to the boundary between supersonic and subsonic flow (Mach=1) - throat_temperature = chamber_temperature*(2/(specific_heat_ratio+1)); - throat_pressure = chamber_pressure*(2/(specific_heat_ratio+1))^(specific_heat_ratio/(specific_heat_ratio-1)); - throat_velocity = sqrt(specific_heat_ratio*specific_gas_constant*throat_temperature); - % use nozzle areas to find Mach number at the exit via the Newton-Raphson method - exit_mach = solve_mach(exit_area,throat_area,specific_heat_ratio); -% since -% mass_flow_rate = density*speed*cross_sectional_area -% and from the ideal gas law, PV=RT - throat_density = throat_pressure/(specific_gas_constant*throat_temperature); %[kg/m^3] -% and at the throat Mach = 1 so flow rate (speed of flow) is the fluid's speed of sound at the throat - throat_flowrate = sqrt(specific_heat_ratio*specific_gas_constant*throat_temperature); - mass_flowrate = throat_density*throat_flowrate*throat_area; - mass_flowrate_units = '[kg/s]'; + % Throat conditions - %% Exit conditions -% it's much easier to define the flow characteristics at the exit in -% terms of ratios for now. We'll resolve them into actual useful values -% later. -% these are ratios between the throat parameter (numerator) and the exit -% parameter (denominator). -% For example, area_ratio = throat_area/exit_area - -%(already know this one) area_ratio = exit_mach(((specific_heat_ratio+1)/2)/(1+(specific_heat_ratio-1)/2*exit_mach^2))^((specific_heat_ratio+1)/(2*(specific_heat_ratio-1))); - temperature_ratio = 1+((specific_heat_ratio-1)/2)*exit_mach^2; - pressure_ratio = temperature_ratio^(specific_heat_ratio/(specific_heat_ratio-1)); - - exit_temperature = chamber_temperature/temperature_ratio; - exit_pressure = chamber_pressure/pressure_ratio; - exit_velocity = sqrt((2*specific_heat_ratio*specific_gas_constant*chamber_temperature)/(specific_heat_ratio-1)*(1-1/(1+(specific_heat_ratio-1)/2*exit_mach^2))); - velocity_units = '[m/s]'; - - %% Thrust and Specific Impulse - thrust = mass_flowrate*exit_velocity + exit_pressure*exit_area; - force_units = '[N]'; - specific_impulse = thrust/(mass_flowrate*standard_gravity); - isp_units='[s]'; - - %% format & display outputs + % format & display outputs linedivider='------------'; - result = {'Propellant','',propellant_name; + result = {'Propellant','',prop_name; linedivider,'',''; - 'Specific heat ratio', specific_heat_ratio, unitless; - 'Molar mass', molecular_weight, molar_mass_units; - 'Specific gas constant',specific_gas_constant,specific_gas_constant_units; + 'Specific heat ratio', k, unitless; + 'Molar mass', mw, mw_units; + 'Specific gas constant',R,R_units; linedivider,'',''; - 'Chamber temperature', chamber_temperature, temperature_units; - 'Chamber pressure', chamber_pressure, pressure_units; + 'Chamber temperature', T_0, temperature_units; + 'Chamber pressure', P_0, pressure_units; 'Exit radius', exit_radius, length_units; 'Throat radius', throat_radius, length_units; - 'Exit area',exit_area,'[m^2]'; - 'Throat area',throat_area,'[m^2]'; + 'Exit area',A_e,'[m^2]'; + 'Throat area',A_t,'[m^2]'; 'Half-angle',conical_half_angle, angle_units; linedivider,'',''; 'Length',length,length_units; - 'Exit area',exit_area,area_units; - 'Throat area',throat_area,area_units; + 'Exit area',A_e,area_units; + 'Throat area',A_t,area_units; 'Throat temperature',throat_temperature,temperature_units; 'Throat pressure',throat_pressure,pressure_units; 'Mass flow rate',mass_flowrate,mass_flowrate_units; @@ -124,12 +89,11 @@ 'Specific impulse',specific_impulse,isp_units; linedivider,'',''; 'Exit Mach',exit_mach,unitless; - 'A/At',exit_area/throat_area,unitless; - 'T/Tc',exit_temperature/chamber_temperature,unitless; - 'P/Pc',exit_pressure/chamber_pressure,unitless; + 'A/At',A_e/A_t,unitless; + 'T/Tc',exit_temperature/T_0,unitless; + 'P/Pc',exit_pressure/P_0,unitless; 'v/at',exit_velocity/throat_velocity,unitless; }; - display(result); end @@ -153,7 +117,31 @@ end Mach = 1/sqrt(X); end - +function loadfromfile(filein) +% read data from input file + fid = fopen(filein,'r'); + %%fprintf('reading data from input file %s...\n',filein); + header = fscanf(fid,'%*s',[1 1]) % descriptive header (no quotes, no spaces) + prop_params = fscanf(fid,'%s\t%g\t%g\t%s',[1 4]) % scan propellant parameters + prop_name = prop_params(1,1) % name of propellant + k = prop_params(1,2) % specific heat ratio + mw = prop_params(1,3) % molecular weight +% mw_units +% +% g = fscanf(ini,'%g',[1 1]); % body force acceleration (gravity) +% nodes=zeros(numnod,4); +% for i=1:numnod +% nodes(i,:) = fscanf(ini,'%g',[1 4]); % node xcoord bcflag bcvalue +% end +% %%fprintf('\tnodes imported.\n'); +% elems=zeros(numel,6); +% for i=1:numel +% elems(i,:) = fscanf(ini,'%g',[1 6]); % elem node1 node2 dia E dens +% end + %%fprintf('\telements imported.\n'); + fclose('all'); %close input file + %fprintf('\tinput file closed.\n'); +end function display(result) [n,~]=size(result); for i = 1:n @@ -161,3 +149,4 @@ function display(result) end fprintf('\n') end +[ k mw mw_units] = \ No newline at end of file diff --git a/Matlab/tduOBSOLETE.m b/Matlab/tduOBSOLETE.m new file mode 100644 index 0000000..ba9fa92 --- /dev/null +++ b/Matlab/tduOBSOLETE.m @@ -0,0 +1,164 @@ +function varargout = tdu % main function + clear all; + clc; + format long + %% universal constants + R_0= 8.3144598; % [J/(mol*K)] universal gas constant + g_0 = 9.81; % [m/s^2] standard gravity + unitless='[-]'; + +%% === EDIT THESE === + % gas properties of propellant + propellant_name = 'Air'; + k = 1.4; % 1.4 for air + molecular_weight = .0289645; % .0289645 for air + molar_mass_units = '[kg/mol]'; + + % chamber conditions + T_0 = 273; % stagnation temperature + temperature_units = '[K]'; + P_0 = 101325; % stagnation pressure + pressure_units = '[Pa]'; + + % nozzle geometry + inlet_radius = .0075; % radius at inlet of converging section + exit_radius = .00708; % radius at exit of diverging section + throat_radius = .005; % radius at throat +% exit_radius = .00708; % radius at nozzle exit +% throat_radius = .005; % radius at nozzle throat + length_units = '[m]'; + conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal + angle_units = '[deg]'; +% ===------------=== + +%% Math (the fun part) + +% NOZZLE NOMENCLATURE +% ******************************** +% +% /- +% / +% ===\---/ c = chamber +% (c) (t) (e) t = throat +% ===/---\ e = exit +% \ +% \- +% ******************************** + + R = R_0/molecular_weight; % gas constant + R_units = '[J/kg K]'; + A_e = pi*exit_radius^2; % exit area + A_t = pi*throat_radius^2; % throat area + area_units = '[m^2]'; + length = (exit_radius-throat_radius)/tan(deg2rad(conical_half_angle)); % [m] + + %% Throat conditions + % We are designing a nozzle that "chokes" the fluid flow at the throat, + % so we assume that at the throat, Mach number is isentropically + % brought to the boundary between supersonic and subsonic flow (Mach=1) + throat_temperature = T_0*(2/(k+1)); + throat_pressure = P_0*(2/(k+1))^(k/(k-1)); + throat_velocity = sqrt(k*R*throat_temperature); + % use nozzle areas to find Mach number at the exit via the Newton-Raphson method + exit_mach = solve_mach(A_e,A_t,k); +% since +% mass_flow_rate = density*speed*cross_sectional_area +% and from the ideal gas law, PV=RT + throat_density = throat_pressure/(R*throat_temperature); %[kg/m^3] +% and at the throat Mach = 1 so flow rate (speed of flow) is the fluid's speed of sound at the throat + throat_flowrate = sqrt(k*R*throat_temperature); + mass_flowrate = throat_density*throat_flowrate*A_t; + mass_flowrate_units = '[kg/s]'; + + %% Exit conditions +% it's much easier to define the flow characteristics at the exit in +% terms of ratios for now. We'll resolve them into actual useful values +% later. +% these are ratios between the throat parameter (numerator) and the exit +% parameter (denominator). +% For example, area_ratio = A_t/A_e + +%(already know this one) area_ratio = exit_mach(((k+1)/2)/(1+(k-1)/2*exit_mach^2))^((k+1)/(2*(k-1))); + temperature_ratio = 1+((k-1)/2)*exit_mach^2; + pressure_ratio = temperature_ratio^(k/(k-1)); + + exit_temperature = T_0/temperature_ratio; + exit_pressure = P_0/pressure_ratio; + + exit_velocity = sqrt((2*k*R*T_0)/(k-1)*(1-1/(1+(k-1)/2*exit_mach^2))); + velocity_units = '[m/s]'; + + %% Thrust and Specific Impulse + thrust = mass_flowrate*exit_velocity + exit_pressure*A_e; + force_units = '[N]'; + specific_impulse = thrust/(mass_flowrate*g_0); + isp_units='[s]'; + + %% format & display outputs + linedivider='------------'; + result = {'Propellant','',propellant_name; + linedivider,'',''; + 'Specific heat ratio', k, unitless; + 'Molar mass', molecular_weight, molar_mass_units; + 'Specific gas constant',R,R_units; + linedivider,'',''; + 'Chamber temperature', T_0, temperature_units; + 'Chamber pressure', P_0, pressure_units; + 'Exit radius', exit_radius, length_units; + 'Throat radius', throat_radius, length_units; + 'Exit area',A_e,'[m^2]'; + 'Throat area',A_t,'[m^2]'; + 'Half-angle',conical_half_angle, angle_units; + linedivider,'',''; + 'Length',length,length_units; + 'Exit area',A_e,area_units; + 'Throat area',A_t,area_units; + 'Throat temperature',throat_temperature,temperature_units; + 'Throat pressure',throat_pressure,pressure_units; + 'Mass flow rate',mass_flowrate,mass_flowrate_units; + linedivider,'',''; + 'Exit temperature',exit_temperature,temperature_units; + 'Exit pressure',exit_pressure,pressure_units; + linedivider,'',''; + 'Exhaust velocity',exit_velocity,velocity_units; + 'Thrust',thrust,force_units; + 'Specific impulse',specific_impulse,isp_units; + linedivider,'',''; + 'Exit Mach',exit_mach,unitless; + 'A/At',A_e/A_t,unitless; + 'T/Tc',exit_temperature/T_0,unitless; + 'P/Pc',exit_pressure/P_0,unitless; + 'v/at',exit_velocity/throat_velocity,unitless; + }; + + display(result); +end + +function Mach = solve_mach(A,At,k) +% solve Mach number from area ratio by Newton-Raphson Method. (assume +% supersonic) +% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html + P = 2/(k+1); + Q = 1-P; + R = (A/At).^((2*Q)/P); + a = Q.^(1/P); + r = (R-1)/(2*a); + X = 1/((1+r)+sqrt(r*(r+2))); % initial guess + diff = 1; % initalize termination criteria + while abs(diff) > .0001 + F = (P*X+Q).^(1/P)-R*X; + dF = (P*X+Q).^((1/P)-1)-R; + Xnew = X - F/dF; + diff = Xnew - X; + X = Xnew; + end + Mach = 1/sqrt(X); +end + +function display(result) + [n,~]=size(result); + for i = 1:n + fprintf('\n%24s\t%15.8f\t%s',result{i,:}); + end + fprintf('\n') +end From ba49ab95f39aff0d92925e147e0214bbef692bd1 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Mon, 27 Mar 2017 07:49:21 -0400 Subject: [PATCH 02/21] import data from .tdu file --- Matlab/sample.tdu | 11 ++++---- Matlab/tdu.asv | 63 ++++++++++++++++++++++---------------------- Matlab/tdu.m | 66 +++++++++++++++++++++++++---------------------- 3 files changed, 72 insertions(+), 68 deletions(-) diff --git a/Matlab/sample.tdu b/Matlab/sample.tdu index 70f46ad..6406f0a 100644 --- a/Matlab/sample.tdu +++ b/Matlab/sample.tdu @@ -1,7 +1,6 @@ -sample -Air 1.4 .0289645 kg/mol -273 K -101325 Pa -.0075 .00708 .005 m -15 deg +Air +1.4 .0289645 +273 101325 +.0075 .005 .00708 +15 diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv index bcf8df9..e0cbbb6 100644 --- a/Matlab/tdu.asv +++ b/Matlab/tdu.asv @@ -2,6 +2,9 @@ function varargout = tdu % main function clear all; clc; format long + + filein = 'sample.tdu'; % input file + % universal constants R_0 = 8.3144598; % [J/(mol*K)] universal gas constant g_0 = 9.81; % [m/s^2] standard gravity @@ -9,8 +12,31 @@ function varargout = tdu % main function global debug; debug = true; % === THRUSTER PARAMETERS === - +% read data from input file +fid = fopen(filein,'r'); +if debug;fprintf('reading data from input file %s...\n',filein);end; +prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) + if debug;fprintf('Propellant:\t%s',prop_name);end; +prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters + k = prop_params(1,1); % specific heat ratio + mw = prop_params(1,2); % molecular weight + if debug;fprintf('k:\t%g\nmw:',k);end; +total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters + T_0 = total_params(1,1); % total temperature + P_0 = total_params(1,2); % total pressure +geom = fscanf(fid,'%g',[1 3]); % scan nozzle geometry + inlet_radius= geom(1,1); % radius at inlet of converging section + exit_radius = geom(1,2); % radius at exit of diverging section + throat_radius= geom(1,3); % radius at throat +alpha = fscanf(fid,'%g',[1,1]); % conical half angle +fclose('all'); %close input file +if debug;fprintf('\tinput file closed.\n');end; % % MANUAL ENTRY +mw_units = '[kg/mol]'; +temperature_units = '[K]'; +pressure_units = '[Pa]'; +length_units = '[m]'; +angle_units = '[deg]'; % % % gas properties of propellant % % prop_name = 'Air'; % % k = 1.4; % 1.4 for air @@ -30,7 +56,7 @@ function varargout = tdu % main function % % % exit_radius = .00708; % radius at nozzle exit % % % throat_radius = .005; % radius at nozzle throat % % length_units = '[m]'; -% % conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal +% % alpha = 15; % half angle of conical nozzle, 15 degrees is optimal % % angle_units = '[deg]'; % ===------------=== @@ -53,7 +79,7 @@ function varargout = tdu % main function A_e = pi*exit_radius^2; % exit area A_t = pi*throat_radius^2; % throat area area_units = '[m^2]'; - length = (exit_radius-throat_radius)/tan(deg2rad(conical_half_angle)); % [m] + length = (exit_radius-throat_radius)/tan(deg2rad(alpha)); % [m] % Throat conditions @@ -72,7 +98,7 @@ function varargout = tdu % main function 'Throat radius', throat_radius, length_units; 'Exit area',A_e,'[m^2]'; 'Throat area',A_t,'[m^2]'; - 'Half-angle',conical_half_angle, angle_units; + 'Half-angle',alpha, angle_units; linedivider,'',''; 'Length',length,length_units; 'Exit area',A_e,area_units; @@ -117,36 +143,11 @@ function Mach = solve_mach(A,At,k) end Mach = 1/sqrt(X); end -function loadfromfile(filein) -% read data from input file - ini = fopen(filein,'r'); - %%fprintf('reading data from input file %s...\n',filein); - header = fscanf(ini,'%s',[1 1]); % descriptive header (no quotes, no spaces) - prop_params = fscanf(ini,'%s%g',[1 4]); % scan propellant parameters - prop_name = prop_params(1,1); % name of propellant - k = prop_params(1,2); % specific heat ratio - mw = prop_params(1,3); % molecular weight - mw_units - - g = fscanf(ini,'%g',[1 1]); % body force acceleration (gravity) - nodes=zeros(numnod,4); - for i=1:numnod - nodes(i,:) = fscanf(ini,'%g',[1 4]); % node xcoord bcflag bcvalue - end - %%fprintf('\tnodes imported.\n'); - elems=zeros(numel,6); - for i=1:numel - elems(i,:) = fscanf(ini,'%g',[1 6]); % elem node1 node2 dia E dens - end - %%fprintf('\telements imported.\n'); - fclose('all'); %close input file - %fprintf('\tinput file closed.\n'); -end + function display(result) [n,~]=size(result); for i = 1:n fprintf('\n%24s\t%15.8f\t%s',result{i,:}); end fprintf('\n') -end -[ k mw mw_units] = \ No newline at end of file +end \ No newline at end of file diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 9225c53..2fdd21a 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -2,6 +2,9 @@ clear all; clc; format long + + filein = 'sample.tdu'; % input file + % universal constants R_0 = 8.3144598; % [J/(mol*K)] universal gas constant g_0 = 9.81; % [m/s^2] standard gravity @@ -9,8 +12,34 @@ global debug; debug = true; % === THRUSTER PARAMETERS === - +% read data from input file +fid = fopen(filein,'r'); +if debug;fprintf('reading data from input file (%s)...\n',filein);end; +prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) + if debug;fprintf('\tPropellant:\t%s\n',prop_name);end; +prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters + k = prop_params(1,1); % specific heat ratio + mw = prop_params(1,2); % molecular weight + if debug;fprintf('\tk:\t%g\n\tmw:\t%g\n',k,mw);end; +total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters + T_0 = total_params(1,1); % total temperature + P_0 = total_params(1,2); % total pressure + if debug;fprintf('\tT_0:\t%g\n\tP_0:%g\n',T_0,P_0);end; +geom = fscanf(fid,'%g',[1 3]); % scan nozzle geometry + inlet_radius= geom(1,1); % radius at inlet of converging section + throat_radius= geom(1,2); % radius at throat + exit_radius = geom(1,3); % radius at exit of diverging section + if debug;fprintf('\tinlet radius:\t%g\n\tthroat radius:\t%g\n\texit radius:\t%g\n',inlet_radius,throat_radius,exit_radius);end; +alpha = fscanf(fid,'%g',[1,1]); % conical half angle + if debug;fprintf('\talpha:\t%g\n',alpha);end; +fclose('all'); %close input file +if debug;fprintf('input file closed.\n');end; % % MANUAL ENTRY +mw_units = '[kg/mol]'; +temperature_units = '[K]'; +pressure_units = '[Pa]'; +length_units = '[m]'; +angle_units = '[deg]'; % % % gas properties of propellant % % prop_name = 'Air'; % % k = 1.4; % 1.4 for air @@ -30,7 +59,7 @@ % % % exit_radius = .00708; % radius at nozzle exit % % % throat_radius = .005; % radius at nozzle throat % % length_units = '[m]'; -% % conical_half_angle = 15; % half angle of conical nozzle, 15 degrees is optimal +% % alpha = 15; % half angle of conical nozzle, 15 degrees is optimal % % angle_units = '[deg]'; % ===------------=== @@ -53,7 +82,7 @@ A_e = pi*exit_radius^2; % exit area A_t = pi*throat_radius^2; % throat area area_units = '[m^2]'; - length = (exit_radius-throat_radius)/tan(deg2rad(conical_half_angle)); % [m] + length = (exit_radius-throat_radius)/tan(deg2rad(alpha)); % [m] % Throat conditions @@ -72,7 +101,7 @@ 'Throat radius', throat_radius, length_units; 'Exit area',A_e,'[m^2]'; 'Throat area',A_t,'[m^2]'; - 'Half-angle',conical_half_angle, angle_units; + 'Half-angle',alpha, angle_units; linedivider,'',''; 'Length',length,length_units; 'Exit area',A_e,area_units; @@ -117,36 +146,11 @@ end Mach = 1/sqrt(X); end -function loadfromfile(filein) -% read data from input file - fid = fopen(filein,'r'); - %%fprintf('reading data from input file %s...\n',filein); - header = fscanf(fid,'%*s',[1 1]) % descriptive header (no quotes, no spaces) - prop_params = fscanf(fid,'%s\t%g\t%g\t%s',[1 4]) % scan propellant parameters - prop_name = prop_params(1,1) % name of propellant - k = prop_params(1,2) % specific heat ratio - mw = prop_params(1,3) % molecular weight -% mw_units -% -% g = fscanf(ini,'%g',[1 1]); % body force acceleration (gravity) -% nodes=zeros(numnod,4); -% for i=1:numnod -% nodes(i,:) = fscanf(ini,'%g',[1 4]); % node xcoord bcflag bcvalue -% end -% %%fprintf('\tnodes imported.\n'); -% elems=zeros(numel,6); -% for i=1:numel -% elems(i,:) = fscanf(ini,'%g',[1 6]); % elem node1 node2 dia E dens -% end - %%fprintf('\telements imported.\n'); - fclose('all'); %close input file - %fprintf('\tinput file closed.\n'); -end + function display(result) [n,~]=size(result); for i = 1:n fprintf('\n%24s\t%15.8f\t%s',result{i,:}); end fprintf('\n') -end -[ k mw mw_units] = \ No newline at end of file +end \ No newline at end of file From ee5f6372509ea98feb3b07be599afdc473bd1007 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Mon, 27 Mar 2017 07:50:10 -0400 Subject: [PATCH 03/21] cleanup --- .gitignore | 1 + Matlab/sample.asv | 19 ------ Matlab/tdu-air.asv | 10 --- Matlab/tdu.asv | 153 --------------------------------------------- 4 files changed, 1 insertion(+), 182 deletions(-) delete mode 100644 Matlab/sample.asv delete mode 100644 Matlab/tdu-air.asv delete mode 100644 Matlab/tdu.asv diff --git a/.gitignore b/.gitignore index ed05ea3..e1d5b53 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,4 @@ dist/ .doctrees/ *.egg-info/ __pycache__/ +.asv/ diff --git a/Matlab/sample.asv b/Matlab/sample.asv deleted file mode 100644 index bcc7de1..0000000 --- a/Matlab/sample.asv +++ /dev/null @@ -1,19 +0,0 @@ -[Prop] -Name = Air -k = 1.4 -MW = .0289645 -MW_unit = kg/mol -[Chamber] -T_0 = 273 -T_unit = K -P_0 = 101325 -P_unit = Pa -[Nozzle] -Inlet_radius= .0075 -Exit_radius = .00708 -Throat_radius = .005 -Radius_unit =m -15 deg - - - diff --git a/Matlab/tdu-air.asv b/Matlab/tdu-air.asv deleted file mode 100644 index 872e05c..0000000 --- a/Matlab/tdu-air.asv +++ /dev/null @@ -1,10 +0,0 @@ -[Prop] -Name=Air -k=1.4 .0289645 kg/mol -273 K -101325 Pa -.0075 .00708 .005 m -15 deg - - - diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv deleted file mode 100644 index e0cbbb6..0000000 --- a/Matlab/tdu.asv +++ /dev/null @@ -1,153 +0,0 @@ -function varargout = tdu % main function - clear all; - clc; - format long - - filein = 'sample.tdu'; % input file - - % universal constants - R_0 = 8.3144598; % [J/(mol*K)] universal gas constant - g_0 = 9.81; % [m/s^2] standard gravity - unitless='[-]'; - global debug; - debug = true; -% === THRUSTER PARAMETERS === -% read data from input file -fid = fopen(filein,'r'); -if debug;fprintf('reading data from input file %s...\n',filein);end; -prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) - if debug;fprintf('Propellant:\t%s',prop_name);end; -prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters - k = prop_params(1,1); % specific heat ratio - mw = prop_params(1,2); % molecular weight - if debug;fprintf('k:\t%g\nmw:',k);end; -total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters - T_0 = total_params(1,1); % total temperature - P_0 = total_params(1,2); % total pressure -geom = fscanf(fid,'%g',[1 3]); % scan nozzle geometry - inlet_radius= geom(1,1); % radius at inlet of converging section - exit_radius = geom(1,2); % radius at exit of diverging section - throat_radius= geom(1,3); % radius at throat -alpha = fscanf(fid,'%g',[1,1]); % conical half angle -fclose('all'); %close input file -if debug;fprintf('\tinput file closed.\n');end; -% % MANUAL ENTRY -mw_units = '[kg/mol]'; -temperature_units = '[K]'; -pressure_units = '[Pa]'; -length_units = '[m]'; -angle_units = '[deg]'; -% % % gas properties of propellant -% % prop_name = 'Air'; -% % k = 1.4; % 1.4 for air -% % mw = .0289645; % .0289645 for air -% % mw_units = '[kg/mol]'; -% % -% % % chamber conditions -% % T_0 = 273; % stagnation temperature -% % temperature_units = '[K]'; -% % P_0 = 101325; % stagnation pressure -% % pressure_units = '[Pa]'; -% % -% % % nozzle geometry -% % inlet_radius = .0075; % radius at inlet of converging section -% % exit_radius = .00708; % radius at exit of diverging section -% % throat_radius = .005; % radius at throat -% % % exit_radius = .00708; % radius at nozzle exit -% % % throat_radius = .005; % radius at nozzle throat -% % length_units = '[m]'; -% % alpha = 15; % half angle of conical nozzle, 15 degrees is optimal -% % angle_units = '[deg]'; -% ===------------=== - -% Math (the fun part) - -% NOZZLE NOMENCLATURE -% ******************************** -% -% /- -% / -% ===\---/ c = chamber -% (c) (t) (e) t = throat -% ===/---\ e = exit -% \ -% \- -% ******************************** - - R = R_0/mw; % gas constant - R_units = '[J/kg K]'; - A_e = pi*exit_radius^2; % exit area - A_t = pi*throat_radius^2; % throat area - area_units = '[m^2]'; - length = (exit_radius-throat_radius)/tan(deg2rad(alpha)); % [m] - - % Throat conditions - - - % format & display outputs - linedivider='------------'; - result = {'Propellant','',prop_name; - linedivider,'',''; - 'Specific heat ratio', k, unitless; - 'Molar mass', mw, mw_units; - 'Specific gas constant',R,R_units; - linedivider,'',''; - 'Chamber temperature', T_0, temperature_units; - 'Chamber pressure', P_0, pressure_units; - 'Exit radius', exit_radius, length_units; - 'Throat radius', throat_radius, length_units; - 'Exit area',A_e,'[m^2]'; - 'Throat area',A_t,'[m^2]'; - 'Half-angle',alpha, angle_units; - linedivider,'',''; - 'Length',length,length_units; - 'Exit area',A_e,area_units; - 'Throat area',A_t,area_units; - 'Throat temperature',throat_temperature,temperature_units; - 'Throat pressure',throat_pressure,pressure_units; - 'Mass flow rate',mass_flowrate,mass_flowrate_units; - linedivider,'',''; - 'Exit temperature',exit_temperature,temperature_units; - 'Exit pressure',exit_pressure,pressure_units; - linedivider,'',''; - 'Exhaust velocity',exit_velocity,velocity_units; - 'Thrust',thrust,force_units; - 'Specific impulse',specific_impulse,isp_units; - linedivider,'',''; - 'Exit Mach',exit_mach,unitless; - 'A/At',A_e/A_t,unitless; - 'T/Tc',exit_temperature/T_0,unitless; - 'P/Pc',exit_pressure/P_0,unitless; - 'v/at',exit_velocity/throat_velocity,unitless; - }; - display(result); -end - -function Mach = solve_mach(A,At,k) -% solve Mach number from area ratio by Newton-Raphson Method. (assume -% supersonic) -% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html - P = 2/(k+1); - Q = 1-P; - R = (A/At).^((2*Q)/P); - a = Q.^(1/P); - r = (R-1)/(2*a); - X = 1/((1+r)+sqrt(r*(r+2))); % initial guess - diff = 1; % initalize termination criteria - while abs(diff) > .0001 - F = (P*X+Q).^(1/P)-R*X; - dF = (P*X+Q).^((1/P)-1)-R; - Xnew = X - F/dF; - diff = Xnew - X; - X = Xnew; - end - Mach = 1/sqrt(X); -end - -function display(result) - [n,~]=size(result); - for i = 1:n - fprintf('\n%24s\t%15.8f\t%s',result{i,:}); - end - fprintf('\n') -end \ No newline at end of file From 2d441e4a9500201cd0d792ea83d92da347f08199 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Mon, 27 Mar 2017 21:14:58 -0400 Subject: [PATCH 04/21] reformat input file and data import --- Matlab/sample.tdu | 16 ++++++- Matlab/tdu.m | 106 ++++++++++++++++++++-------------------------- 2 files changed, 60 insertions(+), 62 deletions(-) diff --git a/Matlab/sample.tdu b/Matlab/sample.tdu index 6406f0a..d1d79e0 100644 --- a/Matlab/sample.tdu +++ b/Matlab/sample.tdu @@ -1,6 +1,18 @@ Air 1.4 .0289645 273 101325 -.0075 .005 .00708 -15 +11 +0.00 0.07 +0.01 0.06 +0.02 0.05 +0.03 0.05258819 +0.04 0.055176381 +0.05 0.057764571 +0.06 0.060352762 +0.07 0.062940952 +0.08 0.065529143 +0.09 0.068117333 +0.10 0.070705524 + + diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 2fdd21a..d279ae7 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -12,80 +12,66 @@ global debug; debug = true; % === THRUSTER PARAMETERS === -% read data from input file -fid = fopen(filein,'r'); -if debug;fprintf('reading data from input file (%s)...\n',filein);end; -prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) - if debug;fprintf('\tPropellant:\t%s\n',prop_name);end; -prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters - k = prop_params(1,1); % specific heat ratio - mw = prop_params(1,2); % molecular weight - if debug;fprintf('\tk:\t%g\n\tmw:\t%g\n',k,mw);end; -total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters - T_0 = total_params(1,1); % total temperature - P_0 = total_params(1,2); % total pressure - if debug;fprintf('\tT_0:\t%g\n\tP_0:%g\n',T_0,P_0);end; -geom = fscanf(fid,'%g',[1 3]); % scan nozzle geometry - inlet_radius= geom(1,1); % radius at inlet of converging section - throat_radius= geom(1,2); % radius at throat - exit_radius = geom(1,3); % radius at exit of diverging section - if debug;fprintf('\tinlet radius:\t%g\n\tthroat radius:\t%g\n\texit radius:\t%g\n',inlet_radius,throat_radius,exit_radius);end; -alpha = fscanf(fid,'%g',[1,1]); % conical half angle - if debug;fprintf('\talpha:\t%g\n',alpha);end; -fclose('all'); %close input file -if debug;fprintf('input file closed.\n');end; -% % MANUAL ENTRY -mw_units = '[kg/mol]'; -temperature_units = '[K]'; -pressure_units = '[Pa]'; -length_units = '[m]'; -angle_units = '[deg]'; -% % % gas properties of propellant -% % prop_name = 'Air'; -% % k = 1.4; % 1.4 for air -% % mw = .0289645; % .0289645 for air -% % mw_units = '[kg/mol]'; -% % -% % % chamber conditions -% % T_0 = 273; % stagnation temperature -% % temperature_units = '[K]'; -% % P_0 = 101325; % stagnation pressure -% % pressure_units = '[Pa]'; -% % -% % % nozzle geometry -% % inlet_radius = .0075; % radius at inlet of converging section -% % exit_radius = .00708; % radius at exit of diverging section -% % throat_radius = .005; % radius at throat -% % % exit_radius = .00708; % radius at nozzle exit -% % % throat_radius = .005; % radius at nozzle throat -% % length_units = '[m]'; -% % alpha = 15; % half angle of conical nozzle, 15 degrees is optimal -% % angle_units = '[deg]'; + % IMPORT FROM FILE + fid = fopen(filein,'r'); + if debug;fprintf('reading data from input file (%s)...\n',filein);end + prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) + if debug;fprintf('\tPropellant:\t\t%8s\n',prop_name);end + prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters + k = prop_params(1,1); % specific heat ratio + mw = prop_params(1,2); % molecular weight + if debug;fprintf('\tk:\t\t%16g\n\tmw:\t\t%16g\n',k,mw);end + total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters + T_0 = total_params(1,1); % total temperature + P_0 = total_params(1,2); % total pressure + if debug;fprintf('\tT_0:\t%16g\n\tP_0:\t%16g\n',T_0,P_0);end + geom_size = fscanf(fid,'%g',[1 1]); % number of geometry nodes + xcoord = zeros(geom_size,1); radius = zeros(geom_size,1); + for i=1:geom_size + geom = fscanf(fid,'%g',[1 2]); + xcoord(i) = geom(1,1); % x coordinate of geometry node + radius(i) = geom(1,2); % radius at xcoord + end + if debug + fprintf('\tinlet radius:\t%8f\n\tthroat radius:\t%8f\n\texit radius:\t%8f\n',radius(1),min(radius),radius(end)); + fprintf('\tlength:\t%16f\n\tgeometry nodes:\t%8i\n',xcoord(end),geom_size); + end + fclose('all'); %close input file + if debug;fprintf('input file closed.\n');end + % MANUAL ENTRY + mw_units = '[kg/mol]'; + temperature_units = '[K]'; + pressure_units = '[Pa]'; + length_units = '[m]'; + angle_units = '[deg]'; % ===------------=== -% Math (the fun part) - % NOZZLE NOMENCLATURE % ******************************** % % /- -% / -% ===\---/ c = chamber -% (c) (t) (e) t = throat -% ===/---\ e = exit +% / 0 = total parameters +% ===\---/ 1 = chamber +% (1) (2) (3) 2 = throat +% ===/---\ 3 = exit % \ % \- % ******************************** +% % ASSUMPTIONS +% - Isentropic flow +% - Initial temperature and pressure are total parameters + R = R_0/mw; % gas constant R_units = '[J/kg K]'; A_e = pi*exit_radius^2; % exit area A_t = pi*throat_radius^2; % throat area area_units = '[m^2]'; - length = (exit_radius-throat_radius)/tan(deg2rad(alpha)); % [m] - % Throat conditions - + for index = 2:geom_size + + end + % format & display outputs linedivider='------------'; @@ -126,13 +112,13 @@ display(result); end -function Mach = solve_mach(A,At,k) +function Mach = arearatio2mach(A,Astar,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume % supersonic) % https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html P = 2/(k+1); Q = 1-P; - R = (A/At).^((2*Q)/P); + R = (A/Astar).^((2*Q)/P); a = Q.^(1/P); r = (R-1)/(2*a); X = 1/((1+r)+sqrt(r*(r+2))); % initial guess From d5ed2c9b3955466fc596a0084d2be9b00a225cff Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 31 Mar 2017 00:10:15 -0400 Subject: [PATCH 05/21] estimate super and subsonic, accurate msd nodes --- Matlab/msd.tdu | 155 ++++++++++++++++++++++++++++++++++++++++++++++ Matlab/sample.tdu | 97 +++++++++++++++++++++++++++-- Matlab/tdu.m | 142 ++++++++++++++++++++++++++++++++---------- 3 files changed, 357 insertions(+), 37 deletions(-) create mode 100644 Matlab/msd.tdu diff --git a/Matlab/msd.tdu b/Matlab/msd.tdu new file mode 100644 index 0000000..c3157bb --- /dev/null +++ b/Matlab/msd.tdu @@ -0,0 +1,155 @@ +Nitrogen +1.4 .0280134 +300 23922.146 +151 +0 0.009525 +0.01 0.009254657 +0.02 0.008984314 +0.03 0.008713971 +0.04 0.008443629 +0.05 0.008173286 +0.06 0.007902943 +0.07 0.0076326 +0.08 0.007362257 +0.09 0.007091914 +0.1 0.006821571 +0.11 0.006551229 +0.12 0.006280886 +0.13 0.006010543 +0.14 0.0057402 +0.15 0.005469857 +0.16 0.005199514 +0.17 0.004929171 +0.18 0.004658829 +0.19 0.004388486 +0.2 0.004118143 +0.21 0.0038478 +0.22 0.003577457 +0.23 0.003307114 +0.24 0.003036771 +0.25 0.002766429 +0.26 0.002496086 +0.27 0.002225743 +0.28 0.0019554 +0.29 0.001685057 +0.3 0.001414714 +0.31 0.001144371 +0.32 0.000874029 +0.33 0.000603686 +0.34 0.000333343 +0.35 0.000063 +0.36 0.000063 +0.37 0.000127754 +0.38 0.000192509 +0.39 0.000257263 +0.4 0.000322018 +0.41 0.000386772 +0.42 0.000451526 +0.43 0.000516281 +0.44 0.000581035 +0.45 0.000645789 +0.46 0.000710544 +0.47 0.000775298 +0.48 0.000840053 +0.49 0.000904807 +0.5 0.000969561 +0.51 0.001034316 +0.52 0.00109907 +0.53 0.001163825 +0.54 0.001228579 +0.55 0.001293333 +0.56 0.001358088 +0.57 0.001422842 +0.58 0.001487596 +0.59 0.001552351 +0.6 0.001617105 +0.61 0.00168186 +0.62 0.001746614 +0.63 0.001811368 +0.64 0.001876123 +0.65 0.001940877 +0.66 0.002005632 +0.67 0.002070386 +0.68 0.00213514 +0.69 0.002199895 +0.7 0.002264649 +0.71 0.002329404 +0.72 0.002394158 +0.73 0.002458912 +0.74 0.002523667 +0.75 0.002588421 +0.76 0.002653175 +0.77 0.00271793 +0.78 0.002782684 +0.79 0.002847439 +0.8 0.002912193 +0.81 0.002976947 +0.82 0.003041702 +0.83 0.003106456 +0.84 0.003171211 +0.85 0.003235965 +0.86 0.003300719 +0.87 0.003365474 +0.88 0.003430228 +0.89 0.003494982 +0.9 0.003559737 +0.91 0.003624491 +0.92 0.003689246 +0.93 0.003754 +0.94 0.003818754 +0.95 0.003883509 +0.96 0.003948263 +0.97 0.004013018 +0.98 0.004077772 +0.99 0.004142526 +1 0.004207281 +1.01 0.004272035 +1.02 0.004336789 +1.03 0.004401544 +1.04 0.004466298 +1.05 0.004531053 +1.06 0.004595807 +1.07 0.004660561 +1.08 0.004725316 +1.09 0.00479007 +1.1 0.004854825 +1.11 0.004919579 +1.12 0.004984333 +1.13 0.005049088 +1.14 0.005113842 +1.15 0.005178596 +1.16 0.005243351 +1.17 0.005308105 +1.18 0.00537286 +1.19 0.005437614 +1.2 0.005502368 +1.21 0.005567123 +1.22 0.005631877 +1.23 0.005696632 +1.24 0.005761386 +1.25 0.00582614 +1.26 0.005890895 +1.27 0.005955649 +1.28 0.006020404 +1.29 0.006085158 +1.3 0.006149912 +1.31 0.006214667 +1.32 0.006279421 +1.33 0.006344175 +1.34 0.00640893 +1.35 0.006473684 +1.36 0.006538439 +1.37 0.006603193 +1.38 0.006667947 +1.39 0.006732702 +1.4 0.006797456 +1.41 0.006862211 +1.42 0.006926965 +1.43 0.006991719 +1.44 0.007056474 +1.45 0.007121228 +1.46 0.007185982 +1.47 0.007250737 +1.48 0.007315491 +1.49 0.007380246 +1.5 0.007445 \ No newline at end of file diff --git a/Matlab/sample.tdu b/Matlab/sample.tdu index d1d79e0..c9d9e96 100644 --- a/Matlab/sample.tdu +++ b/Matlab/sample.tdu @@ -1,18 +1,107 @@ Air 1.4 .0289645 273 101325 -11 -0.00 0.07 +101 +0 0.07 +0.001 0.069 +0.002 0.068 +0.003 0.067 +0.004 0.066 +0.005 0.065 +0.006 0.064 +0.007 0.063 +0.008 0.062 +0.009 0.061 0.01 0.06 +0.011 0.059 +0.012 0.058 +0.013 0.057 +0.014 0.056 +0.015 0.055 +0.016 0.054 +0.017 0.053 +0.018 0.052 +0.019 0.051 0.02 0.05 +0.021 0.050258819 +0.022 0.050517638 +0.023 0.050776457 +0.024 0.051035276 +0.025 0.051294095 +0.026 0.051552914 +0.027 0.051811733 +0.028 0.052070552 +0.029 0.052329371 0.03 0.05258819 +0.031 0.052847009 +0.032 0.053105829 +0.033 0.053364648 +0.034 0.053623467 +0.035 0.053882286 +0.036 0.054141105 +0.037 0.054399924 +0.038 0.054658743 +0.039 0.054917562 0.04 0.055176381 +0.041 0.0554352 +0.042 0.055694019 +0.043 0.055952838 +0.044 0.056211657 +0.045 0.056470476 +0.046 0.056729295 +0.047 0.056988114 +0.048 0.057246933 +0.049 0.057505752 0.05 0.057764571 +0.051 0.05802339 +0.052 0.058282209 +0.053 0.058541028 +0.054 0.058799848 +0.055 0.059058667 +0.056 0.059317486 +0.057 0.059576305 +0.058 0.059835124 +0.059 0.060093943 0.06 0.060352762 +0.061 0.060611581 +0.062 0.0608704 +0.063 0.061129219 +0.064 0.061388038 +0.065 0.061646857 +0.066 0.061905676 +0.067 0.062164495 +0.068 0.062423314 +0.069 0.062682133 0.07 0.062940952 +0.071 0.063199771 +0.072 0.06345859 +0.073 0.063717409 +0.074 0.063976228 +0.075 0.064235047 +0.076 0.064493867 +0.077 0.064752686 +0.078 0.065011505 +0.079 0.065270324 0.08 0.065529143 +0.081 0.065787962 +0.082 0.066046781 +0.083 0.0663056 +0.084 0.066564419 +0.085 0.066823238 +0.086 0.067082057 +0.087 0.067340876 +0.088 0.067599695 +0.089 0.067858514 0.09 0.068117333 -0.10 0.070705524 - +0.091 0.068376152 +0.092 0.068634971 +0.093 0.06889379 +0.094 0.069152609 +0.095 0.069411428 +0.096 0.069670247 +0.097 0.069929066 +0.098 0.070187886 +0.099 0.070446705 +0.1 0.070705524 diff --git a/Matlab/tdu.m b/Matlab/tdu.m index d279ae7..031f746 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -1,6 +1,6 @@ function varargout = tdu % main function - clear all; clc; + close all format long filein = 'sample.tdu'; % input file @@ -64,16 +64,72 @@ R = R_0/mw; % gas constant R_units = '[J/kg K]'; - A_e = pi*exit_radius^2; % exit area - A_t = pi*throat_radius^2; % throat area +% A_e = pi*exit_radius^2; % exit area +% A_t = pi*throat_radius^2; % throat area + A=pi.*radius.^2; area_units = '[m^2]'; - for index = 2:geom_size - - end - + A_t=min(A); + T_star=T_0*(2/(k+1)); %K + P_star=(2/(k+1))^(k/(k-1)); %Pa + rho_star=P_star/(R*T_star); % kg/m^3 + mdot=rho_star*sqrt(k*R*T_star)*A_t; %choked - % format & display outputs + M_idx=linspace(.1,4,length(xcoord)); + k1=(2/(k+1)); + k2=((k-1)/2); + k3=(.5*(k+1)/(k-1)); + area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); + + M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; + M(1)=0; + choked=false; + for x=2:length(xcoord) + if A(x)=1 + choked=true + else + disp('flag') + end + else + M_sup(x)=arearatio2mach_sup(A(x),A_t,k); + M_sub(x)=arearatio2mach_sub(A(x),A_t,k); + end + temp_ratio(x)=(1+((k-1)/2)*M(x)^2); + T(x)=T_0/temp_ratio(x); + pres_ratio(x)=temp_ratio(x)^(k/(k-1)); + P(x)=P_0/pres_ratio(x); + + temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); + T_sub(x)=T_0/temp_ratio_sub(x); + pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); + P_sub(x)=P_0/pres_ratio_sub(x); + temp_ratio_sup(x)=(1+((k-1)/2)*M_sup(x)^2); + T_sup(x)=T_0/temp_ratio_sup(x); + pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); + P_sup(x)=P_0/pres_ratio_sup(x); + end + if debug + figure + numplots=4;plotcounter=1; + subplot(numplots,1,plotcounter) + plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 radius(end) 0 xcoord(end)]); +% subplot(numplots,1,plotcounter) +% semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; + subplot(numplots,1,plotcounter) + plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1; + subplot(numplots,1,plotcounter) + plot(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1; + subplot(numplots,1,plotcounter) + plot(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1; +% figure +% semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); +% figure +% semilogy(M_idx,area_ratio,M,A./A_t,'--') + end +% % format & display outputs linedivider='------------'; result = {'Propellant','',prop_name; linedivider,'',''; @@ -81,44 +137,65 @@ 'Molar mass', mw, mw_units; 'Specific gas constant',R,R_units; linedivider,'',''; - 'Chamber temperature', T_0, temperature_units; - 'Chamber pressure', P_0, pressure_units; - 'Exit radius', exit_radius, length_units; - 'Throat radius', throat_radius, length_units; - 'Exit area',A_e,'[m^2]'; - 'Throat area',A_t,'[m^2]'; - 'Half-angle',alpha, angle_units; + 'Total temperature', T_0, temperature_units; + 'Total pressure', P_0, pressure_units; linedivider,'',''; - 'Length',length,length_units; - 'Exit area',A_e,area_units; + 'Length',xcoord(end),length_units; + 'Inlet radius',radius(1),length_units; + 'Throat radius', min(radius), length_units; + 'Exit radius', radius(end), length_units; + 'Inlet area',A(1),area_units; 'Throat area',A_t,area_units; - 'Throat temperature',throat_temperature,temperature_units; - 'Throat pressure',throat_pressure,pressure_units; - 'Mass flow rate',mass_flowrate,mass_flowrate_units; + 'Exit area',A(end),area_units; +% 'Half-angle',alpha, angle_units; + linedivider,'',''; + 'Throat temperature',T(find(A==A_t)),temperature_units; + 'Throat pressure',P(find(A==A_t)),pressure_units; +% 'Mass flow rate',mass_flowrate,mass_flowrate_units; linedivider,'',''; - 'Exit temperature',exit_temperature,temperature_units; - 'Exit pressure',exit_pressure,pressure_units; + 'Exit temperature',T(end),temperature_units; + 'Exit pressure',P(end),pressure_units; linedivider,'',''; - 'Exhaust velocity',exit_velocity,velocity_units; - 'Thrust',thrust,force_units; - 'Specific impulse',specific_impulse,isp_units; +% 'Exhaust velocity',exit_velocity,velocity_units; +% 'Thrust',thrust,force_units; +% 'Specific impulse',specific_impulse,isp_units; linedivider,'',''; - 'Exit Mach',exit_mach,unitless; - 'A/At',A_e/A_t,unitless; - 'T/Tc',exit_temperature/T_0,unitless; - 'P/Pc',exit_pressure/P_0,unitless; - 'v/at',exit_velocity/throat_velocity,unitless; + 'Exit Mach',M(end),unitless; + 'A/At',A(end)/A_t,unitless; + 'T/T0',T(end)/T_0,unitless; + 'P/P0',P(end)/P_0,unitless; +% 'v/at',exit_velocity/throat_velocity,unitless; }; display(result); end -function Mach = arearatio2mach(A,Astar,k) +function Mach = arearatio2mach_sub(A,A_t,k) +% solve Mach number from area ratio by Newton-Raphson Method. (assume +% subsonic) +% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html + P = 2/(k+1); + Q = 1-P; + R = (A/A_t).^2; + a = P.^(1/Q); + r = (R-1)/(2*a); + X = 1/((1+r)+sqrt(r*(r+2))); % initial guess + diff = 1; % initalize termination criteria + while abs(diff) > .0001 + F = (P*X+Q).^(1/P)-R*X; + dF = (P*X+Q).^((1/P)-1)-R; + Xnew = X - F/dF; + diff = Xnew - X; + X = Xnew; + end + Mach = sqrt(X); +end +function Mach = arearatio2mach_sup(A,A_t,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume % supersonic) % https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html P = 2/(k+1); Q = 1-P; - R = (A/Astar).^((2*Q)/P); + R = (A/A_t).^((2*Q)/P); a = Q.^(1/P); r = (R-1)/(2*a); X = 1/((1+r)+sqrt(r*(r+2))); % initial guess @@ -132,7 +209,6 @@ end Mach = 1/sqrt(X); end - function display(result) [n,~]=size(result); for i = 1:n From f5d6fd9525f7fb926ecc86f3ac51abfd94d37f24 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Tue, 11 Apr 2017 10:51:33 -0400 Subject: [PATCH 06/21] fix plot axes --- Matlab/tdu.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 031f746..4f6c32f 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -115,15 +115,15 @@ figure numplots=4;plotcounter=1; subplot(numplots,1,plotcounter) - plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 radius(end) 0 xcoord(end)]); + plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % subplot(numplots,1,plotcounter) % semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; subplot(numplots,1,plotcounter) - plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1; + plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); subplot(numplots,1,plotcounter) - plot(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1; + plot(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); subplot(numplots,1,plotcounter) - plot(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1; + plot(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % figure % semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); % figure From cb02f56fddb4a0587748c805ed4b6a43ca47331d Mon Sep 17 00:00:00 2001 From: runphilrun Date: Thu, 4 May 2017 18:35:09 -0400 Subject: [PATCH 07/21] update nozzle profiles --- Matlab/msd-025.tdu | 782 +++++++++++++++++++++++++++++++++++++++++ Matlab/msd-050.tdu | 855 +++++++++++++++++++++++++++++++++++++++++++++ Matlab/msd.tdu | 155 -------- Matlab/tdu.m | 11 +- 4 files changed, 1646 insertions(+), 157 deletions(-) create mode 100644 Matlab/msd-025.tdu create mode 100644 Matlab/msd-050.tdu delete mode 100644 Matlab/msd.tdu diff --git a/Matlab/msd-025.tdu b/Matlab/msd-025.tdu new file mode 100644 index 0000000..8d62771 --- /dev/null +++ b/Matlab/msd-025.tdu @@ -0,0 +1,782 @@ +Nitrogen +1.4 0.0280134 +300 239220 +778 +0.0000 0.019050 +0.0001 0.018979 +0.0002 0.018909 +0.0003 0.018838 +0.0004 0.018767 +0.0005 0.018696 +0.0006 0.018626 +0.0007 0.018555 +0.0008 0.018484 +0.0009 0.018414 +0.0010 0.018343 +0.0011 0.018272 +0.0012 0.018201 +0.0013 0.018131 +0.0014 0.018060 +0.0015 0.017989 +0.0016 0.017919 +0.0017 0.017848 +0.0018 0.017777 +0.0019 0.017706 +0.0020 0.017636 +0.0021 0.017565 +0.0022 0.017494 +0.0023 0.017424 +0.0024 0.017353 +0.0025 0.017282 +0.0026 0.017212 +0.0027 0.017141 +0.0028 0.017070 +0.0029 0.016999 +0.0030 0.016929 +0.0031 0.016858 +0.0032 0.016787 +0.0033 0.016717 +0.0034 0.016646 +0.0035 0.016575 +0.0036 0.016504 +0.0037 0.016434 +0.0038 0.016363 +0.0039 0.016292 +0.0040 0.016222 +0.0041 0.016151 +0.0042 0.016080 +0.0043 0.016009 +0.0044 0.015939 +0.0045 0.015868 +0.0046 0.015797 +0.0047 0.015727 +0.0048 0.015656 +0.0049 0.015585 +0.0050 0.015514 +0.0051 0.015444 +0.0052 0.015373 +0.0053 0.015302 +0.0054 0.015232 +0.0055 0.015161 +0.0056 0.015090 +0.0057 0.015019 +0.0058 0.014949 +0.0059 0.014878 +0.0060 0.014807 +0.0061 0.014737 +0.0062 0.014666 +0.0063 0.014595 +0.0064 0.014525 +0.0065 0.014454 +0.0066 0.014383 +0.0067 0.014312 +0.0068 0.014242 +0.0069 0.014171 +0.0070 0.014100 +0.0071 0.014030 +0.0072 0.013959 +0.0073 0.013888 +0.0074 0.013817 +0.0075 0.013747 +0.0076 0.013676 +0.0077 0.013605 +0.0078 0.013535 +0.0079 0.013464 +0.0080 0.013393 +0.0081 0.013322 +0.0082 0.013252 +0.0083 0.013181 +0.0084 0.013110 +0.0085 0.013040 +0.0086 0.012969 +0.0087 0.012898 +0.0088 0.012827 +0.0089 0.012757 +0.0090 0.012686 +0.0091 0.012615 +0.0092 0.012545 +0.0093 0.012474 +0.0094 0.012403 +0.0095 0.012332 +0.0096 0.012262 +0.0097 0.012191 +0.0098 0.012120 +0.0099 0.012050 +0.0100 0.011979 +0.0101 0.011908 +0.0102 0.011838 +0.0103 0.011767 +0.0104 0.011696 +0.0105 0.011625 +0.0106 0.011555 +0.0107 0.011484 +0.0108 0.011413 +0.0109 0.011343 +0.0110 0.011272 +0.0111 0.011201 +0.0112 0.011130 +0.0113 0.011060 +0.0114 0.010989 +0.0115 0.010918 +0.0116 0.010848 +0.0117 0.010777 +0.0118 0.010706 +0.0119 0.010635 +0.0120 0.010565 +0.0121 0.010494 +0.0122 0.010423 +0.0123 0.010353 +0.0124 0.010282 +0.0125 0.010211 +0.0126 0.010140 +0.0127 0.010070 +0.0128 0.009999 +0.0129 0.009928 +0.0130 0.009858 +0.0131 0.009787 +0.0132 0.009716 +0.0133 0.009645 +0.0134 0.009575 +0.0135 0.009504 +0.0136 0.009433 +0.0137 0.009363 +0.0138 0.009292 +0.0139 0.009221 +0.0140 0.009151 +0.0141 0.009080 +0.0142 0.009009 +0.0143 0.008938 +0.0144 0.008868 +0.0145 0.008797 +0.0146 0.008726 +0.0147 0.008656 +0.0148 0.008585 +0.0149 0.008514 +0.0150 0.008443 +0.0151 0.008373 +0.0152 0.008302 +0.0153 0.008231 +0.0154 0.008161 +0.0155 0.008090 +0.0156 0.008019 +0.0157 0.007948 +0.0158 0.007878 +0.0159 0.007807 +0.0160 0.007736 +0.0161 0.007666 +0.0162 0.007595 +0.0163 0.007524 +0.0164 0.007453 +0.0165 0.007383 +0.0166 0.007312 +0.0167 0.007241 +0.0168 0.007171 +0.0169 0.007100 +0.0170 0.007029 +0.0171 0.006958 +0.0172 0.006888 +0.0173 0.006817 +0.0174 0.006746 +0.0175 0.006676 +0.0176 0.006605 +0.0177 0.006534 +0.0178 0.006463 +0.0179 0.006393 +0.0180 0.006322 +0.0181 0.006251 +0.0182 0.006181 +0.0183 0.006110 +0.0184 0.006039 +0.0185 0.005969 +0.0186 0.005898 +0.0187 0.005827 +0.0188 0.005756 +0.0189 0.005686 +0.0190 0.005615 +0.0191 0.005544 +0.0192 0.005474 +0.0193 0.005403 +0.0194 0.005332 +0.0195 0.005261 +0.0196 0.005191 +0.0197 0.005120 +0.0198 0.005049 +0.0199 0.004979 +0.0200 0.004908 +0.0201 0.004837 +0.0202 0.004766 +0.0203 0.004696 +0.0204 0.004625 +0.0205 0.004554 +0.0206 0.004484 +0.0207 0.004413 +0.0208 0.004342 +0.0209 0.004271 +0.0210 0.004201 +0.0211 0.004130 +0.0212 0.004059 +0.0213 0.003989 +0.0214 0.003918 +0.0215 0.003847 +0.0216 0.003776 +0.0217 0.003706 +0.0218 0.003635 +0.0219 0.003564 +0.0220 0.003494 +0.0221 0.003423 +0.0222 0.003352 +0.0223 0.003282 +0.0224 0.003211 +0.0225 0.003140 +0.0226 0.003069 +0.0227 0.002999 +0.0228 0.002928 +0.0229 0.002857 +0.0230 0.002787 +0.0231 0.002716 +0.0232 0.002645 +0.0233 0.002574 +0.0234 0.002504 +0.0235 0.002433 +0.0236 0.002362 +0.0237 0.002292 +0.0238 0.002221 +0.0239 0.002150 +0.0240 0.002079 +0.0241 0.002009 +0.0242 0.001938 +0.0243 0.001867 +0.0244 0.001797 +0.0245 0.001726 +0.0246 0.001655 +0.0247 0.001584 +0.0248 0.001514 +0.0249 0.001443 +0.0250 0.001372 +0.0251 0.001302 +0.0252 0.001231 +0.0253 0.001160 +0.0254 0.001089 +0.0255 0.001019 +0.0256 0.000948 +0.0257 0.000877 +0.0258 0.000807 +0.0259 0.000736 +0.0260 0.000665 +0.0261 0.000635 +0.0262 0.000635 +0.0263 0.000635 +0.0264 0.000661 +0.0265 0.000687 +0.0266 0.000713 +0.0267 0.000739 +0.0268 0.000764 +0.0269 0.000790 +0.0270 0.000816 +0.0271 0.000842 +0.0272 0.000868 +0.0273 0.000894 +0.0274 0.000920 +0.0275 0.000946 +0.0276 0.000971 +0.0277 0.000997 +0.0278 0.001023 +0.0279 0.001049 +0.0280 0.001075 +0.0281 0.001101 +0.0282 0.001127 +0.0283 0.001153 +0.0284 0.001179 +0.0285 0.001204 +0.0286 0.001230 +0.0287 0.001256 +0.0288 0.001282 +0.0289 0.001308 +0.0290 0.001334 +0.0291 0.001360 +0.0292 0.001386 +0.0293 0.001411 +0.0294 0.001437 +0.0295 0.001463 +0.0296 0.001489 +0.0297 0.001515 +0.0298 0.001541 +0.0299 0.001567 +0.0300 0.001593 +0.0301 0.001619 +0.0302 0.001644 +0.0303 0.001670 +0.0304 0.001696 +0.0305 0.001722 +0.0306 0.001748 +0.0307 0.001774 +0.0308 0.001800 +0.0309 0.001826 +0.0310 0.001851 +0.0311 0.001877 +0.0312 0.001903 +0.0313 0.001929 +0.0314 0.001955 +0.0315 0.001981 +0.0316 0.002007 +0.0317 0.002033 +0.0318 0.002059 +0.0319 0.002084 +0.0320 0.002110 +0.0321 0.002136 +0.0322 0.002162 +0.0323 0.002188 +0.0324 0.002214 +0.0325 0.002240 +0.0326 0.002266 +0.0327 0.002291 +0.0328 0.002317 +0.0329 0.002343 +0.0330 0.002369 +0.0331 0.002395 +0.0332 0.002421 +0.0333 0.002447 +0.0334 0.002473 +0.0335 0.002498 +0.0336 0.002524 +0.0337 0.002550 +0.0338 0.002576 +0.0339 0.002602 +0.0340 0.002628 +0.0341 0.002654 +0.0342 0.002680 +0.0343 0.002706 +0.0344 0.002731 +0.0345 0.002757 +0.0346 0.002783 +0.0347 0.002809 +0.0348 0.002835 +0.0349 0.002861 +0.0350 0.002887 +0.0351 0.002913 +0.0352 0.002938 +0.0353 0.002964 +0.0354 0.002990 +0.0355 0.003016 +0.0356 0.003042 +0.0357 0.003068 +0.0358 0.003094 +0.0359 0.003120 +0.0360 0.003146 +0.0361 0.003171 +0.0362 0.003197 +0.0363 0.003223 +0.0364 0.003249 +0.0365 0.003275 +0.0366 0.003301 +0.0367 0.003327 +0.0368 0.003353 +0.0369 0.003378 +0.0370 0.003404 +0.0371 0.003430 +0.0372 0.003456 +0.0373 0.003482 +0.0374 0.003508 +0.0375 0.003534 +0.0376 0.003560 +0.0377 0.003586 +0.0378 0.003611 +0.0379 0.003637 +0.0380 0.003663 +0.0381 0.003689 +0.0382 0.003715 +0.0383 0.003741 +0.0384 0.003767 +0.0385 0.003793 +0.0386 0.003818 +0.0387 0.003844 +0.0388 0.003870 +0.0389 0.003896 +0.0390 0.003922 +0.0391 0.003948 +0.0392 0.003974 +0.0393 0.004000 +0.0394 0.004026 +0.0395 0.004051 +0.0396 0.004077 +0.0397 0.004103 +0.0398 0.004129 +0.0399 0.004155 +0.0400 0.004181 +0.0401 0.004207 +0.0402 0.004233 +0.0403 0.004258 +0.0404 0.004284 +0.0405 0.004310 +0.0406 0.004336 +0.0407 0.004362 +0.0408 0.004388 +0.0409 0.004414 +0.0410 0.004440 +0.0411 0.004466 +0.0412 0.004491 +0.0413 0.004517 +0.0414 0.004543 +0.0415 0.004569 +0.0416 0.004595 +0.0417 0.004621 +0.0418 0.004647 +0.0419 0.004673 +0.0420 0.004698 +0.0421 0.004724 +0.0422 0.004750 +0.0423 0.004776 +0.0424 0.004802 +0.0425 0.004828 +0.0426 0.004854 +0.0427 0.004880 +0.0428 0.004906 +0.0429 0.004931 +0.0430 0.004957 +0.0431 0.004983 +0.0432 0.005009 +0.0433 0.005035 +0.0434 0.005061 +0.0435 0.005087 +0.0436 0.005113 +0.0437 0.005138 +0.0438 0.005164 +0.0439 0.005190 +0.0440 0.005216 +0.0441 0.005242 +0.0442 0.005268 +0.0443 0.005294 +0.0444 0.005320 +0.0445 0.005346 +0.0446 0.005371 +0.0447 0.005397 +0.0448 0.005423 +0.0449 0.005449 +0.0450 0.005475 +0.0451 0.005501 +0.0452 0.005527 +0.0453 0.005553 +0.0454 0.005578 +0.0455 0.005604 +0.0456 0.005630 +0.0457 0.005656 +0.0458 0.005682 +0.0459 0.005708 +0.0460 0.005734 +0.0461 0.005760 +0.0462 0.005785 +0.0463 0.005811 +0.0464 0.005837 +0.0465 0.005863 +0.0466 0.005889 +0.0467 0.005915 +0.0468 0.005941 +0.0469 0.005967 +0.0470 0.005993 +0.0471 0.006018 +0.0472 0.006044 +0.0473 0.006070 +0.0474 0.006096 +0.0475 0.006122 +0.0476 0.006148 +0.0477 0.006174 +0.0478 0.006200 +0.0479 0.006225 +0.0480 0.006251 +0.0481 0.006277 +0.0482 0.006303 +0.0483 0.006329 +0.0484 0.006355 +0.0485 0.006381 +0.0486 0.006407 +0.0487 0.006433 +0.0488 0.006458 +0.0489 0.006484 +0.0490 0.006510 +0.0491 0.006536 +0.0492 0.006562 +0.0493 0.006588 +0.0494 0.006614 +0.0495 0.006640 +0.0496 0.006665 +0.0497 0.006691 +0.0498 0.006717 +0.0499 0.006743 +0.0500 0.006769 +0.0501 0.006795 +0.0502 0.006821 +0.0503 0.006847 +0.0504 0.006873 +0.0505 0.006898 +0.0506 0.006924 +0.0507 0.006950 +0.0508 0.006976 +0.0509 0.007002 +0.0510 0.007028 +0.0511 0.007054 +0.0512 0.007080 +0.0513 0.007105 +0.0514 0.007131 +0.0515 0.007157 +0.0516 0.007183 +0.0517 0.007209 +0.0518 0.007235 +0.0519 0.007261 +0.0520 0.007287 +0.0521 0.007313 +0.0522 0.007338 +0.0523 0.007364 +0.0524 0.007390 +0.0525 0.007416 +0.0526 0.007442 +0.0527 0.007468 +0.0528 0.007494 +0.0529 0.007520 +0.0530 0.007545 +0.0531 0.007571 +0.0532 0.007597 +0.0533 0.007623 +0.0534 0.007649 +0.0535 0.007675 +0.0536 0.007701 +0.0537 0.007727 +0.0538 0.007753 +0.0539 0.007778 +0.0540 0.007804 +0.0541 0.007830 +0.0542 0.007856 +0.0543 0.007882 +0.0544 0.007908 +0.0545 0.007934 +0.0546 0.007960 +0.0547 0.007985 +0.0548 0.008011 +0.0549 0.008037 +0.0550 0.008063 +0.0551 0.008089 +0.0552 0.008115 +0.0553 0.008141 +0.0554 0.008167 +0.0555 0.008193 +0.0556 0.008218 +0.0557 0.008244 +0.0558 0.008270 +0.0559 0.008296 +0.0560 0.008322 +0.0561 0.008348 +0.0562 0.008374 +0.0563 0.008400 +0.0564 0.008425 +0.0565 0.008451 +0.0566 0.008477 +0.0567 0.008503 +0.0568 0.008529 +0.0569 0.008555 +0.0570 0.008581 +0.0571 0.008607 +0.0572 0.008633 +0.0573 0.008658 +0.0574 0.008684 +0.0575 0.008710 +0.0576 0.008736 +0.0577 0.008762 +0.0578 0.008788 +0.0579 0.008814 +0.0580 0.008840 +0.0581 0.008865 +0.0582 0.008891 +0.0583 0.008917 +0.0584 0.008943 +0.0585 0.008969 +0.0586 0.008995 +0.0587 0.009021 +0.0588 0.009047 +0.0589 0.009073 +0.0590 0.009098 +0.0591 0.009124 +0.0592 0.009150 +0.0593 0.009176 +0.0594 0.009202 +0.0595 0.009228 +0.0596 0.009254 +0.0597 0.009280 +0.0598 0.009305 +0.0599 0.009331 +0.0600 0.009357 +0.0601 0.009383 +0.0602 0.009409 +0.0603 0.009435 +0.0604 0.009461 +0.0605 0.009487 +0.0606 0.009512 +0.0607 0.009538 +0.0608 0.009564 +0.0609 0.009590 +0.0610 0.009616 +0.0611 0.009642 +0.0612 0.009668 +0.0613 0.009694 +0.0614 0.009720 +0.0615 0.009745 +0.0616 0.009771 +0.0617 0.009797 +0.0618 0.009823 +0.0619 0.009849 +0.0620 0.009875 +0.0621 0.009901 +0.0622 0.009927 +0.0623 0.009952 +0.0624 0.009978 +0.0625 0.010004 +0.0626 0.010030 +0.0627 0.010056 +0.0628 0.010082 +0.0629 0.010108 +0.0630 0.010134 +0.0631 0.010160 +0.0632 0.010185 +0.0633 0.010211 +0.0634 0.010237 +0.0635 0.010263 +0.0636 0.010289 +0.0637 0.010315 +0.0638 0.010341 +0.0639 0.010367 +0.0640 0.010392 +0.0641 0.010418 +0.0642 0.010444 +0.0643 0.010470 +0.0644 0.010496 +0.0645 0.010522 +0.0646 0.010548 +0.0647 0.010574 +0.0648 0.010600 +0.0649 0.010625 +0.0650 0.010651 +0.0651 0.010677 +0.0652 0.010703 +0.0653 0.010729 +0.0654 0.010755 +0.0655 0.010781 +0.0656 0.010807 +0.0657 0.010832 +0.0658 0.010858 +0.0659 0.010884 +0.0660 0.010910 +0.0661 0.010936 +0.0662 0.010962 +0.0663 0.010988 +0.0664 0.011014 +0.0665 0.011040 +0.0666 0.011065 +0.0667 0.011091 +0.0668 0.011117 +0.0669 0.011143 +0.0670 0.011169 +0.0671 0.011195 +0.0672 0.011221 +0.0673 0.011247 +0.0674 0.011272 +0.0675 0.011298 +0.0676 0.011324 +0.0677 0.011350 +0.0678 0.011376 +0.0679 0.011402 +0.0680 0.011428 +0.0681 0.011454 +0.0682 0.011480 +0.0683 0.011505 +0.0684 0.011531 +0.0685 0.011557 +0.0686 0.011583 +0.0687 0.011609 +0.0688 0.011635 +0.0689 0.011661 +0.0690 0.011687 +0.0691 0.011712 +0.0692 0.011738 +0.0693 0.011764 +0.0694 0.011790 +0.0695 0.011816 +0.0696 0.011842 +0.0697 0.011868 +0.0698 0.011894 +0.0699 0.011920 +0.0700 0.011945 +0.0701 0.011971 +0.0702 0.011997 +0.0703 0.012023 +0.0704 0.012049 +0.0705 0.012075 +0.0706 0.012101 +0.0707 0.012127 +0.0708 0.012152 +0.0709 0.012178 +0.0710 0.012204 +0.0711 0.012230 +0.0712 0.012256 +0.0713 0.012282 +0.0714 0.012308 +0.0715 0.012334 +0.0716 0.012360 +0.0717 0.012385 +0.0718 0.012411 +0.0719 0.012437 +0.0720 0.012463 +0.0721 0.012489 +0.0722 0.012515 +0.0723 0.012541 +0.0724 0.012567 +0.0725 0.012592 +0.0726 0.012618 +0.0727 0.012644 +0.0728 0.012670 +0.0729 0.012696 +0.0730 0.012722 +0.0731 0.012748 +0.0732 0.012774 +0.0733 0.012799 +0.0734 0.012825 +0.0735 0.012851 +0.0736 0.012877 +0.0737 0.012903 +0.0738 0.012929 +0.0739 0.012955 +0.0740 0.012981 +0.0741 0.013007 +0.0742 0.013032 +0.0743 0.013058 +0.0744 0.013084 +0.0745 0.013110 +0.0746 0.013136 +0.0747 0.013162 +0.0748 0.013188 +0.0749 0.013214 +0.0750 0.013239 +0.0751 0.013265 +0.0752 0.013291 +0.0753 0.013317 +0.0754 0.013343 +0.0755 0.013369 +0.0756 0.013395 +0.0757 0.013421 +0.0758 0.013447 +0.0759 0.013472 +0.0760 0.013498 +0.0761 0.013524 +0.0762 0.013550 +0.0763 0.013576 +0.0764 0.013602 +0.0765 0.013628 +0.0766 0.013654 +0.0767 0.013679 +0.0768 0.013705 +0.0769 0.013731 +0.0770 0.013757 +0.0771 0.013783 +0.0772 0.013809 +0.0773 0.013835 +0.0774 0.013861 +0.0775 0.013887 +0.0776 0.013912 +0.0777 0.013938 diff --git a/Matlab/msd-050.tdu b/Matlab/msd-050.tdu new file mode 100644 index 0000000..76ccf54 --- /dev/null +++ b/Matlab/msd-050.tdu @@ -0,0 +1,855 @@ +Nitrogen +1.4 0.0280134 +300 239220 +778 +0.0000 0.019050 +0.0001 0.018979 +0.0002 0.018909 +0.0003 0.018838 +0.0004 0.018767 +0.0005 0.018696 +0.0006 0.018626 +0.0007 0.018555 +0.0008 0.018484 +0.0009 0.018414 +0.0010 0.018343 +0.0011 0.018272 +0.0012 0.018201 +0.0013 0.018131 +0.0014 0.018060 +0.0015 0.017989 +0.0016 0.017919 +0.0017 0.017848 +0.0018 0.017777 +0.0019 0.017706 +0.0020 0.017636 +0.0021 0.017565 +0.0022 0.017494 +0.0023 0.017424 +0.0024 0.017353 +0.0025 0.017282 +0.0026 0.017212 +0.0027 0.017141 +0.0028 0.017070 +0.0029 0.016999 +0.0030 0.016929 +0.0031 0.016858 +0.0032 0.016787 +0.0033 0.016717 +0.0034 0.016646 +0.0035 0.016575 +0.0036 0.016504 +0.0037 0.016434 +0.0038 0.016363 +0.0039 0.016292 +0.0040 0.016222 +0.0041 0.016151 +0.0042 0.016080 +0.0043 0.016009 +0.0044 0.015939 +0.0045 0.015868 +0.0046 0.015797 +0.0047 0.015727 +0.0048 0.015656 +0.0049 0.015585 +0.0050 0.015514 +0.0051 0.015444 +0.0052 0.015373 +0.0053 0.015302 +0.0054 0.015232 +0.0055 0.015161 +0.0056 0.015090 +0.0057 0.015019 +0.0058 0.014949 +0.0059 0.014878 +0.0060 0.014807 +0.0061 0.014737 +0.0062 0.014666 +0.0063 0.014595 +0.0064 0.014525 +0.0065 0.014454 +0.0066 0.014383 +0.0067 0.014312 +0.0068 0.014242 +0.0069 0.014171 +0.0070 0.014100 +0.0071 0.014030 +0.0072 0.013959 +0.0073 0.013888 +0.0074 0.013817 +0.0075 0.013747 +0.0076 0.013676 +0.0077 0.013605 +0.0078 0.013535 +0.0079 0.013464 +0.0080 0.013393 +0.0081 0.013322 +0.0082 0.013252 +0.0083 0.013181 +0.0084 0.013110 +0.0085 0.013040 +0.0086 0.012969 +0.0087 0.012898 +0.0088 0.012827 +0.0089 0.012757 +0.0090 0.012686 +0.0091 0.012615 +0.0092 0.012545 +0.0093 0.012474 +0.0094 0.012403 +0.0095 0.012332 +0.0096 0.012262 +0.0097 0.012191 +0.0098 0.012120 +0.0099 0.012050 +0.0100 0.011979 +0.0101 0.011908 +0.0102 0.011838 +0.0103 0.011767 +0.0104 0.011696 +0.0105 0.011625 +0.0106 0.011555 +0.0107 0.011484 +0.0108 0.011413 +0.0109 0.011343 +0.0110 0.011272 +0.0111 0.011201 +0.0112 0.011130 +0.0113 0.011060 +0.0114 0.010989 +0.0115 0.010918 +0.0116 0.010848 +0.0117 0.010777 +0.0118 0.010706 +0.0119 0.010635 +0.0120 0.010565 +0.0121 0.010494 +0.0122 0.010423 +0.0123 0.010353 +0.0124 0.010282 +0.0125 0.010211 +0.0126 0.010140 +0.0127 0.010070 +0.0128 0.009999 +0.0129 0.009928 +0.0130 0.009858 +0.0131 0.009787 +0.0132 0.009716 +0.0133 0.009645 +0.0134 0.009575 +0.0135 0.009504 +0.0136 0.009433 +0.0137 0.009363 +0.0138 0.009292 +0.0139 0.009221 +0.0140 0.009151 +0.0141 0.009080 +0.0142 0.009009 +0.0143 0.008938 +0.0144 0.008868 +0.0145 0.008797 +0.0146 0.008726 +0.0147 0.008656 +0.0148 0.008585 +0.0149 0.008514 +0.0150 0.008443 +0.0151 0.008373 +0.0152 0.008302 +0.0153 0.008231 +0.0154 0.008161 +0.0155 0.008090 +0.0156 0.008019 +0.0157 0.007948 +0.0158 0.007878 +0.0159 0.007807 +0.0160 0.007736 +0.0161 0.007666 +0.0162 0.007595 +0.0163 0.007524 +0.0164 0.007453 +0.0165 0.007383 +0.0166 0.007312 +0.0167 0.007241 +0.0168 0.007171 +0.0169 0.007100 +0.0170 0.007029 +0.0171 0.006958 +0.0172 0.006888 +0.0173 0.006817 +0.0174 0.006746 +0.0175 0.006676 +0.0176 0.006605 +0.0177 0.006534 +0.0178 0.006463 +0.0179 0.006393 +0.0180 0.006322 +0.0181 0.006251 +0.0182 0.006181 +0.0183 0.006110 +0.0184 0.006039 +0.0185 0.005969 +0.0186 0.005898 +0.0187 0.005827 +0.0188 0.005756 +0.0189 0.005686 +0.0190 0.005615 +0.0191 0.005544 +0.0192 0.005474 +0.0193 0.005403 +0.0194 0.005332 +0.0195 0.005261 +0.0196 0.005191 +0.0197 0.005120 +0.0198 0.005049 +0.0199 0.004979 +0.0200 0.004908 +0.0201 0.004837 +0.0202 0.004766 +0.0203 0.004696 +0.0204 0.004625 +0.0205 0.004554 +0.0206 0.004484 +0.0207 0.004413 +0.0208 0.004342 +0.0209 0.004271 +0.0210 0.004201 +0.0211 0.004130 +0.0212 0.004059 +0.0213 0.003989 +0.0214 0.003918 +0.0215 0.003847 +0.0216 0.003776 +0.0217 0.003706 +0.0218 0.003635 +0.0219 0.003564 +0.0220 0.003494 +0.0221 0.003423 +0.0222 0.003352 +0.0223 0.003282 +0.0224 0.003211 +0.0225 0.003140 +0.0226 0.003069 +0.0227 0.002999 +0.0228 0.002928 +0.0229 0.002857 +0.0230 0.002787 +0.0231 0.002716 +0.0232 0.002645 +0.0233 0.002574 +0.0234 0.002504 +0.0235 0.002433 +0.0236 0.002362 +0.0237 0.002292 +0.0238 0.002221 +0.0239 0.002150 +0.0240 0.002079 +0.0241 0.002009 +0.0242 0.001938 +0.0243 0.001867 +0.0244 0.001797 +0.0245 0.001726 +0.0246 0.001655 +0.0247 0.001584 +0.0248 0.001514 +0.0249 0.001443 +0.0250 0.001372 +0.0251 0.001302 +0.0252 0.001270 +0.0253 0.001270 +0.0254 0.001270 +0.0255 0.001296 +0.0256 0.001322 +0.0257 0.001348 +0.0258 0.001374 +0.0259 0.001399 +0.0260 0.001425 +0.0261 0.001451 +0.0262 0.001477 +0.0263 0.001503 +0.0264 0.001529 +0.0265 0.001555 +0.0266 0.001581 +0.0267 0.001606 +0.0268 0.001632 +0.0269 0.001658 +0.0270 0.001684 +0.0271 0.001710 +0.0272 0.001736 +0.0273 0.001762 +0.0274 0.001788 +0.0275 0.001814 +0.0276 0.001839 +0.0277 0.001865 +0.0278 0.001891 +0.0279 0.001917 +0.0280 0.001943 +0.0281 0.001969 +0.0282 0.001995 +0.0283 0.002021 +0.0284 0.002046 +0.0285 0.002072 +0.0286 0.002098 +0.0287 0.002124 +0.0288 0.002150 +0.0289 0.002176 +0.0290 0.002202 +0.0291 0.002228 +0.0292 0.002254 +0.0293 0.002279 +0.0294 0.002305 +0.0295 0.002331 +0.0296 0.002357 +0.0297 0.002383 +0.0298 0.002409 +0.0299 0.002435 +0.0300 0.002461 +0.0301 0.002486 +0.0302 0.002512 +0.0303 0.002538 +0.0304 0.002564 +0.0305 0.002590 +0.0306 0.002616 +0.0307 0.002642 +0.0308 0.002668 +0.0309 0.002694 +0.0310 0.002719 +0.0311 0.002745 +0.0312 0.002771 +0.0313 0.002797 +0.0314 0.002823 +0.0315 0.002849 +0.0316 0.002875 +0.0317 0.002901 +0.0318 0.002926 +0.0319 0.002952 +0.0320 0.002978 +0.0321 0.003004 +0.0322 0.003030 +0.0323 0.003056 +0.0324 0.003082 +0.0325 0.003108 +0.0326 0.003133 +0.0327 0.003159 +0.0328 0.003185 +0.0329 0.003211 +0.0330 0.003237 +0.0331 0.003263 +0.0332 0.003289 +0.0333 0.003315 +0.0334 0.003341 +0.0335 0.003366 +0.0336 0.003392 +0.0337 0.003418 +0.0338 0.003444 +0.0339 0.003470 +0.0340 0.003496 +0.0341 0.003522 +0.0342 0.003548 +0.0343 0.003573 +0.0344 0.003599 +0.0345 0.003625 +0.0346 0.003651 +0.0347 0.003677 +0.0348 0.003703 +0.0349 0.003729 +0.0350 0.003755 +0.0351 0.003781 +0.0352 0.003806 +0.0353 0.003832 +0.0354 0.003858 +0.0355 0.003884 +0.0356 0.003910 +0.0357 0.003936 +0.0358 0.003962 +0.0359 0.003988 +0.0360 0.004013 +0.0361 0.004039 +0.0362 0.004065 +0.0363 0.004091 +0.0364 0.004117 +0.0365 0.004143 +0.0366 0.004169 +0.0367 0.004195 +0.0368 0.004221 +0.0369 0.004246 +0.0370 0.004272 +0.0371 0.004298 +0.0372 0.004324 +0.0373 0.004350 +0.0374 0.004376 +0.0375 0.004402 +0.0376 0.004428 +0.0377 0.004453 +0.0378 0.004479 +0.0379 0.004505 +0.0380 0.004531 +0.0381 0.004557 +0.0382 0.004583 +0.0383 0.004609 +0.0384 0.004635 +0.0385 0.004661 +0.0386 0.004686 +0.0387 0.004712 +0.0388 0.004738 +0.0389 0.004764 +0.0390 0.004790 +0.0391 0.004816 +0.0392 0.004842 +0.0393 0.004868 +0.0394 0.004893 +0.0395 0.004919 +0.0396 0.004945 +0.0397 0.004971 +0.0398 0.004997 +0.0399 0.005023 +0.0400 0.005049 +0.0401 0.005075 +0.0402 0.005101 +0.0403 0.005126 +0.0404 0.005152 +0.0405 0.005178 +0.0406 0.005204 +0.0407 0.005230 +0.0408 0.005256 +0.0409 0.005282 +0.0410 0.005308 +0.0411 0.005333 +0.0412 0.005359 +0.0413 0.005385 +0.0414 0.005411 +0.0415 0.005437 +0.0416 0.005463 +0.0417 0.005489 +0.0418 0.005515 +0.0419 0.005541 +0.0420 0.005566 +0.0421 0.005592 +0.0422 0.005618 +0.0423 0.005644 +0.0424 0.005670 +0.0425 0.005696 +0.0426 0.005722 +0.0427 0.005748 +0.0428 0.005773 +0.0429 0.005799 +0.0430 0.005825 +0.0431 0.005851 +0.0432 0.005877 +0.0433 0.005903 +0.0434 0.005929 +0.0435 0.005955 +0.0436 0.005981 +0.0437 0.006006 +0.0438 0.006032 +0.0439 0.006058 +0.0440 0.006084 +0.0441 0.006110 +0.0442 0.006136 +0.0443 0.006162 +0.0444 0.006188 +0.0445 0.006213 +0.0446 0.006239 +0.0447 0.006265 +0.0448 0.006291 +0.0449 0.006317 +0.0450 0.006343 +0.0451 0.006369 +0.0452 0.006395 +0.0453 0.006420 +0.0454 0.006446 +0.0455 0.006472 +0.0456 0.006498 +0.0457 0.006524 +0.0458 0.006550 +0.0459 0.006576 +0.0460 0.006602 +0.0461 0.006628 +0.0462 0.006653 +0.0463 0.006679 +0.0464 0.006705 +0.0465 0.006731 +0.0466 0.006757 +0.0467 0.006783 +0.0468 0.006809 +0.0469 0.006835 +0.0470 0.006860 +0.0471 0.006886 +0.0472 0.006912 +0.0473 0.006938 +0.0474 0.006964 +0.0475 0.006990 +0.0476 0.007016 +0.0477 0.007042 +0.0478 0.007068 +0.0479 0.007093 +0.0480 0.007119 +0.0481 0.007145 +0.0482 0.007171 +0.0483 0.007197 +0.0484 0.007223 +0.0485 0.007249 +0.0486 0.007275 +0.0487 0.007300 +0.0488 0.007326 +0.0489 0.007352 +0.0490 0.007378 +0.0491 0.007404 +0.0492 0.007430 +0.0493 0.007456 +0.0494 0.007482 +0.0495 0.007508 +0.0496 0.007533 +0.0497 0.007559 +0.0498 0.007585 +0.0499 0.007611 +0.0500 0.007637 +0.0501 0.007663 +0.0502 0.007689 +0.0503 0.007715 +0.0504 0.007740 +0.0505 0.007766 +0.0506 0.007792 +0.0507 0.007818 +0.0508 0.007844 +0.0509 0.007870 +0.0510 0.007896 +0.0511 0.007922 +0.0512 0.007948 +0.0513 0.007973 +0.0514 0.007999 +0.0515 0.008025 +0.0516 0.008051 +0.0517 0.008077 +0.0518 0.008103 +0.0519 0.008129 +0.0520 0.008155 +0.0521 0.008180 +0.0522 0.008206 +0.0523 0.008232 +0.0524 0.008258 +0.0525 0.008284 +0.0526 0.008310 +0.0527 0.008336 +0.0528 0.008362 +0.0529 0.008388 +0.0530 0.008413 +0.0531 0.008439 +0.0532 0.008465 +0.0533 0.008491 +0.0534 0.008517 +0.0535 0.008543 +0.0536 0.008569 +0.0537 0.008595 +0.0538 0.008620 +0.0539 0.008646 +0.0540 0.008672 +0.0541 0.008698 +0.0542 0.008724 +0.0543 0.008750 +0.0544 0.008776 +0.0545 0.008802 +0.0546 0.008828 +0.0547 0.008853 +0.0548 0.008879 +0.0549 0.008905 +0.0550 0.008931 +0.0551 0.008957 +0.0552 0.008983 +0.0553 0.009009 +0.0554 0.009035 +0.0555 0.009060 +0.0556 0.009086 +0.0557 0.009112 +0.0558 0.009138 +0.0559 0.009164 +0.0560 0.009190 +0.0561 0.009216 +0.0562 0.009242 +0.0563 0.009268 +0.0564 0.009293 +0.0565 0.009319 +0.0566 0.009345 +0.0567 0.009371 +0.0568 0.009397 +0.0569 0.009423 +0.0570 0.009449 +0.0571 0.009475 +0.0572 0.009500 +0.0573 0.009526 +0.0574 0.009552 +0.0575 0.009578 +0.0576 0.009604 +0.0577 0.009630 +0.0578 0.009656 +0.0579 0.009682 +0.0580 0.009708 +0.0581 0.009733 +0.0582 0.009759 +0.0583 0.009785 +0.0584 0.009811 +0.0585 0.009837 +0.0586 0.009863 +0.0587 0.009889 +0.0588 0.009915 +0.0589 0.009940 +0.0590 0.009966 +0.0591 0.009992 +0.0592 0.010018 +0.0593 0.010044 +0.0594 0.010070 +0.0595 0.010096 +0.0596 0.010122 +0.0597 0.010147 +0.0598 0.010173 +0.0599 0.010199 +0.0600 0.010225 +0.0601 0.010251 +0.0602 0.010277 +0.0603 0.010303 +0.0604 0.010329 +0.0605 0.010355 +0.0606 0.010380 +0.0607 0.010406 +0.0608 0.010432 +0.0609 0.010458 +0.0610 0.010484 +0.0611 0.010510 +0.0612 0.010536 +0.0613 0.010562 +0.0614 0.010587 +0.0615 0.010613 +0.0616 0.010639 +0.0617 0.010665 +0.0618 0.010691 +0.0619 0.010717 +0.0620 0.010743 +0.0621 0.010769 +0.0622 0.010795 +0.0623 0.010820 +0.0624 0.010846 +0.0625 0.010872 +0.0626 0.010898 +0.0627 0.010924 +0.0628 0.010950 +0.0629 0.010976 +0.0630 0.011002 +0.0631 0.011027 +0.0632 0.011053 +0.0633 0.011079 +0.0634 0.011105 +0.0635 0.011131 +0.0636 0.011157 +0.0637 0.011183 +0.0638 0.011209 +0.0639 0.011235 +0.0640 0.011260 +0.0641 0.011286 +0.0642 0.011312 +0.0643 0.011338 +0.0644 0.011364 +0.0645 0.011390 +0.0646 0.011416 +0.0647 0.011442 +0.0648 0.011467 +0.0649 0.011493 +0.0650 0.011519 +0.0651 0.011545 +0.0652 0.011571 +0.0653 0.011597 +0.0654 0.011623 +0.0655 0.011649 +0.0656 0.011675 +0.0657 0.011700 +0.0658 0.011726 +0.0659 0.011752 +0.0660 0.011778 +0.0661 0.011804 +0.0662 0.011830 +0.0663 0.011856 +0.0664 0.011882 +0.0665 0.011907 +0.0666 0.011933 +0.0667 0.011959 +0.0668 0.011985 +0.0669 0.012011 +0.0670 0.012037 +0.0671 0.012063 +0.0672 0.012089 +0.0673 0.012115 +0.0674 0.012140 +0.0675 0.012166 +0.0676 0.012192 +0.0677 0.012218 +0.0678 0.012244 +0.0679 0.012270 +0.0680 0.012296 +0.0681 0.012322 +0.0682 0.012347 +0.0683 0.012373 +0.0684 0.012399 +0.0685 0.012425 +0.0686 0.012451 +0.0687 0.012477 +0.0688 0.012503 +0.0689 0.012529 +0.0690 0.012555 +0.0691 0.012580 +0.0692 0.012606 +0.0693 0.012632 +0.0694 0.012658 +0.0695 0.012684 +0.0696 0.012710 +0.0697 0.012736 +0.0698 0.012762 +0.0699 0.012787 +0.0700 0.012813 +0.0701 0.012839 +0.0702 0.012865 +0.0703 0.012891 +0.0704 0.012917 +0.0705 0.012943 +0.0706 0.012969 +0.0707 0.012995 +0.0708 0.013020 +0.0709 0.013046 +0.0710 0.013072 +0.0711 0.013098 +0.0712 0.013124 +0.0713 0.013150 +0.0714 0.013176 +0.0715 0.013202 +0.0716 0.013227 +0.0717 0.013253 +0.0718 0.013279 +0.0719 0.013305 +0.0720 0.013331 +0.0721 0.013357 +0.0722 0.013383 +0.0723 0.013409 +0.0724 0.013434 +0.0725 0.013460 +0.0726 0.013486 +0.0727 0.013512 +0.0728 0.013538 +0.0729 0.013564 +0.0730 0.013590 +0.0731 0.013616 +0.0732 0.013642 +0.0733 0.013667 +0.0734 0.013693 +0.0735 0.013719 +0.0736 0.013745 +0.0737 0.013771 +0.0738 0.013797 +0.0739 0.013823 +0.0740 0.013849 +0.0741 0.013874 +0.0742 0.013900 +0.0743 0.013926 +0.0744 0.013952 +0.0745 0.013978 +0.0746 0.014004 +0.0747 0.014030 +0.0748 0.014056 +0.0749 0.014082 +0.0750 0.014107 +0.0751 0.014133 +0.0752 0.014159 +0.0753 0.014185 +0.0754 0.014211 +0.0755 0.014237 +0.0756 0.014263 +0.0757 0.014289 +0.0758 0.014314 +0.0759 0.014340 +0.0760 0.014366 +0.0761 0.014392 +0.0762 0.014418 +0.0763 0.014444 +0.0764 0.014470 +0.0765 0.014496 +0.0766 0.014522 +0.0767 0.014547 +0.0768 0.014573 +0.0769 0.014599 +0.0770 0.014625 +0.0771 0.014651 +0.0772 0.014677 +0.0773 0.014703 +0.0774 0.014729 +0.0775 0.014754 +0.0776 0.014780 +0.0777 0.014806 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Matlab/msd.tdu b/Matlab/msd.tdu deleted file mode 100644 index c3157bb..0000000 --- a/Matlab/msd.tdu +++ /dev/null @@ -1,155 +0,0 @@ -Nitrogen -1.4 .0280134 -300 23922.146 -151 -0 0.009525 -0.01 0.009254657 -0.02 0.008984314 -0.03 0.008713971 -0.04 0.008443629 -0.05 0.008173286 -0.06 0.007902943 -0.07 0.0076326 -0.08 0.007362257 -0.09 0.007091914 -0.1 0.006821571 -0.11 0.006551229 -0.12 0.006280886 -0.13 0.006010543 -0.14 0.0057402 -0.15 0.005469857 -0.16 0.005199514 -0.17 0.004929171 -0.18 0.004658829 -0.19 0.004388486 -0.2 0.004118143 -0.21 0.0038478 -0.22 0.003577457 -0.23 0.003307114 -0.24 0.003036771 -0.25 0.002766429 -0.26 0.002496086 -0.27 0.002225743 -0.28 0.0019554 -0.29 0.001685057 -0.3 0.001414714 -0.31 0.001144371 -0.32 0.000874029 -0.33 0.000603686 -0.34 0.000333343 -0.35 0.000063 -0.36 0.000063 -0.37 0.000127754 -0.38 0.000192509 -0.39 0.000257263 -0.4 0.000322018 -0.41 0.000386772 -0.42 0.000451526 -0.43 0.000516281 -0.44 0.000581035 -0.45 0.000645789 -0.46 0.000710544 -0.47 0.000775298 -0.48 0.000840053 -0.49 0.000904807 -0.5 0.000969561 -0.51 0.001034316 -0.52 0.00109907 -0.53 0.001163825 -0.54 0.001228579 -0.55 0.001293333 -0.56 0.001358088 -0.57 0.001422842 -0.58 0.001487596 -0.59 0.001552351 -0.6 0.001617105 -0.61 0.00168186 -0.62 0.001746614 -0.63 0.001811368 -0.64 0.001876123 -0.65 0.001940877 -0.66 0.002005632 -0.67 0.002070386 -0.68 0.00213514 -0.69 0.002199895 -0.7 0.002264649 -0.71 0.002329404 -0.72 0.002394158 -0.73 0.002458912 -0.74 0.002523667 -0.75 0.002588421 -0.76 0.002653175 -0.77 0.00271793 -0.78 0.002782684 -0.79 0.002847439 -0.8 0.002912193 -0.81 0.002976947 -0.82 0.003041702 -0.83 0.003106456 -0.84 0.003171211 -0.85 0.003235965 -0.86 0.003300719 -0.87 0.003365474 -0.88 0.003430228 -0.89 0.003494982 -0.9 0.003559737 -0.91 0.003624491 -0.92 0.003689246 -0.93 0.003754 -0.94 0.003818754 -0.95 0.003883509 -0.96 0.003948263 -0.97 0.004013018 -0.98 0.004077772 -0.99 0.004142526 -1 0.004207281 -1.01 0.004272035 -1.02 0.004336789 -1.03 0.004401544 -1.04 0.004466298 -1.05 0.004531053 -1.06 0.004595807 -1.07 0.004660561 -1.08 0.004725316 -1.09 0.00479007 -1.1 0.004854825 -1.11 0.004919579 -1.12 0.004984333 -1.13 0.005049088 -1.14 0.005113842 -1.15 0.005178596 -1.16 0.005243351 -1.17 0.005308105 -1.18 0.00537286 -1.19 0.005437614 -1.2 0.005502368 -1.21 0.005567123 -1.22 0.005631877 -1.23 0.005696632 -1.24 0.005761386 -1.25 0.00582614 -1.26 0.005890895 -1.27 0.005955649 -1.28 0.006020404 -1.29 0.006085158 -1.3 0.006149912 -1.31 0.006214667 -1.32 0.006279421 -1.33 0.006344175 -1.34 0.00640893 -1.35 0.006473684 -1.36 0.006538439 -1.37 0.006603193 -1.38 0.006667947 -1.39 0.006732702 -1.4 0.006797456 -1.41 0.006862211 -1.42 0.006926965 -1.43 0.006991719 -1.44 0.007056474 -1.45 0.007121228 -1.46 0.007185982 -1.47 0.007250737 -1.48 0.007315491 -1.49 0.007380246 -1.5 0.007445 \ No newline at end of file diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 4f6c32f..54245e9 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -3,7 +3,7 @@ close all format long - filein = 'sample.tdu'; % input file + filein = 'msd-050.tdu'; % input file % universal constants R_0 = 8.3144598; % [J/(mol*K)] universal gas constant @@ -119,11 +119,14 @@ % subplot(numplots,1,plotcounter) % semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; subplot(numplots,1,plotcounter) - plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + semilogy(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); subplot(numplots,1,plotcounter) plot(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% semilogy(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); subplot(numplots,1,plotcounter) plot(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% semilogy(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % figure % semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); % figure @@ -167,6 +170,7 @@ % 'v/at',exit_velocity/throat_velocity,unitless; }; display(result); + end function Mach = arearatio2mach_sub(A,A_t,k) @@ -210,9 +214,12 @@ Mach = 1/sqrt(X); end function display(result) + global debug + if debug;fprintf('displaying results...\n');end; [n,~]=size(result); for i = 1:n fprintf('\n%24s\t%15.8f\t%s',result{i,:}); end fprintf('\n') + if debug;fprintf('done.\n');end; end \ No newline at end of file From c894f86ad0a4b4ec450139c579afefbb6aeeb753 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 13:29:47 -0400 Subject: [PATCH 08/21] refine plots --- Matlab/msd-025.tdu | 2 +- Matlab/msd-050.tdu | 2 +- Matlab/tdu.asv | 250 +++++++++++++++++++++++++++++++++++++ Matlab/tdu.m | 303 ++++++++++++++++++++++++--------------------- 4 files changed, 416 insertions(+), 141 deletions(-) create mode 100644 Matlab/tdu.asv diff --git a/Matlab/msd-025.tdu b/Matlab/msd-025.tdu index 8d62771..f1d9acd 100644 --- a/Matlab/msd-025.tdu +++ b/Matlab/msd-025.tdu @@ -1,6 +1,6 @@ Nitrogen 1.4 0.0280134 -300 239220 +300 239220 101325 778 0.0000 0.019050 0.0001 0.018979 diff --git a/Matlab/msd-050.tdu b/Matlab/msd-050.tdu index 76ccf54..7f0f4ca 100644 --- a/Matlab/msd-050.tdu +++ b/Matlab/msd-050.tdu @@ -1,6 +1,6 @@ Nitrogen 1.4 0.0280134 -300 239220 +300 239220 101325 778 0.0000 0.019050 0.0001 0.018979 diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv new file mode 100644 index 0000000..841a704 --- /dev/null +++ b/Matlab/tdu.asv @@ -0,0 +1,250 @@ +% HEADER +clc; +close all +format long + +filein = 'msd-050.tdu'; % input file + +% universal constants +R_0 = 8.3144598; % [J/(mol*K)] universal gas constant +g_0 = 9.81; % [m/s^2] standard gravity +unitless='[-]'; +global debug; +debug = true; +% === THRUSTER PARAMETERS === +% IMPORT FROM FILE +fid = fopen(filein,'r'); +if debug;fprintf('reading data from input file (%s)...\n',filein);end +prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) + if debug;fprintf('\tPropellant:\t\t%8s\n',prop_name);end +prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters + k = prop_params(1,1); % specific heat ratio + mw = prop_params(1,2); % molecular weight + if debug;fprintf('\tk:\t\t%16g\n\tmw:\t\t%16g\n',k,mw);end +total_params= fscanf(fid,'%g',[1 3]); % scan total/stagnation parameters + T_0 = total_params(1,1); % total temperature + P_0 = total_params(1,2); % total pressure + P_b = total_params(1,3); % back pressure + if debug;fprintf('\tT_0:\t%16g\n\tP_0:\t%16g\n',T_0,P_0);end +geom_size = fscanf(fid,'%g',[1 1]); % number of geometry nodes +xcoord = zeros(geom_size,1); radius = zeros(geom_size,1); + for i=1:geom_size + geom = fscanf(fid,'%g',[1 2]); + xcoord(i) = geom(1,1); % x coordinate of geometry node + radius(i) = geom(1,2); % radius at xcoord + end + if debug + fprintf('\tinlet radius:\t%8f\n\tthroat radius:\t%8f\n\texit radius:\t%8f\n',radius(1),min(radius),radius(end)); + fprintf('\tlength:\t%16f\n\tgeometry nodes:\t%8i\n',xcoord(end),geom_size); + end +fclose('all'); %close input file +if debug;fprintf('input file closed.\n');end +% MANUAL ENTRY +mw_units = '[kg/mol]'; +temperature_units = '[K]'; +pressure_units = '[Pa]'; +length_units = '[m]'; +angle_units = '[deg]'; +% ===------------=== + +% NOZZLE NOMENCLATURE +% ******************************** +% +% /- +% / 0 = total parameters +% ===\---/ 1 = chamber +% (1) (2) (3) 2 = throat +% ===/---\ 3 = exit +% \ +% \- +% ******************************** + +% % ASSUMPTIONS +% - Isentropic flow +% - Initial temperature and pressure are total parameters + +R = R_0/mw; % gas constant +R_units = '[J/kg K]'; +% A_e = pi*exit_radius^2; % exit area +% A_t = pi*throat_radius^2; % throat area +A=pi.*radius.^2; +area_units = '[m^2]'; + +[A_t,A_t_idx]=min(A); +T_star=T_0*(2/(k+1)); %K +P_star=(2/(k+1))^(k/(k-1)); %Pa +rho_star=P_star/(R*T_star); % kg/m^3 +mdot=rho_star*sqrt(k*R*T_star)*A_t; %choked + +M_idx=linspace(.1,4,length(xcoord)); +k1=(2/(k+1)); +k2=((k-1)/2); +k3=(.5*(k+1)/(k-1)); +area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); + +M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; +M(1)=1; +choked=false; +if debug;fprintf('Computing flow conditions along x-axis...\n');end +for x=2:length(xcoord) + M_sub(x)=arearatio2mach_sub(A(x),A_t,k); + if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end; + if M(x)<1 + M(x)=M_sub(x); + elseif M(x)==1 + choked=true; + if debug;fprintf('\tFlow is choked at A = %f\n',A(x));end + M(x)=M_sub(x); + elseif M(x)>1 + M(x)=M_sup(x); + end + temp_ratio(x)=(1+((k-1)/2)*M(x)^2); + T(x)=T_0/temp_ratio(x); + pres_ratio(x)=temp_ratio(x)^(k/(k-1)); + P(x)=P_0/pres_ratio(x); + + temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); + T_sub(x)=T_0/temp_ratio_sub(x); + pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); + P_sub(x)=P_0/pres_ratio_sub(x); + temp_ratio_sup(x)=(1+((k-1)/2)*M_sup(x)^2); + T_sup(x)=T_0/temp_ratio_sup(x); + pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); + P_sup(x)=P_0/pres_ratio_sup(x); + +end +if debug;fprintf('done.\n');end +if debug + figure + numplots=4;plotcounter=1; + subplot(numplots,1,plotcounter) + plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% subplot(numplots,1,plotcounter) +% semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; + subplot(numplots,1,plotcounter) +% plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--');legend('M','M<1','M>1');ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + subplot(numplots,1,plotcounter) + plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end));ylabel('T (K)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% semilogy(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + subplot(numplots,1,plotcounter) + plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--');legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% semilogy(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +% figure +% semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); +% figure +% semilogy(M_idx,area_ratio,M,A./A_t,'--') + +figure + plot(xcoord,radius);ylabel('radius');axis([0 xcoord(end) 0 inf]); +figure + semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--'); + legend('M','M<1','M>1');ylabel('M');axis([0 xcoord(end) 0 inf]); +figure + plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); + ylabel('T (K)');axis([0 xcoord(end) 0 inf]); +figure + plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); + legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); +end +% % format & display outputs +linedivider='------------'; +result = {'Propellant','',prop_name; + linedivider,'',''; + 'Specific heat ratio', k, unitless; + 'Molar mass', mw, mw_units; + 'Specific gas constant',R,R_units; + linedivider,'',''; + 'Total temperature', T_0, temperature_units; + 'Total pressure', P_0, pressure_units; + linedivider,'',''; + 'Length',xcoord(end),length_units; + 'Inlet radius',radius(1),length_units; + 'Throat radius', min(radius), length_units; + 'Exit radius', radius(end), length_units; + 'Inlet area',A(1),area_units; + 'Throat area',A_t,area_units; + 'Exit area',A(end),area_units; +% 'Half-angle',alpha, angle_units; + linedivider,'',''; + 'Throat temperature',T(A_t_idx),temperature_units; + 'Throat pressure',P(A_t_idx),pressure_units; +% 'Mass flow rate',mass_flowrate,mass_flowrate_units; + linedivider,'',''; + 'Exit temperature',T(end),temperature_units; + 'Exit pressure',P(end),pressure_units; + linedivider,'',''; +% 'Exhaust velocity',exit_velocity,velocity_units; +% 'Thrust',thrust,force_units; +% 'Specific impulse',specific_impulse,isp_units; + linedivider,'',''; + 'Exit Mach',M(end),unitless; + 'A/At',A(end)/A_t,unitless; + 'T/T0',T(end)/T_0,unitless; + 'P/P0',P(end)/P_0,unitless; +% 'v/at',exit_velocity/throat_velocity,unitless; + }; +display(result); +figure + g=1.4; + g1=2/(g+1); + g2=(g-1)/2; + g3=0.5*(g+1)/(g-1); + M=0.1:0.1:4; + M=linspace(0.1,4); + A_theoretical=(1./M).*(g1*(1+g2*M.^2)).^g3; + A_trunc = round(A_theoretical.*1000)./1000; + semilogx(M,A_theoretical,M(find(A_trunc==1)),,'o','linewidth',2) + grid on + xlabel('Mach Number') + ylabel('Area Ratio') +function Mach = arearatio2mach_sub(A,A_t,k) +% solve Mach number from area ratio by Newton-Raphson Method. (assume +% subsonic) +% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html + P = 2/(k+1); + Q = 1-P; + R = (A/A_t).^2; + a = P.^(1/Q); + r = (R-1)/(2*a); + X = 1/((1+r)+sqrt(r*(r+2))); % initial guess + diff = 1; % initalize termination criteria + while abs(diff) > .00001 + F = (P*X+Q).^(1/P)-R*X; + dF = (P*X+Q).^((1/P)-1)-R; + Xnew = X - F/dF; + diff = Xnew - X; + X = Xnew; + end + Mach = sqrt(X); +end +function Mach = arearatio2mach_sup(A,A_t,k) +% solve Mach number from area ratio by Newton-Raphson Method. (assume +% supersonic) +% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html + P = 2/(k+1); + Q = 1-P; + R = (A/A_t).^((2*Q)/P); + a = Q.^(1/P); + r = (R-1)/(2*a); + X = 1/((1+r)+sqrt(r*(r+2))); % initial guess + diff = 1; % initalize termination criteria + while abs(diff) > .00001 + F = (P*X+Q).^(1/P)-R*X; + dF = (P*X+Q).^((1/P)-1)-R; + Xnew = X - F/dF; + diff = Xnew - X; + X = Xnew; + end + Mach = 1/sqrt(X); +end +function display(result) + global debug + if debug;fprintf('displaying results...\n');end; + [n,~]=size(result); + for i = 1:n + fprintf('\n%24s\t%15.8f\t%s',result{i,:}); + end + fprintf('\n') + if debug;fprintf('done.\n');end; +end \ No newline at end of file diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 54245e9..0007af1 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -1,49 +1,50 @@ -function varargout = tdu % main function - clc; - close all - format long - - filein = 'msd-050.tdu'; % input file - - % universal constants - R_0 = 8.3144598; % [J/(mol*K)] universal gas constant - g_0 = 9.81; % [m/s^2] standard gravity - unitless='[-]'; - global debug; - debug = true; +% HEADER +clc; +close all +format long + +filein = 'msd-050.tdu'; % input file + +% universal constants +R_0 = 8.3144598; % [J/(mol*K)] universal gas constant +g_0 = 9.81; % [m/s^2] standard gravity +unitless='[-]'; +global debug; +debug = true; % === THRUSTER PARAMETERS === - % IMPORT FROM FILE - fid = fopen(filein,'r'); - if debug;fprintf('reading data from input file (%s)...\n',filein);end - prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) - if debug;fprintf('\tPropellant:\t\t%8s\n',prop_name);end - prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters - k = prop_params(1,1); % specific heat ratio - mw = prop_params(1,2); % molecular weight - if debug;fprintf('\tk:\t\t%16g\n\tmw:\t\t%16g\n',k,mw);end - total_params= fscanf(fid,'%g',[1 2]); % scan total/stagnation parameters - T_0 = total_params(1,1); % total temperature - P_0 = total_params(1,2); % total pressure - if debug;fprintf('\tT_0:\t%16g\n\tP_0:\t%16g\n',T_0,P_0);end - geom_size = fscanf(fid,'%g',[1 1]); % number of geometry nodes - xcoord = zeros(geom_size,1); radius = zeros(geom_size,1); - for i=1:geom_size - geom = fscanf(fid,'%g',[1 2]); - xcoord(i) = geom(1,1); % x coordinate of geometry node - radius(i) = geom(1,2); % radius at xcoord - end - if debug - fprintf('\tinlet radius:\t%8f\n\tthroat radius:\t%8f\n\texit radius:\t%8f\n',radius(1),min(radius),radius(end)); - fprintf('\tlength:\t%16f\n\tgeometry nodes:\t%8i\n',xcoord(end),geom_size); - end - fclose('all'); %close input file - if debug;fprintf('input file closed.\n');end - % MANUAL ENTRY - mw_units = '[kg/mol]'; - temperature_units = '[K]'; - pressure_units = '[Pa]'; - length_units = '[m]'; - angle_units = '[deg]'; +% IMPORT FROM FILE +fid = fopen(filein,'r'); +if debug;fprintf('reading data from input file (%s)...\n',filein);end +prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) + if debug;fprintf('\tPropellant:\t\t%8s\n',prop_name);end +prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters + k = prop_params(1,1); % specific heat ratio + mw = prop_params(1,2); % molecular weight + if debug;fprintf('\tk:\t\t%16g\n\tmw:\t\t%16g\n',k,mw);end +total_params= fscanf(fid,'%g',[1 3]); % scan total/stagnation parameters + T_0 = total_params(1,1); % total temperature + P_0 = total_params(1,2); % total pressure + P_b = total_params(1,3); % back pressure + if debug;fprintf('\tT_0:\t%16g\n\tP_0:\t%16g\n',T_0,P_0);end +geom_size = fscanf(fid,'%g',[1 1]); % number of geometry nodes +xcoord = zeros(geom_size,1); radius = zeros(geom_size,1); + for i=1:geom_size + geom = fscanf(fid,'%g',[1 2]); + xcoord(i) = geom(1,1); % x coordinate of geometry node + radius(i) = geom(1,2); % radius at xcoord + end + if debug + fprintf('\tinlet radius:\t%8f\n\tthroat radius:\t%8f\n\texit radius:\t%8f\n',radius(1),min(radius),radius(end)); + fprintf('\tlength:\t%16f\n\tgeometry nodes:\t%8i\n',xcoord(end),geom_size); + end +fclose('all'); %close input file +if debug;fprintf('input file closed.\n');end +% MANUAL ENTRY +mw_units = '[kg/mol]'; +temperature_units = '[K]'; +pressure_units = '[Pa]'; +length_units = '[m]'; +angle_units = '[deg]'; % ===------------=== % NOZZLE NOMENCLATURE @@ -62,117 +63,141 @@ % - Isentropic flow % - Initial temperature and pressure are total parameters - R = R_0/mw; % gas constant - R_units = '[J/kg K]'; +R = R_0/mw; % gas constant +R_units = '[J/kg K]'; % A_e = pi*exit_radius^2; % exit area % A_t = pi*throat_radius^2; % throat area - A=pi.*radius.^2; - area_units = '[m^2]'; - - A_t=min(A); - T_star=T_0*(2/(k+1)); %K - P_star=(2/(k+1))^(k/(k-1)); %Pa - rho_star=P_star/(R*T_star); % kg/m^3 - mdot=rho_star*sqrt(k*R*T_star)*A_t; %choked - - M_idx=linspace(.1,4,length(xcoord)); - k1=(2/(k+1)); - k2=((k-1)/2); - k3=(.5*(k+1)/(k-1)); - area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); +A=pi.*radius.^2; +area_units = '[m^2]'; - M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; - M(1)=0; - choked=false; - for x=2:length(xcoord) - if A(x)=1 - choked=true - else - disp('flag') - end - else - M_sup(x)=arearatio2mach_sup(A(x),A_t,k); - M_sub(x)=arearatio2mach_sub(A(x),A_t,k); - end - temp_ratio(x)=(1+((k-1)/2)*M(x)^2); - T(x)=T_0/temp_ratio(x); - pres_ratio(x)=temp_ratio(x)^(k/(k-1)); - P(x)=P_0/pres_ratio(x); - - temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); - T_sub(x)=T_0/temp_ratio_sub(x); - pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); - P_sub(x)=P_0/pres_ratio_sub(x); - temp_ratio_sup(x)=(1+((k-1)/2)*M_sup(x)^2); - T_sup(x)=T_0/temp_ratio_sup(x); - pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); - P_sup(x)=P_0/pres_ratio_sup(x); - end - if debug - figure - numplots=4;plotcounter=1; - subplot(numplots,1,plotcounter) - plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +[A_t,A_t_idx]=min(A); +T_star=T_0*(2/(k+1)); %K +P_star=(2/(k+1))^(k/(k-1)); %Pa +rho_star=P_star/(R*T_star); % kg/m^3 +mdot=rho_star*sqrt(k*R*T_star)*A_t; %choked + +M_idx=linspace(.1,4,length(xcoord)); +k1=(2/(k+1)); +k2=((k-1)/2); +k3=(.5*(k+1)/(k-1)); +area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); + +M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; +M(1)=1; +choked=false; +if debug;fprintf('Computing flow conditions along x-axis...\n');end +for x=2:length(xcoord) + M_sub(x)=arearatio2mach_sub(A(x),A_t,k); + if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end; + if M(x)<1 + M(x)=M_sub(x); + elseif M(x)==1 + choked=true; + if debug;fprintf('\tFlow is choked at A = %f\n',A(x));end + M(x)=M_sub(x); + elseif M(x)>1 + M(x)=M_sup(x); + end + temp_ratio(x)=(1+((k-1)/2)*M(x)^2); + T(x)=T_0/temp_ratio(x); + pres_ratio(x)=temp_ratio(x)^(k/(k-1)); + P(x)=P_0/pres_ratio(x); + + temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); + T_sub(x)=T_0/temp_ratio_sub(x); + pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); + P_sub(x)=P_0/pres_ratio_sub(x); + temp_ratio_sup(x)=(1+((k-1)/2)*M_sup(x)^2); + T_sup(x)=T_0/temp_ratio_sup(x); + pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); + P_sup(x)=P_0/pres_ratio_sup(x); + +end +if debug;fprintf('done.\n');end +if debug + figure + numplots=4;plotcounter=1; + subplot(numplots,1,plotcounter) + plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % subplot(numplots,1,plotcounter) % semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; - subplot(numplots,1,plotcounter) + subplot(numplots,1,plotcounter) % plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - semilogy(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--');legend('M','M<1','M>1');ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + subplot(numplots,1,plotcounter) + plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end));ylabel('T (K)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % semilogy(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); + subplot(numplots,1,plotcounter) + plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--');legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % semilogy(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % figure % semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); % figure % semilogy(M_idx,area_ratio,M,A./A_t,'--') - end + +figure + plot(xcoord,radius);ylabel('radius');axis([0 xcoord(end) 0 inf]); +figure + semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--'); + legend('M','M<1','M>1');ylabel('M');axis([0 xcoord(end) 0 inf]); +figure + plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); + ylabel('T (K)');axis([0 xcoord(end) 0 inf]); +figure + plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); + legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); +end % % format & display outputs - linedivider='------------'; - result = {'Propellant','',prop_name; - linedivider,'',''; - 'Specific heat ratio', k, unitless; - 'Molar mass', mw, mw_units; - 'Specific gas constant',R,R_units; - linedivider,'',''; - 'Total temperature', T_0, temperature_units; - 'Total pressure', P_0, pressure_units; - linedivider,'',''; - 'Length',xcoord(end),length_units; - 'Inlet radius',radius(1),length_units; - 'Throat radius', min(radius), length_units; - 'Exit radius', radius(end), length_units; - 'Inlet area',A(1),area_units; - 'Throat area',A_t,area_units; - 'Exit area',A(end),area_units; +linedivider='------------'; +result = {'Propellant','',prop_name; + linedivider,'',''; + 'Specific heat ratio', k, unitless; + 'Molar mass', mw, mw_units; + 'Specific gas constant',R,R_units; + linedivider,'',''; + 'Total temperature', T_0, temperature_units; + 'Total pressure', P_0, pressure_units; + linedivider,'',''; + 'Length',xcoord(end),length_units; + 'Inlet radius',radius(1),length_units; + 'Throat radius', min(radius), length_units; + 'Exit radius', radius(end), length_units; + 'Inlet area',A(1),area_units; + 'Throat area',A_t,area_units; + 'Exit area',A(end),area_units; % 'Half-angle',alpha, angle_units; - linedivider,'',''; - 'Throat temperature',T(find(A==A_t)),temperature_units; - 'Throat pressure',P(find(A==A_t)),pressure_units; + linedivider,'',''; + 'Throat temperature',T(A_t_idx),temperature_units; + 'Throat pressure',P(A_t_idx),pressure_units; % 'Mass flow rate',mass_flowrate,mass_flowrate_units; - linedivider,'',''; - 'Exit temperature',T(end),temperature_units; - 'Exit pressure',P(end),pressure_units; - linedivider,'',''; + linedivider,'',''; + 'Exit temperature',T(end),temperature_units; + 'Exit pressure',P(end),pressure_units; + linedivider,'',''; % 'Exhaust velocity',exit_velocity,velocity_units; % 'Thrust',thrust,force_units; % 'Specific impulse',specific_impulse,isp_units; - linedivider,'',''; - 'Exit Mach',M(end),unitless; - 'A/At',A(end)/A_t,unitless; - 'T/T0',T(end)/T_0,unitless; - 'P/P0',P(end)/P_0,unitless; + linedivider,'',''; + 'Exit Mach',M(end),unitless; + 'A/At',A(end)/A_t,unitless; + 'T/T0',T(end)/T_0,unitless; + 'P/P0',P(end)/P_0,unitless; % 'v/at',exit_velocity/throat_velocity,unitless; - }; - display(result); - -end - + }; +display(result); +figure + g=1.4; + g1=2/(g+1); + g2=(g-1)/2; + g3=0.5*(g+1)/(g-1); + M=0.1:0.1:4; + M=linspace(0.1,4); + A_theoretical=(1./M).*(g1*(1+g2*M.^2)).^g3; + A_trunc = round(A_theoretical.*1000)./1000; + semilogx(M,A_theoretical,M(find(A_trunc==1)),1,'o','linewidth',2) + grid on + xlabel('Mach Number') + ylabel('Area Ratio') function Mach = arearatio2mach_sub(A,A_t,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume % subsonic) @@ -184,7 +209,7 @@ r = (R-1)/(2*a); X = 1/((1+r)+sqrt(r*(r+2))); % initial guess diff = 1; % initalize termination criteria - while abs(diff) > .0001 + while abs(diff) > .00001 F = (P*X+Q).^(1/P)-R*X; dF = (P*X+Q).^((1/P)-1)-R; Xnew = X - F/dF; @@ -204,7 +229,7 @@ r = (R-1)/(2*a); X = 1/((1+r)+sqrt(r*(r+2))); % initial guess diff = 1; % initalize termination criteria - while abs(diff) > .0001 + while abs(diff) > .00001 F = (P*X+Q).^(1/P)-R*X; dF = (P*X+Q).^((1/P)-1)-R; Xnew = X - F/dF; From 3dc14c90b8c3ffc7e93473b8abc77b8a4d7af6d8 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 13:48:17 -0400 Subject: [PATCH 09/21] extra plots --- Matlab/tdu.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 0007af1..8c90c0e 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -88,7 +88,7 @@ if debug;fprintf('Computing flow conditions along x-axis...\n');end for x=2:length(xcoord) M_sub(x)=arearatio2mach_sub(A(x),A_t,k); - if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end; + if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end if M(x)<1 M(x)=M_sub(x); elseif M(x)==1 @@ -144,7 +144,7 @@ plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); ylabel('T (K)');axis([0 xcoord(end) 0 inf]); figure - plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); + semilogy(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); end % % format & display outputs From 559d8051d2e031d8e27e039904e1c8ff3d4b2a76 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 15:24:56 -0400 Subject: [PATCH 10/21] area vs m plot --- Matlab/tdu.asv | 41 +++++++++++++++++++++++++---------------- Matlab/tdu.m | 38 ++++++++++++++++++++++++-------------- 2 files changed, 49 insertions(+), 30 deletions(-) diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv index 841a704..3c48d0e 100644 --- a/Matlab/tdu.asv +++ b/Matlab/tdu.asv @@ -60,8 +60,11 @@ angle_units = '[deg]'; % ******************************** % % ASSUMPTIONS -% - Isentropic flow +% - Uniform flow across nozzle cross sections +% - Isentropic flow (except when normal shock occurs) % - Initial temperature and pressure are total parameters +% - Inlet flow is unmoving +% R = R_0/mw; % gas constant R_units = '[J/kg K]'; @@ -88,7 +91,7 @@ choked=false; if debug;fprintf('Computing flow conditions along x-axis...\n');end for x=2:length(xcoord) M_sub(x)=arearatio2mach_sub(A(x),A_t,k); - if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end; + if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end if M(x)<1 M(x)=M_sub(x); elseif M(x)==1 @@ -144,7 +147,7 @@ figure plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); ylabel('T (K)');axis([0 xcoord(end) 0 inf]); figure - plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); + semilogy(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); end % % format & display outputs @@ -185,19 +188,25 @@ result = {'Propellant','',prop_name; % 'v/at',exit_velocity/throat_velocity,unitless; }; display(result); -figure - g=1.4; - g1=2/(g+1); - g2=(g-1)/2; - g3=0.5*(g+1)/(g-1); - M=0.1:0.1:4; - M=linspace(0.1,4); - A_theoretical=(1./M).*(g1*(1+g2*M.^2)).^g3; - A_trunc = round(A_theoretical.*1000)./1000; - semilogx(M,A_theoretical,M(find(A_trunc==1)),,'o','linewidth',2) - grid on - xlabel('Mach Number') - ylabel('Area Ratio') +if debug + figure + g=1.4; + g1=2/(g+1); + g2=(g-1)/2; + g3=0.5*(g+1)/(g-1); + M_span=.1:.01:12; + Aratio =((1./M_span).*(g1*(1+g2*M_span.^2)).^g3); + Aratio_sonic =((1./1).*(g1*(1+g2*1.^2)).^g3); + A_star = A/Aratio_sonic; + A_theoretical=((1./M_span).*(g1*(1+g2*M_span.^2)).^g3).*A_star; + A_inlet = ones(length(M_span))*A(1)./A_t; + A_exit = ones(length(M_span))*A(end)./A_t; + A_throat = ones(length(M_span))*A_t; + loglog(M_span,A_theoretical,M_span,A_inlet,'--',M_span,A_exit,'--',M_span,A_throat,'--',M_sub,A,'k',M_sup,A,'k') + grid on + xlabel('Mach Number') + ylabel('Area') +end function Mach = arearatio2mach_sub(A,A_t,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume % subsonic) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 8c90c0e..3d9f5ea 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -60,8 +60,11 @@ % ******************************** % % ASSUMPTIONS -% - Isentropic flow +% - Uniform flow across nozzle cross sections +% - Isentropic flow (except when normal shock occurs) % - Initial temperature and pressure are total parameters +% - Inlet flow is unmoving +% R = R_0/mw; % gas constant R_units = '[J/kg K]'; @@ -185,19 +188,26 @@ % 'v/at',exit_velocity/throat_velocity,unitless; }; display(result); -figure - g=1.4; - g1=2/(g+1); - g2=(g-1)/2; - g3=0.5*(g+1)/(g-1); - M=0.1:0.1:4; - M=linspace(0.1,4); - A_theoretical=(1./M).*(g1*(1+g2*M.^2)).^g3; - A_trunc = round(A_theoretical.*1000)./1000; - semilogx(M,A_theoretical,M(find(A_trunc==1)),1,'o','linewidth',2) - grid on - xlabel('Mach Number') - ylabel('Area Ratio') +if debug + figure + g=1.4; + g1=2/(g+1); + g2=(g-1)/2; + g3=0.5*(g+1)/(g-1); + M_span=.0015:.001:7.5; + Aratio =((1./M_span).*(g1*(1+g2*M_span.^2)).^g3); + Aratio_sonic =((1./1).*(g1*(1+g2*1.^2)).^g3); + A_star = A_t/Aratio_sonic; + A_theoretical=((1./M_span).*(g1*(1+g2*M_span.^2)).^g3).*A_star; + A_inlet = ones(length(M_span))*A(1); + A_exit = ones(length(M_span))*A(end); + A_throat = ones(length(M_span))*A_t; + loglog(M_span,A_theoretical,M_span,A_inlet,'--',M_span,A_exit,'--',M_span,A_throat,'--',M_sub,A,'k',M_sup,A,'k') + legend('theoretical area from F_i_s(M)','inlet area','exit area','throat area','area as function of M found numerically',location,'southeastoutside') + grid on + xlabel('Mach') + ylabel('Area (m^2)') +end function Mach = arearatio2mach_sub(A,A_t,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume % subsonic) From 9bd4a33f23b6d54915cf80003ba2ebc2d4e04b2d Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 15:27:08 -0400 Subject: [PATCH 11/21] area ratio vs m plot with debug notes --- Matlab/tdu.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 3d9f5ea..a6583e5 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -189,6 +189,7 @@ }; display(result); if debug + fprintf('plotting area ratio...'); figure g=1.4; g1=2/(g+1); @@ -203,10 +204,11 @@ A_exit = ones(length(M_span))*A(end); A_throat = ones(length(M_span))*A_t; loglog(M_span,A_theoretical,M_span,A_inlet,'--',M_span,A_exit,'--',M_span,A_throat,'--',M_sub,A,'k',M_sup,A,'k') - legend('theoretical area from F_i_s(M)','inlet area','exit area','throat area','area as function of M found numerically',location,'southeastoutside') + legend('theoretical area from F_i_s(M)','inlet area','exit area','throat area','area as function of M found numerically','location','southeastoutside') grid on xlabel('Mach') ylabel('Area (m^2)') + fprintf('done'); end function Mach = arearatio2mach_sub(A,A_t,k) % solve Mach number from area ratio by Newton-Raphson Method. (assume From 295ad088396e55f98dfa0b369493d6d03148600f Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 15:31:07 -0400 Subject: [PATCH 12/21] tweak --- Matlab/tdu.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index a6583e5..c0bdd2c 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -204,7 +204,7 @@ A_exit = ones(length(M_span))*A(end); A_throat = ones(length(M_span))*A_t; loglog(M_span,A_theoretical,M_span,A_inlet,'--',M_span,A_exit,'--',M_span,A_throat,'--',M_sub,A,'k',M_sup,A,'k') - legend('theoretical area from F_i_s(M)','inlet area','exit area','throat area','area as function of M found numerically','location','southeastoutside') + legend('theoretical area from F_i_s(M)','inlet area','exit area','throat area','area as function of M found numerically','location','eastoutside') grid on xlabel('Mach') ylabel('Area (m^2)') From 766c055aca91e5449a98f76d3873767a83518d30 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Fri, 5 May 2017 15:33:19 -0400 Subject: [PATCH 13/21] coarser area ratio plot --- Matlab/tdu.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index c0bdd2c..fc5bedf 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -195,7 +195,7 @@ g1=2/(g+1); g2=(g-1)/2; g3=0.5*(g+1)/(g-1); - M_span=.0015:.001:7.5; + M_span=.0015:.01:7.5; Aratio =((1./M_span).*(g1*(1+g2*M_span.^2)).^g3); Aratio_sonic =((1./1).*(g1*(1+g2*1.^2)).^g3); A_star = A_t/Aratio_sonic; From dfa93122341c3ccdf71dc90b99fda339f52d2917 Mon Sep 17 00:00:00 2001 From: runphilrun Date: Sat, 6 May 2017 00:41:08 -0400 Subject: [PATCH 14/21] remove newton raphson --- Matlab/tdu.m | 91 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 36 deletions(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index fc5bedf..43c1ca2 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -210,45 +210,64 @@ ylabel('Area (m^2)') fprintf('done'); end + function Mach = arearatio2mach_sub(A,A_t,k) -% solve Mach number from area ratio by Newton-Raphson Method. (assume -% subsonic) -% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html - P = 2/(k+1); - Q = 1-P; - R = (A/A_t).^2; - a = P.^(1/Q); - r = (R-1)/(2*a); - X = 1/((1+r)+sqrt(r*(r+2))); % initial guess - diff = 1; % initalize termination criteria - while abs(diff) > .00001 - F = (P*X+Q).^(1/P)-R*X; - dF = (P*X+Q).^((1/P)-1)-R; - Xnew = X - F/dF; - diff = Xnew - X; - X = Xnew; - end - Mach = sqrt(X); +% % solve Mach number from area ratio by Newton-Raphson Method. (assume +% % subsonic) +% % https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html +% P = 2/(k+1); +% Q = 1-P; +% R = (A/A_t).^2; +% a = P.^(1/Q); +% r = (R-1)/(2*a); +% X = 1/((1+r)+sqrt(r*(r+2))); % initial guess +% diff = 1; % initalize termination criteria +% while abs(diff) > .00001 +% F = (P*X+Q).^(1/P)-R*X; +% dF = (P*X+Q).^((1/P)-1)-R; +% Xnew = X - F/dF; +% diff = Xnew - X; +% X = Xnew; +% end +% Mach = sqrt(X); +M=.001:.001:1; +g=k; +g1=2/(g+1); +g2=(g-1)/2; +g3=0.5*(g+1)/(g-1); +arearatio=1./M.*(g1*(1+g2*M.^2)).^g3; +diff=abs(arearatio-A/A_t); +[res,i]=min(diff); +Mach=M(i); end function Mach = arearatio2mach_sup(A,A_t,k) -% solve Mach number from area ratio by Newton-Raphson Method. (assume -% supersonic) -% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html - P = 2/(k+1); - Q = 1-P; - R = (A/A_t).^((2*Q)/P); - a = Q.^(1/P); - r = (R-1)/(2*a); - X = 1/((1+r)+sqrt(r*(r+2))); % initial guess - diff = 1; % initalize termination criteria - while abs(diff) > .00001 - F = (P*X+Q).^(1/P)-R*X; - dF = (P*X+Q).^((1/P)-1)-R; - Xnew = X - F/dF; - diff = Xnew - X; - X = Xnew; - end - Mach = 1/sqrt(X); +% % solve Mach number from area ratio by Newton-Raphson Method. (assume +% % supersonic) +% % https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html +% P = 2/(k+1); +% Q = 1-P; +% R = (A/A_t).^((2*Q)/P); +% a = Q.^(1/P); +% r = (R-1)/(2*a); +% X = 1/((1+r)+sqrt(r*(r+2))); % initial guess +% diff = 1; % initalize termination criteria +% while abs(diff) > .00001 +% F = (P*X+Q).^(1/P)-R*X; +% dF = (P*X+Q).^((1/P)-1)-R; +% Xnew = X - F/dF; +% diff = Xnew - X; +% X = Xnew; +% end +% Mach = 1/sqrt(X); + M=1:.01:10; + g=k; + g1=2/(g+1); + g2=(g-1)/2; + g3=0.5*(g+1)/(g-1); + arearatio=1./M.*(g1*(1+g2*M.^2)).^g3; + diff=abs(arearatio-A/A_t); + [res,i]=min(diff); + Mach=M(i); end function display(result) global debug From 7cabbff4ab169caeb0972b7b06232db03339c18d Mon Sep 17 00:00:00 2001 From: runphilrun Date: Sat, 6 May 2017 01:32:08 -0400 Subject: [PATCH 15/21] m from area ratio using graph --- .gitignore | 2 +- Matlab/tdu.m | 71 +++++++++++++++++++++++++++------------------------- 2 files changed, 38 insertions(+), 35 deletions(-) diff --git a/.gitignore b/.gitignore index e1d5b53..d7c7dc8 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,4 @@ dist/ .doctrees/ *.egg-info/ __pycache__/ -.asv/ +*.asv/ diff --git a/Matlab/tdu.m b/Matlab/tdu.m index 43c1ca2..ac3a0f5 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -83,22 +83,23 @@ k1=(2/(k+1)); k2=((k-1)/2); k3=(.5*(k+1)/(k-1)); -area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); +% area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; -M(1)=1; +M(1)=0; choked=false; +shocks=false; if debug;fprintf('Computing flow conditions along x-axis...\n');end for x=2:length(xcoord) M_sub(x)=arearatio2mach_sub(A(x),A_t,k); - if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end - if M(x)<1 + M_sup(x)=arearatio2mach_sup(A(x),A_t,k); + if M(x-1)<1 M(x)=M_sub(x); - elseif M(x)==1 + elseif M(x-1)==1 choked=true; if debug;fprintf('\tFlow is choked at A = %f\n',A(x));end - M(x)=M_sub(x); - elseif M(x)>1 + M(x)=M_sup(x); + elseif M(x-1)>1 M(x)=M_sup(x); end temp_ratio(x)=(1+((k-1)/2)*M(x)^2); @@ -114,41 +115,43 @@ T_sup(x)=T_0/temp_ratio_sup(x); pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); P_sup(x)=P_0/pres_ratio_sup(x); - end +if debug&¬(choked);fprintf('\tFlow is NOT choked.\n');end if debug;fprintf('done.\n');end -if debug - figure - numplots=4;plotcounter=1; - subplot(numplots,1,plotcounter) - plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +figure +numplots=4;plotcounter=1; +subplot(numplots,1,plotcounter) +plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % subplot(numplots,1,plotcounter) % semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; - subplot(numplots,1,plotcounter) +subplot(numplots,1,plotcounter) % plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--');legend('M','M<1','M>1');ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end));ylabel('T (K)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +semilogy(xcoord,M,'k',xcoord,M_sub,xcoord(A_t_idx:end),M_sup(A_t_idx:end)); + legend('M','M<1','M>1');ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +subplot(numplots,1,plotcounter) +plot(xcoord,T,'k',xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); + legend('T','T_M_<_1','T_M_>_1');ylabel('T (K)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % semilogy(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--');legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +subplot(numplots,1,plotcounter) +plot(xcoord,P,'k',xcoord,P_sub,xcoord(A_t_idx:end),P_sup(A_t_idx:end)); + legend('P','P_M_<_1','P_M_>_1');ylabel('P (Pa)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % semilogy(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % figure % semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); % figure % semilogy(M_idx,area_ratio,M,A./A_t,'--') - -figure - plot(xcoord,radius);ylabel('radius');axis([0 xcoord(end) 0 inf]); -figure - semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--'); - legend('M','M<1','M>1');ylabel('M');axis([0 xcoord(end) 0 inf]); -figure - plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); - ylabel('T (K)');axis([0 xcoord(end) 0 inf]); -figure - semilogy(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); - legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); +if debug + figure + plot(xcoord,radius);ylabel('radius');axis([0 xcoord(end) 0 inf]); + figure + semilogy(xcoord,M,'k',xcoord,M_sub,xcoord(A_t_idx:end),M_sup(A_t_idx:end)); + legend('M','M<1','M>1');ylabel('M');axis([0 xcoord(end) 0 inf]); + figure + plot(xcoord,T,'k',xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); + legend('T','T_M_<_1','T_M_>_1');ylabel('T (K)');axis([0 xcoord(end) 0 inf]); + figure + semilogy(xcoord,P,'k',xcoord,P_sub,xcoord(A_t_idx:end),P_sup(A_t_idx:end),xcoord,ones(length(xcoord))*P_b,'k--'); + legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); end % % format & display outputs linedivider='------------'; @@ -208,7 +211,7 @@ grid on xlabel('Mach') ylabel('Area (m^2)') - fprintf('done'); + fprintf('done.\n'); end function Mach = arearatio2mach_sub(A,A_t,k) @@ -230,7 +233,7 @@ % X = Xnew; % end % Mach = sqrt(X); -M=.001:.001:1; +M=.001:.0005:1; g=k; g1=2/(g+1); g2=(g-1)/2; @@ -259,7 +262,7 @@ % X = Xnew; % end % Mach = 1/sqrt(X); - M=1:.01:10; + M=1:.005:12; g=k; g1=2/(g+1); g2=(g-1)/2; From 6d763868a615ed7356468d14450133bce3d5199a Mon Sep 17 00:00:00 2001 From: runphilrun Date: Sat, 6 May 2017 01:44:01 -0400 Subject: [PATCH 16/21] set up shock calcs --- Matlab/tdu.m | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/Matlab/tdu.m b/Matlab/tdu.m index ac3a0f5..f500f18 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -85,11 +85,27 @@ k3=(.5*(k+1)/(k-1)); % area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); +if debug;fprintf('Checking for normal shocks...\n');end +M_sub=arearatio2mach_sub(A(end),A_t,k); +M_sup=arearatio2mach_sup(A(end),A_t,k); +pres_ratio_sub=(1+((k-1)/2)*M_sub^2)^(k/(k-1)); +pres_ratio_sup=(1+((k-1)/2)*M_sup^2)^(k/(k-1)); +Pe1=P_0/pres_ratio_sub; +Pe2=P_0/pres_ratio_sup; +if P_b >= Pe1 + if debug;fprintf('\tAlways subsonic.\n');end +elseif P_b < Pe1 && P_b > Pe2 + if debug;fprintf('\tNormal shock exists.\n');end +elseif P_b == Pe2 + if debug;fprintf('\tIsentropically supersonic through entire nozzle.\n');end +else + if debug;fprintf('\tERROR!\n');end +end + +if debug;fprintf('Computing flow conditions along x-axis...\n');end M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; M(1)=0; choked=false; -shocks=false; -if debug;fprintf('Computing flow conditions along x-axis...\n');end for x=2:length(xcoord) M_sub(x)=arearatio2mach_sub(A(x),A_t,k); M_sup(x)=arearatio2mach_sup(A(x),A_t,k); From a43b86b2e49aa367e3df27ff68c46e13e09200c1 Mon Sep 17 00:00:00 2001 From: runphilrun Date: Sun, 7 May 2017 10:54:08 -0400 Subject: [PATCH 17/21] add shock calcs --- .gitignore | 1 + Matlab/tdu.m | 62 +++++++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 58 insertions(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index d7c7dc8..8e56b46 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ dist/ *.egg-info/ __pycache__/ *.asv/ +*.zip diff --git a/Matlab/tdu.m b/Matlab/tdu.m index f500f18..afa4444 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -45,6 +45,7 @@ pressure_units = '[Pa]'; length_units = '[m]'; angle_units = '[deg]'; +mdot_units = '[kg/s]'; % ===------------=== % NOZZLE NOMENCLATURE @@ -96,6 +97,9 @@ if debug;fprintf('\tAlways subsonic.\n');end elseif P_b < Pe1 && P_b > Pe2 if debug;fprintf('\tNormal shock exists.\n');end + [ashock, Mminus]=ShockPosition2(A(end),A_t,k,P_0,P_b); +% xshock=A(find(A==round(ashock*1000)/1000)) +% xshock2=A(find(M==Mminus)) elseif P_b == Pe2 if debug;fprintf('\tIsentropically supersonic through entire nozzle.\n');end else @@ -121,8 +125,8 @@ temp_ratio(x)=(1+((k-1)/2)*M(x)^2); T(x)=T_0/temp_ratio(x); pres_ratio(x)=temp_ratio(x)^(k/(k-1)); - P(x)=P_0/pres_ratio(x); - + P(x)=P_0/pres_ratio(x); + temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); T_sub(x)=T_0/temp_ratio_sub(x); pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); @@ -132,6 +136,8 @@ pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); P_sup(x)=P_0/pres_ratio_sup(x); end +thrust = mdot*M(end)*sqrt(k*R*T(end))-(P(end)-P_b)*A(end); +force_units='[N]'; if debug&¬(choked);fprintf('\tFlow is NOT choked.\n');end if debug;fprintf('done.\n');end figure @@ -191,13 +197,13 @@ linedivider,'',''; 'Throat temperature',T(A_t_idx),temperature_units; 'Throat pressure',P(A_t_idx),pressure_units; -% 'Mass flow rate',mass_flowrate,mass_flowrate_units; + 'Mass flow rate',mdot,mdot_units; linedivider,'',''; 'Exit temperature',T(end),temperature_units; 'Exit pressure',P(end),pressure_units; linedivider,'',''; % 'Exhaust velocity',exit_velocity,velocity_units; -% 'Thrust',thrust,force_units; + 'Thrust',thrust,force_units; % 'Specific impulse',specific_impulse,isp_units; linedivider,'',''; 'Exit Mach',M(end),unitless; @@ -288,12 +294,58 @@ [res,i]=min(diff); Mach=M(i); end +function [shockArea,Mminus]=ShockPosition2(Ae,At,g,P0,Pe) +global debug +% Ae=3; %exit area +% At=1; %Throat area +%Me=0.4 %Exit Mach +% P0=1; %Total P at the inlet +% Pe=0.6; %Static P at exit +% g=1.4 %isentropic constant +g1=(g-1)/2; +g2=2/(g+1); +g3=0.5*(g+1)/(g-1); + +FUN=@(M) 1+g1*M.^2; +FA=@(M) (1./M).*(g2*FUN(M)).^g3; %isentropic Area/Areat +MSHOCK=@(M) FUN(M)./(g*M.^2-g1); %shock for M +Pis=@(M) FUN(M).^(g/g-1); +PSHOCK=@(M) 1+2*g/(g+1)*(M.^2-1); + +%Arrays of shock areas satisfying upstream and downstream conditions +MM=1:0.01:10; +MMPLUS=MSHOCK(MM); +ASH1=At*FA(MM); + +P0SHOCK=(PSHOCK(MM).*Pis(MMPLUS))./Pis(MM); +P02=P0*P0SHOCK; +Me=2/(g-1)*((P02/Pe).^((g-1)/g)-1); Me=sqrt(Me); +ASH2=Ae*FA(Me)./FA(MMPLUS); +if debug + figure + plot(MM,ASH1,'--r','LineWidth',3) + hold on + plot(MM,ASH2,'--g','LineWidth',3) + grid on + xlabel('Mach Number') + ylabel('A - Shock Area') + title('Shock Cross Section Area') +end + +%Solution as a min of residual +[RES,i]=min((ASH1-ASH2).^2); +% RES +if debug + shockArea=ASH1(i) + Mminus=MM(i) +end +end function display(result) global debug if debug;fprintf('displaying results...\n');end; [n,~]=size(result); for i = 1:n - fprintf('\n%24s\t%15.8f\t%s',result{i,:}); + fprintf('\n%24s\t%15.4g\t%s',result{i,:}); end fprintf('\n') if debug;fprintf('done.\n');end; From edb7d0f73f9fb84d7ab1ebb923e759d800909544 Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Tue, 9 May 2017 14:00:15 -0400 Subject: [PATCH 18/21] update for v1.0 release -notdone --- README.md | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index 4c52126..b79af2a 100644 --- a/README.md +++ b/README.md @@ -10,8 +10,8 @@ > Tool to aid in design of small monopropellant thrusters. *[View this project on STEMN](http://stemn.com/projects/thruster-design-tool)* -## Abstract -The thruster's nozzle, propellant, and chamber conditions each have a huge impact on performance. The fluid mechanics that model engine effectiveness are pretty complicated, unfortunately. The purpose of this script is to make it easy to see how tweaks in nozzle geometry, propellants, and chamber conditions affect performance in order to find the optimal design solution. +## Scope +The purpose of this script is to make it easy to see how tweaks in nozzle geometry, propellants, chamber conditions, and ambient conditions affect nozzle performance in order to find the optimal design solution. This script does not perform any optimization. ## Features * Simulate the performance of a thruster in space for a given set of parameters and output performance metrics. @@ -19,22 +19,16 @@ The thruster's nozzle, propellant, and chamber conditions each have a huge impac * Allow the user to easily tweak parameters. ## Usage -### Matlab -Open `Matlab/tdu.m` in Matlab 2014 or newer. -Edit fields as indicated to specify the propellant gas properties and nozzle dimensions of the engine, then run the script. -Mach number, temperature ratio, and pressure ratio at the exit of the nozzle agree with [NASA Report 1135](http://www.nasa.gov/sites/default/files/734673main_Equations-Tables-Charts-CompressibleFlow-Report-1135.pdf) for air at 1 atm with an area ratio of 2.005. - -### Python +### Generating an input file +TDU loads propellant properties, inlet conditions, and nozzle geometry from a specially formatted tab-delimited text file with the extension `*.tdu`. -```bash -$ pip install tdu -$ tdu --help -``` -Set propellant gas properties, nozzle geometry and chamber conditions in `config.ini`. -Run `tdu.py`. +### Running the script +* Open `tdu.m` in Matlab 2014 or newer. +* Specify the desired input file as the value of `filein`. (For example, `filein='sample.tdu';`) +* To show additional plots and print verbose actions and data to the command line, set `debug=true;`. +* Run `tdu` +Mach number, temperature ratio, and pressure ratio at the exit of the nozzle agree with [NASA Report 1135](http://www.nasa.gov/sites/default/files/734673main_Equations-Tables-Charts-CompressibleFlow-Report-1135.pdf) for air at 1 atm with an area ratio of 2.005. -## About & Documentation -Please refer to the [STEMN project page](http://stemn.com/projects/thruster-design-tool) for detailed information about the project and general theory. > *If you encounter any bugs, please report them in the [Issue Tracker](https://github.com/runphilrun/ThrusterDesign/issues)!* From 80a92c789e11ef31574c8eac395ca6c61f58eaee Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Tue, 9 May 2017 14:11:15 -0400 Subject: [PATCH 19/21] format of input file --- README.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/README.md b/README.md index b79af2a..9d1e348 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,19 @@ The purpose of this script is to make it easy to see how tweaks in nozzle geomet ### Generating an input file TDU loads propellant properties, inlet conditions, and nozzle geometry from a specially formatted tab-delimited text file with the extension `*.tdu`. +In general, the format of an input file is as follows: +``` +PropellantNameString +Gamma MolecularWeight +T0 P0 +NumberofNodes +xLocation0 radius0 +xlocation1 radius1 +. . +. . +. . +xlocationN radiusN +``` ### Running the script * Open `tdu.m` in Matlab 2014 or newer. From 84d1995ff7714a55ee4196b8769e0434a1a9485e Mon Sep 17 00:00:00 2001 From: Philip Linden Date: Tue, 9 May 2017 14:12:18 -0400 Subject: [PATCH 20/21] format of input file --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 9d1e348..32dbe26 100644 --- a/README.md +++ b/README.md @@ -30,9 +30,9 @@ T0 P0 NumberofNodes xLocation0 radius0 xlocation1 radius1 -. . -. . -. . + . . + . . + . . xlocationN radiusN ``` From d0d7651002c72726f217178798d0a556379d88ad Mon Sep 17 00:00:00 2001 From: runphilrun Date: Wed, 10 May 2017 17:12:21 -0400 Subject: [PATCH 21/21] minor output tweaks --- Matlab/tdu.asv | 259 ------------------------------------------------- Matlab/tdu.m | 10 +- 2 files changed, 4 insertions(+), 265 deletions(-) delete mode 100644 Matlab/tdu.asv diff --git a/Matlab/tdu.asv b/Matlab/tdu.asv deleted file mode 100644 index 3c48d0e..0000000 --- a/Matlab/tdu.asv +++ /dev/null @@ -1,259 +0,0 @@ -% HEADER -clc; -close all -format long - -filein = 'msd-050.tdu'; % input file - -% universal constants -R_0 = 8.3144598; % [J/(mol*K)] universal gas constant -g_0 = 9.81; % [m/s^2] standard gravity -unitless='[-]'; -global debug; -debug = true; -% === THRUSTER PARAMETERS === -% IMPORT FROM FILE -fid = fopen(filein,'r'); -if debug;fprintf('reading data from input file (%s)...\n',filein);end -prop_name = fscanf(fid,'%s',[1,1]); % descriptive header (no quotes, no spaces) - if debug;fprintf('\tPropellant:\t\t%8s\n',prop_name);end -prop_params = fscanf(fid,'%g',[1 2]); % scan propellant parameters - k = prop_params(1,1); % specific heat ratio - mw = prop_params(1,2); % molecular weight - if debug;fprintf('\tk:\t\t%16g\n\tmw:\t\t%16g\n',k,mw);end -total_params= fscanf(fid,'%g',[1 3]); % scan total/stagnation parameters - T_0 = total_params(1,1); % total temperature - P_0 = total_params(1,2); % total pressure - P_b = total_params(1,3); % back pressure - if debug;fprintf('\tT_0:\t%16g\n\tP_0:\t%16g\n',T_0,P_0);end -geom_size = fscanf(fid,'%g',[1 1]); % number of geometry nodes -xcoord = zeros(geom_size,1); radius = zeros(geom_size,1); - for i=1:geom_size - geom = fscanf(fid,'%g',[1 2]); - xcoord(i) = geom(1,1); % x coordinate of geometry node - radius(i) = geom(1,2); % radius at xcoord - end - if debug - fprintf('\tinlet radius:\t%8f\n\tthroat radius:\t%8f\n\texit radius:\t%8f\n',radius(1),min(radius),radius(end)); - fprintf('\tlength:\t%16f\n\tgeometry nodes:\t%8i\n',xcoord(end),geom_size); - end -fclose('all'); %close input file -if debug;fprintf('input file closed.\n');end -% MANUAL ENTRY -mw_units = '[kg/mol]'; -temperature_units = '[K]'; -pressure_units = '[Pa]'; -length_units = '[m]'; -angle_units = '[deg]'; -% ===------------=== - -% NOZZLE NOMENCLATURE -% ******************************** -% -% /- -% / 0 = total parameters -% ===\---/ 1 = chamber -% (1) (2) (3) 2 = throat -% ===/---\ 3 = exit -% \ -% \- -% ******************************** - -% % ASSUMPTIONS -% - Uniform flow across nozzle cross sections -% - Isentropic flow (except when normal shock occurs) -% - Initial temperature and pressure are total parameters -% - Inlet flow is unmoving -% - -R = R_0/mw; % gas constant -R_units = '[J/kg K]'; -% A_e = pi*exit_radius^2; % exit area -% A_t = pi*throat_radius^2; % throat area -A=pi.*radius.^2; -area_units = '[m^2]'; - -[A_t,A_t_idx]=min(A); -T_star=T_0*(2/(k+1)); %K -P_star=(2/(k+1))^(k/(k-1)); %Pa -rho_star=P_star/(R*T_star); % kg/m^3 -mdot=rho_star*sqrt(k*R*T_star)*A_t; %choked - -M_idx=linspace(.1,4,length(xcoord)); -k1=(2/(k+1)); -k2=((k-1)/2); -k3=(.5*(k+1)/(k-1)); -area_ratio=(1./M_idx).*((2/(k+1))*(1+((k-1)/2)*M_idx.^2)).^(.5*(k+1)/(k-1)); - -M=zeros(length(xcoord),1);M_sub=M;M_sup=M;temp_ratio=M;T=M;pres_ratio=M;P=M; -M(1)=1; -choked=false; -if debug;fprintf('Computing flow conditions along x-axis...\n');end -for x=2:length(xcoord) - M_sub(x)=arearatio2mach_sub(A(x),A_t,k); - if x >=A_t_idx;M_sup(x)=arearatio2mach_sup(A(x),A_t,k);end - if M(x)<1 - M(x)=M_sub(x); - elseif M(x)==1 - choked=true; - if debug;fprintf('\tFlow is choked at A = %f\n',A(x));end - M(x)=M_sub(x); - elseif M(x)>1 - M(x)=M_sup(x); - end - temp_ratio(x)=(1+((k-1)/2)*M(x)^2); - T(x)=T_0/temp_ratio(x); - pres_ratio(x)=temp_ratio(x)^(k/(k-1)); - P(x)=P_0/pres_ratio(x); - - temp_ratio_sub(x)=(1+((k-1)/2)*M_sub(x)^2); - T_sub(x)=T_0/temp_ratio_sub(x); - pres_ratio_sub(x)=temp_ratio_sub(x)^(k/(k-1)); - P_sub(x)=P_0/pres_ratio_sub(x); - temp_ratio_sup(x)=(1+((k-1)/2)*M_sup(x)^2); - T_sup(x)=T_0/temp_ratio_sup(x); - pres_ratio_sup(x)=temp_ratio_sup(x)^(k/(k-1)); - P_sup(x)=P_0/pres_ratio_sup(x); - -end -if debug;fprintf('done.\n');end -if debug - figure - numplots=4;plotcounter=1; - subplot(numplots,1,plotcounter) - plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); -% subplot(numplots,1,plotcounter) -% semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; - subplot(numplots,1,plotcounter) -% plot(xcoord,M_sub,xcoord,M_sup);ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--');legend('M','M<1','M>1');ylabel('M');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end));ylabel('T (K)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); -% semilogy(xcoord,T_sub,xcoord,T_sup);ylabel('T');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); - subplot(numplots,1,plotcounter) - plot(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--');legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); -% semilogy(xcoord,P_sub/10^3,xcoord,P_sup/10^3);ylabel('P');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); -% figure -% semilogy(M(1:find(A==A_t)),A(1:find(A==A_t))./A_t,M(find(A==A_t)+1:end),A(find(A==A_t)+1:end)./A_t,'--');xlabel('M');ylabel('A/A_t'); -% figure -% semilogy(M_idx,area_ratio,M,A./A_t,'--') - -figure - plot(xcoord,radius);ylabel('radius');axis([0 xcoord(end) 0 inf]); -figure - semilogy(xcoord,M,'k',xcoord,M_sub,'--',xcoord(A_t_idx:end),M_sup(A_t_idx:end),'--'); - legend('M','M<1','M>1');ylabel('M');axis([0 xcoord(end) 0 inf]); -figure - plot(xcoord,T_sub,xcoord(A_t_idx:end),T_sup(A_t_idx:end)); - ylabel('T (K)');axis([0 xcoord(end) 0 inf]); -figure - semilogy(xcoord,P,'k',xcoord,P_sub,'--',xcoord(A_t_idx:end),P_sup(A_t_idx:end),'--',xcoord,ones(length(xcoord))*P_b,'k--'); - legend('P','P_M_<_1','P_M_>_1','P_a_m_b');ylabel('P (Pa)');axis([0 xcoord(end) 0 inf]); -end -% % format & display outputs -linedivider='------------'; -result = {'Propellant','',prop_name; - linedivider,'',''; - 'Specific heat ratio', k, unitless; - 'Molar mass', mw, mw_units; - 'Specific gas constant',R,R_units; - linedivider,'',''; - 'Total temperature', T_0, temperature_units; - 'Total pressure', P_0, pressure_units; - linedivider,'',''; - 'Length',xcoord(end),length_units; - 'Inlet radius',radius(1),length_units; - 'Throat radius', min(radius), length_units; - 'Exit radius', radius(end), length_units; - 'Inlet area',A(1),area_units; - 'Throat area',A_t,area_units; - 'Exit area',A(end),area_units; -% 'Half-angle',alpha, angle_units; - linedivider,'',''; - 'Throat temperature',T(A_t_idx),temperature_units; - 'Throat pressure',P(A_t_idx),pressure_units; -% 'Mass flow rate',mass_flowrate,mass_flowrate_units; - linedivider,'',''; - 'Exit temperature',T(end),temperature_units; - 'Exit pressure',P(end),pressure_units; - linedivider,'',''; -% 'Exhaust velocity',exit_velocity,velocity_units; -% 'Thrust',thrust,force_units; -% 'Specific impulse',specific_impulse,isp_units; - linedivider,'',''; - 'Exit Mach',M(end),unitless; - 'A/At',A(end)/A_t,unitless; - 'T/T0',T(end)/T_0,unitless; - 'P/P0',P(end)/P_0,unitless; -% 'v/at',exit_velocity/throat_velocity,unitless; - }; -display(result); -if debug - figure - g=1.4; - g1=2/(g+1); - g2=(g-1)/2; - g3=0.5*(g+1)/(g-1); - M_span=.1:.01:12; - Aratio =((1./M_span).*(g1*(1+g2*M_span.^2)).^g3); - Aratio_sonic =((1./1).*(g1*(1+g2*1.^2)).^g3); - A_star = A/Aratio_sonic; - A_theoretical=((1./M_span).*(g1*(1+g2*M_span.^2)).^g3).*A_star; - A_inlet = ones(length(M_span))*A(1)./A_t; - A_exit = ones(length(M_span))*A(end)./A_t; - A_throat = ones(length(M_span))*A_t; - loglog(M_span,A_theoretical,M_span,A_inlet,'--',M_span,A_exit,'--',M_span,A_throat,'--',M_sub,A,'k',M_sup,A,'k') - grid on - xlabel('Mach Number') - ylabel('Area') -end -function Mach = arearatio2mach_sub(A,A_t,k) -% solve Mach number from area ratio by Newton-Raphson Method. (assume -% subsonic) -% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html - P = 2/(k+1); - Q = 1-P; - R = (A/A_t).^2; - a = P.^(1/Q); - r = (R-1)/(2*a); - X = 1/((1+r)+sqrt(r*(r+2))); % initial guess - diff = 1; % initalize termination criteria - while abs(diff) > .00001 - F = (P*X+Q).^(1/P)-R*X; - dF = (P*X+Q).^((1/P)-1)-R; - Xnew = X - F/dF; - diff = Xnew - X; - X = Xnew; - end - Mach = sqrt(X); -end -function Mach = arearatio2mach_sup(A,A_t,k) -% solve Mach number from area ratio by Newton-Raphson Method. (assume -% supersonic) -% https://www.grc.nasa.gov/WWW/winddocs/utilities/b4wind_guide/mach.html - P = 2/(k+1); - Q = 1-P; - R = (A/A_t).^((2*Q)/P); - a = Q.^(1/P); - r = (R-1)/(2*a); - X = 1/((1+r)+sqrt(r*(r+2))); % initial guess - diff = 1; % initalize termination criteria - while abs(diff) > .00001 - F = (P*X+Q).^(1/P)-R*X; - dF = (P*X+Q).^((1/P)-1)-R; - Xnew = X - F/dF; - diff = Xnew - X; - X = Xnew; - end - Mach = 1/sqrt(X); -end -function display(result) - global debug - if debug;fprintf('displaying results...\n');end; - [n,~]=size(result); - for i = 1:n - fprintf('\n%24s\t%15.8f\t%s',result{i,:}); - end - fprintf('\n') - if debug;fprintf('done.\n');end; -end \ No newline at end of file diff --git a/Matlab/tdu.m b/Matlab/tdu.m index afa4444..19250c6 100644 --- a/Matlab/tdu.m +++ b/Matlab/tdu.m @@ -10,7 +10,7 @@ g_0 = 9.81; % [m/s^2] standard gravity unitless='[-]'; global debug; -debug = true; +debug = false; % === THRUSTER PARAMETERS === % IMPORT FROM FILE fid = fopen(filein,'r'); @@ -143,7 +143,7 @@ figure numplots=4;plotcounter=1; subplot(numplots,1,plotcounter) -plot(xcoord,radius);ylabel('radius');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); +plot(xcoord,radius);ylabel('Radius (m)');plotcounter=plotcounter+1;axis([0 xcoord(end) 0 inf]); % subplot(numplots,1,plotcounter) % semilogy(xcoord,A./A_t,xcoord(mark),A(mark)./A_t,'x');ylabel('A/A_t');plotcounter=plotcounter+1; subplot(numplots,1,plotcounter) @@ -335,10 +335,8 @@ %Solution as a min of residual [RES,i]=min((ASH1-ASH2).^2); % RES -if debug - shockArea=ASH1(i) - Mminus=MM(i) -end +shockArea=ASH1(i); +Mminus=MM(i); end function display(result) global debug