六自由度外弹道计算程序一步到位版（含气动参数计算）
发布时间
阅读量：
阅读量
声明：
本程序仅用于PYTHON编程学习交流用途。
涉及的计算公式主要来自《弹箭外弹道学 第二版》一书；而探讨六自由度外弹道计算的知识，则需要进一步了解相关的理论基础与应用方法，在具体资料中可找到详细的说明。
该系统采用的弹道模型具有6自由度刚体运动学描述，并在《弹箭外弹道学 第二版》一书中有所体现
1. 程序介绍：
该程序由三个data.csv文件组成，负责接收和转导计算数据；此外又有一个名为SUPPORT的文件夹用来存储无界面版的DATCOM软件；其具体内容如图所示。
其中：
001.csv用于存放六月自由度外弹道相关计算基础输入参数，如下图所示：
002.csv用于存放 DATCOM计算出的气动参数结果，如下图所示：
003.csv用于存放 DATCOM气动计算的输入参数，如下图所示：
程序计算过程：
用户基于外弹道计算的需求对包含相关参数的两个CSV文件（编号分别为001和003）进行了设置，并随后在程序中启动了六自由度的外弹道计算流程。
计算结束后，程序自动绘制出基本的弹丸飞行轨迹曲线，如下图所示：
在程序目录中同时生成并存储一个名为result.csv的文件用于存储计算过程中的所有相关变量以便于后续的数据分析和可视化处理该文件的内容如图所示
2. 程序代码实现：
2.1 基础内容设置
首先引入程序相关的库：
# 引入必要的库
#..........................................................................................
import math #数学计算库
from matplotlib import pyplot as plt #绘图库
import numpy as np #数值计算库
from scipy.integrate import solve_ivp #ODE求解库
from scipy.interpolate import interp1d, RectBivariateSpline #插值库
import pandas as pd #数据处理库
import chardet #编码检测库
import subprocess #子进程库
import os #操作系统库
import re #正则表达式库
#..........................................................................................
python
并定义用于检测文件编码的函数：
# 定义函数，用于检测文件编码
def detect_encoding(file_path):
try:
with open(file_path, 'rb') as f:
rawdata = f.read()
result = chardet.detect(rawdata)
return result['encoding']
except Exception as e:
print(f"Error detecting encoding: {e}")
return None
python
2.2 相关参数读取及气动参数计算
请获取003.csv文件内容，并提取其中用于气动计算的相关参数信息；随后将这些数据以指定的方式写入至DATCOM的运算输入文件中。
# 指定弹丸气动参数计算DATCOM软件的输入文件路径
csv_file_path_datcom_params = '003.csv'
# 尝试检测弹丸气动参数计算DATCOM软件的输入文件编码
detected_encoding = detect_encoding(csv_file_path_datcom_params)
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
exit()
else:
try:
# 读取CSV文件
csv_reader = pd.read_csv(csv_file_path_datcom_params, encoding=detected_encoding, header=0)
# 创建空列表来存储数据
first_column_data = []
second_column_data = []
third_column_data = []
# 遍历CSV文件的每一行
for index, row in csv_reader.iterrows():
first_column_data.append(row[0])
second_column_data.append(row[1])
third_column_data.append(row[2])
data = list(zip(first_column_data, second_column_data, third_column_data))
# 打印读取结果
#print(data)
# 赋值给变量
XCG = second_column_data[0]
LNOSE = second_column_data[1]
DNOSE = second_column_data[2]
BNOSE = second_column_data[3]
LCENTR = second_column_data[4]
DCENTR = second_column_data[5]
LAFT = second_column_data[6]
DAFT = second_column_data[7]
# 打印弹丸气动参数
print(f"XCG = {XCG}")
print(f"LNOSE = {LNOSE}")
print(f"DNOSE = {DNOSE}")
print(f"BNOSE = {BNOSE}")
print(f"LCENTR = {LCENTR}")
print(f"DCENTR = {DCENTR}")
print(f"LAFT = {LAFT}")
print(f"DAFT = {DAFT}")
except Exception as e:
print(f"Error processing CSV file 003.csv: {e}")
exit()
python
DATCOM程序并读取相关计算结果转存至002.csv中。
读取外弹道计算基础参数：
# 指定六自由度外弹道计算基础参数文件路径
basic_params_csv_file_path = '001.csv'
detected_encoding = detect_encoding(basic_params_csv_file_path)
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
else:
try:
# 读取CSV文件
csv_reader = pd.read_csv(basic_params_csv_file_path, encoding=detected_encoding, header=0)
# 创建空列表来存储数据
first_column_data = []
second_column_data = []
third_column_data = []
# 遍历CSV文件的每一行
for index, row in csv_reader.iterrows():
first_column_data.append(row[0])
second_column_data.append(row[1])
third_column_data.append(row[2])
data = list(zip(first_column_data, second_column_data, third_column_data))
# 打印读取结果
print(data)
except Exception as e:
print(f"Error processing CSV file 001.csv: {e}")
python
读取外弹道计算气动参数：
# 指定六自由度外弹道气动基础参数文件路径
aero_params_csv_file_path = '002.csv'
detected_encoding = detect_encoding(aero_params_csv_file_path)
# 读取CSV文件
if detected_encoding is None:
print("Failed to detect encoding. Exiting.")
else:
try:
# 读取CSV文件
df = pd.read_csv(aero_params_csv_file_path, encoding=detected_encoding, header=0)
# 创建一个字典来存储二维数组
arrays = {}
# 遍历DataFrame，提取每个二维数组
current_array_name = None
current_array = []
for index, row in df.iterrows():
array_name = row[0]
if array_name != current_array_name: # 判断二维数组的名字是否改变
if current_array_name is not None: # 如果当前二维数组的名字不为空，则将上一个二维数组添加到字典中
arrays[current_array_name] = current_array # 添加到字典中
current_array_name = array_name # 更新当前二维数组的名字
current_array = [] # 重置当前二维数组
current_array.append(row[1:].tolist()) # 添加到当前二维数组中
# 添加最后一个数组
if current_array_name is not None:
arrays[current_array_name] = current_array
# 打印结果
for array_name, array in arrays.items():
print(f"{array_name}:")
for row in array:
print(row)
print()
except Exception as e:
print(f"Error processing CSV file 002.csv: {e}")
python
2.3 六自由度外弹道计算
调用ode求解器，使用RK45方法求解微分方程组。
# 调用ode求解器，使用RK45方法求解微分方程组
def wrap_sdofmotion(t, mp, global_params):
# 解包全局参数
X_cm, tau, a, C, A, m, d, l, S, v_0, theta_a_0, psi_2_0, phi_a_0, phi_2_0, gamma_0, omiga_xi_0, omiga_eta_0, omiga_zeta_0, x_0, y_0, z_0, w, a_W, g, rou, beta_D_1, beta_D_2, Lambda, delta_f, delta_N, delta_M, gamma_0_1, gamma_0_2, Omiga_E, a_N, ma1, Delta, CA, CN, X_CP, CNA, CYPA, Cnp1, Cnp3, Cnp5= global_params
return sdofmotion(t, mp, global_params)
# 设置初始条件和参数
t0 = 0.0
mp0 = np.array([v_0,theta_a_0,psi_2_0,omiga_xi_0,omiga_eta_0,omiga_zeta_0,phi_a_0,phi_2_0,gamma_0,x_0,y_0,z_0]) # 设置初始状态向量
t_end = t_estimation # 设置积分终止时间
global_params = [X_cm, tau, a, C, A, m, d, l, S, v_0, theta_a_0, psi_2_0, phi_a_0, phi_2_0, gamma_0,
omiga_xi_0, omiga_eta_0, omiga_zeta_0, x_0, y_0, z_0, w, a_W, g, rou, beta_D_1, beta_D_2,
Lambda, delta_f, delta_N, delta_M, gamma_0_1, gamma_0_2, Omiga_E, a_N,
ma1, Delta, CA, CN, X_CP, CNA, CYPA, Cnp1, Cnp3, Cnp5] # 设置全局参数列表
# 使用solve_ivp进行积分，指定方法为RK45并设置容差
sol = solve_ivp(wrap_sdofmotion, (t0, t_end), mp0, args=(global_params,), method='RK45', rtol=1e-7, atol=1e-7)
python
2.4 结果处理
将结果输出值result.csv文件，并绘制出基本外弹道飞行轨迹曲线。
# 处理结果
ts = sol.t # 时间序列
mps = sol.y.T # 状态向量序列
#..............................................................................................................................................
# 结果处理及导出
#..............................................................................................................................................
# 导出计算结果至外部csv文件中
df = pd.DataFrame(mps, columns=['v', 'theta_a', 'psi_2', 'omiga_xi', 'omiga_eta', 'omiga_zeta', 'phi_a', 'phi_2', 'gamma', 'x', 'y', 'z'])
df['t'] = ts
df.to_csv('result.csv', index=False)
print('The calculation result has been exported to the external csv file successfully!')
# 绘制结果图
fig, ax = plt.subplots(2, 2, figsize=(12, 10))
ax[0, 0].plot(ts, mps[:, 0], label='v')
ax[0, 0].set_xlabel('t')
ax[0, 0].set_ylabel('v')
ax[0, 0].legend()
ax[0, 0].set_title('Velocity')
ax[0, 1].plot(ts, mps[:, 9], label='X')
ax[0, 1].set_xlabel('t')
ax[0, 1].set_ylabel('X')
ax[0, 1].legend()
ax[0, 1].set_title('Position X')
ax[1, 0].plot(ts, mps[:, 10], label='Y')
ax[1, 0].set_xlabel('t')
ax[1, 0].set_ylabel('Y')
ax[1, 0].legend()
ax[1, 0].set_title('Position Y')
ax[1, 1].plot(ts, mps[:, 11], label='Z')
ax[1, 1].set_xlabel('t')
ax[1, 1].set_ylabel('Z')
ax[1, 1].legend()
ax[1, 1].set_title('Position Z')
fig.suptitle('Time Series Analysis of Motion Data')
plt.tight_layout() # 调整布局以防止标题重叠
plt.show()
python
3. 程序打包
利用PyInstaller工具将程序内容打包生成可以直接运行的.exe文件；配置好相关输入参数后，用户可以通过运行该.exe文件来进行外弹道解算。
利用PyInstaller工具将程序内容打包生成可以直接运行的.exe文件；配置好相关输入参数后, 用户可以通过运行该.exe文件来进行外弹道解算.
4. 总结
基于PYTHON技术的基础上开发出了一种简便的方法实现了将第三-party aerodynamic parameter calculation software（此处指DATCOM无界面版）集成到本程序中。这种方法大幅降低了整个计算流程所需的资源投入（相比之前需从006.dat文件中逐一查找数据并进行复制粘贴的操作）。此外还利用pyinstller工具实现了程序的打包成独立运行的.exe文件这一功能使用户无需额外准备相关环境即可轻松运行该程序在不同设备上使用更加便捷。关于外弹道相关的知识以及与DATCOM气动参数计算软件相关的详细手册和技术资料还有编程过程中积累的经验与心得我们热忱欢迎对此感兴趣的所有人进行交流讨论具体联系方式为QQ：3383016937
全部评论 (0)
还没有任何评论哟~