CN109740182A - 一种基于再生核粒子的无网格物理变形仿真方法 - Google Patents

Abstract

本发明公开了一种基于再生核粒子的无网格物理变形仿真方法，包括以下步骤：按照粒子密度对仿真几何模型进行粒子离散，得到离散颗粒；选用修正函数和窗函数构造形函数；根据影响域内离散颗粒的位移和形函数构造全场位移场函数；利用边界转化方法进行边界处理；联立动量方程和几何方程、物理方程、边界条件构造等效积分弱形式方程；将全场位移场函数带入等效积分弱形式方程求解，得到刚度矩阵、等效内力和外力矩阵；利用中心差分法进行代数求解，得到不同时间步的形变和应力。通过对比有限元方法，该方法能够快速有效地模拟大变形等工程问题，解决了传统方法不收敛或需要人工网格重画分计算的问题，有利于工程实际问题结构的设计和优化。

Description

一种基于再生核粒子的无网格物理变形仿真方法

技术领域

本发明涉及计算机辅助工程中的仿真计算领域，尤其涉及一种基于再生核粒子的无网格物理变形仿真方法。

背景技术

结构仿真是通过虚拟仿真对工程实际结构进行相关验证的过程，目前主要是通过有限元方法(FEM)实现。有限元方法在早期是以变分原理为基础发展起来的，广泛地应用于以拉普拉斯方程和泊松方程所描述的各类物理场中，其基本思想是通过将无限自由度的结构近似离散成有个自由度，通过构建自由度之间的关系函数，求解泛函的极值。

尽管有限元方法在结构仿真中已经取得了很多的成果，此项技术领域仍然具有很多有待发展的方面，比如大变形等易引起网格畸变的工程问题。由于有限元单元方法对网格的敏感性和依赖性，结构仿真结果很大程度依赖于结构网格处理的好坏。对于大变形等工程问题，由于结构发生了较大变形，网格畸变，会造成刚度矩阵奇异从而导致求解不收敛。为了解决此问题，需要网格重划技术，但网格重划会使单元数量增加，计算时间加长，计算精度降低，当不均匀变形继续增大，网格畸变严重甚至使得模拟无法继续进行。

橡胶材料是一种典型的超弹性材料，具有明显的大变形、大应变及高度非线性的力学特性，因此在对橡胶元件承载过程的模拟分析中，往往会因大变形导致橡胶单元网格出现严重扭曲，从而导致程序收敛失败使计算模拟过程无法进行。网格法不再适应于处理大变形问题是一直是困扰工程师们的问题。

近年来无网格法的思想被提了出来，并得到了迅速的发展。与有限元方法不同，无网格法是在对于一个影响域建立离散的系统方程时不用事先定义好网格的一种数值方法。与有限元方法相比，无网格法的优点包括：不需要网格；容易构造高阶形状函数；容易进行自适应分析。无网格法的这些优点使其具有相当的潜力解决上述提及的传统有限元方法所不能处理的极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题。

再生核粒子算法(RKPM,reproducing kernel particle method)作为先进的无网格法被引入到计算机辅助工程(Computer Aided Engineering，CAE)领域，其具有优于基于网格的有限元方法所不具备的优点。由于再生核粒子算法是一种无网格法，只需要粒子离散结构，克服了结构网格离散的不足，非常适合于大变形结构工程问题计算。与此同时，模型前处理简单，可以快速处理复杂结构模型。目前这一技术题目成为了当前前沿的方向，受到技术人员的广泛关注。因此，研究一种适合大变形问题的无网格物理变形仿真方法具有很重要的意义。

发明内容

本发明的目的是提供一种基于再生核粒子的无网格物理变形仿真方法，可以克服传统方法中对网格的依赖，适用于大变形问题的结构仿真。

为实现上述目的，本发明采用以下技术方案：

一种基于再生核粒子的无网格物理变形仿真方法，包括如下步骤：

(1)按照给定的粒子密度对仿真几何模型进行粒子离散，得到离散颗粒，确定所述离散颗粒的初始位置、影响域Ω和影响域半径ρ；

(2)基于再生核粒子算法，选用修正函数和窗函数的类型，构造所述影响域内所述离散颗粒的形函数；

(3)根据所述影响域Ω内的所述离散颗粒的位移和所述形函数构造全场位移场函数；

(4)利用边界转化方法进行边界处理，消除边界不一致性引起的误差；

(5)联立动量方程和几何方程、物理方程、边界条件构造等效积分弱形式方程；

(6)将所述全场位移场函数带入所述等效积分弱形式方程求解，得到刚度矩阵、等效内力和外力矩阵；

(7)利用中心差分法对等效积分弱形式方程进行代数求解，得到不同时间步的形变和应力。

优选地，步骤(1)中，所述离散颗粒的几何形状包括、但不限于圆形、球形、矩形、立方体。

优选地，步骤(1)中，所述影响域Ω由每个所述离散颗粒指定。

优选地，步骤(1)中，所述影响域半径ρ对于所有所述离散颗粒为恒定值或者对于每个单独的所述离散颗粒为不同的。

优选地，步骤(2)中，基于再生核粒子算法，所述形函数Φ_ρ(x；x-y)为：

Φ_ρ(x；x-y)＝C_ρ(x；x-y)k_ρ(x-y)

式中，C_ρ(x；x-y)为满足再生条件的修正函数，k_ρ(x-y)为满足再生条件的窗函数；

所述修正函数C_ρ(x；x-y)是由多项式基函数的线性组合表示而成，如下：

C_ρ(x；x-y)＝C₀(x)+C₁(x)(x-y)+C₂(x)(x-y)²+…+C_N(x)(x-y)^N＝H^T(x-y)C(x)

其中，N是满足再生条件的次数，H^T(x-y)是N次多项式基函数的向量；C(x) 为系数向量，通过施加N次再生条件来确定；C_i(x)为系数(i＝1,2,3……N)；

H^T(x-y)＝[1,x-y,(x-y)²,…,(x-y)^N]

C(x)＝[C₀(x),C₁(x),…,C_N(x)]

所述窗函数κ_ρ(x-y)是二维窗函数：

其中，为中间变量。

更优选地，所述窗函数选用三次样条函数，如下：

式中，

更优选地，所述窗函数包括、但不限于三次样条函数、B样条函数、高斯函数。

优选地，步骤(3)中，所述全场位移场函数u^R(x)为：

式中，u(y)表示所述影响域Ω内的所述离散颗粒的位移；

将影响域Ω用NP个所述离散粒子离散，获得所述全场位移场函数的离散近似式：

式中，NP表示所述影响域Ω内的所述离散粒子的总点数，u_I对应所述离散粒子I处的位移,ΔV_I为离散粒子的体积权重，N_I(x)为所述离散粒子I处的形函数，I＝1，2，……NP。

优选地，所述步骤(4)中，所述边界转化方法将所述全场位移场函数u^R(x) 分为两部分

式中，b表示位移边界上的离散颗粒，nb表示除去位移边界上的离散颗粒剩下的所有离散颗粒；位移边界上的离散颗粒总数为N_b，N_nb＝NP-N_b；N^b(x)为边界插值点上的形函数，N^nb(x)为非边界插值点的形函数，u^b为边界插值点位移， u^nb为非边界插值点位移。

通过强制位移边界条件消除误差，假定

u^b＝(D^b)^-1g

式中：g表示边界位移，D^b通过下式构造

优选地，所述步骤(5)中，

动量方程：

几何方程：ε_ij＝Lu_ij (2)

物理方程：σ_ij＝Dε_ij (3)

边界条件：

式中：σ_ij表示应力分量，ρ为密度，b_i为体积力，ε_ij表示应变分量，L为计算微分算子，u_ij为位移，D为弹性矩阵，n_j表示边界法向量，T_i为边界载荷力， S_σ为载荷边界域，S_u为位移边界域，ν_i为边界位移，为边界位移；

把式(2)、式(3)、式(4)带入式(1)，联立上述方程得到动量方程的等效积分弱形式方程：

上式中，δu_i,j为i方向位移对x_j坐标求偏导，V为求解体积域，s为边界域， A为面积域，u_ij为ij方向的位移，u_i为i方向位移，为i方向的加速度，为 S_σ边界域上的外力。

优选地，所述步骤(6)中，将重构的位移场函数u^R(x)带入等效积分弱形式方程，并将该方程写成矩阵形式得到下式：

式中，M为质量矩阵，为加速度场，f^int为等效内力矩阵，f^ext为等效外力矩阵。

优选地，所述步骤(7)中，所述中心差分法将时间平均分为n个时间间隔， 0，Δt，…t时刻的位移、速度、加速度已知，求解t+Δt时刻对应的解；

t时刻的加速度和速度利用中心差分法表示为

式中，u_t-Δt(x)为离散颗粒在t-Δt时刻的位移，u_t(x)为离散颗粒在t时刻的位移，u_t+Δt(x)为离散颗粒在t+Δt时刻的位移；

t+Δt时刻的位移解答可以根据时间t时刻的运动方程得到，即为下式

式中，M、C和K分别代表影响域的质量矩阵、阻尼矩阵和刚度矩阵；F_t是 t时刻外力和；

由此可得到中心差分法的递推公式

若t与t-Δt时刻的位移已知，则可以进一步求解出t+Δt时刻的位移；利用步骤(5)中的几何方程，将不同时刻位移带入计算，得到形变，即应变，再利用步骤(5)中的物理方程求解，将应变带入计算得到应力。

本发明提供了一种全新的算法应用手段，与现有技术相比，有益效果在于：本发明针对传统FEM方法对大变形工程问题难以求解的现状，基于再生核粒子算法，依据离散颗粒的初始位置，结合选用的修正函数和窗函数构造形函数，联合动量微分方程进行求解。通过对比表明，再生核粒子算法能够更加快速有效地模拟大变形等工程问题，解决了传统计算方法不收敛或者需要人工网格重画分计算的问题，极大地提高了计算的速度与结果的准确性，从而有利于工程实际问题结构的设计和优化。

附图说明

构成本申请的一部分附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1是本发明基于再生核粒子的无网格物理变形仿真方法的流程图。

具体实施方式

本发明提供一种基于再生核粒子的无网格物理变形仿真方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下参照附图并举实例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

图1是本发明基于再生核粒子的无网格物理变形仿真方法的流程图。

结合图1，本发明所述的一种基于再生核粒子的无网格物理变形仿真方法，包括如下步骤：

步骤1：按照给定的粒子密度对仿真几何模型进行粒子离散，得到离散颗粒，确定所述离散颗粒的初始位置、影响域Ω和影响域半径ρ。

其中，所述离散颗粒可以为圆形或球形，也可以为其他不同的几何形状，例如，矩形、立方体等。粒子离散后的所述几何模型由多个所述离散颗粒代表。影响域是指一个离散颗粒的影响区域，所述影响域由每个离散颗粒指定，这通过为每个离散颗粒指定影响域半径而完成。

每个离散颗粒具有它自身的影响域。所述影响域半径ρ用来控制所述离散颗粒的影响域的尺寸，对于所有离散颗粒为恒定值或者对于每个单独的所述离散颗粒为不同的。

步骤2：基于再生核粒子算法，选用修正函数和窗函数的类型，构造所述影响域内所述离散颗粒的形函数。

基于再生核粒子算法，所述形函数Φ_ρ(x,x-y)为：

Φ_ρ(x；x-y)＝C_ρ(x；x-y)k_ρ(x-y)

式中，C_ρ(x,x-y)为满足再生条件的修正函数或改进函数，它决定近似的完备性与一致性；k_ρ(x-y)为满足再生条件的窗函数，它是用来表征影响域内离散颗粒的光滑性和紧支性的函数，用于控制计算的连续性。

所述修正函数C_ρ(x,x-y)是由多项式基函数的线性组合表示而成，如下：

C_ρ(x；x-y)＝C₀(x)+C₁(x)(x-y)+C₂(x)(x-y)²+…+C_N(x)(x-y)^N＝H^T(x-y)C(x)

其中，N是满足再生条件的次数，H^T(x-y)是N次多项式基函数的向量；C(x) 为系数向量，通过施加N次再生条件来确定；C_i(x)为系数(i＝1,2,3……N)；

H^T(x-y)＝[1,x-y,(x-y)²,…,(x-y)^N]

C(x)＝[C₀(x),C₁(x),…,C_N(x)]

所述窗函数κ_ρ(x-y)是二维窗函数：

其中，为中间变量。

所述窗函数包括、但不限于三次样条函数、B样条函数、高斯函数。

由于样条型核函数具有紧支撑特性，因此在大多数数值计算中经常采用该类型的核函数以节省计算时间。本实例中，所述窗函数选用三次样条函数，具体如下：

式中，

步骤3：根据所述影响域Ω内的所述离散颗粒的位移和所述形函数构造全场位移场函数。

再生核粒子算法通过构造形函数建立位移场函数和离散颗粒处函数的关系，从而得到全场位移场函数，进一步进行求解。

所述全场位移场函数u^R(x)具体构造如下所示：

式中，u(y)表示所述影响域Ω内的所述离散颗粒的位移；

将影响域Ω用NP个所述离散粒子离散，获得所述全场位移场函数的离散近似式：

式中，NP表示所述影响域Ω内的所述离散粒子的总点数，u_I对应所述离散粒子I处的位移,ΔV_I为离散粒子的体积权重，N_I(x)为所述离散粒子I处的形函数，I＝1，2，……NP。

步骤4：利用边界转化方法进行边界处理，消除边界不一致性引起的误差。

无网格法均不能自然满足位移边界条件，因此在计算过程中必须采用适当的方法来强制满足边界条件，常用的方法有：拉格朗日乘子法、再生核粒子插值技术、修正变分原理法、坐标变换法等。这些方法大多数致力于构建具有Kronecker dalta特性的近似函数，各有其特点，适用于不同场合。例如拉格朗日乘子法其系数矩阵即不正定又不呈带状，所以运算量很大，但这一方法的精确度很高，适用于求解小规模二维问题；修正变原理法的系数矩阵虽然保持正定和带状性，但精度不高。但总的说来，和其他无网络方法一样，还没有一种简单而又完美的方法来处理再生核粒子的边界条件。因此，边界条件的处理是目前无网格方法研究的热点之一，直接关系到无网格方法的实用性。

本实施例中，采用一种对再生核粒子方法适用的一致方法。所述一致方法将所述全场位移场函数u^R(x)分为两部分

式中，b表示位移边界上的离散颗粒，nb表示除去位移边界上的离散颗粒剩下的所有离散颗粒；位移边界上的离散颗粒总数为N_b，N_nb＝NP-N_b；N^b(x)为边界插值点上的形函数，N^nb(x)为非边界插值点的形函数，u^b为边界插值点位移， u^nb为非边界插值点位移。

通过强制位移边界条件消除误差，假定

u^b＝(D^b)^-1g

式中：g表示边界位移，D^b通过下式构造

步骤5：联立动量方程和几何方程、物理方程、边界条件构造等效积分弱形式方程。

动量方程：

几何方程：ε_ij＝Lu_ij (2)

物理方程：σ_ij＝Dε_ij (3)

边界条件：

式中：σ_ij表示应力分量，ε_ij表示应变分量，ρ为密度，b_i为体积力，L为计算微分算子，D为弹性矩阵，n_j表示边界法向量，T_i为边界载荷力，为边界位移。

把式(2)、式(3)、式(4)带入式(1)，联立上述方程得到动量方程的等效积分弱形式方程：

上式中，δu_i,j为i方向位移对x_j坐标求偏导，V为求解体积域，s为边界域， A为面积域，u_ij为ij方向的位移，u_i为i方向位移，为i方向的加速度，为 S_σ边界域上的外力。

步骤6：将所述全场位移场函数带入所述等效积分弱形式方程求解，得到刚度矩阵、等效内力和外力矩阵。将重构的位移场函数u^R(x)带入等效积分弱形式方程，并将该方程写成矩阵形式得到下式：

式中：M为质量矩阵，为加速度场，f^int为等效内力矩阵，f^ext为等效外力矩阵。

步骤7：利用中心差分法对等效积分弱形式方程进行代数求解，得到不同时间步的形变和应力。

所述中心差分法将时间平均分为n个时间间隔，0，Δt，…t时刻的位移、速度、加速度已知，求解t+Δt时刻对应的解；

t时刻的加速度和速度利用中心差分法表示为

式中，u_t-Δt(x)为离散颗粒在t-Δt时刻的位移，u_t(x)为离散颗粒在t时刻的位移，u_t+Δt(x)为离散颗粒在t+Δt时刻的位移；

t+Δt时刻的位移解答可以根据时间t时刻的运动方程得到，即为下式

式中，M、C和K分别代表影响域的质量矩阵、阻尼矩阵和刚度矩阵；F_t是 t时刻外力和；

由此可得到中心差分法的递推公式

若t与t-Δt时刻的位移已知，则可以进一步求解出t+Δt时刻的位移。

需要注意的是，该算法开始计算时，即t＝0时，需要知道u_-Δt(x)，为此利用和的计算公式可以得到

其中，u₀(x)和可以在初始条件中设置，而可利用初始时运动方程得到，即为

利用步骤5中的几何方程，即式(2)，将不同时刻位移带入计算，得到形变，即应变，在利用步骤5中的物理方程，即式(3)求解，将应变带入计算得到应力。

实施例：

再生核粒子方法为理解大变形机理提供了有效的手段。下面以橡胶大变形为例说明再生核粒子在解决大变形问题的优越性。

以2D方形橡胶块的压缩为例，橡胶块尺寸：长a＝1inch，宽b＝2inch；橡胶块底端固定约束，顶端受到垂向方向的压缩载荷，加载速度为V0＝50inch/s。橡胶块材料为超弹性材料，材料本构选择Mooney-Rivilin本构，具体材料属性为：弹性模量E＝1.384MPa，泊松比nu＝0.5。

采用有限元方法对大变形压缩问题进行仿真，随着橡胶块被压缩，网格发生严重畸变，长宽比急剧增大，当压缩量达到60％时，在大变形区由于网格畸变导致计算无法进行；而基于再生核粒子的无网格物理变形仿真方法不依赖于网格，即使压缩量达到90％，依然未出现畸变现象，且仿真结果是光滑连续的。由此，验证了再生核粒子算法可以进行大变形压缩情况的收敛计算，且采用再生核粒子算法计算大变形问题较有限元方法有很大的优势。

以上对本发明的具体实施例进行了详细描述，但其只是作为范例，本发明并不限制于以上描述的具体实施例。对于本领域技术人员而言，任何对本发明进行的等同修改和替代也都在本发明的范畴之中。因此，在不脱离本发明的精神和范围下所作的均等变换和修改，都应涵盖在本发明的范围内。

Claims (10)

1.一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：包括如下步骤：

(1)按照给定的粒子密度对仿真几何模型进行粒子离散，得到离散颗粒，确定所述离散颗粒的初始位置、影响域Ω和影响域半径ρ；

(2)基于再生核粒子算法，选用修正函数和窗函数的类型，构造所述影响域内所述离散颗粒的形函数；

(3)根据所述影响域Ω内的所述离散颗粒的位移和所述形函数构造全场位移场函数；

(4)利用边界转化方法进行边界处理，消除边界不一致性引起的误差；

(5)联立动量方程和几何方程、物理方程、边界条件构造等效积分弱形式方程；

(6)将所述全场位移场函数带入所述等效积分弱形式方程求解，得到刚度矩阵、等效内力和外力矩阵；

(7)利用中心差分法对等效积分弱形式方程进行代数求解，得到不同时间步的形变和应力。

2.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：所述离散颗粒的几何形状包括、但不限于圆形、球形、矩形、立方体。

3.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于，所述步骤(2)中，基于再生核粒子法，所述形函数Φ_ρ(x,x-y)为：

Φ_ρ(x；x-y)＝C_ρ(x；x-y)k_ρ(x-y)

式中，C_ρ(x,x-y)为满足再生条件的修正函数，k_ρ(x-y)为满足再生条件的窗函数；

所述修正函数C_ρ(x,x-y)是由多项式基函数的线性组合表示而成，如下：

C_ρ(x；x-y)＝C₀(x)+C₁(x)(x-y)+C₂(x)(x-y)²+…+C_N(x)(x-y)^N

＝H^T(x-y)C(x)

其中，N是满足再生条件的次数，H^T(x-y)是N次多项式基函数的向量；C(x)为系数向量，通过施加N次再生条件来确定；C_i(x)为系数(i＝1,2,3……N)；

H^T(x-y)＝[1,x-y,(x-y)²,…,(x-y)^N]

C(x)＝[C₀(x),C₁(x),…,C_N(x)]

所述窗函数κ_ρ(x-y)是二维窗函数：

其中，为中间变量。

4.根据权利要求3所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于，所述窗函数选用三次样条函数，如下：

式中，

5.根据权利要求3所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：所述窗函数包括、但不限于三次样条函数、B样条函数、高斯函数。

6.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：所述步骤(3)中，所述全场位移场函数u^R(x)为：

式中，u(y)表示所述影响域Ω内的所述离散颗粒的位移；

将影响域Ω用NP个所述离散粒子离散，获得所述全场位移场函数的离散近似式：

式中，NP表示所述影响域Ω内的所述离散粒子的总点数，u_I对应所述离散粒子I处的位移,ΔV_I为离散粒子的体积权重，N_I(x)为所述离散粒子I处的形函数，I＝1，2，……NP。

7.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：所述步骤(4)中，所述边界转化方法将所述全场位移场函数u^R(x)分为两部分

式中，b表示位移边界上的离散颗粒，nb表示除去位移边界上的离散颗粒剩下的所有离散颗粒；位移边界上的离散颗粒总数为N_b，N_nb＝NP-N_b；N^b(x)为边界插值点上的形函数，N^nb(x)为非边界插值点的形函数，u^b为边界插值点位移，u^nb为非边界插值点位移；

通过强制位移边界条件消除误差，假定

u^b＝(D^b)^-1g

式中：g表示边界位移，D^b通过下式构造

8.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于，所述步骤(5)中：

动量方程：

几何方程：ε_ij＝Lu_ij (2)

物理方程：σ_ij＝Dε_ij (3)

边界条件：

式中：σ_ij表示应力分量，ρ为密度，b_i为体积力，ε_ij表示应变分量，L为计算微分算子，u_ij为位移，D为弹性矩阵，n_j表示边界法向量，T_i为边界载荷力，S_σ为载荷边界域，S_u为位移边界域，ν_i为边界位移，为边界位移；

把式(2)、式(3)、式(4)带入式(1)，联立上述方程得到动量方程的等效积分弱形式方程：

式中，δu_i,j为i方向位移对x_j坐标求偏导，V为求解体积域，s为边界域，A为面积域，u_ij为ij方向的位移，u_i为i方向位移，为i方向的加速度，为S_σ边界域上的外力。

9.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于，所述步骤(6)中，将重构的位移场函数u^R(x)带入等效积分弱形式方程，并将该方程写成矩阵形式得到下式：

式中，M为质量矩阵，为加速度场，f^int为等效内力矩阵，f^ext为等效外力矩阵。

10.根据权利要求1所述的一种基于再生核粒子的无网格物理变形仿真方法，其特征在于：所述步骤(7)中，所述中心差分法将时间平均分为n个时间间隔，0，Δt，…t时刻的位移、速度、加速度已知，求解t+Δt时刻对应的解；

t时刻的加速度和速度利用中心差分法表示为

式中，u_t-Δt(x)为离散颗粒在t-Δt时刻的位移，u_t(x)为离散颗粒在t时刻的位移，u_t+Δt(x)为离散颗粒在t+Δt时刻的位移；

t+Δt时刻的位移解答可以根据时间t时刻的运动方程得到，即为下式

式中，M、C和K分别代表影响域的质量矩阵、阻尼矩阵和刚度矩阵；F_t是t时刻外力和；

由此可得到中心差分法的递推公式

若t与t-Δt时刻的位移已知，则可以进一步求解出t+Δt时刻的位移；利用步骤(5)中的几何方程，将不同时刻位移带入计算，得到形变，即应变，再利用步骤(5)中的物理方程求解，将应变带入计算得到应力。