CN112580223B - 一种基于选区应力判据的双金属界面自适应加载模拟方法 - Google Patents

一种基于选区应力判据的双金属界面自适应加载模拟方法 Download PDF

Info

Publication number
CN112580223B
CN112580223B CN202011587163.6A CN202011587163A CN112580223B CN 112580223 B CN112580223 B CN 112580223B CN 202011587163 A CN202011587163 A CN 202011587163A CN 112580223 B CN112580223 B CN 112580223B
Authority
CN
China
Prior art keywords
stress
region
atom
criterion
iteration
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
CN202011587163.6A
Other languages
English (en)
Other versions
CN112580223A (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Publication of CN112580223A publication Critical patent/CN112580223A/zh
Application granted granted Critical
Publication of CN112580223B publication Critical patent/CN112580223B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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

Abstract

本发明公开了一种基于选区应力判据的双金属界面自适应加载模拟方法，属于分子静力学模拟领域。首先模拟双金属界面模型并进行区域划分，通过周期性边界条件得到扩展盒子，建立扩展盒子内各原子的近邻集合，并计算各原子的体积，统计特征区域内的应力。如果特征区域内上下层的应力未达到应力平衡，则针对不同加载模式，对双金属界面模型的可变区域施加局部应变，直至特征区域内的应力达到平衡。然后统计施加局部应变的迭代失败次数，直至总迭代失败次数大于用户预设的迭代次数，输出当前特征区域应力，通过CG和FIRE方法计算双金属界面模型最小化能量，更新原子坐标。本发明实现了多种体系状态的内应力平衡，得到双金属界面加载模拟的精确结果。

Description

一种基于选区应力判据的双金属界面自适应加载模拟方法
技术领域
本发明属于材料学中的分子静力学模拟领域，具体为一种基于选区应力判据的双金属界面自适应加载模拟方法。
背景技术
界面是晶体中的面缺陷，它对于晶体材料的性质和发生的转变有着极其重要的影响。对于位错，界面一方面能够阻碍滑移位错的运动，也能够分解、吸收位错；另一方面，界面也能作为位错等缺陷的源头，从而引起界面强化，提高材料的强度。
近年来，大量的计算模拟和实验研究发现：界面结构会直接影响位错的形核方式和界面的滑移方式，从而改变纳米复合物材料和异质结构材料的力学强度和塑性等。在对界面位错和模型变形进行分析之前，重中之重是构建出符合力学性质匹配的模型，对于层状复合材料一般认为界面法方向符合并联(Ruess)模型，如文献1：Reuss A.Berechnungder Flieβgrenze von Mischkristallen auf Grund der
Figure BDA0002867487250000011
fürEinkristalle[J].ZAMM Journal of applied mathematics and mechanics:Zeitschriftfür angewandte Mathematik und Mechanik,1929,9(1):49-58，公开了构成复合材料的两层应力相等；同时平行界面方向符合串联(Voigt)模型，如文献2：Voigt W.Ueber dieBeziehung Zwischen den beiden Elastizitatsconstanten isotroper Korper[J].Annalen der Physik und Chemie,1989,38.公开了构成符合材料两层应变相等。另外，为了绘制完整的应力应变图，一般需要让初始应力为零。
对于含有界面的异质模型而言，如何控制界面两侧的原子从而使模型达到平衡状态，并与实际界面具有相近的结构性质是一切后续研究的前提，但是在现有的相关模拟研究中一般采用均匀加载的方式，使得双金属界面材料在计算模拟准静态加载过程中，内应力无法达到对应加载模式应力分布的假设。
发明内容
本发明为了从根本上解决界面位错形核机理和变形机制前的模型构建问题，研究出了一种可应用于广泛的体系范围和模型结构的内应力平衡方法，具体为一种基于选区应力判据的双金属界面自适应加载模拟方法。
所述的基于选区应力判据的双金属界面自适应加载模拟方法，包括如下步骤：
步骤一，通过计算机模拟双金属界面模型，对双金属界面模型进行区域划分，具体划分为可变区域、特征区域与固定区域。
具体为：
首先，计算机预制双金属界面模型，以及与界面模型对应的势函数；
所述的双金属界面模型包括盒子信息、盒子内原子类型和所有原子的坐标。
盒子信息包括模拟正交盒子的对角线两个端点坐标以及盒子内原子总数；
对角线两个端点坐标分别为起点坐标(xs,ys,zs)和终点坐标(xe,ye,ze)；
盒子内第i个原子坐标为(xi,yi,zi)，用户指定界面法方向，界面法方向必须平行于直角坐标系中的其中一个轴，且模型界面在z＝0面上，界面层厚度大于
Figure BDA0002867487250000021
。
所述的势函数为分子动力学或静力学可用势函数，不同的势函数对应的截断半径不同。
然后，利用模型界面和势函数对双金属界面模型进行区域划分，得到可变区域、特征区域与固定区域：
可变区域分上下层，上层称为可变A1区，下层称为可变B1区，可变A1区为模型界面z>0的区域，可变B1区为z<0的区域。
特征区域分上下层，上层称为特征A2区，下层称为特征B2区，特征A2区为模型界面z>h1∩z<h2的区域，特征B2区为模型界面z>-h2∩z<-h1的区域；其中h1和h2分别为用户指定参数。
固定区域为距离盒子z方向边界h3以内的区域；h3为用户指定参数，且h3大于势函数的截断半径。
步骤二，通过周期性边界条件得到扩展盒子，建立扩展盒子内各原子的近邻集合，并计算各原子的体积，最后利用各原子的体积统计特征区域内上下层的应力：
具体为：
步骤201，针对模拟正交盒子内的原子，应用周期性边界条件选择与盒子边界的距离小于截断半径的原子复制到扩展盒子中；
具体方法是：
首先，对模拟正交盒子外的原子重整，将位于周期性边界外的原子移入模拟正交盒子中。
然后，将原有的模拟正交盒子沿边界向外扩展，形成包括原正交盒子在内的新的扩展盒子；
扩展长度为原有的模拟正交盒子的势函数对应的截断半径的长度。
接着，对位于原正交盒子内的原子逐一进行判断：将某个原子a在某一方向平移一个周期性边界对应的周期后，判断平移后的原子a是否到了扩展部分中。如果是，则将原子a复制到扩展部分中，并继承原子a平移之前的原子类型和原子坐标。否则，原子a仍在原有的模拟正交盒子中，不做处理。
步骤202，对扩展盒子内的各原子分别建立近邻集合；
具体为：
首先，利用各原子坐标分别除以截断半径并向下取整，作为该原子所在像素区域号；
然后，针对某个原子b，判断相邻像素区域内的原子c与原子b的距离是否小于原子b的截断半径，如果是，则将原子c添加到原子b的近邻集合中；否则，对原子c不做处理；依次对原子b的相邻像素区域内的原子进行判断，得到的原子b的近邻集合；
同理，得到扩展盒子内其余每个原子的近邻集合。
步骤203，针对扩展盒子内的第i个原子，将该原子与其近邻集合内各原子连线的垂直平分线构成闭包，通过voronoi方法计算闭包体积作为原子i的原子体积Vi。
同理，得到扩展盒子内每个原子的原子体积。
步骤204，通过扩展盒子内每个原子的体积和应力张量，计算特征区域的统计应力：
首先，计算每个原子的原子应力：
以坐标轴方向为参考，采取6个应力张量来表示，分别为：Sxx、Syy、Szz、Sxy、Syz、Sxz。
对于每一个应力张量，表达为：
Figure BDA0002867487250000031
其中，a、b分别代表x、y、z分量中的任意一个，m代表原子质量，v代表原子此刻的速度；Np代表原子的近邻数，ria表示第i个原子由于成键导致的原子力在a方向上的分量差，Fib表示第i个原子由于成键导致的原子力在b方向上的分量，i＝1,2,3,4；Nb代表原子所参与形成的键的数目，Na代表原子所参与形成的键角的个数，Nd代表原子所参与形成的二面角的个数，Ni代表第i个原子所参与形成的非正常二面角的个数，Nf代表施加在原子上的内约束力作用数；Kspace代表长程库伦相互作用。
然后，对于选取的模拟双金属界面模型的特征区域的上下层，分别利用各原子的体积计算特征区域上下层的统计应力；
即：
Figure BDA0002867487250000032
σ表示特征区域应力，N代表待测特征区域的原子数量。
步骤三，分别判断特征区域内上下层的应力是否达到应力平衡，如果是，进入步骤四；
否则，利用不同加载模式对可变区域施加局部应变，直至特征区域内的应力达到平衡，进入步骤四。
加载模式包括：初始加载模式、平行界面单轴加载模式、平行界面双轴加载模式和界面法向加载模式。
(1)对初始加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1011、计算初始加载模式下可变区域界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000041
Figure BDA0002867487250000042
Figure BDA0002867487250000043
其中，
Figure BDA0002867487250000044
表示特征A2区统计应力xx方向上的分量，
Figure BDA0002867487250000045
表示特征A2区统计应力yy方向上的分量，
Figure BDA0002867487250000046
表示特征A2区统计应力zz方向上的分量，
Figure BDA0002867487250000047
表示特征B2区统计应力xx方向上的分量，
Figure BDA0002867487250000048
表示特征B2区统计应力yy方向上的分量，
Figure BDA0002867487250000049
表示特征B2区统计应力zz方向上的分量。
步骤1012、利用特征A2区和特征B2区的目标应力状态计算特征区域的应力平衡判据；计算公式为：
Figure BDA00028674872500000410
步骤1013、用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1014、判断当前次迭代的应力平衡判据是否满足Dk,0＞Dmax，如果是，进入步骤1015；否则，特征区域内的应力达到平衡，结束算法；
步骤1015、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk，f＜Dk，0，将失败次数f置0，迭代次数k自增1，输出Dk，f作为应力平衡判据Dk+1，0，进入步骤1014重复判断；
具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA00028674872500000411
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000051
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000052
其中λ为系数常数，E为材料平均弹性模量的常数。
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1016、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(2)对平行界面单轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1021、以平行于x轴加载的可变区域为例，计算其两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000053
Figure BDA0002867487250000054
步骤1022、利用目标应力状态计算特征区域的应力平衡判据；
计算公式为：
Figure BDA0002867487250000055
步骤1023、用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1024、判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1025；否则，特征区域内的应力达到平衡，结束算法；
步骤1025、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1024重复判断；
对平行于x轴记载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000061
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000062
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000063
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1026、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(3)对平行界面双轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1031，以平行于xy轴加载模式的可变区域为例，计算界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000064
步骤1032，利用目标应力状态计算特征区域的应力平衡判据；
Figure BDA0002867487250000071
步骤1033，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1034，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1035；否则，特征区域内的应力达到平衡，结束算法；
步骤1035、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1034重复判断；
对平行于xy轴加载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000072
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000073
对盒子无需修正。
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1036、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(4)对界面法向加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1041，以平行于x轴界面法向加载模式的可变区域为例，其界面两侧的特征区域A2和B2的目标应力状态为：
Figure BDA0002867487250000081
Figure BDA0002867487250000082
Figure BDA0002867487250000083
步骤1042，利用目标应力状态计算特征区域的应力平衡判据：
Figure BDA0002867487250000084
步骤1043，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1044，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1045；否则，特征区域内的应力达到平衡，结束算法；
步骤1045、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1044重复判断；
具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000085
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000086
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000087
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1046、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
步骤四，输出当前特征区域应力平衡判据，通过共轭梯度算法以及快速松弛引擎算法对达到应力平衡的双金属界面模型计算最小化能量，并更新原子坐标，获得双金属界面在内应力平衡条件下的加载模拟结果。
本发明的优点在于：
(1)本发明一种基于选区应力判据的双金属界面自适应加载模拟方法，实现了多种体系状态下内应力平衡模型的建立，实现了指定加载方式特定应力条件下双金属界面的加载；
(2)本发明一种基于选区应力判据的双金属界面自适应加载模拟方法，实现了准静态加载的内应力平衡条件，能得到一般分子动力学模拟难以得到的精确结果。
附图说明
图1为本发明基于选区应力判据的双金属界面自适应加载模拟方法的原理图；
图2为本发明基于选区应力判据的双金属界面自适应加载模拟方法流程图。
具体实施方式
以下结合附图和实施例对本发明进行详细说明。
本发明提供的基于选区应力判据的双金属界面自适应加载模拟方法，如图1所示，首先读取盒子大小内的原子坐标，对模型划分出特征区域，统计特征区域内的原子应力，根据应力控制判据判断模型的各部分是否达到平衡，对未达到平衡状态的特征区域，通过对可变区域进行拉伸、压缩或剪切变形，循环调整直到模型达到特征区域的内应力平衡。
内应力平衡状态包括：未加载条件界面法向上的无正应力并联模型，界面方向上的无正应力串联模型；平行界面加载条件下界面法向上的无正应力并联模型，界面方向上串联模型；界面法向加载条件下界面法向上的并联模型，界面方向上的无正应力串联模型；剪切加载下剪切方向与界面法向上的并联模型，界面方向上的无正应力串联模型。
满足内应力平衡后，使当前模型达到符合力学性质匹配条件的状态，之后通过能量最小化原则进行原子坐标的调整(即材料驰豫过程)，便可以输出精确可靠的模拟结果，包括应力应变和加载后的原子坐标。
本发明不但可应用于初始模型构建过程中通过调整模型内应力，达到内应力平衡，还可以应用于准静态加载过程中，在模拟计算中得到精细准确的结果，且该结果可以修正常用模拟方法导致模拟过程中产生的误差，对后续的计算和分析有着不可替代的作用。
一种基于选区应力判据的双金属界面自适应加载模拟方法，如图2所示，包括如下步骤：
步骤一，通过计算机模拟双金属界面模型，对双金属界面模型进行区域划分，具体划分为可变区域、特征区域与固定区域。
具体为：
首先，从计算机获取预制的双金属界面模型，以及与预制的双金属界面模型对应的势函数；
所述的双金属界面模型包括盒子信息、原子类型和所有原子的原子坐标。
盒子信息包括模拟正交盒子的对角线两个端点坐标以及盒子内原子总数；
对角线两个端点坐标分别为起点坐标(xs,ys,zs)和终点坐标(xe,ye,ze)；起点坐标(xs,ys,zs)取所有原子坐标的最小值，终点坐标(xe,ye,ze)取所有原子坐标中的最大值；
第i个原子坐标为(xi,yi,zi)，由于是双金属界面模型，用户指定界面法方向，界面法方向必须平行于直角坐标系中的其中一个轴，以z轴为例，z轴可替换为x轴或y轴，且模型界面在z＝0面上，要求界面层够厚，界面层厚度大于
Figure BDA0002867487250000101
。
所述的势函数为分子动力学或静力学可用势函数，可以为EAM，MEAM和BOP等LAMMPS常见格式，不同的势函数对应的截断半径不同。
然后，利用模型界面和势函数对双金属界面模型进行区域划分，得到可变区域、特征区域与固定区域：
可变区域分上下层，上层称为可变A1区，下层称为可变B1区，可变A1区为模型界面z>0的区域，可变B1区为z<0的区域。
特征区域分上下层，上层称为特征A2区，下层称为特征B2区，特征A2区为模型界面z>h1∩z<h2的区域，特征B2区为模型界面z>-h2∩z<-h1的区域；
其中h1和h2分别为用户指定参数，一般要求
Figure BDA0002867487250000102
标记特征区域中的原子为特征区域原子，且特征区域原子不会随着模拟特征区域变化而变化。
固定区域为距离盒子z方向边界h3以内的区域；
h3为用户指定参数，且h3大于势函数的截断半径。固定区域不与特征区域重叠，固定区域内原子不会随着模拟过程中固定区域变化而变化，并且固定区域原子不参与共轭梯度算法(CG)以及快速松弛引擎算法(FIRE)能量最小化计算。
步骤二，通过周期性边界条件得到扩展盒子，建立扩展盒子内各原子的近邻集合，并计算各原子的体积，最后利用各原子的体积统计特征区域内的应力：
具体为：
步骤201，针对模拟正交盒子内的原子，应用周期性边界条件选择与盒子边界的距离小于截断半径的原子复制到扩展盒子中；
具体方法是：
首先，对模拟正交盒子外的原子重整，将位于周期性边界外的原子移入模拟正交盒子中。
然后，将原有的模拟正交盒子沿边界外扩展，形成包括原正交盒子在内的新的扩展盒子；
扩展长度为原有的模拟正交盒子的势函数对应的截断半径的长度。
接着，利用扩展盒子对位于原正交盒子内的原子逐一进行判断：将某个原子a在某一方向平移一个周期性边界对应的周期后，判断平移后的原子a是否到了扩展部分的盒子中。如果是，则将原子a复制到扩展部分的盒子中，并继承原子a平移之前的原子类型和原子坐标。否则，原子a仍在原有的模拟正交盒子中，不做处理。
步骤202，对扩展盒子内的各原子分别建立近邻集合；
具体为：
首先，利用扩展盒子内的各原子坐标分别除以截断半径并向下取整，作为该原子所在像素区域号；
然后，针对某个原子b，判断相邻像素区域内的原子c与原子b的距离是否小于原子b的截断半径，如果是，则将原子c添加到原子b的近邻集合中；否则，对原子c不做处理；依次对原子b的相邻像素区域内的原子进行判断，得到的原子b的近邻集合；
同理，得到扩展盒子内其余每个原子的近邻集合。
步骤203，针对扩展盒子内的第i个原子，将该原子与其近邻集合内各原子连线的垂直平分线构成闭包，通过voronoi方法计算闭包体积作为原子i的原子体积Vi。
同理，得到扩展盒子内每个原子的原子体积。
步骤204，通过扩展盒子内每个原子的体积和原子应力张量，计算特征区域的统计应力：
首先，计算每个原子的原子应力：
以坐标轴方向为参考，采取6个应力张量来表示，分别为：Sxx、Syy、Szz、Sxy、Syz、Sxz。
对于每一个应力张量，表达为：
Figure BDA0002867487250000111
其中，a、b分别代表x、y、z分量中的任意一个，m代表原子质量，v代表原子此刻的速度；Np代表原子的近邻数，ria表示构成原子成键关系的第i个原子由于成键导致的原子力在a方向上的分量差，Fib表示构成原子成键关系的第i个原子由于成键导致的原子力在b方向上的分量，i＝1,2,3,4；Nb代表原子所参与形成的键的数目，Na代表原子所参与形成的键角的个数，Nd代表原子所参与形成的二面角的个数，Ni代表第i个原子所参与形成的非正常二面角的个数，Nf代表施加在原子上的内约束力作用数；Kspace代表长程库伦相互作用。
这些基本量，在存在近邻集合的条件下，通过匹配势函数，可以直接从系统中获取。
然后，对于选取的模拟双金属界面模型的特征区域的上下层，分别利用各原子的体积计算特征区域上下层的统计应力；
即：
Figure BDA0002867487250000121
Sab表示原子应力张量，a、b分别代表x、y、z分量中的任意一个，σ表示特征区域应力，N代表待测特征区域内的所有原子数量。
步骤三，分别判断特征区域内上下层的应力是否达到应力平衡，如果是，则执行步骤五；否则，进入步骤四；
步骤四，针对不同加载模式，对双金属界面模型的可变区域施加局部应变，直至特征区域内的应力达到平衡，进入步骤五。
加载模式包括：初始加载模式、平行界面单轴加载模式、平行界面双轴加载模式和界面法向加载模式。
(1)针对初始加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡的具体过程如下：
步骤1011，计算初始加载模式的可变区域界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000122
Figure BDA0002867487250000123
Figure BDA0002867487250000124
其中，
Figure BDA0002867487250000125
表示特征A2区统计应力xx方向上的分量，
Figure BDA0002867487250000126
表示特征A2区统计应力yy方向上的分量，
Figure BDA0002867487250000127
表示特征A2区统计应力zz方向上的分量，
Figure BDA0002867487250000128
表示特征B2区统计应力xx方向上的分量，
Figure BDA0002867487250000129
表示特征B2区统计应力yy方向上的分量，
Figure BDA00028674872500001210
表示特征B2区统计应力zz方向上的分量。
步骤1012，利用目标应力状态计算特征区域的应力平衡判据；
计算公式为：
Figure BDA00028674872500001211
步骤1013，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1014、判断当前次迭代的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1015；否则，特征区域内的应力达到平衡，结束算法；
步骤1015、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1014重复判断；
具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000131
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000132
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000133
其中λ为系数常数。E为常数，通常为材料平均弹性模量。
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1016、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(2)针对平行界面单轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡的具体过程如下：
步骤1021，以平行于x轴加载的可变区域为例，计算其界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000141
Figure BDA0002867487250000142
步骤1022，利用目标应力状态计算特征区域的应力平衡判据；
计算公式为：
Figure BDA0002867487250000143
步骤1023，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1024、判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1025；否则，特征区域内的应力达到平衡，结束算法；
步骤1025、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1024重复判断；
对平行于x轴记载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000144
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000145
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000151
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1026、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(3)针对平行界面双轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡的具体过程如下：
步骤1031，以平行于xy轴加载模式的可变区域为例，计算界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure BDA0002867487250000152
步骤1032，利用目标应力状态计算特征区域的应力平衡判据；
Figure BDA0002867487250000153
步骤1033，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1034，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1035；否则，特征区域内的应力达到平衡，结束算法；
步骤1035、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1034重复判断；
对平行于xy轴加载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000154
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000161
对盒子无需修正。
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1036、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
(4)针对界面法向加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡的具体过程如下：
步骤1041，以平行于x轴界面法向加载模式的可变区域为例，其界面两侧的特征区域A2和B2的目标应力状态为：
Figure BDA0002867487250000162
Figure BDA0002867487250000163
Figure BDA0002867487250000164
步骤1042，利用目标应力状态计算特征区域的应力平衡判据：
Figure BDA0002867487250000165
步骤1043，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1044，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1045；否则，特征区域内的应力达到平衡，结束算法；
步骤1045、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1044重复判断；
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，对可变A1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000171
对可变B1区内原子坐标乘以应变矩阵：
Figure BDA0002867487250000172
对盒子对角线坐标乘以应变矩阵：
Figure BDA0002867487250000173
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1046、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法。
步骤五，输出双金属界面模型达到平衡的特征区域应力判据，通过共轭梯度算法(CG)以及快速松弛引擎算法(FIRE)对计算最小化能量，并更新原子坐标，获得双金属界面在内应力平衡条件下的加载模拟结果。
具体为：
步骤501，步骤四输出的平衡特征区域应力判据D0，通过最小化能量更新了原子坐标后，对应的各原子的体积和应力张量也发生变化，得到更新后的统计应力和应力平衡判据D1；
初始设定迭代次数m＝1，
步骤502，对应力平衡判据D1进行判断，是否满足D1＜Dmax，如果是，应力平衡判据D1为双金属界面在内应力平衡条件下的加载模拟结果。否则，进入步骤503；
步骤503，迭代次数m自增1，利用步骤四对应力平衡判据D1加载不同模式，得到新的应力平衡判据D1'；
步骤504，新的应力平衡判据D1'通过最小化能量更新原子坐标，重复更新应力平衡判据并返回步骤502进行判断；
步骤505，针对m达到最大迭代次数时，仍没有应力平衡判据满足Dm＜Dmax，输出最大m对应的应力平衡判据作为最终结果输出，结束算法。

Claims (5)

1.一种基于选区应力判据的双金属界面自适应加载模拟方法，其特征在于，包括如下步骤：
步骤一，通过计算机模拟双金属界面模型，对双金属界面模型进行区域划分，具体划分为可变区域、特征区域与固定区域；
具体为：
首先，计算机预制双金属界面模型，以及与界面模型对应的势函数；
然后，利用模型界面和势函数对双金属界面模型进行区域划分，得到可变区域、特征区域与固定区域：
可变区域分上下层，上层称为可变A1区，下层称为可变B1区，可变A1区为模型界面z>0的区域，可变B1区为z<0的区域；
特征区域分上下层，上层称为特征A2区，下层称为特征B2区，特征A2区为模型界面z>h1∩z<h2的区域，特征B2区为模型界面z>-h2∩z<-h1的区域；其中h1和h2分别为用户指定参数；
固定区域为距离盒子z方向边界h3以内的区域；h3为用户指定参数，且h3大于势函数的截断半径；
步骤二，通过周期性边界条件得到扩展盒子，建立扩展盒子内各原子的近邻集合，并计算各原子的体积，最后利用各原子的体积统计特征区域内上下层的应力：
步骤三，分别判断特征区域内上下层的应力是否达到应力平衡，如果是，进入步骤五；否则，进入步骤四；
步骤四，利用不同加载模式对可变区域施加局部应变，直至特征区域内的应力达到平衡，进入步骤五；
加载模式包括：初始加载模式、平行界面单轴加载模式、平行界面双轴加载模式和界面法向加载模式；
(1)对初始加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1011、计算初始加载模式下可变区域界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure FDA0003641201820000011
Figure FDA0003641201820000012
Figure FDA0003641201820000013
其中，
Figure FDA0003641201820000014
表示特征A2区统计应力xx方向上的分量，
Figure FDA0003641201820000015
表示特征A2区统计应力yy方向上的分量，
Figure FDA0003641201820000016
表示特征A2区统计应力zz方向上的分量，
Figure FDA0003641201820000017
表示特征B2区统计应力xx方向上的分量，
Figure FDA0003641201820000021
表示特征B2区统计应力yy方向上的分量，
Figure FDA0003641201820000022
表示特征B2区统计应力zz方向上的分量；
步骤1012、利用特征A2区和特征B2区的目标应力状态计算特征区域的应力平衡判据；
计算公式为：
Figure FDA0003641201820000023
步骤1013、用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1014、判断当前次迭代的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1015；否则，特征区域内的应力达到平衡，结束算法；
步骤1015、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0，进入步骤1014重复判断；
具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000024
对可变B1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000025
对盒子对角线坐标乘以应变矩阵：
Figure FDA0003641201820000026
其中λ为系数常数，E为材料平均弹性模量的常数；
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1016、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法；
(2)对平行界面单轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1021、以平行于x轴加载的可变区域为例，计算其两侧的特征A2区和特征B2区的目标应力状态为：
Figure FDA0003641201820000031
Figure FDA0003641201820000032
步骤1022、利用目标应力状态计算特征区域的应力平衡判据；
计算公式为：
Figure FDA0003641201820000033
步骤1023、用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1024、判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1025；否则，特征区域内的应力达到平衡，结束算法；
步骤1025、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1024重复判断；
对平行于x轴记载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000034
对可变B1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000041
对盒子对角线坐标乘以应变矩阵：
Figure FDA0003641201820000042
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1026、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法；
(3)对平行界面双轴加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1031，以平行于xy轴加载模式的可变区域为例，计算界面两侧的特征A2区和特征B2区的目标应力状态为：
Figure FDA0003641201820000043
步骤1032，利用目标应力状态计算特征区域的应力平衡判据；
Figure FDA0003641201820000044
步骤1033，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1034，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1035；否则，特征区域内的应力达到平衡，结束算法；
步骤1035、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1034重复判断；
对平行于xy轴加载模式的可变区域施加局部应变，具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，具体为：对可变A1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000051
对可变B1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000052
对盒子无需修正；
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk，f＜Dk,0，结束迭代过程；
步骤1036、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法；
(4)对界面法向加载模式的可变区域施加局部应变，使特征区域内的应力达到平衡；
具体过程如下：
步骤1041，以平行于x轴界面法向加载模式的可变区域为例，其界面两侧的特征区域A2和B2的目标应力状态为：
Figure FDA0003641201820000053
Figure FDA0003641201820000054
Figure FDA0003641201820000055
步骤1042，利用目标应力状态计算特征区域的应力平衡判据：
Figure FDA0003641201820000056
步骤1043，用户设定迭代精度Dmax，定义迭代失败次数为f，f赋值为0，定义迭代次数为k，初始k＝0；
步骤1044，判断当前次的应力平衡判据是否满足Dk，0＞Dmax，如果是，进入步骤1045；否则，特征区域内的应力达到平衡，结束算法；
步骤1045、通过不断增加失败次数f的值，改变施加的局部应变矩阵，得到不同的应力平衡判据，直至某次应力平衡判据满足Dk,f＜Dk,0，将失败次数f置0，迭代次数k自增1，输出Dk,f作为应力平衡判据Dk+1，0进入步骤1044重复判断；
具体为：
首先，当前特征区域判据Dk，0未达到应力平衡，对可变区域施加f＝0的局部应变，对可变A1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000061
对可变B1区内原子坐标乘以应变矩阵：
Figure FDA0003641201820000062
对盒子对角线坐标乘以应变矩阵：
Figure FDA0003641201820000063
然后、对可变区域施加f＝0的局部应变后，特征区域内的原子坐标发生改变，导致各原子的应力张量和体积对应改变，重新计算应力平衡判据Dk，1；
判断是否满足Dk，1＞Dk，0；如果是，则f自增1，并撤销已施加的原f＝0对应的应变，继续对可变区域施加f＝1的局部应变，并判断得到的应力平衡判据是否满足Dk，2＞Dk，0；如果是，f继续自增1，撤销上次施加的应变，重复迭代直至Dk,f＜Dk,0，结束迭代过程；
步骤1046、当迭代次数k达到用户预设的总迭代次数时，仍不满足Dk，0＜Dmax，输出最大k对应的应力平衡判据作为最终结果输出，结束算法；
步骤五，输出当前特征区域应力平衡判据，通过共轭梯度算法以及快速松弛引擎算法对达到应力平衡的双金属界面模型计算最小化能量，并更新原子坐标，获得双金属界面在内应力平衡条件下的加载模拟结果。
2.如权利要求1所述的一种基于选区应力判据的双金属界面自适应加载模拟方法，其特征在于，步骤一中所述的双金属界面模型包括盒子信息、盒子内原子类型和所有原子的坐标；
盒子信息包括模拟正交盒子的对角线两个端点坐标以及盒子内原子总数；
盒子内第i个原子坐标为(xi,yi,zi)，每个原子坐标均包含界面法方向，界面法方向平行于直角坐标系中的z轴，且模型界面在z＝0面上，界面层厚度大于
Figure FDA0003641201820000071
其中，z轴能替换为x轴或y轴；
所述的势函数为分子动力学或静力学可用势函数，不同的势函数对应的截断半径不同。
3.如权利要求1所述的一种基于选区应力判据的双金属界面自适应加载模拟方法，其特征在于，所述的步骤二具体是：
步骤201，针对模拟正交盒子内的原子，应用周期性边界条件选择与盒子边界的距离小于截断半径的原子复制到扩展盒子中；
首先，对模拟正交盒子外的原子重整，将位于周期性边界外的原子移入模拟正交盒子中；
然后，将原有的模拟正交盒子沿边界向外扩展，形成包括原正交盒子在内的新的扩展盒子；
接着，对位于原正交盒子内的原子逐一进行判断：将某个原子a在某一方向平移一个周期性边界对应的周期后，判断平移后的原子a是否到了扩展部分中；如果是，则将原子a复制到扩展部分中，并继承原子a平移之前的原子类型和原子坐标；否则，原子a仍在原有的模拟正交盒子中，不做处理；
步骤202，对扩展盒子内的各原子分别建立近邻集合；
具体为：
首先，利用各原子坐标分别除以截断半径并向下取整，作为该原子所在像素区域号；
然后，针对某个原子b，判断相邻像素区域内的原子c与原子b的距离是否小于原子b的截断半径，如果是，则将原子c添加到原子b的近邻集合中；否则，对原子c不做处理；依次对原子b的相邻像素区域内的原子进行判断，得到的原子b的近邻集合；
同理，得到扩展盒子内其余每个原子的近邻集合；
步骤203，针对扩展盒子内的第i个原子，将该原子与其近邻集合内各原子连线的垂直平分线构成闭包，通过voronoi方法计算闭包体积作为原子i的原子体积Vi；
同理，得到扩展盒子内每个原子的原子体积；
步骤204，通过扩展盒子内每个原子的体积和应力张量，计算特征区域的统计应力：
首先，计算每个原子的原子应力：
以坐标轴方向为参考，采取6个应力张量来表示，分别为：Sxx、Syy、Szz、Sxy、Syz、Sxz；对于每一个应力张量，表达为：
Figure FDA0003641201820000081
其中，a、b分别代表x、y、z分量中的任意一个，m代表原子质量，v代表原子此刻的速度；Np代表原子的近邻数，ria表示第i个原子由于成键导致的原子力在a方向上的分量差，Fib表示第i个原子由于成键导致的原子力在b方向上的分量，i＝1,2,3,4；Nb代表原子所参与形成的键的数目，Na代表原子所参与形成的键角的个数，Nd代表原子所参与形成的二面角的个数，Ni代表第i个原子所参与形成的非正常二面角的个数，Nf代表施加在原子上的内约束力作用数；Kspace代表长程库伦相互作用；
然后，对于选取的模拟双金属界面模型的特征区域的上下层，分别利用各原子的体积计算特征区域上下层的统计应力；
即：
Figure FDA0003641201820000082
σ表示特征区域应力，N代表扩展盒子内的所有原子数量。
4.如权利要求3所述的一种基于选区应力判据的双金属界面自适应加载模拟方法，其特征在于，所述的步骤201中，扩展长度为原有的模拟正交盒子的势函数对应的截断半径的长度。
5.如权利要求1所述的一种基于选区应力判据的双金属界面自适应加载模拟方法，其特征在于，所述步骤五具体为：
步骤501，当前输出的平衡特征区域应力判据D0，通过最小化能量更新了原子坐标后，对应的各原子的体积和应力张量也发生变化，得到更新后的统计应力和应力平衡判据D1；
初始设定迭代次数m＝1，
步骤502，对应力平衡判据D1进行判断，是否满足D1＜Dmax，如果是，应力平衡判据D1为双金属界面在内应力平衡条件下的加载模拟结果；否则，进入步骤503；
步骤503，迭代次数m自增1，对应力平衡判据D1加载不同模式，得到新的应力平衡判据D1'；
步骤504，新的应力平衡判据D1'通过最小化能量更新原子坐标，重复更新应力平衡判据并进入步骤502进行判断；
步骤505，针对m达到最大迭代次数时，仍没有应力平衡判据满足Dm＜Dmax，输出最大m对应的应力平衡判据作为最终结果输出，结束算法。
CN202011587163.6A 2019-12-30 2020-12-29 一种基于选区应力判据的双金属界面自适应加载模拟方法 Active CN112580223B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2019113915603 2019-12-30
CN201911391560 2019-12-30

Publications (2)

Publication Number Publication Date
CN112580223A CN112580223A (zh) 2021-03-30
CN112580223B true CN112580223B (zh) 2022-08-02

Family

ID=75143801

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011587163.6A Active CN112580223B (zh) 2019-12-30 2020-12-29 一种基于选区应力判据的双金属界面自适应加载模拟方法

Country Status (1)

Country Link
CN (1) CN112580223B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113421616B (zh) * 2021-06-23 2022-07-19 北京航空航天大学 一种基于自由体积分布的溶质低能偏聚界面模型建模方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105598178A (zh) * 2015-12-28 2016-05-25 北京科技大学 基于数值模拟的复合板界面结合强度的工艺参数控制方法
CN105930619A (zh) * 2016-05-17 2016-09-07 上海交通大学 纤维增强复合材料物理非线性模拟的态型近场动力学方法
CN110614275A (zh) * 2019-11-12 2019-12-27 太原理工大学 一种强变形轧制双金属复合板的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9524555B2 (en) * 2011-12-12 2016-12-20 Beihang University Method and computer program product of the simultaneous pose and points-correspondences determination from a planar model

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105598178A (zh) * 2015-12-28 2016-05-25 北京科技大学 基于数值模拟的复合板界面结合强度的工艺参数控制方法
CN105930619A (zh) * 2016-05-17 2016-09-07 上海交通大学 纤维增强复合材料物理非线性模拟的态型近场动力学方法
CN110614275A (zh) * 2019-11-12 2019-12-27 太原理工大学 一种强变形轧制双金属复合板的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
内爆加载金属界面不稳定性的数值分析;郝鹏程;《爆炸与冲击》;北京应用物理与计算数学研究所;20161130;第36卷(第6期);第740页-第743页 *

Also Published As

Publication number Publication date
CN112580223A (zh) 2021-03-30

Similar Documents

Publication Publication Date Title
Keshavarz et al. Hierarchical crystal plasticity FE model for nickel-based superalloys: Sub-grain microstructures to polycrystalline aggregates
CN113051831B (zh) 机床热误差自学习预测模型建模方法及热误差控制方法
Curreli et al. Application of the finite element submodeling technique in a single point contact and wear problem
CN110515301B (zh) 一种结合dmpc的改进的admm算法
CN112580223B (zh) 一种基于选区应力判据的双金属界面自适应加载模拟方法
CN104539601B (zh) 动态网络攻击过程可靠性分析方法及系统
CN112800675A (zh) 一种基于kpca和elm的时空分离分布参数系统建模方法
CN108762072B (zh) 基于核范数子空间法和增广向量法的预测控制方法
CN115997214A (zh) 计算机辅助设计和制造的载荷循环内防止损伤的生成设计形状优化
CN116018595A (zh) 计算机辅助设计和制造的使用构建材料强度模型的生成设计形状优化
WO2009067952A1 (fr) Procédé de commande de technique et son dispositif
Bonnet et al. Evaluation of numerical time‐integration schemes for real‐time hybrid testing
CN116306248A (zh) 基于时空非线性误差补偿模型的锂电池温度场预测方法
CN106919759B (zh) 基于拟合灵敏度的航空发动机性能的建模方法及模型应用
CN106067075B (zh) 一种建筑用能负荷预测模型建立、负荷预测方法及其装置
CN111625995B (zh) 一种集成遗忘机制和双超限学习机的在线时空建模方法
Nair et al. Crack tip enhanced phase-field model for crack evolution in crystalline Ti6Al from concurrent crystal plasticity FE-molecular dynamics simulations
CN113159318A (zh) 一种神经网络的量化方法、装置、电子设备及存储介质
CN111103789B (zh) 源网荷综合能源调度分析方法、系统及终端设备
Pedro et al. NARMA-L2 control of a nonlinear half-car servo-hydraulic vehicle suspension system
CN114925584B (zh) 变向锻造成形的金属晶粒尺寸和取向的智能协同调控方法
CN113032892B (zh) 一种基于稳定性算法的蒙皮壁板参数优化方法
CN116432411B (zh) 一种纯金属高压原子间相互作用势函数的构建方法
WO2023143110A1 (zh) 一种偏微分方程求解方法及其相关设备
CN117648533A (zh) 一种基于迁移学习的薄壁件加工变形预测方法

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