杂谈21深圳杯A题

目录

一、前言

二、模型建立（设想）

三、源码解读

从源码中获得的知识点

---

一、前言

这次和我的小伙伴们合作了深圳杯的数模比赛，我们选择的是A题，题目是主要是建立火星登陆器着陆的问题，我在团队中是负责编程实现的（虽然，我是matlab小白）。在建立整个模型之前，我参考了别人的一个源码https://ww2.mathworks.cn/matlabcentral/fileexchange/58135-mars-cubesat-descent-and-impact-simulation

这是链接，有需要的可以下载（如果涉及侵权、请告知作者删除）

以下主要是对题目模型的设想与对该源码的部分函数和代码块的解读

二、模型建立（设想）

该设想只是我从一些论文获得之后针对该题目的模型假设

分析题目很容易知道，题目是要求我们建立火星着陆器从火星大气层到陆地的过程模型，并不需要考虑卫星绕火星运动等问题。查阅相关资料、尤其是我国天问卫星的资料之后对这一过程模型可以依时间划分 从进入大气层加速向下到离地一定高度时展开降落伞、再到降落伞脱离采用动力减速至悬停最后着陆。这些过程涉及到的问题方向很多，比如一开始的进入大气层时的入射角度、本身的速度……

而我所参考的源码、其实只涵盖了进入大气层加速以及降落伞打开的着陆是否安全这几个方面展开、且大气加速时属于垂直向下的。

三、源码解读

以下是加入了我个人的注释的matlab源码

% Mars Descent Simulaltion
% 火星下降模拟

% This code simulates, animates, analizes impact force and saves a video file
% of a custom descent of a lander on Mars. In the animation, blue represents
% parachute deployment, green represents impact survival and red represents
% impact accident.
% 改代码模拟，设置动画，分析冲击力，并保存火星着陆器自定义设定的视频文件。在动画中，蓝色表示降落伞展开，
% 绿色表示碰撞幸存以及红色表示碰撞事故。

% Rodrigo Santos Valente da Costa
% 罗德里戈·桑托斯·瓦伦特·达·科斯塔
% RodrigoS.V.daCosta@gmail.com
% 邮箱下述是时间
% June 2016

% clear 从工作区删除项目，释放系统内存, (就是清楚工作空间的所有变量、函数和MEX文件)
% close all 关闭所有的figuer（动画、图象）窗口
clear; close all;

% Lander and Planet Settings
% 着陆器和行星设置

% 着陆器在动画中的尺寸，单位为m
L1= 1000; % Size of the lander in the animation (does not affect calculations), m

% 重力加速度，单位为 m/s^2
g= 3.711; % Acceleratiojn of gravity, m/s^2

% 立方体卫星的投影面积，单位为 m^2
A= 1; % Projected area of the cubesat, m^2

% 阻力系数（每立方体 1.05），无单位
Cd= 1.05; % Coefficient of drag (1.05 for a cube), unitless

% 着陆器质量，千克
m= 5; % Mass of the lander, kg

% 大气初始密度， kg/m^3
rho= 0.0044; % Initial density of the atmosphere, kg/m^3

% 初始阻力，牛顿
drag= 0; % Initial drag force, N

%% 分隔线————

% 火星质量，千克
mars_mass= 6.39*10^23; % Mass of mars, kg

% 引力常数
G= 6.67*10^(-11); % Gravitational constant

% 火星半径，米
mars_radius= 3.39 * 10^6; % Mars radius, m

% 降落伞展开高度，米
para_altitude= 12000; % Altitude of parachute deployment, m

%% 分隔线————

% 初始条件
% Initial Conditions

% 初始高度，米
y0= 15000; % Initial altitude, m

% 初始垂直速度，m/s
v0= -975; % Initial vertical velocity, m/s

% 初始加速度，m/s^2
yAcc= 0; % Initial acceleration, m/s^2

% ？垂直速度赋值？ m/s
yVel= v0; % Vertical velocity, m/s

% ？初始高度赋值？ 米
y= y0; % Altitude, m

%% 分隔线————

% 声明了一系列的数组变量，（应该是用来记录相关数据）？

% 高度数据
ydata= [];
% 垂直速度数据
yVeldata= [];
% （垂直）加速度
yAccdata= [];
% 时间
time= [];
% 密度
density= [];
% 重力
gravity= [];
% 阻力数据
dragdata= [];

%% 分隔线————

% 动画的外部框架
% 初始高度，和着陆器在动画中的尺寸 之和
frame_size= y0 + L1; % External frame of the animation

%% 分隔线————

% 时间设定
% Time Settings
% 模拟持续的时间，秒
duration= 690; % Duration of simluation, s
% 模拟的增量时间（延迟时间），秒
delta_t= 0.01; % Delta time of the simulation, s
% 时间步长的数字（时频）
% round() 函数，取整（四舍五入）时间步数
time_steps= round(duration/delta_t); % Number of time steps
% 帧率，2帧每秒
frameRate= 2; % Animation frame rate, frames/s
% 帧时，一帧的时间
framePeriod= 1./frameRate; % Period of one frame, s

%% 分隔线————
% 绘制初始图形
% Draw initial figure

% figure() 函数，建立图像，其参数 1 表示第一幅图象
fig= figure(1);
% axes() 函数，创建坐标系图形图象，其中第一个参数表示图象文件名字，关键字fig表示文件格式为fig（figure）格式
axs=axes('Parent',fig);

%% 分隔线————
% 打开移动文件
% Open movie file

% VideoWriter()函数，建立视频类文件，其参数为文件名
vidObj= VideoWriter('mars_descent_simulation.avi');
% 设置动画视频文件的帧率，1秒2帧
vidObj.FrameRate= frameRate;
% 打开动画视频文件
open(vidObj);

%% 分隔线————
% 绘制登陆器
% Draw Lander	
% rectangle()函数，创建二维矩阵对象，第一个参数‘Position’表示位置，其后的[x,y,w,h]指从(x,y)开始绘制宽w高h的矩形
% ‘Curvature’表示曲率，默认为[0,0]
cubesat= rectangle('Position', [0, -(frame_size - 1)*L1, L1, L1], 'Curvature', [0,0]);
% line()函数，画线函数，使用其画边界border
border= line([-frame_size/2, -frame_size/2, frame_size/2, frame_size/2, -frame_size/2],...
    [0, frame_size, frame_size, 0, 0]);
%axis()函数，主要对坐标轴进行缩放的操作，[xmin, xmax, ymin, ymax], 设置坐标轴x轴和y轴的限制范围
axis(axs, [-frame_size/2, frame_size/2, 0, frame_size]);
% ‘square’表示将坐标图设置为方形。横轴和纵轴的比例为1:1
axis(axs, 'square');
% 打开轴值、坐标轴标签、刻度、背景
axis(axs, 'on'); % Turns axes values on
% 关闭x轴
set(axs,'XTick',[]); % Turns x axis off
% label指标签，ylabel()函数，就是对y轴标签 ‘Altitude(m)’表示标签为‘海拔高度（米）’
ylabel('Altitude (m)');
% 声明压力改变变量，用了画线函数，其中参数‘Color’和‘red’表示用红线表示
pressure_change= line([-frame_size, frame_size], [7000,7000], 'Color', 'red'); ...
	% 在大气密度变化的位置绘制一条线
    % Draws a line where the atmospheric density changes

%% 分隔线————
% 视频生成
% Video generation
% 一个for语句，在时间步数中循环
for i= 1:time_steps
    t= (i-1)*delta_t;
    % 根据高度计算重力加速度
    % Calculation of the gravity acceleration depending on the altitude
	% 用了万有引力方程
    g= G*mars_mass/(y + mars_radius)^2;
    
	% 当离地高度小于等于0时，停止（就是触地停止）
    if y<=0 % Stops when it hits the ground
        break;
    end	%这里的end对应的是if语句
    
	% 当离地高度大于降落伞展开高度时
    if y>para_altitude
		% 设置卫星的位置和颜色
        set(cubesat, 'Position', [-L1/2, y, L1, L1], 'FaceColor', 'black');
		% 立方体卫星的阻力系数，无单位
        Cd= 1.05; % Coefficient of drag of Cubesat, unitless
		% 仅立方体的表面面积
        A= 0.01; % Surface area of the Cubesat only
    end 
	%当离地高度小于等于降落伞展开高度时
    if y<=para_altitude
		% 设置卫星的位置和颜色
        set(cubesat, 'Position', [-L1/2, y, L1, L1], 'FaceColor', 'blue');
		% 此时降落伞打开，阻力系数改变（其实就是降落伞的阻力系数）
        Cd= 1.3; % Coefficient of drag of parachute, unitless
		% 降落伞的表面积，平方米
        A= 2; % Surface area of the parachute, m^2
    end     
    
	% 根据海拔高度计算大气密度
    % Atmospheric density calculation depending on the altitude
    if y>7000
		% temp表示一个暂态变量，具体公式我也不知道
        temp=-23.4 - 0.00222*y;
		% exp()函数，表示对其内的参数，求其的e指数
        pressure= 0.699*exp(-0.00009*y);
    end
    if y<7000
        temp= -31 - 0.000998*y;
        pressure= 0.699*exp(-0.00009*y);
    end
	% 计算此时的大气密度
    rho= pressure/(.1921*(temp+273.1));
    
	% 重新设置视频的数据并且获取位置/速度数据
    % Resettig video data and getting position/velocity
	
	% 如果获取时间（带有延迟时间）与帧时取模 为0
    if (mod(t,framePeriod) == 0)
		% 当前节点帧赋值为获取帧
        currFrame= getframe;
		% 建立视频文件
        writeVideo(vidObj, currFrame);
		% 数组赋值操作，time 后跟t相当于把time数组原有的元素后面加上t（t也可以是数组）
        time= [time t];
        ydata= [ydata y];
        yVeldata= [yVeldata yVel];
        yAccdata= [yAccdata yAcc];
        density= [density rho];
        gravity= [gravity g];
    end
        
	% 拉力（升力、下降阻力）
    drag= A/2*Cd*rho*yVel^2; % Drag force, N
    
	% 加速度、方向取拉力方向
    yAcc= (drag - m*g)/m; % Acceleration considering drag, m/s^2
    
	% 对应的速度
    yVel= yVel + yAcc*delta_t; % Velocity considering drag, m/s
    
	% 对应的竖直高度
    y= y + yVel*delta_t; % Altitude considering drag, m
end	%循环终止

% 计算的影响因素
% Impact calculation

% 碰撞后减速距离？？？
s= 0.001; % Slow down distance after impact, m
% 冲击力？？
impactForce= 1/2*m*yVel^2/s; % Impact force, N
% 冲击压强？？
impactStress= impactForce/0.01; % Impact stress, Pa

% 用冲击压强来考虑是否能够安全着陆（撞击地面）
if impactStress<276*10^6
    set(cubesat, 'Position', [-L1/2, y, L1, L1], 'FaceColor', 'green');
	% disp() 类似c中的printf，作输出打印
    disp('Cubesat survived the land.');
end
if impactStress>=276*10^6
    set(cubesat, 'Position', [-L1/2, y, L1, L1], 'FaceColor', 'red');
    disp('Cubesat did not survive the land.');
end

% elapsedTime表示拖延时间（全程总时间）
elapsedTime= t;
% 输出下降和着陆耗时
disp('The descent and landing process took (in seconds):');
disp(elapsedTime);

% 关闭动画文件
%Close the movie file.
close(vidObj);

% 绘制时间与高度图象
% Plot Time vs. Altitude
figure(2);
% plot函数，绘制二维图像
plot(time, ydata);
xlabel('Time (s)');
ylabel('Altitude (m)');
% 横纵轴范围，颜色
line([0, t], [7000,7000], 'Color', 'red');

% 绘制时间和速度图象
% Plot Time vs. Velocity
figure(3);
plot(time, yVeldata);
xlabel('Time (s)');
ylabel('Velocity (m/s)');

% 绘制时间和加速度图象
% Plot Time vs. Acceleration
figure(4)
plot(time, yAccdata);
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');

% 绘制时间和大气密度图象
% Plot Time vs. Atmospheric Density
figure(5)
plot(time, density);
xlabel('Time (s)');
ylabel('Air density (kg/m^3)');

% 绘制时间和重力图象
% Plot Time vs. Gravity
figure(6)
plot(time, gravity);
xlabel('Time (s)');
ylabel('Gravity Acceleration (m/s^2)');

从源码中获得的知识点

前面提到我自己是一个matlab小白，那么从这份源码中其实学到了不少东西的

1、相关的动画和制图函数：VideoWriter()、axes()、figure()、（figure用来输出图片）、plot（）

2、设定函数: set（（x,y）...）、边界函数border()、画线函数line（）、标签函数labe();

3、数组类变量的设定和操作

        1）matlab中数组类的设定比较简单、只需 "变量名 = []" 即可、有点类似于 C++中的new的作用

        2）数组操作可以用 “数组名 = [数组名 变量]”这样子来填写数组元素、注意这个变量既可以是元素、也可以是一个数组。