Small Satellites

Home » 2013 » May (Page 2)

Monthly Archives: May 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
load matlab2
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;project_01

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

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

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

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]');
    ylabel('Normalized radiation pattern');
    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

PhasedArray_01 PhasedArray_02 PhasedArray_03 PhasedArray_04 PhasedArray_05

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;
    % Ed - Array Factor, Radiation patern for isotropic radiators
    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]');
    ylabel('Normalized radiation pattern');
    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

PhasedArray_06 PhasedArray_07 PhasedArray_08 PhasedArray_09 PhasedArray_10

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]');
    ylabel('Normalized array factor,Radiation pattern');
    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

PhasedArray_11 PhasedArray_12 PhasedArray_13 PhasedArray_14 PhasedArray_15

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]');
    ylabel('Normalized array factor,Radiation pattern');
    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]');
    ylabel('Radiation pattern dB');
    title(['Number of elements = ',num2str(N),', d/\lambda = 0.5']);
PhasedArray_16

PhasedArray_17

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');
    ylabel('Radiation pattern dB');
    title(['Number of elements = ',num2str(N),', d/\lambda = ',num2str(d_l)]);
    subplot(3,1,2);
    plot(alfa,Es,'r');
    axis([-90 90 -100 0]);
    ylabel('Radiation pattern dB');
    grid on;
    legend('Weighted Sin');
    subplot(3,1,3);
    plot(alfa,E,'g');
    axis([-90 90 -100 0]);
    grid on;
    xlabel('Zenith angle [deg]');
    ylabel('Radiation pattern dB');
    legend('Weighted Poly');
    hold off;

PhasedArray_18 PhasedArray_19

 

%d bloggers like this: