CN115906583A - 一种基于虚单元法的药柱结构完整性仿真分析方法及系统 - Google Patents

Abstract

本发明公开了一种基于虚单元法的药柱结构完整性仿真分析方法及系统，其中方法包括：获取多边形网格数据、药柱材料数据和载荷数据；基于多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；基于多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；根据单元内力计算单元残余力，当单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；计算时间满足设定要求后输出节点位移、单元应力以及应变信息。本发明提出的基于虚单元法的药柱结构完整性仿真分析方法能够应对复杂的药柱结构并开展精确且快速的药柱结构完整性分析计算。

Description

一种基于虚单元法的药柱结构完整性仿真分析方法及系统

技术领域

本发明属于药柱结构完整性分析技术领域，具体涉及一种基于虚单元法的药柱结构完整性仿真分析方法及系统。

背景技术

药柱是具有一定几何形状和尺寸的固体推进剂。安放于固体火箭发动机燃烧室中的药柱，几何形状和尺寸的选择与发动机的工作时间、燃烧室压力和推力有关，同时也影响药柱的结构完整性。为了满足总体指标，药柱的构型往往会比较复杂，这就对仅仅支持三角形、四边形、四面体以及六面体等固定单元形状的有限元方法带来极大的挑战。因此，在利用有限元进行药柱结构完整性分析过程中，几何简化以及网格划分处理几乎占据了绝大部分时间。另一方面，在运输以及战备服役阶段，药柱内部会发生裂纹等缺陷，而有限元方法在处理裂纹扩展问题时面临的最大挑战就是裂纹扩展后单元的重构问题。因此，有必要在保证计算精度的前提条件下，寻求一种新的快捷高效的药柱结构完整性分析手段。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种基于虚单元法的药柱结构完整性仿真分析方法及系统，能够应对复杂的药柱结构并开展精确且快速的药柱结构完整性分析计算。

本发明提供了如下的技术方案：

第一方面，提供一种基于虚单元法的药柱结构完整性仿真分析方法，包括：

获取多边形网格数据、药柱材料数据和载荷数据；

基于所述多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；

基于所述多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；

根据所述单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；

计算时间满足设定要求后输出节点位移、单元应力以及应变信息。

进一步的，所述多边形网格数据包括所有节点的坐标参数和所有单元的节点组成，所述药柱材料数据包括松弛模量和泊松比参数，所述载荷数据包括温度载荷参数和压力载荷参数。

进一步的，所述多边形单元刚度的计算方法包括：

按照逆时针方向标记多边形单元顶点分别为V_i，标记多边形单元边分别为e_i，i＝1,...,N^V，则单元的形函数矩阵W_c的表达式为：

其中，

分别表示单元的形函数矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_c)_a的表达式为：

其中，q_1a表示单元的形函数矩阵第a个分量的第1分量，q_2a表示单元的形函数矩阵第a个分量的第2分量，q_1a和q_2a的表达式为：

其中，|e_a-1|、|e_a|分别表示单元边e_a-1、e_a的长度，|E|表示单元面积，{n₁}_a-1、{n₁}_a分别表示单元边e_a-1、e_a外法线方向的第1个分量，{n₂}_a-1、{n₂}_a分别表示单元边e_a-1、e_a外法线方向的第2个分量；

药柱的松弛模量E(t)的表达式为：

其中，E_n和

分别表示第n项松弛模量的两个参数，N_E表示松弛模量Prony级数的项数，E₀表示初始松弛模量，t表示加载时间；

单元的材料矩阵D的表达式为：

其中，

其中，ν表示药柱材料的泊松比，γ^E()、

分别表示应力张量辅助第一变量、应力张量辅助第二变量、应力张量辅助第三变量，Δt表示当前步的时间增量；

单元刚度矩阵连续项K_c的表达式为：

单元刚度矩阵稳定项K_s的表达式为：

其中，

表示2N^V×2N^V的单位矩阵，P表示刚度矩阵的稳定项辅助第一矩阵，其表达式为：

其中，H_R表示刚度矩阵的稳定项辅助第二矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第二矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_R)_a的表达式为：

其中，x_1a和x_2a分别表示顶点V_a的横纵坐标，

和

表示多边形的中心点的横纵坐标；

H_C表示刚度矩阵的稳定项辅助第三矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第三矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_C)_a的表达式为：

W_R表示刚度矩阵的稳定项辅助第四矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第四矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_R)_a的表达式为：

S_E表示刚度矩阵的稳定项辅助第五矩阵，表达式为：

其中，α_E表示刚度矩阵的稳定项辅助第五矩阵的系数，其表达式为：

则多边形单元刚度K的表达式为：

进一步的，所述单元应力与应力增量的计算方法包括：

所述应力增量包括偏应力增量和球应力增量，t_m+1时刻的偏应力增量ΔS_ij(t_m+1)的表达式为：

其中，γ^E(Δt_m+1)表示应力张量辅助第一变量对应Δt_m+1时的函数值，Δt_m+1表示t_m+1时刻的时间增量，Δe_ij(t_m+1)表示t_m+1时刻的偏应变增量，

表示单元的第一状态变量，其满足递推关系：

其中，

表示应力张量辅助第二变量对应Δt_m时的函数值，Δt_m表示t_m时刻的时间增量，Δe_ij(t_m)表示t_m时刻的偏应变增量，

表示单元的第一状态变量对应t_m时的函数值；

t_m+1时刻的球应力增量Δσ_kk(t_m+1)的表达式为：

其中，Δe_kk(t_m+1)表示t_m+1时刻的球应变增量，

表示单元的第二状态变量，其满足递推关系：

其中，Δe_kk(t_m)表示t_m时刻的球应变增量，

表示单元的第二状态变量对应t_m时的函数值；

单元t_m+1时刻的应力增量Δσ_ij(t_m+1)的表达式为：

其中，δ_ij表示克罗内克尔符号；

单元t_m+1时刻的应力表达式为：

其中，σ_ij(t_m+1)为单元t_m+1时刻的应力，σ_ij(t_m)为单元t_m时刻的应力。

进一步的，所述单元内力的计算方法包括：

所述单元内力包括单元内力连续项和单元内力稳定项，所述单元内力连续项T_c的表达式为：

T_c＝K_cu (42)

其中，u为单元的位移矩阵；

所述单元内力稳定项T_s的表达式为：

其中，σ为单元的应力矩阵；

t_m+1时刻单元内力T总的表达式为：

T＝T_c+T_s (44)。

进一步的，所述单元残余力的计算方法包括：

当前增量步内的残余力R_j的表达式为：

R_j＝F_ct-T_j                     (45)

其中，F_ct为当前载荷步结构上的外力，T_j为当前载荷步第j个增量步上的单元总的内力。

进一步的，所述节点位移的计算方法包括：

当前增量步上节点位移u_j的表达式为：

u_j＝u_j-1+Δu_j (46)

其中，u_j-1为上一增量步上的节点位移，Δu_j表示当前增量步上的节点位移增量，其表达式为：

其中，K_j为当前载荷步第j个增量步上的单元总的刚度。

第二方面，提供一种基于虚单元法的药柱结构完整性仿真分析系统，包括：

数据获取模块，用于获取多边形网格数据、药柱材料数据和载荷数据；

第一数据处理模块，用于基于所述多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；

第二数据处理模块，用于基于所述多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；

第三数据处理模块，用于根据所述单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；

数据输出模块，用于计算时间满足设定要求后输出节点位移、单元应力以及应变信息。

第三方面，提供一种基于虚单元法的药柱结构完整性仿真分析装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行第一方面所述方法的步骤。

第四方面，提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现第一方面所述方法的步骤。

与现有技术相比，本发明的有益效果是：

本发明基于多边形网格数据、药柱材料数据和载荷数据计算多边形单元刚度并更新单元应力与应力增量，然后计算单元内力，计算单元残余力并迭代至单元残余力达到设定标准后开展下一载荷步的计算，当计算时间满足设定要求后输出节点位移、单元应力以及应变信息；相比于传统的有限元分析方法，本发明提出的基于虚单元法的药柱结构完整性仿真分析方法，能够应对十分复杂的药柱结构并开展精确且快速的药柱结构完整性分析计算。

附图说明

图1是本发明实施例中基于虚单元法的药柱结构完整性仿真分析的流程图；

图2是本发明实施例中典型多边形E的标记示意图；

图3是本发明实施例中基于虚单元法计算出的药柱结构在内压载荷下的径向应变云图；

图4是本发明实施例中基于虚单元法计算出的药柱结构在内压载荷下的环向应变云图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例1

如图1所示，本实施例提供一种基于虚单元法的药柱结构完整性仿真分析方法，步骤如下：

步骤1、获取多边形网格数据，多边形网格数据包括所有节点的坐标参数和所有单元的节点组成。

步骤2、获取药柱材料数据和载荷数据，其中，药柱材料数据包括松弛模量和泊松比参数，载荷数据包括温度载荷参数和压力载荷参数。

步骤3、基于步骤1和步骤2中的数据计算多边形单元刚度，具体计算方法包括：

按照逆时针方向标记多边形单元顶点分别为V_i，标记多边形单元边分别为e_i，i＝1,...,N^V，图2所示为典型多边形E的标记示意图，则单元的形函数矩阵W_c的表达式为：

其中，

分别表示单元的形函数矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_c)_a的表达式为：

其中，q_1a表示单元的形函数矩阵第a个分量的第1分量，q_2a表示单元的形函数矩阵第a个分量的第2分量，q_1a和q_2a的表达式为：

其中，|e_a-1|、|e_a|分别表示单元边e_a-1、e_a的长度，|E|表示单元面积，{n₁}_a-1、{n₁}_a分别表示单元边e_a-1、e_a外法线方向的第1个分量，{n₂}_a-1、{n₂}_a分别表示单元边e_a-1、e_a外法线方向的第2个分量；

药柱的松弛模量E(t)的表达式为：

其中，E_n和

分别表示第n项松弛模量的两个参数，N_E表示松弛模量Prony级数的项数，E₀表示初始松弛模量，t表示加载时间；

单元的材料矩阵D的表达式为：

其中，

其中，ν表示药柱材料的泊松比，γ^E()、

分别表示应力张量辅助第一变量、应力张量辅助第二变量、应力张量辅助第三变量，Δt表示当前步的时间增量；

单元刚度矩阵连续项K_c的表达式为：

K_c＝(W_c)^TD(W_c) (47)

单元刚度矩阵稳定项K_s的表达式为：

其中，

表示2N^V×2N^V的单位矩阵，P表示刚度矩阵的稳定项辅助第一矩阵，其表达式为：

其中，H_R表示刚度矩阵的稳定项辅助第二矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第二矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_R)_a的表达式为：

其中，x_1a和x_2a分别表示顶点V_a的横纵坐标，

和

表示多边形的中心点的横纵坐标；

H_C表示刚度矩阵的稳定项辅助第三矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第三矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_C)_a的表达式为：

W_R表示刚度矩阵的稳定项辅助第四矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第四矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_R)_a的表达式为：

S_E表示刚度矩阵的稳定项辅助第五矩阵，表达式为：

其中，α_E表示刚度矩阵的稳定项辅助第五矩阵的系数，其表达式为：

则多边形单元刚度K的表达式为：

K＝K_c+K_s                  (58)。

步骤4、基于步骤1和步骤2中的数据更新单元应力与应力增量，单元应力与应力增量的具体计算方法包括：

应力增量包括偏应力增量和球应力增量，t_m+1时刻的偏应力增量ΔS_ij(t_m+1)的表达式为：

其中，γ^E(Δt_m+1)表示应力张量辅助第一变量对应Δt_m+1时的函数值，Δt_m+1表示t_m+1时刻的时间增量，Δe_ij(t_m+1)表示t_m+1时刻的偏应变增量，

表示单元的第一状态变量，其满足递推关系：

其中，

表示应力张量辅助第二变量对应Δt_m时的函数值，Δt_m表示t_m时刻的时间增量，Δe_ij(t_m)表示t_m时刻的偏应变增量，

表示单元的第一状态变量对应t_m时的函数值；

t_m+1时刻的球应力增量ΔS_kk(t_m+1)的表达式为：

其中，Δe_kk(t_m+1)表示t_m+1时刻的球应变增量，

表示单元的第二状态变量，其满足递推关系：

其中，Δe_kk(t_m)表示t_m时刻的球应变增量，

表示单元的第二状态变量对应t_m时的函数值；

单元t_m+1时刻的应力增量Δσ_ij(t_m+1)的表达式为：

其中，δ_ij表示克罗内克尔符号；

单元t_m+1时刻的应力表达式为：

σ_ij(t_m+1)＝Δσ_ij(t_m+1)+σ_ij(t_m) (64)

其中，σ_ij(t_m+1)为单元t_m+1时刻的应力，σ_ij(t_m)为单元t_m时刻的应力。

步骤5、基于步骤3得到的多边形单元刚度以及步骤4更新的单元应力与应力增量，计算单元内力，具体计算方法包括：

所述单元内力包括单元内力连续项和单元内力稳定项，所述单元内力连续项T_c的表达式为：

T_c＝K_cu (65)

其中，u为单元的位移矩阵；

所述单元内力稳定项T_s的表达式为：

其中，σ为单元的应力矩阵；

t_m+1时刻单元内力T总的表达式为：

T＝T_c+T_s (67)。

步骤6、根据步骤5得到的单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算。

单元残余力的具体计算方法包括：

当前增量步内的残余力R_j的表达式为：

R_j＝F_ct-T_j                     (68)

其中，F_ct为当前载荷步结构上的外力，T_j为当前载荷步第j个增量步上的单元总的内力。

节点位移的具体计算方法包括：

当前增量步上节点位移u_j的表达式为：

u_j＝u_j-1+Δu_j (69)

其中，u_j-1为上一增量步上的节点位移，Δu_j表示当前增量步上的节点位移增量，其表达式为：

其中，K_j为当前载荷步第j个增量步上的单元总的刚度。

步骤7、计算时间满足设定要求后输出步骤6得到的节点位移、步骤4得到的单元应力以及应变信息，其中应变ε的表达式为：

其中，

表示位移的梯度。

实施例2

本实施例采用实施例1提供的基于虚单元法的药柱结构完整性仿真分析方法对圆管形式的药柱结构进行完整性分析。

图3为基于虚单元法计算出的药柱结构在内压载荷下的径向应变云图，图4为基于虚单元法计算出的药柱结构在内压载荷下的环向应变云图。

实施例3

本实施例提供一种基于虚单元法的药柱结构完整性仿真分析系统，包括：

数据获取模块，用于获取多边形网格数据、药柱材料数据和载荷数据；

第一数据处理模块，用于基于所述多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；

第二数据处理模块，用于基于所述多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；

第三数据处理模块，用于根据所述单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；

数据输出模块，用于计算时间满足设定要求后输出节点位移、单元应力以及应变信息。

实施例4

本实施例提供一种基于虚单元法的药柱结构完整性仿真分析装置，包括处理器及存储介质；所述存储介质用于存储指令；所述处理器用于根据所述指令进行操作以执行实施例1所述方法的步骤。

实施例5

本实施例提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现实施例1所述方法的步骤。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (10)

1.一种基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，包括：

获取多边形网格数据、药柱材料数据和载荷数据；

基于所述多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；

基于所述多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；

根据所述单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；

计算时间满足设定要求后输出节点位移、单元应力以及应变信息。

2.根据权利要求1所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述多边形网格数据包括所有节点的坐标参数和所有单元的节点组成，所述药柱材料数据包括松弛模量和泊松比参数，所述载荷数据包括温度载荷参数和压力载荷参数。

3.根据权利要求1所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述多边形单元刚度的计算方法包括：

按照逆时针方向标记多边形单元顶点分别为V_i，标记多边形单元边分别为e_i，i＝1,...,N^V，则单元的形函数矩阵W_c的表达式为：

其中，

分别表示单元的形函数矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_c)_a的表达式为：

其中，q_1a表示单元的形函数矩阵第a个分量的第1分量，q_2a表示单元的形函数矩阵第a个分量的第2分量，q_1a和q_2a的表达式为：

其中，|e_a-1|、|e_a|分别表示单元边e_a-1、e_a的长度，|E|表示单元面积，{n₁}_a-1、{n₁}_a分别表示单元边e_a-1、e_a外法线方向的第1个分量，{n₂}_a-1、{n₂}_a分别表示单元边e_a-1、e_a外法线方向的第2个分量；

药柱的松弛模量E(t)的表达式为：

其中，E_n和

分别表示第n项松弛模量的两个参数，N_E表示松弛模量Prony级数的项数，E₀表示初始松弛模量，t表示加载时间；

单元的材料矩阵D的表达式为：

其中，

其中，ν表示药柱材料的泊松比，γ^E()、

分别表示应力张量辅助第一变量、应力张量辅助第二变量、应力张量辅助第三变量，Δt表示当前步的时间增量；

单元刚度矩阵连续项K_c的表达式为：

K_c＝(W_c)^TD(W_c) (1)

单元刚度矩阵稳定项K_s的表达式为：

其中，

表示2N^V×2N^V的单位矩阵，P表示刚度矩阵的稳定项辅助第一矩阵，其表达式为：

其中，H_R表示刚度矩阵的稳定项辅助第二矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第二矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_R)_a的表达式为：

其中，x_1a和x_2a分别表示顶点V_a的横纵坐标，

和

表示多边形的中心点的横纵坐标；

H_C表示刚度矩阵的稳定项辅助第三矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第三矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(H_C)_a的表达式为：

W_R表示刚度矩阵的稳定项辅助第四矩阵，表达式为：

其中，

分别表示刚度矩阵的稳定项辅助第四矩阵的第1个、第a个以及第N^V个分量，顶点V_a对应的(W_R)_a的表达式为：

S_E表示刚度矩阵的稳定项辅助第五矩阵，表达式为：

其中，α_E表示刚度矩阵的稳定项辅助第五矩阵的系数，其表达式为：

则多边形单元刚度K的表达式为：

K＝K_c+K_s (12)。

4.根据权利要求3所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述单元应力与应力增量的计算方法包括：

所述应力增量包括偏应力增量和球应力增量，t_m+1时刻的偏应力增量ΔS_ij(t_m+1)的表达式为：

其中，γ^E(Δt_m+1)表示应力张量辅助第一变量对应Δt_m+1时的函数值，Δt_m+1表示t_m+1时刻的时间增量，Δe_ij(t_m+1)表示t_m+1时刻的偏应变增量，

表示单元的第一状态变量，其满足递推关系：

其中，

表示应力张量辅助第二变量对应Δt_m时的函数值，Δt_m表示t_m时刻的时间增量，Δe_ij(t_m)表示t_m时刻的偏应变增量，

表示单元的第一状态变量对应t_m时的函数值；

t_m+1时刻的球应力增量Δσ_kk(t_m+1)的表达式为：

其中，Δe_kk(t_m+1)表示t_m+1时刻的球应变增量，

表示单元的第二状态变量，其满足递推关系：

其中，Δe_kk(t_m)表示t_m时刻的球应变增量，

表示单元的第二状态变量对应t_m时的函数值；

单元t_m+1时刻的应力增量Δσ_ij(t_m+1)的表达式为：

其中，δ_ij表示克罗内克尔符号；

单元t_m+1时刻的应力表达式为：

σ_ij(t_m+1)＝Δσ_ij(t_m+1)+σ_ij(t_m) (18)

其中，σ_ij(t_m+1)为单元t_m+1时刻的应力，σ_ij(t_m)为单元t_m时刻的应力。

5.根据权利要求4所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述单元内力的计算方法包括：

所述单元内力包括单元内力连续项和单元内力稳定项，所述单元内力连续项T_c的表达式为：

T_c＝K_cu (19)

其中，u为单元的位移矩阵；

所述单元内力稳定项T_s的表达式为：

其中，σ为单元的应力矩阵；

t_m+1时刻单元内力T总的表达式为：

T＝T_c+T_s   (21)。

6.根据权利要求1所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述单元残余力的计算方法包括：

当前增量步内的残余力R_j的表达式为：

R_j＝F_ct-T_j                     (22)

其中，F_ct为当前载荷步结构上的外力，T_j为当前载荷步第j个增量步上的单元总的内力。

7.根据权利要求6所述的基于虚单元法的药柱结构完整性仿真分析方法，其特征在于，所述节点位移的计算方法包括：

当前增量步上节点位移u_j的表达式为：

u_j＝u_j-1+Δu_j (23)

其中，u_j-1为上一增量步上的节点位移，Δu_j表示当前增量步上的节点位移增量，其表达式为：

其中，K_j为当前载荷步第j个增量步上的单元总的刚度。

8.一种基于虚单元法的药柱结构完整性仿真分析系统，其特征在于，包括：

数据获取模块，用于获取多边形网格数据、药柱材料数据和载荷数据；

第一数据处理模块，用于基于所述多边形网格数据、药柱材料数据和载荷数据，计算多边形单元刚度，并更新单元应力与应力增量；

第二数据处理模块，用于基于所述多边形单元刚度以及更新的单元应力与应力增量，计算单元内力；

第三数据处理模块，用于根据所述单元内力计算单元残余力，当所述单元残余力大于设定标准时，更新节点位移并重新计算单元残余力，直至单元残余力达到设定标准后，开展下一载荷步的计算；

数据输出模块，用于计算时间满足设定要求后输出节点位移、单元应力以及应变信息。

9.一种基于虚单元法的药柱结构完整性仿真分析装置，其特征在于，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行权利要求1～7任一项所述方法的步骤。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1～7任一项所述方法的步骤。

Images