Trajectory of a 3d body

12 次查看（过去 30 天）
Rakshit Patil
Rakshit Patil 2024-9-20
编辑： Ashok 2024-10-11
Problem Statement: -
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.
the code: -
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)
% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)
% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)
% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)
% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)
% Distance between center of mass and center of pressure
d_cp_cm = cp_z - cm_z;
% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)
% Calculate initial velocity components
v0x = v0_aircraft - v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);
% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];
% Time settings
tspan = [0 300]; % Time span for integration
% ODE solver options
options = odeset('Events', @ground_impact, 'RelTol', 1e-6, 'AbsTol', 1e-6);
% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);
% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);
% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);
% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel('Horizontal distance (km)');
ylabel('Altitude (km)');
title('2D Trajectory of Body (Advanced 6DOF Simulation)');
grid on;
% Plot velocity
figure;
plot(t, v);
xlabel('Time (s)');
ylabel('Velocity (m/s)');
title('Velocity vs Time');
grid on;
% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel('Angular Position (degrees)');
title('Angular Motion');
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel('Time (s)');
ylabel('Angular Velocity (deg/s)');
grid on;
% Display results
disp(['Time of flight: ', num2str(t(end)), ' s']);
disp(['Horizontal distance traveled: ', num2str(x(end)/1000), ' km']);
disp(['Impact velocity: ', num2str(v(end)), ' m/s']);
disp(['Final angle: ', num2str(rad2deg(theta(end))), ' degrees']);
disp(['Max altitude: ', num2str(max(h)/1000), ' km']);
disp(['Max velocity: ', num2str(max(v)), ' m/s']);
% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp(['State at t = ', num2str(t(idx)), ' s:']);
disp([' Position: (', num2str(x(idx)/1000), ' km, ', num2str(h(idx)/1000), ' km)']);
disp([' Velocity: ', num2str(v(idx)), ' m/s']);
disp([' Angle: ', num2str(rad2deg(theta(idx))), ' degrees']);
end
end
end
function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);
% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) - theta;
% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);
% Calculate Mach number
M = v / a;
% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);
% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment
% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;
% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end
% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) - Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration
% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end
function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air
% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)
if h < 11000 % Troposphere
T = T0 - 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h - 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h - 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end
rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end
function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models
% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);
% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;
% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 - M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility
% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;
% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end
function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
end

请先登录，再回答此问题。

采纳的回答

Ashok
Ashok 2024-10-7
编辑：Ashok 2024-10-11
Hi Rakshit,
In the code shared, I am assuming that 'theta' is the pitch angle of the projectile,'alpha' is the angle of attack, the positive angles indicate pitch up and negative angles indicate pitch down. Based on this, you should check the correctness of the following relationships:
1. The initial velocity components in the ‘aircraft_launch_trajectory_6DOF_advanced’ function.
v0x = v0_aircraft - v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);
2. Relations in the ‘equations_of_motion_6DOF_advanced’ function for 'alpha', total moment, and acceleration.
alpha = atan2(vy, vx) - theta;
M_total = M_aero + Fl * d_cp_cm;
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) - Fl * cos(theta) + boost_fy) / m;
Additionally, the drag coefficient calculated in the ‘aero_coefficients’ function and the characteristic length in the moment relation, ‘M_aero = q * Cm * A * length;’ seem quite high.
Since the ‘aero_coefficients’ function doesn’t account for the stall region, it might be better to use a smaller launch angle. Instead of ‘angle = -45;’ consider using a launch angle in the range -15 to 15 degrees.
To improve the visualization of the projectile’s trajectory, you can make the X and Y axes have the same scale by using the command ‘axis equal’.
You can find more details about this command in the link below:
https://www.mathworks.com/help/matlab/ref/axis.html#buk989s-1-style
I believe this will assist you!

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Aerospace Blockset 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by