CN116894302B - 一种基于精算迎风面积的leo航天器大气阻力算法 - Google Patents

一种基于精算迎风面积的leo航天器大气阻力算法 Download PDF

Info

Publication number
CN116894302B
CN116894302B CN202311161851.XA CN202311161851A CN116894302B CN 116894302 B CN116894302 B CN 116894302B CN 202311161851 A CN202311161851 A CN 202311161851A CN 116894302 B CN116894302 B CN 116894302B
Authority
CN
China
Prior art keywords
spacecraft
grid
area
windward
projection
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202311161851.XA
Other languages
English (en)
Other versions
CN116894302A (zh
Inventor
唐小伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Original Assignee
Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center filed Critical Ultra High Speed Aerodynamics Institute China Aerodynamics Research and Development Center
Priority to CN202311161851.XA priority Critical patent/CN116894302B/zh
Publication of CN116894302A publication Critical patent/CN116894302A/zh
Application granted granted Critical
Publication of CN116894302B publication Critical patent/CN116894302B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种基于精算迎风面积的LEO航天器大气阻力算法，涉及风洞试验领域，包括：S1、在通过航天器数模生成表面网格后，基于网格面元与迎风截面格子投射的原理对航天器迎风面积S W 进行精算；S2、将S1中得到的航天器迎风面积S W 引入至LEO航天器大气阻力计算中。本发明提供一种基于精算迎风面积的LEO航天器大气阻力算法，可对航天器飞行过程中任何姿态精确计算迎风面积，并应用于LEO航天器大气阻力中，计算精度更高。

Description

一种基于精算迎风面积的LEO航天器大气阻力算法
技术领域
本发明涉及风洞试验领域，具体涉及在风洞试验中集成空气动力学、计算机辅助设计及数模网格解析几何方面工程技术，共同协作解决相关试验问题的技术应用领域。更具体地说，本发明涉及一种基于精算迎风面积的LEO航天器大气阻力算法。
背景技术
低地球轨道（LEO：Low Earth Orbit）上运行的航天器，一般运行于距离地球表面高度小于2000km的范围。在LEO高度范围内，稀薄大气环境对航天器的作用虽然微弱但产生的影响却是实实在在的。对于高度范围在200km至500km内的超低轨道航天器，如果没有能量补充，一般都会在有限的时间内由于轨道高度衰减而再入稠密大气层陨落。诸如运行在高度400km左右的空间站等航天器，其轨道衰降至陨落的时间一般在几年以内。
绕地球运行的航天器轨道参数的变化受到许多方面因素影响，但是其高度的衰减主要由于大气环境的作用而导致。因此，当航天器运行高度不是特别高时是不能忽略大气环境作用的。在轨道动力学中，大气环境作用的基本参量即为航天器受到的气动阻力，在高度120km以上，空气动力学中认为此时处于自由分子流状态，因此，LEO航天器大气阻力的确定是非常重要的一个技术问题。
学术界对航天器大气阻力的各种计算方法手段统称为大气阻力模型，大气阻力模型是影响轨道确定及预报精度的主要因素之一。对于LEO航天器，大气阻力模型精度主要取决于大气密度、阻力系数以及有效迎风面积的精度。现有大气密度模型已比较精确并得到了普遍应用，如CIRA系列、DTM系列、Jacchia系列、MSIS系列等。另外，由于轨道环境中自由分子流假设成立，故可认为在轨航天器的阻力系数为常数。因此，有效迎风面积就成为确定大气阻力模型精度的一个关键因素。
目前，大气阻力建模时，迎风面积通过两种方式获得，一种方式认为迎风面积为常值；另一种是将迎风面积作为参数和轨道参数一起估计。对于复杂外形且姿轨运动形式变化多样的航天器，前一种精度太差，后一种不易求解。国内外航空领域学者曾提出了一些方法计算飞行器迎风面积，例如基于3DMAX和数据拟合技术、B样条曲线曲面造型和包含互斥理论、各类几何及拓扑变换模型等方法，均可有效解决计算精度不高的问题，但这些传统方法以飞机为对象，对于航天器其适用性有待于进一步检验。另外，也有学者针对航天器研究了考虑遮挡的投影面积求解问题，对于较简单的空间任务，可在一定程度上提高迎风面积的计算精度。但是针对未来复杂操作任务轨道的精确确定及预报，其精度还有待于进一步提高。
发明内容
本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。
为了实现本发明的这些目的和其它优点，提供了一种基于精算迎风面积的LEO航天器大气阻力算法，其特征在于，包括：
S1、在通过航天器数模生成表面网格后，基于网格面元与迎风截面格子投射的原理对航天器迎风面积S W 进行精算；
S2、将S1中得到的航天器迎风面积S W 引入至下式的LEO航天器大气阻力F Free 计算中：
其中，ρ为航天器飞行环境处的大气密度，根据不同大气模式计算或查询获得；V为来流速度亦即航天器飞行速度；C D0 为单位面积自由分子流阻力系数。
优选的是，在S1中，所述航天器迎风面积S W 的精算过程为：
S11、基于航天器数模生成表面网格，建立基准状态坐标系，计算航天器网格面元中心位置、面积及单位外法向量；
S12、计算航天器迎风截面上面元投射格子面积和尺度；
S13、根据航天器姿态角计算迎风截面投影范围，并根据面元投射尺度在迎风截面投影范围内划分规范化格子；
S14、遍历投影范围内格子及航天器网格面元记录投射信息，并据此统计计算迎风面积。
优选的是，在S11中，所述表面网格采用网格生成软件生成，且表面网格包括四边形结构网格和三角形非结构网格，航天器基准状态坐标系基于空气动力学建立；
所述航天器表面网格的面元中心位置（x,y,z）、面积ds及单位法向量基于平面解析几何算法得到，且所述航天器基准坐标系的纵向体轴和气流方向平行。
优选的是，在S12中，航天器迎风截面上面元投射格子面积和尺度的获取方法为：
S121、航天器迎风截面上面元投射格子面积为的计算公式如下：
其中，c为航天器迎风截面上面元投射格子尺度调节系数，N为航天器面元总数量，ds i 为序号i的面元面积；
S122、航天器迎风截面上面元投射格子尺度L C 的计算公式如下：
。
优选的是，在S13中，根据航天器姿态角计算迎风截面投影范围，并根据面元投射尺度在迎风截面投影范围内划分规范化格子的方法如下：
S131、在飞行器处于攻角α和侧滑角β时，任意基准面元纵向位置角变为α 0 +α，侧向位置角变为β 0 +β，则面元中心位置（x,y,z）的计算公式为：
其中，r 0 为面元距坐标原点或重心的距离，其表达式为：
（x 0 ,y 0 ,z 0 ）为基准状态下航天器面元中心位置坐标；
α 0 和β 0 分别为面元基准状态时的纵向位置角和侧向位置角，其表达式为：
其中，面元中心坐标Z分量为负时纵向位置角为正；当面元中心坐标X分量为零，则Y分量为正时侧向位置角为正90度，Y分量为负时侧向位置角为负90度；
S132、航天器迎风截面投影范围坐标值域分别用(Y max ,Y min )和(Z max ,Z min )表示，则其计算方法如下：
其中，下标j代表所有全部面元对应的坐标值，MAXVAL为取最大值，MINVAL为取最小值；
则航天器迎风截面所在平面等价于基准状态体轴系坐标系的YZ平面，迎风截面投影范围坐标值域等价于面元中心位置YZ坐标的极大值和极小值；
S133、定义坐标投影范围内在横轴和竖轴上的规范化坐标为(I Y ,I Z )，则一组(I Y ,I Z )坐标对应一个边长为L C 的正方形格子，所述(I Y ,I Z )的计算公式如下：
其中，INT（）为最小整数，一组面元坐标(y,z)对应一组规范化整数坐标(I Y ,I Z )，而一组(I Y ,I Z )能对应多组(y,z)；
S134、规范化坐标(I Y ,I Z )的值域范围I Y,min、I Y,max、I Z,min、I Z,max由投影范围边界确定，其计算方法如下：
S135、投影范围内规范化格子的数量N YZ 的计算公式如下：
。
优选的是，在S14中，所述迎风面积的计算方式为：
S141、将投影范围内每个格子的面元投影信息通过下述数据结构来表示：
其中，I X ，I Y 为格子规范化坐标，N jXY 为某个格子内投影的面元数量，
I J （j）和K J （j）分别为格子内投影的面元序号和单位外法向量状态，i XY 为格子序号；
S142、在进行投射信息记录前，遍历所有格子并初始化，以使所有N jXY 为0，I J （j）和K J （j）为空数组；
S143、遍历所有格子和面元，记录每个格子中投射的面元数量N jXY 、面元序号I J （j）、单位外法向量状态K J （j），则迎风面积为S W 的计算公式为：
其中，N XY 为投影范围内N jXY 不为0的所有格子数量，为每个格子的面积。
优选的是，在航天器的在轨动态飞行中，对任意资态下的气动阻力进行计算时，需基于圆球基准模型网格尺度变化和理论值标定，对不同精度指标的迎风面积计算进行如下的网格尺度划分：
精度指标96%时，L grid ≈L ref/50；
精度指标98%时，L grid ≈L ref/100；
精度指标99%时，L grid ≈L ref/200；
精度指标99.7%时；L grid ≈L ref/500；
精度指标99.8%时；L grid ≈L ref/1000；
其中，特征尺度L ref为航天器长度、展宽、直径的最大值，L grid 为面元网格尺度；
在任意姿态中，攻角范围为-180°~180°，侧滑角范围为0°~180°。
本发明至少包括以下有益效果：
其一，利用现代航天工程中的航天器构型数模采用完全的数字化方法快速高效且精准地获得航天器飞行过程中的迎风面积参数。
其二，针对航天器飞行过程中，无论处于任何姿态均可以采用本发明算法精确计算迎风面积。
其三，根据精算迎风面积从而较为快速准确地计算LEO航天器大气阻力。
本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
图1为基于精算迎风面积的LEO航天器大气阻力算法流程图
图2为航天器SCSS数模及基准坐标系；
图3为航天器SCSS表面三角形网格面元示意图；
图4为航天器SCSS在纵向0°、横向0°情况下，其迎风截面投影示意图；
图5为航天器SCSS在纵向90°、横向0°情况下，其迎风截面投影示意图；
图6为航天器SCSS在纵向0°、横向90°情况下，其迎风截面投影示意图；
图7为航天器SCSS在纵向45°、横向45°情况下，其迎风截面投影示意图；
图8为典型高度及姿态角条件下航天器SCSS大气阻力变化曲线。
具体实施方式
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
航天器往往外形复杂、姿态变化范围大，运行过程中其有效迎风面积精确求解困难，会影响大气阻力建模精度。为此，本发明基于现代航天器工程中不可或缺的数模为基础，采用物面网格面元和迎风截面格子投射技术精算航天器迎风面积，并以此为基础构建一种航天器大气阻力的快速高效算法。具体来说，本发明提供的一种基于精算迎风面积的LEO航天器大气阻力算法，在航天工程设计加工高度数字化条件下，基于航天器数模生成表面网格，应用网格面元-迎风截面格子投射技术，利用构思精巧的纯粹算法实现对迎风面积的精准计算，从而快速准确获得LEO航天器大气阻力。本发明对LEO航天器大气阻力的计算可以是航天器动态飞行处于任意姿态下的算法结果，是一种快速高效且具有高精度的算法。
本发明提供的一种基于精算迎风面积的LEO航天器大气阻力算法，如图1流程图所示，具有以下五个步骤：
步骤一、构建或应用航天器数模生成表面网格并建立航天器纵向体轴和气流方向平行的基准坐标系，计算航天器网格面元在基准坐标系的位置、面积及外法向量
航天器数模即航天器设计、加工或运行中常用的外形构型的三维计算机辅助设计（CAD：Computer Aided Design）电子文件，这些电子文件的来源建模软件包括但不限于ProE、UG、Solidworks、CATIA等。
采用网格生成软件生成航天器表面网格，这些网格生成软件包括但不限于HyperMesh、Pointwise、GridStar、Gmsh等。
针对航天器数模生成的表面网格包括四边形结构网格或三角形非结构网格两类；四边形结构网格文件采用PLOT3D格式（行业通用结构化数据格式，源于NASA），三角形非结构网格文件采用UCD格式（Unstructured Cell Data，行业通用非结构网格文件格式）。
建立航天器纵向体轴和气流方向平行的基准坐标系，且航天器基准坐标系的纵向体轴和气流方向平行，具体方式如下：
坐标原点位于航天器姿态旋转中心位置（一般定为重心），记为原点O，坐标为(0,0,0)；根据几何关系易知，如旋转中心变化，并不影响迎风面积量值，影响的是迎风截面投影区域在投影平面上的整体平移位置。故本发明算法中设定旋转中心不变，不影响最终算法结果。
坐标X轴为航天器主纵轴，其正方向与气流方向相反，由重心指向航天器头部或前部；
坐标Y轴为航天器横轴，其正向由重心指向航天器右侧；
坐标Z轴根据坐标X轴和坐标Y轴按右手法则确定，具体其正向由重心指向航天器下方。
基准坐标系建立后不再变化，航天器姿态变化后的面元坐标等信息均在这个基准坐标系中进行描述；航天器的迎风截面参数等价于与气流方向垂直的基准坐标系的YZ平面内参数，故本发明中采用航天器气流下游任意位置YZ平面内参数作为迎风面积计算的重要基础。
对于四边形面元，设4个角点M 1 M 2 M 3 M 4的坐标为 (x i ,y i ,z i ,),i=1,2,3,4于是四边形面元中心位置、面积及单位外法向量的计算方法如下：
四边形面元的中心点坐标设为(x,y,z,)，其坐标值的计算公式为：
其中，面元中心各坐标值即为四边形4个角点相应坐标的平均值；
设四边形面元的单位外法向量为，其表达式如下：
其中，分别是3个坐标方向上的单位向量；n x 、n y 、n z 是3个坐标方向上的单位外法向量分量，其表达式如下：
其中，N x 、N y 、N z 为四边形对角线构成的两个向量叉乘得到的结果向量在3个坐标方向上分量，其表达式如下：
设四边形面元的面积为ds，其计算公式如下：
其中，a 1为点M 1到M 2的距离，a 2为点M 1到M 4的距离，b 1为点M 2到M 3的距离，b 2为点M 4到M 3的距离，c为点M 1到M 3的距离。
其中，两个中间变量s 1和s 2的计算公式如下：
设三角形面元ΔABC的3个角点坐标分别为，于是三角形面元中心位置、面积及单位外法向量的计算方法如下：
三角形面元的中心点坐标设为，其坐标值的计算公式为：
其中，三角形面元中心各坐标值为三角形3个角点相应坐标的平均值；
设三角形面元的单位外法向量为，其计算公式如下：
其中，分别是3个坐标方向上的单位向量；n x 、n y 、n z 是3个坐标方向上的单位外法向量分量；其中，三角形两个邻边构成的向量/>的表达式为：
设三角形面元的面积为ds，其计算公式如下：
步骤二、根据航天器网格面元面积统计信息计算迎风截面投影格子面积和格子尺度
定义航天器迎风截面投影格子面积，其计算公式如下：
其中，c定义为格子尺度调节系数，其值大小决定了一个格子可能投射的面元数量，c值过小可能导致格子无投射面元而丢失应计入面积，c值过大则显得粗糙可能导致计算精度降低；根据经验推荐1≤c≤3；
其中，N为航天器面元总数量，ds i 为序号i的面元面积。
定义航天器迎风截面投影格子尺度L C ，其计算公式如下：
其中，航天器迎风截面投影格子尺度即迎风截面投影格子面积相应的正方形边长。
步骤三、计算航天器任意姿态角的迎风截面投影范围坐标极限，并采用格子尺度在迎风截面投影范围内划分规范化格子
设航天器处于任意攻角a和侧滑角β时，面元纵向位置角变为a 0 +a，侧向位置角变为β 0 +β。此时面元中心位置(x,y,z)的计算公式为：
其中，a和β分别为航天器的攻角和侧滑角。
上述式子中，r 0为面元距坐标原点即旋转中心的距离，即面元矢径，在前述基准坐标系下任意一个面元在航天器姿态变化过程中保持不变，其表达式为：
其中，(x 0,y0,z 0)为基准状态下航天器表面网格面元中心位置坐标。
上述式子中，a 0和β 0分别为面元基准状态时（攻角和侧滑角均为零）的纵向位置角和侧向位置角，其表达式如下：
其中，面元中心坐标Z分量为负时纵向位置角为正；当面元中心坐标X分量为零，则Y分量为正时侧向位置角为正90度，Y分量为负时侧向位置角为负90度。上述式子同样适用于面元处于任意位置处。
设航天器迎风截面投影区域坐标值域分别用(Y max ,Y min )和(Z max ,Z min )表示，则其计算方法如下：
其中，下标j代表所有全部面元对应的坐标值；
上述式子中，表明航天器迎风截面所在平面即等价于基准状态体轴系坐标系的YZ平面，迎风截面投影区域坐标值域即等价于面元中心位置YZ坐标的极大值和极小值。
根据格子尺度在迎风截面投影极限范围值域内划分规范化格子的方法如下：
定义坐标投影范围内在横轴和竖轴上的规范化坐标为(I Y ,I Z )，则一组(I Y ,I Z )坐标对应一个边长为L C 的正方形格子。
面元投影对应的(I Y ,I Z )由面元中心位置在横轴和竖轴上的坐标分量(y,z)转化而来，其计算方法如下：
其中， INT（）函数的意义为对括号中取最小整数。一组面元坐标(y,z)对应一组规范化整数坐标(I Y ,I Z )，一组(I Y ,I Z )可能对应多组(y,z)。
投影值域内规范化坐标(I Y ,I Z )的值域范围由投影区域边界确定，计算方法如下：
设计算投影值域内规范化格子的数量为N YZ ，计算公式如下：
由此即完成采用格子尺度在迎风截面投影极限范围值域内划分规范化格子的过程。
步骤四、历迎风截面投影范围内格子及航天器网格面元记录投射信息，根据面元-格子投射信息统计计算迎风面积
首先定义投影值域内每个格子的面元投影信息以数据结构表示，I X ，I Y 为格子规范化坐标，N jXY 为某个格子内投影的面元数量，I J （j）和K J （j）分别为格子内投影的面元序号和单位外法向量状态，i XY 为格子序号。
进行投射信息记录前，遍历所有格子初始化，所有N jXY 为0，I J （j）和K J （j）为空数组。
遍历所有格子和面元，记录每个格子中投射的面元数量N jXY 、面元序号I J （j）、单位外法向量状态K J （j），其中，单位外法向量状态取值的规则如下：
若n x ＜0，则K J （j）=1，推断航天器面元处于迎风状态；
若n x ＞0，则K J （j）=-1，推断航天器面元处于背风状态；
若n x =0，则K J （j）=0，推断航天器面元处于顺风状态。
设航天器迎风面积为S W ，其计算公式为：
其中，N XY 为投影范围内N jXY 不为0的所有格子数量，其值通过数据结构进行逐一判识计数N jXY 不为0状态而得；/>为每个格子的面积。
步骤五、根据精算的航天器迎风面积计算航天器大气阻力
设航天器大气阻力为F Free ，其计算公式如下：
其中，ρ为航天器飞行环境处的大气密度，根据不同大气模式计算或查询获得；V为来流速度亦即航天器飞行速度；C D0 为单位面积自由分子流阻力系数，根据分子碰撞理论计算获得，一般取值为2.2左右。
在上述航天器大气阻力F Free 计算公式中，在其它参数都确定情况下，航天器迎风面积S W 的精确计算成为决定性因素；航天器飞行姿态可以通过观测手段（包括内测或外测）进行确定，从而获得精算迎风面积的前提条件。
本发明是一种快速高效的LEO航天器在轨动态运行气动阻力的高精度算法，其中，所述在轨动态运行气动阻力可以是任意姿态下的算法结果；根据圆球基准模型网格尺度变化和理论值标定给出了迎风面积计算不同精度要求的具体网格尺度划分要求如下：
航天器自由分子流阻力计算的动态飞行，可以是任意姿态，具体为：攻角范围-180°~180°，侧滑角范围0°~180°。
设航天器特征尺度L ref为其长度、展宽、直径的最大值，航天器迎风面积的计算精度由主要面元网格尺度L grid 决定，根据圆球不同网格尺度面元划分迎风面积计算结果及理论值比较标定，具体的精度指标要求下的网格尺度参考取值如下：
精度指标96%时，L grid ≈L ref/50；
精度指标98%时，L grid ≈L ref/100；
精度指标99%时，L grid ≈L ref/200；
精度指标99.7%时；L grid ≈L ref/500；
精度指标99.8%时；L grid ≈L ref/1000。
实施例：
本发明的一种基于精算迎风面积的LEO航天器大气阻力算法，为了更清楚地说明本发明的技术方案，通过实施例进行说明。实施例采用一个类似中国空间站的简化模型（后续简称为SCSS：Simple China Space Station）作为例子，针对典型的运行高度和不同姿态情况，应用本发明描述的方法对其自由分子流阻力的计算过程进行详细描述。具体算法过程和结果如下。
步骤一、构建或应用航天器数模生成表面网格并建立航天器纵向体轴和气流方向平行的基准坐标系，计算航天器网格面元在基准坐标系的位置、面积及外法向量
采用建模软件生成的SCSS三维数模及基准坐标系如图2所示，该航天器为一种模块化组合体，在基准坐标系基准姿态条件下其纵向尺度为19米，横向尺度约20.5米，垂向尺度约16.7米；模型表面总面积约为700平方米。
在实施例中，SCSS基准坐标系的坐标原点为中部六面体核心舱的中心位置，此即航天器姿态变化的旋转中心。实际运行过程中改变姿态时，航天器的旋转中心可能会有明显变化，但是根据几何关系易知，即使旋转中心变化，并不影响迎风面积量值，影响的是迎风截面投影区域在投影平面上的整体平移位置。故本发明算法中设定旋转中心不变，不影响最终算法结果。基准坐标系实际即为航天器姿态归零情况下的体轴系，因此称之为基准坐标系。该基准状态坐标系建立后不再变化，后续航天器姿态变化后的面元坐标法向量等信息均在这个基准坐标系中进行描述。
采用网格软件生成的SCSS航天器表面非结构网格如图3所示，图中给出了一个局部放大图。根据统计，本实例划分的结构化四边形网格面元数量为864280个。
步骤二、根据航天器网格面元面积统计信息计算迎风截面投影格子面积和格子尺度
对实施例SCSS表面网格的上述864280个面元逐个计算其中心位置、面积和单位外法向量。需要注意的是，除了面元面积保持不变，随着航天器姿态变化，其每个面元的中心位置和单位外法向量都是动态变化的，这也是导致迎风面积随姿态不断变化的几何原因。
根据上述面元的基本信息，计算得出实施例航天器SCSS的迎风投影格子面积约为4051平方毫米，对应的迎风投影格子尺度约为63.6毫米。在本例计算中，格子尺度调节系数c=2；一般情况下要求c≥1，建议1≤c≤3；取值过小可能丢失格子，但取值过大则会影响计算精度。
步骤三、计算航天器任意姿态角的迎风截面投影范围坐标极限，并采用格子尺度在迎风截面投影范围内划分规范化格子
针对航天器SCSS任意姿态角计算迎风面积的中间步骤之一是获得投影范围的坐标值域，基于空气动力学中攻角（绕Z轴转角）和侧滑角（绕Y轴转角）的对应关系，实施例给出了几个典型姿态角组合下的投影区域坐标范围实际计算结果。具体计算结果如表1所示。
表1
根据上述投影范围值域，结合航天器SCSS迎风投影格子尺度在迎风截面投影区域坐标范围内划分规范化格子，这些格子都是规范化的正方形，它们的规范化坐标可以用整数进行表示。具体计算结果如表2所示，表中给出了规范化格子的值域范围，同时对格子数量进行了统计，需要了解的是这些格子并不一定都会被面元进行投射。
表2
步骤四、遍历迎风截面投影范围内格子及航天器网格面元记录投射信息，根据面元-格子投射信息统计计算迎风面积
面元-格子投射信息的记录，需要遍历投影范围内的所有格子。首先对属于每个格子的投射信息数据结构中的投影面元数量、序号数组及法向量状态数组初始化，然后再同时遍历格子和面元进行投射匹配同时更新投射信息。
图4中给出了表2所示几个典型姿态下航天器SCSS及其迎风截面投影示意图。
本实施例航天器SCSS的迎风面积在几个典型姿态条件下的统计计算结果如表3所示。本实施例航天器SCSS具有近百万量级的面元数量条件下，用主流PC每个状态的计算时间约为几秒钟，算是非常快速了。此外，表中还给出了SCSS在不同姿态下在基准坐标系中三个维度上的占位尺度，这些都是一些较为重要参数。
表3
步骤五、根据精算的航天器迎风面积计算航天器大气阻力
在航天器SCSS自由分子流大气阻力计算时，采用单位球体构型自由分子流阻力系数的面元法理论值2.23；根据上述航天器SCSS典型姿态下的迎风面积，在载人航天飞行器的运行和轨道衰降过程的典型高度范围内，应用标准大气模式给出的大气密度，设航天器飞行速度约为7500米/秒，计算出航天器SCSS的自由分子流大气阻力随高度的变化情况如表4-7所示，表中自由分子流大气阻力栏的序号和前述表1至表3中第一列的序号对应，表示不同的运行姿态状态。
表4
同时，在图8中绘制了相应的SCSS大气阻力变化曲线，图中同时还给出了对应高度的环境大气密度参数。从航天器SCSS自由分子流大气阻力随高度变化情况可知，它们的数值与大气密度和迎风面积同步变化；在大气模式和单位阻力系数确定的情况下，对迎风面积的精准计算是获得自由分子流大气阻力的重要前提条件。
实例中SCSS的表面网格平均尺度毫米，特征尺度取为L ref=20.5米，两者比值约为0.2‰，对应的迎风面积计算精度超过本发明用圆球理论值标定的99.8%，具有很高精度。/>
以上方案只是一种较佳实例的说明，但并不局限于此。在实施本发明时，可以根据使用者需求进行适当的替换和/或修改。
这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言，可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (3)

1.一种基于精算迎风面积的LEO航天器大气阻力算法，其特征在于，包括：
S1、在通过航天器数模生成表面网格后，基于网格面元与迎风截面格子投射的原理对航天器迎风面积S W 进行精算；
S2、将S1中得到的航天器迎风面积S W 引入至下式的LEO航天器大气阻力F Free 计算中：
其中，ρ为航天器飞行环境处的大气密度，根据不同大气模式计算或查询获得；V为来流速度亦即航天器飞行速度；C D0 为单位面积自由分子流阻力系数；
在S1中，所述航天器迎风面积S W 的精算过程为：
S11、基于航天器数模生成表面网格，建立基准状态坐标系，计算航天器网格面元中心位置、面积及单位外法向量；
S12、计算航天器迎风截面上面元投射格子面积和尺度；
S13、根据航天器姿态角计算迎风截面投影范围，并根据面元投射尺度在迎风截面投影范围内划分规范化格子；
S14、遍历投影范围内格子及航天器网格面元记录投射信息，并据此统计计算迎风面积；
在S11中，所述表面网格采用网格生成软件生成，且表面网格包括四边形结构网格和三角形非结构网格，航天器基准状态坐标系基于空气动力学建立；
所述航天器表面网格的面元中心位置（x,y,z）、面积ds及单位法向量基于平面解析几何算法得到，且所述航天器基准坐标系的纵向体轴和气流方向平行；
在S12中，航天器迎风截面上面元投射格子面积和尺度的获取方法为：
S121、航天器迎风截面上面元投射格子面积为的计算公式如下：
其中，c为航天器迎风截面上面元投射格子尺度调节系数，N为航天器面元总数量，ds i 为序号i的面元面积；
S122、航天器迎风截面上面元投射格子尺度L C 的计算公式如下：
；
在S14中，所述迎风面积的计算方式为：
S141、将投影范围内每个格子的面元投影信息通过下述数据结构来表示：
其中，I X ，I Y 为格子规范化坐标，N jXY 为某个格子内投影的面元数量，
I J （j）和K J （j）分别为格子内投影的面元序号和单位外法向量状态，i XY 为格子序号；
S142、在进行投射信息记录前，遍历所有格子并初始化，以使所有N jXY 为0，I J （j）和K J （j）为空数组；
S143、遍历所有格子和面元，记录每个格子中投射的面元数量N jXY 、面元序号I J （j）、单位外法向量状态K J （j），则迎风面积为S W 的计算公式为：
其中，N XY 为投影范围内N jXY 不为0的所有格子数量，为每个格子的面积。
2.如权利要求1所述的基于精算迎风面积的LEO航天器大气阻力算法，其特征在于，在S13中，根据航天器姿态角计算迎风截面投影范围，并根据面元投射尺度在迎风截面投影范围内划分规范化格子的方法如下：
S131、在飞行器处于攻角α和侧滑角β时，任意基准面元纵向位置角变为α 0 +α，侧向位置角变为β 0 +β，则面元中心位置（x,y,z）的计算公式为：
其中，r 0 为面元距坐标原点或重心的距离，其表达式为：
（x 0 ,y 0 ,z 0 ）为基准状态下航天器面元中心位置坐标；
α 0 和β 0 分别为面元基准状态时的纵向位置角和侧向位置角，其表达式为：
其中，面元中心坐标Z分量为负时纵向位置角为正；当面元中心坐标X分量为零，则Y分量为正时侧向位置角为正90度，Y分量为负时侧向位置角为负90度；
S132、航天器迎风截面投影范围坐标值域分别用(Y max ,Y min )和(Z max ,Z min )表示，则其计算方法如下：
其中，下标j代表所有全部面元对应的坐标值，MAXVAL为取最大值，MINVAL为取最小值；
则航天器迎风截面所在平面等价于基准状态体轴系坐标系的YZ平面，迎风截面投影范围坐标值域等价于面元中心位置YZ坐标的极大值和极小值；
S133、定义坐标投影范围内在横轴和竖轴上的规范化坐标为(I Y ,I Z )，则一组(I Y ,I Z )坐标对应一个边长为L C 的正方形格子，所述(I Y ,I Z )的计算公式如下：
其中，INT（）为最小整数，一组面元坐标(y,z)对应一组规范化整数坐标(I Y ,I Z )，而一组(I Y ,I Z )能对应多组(y,z)；
S134、规范化坐标(I Y ,I Z )的值域范围I Y,min、I Y,max、I Z,min、I Z,max由投影范围边界确定，其计算方法如下：
S135、投影范围内规范化格子的数量N YZ 的计算公式如下：
。
3.如权利要求1所述的基于精算迎风面积的LEO航天器大气阻力算法，其特征在于，在航天器的在轨动态飞行中，对任意资态下的气动阻力进行计算时，需基于圆球基准模型网格尺度变化和理论值标定，对不同精度指标的迎风面积计算进行如下的网格尺度划分：
精度指标96%时，L grid ≈L ref/50；
精度指标98%时，L grid ≈L ref/100；
精度指标99%时，L grid ≈L ref/200；
精度指标99.7%时；L grid ≈L ref/500；
精度指标99.8%时；L grid ≈L ref/1000；
其中，特征尺度L ref为航天器长度、展宽、直径的最大值，L grid 为面元网格尺度；
在任意姿态中，攻角范围为-180°~180°，侧滑角范围为0°~180°。
CN202311161851.XA 2023-09-11 2023-09-11 一种基于精算迎风面积的leo航天器大气阻力算法 Active CN116894302B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311161851.XA CN116894302B (zh) 2023-09-11 2023-09-11 一种基于精算迎风面积的leo航天器大气阻力算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311161851.XA CN116894302B (zh) 2023-09-11 2023-09-11 一种基于精算迎风面积的leo航天器大气阻力算法

Publications (2)

Publication Number Publication Date
CN116894302A CN116894302A (zh) 2023-10-17
CN116894302B true CN116894302B (zh) 2023-12-12

Family

ID=88315294

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311161851.XA Active CN116894302B (zh) 2023-09-11 2023-09-11 一种基于精算迎风面积的leo航天器大气阻力算法

Country Status (1)

Country Link
CN (1) CN116894302B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108820260A (zh) * 2018-05-04 2018-11-16 中国人民解放军63920部队 低轨航天器的中期轨道预报方法、装置、存储介质
CN115258197A (zh) * 2022-08-29 2022-11-01 北京航天飞行控制中心 航天器轨道终点的预测方法和装置、处理器及电子设备

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108820260A (zh) * 2018-05-04 2018-11-16 中国人民解放军63920部队 低轨航天器的中期轨道预报方法、装置、存储介质
CN115258197A (zh) * 2022-08-29 2022-11-01 北京航天飞行控制中心 航天器轨道终点的预测方法和装置、处理器及电子设备

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
应用阴影图的航天器迎风面积计算方法;杨成;唐歌实;李勰;陈光明;;计算机辅助设计与图形学学报(第11期);2155-2160 *
类天宫航天器迎风面积建模及变化特性分析;杨成;李勰;孙军;李志辉;;载人航天(第04期);430-435 *

Also Published As

Publication number Publication date
CN116894302A (zh) 2023-10-17

Similar Documents

Publication Publication Date Title
CN108052772A (zh) 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN107544067A (zh) 一种基于高斯混合近似的高超声速再入飞行器跟踪方法
Zhao et al. A viscous vortex particle model for rotor wake and interference analysis
CN108363843A (zh) 一种基于结构降阶模型的几何非线性静气动弹性全机配平方法
CN113761646B (zh) 一种用于移动风场环境中飞行器动响应的确定方法
CN111122899B (zh) 一种用于大气扰动中飞行的迎角侧滑角估计方法
CN106844887A (zh) 旋翼无人机的动力学建模方法及装置
CN104504255A (zh) 一种螺旋翼升力和阻力力矩的确定方法
CN114444216A (zh) 基于数值模拟的高空条件下飞行器姿态控制方法及系统
CN104931970B (zh) 一种面向机载气象雷达仿真的三维云场生成方法
Leng et al. Parameterized modeling and optimization of reusable launch vehicles based on reverse design approach
CN113468828B (zh) 一种飞机空中飞行颠簸强度指数计算方法
CN113408215B (zh) 一种用于移动风场环境中飞行器气动载荷的确定方法
CN116894302B (zh) 一种基于精算迎风面积的leo航天器大气阻力算法
Lobbia Optimization of waverider-derived crew reentry vehicles using a rapid aerodynamics analysis approach
Mocsányi et al. Grid and polytopic LPV modeling of aeroelastic aircraft for co-design
Olejnik et al. Aerodynamic modeling process using reverse engineering and computational fluid dynamics
Iljaszewicz et al. Aerodynamic analysis of the aircraft model made with the 3D printing method
McClain et al. Ice Accretion Roughness Variations on a Hybrid CRM65-Midspan Wing Model
Pueyo Efficient 3D artificial ice shapes simulations with 2D ice accretion codes using a 3-level correction
CN116894408B (zh) 采用数字化对风洞试验模型堵塞度进行计算的方法
CN116894301B (zh) 基于面元和格子投射的航天器迎风面积数字化获取方法
Quanwei et al. Dynamic model for real-time estimation of aerodynamic characteristics
CN112597588B (zh) 一种飞行器动态电磁散射特征快速建模方法
Ruiz Garcia et al. Aerodynamic Model Identification of the Flying V from Sub-Scale Flight Test Data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant