本发明涉及核材料辐照损伤模拟技术领域，具体涉及一种模拟粒子辐照损伤的时间步长优化方法。
背景技术：
材料遭受高能粒子(如中子、离子)辐照时，会造成“辐照损伤”。人们经常使用分子动力学(moleculardynamics,md)模拟研究辐照损伤过程。即给定初始构型和速度，采用一定的积分算法离散求解牛顿运动方程。
对于模拟低能(kev量级)粒子辐照损伤，因模型较小，采用较小的时间步长如一个飞秒即可快速模拟损伤过程。然而，模拟高能(10kev以上)粒子辐照损伤时，通常建立的模型原子数较多(几十万-上千万之间)，采用md模拟方法因计算较多的近邻原子对相互作用而较为耗时，加大计算时间步长虽能提高模拟的时间尺度，但是步长太大可能导致出现受力异常、速度异常等非物理意义现象的出现；降低模拟步长虽能保证物理过程正确，但要达到热平衡，须消耗大量计算资源。因此，需要开发一种根据粒子能量自动优化时间步长的算法，以高效模拟粒子辐照损伤过程。
技术实现要素：
为了解决以上问题，本发明提出一种模拟粒子辐照损伤的时间步长优化方法，以解决模拟粒子辐照损伤时因步长太大导致的受力异常、速度异常等非物理意义现象的出现，以及现有降低模拟步长方法计算量巨大的技术问题。
本发明是通过以下技术方案实现的：
本发明提供了一种模拟粒子辐照损伤的时间步长优化方法，适用于md模拟研究辐照损伤过程中牛顿积分方程中的时间步长的优化，包括以下步骤：
步骤s1：初始化，已知积分方程中的步长待优化体系中的粒子的速度和受力，设定初始时间步长、最大时间步长、相邻两步最大位移和动能最大变化量；
步骤s2：根据步骤s1的已知参数和设定参数，计算待优化体系中各粒子的位移改变量和动能改变量：
若位移改变量不大于设定的相邻两步最大位移，则将时间步长同设定的最大时间步长比较，若大于设定值，则取设定值为时间步长，若小于设定值，则该时间步长即为有效时间步长；若位移改变量大于设定的相邻两步最大位移，则将时间步长减半后重新计算；
若动能改变量不大于设定的动能最大变化量，则将时间步长同设定的最大时间步长比较，若大于设定值，则取设定值为时间步长，若小于设定值，则该时间步长即为有效时间步长；若位移改变量大于设定的动能最大变化量，则将时间步长减半后重新计算；
步骤s3：所有粒子计算完毕后，取步骤s2中获得的有效时间步长中最小者即为优化的步长。
优选地：
所述步骤s2中，位移改变量的计算方法为：
式中，dri为粒子i的位移改变量，vi为粒子i的速度，fi为粒子i的受力，ti为粒子i的时间步长，mi为粒子的质量。
所述步骤s2中，动能改变量的计算方法为：
式中，deki为粒子i的动能改变量，vi为粒子i的速度，fi为粒子i的受力，ti为粒子i的时间步长，mi为粒子的质量。
所述方法可同时适用于kev量级的低能粒子和10kev以上量级的高能粒子的辐照损伤模拟步长选择，尤其适用于后者。
本发明的有益效果是：本发明提出一种模拟粒子辐照损伤的时间步长优化方法，该方法通过求解粒子位移改变量、动能改变量与粒子速度和力的关系方程，能够根据粒子能量自动优化牛顿积分方程中的时间步长，从而实现高能粒子辐照损伤的高效快速模拟。同现有的模拟方法相比，给出了保证粒子在相空间轨迹近似连续的积分时间步长，提高了高能粒子辐照损伤模拟效率。
附图说明
图1为模拟粒子辐照损伤的时间步长优化方法的步骤流程图；
图2为模拟能量为20kev的铁初级碰撞原子诱导的辐照损伤过程中的时间步长、体系动能和势能变化曲线；
图3为优化的时间步长时模拟的铁表面20kev自离子辐照最大级联体积时的快照图。
具体实施方式
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。
实施例1
本实施例提供了一种模拟粒子辐照损伤的时间步长优化方法，包括以下步骤：
步骤s1：初始化，包括：
设定优化算法参数，包括：设定初始时间步长t0、最大时间步长tmax、相邻两步最大位移drmax和动能最大变化量dekmax。
读入md积分过程中某一步体系的粒子i的速度vi和受力fi。
步骤s2：根据位移改变量和动能改变量优化时间步长，使得md积分过程中，体系在相空间轨迹近似准连续变化，包括：
①根据位移改变量优化时间步长
步骤s201：计算粒子i位移改变量dri：
根据速度-verlet积分方法，粒子下一步位置r1与当前位置r、时间步长t、受力f、质量m有如下关系：
可知，位移改变量的计算方法为：
步骤s202：比较dri与drmax：若dri<drmax，再比较ti和tmax，若ti>tmax，则ti＝tmax，若ti<tmax，记录该ti值；若dri>drmax，则将该粒子i的时间步长ti减半，并转至步骤s201重新计算。
②根据位移改变量优化时间步长
步骤s211：计算粒子i动能改变量deki：
根据粒子下一步速度v1与当前速度v、时间步长t、受力f和质量m有如下关系：
动能改变量为：dek＝|mv12-mv2|；
得
步骤s212：比较deki与dekmax：若deki<dekmax，再比较ti和tmax，若ti>tmax，则ti＝tmax，若ti<tmax，记录该ti值；若deki>dekmax，则将该粒子i的时间步长ti减半，转至步骤s211重新计算。
步骤s3：所有粒子计算完毕后，取记录的时间步长值中的最小者为优化的时间步长。
采用本发明方法模拟能量为20kev的铁初级碰撞原子诱导的辐照损伤，测定模拟过程中的时间步长变化曲线、体系动能和势能变化曲线，结果如图2所示。
图2(a)显示了步长在0.001飞秒和5飞秒之间时间步长的动态变化，图中可看出，采用本方法模拟2000步时间尺度可达到1393.3飞秒，若是采用固定较小的时间步长如0.01飞秒，模拟相同的步数，所能达到的时间尺度为20飞秒，说明本算法较为高效。
图2(b)显示了模拟过程中体系动能和势能的变化，二者互为升降的曲线形态特征，说明在模拟中体系总能守恒，侧面地说明了本算法的正确性。
图3为采用本发明优化的时间步长时模拟的铁(100)表面20kev自离子辐照最大级联体积时的快照图，图中，纵横坐标轴长度单位为埃。
以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。