Home » 2013 (Page 2)

# Yearly Archives: 2013

## Composite Rocket Propellant

In this example we are going to build a chemistry table for an AP-composite rocket propellant, and investigate the effects due to increasing metallization of the grain. CEAGUI application from the NASA Glenn Research center was used to calculate molecular weight, ratio of specific heats and combustion temperature for different mixture ratio of AP/HTPB/AL. CEAgui Download for Microsoft Windows available http://www.grc.nasa.gov/WWW/CEAWeb/ceaguiDownload-win.htm here. Those results used to calculate characteristic velocity C* (figure of thermo-chemical merit for a particular propellant) and Ideal Isp.

clc; clear all; close all;

## Input

Combustor pressure – 3000 kPa,

Reactant temperature – 325 K

Ru = 8314;      % Universal Gas Constant
g0 = 9.81;      % [m/s^2] Gravitation Acceleration
Mr = p05(:,1);  % Mixure Ratio
Temp = [p05(:,2),p5(:,2),p10(:,2),p15(:,2),p25(:,2)];
MWv  = [p05(:,4),p5(:,4),p10(:,4),p15(:,4),p25(:,4)];

## Combustion Temperature

figure(1);
hold on
cstring='rgbcmk';              % color string
for i = 1:5
plot(Mr,Temp(:,i),[cstring(mod(i,7)+1),'-*']);
end

title('Combustion Temperature vs Mixture Ratio & Metallized Fuel Fraction ')
legend('Metallized Fuel Fraction Al/HTPB = 0.5%','5%','10%','15%','25%',...
'Location','Southeast');
xlabel('Mixture Ratio  M_O_x_i_d_i_z_e_r/M_F_u_e_l');
ylabel('Combustion Temperature T[K]')
hold off;

## Ratio of Specific Heats,Combustor

Gamma = [p05(:,3),p5(:,3),p10(:,3),p15(:,3),p25(:,3)];
figure(2);
hold on
for i = 1:5
plot(Mr,Gamma(:,i),[cstring(mod(i,7)+1),'-*']);
end

title('Ratio of Specific Heats,Combustor')
legend('Metallized Fuel Fraction Al/HTPB = 0.5%','5%','10%','15%','25%',...
'Location','Northeast');
xlabel('Mixture Ratio  M_O_x_i_d_i_z_e_r/M_F_u_e_l');
ylabel('Gamma')
hold off;

## Characteristic Velocity

Cdv = (Gamma*Ru./MWv.*Temp).^0.5./(Gamma.*((2./(Gamma +1)).^...
((Gamma + 1)./(Gamma -1))).^0.5);

figure(3);
hold on
for i = 1:5
plot(Mr,Cdv(:,i),[cstring(mod(i,7)+1),'-*']);
end
title('Characteristic Velocity')
legend('Metallized Fuel Fraction Al/HTPB = 0.5%','5%','10%','15%','25%',...
'Location','Southeast');
xlabel('Mixture Ratio  M_O_x_i_d_i_z_e_r/M_F_u_e_l');
ylabel('C^*')
hold off;

## Ideal Isp, based on infinite nozzle

Ispv  = Cdv./g0.*Gamma.*(2./(Gamma - 1).*(2./(Gamma +1)).^...
((Gamma + 1)./(Gamma -1))).^0.5;
figure(4);
hold on
for i = 1:5
plot(Mr,Ispv(:,i),[cstring(mod(i,7)+1),'-*']);
end

title('Ideal Isp, based on infinite nozzle')
legend('Metallized Fuel Fraction Al/HTPB = 0.5%','5%','10%','15%','25%',...
'Location','Southeast');
xlabel('Mixture Ratio  M_O_x_i_d_i_z_e_r/M_F_u_e_l');
ylabel('Specific Impulse, Isp[s]')
hold off;

## Orbital Inclination Change

Transfer from a LEO 350 km circular orbit with 53.4 deg inclination to a Geostationary Equatorial Orbit(GEO)

clear all; clc;
close all;

## Input

R_LEO = 6378 + 350;      % km
R_GEO = 42164;           % km
mu    = 398600;          % km^3/s^2
incl  = 53.4;            % deg

## Method 1: Hohman transfer from LEO to GEO and after inclination change

Rp  = R_LEO;
Ra  = R_GEO;
e   = (Ra - Rp)/(Ra + Rp);  % transfer orbit eccentricity
a   = (Ra + Rp)/2;          % transfer orbit semimajor axis
V_LEO   = (mu/R_LEO)^0.5;
Vp      = (2*mu/R_LEO - mu/a)^0.5;
Va      = (2*mu/R_GEO - mu/a)^0.5;
V_GEO   = (mu/R_GEO)^0.5;
dV_LEO  = abs(Vp - V_LEO);
dV_GEO  = abs(V_GEO - Va);
dV_Hoff = dV_GEO + dV_LEO;

% Inclination change at GEO
dV_incl = 2*V_GEO*sind(incl/2);
dV_total1 = dV_incl + dV_Hoff;

fprintf('Method 1: Hohman transfer from LEO to GEO and after inclination change\n');
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_total = %6.2f [km/s] \n',dV_total1);
Method 1: Hohman transfer from LEO to GEO and after inclination change
dV_Hoff  =   3.87 [km/s]
dV_incl  =   2.76 [km/s]
dV_total =   6.64 [km/s]

## Method 2: Inclination change in LEO and after Hohman transfer to GEO

dV_incl   = 2*V_LEO*sind(incl/2);
dV_total2 = dV_incl + dV_Hoff;

fprintf('\nMethod 2: Inclination change in LEO and after Hohman transfer to GEO\n');
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_total = %6.2f [km/s] \n\n',dV_total2);
fprintf('Method 2/Method 1 = %6.2f \n',dV_total2/dV_total1);
Method 2: Inclination change in LEO and after Hohman transfer to GEO
dV_incl  =   6.92 [km/s]
dV_Hoff  =   3.87 [km/s]
dV_total =  10.79 [km/s]

Method 2/Method 1 =   1.63

## Phased Array Antenna, Radiation Pattern and Array Configuration

Phased array consisting of N lined up individual isotropic antennas with the composite main antenna beam into vertical direction. This code based on the governing equation (J. Röttger, The Instrumental Principles of MST Radars, Ch.2.1) and shows how the radiation pattern changes as the parameters are modified.

clc; clear all;
close all;

Ratio between wavelength and distance between individual elements Linear Array of N isotropic radiators

N    = 64;
alfa = -90:0.1:90;       % Zenith angle
d_l  = [0.1,0.5,1,2,5] ; % Ratio
for i = 1:5
Ed = abs(sin(N*pi*d_l(i)*sind(alfa))./sin(pi*d_l(i)*sind(alfa)))/N;
fig = figure(i);
set(fig,'Position', [100, 100, 1049, 400]);
subplot(2,2,[1 3]);
plot(alfa,Ed);
grid on;
axis([-90 90 0 1]);
xlabel('Zenith angle [deg]');
title(['Ratio between distance and wavelength d/\lambda = ',num2str(d_l(i))]);

subplot(2,2,[2 4]);
polar(alfa*pi/180,Ed);
hold on;
polar((alfa+180)*pi/180,Ed);
xlabel(['Polar plot for the radiation pattern d/\lambda = ',num2str(d_l(i))]);
hold off;
end

Distance between individual elements

f_EISCAT =  233E6;   % Centre frequency 233 Hz
c = 3*10^8;          % Light Speed m/s
l = c/f_EISCAT;
d = [0.1,0.5,1,2,3];  % m
for i = 1:5
% Normalizied gain
Ed = abs(sin(N*pi*d(i)/l*sind(alfa))./sin(pi*d(i)/l*sind(alfa)))/N;
fig = figure(i+5);
set(fig,'Position', [100, 100, 1049, 400]);
subplot(2,2,[1 3]);
plot(alfa,Ed);
grid on;
axis([-90 90 0 1]);
xlabel('Zenith angle [deg]');
title([' Distance between individual elements (space weighting) = ',...
num2str(d (i)),' m,  F = 233Mhz']);
subplot(2,2,[2 4]);
polar(alfa*pi/180,Ed);
hold on;
polar((alfa+180)*pi/180,Ed);
xlabel('Polar plot for the radiation pattern');
hold off;
end

Number of antenna elements The distance between individual elements equals to the half wavelength of the signal

d_l = 0.5;
N   = [4 16 64 100 169]; % Number of elements was varied from 4 to 256.
for i = 1:5
Ed = abs(sin(pi*d_l*N(i)*sind(alfa))./sin(pi*d_l*sind(alfa)))/N(i);
fig = figure(i+10);
set(fig,'Position', [100, 100, 1049, 400]);
subplot(2,2,[1 3]);
plot(alfa,Ed);
grid on;
axis([-90 90 0 1]);
xlabel('Zenith angle [deg]');
title(['Number of elements = ',num2str(N(i)),', d/\lambda = 0.5']);
subplot(2,2,[2 4]);
polar(alfa*pi/180,Ed);
hold on;
polar((alfa+180)*pi/180,Ed);
xlabel('Polar plot for the radiation pattern');
hold off;
end

Radiation pattern change when the beam is not pointed vertically

clc;clear all;
N = 8;
d_l = 0.5;
alfa0 = 30;
alfa = -180:0.1:180;       % Zenith angle
Ed = abs(sin(pi*d_l*N*(sind(alfa) - sind(alfa0) ))./...
sin(pi*d_l*(sind(alfa) - sind(alfa0))))/N;
Ed0 = abs(sin(pi*d_l*N*(sind(alfa)))./...
sin(pi*d_l*sind(alfa)))/N;
Edl = 20*log10(Ed);
Edl0= 20*log10(Ed0);
fig = figure(16);
set(fig,'Position', [100, 100, 1049, 400]);
subplot(2,2,[1 3]);
plot(alfa,Ed);
hold on;
plot(alfa,Ed0,'r');
grid on;
axis([-90 90 0 1]);
xlabel('Zenith angle [deg]');
title(['Number of elements = ',num2str(N),', d/\lambda = 0.5']);
legend('\beta = 30 deg','\beta = 0 deg','Location','NorthWest')
subplot(2,2,[2 4]);
polar(alfa*pi/180,Ed);
hold on;
polar(alfa*pi/180,Ed0,'r');
hold on;
xlabel('Polar plot for the radiation pattern');
hold off;
figure(17);
plot(alfa,Edl);
hold on;
plot(alfa,Edl0,'r');
grid on;
legend('\beta = 30 deg','\beta = 0 deg','Location','NorthWest')
axis([-90 90 -100 0]);
xlabel('Zenith angle [deg]');
title(['Number of elements = ',num2str(N),', d/\lambda = 0.5']);


Electrical weighting techniques

N     = 64;
d_l   = 0.1;
Ed = abs(sin(pi*d_l*N*(sind(alfa) ))./...
sin(pi*d_l*sind(alfa)))/N;
Edl = 20*log10(Ed);
% Try to lower the first side lobe relative to main lobe
% To be able to reduce the side lobe levels, we design the array so its radiate
% more power towards the center, and less at the edges.
w  = zeros(1,N);
ws = zeros(1,N);
n = 1:N;
E = 0;Es = 0;
for n = 1:N;
if n <= N/2
w(n) = (n-1)/30.5;
else
w(n) = w(65 - n);
end
ws(n) = sin((n-1)*pi/(N-1));
E  = E + (w(n)*exp(complex(0,((2*pi*(n-1)*d_l)*sind(alfa)))));
Es  = Es + (ws(n)*exp(complex(0,((2*pi*(n-1)*d_l)*sind(alfa)))));
end
figure(19);
hold on; grid on;
plot(w,'g');
plot(ws,'r');
axis([1 64 0 1]);
ylabel('Normalized tapering weight');
xlabel('Array elements index');
legend('polynomial','sin');
hold off;
E = abs(E);
E = 20*log10(E/max(E));
Es = abs(Es);
Es = 20*log10(Es/max(Es));
fig = figure(18);
set(fig,'Position', [100, 100, 600, 700]);
hold on;
subplot(3,1,1);
plot(alfa,Edl);
grid on;
axis([-90 90 -100 0]);
legend('Isotropic');
title(['Number of elements = ',num2str(N),', d/\lambda = ',num2str(d_l)]);
subplot(3,1,2);
plot(alfa,Es,'r');
axis([-90 90 -100 0]);
grid on;
legend('Weighted Sin');
subplot(3,1,3);
plot(alfa,E,'g');
axis([-90 90 -100 0]);
grid on;
xlabel('Zenith angle [deg]');
legend('Weighted Poly');
hold off;

## Kepler Equation Solver Without Transcendental Function Evaluations

clc;clear all;

tic;
j = 1;
for M = 0:0.005:pi
i = 1;
for e = 0:0.005:1
c3 = 5/2 + 560*e;
a = 15/c3*(1-e);
b = -M/c3;
y = (b^2/4+ a^3/27)^0.5;
x = (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
c15 = 3003/14336 + 16384*e;
c13 = 3465/13312 - 61440*e;
c11 = 945/2816 + 92160*e;
c9 = 175/384 -70400*e;
c7 = 75/112 + 28800*e;
c5 = 9/8-6048*e;
x2 = x^2;x3 = x2*x;x4 = x3*x;
x5 = x4*x;x6 = x5*x;x7 = x6*x;
x8 = x7*x;x9 = x8*x;
x10 =x9*x; x11 = x10*x; x12 = x11*x;
x13 = x12*x; x14 = x13*x; x15 = x14*x;

f = c15*x15 + c13*x13 + c11*x11 + c9*x9 + c7*x7 + c5*x5 + c3*x3 +...
15*(1 - e)*x -M;
f1 = 15* c15 *x14 + 13* c13 *x12 + 11* c11 *x10 + 9* c9 *x8 + 7* c7 *x6 + ...
5 *c5 *x4 + 3 *c3 *x2 + 15*(1 - e);
f2 = 210* c15 *x13 + 156* c13 *x11 + 110* c11 *x9 + 72* c9 *x7 + ...
42* c7 *x5 + 20* c5 *x3 + 6* c3 *x;
f3 = 2730*c15 *x12 + 1716* c13 *x10 + 990* c11 *x8 + 504* c9 *x6 +...
210* c7 *x4 + 60* c5 *x2 + 6* c3;
f4 = 32760* c15 *x11 + 17160* c13 *x9 + 7920* c11 *x7 + 3024* c9 *x5 +...
840* c7 *x3 + 120* c5 *x;
f5 = 360360* c15 *x10 + 154440* c13 *x8 + 55440* c11 *x6 + 15120* c9 *x4 ...
+ 2520* c7 *x2 + 120* c5;
f6 = 3603600* c15 *x9 + 1235520* c13 *x7 + 332640* c11 *x5 + 60480* c9 *x3...
+ 5040* c7 *x;
f7 = 32432400* c15 *x8 + 8648640* c13 *x6 + 1663200* c11 *x4 + 181440* c9 *x2 ...
+ 5040* c7;
f8 = 25945920* c15 *x7 + 51891840* c13 *x5 + 6652800* c11 *x3 + ...
362880* c9 *x;
f9 = 1.8162144e9* c15 *x6 + 259459200* c13 *x4 + 19958400* c11 *x2 + ...
362880* c9;
f10 = 1.08972864e10* c15 *x5 + 1.0378368e9*c13 *x3 + 39916800* c11 *x;
f11 = 5.4486432e10* c15 *x4 + 3.1135104e9* c13 *x2 + 39916800* c11;
f12 = 2.17945728e11* c15 *x3 + 6.2270208e9* c13 *x;
f13 = 6.53837184e11* c15 *x2 + 6.2270208e9* c13;
f14 = 1.307674368e13* c15 *x;
f15 = 1.307674368e13* c15;

g1 = 1/2; g2 = 1/6; g3 = 1/24; g4 = 1/120;
g5 = 1/720; g6 = 1/5040; g7 = 1/40320; g8 = 1/362880;
g9 = 1/3628800; g10 = 1/39916800; g11 = 1/479001600;
g12 =1/6.2270208e9 ; g13 = 1/8.71782912e10; g14 = 1/1.307674368e12;

u1 = -f/f1;
h2 = f1 + g1* u1* f2;
u2 = -f/h2;
h3 = f1 + g1* u2* f2 + g2* u2^2*f3;
u3 = -f/h3;
h4 = f1 + g1* u3* f2 + g2* u3^2*f3 + g3* u3^3*f4;
u4 = -f/h4;
h5 = f1 + g1* u4* f2 + g2* u4^2*f3 + g3* u4^3*f4 + g4* u4^4*f5;
u5 = -f/h5;
h6 = f1 + g1* u5* f2 + g2* u5^2*f3 + g3* u5^3*f4 + g4* u5^4*f5 +...
g5* u5^5*f6;
u6 = -f/h6;
h7 = f1 + g1* u6* f2 + g2* u6^2*f3 + g3* u6^3*f4 + g4* u6^4*f5 +...
g5* u6^5*f6 + g6* u6^6*f7;
u7 = -f/h7;
h8 = f1 + g1* u7* f2 + g2* u7^2*f3 + g3* u7^3*f4 + g4* u7^4*f5 +...
g5* u7^5*f6 + g6* u7^6*f7 + g7* u7^7*f8;
u8 = -f/h8;
h9 = f1 + g1* u8* f2 + g2* u8^2*f3 + g3* u8^3*f4 + g4* u8^4*f5 +...
g5* u8^5*f6 + g6* u8^6*f7+ g7* u8^7*f8 + g8* u8^8*f9;
u9 = - f/h9;
h10 = f1 + g1* u9* f2 + g2* u9^2*f3 + g3* u9^3*f4 + g4* u9^4*f5 +...
g5* u9^5*f6 + g6* u9^6*f7+ g7* u9^7*f8 + g8* u9^8*f9 + ...
g9* u9^9*f10;
u10 = -f/h10;
h11 = f1 + g1* u10* f2 + g2* u10^2* f3 + g3* u10^3* f4 + g4* u10^4*f5...
+ g5* u10^5* f6 + g6* u10^6*f7+ g7* u10^7* f8 + g8* u10^8* f9...
+ g9* u10^9* f10 + g10* u10^10* f11;
u11 = -f/h11;
h12 = f1 + g1* u11* f2 + g2* u11^2* f3 + g3* u11^3* f4 + ...
g4* u11^4* f5 + g5* u11^5* f6 + g6* u11^6* f7...
+ g7* u11^7* f8 + g8* u11^8* f9 + g9* u11^9* f10 + g10* u11^10* f11...
+ g11* u11^11* f12;
u12 = -f/h12;
h13 = f1 + g1* u12* f2 + g2* u12^2* f3 + g3* u12^3* f4 +...
g4* u12^4* f5 + g5* u12^5* f6 + g6* u12^6* f7...
+ g7* u12^7* f8 + g8* u12^8* f9 + g9* u12^9* f10 +...
g10* u12^10* f11 + g11* u12^11* f12 + g12* u12^12* f13;
u13 = -f/h13;
h14 = f1 + g1* u13* f2 + g2* u13^2* f3 + g3* u13^3* f4 + ...
g4* u13^4* f5 + g5* u13^5* f6 + g6* u13^6* f7...
+ g7* u13^7* f8 + g8* u13^8* f9 + g9* u13^9* f10 + g10* u13^10*f11...
+ g11* u13^11* f12+ g12* u13^12* f13 + g13* u13^13* f14;
u14 = -f/h14;
h15 = f1 + g1* u14* f2 + g2* u14^2* f3 + g3* u14^3* f4 + ...
g4* u14^4* f5 + g5* u14^5* f6 + g6* u14^6* f7...
+ g7* u14^7* f8 + g8* u14^8* f9 + g9* u14^9* f10 +...
g10* u14^10* f11 + g11*u14^11* f12+ g12* u14^12* f13 + ...
g13* u14^13* f14 + g14* u14^14* f15;
u15 = -f/h15;

x =  x + u15;
w = x - 0.01171875*x^17/(1 + e);
E = M + e*(-16384*w^15 + 61440*w^13 - 92160*w^11+ 70400*w^9 -...
28800*w^7 + 6048*w^5 - 560*w^3 + 15*w);
E_nit(i,j) = E;
i = i + 1;
end
j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 2.591335 [s]


## Newton Iterative Solver

tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.005:pi
i = 1;
for e = 0:0.005:1
En  = M;
Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
while ( abs(Ens-En) > eps )
En = Ens;
Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
end;
E_it(i,j) = Ens;
i = i + 1;
end
j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 3.520123 [s]
Error = log10(abs(E_it - E_nit));
e = 0:0.005:1;
M = (0:0.005:pi);
contourf (M,e,Error);
ylabel('e - Eccentricity');
title('Error in E calculation on Log scale [rad] ');
colorbar;

REFERENCES
Adonis Pimienta-Pe˜nalver and John L. Crassidis†,ACCURATE KEPLER EQUATION
SOLVER WITHOUT TRANSCENDENTAL FUNCTION EVALUATIONS,
University at Buffalo, State University of New York, Amherst, NY, 14260-4400

## Kepler’s equation, Iterative and Non-Iterative Solver Comparison

Kepler’s equation M = E-e*sin(E) M – Mean anomaly [0..pi] rad e – Eccentricity [0..1] E – Eccentric anomaly [rad]

clear all;
clc;

## Non-Iterative Solver

tic;
j = 1;
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
a   =  (1 - e)*3/(4*e+0.5);
b   =  -M/(4*e+0.5);
y   =  (b^2/4 + a^3/27)^0.5;
x   =  (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
w   =  x - 0.078*x^5/(1 + e);
E   =  M + e*(3*w - 4*w^3);
% Two times applying 5th-order Newton correction
for k = 1:2
f   =  (E - e*sin(E)- M);
fd  =  1 - e*cos(E);
f2d =  e*sin(E);
f3d = -e*cos(E);
f4d =  e*sin(E);
E = E -f/fd*(1 + f*f2d/(2*fd^2) + f^2*(3*f2d^2 - fd*f3d)/(6*fd^4)+...
(10*fd*f2d*f3d - 15*f2d^3-fd^2*f4d)*f^3/(24*fd^6));
end
E_nit(i,j) = E;
i = i + 1;
end
j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 0.149860 [s]

## Newton Iterative Solver

tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.01:pi
i = 1;
for e = 0:0.01:1
En  = M;
Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
while ( abs(Ens-En) > eps )
En = Ens;
Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
end;
E_it(i,j) = Ens;
i = i + 1;
end
j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 0.893673 [s]
Error = log10(abs(E_it - E_nit));
e = 0:0.01:1;
M = (0:0.01:pi);
contourf (M,e,Error);
Mechanics and Dynamical Astronomy,Vol. 40, 1987, pp. 329–334.