CN103218490B - 基于数值仿真的卫星防护结构弹道极限自动获取方法 - Google Patents

Abstract

本发明是一种基于数值仿真的卫星防护结构弹道极限自动获取方法，首先伺服程序根据要创建的仿真模型的关键参数和配置，生成建模脚本，调用前处理软件读取建模脚本执行建模操作，输出仿真模型k文件，启动LS-dyna求解器进行计算，并以固定时间间隔输出结果文件，然后利用LS-dyna求解器的重启动分析功能进行仿真终止判断，直到结果稳定，仿真结束，将得到的结果文件中的数据，以图形显示出，若得到临界穿透情况则直接获得极限点，否则采用自动极限直径搜寻得到极限点，最后获得所需要的极限点后，生成弹道极限曲线和方程。本发明方法降低了分析成本，缩短了周期，对仿真结果只需辅助少量实验进行校准，就能提供给工程设计使用。

Description

基于数值仿真的卫星防护结构弹道极限自动获取方法

技术领域

本发明涉及卫星关键部位的防护结构，具体涉及一种基于数值仿真的卫星防护结构弹道极限自动获取方法。

背景技术

为了抵御空间碎片对卫星的撞击，需要为卫星的关键部位设计并安装防护结构。空间碎片击中卫星防护结构会产生“穿透”与“未穿透”两种结果，而这两种状态的临界点就称为弹道极限。弹道极限是评估防护结构防护能力的最重要的指标。

工程上常用弹道极限曲线和弹道极限方程来刻画一种防护结构的防护能力。弹道极限曲线是指极限穿透下的碎片直径和撞击速度之间的关系曲线。而弹道极限方程是针对特定类型的防护结构又引入了其它参考变量，比如：防护板厚度、防护板间距等，通过数值拟合获得的方程。

空间碎片撞击卫星防护结构的过程在力学上被界定为超高速碰撞问题。针对这种问题开展地面试验研究是非常昂贵的。目前超高速碰撞数值仿真使用最广的数值算法是SPH（Smoothed Particle Dynamics，光滑粒子流体动力学）。LS-dyna是目前世界上使用最广的冲击动力学仿真软件之一。利用LS-dyna并结合前后处理器LS-prepost，用户可以建立空间碎片和卫星防护结构的SPH模型，并对碎片撞击防护结构的破碎过程进行SPH仿真计算。LS-dyna的仿真计算需要用户对撞击问题进行具体准确地描述，包括弹丸直径、撞击速度、防护结构具体布局、各部位尺寸和间距等。只有将这些和撞击条件相关的参数全部定义清楚之后才能进行撞击过程的仿真计算。

无论是弹道极限曲线还是弹道极限方程，其基础都是弹道极限点。有了极限点后通过绘图和拟合就能获得曲线和方程。在防护结构的几何尺寸和布局参数全都固定，碎片的撞击速度也固定的条件下，单一的变化碎片直径会获得“穿透”和“未穿透”两种结果。当撞击结果正好处于两者的临界状态时，此时就获得了极限状态下的一组参数组合，包括空间碎片直径和撞击速度等，称为一个极限点。由于空间碎片直径是最后一个被确定下来的量，它也被称为弹道极限直径。弹道极限直径就是指，空间碎片撞击防护结构处于临界点状态时候的空间碎片的直径。

然而通过现有LS-dyna的使用可知，数值仿真软件只能在全部参数都确定下的情况下才能计算，也包含碎片直径。因此现有数值仿真软件只能进行一个特定撞击过程的分析，而无法进行极限点分析。极限点的获取就需要人工建立多种碎片直径下的撞击模型，进行多个仿真计算，并通过一定策略人工寻找极限点。

另外由于人工进行分析时，每次对极限点的判断带有主观性，多次判断之间无法保证标准的客观性，因此分析得到的结论中含有人为误差，这降低了弹道极限分析结果的可信度。

发明内容

本发明的目的是为了解决目前人工分析获得卫星防护结构弹道极限，分析周期长，分析成本高，且还有人为误差的问题，提出一种基于数值仿真的卫星防护结构弹道极限自动获取方法。

本发明提供的一种基于数值仿真的卫星防护结构弹道极限自动获取方法，首先由用户建立工作目录和防护结构模型k文件，并将撞击仿真头文件和防护结构模型k文件拷入工作目录下，然后填写弹道极限自动获取配置文件，开始通过以下步骤进行弹道极限自动获取：

步骤一、自动建模：伺服程序根据当前要创建的仿真模型的关键参数和配置，生成建模脚本，之后调用前处理软件以后台批处理方式读取生成的建模脚本并执行建模操作，最后整理输出仿真模型k文件。所述的关键参数和配置包括：所要建立的空间碎片和防护结构的几何模型以及其物理特性，撞击条件和约束条件，仿真配置信息。所述的空间碎片采用弹丸模型。

步骤二、自动仿真计算：直接在DOS下以批处理方式启动LS-dyna求解器并以命令方式输入配置选项，对仿真模型k文件进行计算，并根据仿真模型k文件中的配置，以固定的时间间隔输出d3plot结果文件。

步骤三、自动仿真结果提取：伺服程序调用后处理软件以后台批处理方式读取d3plot结果文件，得到d3plot结果文件中的结果数据。所述的前处理软件和后处理软件都带有命令模式与后台批处理模式。

步骤四、自动仿真终止判断：对步骤三得到的结果数据进行判断，看结果是否稳定，如果结果稳定结束当前仿真，继续执行步骤五；如果结果还不稳定则继承该结果，通过LS-dyna求解器的重启动技术转步骤二执行。

步骤五、自动极限状态判断：读取当前步骤三得到的结果数据，以图形显示撞击过程和当前状态，如果空间碎片临界穿透防护结构，则直接获得极限点，转步骤七执行，否则继续执行步骤六。

步骤六、自动极限直径搜寻：建立方程：y=f(x)，其中，x是正实数，表示弹丸直径，y表示穿透状态，y的值是0或1，其中0表示未穿透，1表示穿透；用半值法求解函数f，给定初始值x1和x2，且x1<x2，并且f(x1)=0，f(x2)=1，令x=(x1+x2)/2，计算f(x)，如果f(x)=0则令x1=x，如果f(x)=1则令x2=x；再令x=(x1+x2)/2计算f(x)，以此类推，直到|x1-x2|<a，收敛公差a取为有限元网格模型的单元边长或者SPH无网格模型的粒子间距，此时，获得的弹丸的极限直径就是当前的(x1+x2)/2；其中，SPH表示光滑粒子流体动力学。

步骤七、在获得所需要数量的极限点后，调用工程数学分析软件matlab计算引擎进行绘图和拟合操作，生成弹道极限曲线和方程。

本发明的弹道极限自动获取方法，使用数值仿真对超高速碰撞问题进行数值模拟计算，相对传统人工分析，大大降低了分析的成本，缩短了获取的周期，并且在通过数值仿真方法获得了弹道极限特性后，只需要辅助以少量的实验，对仿真结果进行校准，就可以提供给工程设计使用，具有实用价值。

附图说明

图1是本方明弹道极限自动获取方法的步骤流程图；

图2是本发明方法步骤五中在弹丸撞击时防护结构在背面出现鼓包的示意图；

图3是本发明实施例中所建立的单层铝防护板和弹丸的示意图；

图4是本发明实施例采用步骤六的极限直径搜寻方法求解极限点过程的记录；

图5是对应图4中的记录的每一次迭代点的仿真结果图；

图6是本发明实施例得到的弹道极限曲线示意图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明一种基于数值仿真的卫星防护结构弹道极限自动获取方法，首先由用户建立工作目录和防护结构模型k文件，并将撞击仿真头文件和防护结构模型k文件拷入工作目录下，然后填写弹道极限自动获取配置文件，开始通过如图1所示的步骤进行弹道极限自动获取。弹道极限自动获取配置文件中，给定和极限点搜寻相关的参数，包括：工作目录、相关k文件的文件名、自动极限直径搜寻的两个初值x1和x2、收敛误差a；另外由于极限分析程序后台调用前后处理软件LS-prepost以及LS-dyna求解器，因此还需要给出这些软件批处理接口程序的路径和程序名。步骤一至六在弹道极限自动获取中要被迭代和不断重复，所以通过编程由程序自动完成，本发明实施例采用C++语言编程实现。

本发明的弹道极限自动获取方法，如图1所示，具体包括以下步骤：

步骤一、自动建模。伺服程序根据当前要创建的模型的关键参数和配置，动态生成建模脚本，之后调用前处理软件以后台批处理方式读取生成的建模脚本并执行建模操作，最后输出模型k文件。

所述的关键参数和配置，主要包括：建立空间碎片和防护结构的几何模型，如几何尺寸、布局尺寸等，空间碎片和防护结构的物理特性如材料模型、材料参数等，撞击条件和约束条件如撞击速度、防护板边缘的固定等，仿真配置信息。仿真配置信息是指仿真过程如结果输出设置、迭代时间步控制等，和在仿真进行时必要的设置如CPU使用情况、内存分配等。

所述的空间碎片一般以圆球形的铝弹丸来代替，这是在超高速碰撞地面模拟实验中使用最多的方式。铝弹丸一般使用SPH算法进行计算。所述的防护结构则要根据具体要分析的问题进行建模。常见的防护结构有：单板或多板式防护结构、填充式防护结构、带有蜂窝夹层板的防护结构、带有放热层的防护结构等。防护结构要根据具体的部件特征选用SPH、有限元实体或者有限元壳单元等来计算。模型的建立首先要根据每一个部件如弹丸、防护板、填充层等的几何尺寸建立几何体，之后根据选择的算法不同或填充粒子或划分网格。由此也可见如果要更改撞击弹丸尺寸，就需要从建模环节重新开始。

目前人工建模通过人工鼠标键盘操作前处理软件实现。本发明方法中的自动建模就是要让前处理软件在无人工干预的情况下自动执行，包括了两个要点：1、如何替代前处理软件中的人工鼠标操作；2、如何自动启动前处理软件并执行那些操作。

本发明方法中为了实现第1点，通过在前处理软件的选取中选择带有命令模式的软件来实现。软件LS-prepost拥有自己的cmd命令流，并且可以通过SCRIPTO语言进行二次开发，这样就可以用对应的指令来替代人工建模操作。软件ANSYS可以使用APDL（ANSYSParametric Design Language，ANSYS参数化设计语言）进行批处理操作。软件Truegrid也拥有其自身的基于fortran语言。这三种软件都带有命令模式。为了实现第2点，前处理软件还必须具有后台批处理模式，也就是在DOS下，在完全没有图形界面的情况下进行软件启动、仿真模型的创建和软件退出。这三个前处理软件LS-prepost、ANSYS和Truegrid都是具有批处理模式的。所以本发明的方法在具体实现中可以采用这三种前处理软件。本发明在具体说明中是采用软件LS-prepost来说明的。

因为在弹道极限获取过程中，需要改变弹丸直径，改变仿真模型，使得弹丸直径和防护结构模型的范围发生变化，而仿真模型改变了，建模过程就要变化，相应的命令模式下的建模命令也要变化，批处理模式下的建模脚本也要变化，所以为了实现自动建模，就要求在弹丸直径变化后，建模脚本可以自动重写，并根据自动重写的建模脚本建立新的模型，这就需要一个全局的伺服程序。伺服程序就是用来接收当前要创建的仿真模型的关键参数和配置，动态生成建模脚本，之后调用前处理软件以后台批处理方式读取生成的脚本并执行建模操作，最后整合输出下一步计算用的仿真模型k文件。

所述的伺服程序，它实现的功能是根据当前的目标，比如要创建一个弹丸模型，或者提取某条结果数据，动态创建相应的cmd脚本，创建的cmd脚本中就包含了相应的操作命令如创建弹丸或者提取数据相应的命令以及参数，进而调用带有命令模式与后台批处理模式的前处理软件执行cmd脚本中的命令，这样就将原本由人工实现的功能转为自动实现。无论前处理还是后面步骤所述的后处理都是LS-prepost的功能（ANSYS和Truegrid也类似），当人工使用LS-prepost时，由人工点击鼠标来实现相应的目标功能。

所述的仿真模型k文件，包括三部分：第一部分是撞击仿真头文件，撞击仿真头文件是被预先设定的，其内针对仿真冲击碰撞问题提供了设定方法和设定参数。仿真配置信息主要包含在撞击仿真头文件中。在该部分中，还包括弹丸的part（部件）、弹丸模型算法和弹丸模型材料信息，因为这些信息一般是不变的，即使要变动弹丸材料，只在头文件中修改也很方便。所述的弹丸模型算法比如为SPH算法。

第二部分是防护结构模型文件，主要是防护结构信息。由于防护结构的形式有多种，比如有多板式、填充式等，在材料使用上又涉及均质板、编制材料、泡沫材料等。相应的防护结构建模方式也要根据具体情况来确定，因此永远无法寻找一种普适的参数化建模方法适用于任何一种防护结构建模。所以对于防护结构信息，由用户以k文件的方式提供一个防护结构模型k文件。所述的防护结构模型k文件中包含防护结构的几何信息、材料信息和防护结构模型的算法等。所述的防护结构模块的算法比如为SPH、有限元实体或者有限元壳单元等。

第三部分是动态生成参数k文件，记载每次仿真需要更新的参数。为了进行弹道极限分析，需要对同一个防护结构在不同撞击速度下，不同的撞击弹丸直径下进行多次仿真。在这个过程中保持不变的配置参数都存放在头文件中。而对于每次仿真需要更新的参数，比如撞击速度和弹丸直径则需要由伺服程序动态生成。另外，由于弹丸直径的变化要求弹丸模型必须重新建立，因此弹丸模型的节点单元信息也是由伺服程序调用软件LS-prepost动态创建的。这些可变量构成了动态生成参数k文件。

将上述三部分组合就构成了可提供LS-dyna仿真的完整的k文件。

步骤二、自动仿真计算。直接在DOS下以批处理方式启动LS-dyna求解器并以命令方式输入配置选项，对仿真模型k文件进行计算，并根据仿真模型k文件中的配置以固定的时间间隔输出d3plot结果文件。

仿真计算的环节就是要将仿真模型k文件提交给LS-dyna求解器进行冲击动力学数值仿真求解。在计算过程中，LS-dyna求解器会根据用户在仿真模型k文件中做出的配置以固定的时间间隔输出d3plot结果文件。所有的计算结果信息都储存在d3plot结果文件中。

人工方式下，一般通过图形界面设定计算工作路径，并制定计算用的k文件，以及相关的内存和CPU配置选项，最后点击运行，之后LS-dyna求解器就会启动并进行求解。本发明方法中自动仿真计算就要求跳过图形界面，直接在DOS下以批处理方式直接启动LS-dyna求解器，并以命令方式输入各种配置选项。同样的由于每次使用的仿真模型k文件不同，工作路径不同，LS-dyna的启动命令参数也就不一样，这也是需要动态生成的。

比如：用DOS命令启动D盘test文件夹下的仿真模型k文件—exam.k计算，LS-dyna求解器ls971.exe存放在E盘的prog文件夹下，使用2个CPU并行计算，则启动的DOS命令为：

d:

cdtest\

e:\prog\ls971i=exam.k ncpu=2

上面的三行就是启动LS-dyna对仿真模型exam.k进行计算的DOS命令，而存储计算结果的路径就是D:\test\。由于DOS命令本身也是字符串，所以可以由伺服程序以字符变量的方式动态生成。

步骤三、自动仿真结果提取。伺服程序具有动态生成LS-prepost脚本进行后处理的功能。本步骤中，伺服程序自动调用后处理软件以后台批处理方式读取d3plot结果文件，得到d3plot结果文件中的结果数据。

d3plot结果文件是一种二进制文件，由于LS-dyna的生产商未公开d3plot文件内部的数据存储规则，因此用户无法直接提取d3plot结果文件中的有效信息，需要通过后处理软件实现d3plot结果文件的解读。LS-prepost是LS-dyna的标准后处理软件，利用它可以读取d3plot结果文件、提取运动学和变形力学数据并可以对撞击破碎过程进行图形显示。同样，软件ANSYS与Truegrid也可以用来作为后处理软件。

步骤四、自动仿真终止判断。对步骤三得到的结果数据进行判断，看结果是否稳定，如果结果稳定结束当前仿真，执行步骤五；如果结果还不稳定则继承该结果，通过LS-dyna求解器的重启动技术转步骤二执行。

数值仿真是使用计算机数值计算手段对真实过程的模拟，因此在计算前需要用户设定需要模拟真实情况的时间（终止时间），比如由初始状态开始演化20微秒之后的状态，这里的20微秒就是终止时间。对于碰撞问题，这个终止时间的设定一方面应该覆盖整个撞击过程，计算结束时各个撞击体的状态应该趋于稳定，这样才能得到有效的结果；另一方面终止时间不能太长，因为计算量和它成正比，太长会显著影响计算耗时。这就带来一个矛盾：用户使用数值仿真软件就是为了了解撞击过程要持续多久，会发生什么样的破坏，那么在计算之前当然就无从得知终止时间怎么设定是合理的。因此一般情况下，总是先设定一个足够大的仿真终止时间。人工方式下，分析人员通过不断查看当前的状态结果，如果发现力学过程已经趋于稳定，则判断此时的结果已经“有效”，人工强制停止仿真计算。

本发明中利用LS-dyna的重启动技术来实现在仿真过程中抽取当前结果，并根据结果分析确定是否继续仿真的目的。LS-dyna的重启动技术是指：改变模型或加载以后，时间在原来基础上向后顺延，是前一个分析的继续，在本发明中以时间t为步长，在步骤二中针对初始仿真模型进行仿真，在t时刻终止，步骤三对仿真结果文件进行数据读取，在步骤四中对得到的结果数据进行判断，如果结果稳定则当前撞击问题仿真结束；如果结果还不稳定那就修改当前仿真模型，重新设定终止时刻为2t，重新启动LS-dyna求解器以t为初始时刻，继承t时刻的结果，转步骤二继续计算。以此类推直到结果稳定，则当前的撞击仿真结束。

本发明使用加速度准则来进行系统稳定判断。因为撞击过程也就是物体之间力的相互作用，有力就有加速度。因此首先以每一个物体（球、板等）的等效质心加速度为考察量，如果每一个物体的等效质心加速度都趋近于0，则认为撞击过程已经结束。

然而，通过实践，发现仅有加速度指标是不完备的，有几个存在的问题。首先，在撞击过程中，可能出现短暂的加速值趋近于0的状态，具体设每个物体的等效质心加速值小于预设的阈值A时，此时撞击过程并未结束。为了解决这个问题，对加速度指标的考察就不能只对一个时间点进行，需要对一个时间段内的每个物体的等效质心加速度进行评估，一段时间每个物体的等效质心加速度的平均值都小于预设的阈值A，则说明结果稳定，撞击结束。所述的阈值A根据经验可以设定为初始速度值的1/1000。

另外，除了可能出现加速度短暂的小值之外，还可能出现加速度的振荡波动，这时引发一段时间内的加速度矢量叠加结果为小量的情况，这就要求对加速度的方差或者模值进行评估，加速度的方差或者模值趋近于0时，具体是每一个物体的等效质心加速度的方差或者模值小于预设的阈值A时，则说明结果稳定，撞击结束。此外，考虑到模型的对称性，可能会引发在撞击速度的垂直平面内加速度矢量叠加为0的情况出现，为了解决这个问题，对于无论对称还是非对称问题，在考察弹丸的加速度时，总是取其1/4进行考察，而非整个弹丸。

最后，在多板防护结构问题中，有可能出现弹丸穿透了第一块防护板，而未运动到第二块防护板处时，可能出现较长时间的匀速运动的状态。这种情况用加速度准则就无法判断。需要引入速度和位置准则辅助判断。由初始仿真模型，可以通过所有粒子或单元的空间坐标确定模型所占的空间范围，在仿真过程中，如果出现加速度接近0的情况，则对弹丸撞击轴线方向的速度进行判断。如果轴线速度大于预设阈值B，而弹丸粒子当前的坐标还没有超出防护结构的空间范围，说明弹丸还在防护结构内部运动，这意味着撞击过程还没结束。如果轴线速度大于预设阈值B，而弹丸粒子当前坐标已经超出了防护结构范围，说明弹丸已经穿透的防护结构，此时可以终止计算，结果已经稳定。如果弹丸的加速度趋近于0，速度也趋近于0，并且弹丸还在防护结构的空间范围内，这说明弹丸未穿透防护结构，结果已稳定，也可以终止仿真。所述的阈值B可以设定为撞击初始速度的1/100。

步骤五、自动极限状态判断。

在停止仿真后，使用LS-prepost读取d3plot文件，可以图形化的显示撞击过程和当前状态。有图形显示结果可以了解空间碎片（弹丸）是否穿透了防护结构，穿透了多少层，穿孔是多大等直观信息。根据弹道极限的定义，如果弹丸穿透则说明弹丸直径超过了极限直径，如果未穿透则说明未达到极限直径。如果恰好临界穿透则直接获得极限点，当然这种情况在大部分时候不会直接出现。

具体穿透状态的判断方法是：在加速度趋近于0的前提下，如果弹丸速度还很大，大于阈值B，并且弹丸超出了防护结构范围，则已穿透；如果弹丸速度也趋于0，并且弹丸范围在防护结构范围内，则为未穿透。由于弹丸的撞击会引起防护结构的变形，并在背面出现鼓包。因此在防护结构空间范围确定时，除了要考虑初始空间位置外，还要考虑到鼓包引起的范围变化。因为出现鼓包而未穿透的情况是存在的。如图2所示，为若干层防护板，弹丸在撞击最后一层防护板时候，防护板的背面出现鼓包，而最后一层防护板鼓包的出现会扩大原本防护结构占据空间的范围，此时如果按弹丸超出初始防护结构空间范围的准则来判断，则已经出现穿透，而最后一层防护板出现鼓包但是未穿透的情况是存在的，这样前面的判断就错误了。为了防止这种情况出现，在防护结构空间范围计算时应当针对鼓包问题进行修正。考虑到材料的承力极限，鼓包即使出现也不会太大。一般来讲，将防护结构范围修正为初始防护结构空间范围附加上1-1.5倍的最后一层防护板的厚度。

步骤六、自动极限直径搜寻。

如果在步骤五中未能直接得到极限状态，则要根据当前已有的结果，调整弹丸直径重新计算。如何根据当前已有的结果（以穿透点和未穿透点）来搜寻极限点，这可以根据人为经验，也可以结合必要的数学方法。

人工方法得到弹道极限特性，是在获得了弹道极限点后将数据整理成表格，然后使用excel或matlab等工具进行曲线绘制或者方程拟合，最终得到弹道极限特性。

本发明方法中采用极限直径搜寻方法来获得弹丸极限特性，所述的极限直径搜寻方法是用当前的计算结果确定下一步计算用到的弹丸直径的策略。首先把整个建模、计算和穿透状态判断的过程看作一个封闭的系统。这个系统的输入是弹丸的直径，而输出是是否穿透的状态。把“穿透”看作“1”，把“未穿透”看作“0”，则问题可以转化为二值方程问题。也就是y=f(x)，其中，x是正实数，表示弹丸直径，y的值是0或1表示穿透状态，f是没有显式表达式的函数，需要通过数值计算获得。用半值法求解这个问题，给定初始的x1和x2，且x1<x2，并且确保初始的f(x1)=0（未穿透），f(x2)=1（穿透）；令x=(x1+x2)/2，计算f(x)，如果f(x)=0则令x1=x，如果f(x)=1则令x2=x；再令x=(x1+x2)/2……以此类推，直到|x1-x2|<a，a是一个足够小的阈值，一般取为有限元网格模型的单元边长或者SPH无网格模型的粒子间距，表明计算已经收敛。

此时获得的极限直径就是当前的(x1+x2)/2，误差精度就是±a/2。可以通过调整a的大小来平衡极限点精度和计算耗时。

步骤七、更新撞击速度和弹丸直径，重复上述步骤一至步骤六，获得所需要的极限点数量后，调用工程数学分析软件matlab计算引擎进行绘图和拟合操作，生成弹道极限曲线和方程。一般极限点的个数为7个，撞击速度范围是2-8km/s，每间隔1km/s通过上面步骤得到一个极限点。在人工方式下，获得了极限点之后是通过使用其它软件进行绘制曲线和方程拟合的操作。

Matlab是目前使用最广泛的工程数学分析软件之一，有非常完善的绘图和拟合工具。Matlab的常规操作方法是在matlab软件环境下通过命令输入或脚本提交的方式实现功能。本发明中可以通过伺服程序动态生成脚本文件的方式，调用执行matlab。因matlab提供了更直接的方式，就是将matlab的计算引擎（matlab engine）直接集成入C++的开发环境下，所以本发明也可以以一定的方式在C++下直接使用matlab语言编写matlab脚本，并调用matlab引擎进行绘图和拟合操作。

如图3所示，为本发明实施例一个10mm单层铝防护板弹道极限曲线的自动获取。

在该实施例中，以3km/s的弹道极限点获取为例，图4给出了求解极限点过程的记录，图5中对应图4中的记录给出每一个迭代点的仿真结果图。图5左边一列为俯视图，右边一列为侧视图。由图4的结果和图5可以看出，第一次迭代，弹丸穿透靶板，破坏程度很大，穿孔也很大，这是半值法的初始上界x2。第二次迭代给出了半值法初始下界x1的结果，很明显没有穿透，防护板变形也不大。第三次迭代的直径取为x1和x2的平均值，结果已经穿透，弹丸材料并没有超出防护板材料范围，但穿孔明显，用当前的直径替换上界x2。第四次迭代，继续取当前x1和x2的平均值，结果没有穿透，但是变形程度明显大于第一次迭代，用当前直径替换下界x1。继续使用半值法生成新的模型，进行第五次迭代，没有穿透，但是变形较大。

经过五次迭代，当前x2-x1=0.0375。由于提供的防护板模型中的SPH最小粒子间距为0.05cm，因此将弹道极限的收敛公差a定为0.05cm。可见经过五次迭代后已经满足收敛条件。最终取当前上界和下界的均值作为弹道极限点——0.33125cm。

本发明实施例以完全相同的方法对撞击速度从2km/s到8km/s，以1km/s为间隔进行遍历，于是最终获得了每一个特征速度下的弹道极限点。

d＝11.6·v^-0.5258   （1）

图6给出了弹道极限曲线，横轴为撞击速度v，单位km/s，纵轴为弹道极限直径d，单位mm。曲线为用最小二乘法回归得到的弹道极限曲线，其表达式如式（1）所示，为最小二乘方法回归得到的单板防护结构弹道极限方程。

Claims (6)

1.一种基于数值仿真的卫星防护结构弹道极限自动获取方法，其特征在于，首先由用户建立工作目录和防护结构模型k文件，并将撞击仿真头文件和防护结构模型k文件拷入工作目录下，然后填写弹道极限自动获取配置文件，开始通过以下步骤进行弹道极限自动获取：

步骤一、自动建模：伺服程序根据当前要创建的仿真模型的关键参数和配置，生成建模脚本，之后调用前处理软件以后台批处理方式读取生成的建模脚本并执行建模操作，最后整理输出仿真模型k文件；所述的关键参数和配置包括：所要建立的空间碎片和防护结构的几何模型以及其物理特性，撞击条件和约束条件，仿真配置信息；所述的空间碎片采用弹丸模型；

步骤二、自动仿真计算：直接在DOS下以批处理方式启动LS-dyna求解器并以命令方式输入配置选项，对仿真模型k文件进行计算，并根据仿真模型k文件中的配置，以固定的时间间隔输出d3plot结果文件；

步骤三、自动仿真结果提取：伺服程序调用后处理软件以后台批处理方式读取d3plot结果文件，得到d3plot结果文件中的结果数据；所述的前处理软件和后处理软件都带有命令模式与后台批处理模式；

步骤四、自动仿真终止判断：对步骤三得到的结果数据进行判断，看结果是否稳定，如果结果稳定，结束当前仿真，继续执行步骤五；如果结果不稳定，则继承该结果，通过LS-dyna求解器的重启动技术转步骤二执行；

步骤四中所述的结果是否稳定，具体是：(1)如果每一个物体的等效质心加速度在一个时间段内的平均值，都小于阈值A，则结果稳定；(2)若出现加速度振荡波动的情况，对加速度的方差或者模值进行计算，如果每一个物体的等效质心加速度的方差或者模值都小于预设的阈值A，则结果稳定；在撞击速度的垂直平面内加速度矢量叠加为0的情况，计算1/4弹丸的等效加速度，若1/4弹丸的等效加速度小于预设的阈值A，则结果稳定；(3)针对弹丸穿透了一块防护板而未运动到下一块防护板处，出现匀速运动状态的情况，对弹丸撞击轴线方向的速度进行判断，如果轴线速度大于阈值B，而弹丸粒子当前的坐标没有超出防护结构的空间范围，则结果不稳定，撞击过程没结束，如果轴线速度大于阈值B，而弹丸粒子当前坐标超出了防护结构范围，则结果稳定，撞击结束；如果弹丸的加速度趋近于0，速度也趋近于0，并且弹丸还在防护结构的空间范围内，则弹丸未穿透防护结构，但结果稳定，撞击结束；所述的阈值A设定为撞击初始速度的1/1000，所述的阈值B设定为撞击初始速度的1/100；

步骤五、自动极限状态判断：读取当前步骤三得到的结果数据，以图形显示撞击过程和当前状态，如果空间碎片临界穿透防护结构，则直接获得极限点，转步骤七执行，否则继续执行步骤六；

步骤六、自动极限直径搜寻：建立方程：y＝f(x)，其中，x是正实数，表示弹丸直径，y表示穿透状态，y的值是0或1，其中0表示未穿透，1表示穿透；用半值法求解函数f，给定初始值x1和x2，且x1<x2，并且确保初始的f(x1)＝0，f(x2)＝1，令x＝(x1+x2)/2，计算f(x)，如果f(x)＝0则令x1＝x，如果f(x)＝1则令x2＝x；再令x＝(x1+x2)/2计算f(x)，以此类推，直到|x1-x2|<a，收敛公差a取为有限元网格模型的单元边长或者SPH无网格模型的粒子间距，此时，获得的弹丸的极限直径就是当前的(x1+x2)/2；其中，SPH表示光滑粒子流体动力学；

步骤七、更新撞击速度和弹丸直径，重复步骤一到步骤六，获得所需要的极限点数量后，伺服程序调用工程数学分析软件matlab计算引擎进行绘图和拟合操作，生成弹道极限曲线和方程。

2.根据权利要求1所述的卫星防护结构弹道极限自动获取方法，其特征在于，所述的弹道极限自动获取配置文件，包含了与弹道极限点搜寻相关的参数，包括：工作目录、相关k文件的文件名、自动极限直径搜寻的两个初始值x1和x2、收敛误差a、以及前后处理软件以及LS-dyna求解器的批处理接口程序的路径和程序名。

3.根据权利要求1所述的卫星防护结构弹道极限自动获取方法，其特征在于，步骤一中所述的仿真模型k文件，包括三部分：撞击仿真头文件、防护结构模型k文件和动态生成参数k文件；所述的撞击仿真头文件是被预先设定的，针对仿真冲击碰撞问题提供了设定方法和设定参数，包含仿真配置信息以及弹丸的部件、弹丸模型算法和弹丸材料信息；所述的防护结构模型k文件包含防护结构的几何信息、材料信息和防护结构模型的算法；所述的动态生成参数k文件，由伺服程序动态生成，记载每次仿真需要更新的参数，包括撞击速度和弹丸直径。

4.根据权利要求1或2所述的卫星防护结构弹道极限自动获取方法，其特征在于，步骤一所述的前处理软件和步骤三所述的后处理软件为LS-prepost、ANSYS或者Truegrid。

5.根据权利要求1所述的卫星防护结构弹道极限自动获取方法，其特征在于，步骤五中所述的空间碎片临界穿透防护结构穿透的判断方法是：在加速度趋近于0的前提下，如果弹丸速度大于阈值B，并且弹丸超出了防护结构范围，则弹丸已穿透防护结构；如果弹丸速度也趋于0，并且弹丸在防护结构范围内，则弹丸未穿透防护结构；对于由于弹丸的撞击引起防护结构在背面出现鼓包的情况，将防护结构范围修正为初始防护结构范围附加上1-1.5倍的最后一层防护板的厚度。

6.根据权利要求1所述的卫星防护结构弹道极限自动获取方法，其特征在于，步骤七中所述的所需要的极限点数量为7个；撞击速度的范围是2-8km/s，每间隔1km/s得到一个极限点。