CN109145386A - 一种用于工程结构的地震碰撞反应分析方法 - Google Patents

Abstract

本发明公开了一种用于工程结构的地震碰撞反应分析方法，涉及土木工程技术领域，具体包括：结构弹性状态时采用全量方程精细积分公式递推求解，塑性状态采用增量方程精细积分公式递推求解；弹塑性状态转换时，求解精确界点转换时刻，界点转换后的第一个时间步采用全量方程迭代求解；对当前时间步判断结构是否产生接触碰撞，如产生碰撞，求解精确的碰撞时刻，从该时刻点开始引入碰撞力，执行碰撞阶段动力方程的递推计算，直至碰撞结束；结构碰撞状态的积分步长取值较小，分离状态的积分步长取值较大，进一步提高计算效率。与传统的隐式或显式积分方法相比，该方法无条件稳定，计算效率高，精度接近解析解，能精细模拟地震碰撞过程。

Description

一种用于工程结构的地震碰撞反应分析方法

技术领域

本发明涉及土木工程技术领域，特别是涉及一种用于工程结构的地震碰撞反应分析方法。

背景技术

地震是人类面对的最严重的的自然灾害之一，因其巨大的破坏力往往使人类遭受巨大的生命财产损失。工程结构的抗震历来受到人们的重视。随着人们对桥梁、房屋建筑等地震破坏现象和机理的不断认识，抗震设计的理论和实践也在不断深化。在众多工程结构的地震破坏现象中，由于相邻结构位移超过结构预留间隙而产生的碰撞现象以及由于碰撞导致的结构破坏是非常普遍的。

工程结构的碰撞破坏现象引起了国内外众多专家学者的重视，国外从上世纪90年代初开始已有很多学者对结构地震碰撞问题开展了广泛的研究，一些成果已经反映到设计规范中。如我国《公路桥梁抗震设计细则》规定桥梁抗震采用两阶段设计，桥梁结构的动力计算模型应考虑相邻结构和边界条件的影响。

要获得工程结构地震中的反应，最直接的办法就是采用直接积分法进行地震时程反应分析。要考虑相邻结构碰撞的影响，首先要对碰撞作用建立合理的力学模型。碰撞是状态非线性问题，其特点是相邻结构位移差超过结构预留间隙时发生，持续时间很短，瞬时冲击力很大。目前，考虑碰撞的时程反应分析大多是在有限元模型中引入接触单元，建立考虑碰撞力的动力平衡方程，并采用隐式或显式的直接积分法进行求解。

已发展的直接积分算法有中心差分法、Newmark-β法、wilson-θ法、Houbolt法等。在隐式算法中，Newmark-β法使用更为广泛。与Newmark-β法相比，wilson-θ法、Houbolt法在积分步长较大时，会引起高频模态反应的衰减，精度较差。隐式算法是无条件稳定的或其稳定性条件很容易满足，因而积分步长可以取得较大，但是对于材料非线性问题，在每一积分步内都需要对动力平衡方程进行迭代求解，且每次迭代都需要求解联立方程组。这个过程不仅需要占用大量的计算资源和存储资源，也可能会遇到迭代过程不收敛或联立方程组因出现病态而无确定解的问题。

显式算法采用动力学方程的一些差分格式，最常用的为中心差分法，不需要进行平衡迭代，也不需要求解联立方程组，计算速度快，时间步长只要取的足够小，其解是收敛的。Ansys Ls-dyna、Abaqus有限元软件均采用了中心差分格式的显式求解算法，适合对爆炸、冲击、接触碰撞等问题进行求解。但该法也存在缺点，一是该法是条件稳定的，时间步长受到数值积分稳定性的限制，使用时时间步长要取的足够小，使积分点数大大增加；二是该法在质量矩阵为对角阵时才能发挥出优势。

传统的隐式算法或显式算法只具有二阶或三阶精度，算法的理论基础决定其具有“算法阻尼”，导致系统的能量耗散，具体表现为随着积分时间增加的周期延长和振幅衰减。已经从理论上证明，传统数值方法导致能量耗散的根本原因是其具有耗散的截断项。

发明内容

本发明实施例提供了一种用于工程结构的地震碰撞反应分析方法，可以解决现有技术中的问题。

本发明提供了一种用于工程结构的地震碰撞反应分析方法，该方法包括以下步骤：

步骤1，准备阶段：输入工程结构参数、弹塑性模型参数、接触单元模型及参数、地震波文件，以及分离阶段的积分时间步长dt₁和碰撞阶段的积分时间步长dt₂、迭代容差；按照有限元建模，建立弹性状态结构的割线刚度矩阵、质量矩阵和阻尼矩阵；按照弹性阶段的割线刚度矩阵和初始位移状态，确定初始速度和加速度，初始化积分时长t；

步骤2，弹塑性积分计算：以积分时长t是否超出地震波时长为条件判断地震是否结束，如果没有结束则对每一个时间步进行以下计算：

判断当前时间步的结构弹塑性状态：如果结构为弹性状态，则采用全量法计算当前时间步内的结构状态，并在计算完成后判断该时间步内是否产生界点，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3；

如果结构为塑性状态，则判断当前时间步是否为界点后的第一个时间步或切线刚度是否为0，如果是，则采用全量法计算当前时间步内的结构状态，否则按照当前时刻的切线刚度建立结构刚度矩阵，然后采用增量法计算当前时间步的结构状态；计算完成后判断当前时间步内是否产生界点，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3；

步骤3，碰撞阶段积分计算：判断当前时间步内是否产生接触碰撞，若未产生，直接执行下一时间步；若产生接触碰撞，则舍弃当前时间步的计算结果，求解精确的碰撞时刻，从该时刻点开始积分步长变为dt₂，并引入碰撞力，根据碰撞点的位置建立碰撞力向量，然后将碰撞力向量分解成刚度矩阵与位移矩阵乘积+阻尼矩阵与速度向量乘积形式，对总体刚度矩阵、阻尼矩阵进行修正；然后采用全量法精确求解该碰撞时刻的结构反应，直到碰撞点分离；

若当前时间步结构进入分离状态，则采用全量法求解精确的分离时刻和该时刻的结构状态；

步骤4，数据处理阶段：弹塑性和碰撞阶段的积分计算完成后，提取结构的位移反应、速度反应或加速度反应并绘制时程曲线；积分计算过程中，对各接触单元的侵彻位移、碰撞力进行存储，绘制接触单元的碰撞力时程或侵彻位移-碰撞力关系曲线；对结构的恢复力进行存储，绘制弹塑性结构的滞回曲线。

本发明实施例中的一种用于工程结构的地震碰撞反应分析方法，与传统的隐式算法或显示算法相比，首先，精细积分法是无条件稳定的，结果精度接近解析解，确保了地震碰撞时程分析的精确性；其次，通过界点的精确计算和全量法、增量法的联合使用，确保结构滞回曲线完全符合弹塑性模型假定；再次，通过碰撞、分离时刻的精确计算以及分离状态和碰撞状态不同步长设定，不仅能精确模拟地震碰撞过程，而且最大限度的减少积分点增幅，算法仍具有极高的计算效率，解决了传统算法效率与精度间的矛盾。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例中提供的用于工程结构的地震碰撞反应分析方法的流程图；

图2为相邻结构的多点碰撞示意图；

图3为常用的接触单元模型示意图；

图4为常用的折线形弹塑性模型示意图，其中(a)为不考虑刚度退化的双折线模型，(b)为考虑刚度退化的双折线模型，(c)为武田三折线模型。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1，本发明实施例提供了一种用于工程结构的地震碰撞反应分析方法，该方法包括以下步骤：

步骤1，准备阶段。输入工程结构参数、弹塑性模型参数、接触单元模型及参数、地震波文件，以及分离阶段的积分时间步长dt₁和碰撞阶段的积分时间步长dt₂、迭代容差；

按照有限元建模，建立弹性状态结构的割线刚度矩阵、质量矩阵和阻尼矩阵；

按照弹性阶段的割线刚度矩阵和初始位移状态，确定初始速度和加速度；

初始化积分步数s、积分时长t、弹塑性阶段标识数组和碰撞状态标识数组。

步骤2，弹塑性积分计算。以积分时长t是否超出地震波时长为条件判断地震是否结束，如果没有结束则对每一个时间步进行以下计算：

判断当前时间步的结构弹塑性状态：如果结构为弹性状态，则采用全量法计算当前时间步内的结构状态，并在计算完成后判断该时间步内是否产生界点，即弹塑性状态是否发生改变，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻t_a和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3；

如果结构为塑性状态，则判断当前时间步是否为界点后的第一个时间步或切线刚度是否为0，如果是，则采用全量法计算当前时间步内的结构状态，否则按照当前时刻的切线刚度建立结构刚度矩阵，然后采用增量法计算当前时间步的结构状态。计算完成后判断当前时间步内是否产生界点，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻t_a和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3。

该步骤中使用的全量法和增量法的计算过程如下：

全量动力平衡方程：

增量动力平衡方程：

如图2所示，上式中，M表示质量矩阵，U表示位移矩阵，表示位移矩阵的二阶微分，即加速度矩阵，表示位移矩阵的一阶微分，即速度矩阵，C表示阻尼矩阵，K表示割线刚度矩阵，I表示单位阵，表示割线加速度。

式(1)-(2)均为二阶微分方程，通过引入对偶变量，可将二阶微分方程降为一阶，写成如下形式：

式中，H为一维状态方程的系统矩阵，f为一维状态方程的非齐次项。

对于全量动力平衡方程降阶后的一维状态方程系统矩阵，构造形式如下：

增量动力平衡方程降阶后的一维状态方程系统矩阵与式(4)相同，只需把式(4)中的割线刚度矩阵替换为切线刚度矩阵。

方程(3)具有如下形式的通解：

式(5)即表征从已知状态到未知状态的递推关系。求解的关键变成了指数矩阵e^H·Δt的求解和非齐次项的积分计算。指数矩阵求解基于指数函数的加法原理，采用2^N类算法进行。当假定非齐次项在时间步内按照线性变化时，f₁为非齐次项在荷载步内的变化率，式(5)的解析解是：

为避免系统矩阵H求逆，式(5)的非齐次项也可以采用级数解：

采用增量法精确求解界点时刻和界点处的结构状态时，使用折线形模型，如图4所示，折线型模型加载或卸载过程中的开裂点、屈服点以位移作为判定条件，加卸载转换界点以速度为判定条件。假设在t_i～t_i+1时段产生界点，两时刻的结构位移、速度、加速度分别为设界点产生的精确时刻为t_i+τ，界点判定位移为u_p，根据线性加速度假定，可建立如下方程：

而加卸载转换界点，判定条件是t_i+τ时刻速度为0，可建立如下方程：

方程(8)为一元三次方程，运用数值算法求解[0Δt]内的解。方程(9)为一元二次方程，可直接求解[0Δt]内的解。t_i+τ时刻即为精准界点产生时刻，进而舍弃t_i+1时刻求解结果，以精细积分法求解t_i+τ时刻结果。

采用增量法时，由于通过界点后的第一个时间步的增量状态未知，因此需要通过全量法迭代求解该时间步，为增量法的启动创造初始条件。具体步骤如下：

第1步，第i积分点的割线刚度和该积分点的结构反应u_i、为已知，根据第i积分点的割线刚度，按照全量精细积分递推公式递推第i+1积分点的结构反应u_i+1、

第2步，根据弹塑性模型滞回规律确定第i+1积分点的割线刚度矩阵K_i+1；

第3步，按照下式计算第i+1积分点的不平衡力：

第4步，根据第2步确定的第i+1积分点的割线刚度K_i+1和第i积分点的结构反应u_i、按照全量精细积分递推公式递推第i+1积分点的结构反应u_i+1、

第5步，根据第4步计算结果，修正第i+1积分点的割线刚度K_i+1；

第6步，执行第3-5步，直至不平衡力小于设定容差。

步骤3，碰撞阶段积分计算。判断当前时间步内是否产生接触碰撞：若未产生，直接执行下一时间步；

若产生接触碰撞，则舍弃当前时间步的计算结果，求解精确的碰撞时刻t_b，从该时刻点开始积分步长变为dt₂，并引入碰撞力，根据碰撞点的位置建立碰撞力向量，然后将碰撞力向量分解成刚度矩阵与位移矩阵乘积+阻尼矩阵与速度向量乘积形式，对总体刚度矩阵、阻尼矩阵进行修正；然后采用全量法精确求解该碰撞时刻的结构反应，直到碰撞点分离；期间若有新的碰撞点产生或分离，则修正碰撞力向量或刚度矩阵以及阻尼矩阵。

若当前时间步结构进入分离状态，则采用全量法求解精确的分离时刻t_c和该时刻的结构状态，此后转入分离状态求解。

碰撞阶段的全量动力平衡方程为：

其中，F_c表示碰撞力向量，碰撞力的计算按照接触单元模型假定进行。本方法支持多种接触单元模型，包括线性弹簧模型、线性弹簧-阻尼(Kelvin-Voigt)模型、Hertz模型、非线性弹簧阻尼(Hertz-damp)模型、Jan-Hertz-damp模型，如图3所示。几种接触单元模型的碰撞力统一表达成如下形式：

式(12)中，β(t)、c(t)分别表示碰撞弹簧刚度和阻尼系数，g_p为初始间隙，u₁、u₂、为相邻两质点的位移和速度。

对于线性弹簧模型，β(t)为常数，c(t)＝0；对于Kelvin-Voigt模型，β(t)为常数，c(t)为常数，等于0时退化为线性弹簧模型；对于Hertz模型，碰撞弹簧刚度是侵彻位移的函数，而c(t)＝0；对于Hertz-damp模型，是在Hertz模型的基础上计入碰撞阻尼力，该模型阻尼系数为侵彻位移的函数。对于Jan-Hertz-damp模型，恢复阶段：c(t)＝0，接近阶段阻尼系数c(t)为：

碰撞时刻和分离时刻的计算方法为：假定t_i时刻第k个接触单元u_k-u_k+1-g_pk<0(没有接触)，而t_i+1时刻u_k-u_k+1-g_pk>0(已经接触)，说明接触开始时刻t_i+τ介于t_i时刻至t_i+1时刻之间。已解出t_i时刻和t_i+1时刻的位移、速度、加速度，根据线性加速度假定，t_i+τ时刻的位移满足：

方程(13)为一元三次方程，利用数值算法求解τ在区间[0Δt]内的解，t_i+τ时刻即为接触开始时刻。在每一步递推前对所有接触单元循环判定是否有新的接触产生，一旦判定有接触单元从分离状态转换到接触碰撞状态，首先求解精确的碰撞时刻t_i+τ和该时刻的结构反应，舍弃t_i+1的结果。从t_i+τ时刻积分步长变为碰撞阶段的积分步长，进入碰撞阶段求解，接触碰撞状态至分离状态转换的接触结束时刻可以采用同样的原理计算。

步骤4，数据处理阶段。弹塑性和碰撞阶段的积分计算完成后，可提取结构的位移反应、速度反应或加速度反应并绘制时程曲线；积分计算过程中，对各接触单元的侵彻位移、碰撞力进行存储，可绘制接触单元的碰撞力时程或侵彻位移-碰撞力关系曲线；对结构的恢复力进行存储，可绘制弹塑性结构的滞回曲线。

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

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

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

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

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.一种用于工程结构的地震碰撞反应分析方法，其特征在于，该方法包括以下步骤：

步骤1，准备阶段：输入工程结构参数、弹塑性模型参数、接触单元模型及参数、地震波文件，以及分离阶段的积分时间步长dt₁和碰撞阶段的积分时间步长dt₂、迭代容差；按照有限元建模，建立弹性状态结构的割线刚度矩阵、质量矩阵和阻尼矩阵；按照弹性阶段的割线刚度矩阵和初始位移状态，确定初始速度和加速度，初始化积分时长t；

步骤2，弹塑性积分计算：以积分时长t是否超出地震波时长为条件判断地震是否结束，如果没有结束则对每一个时间步进行以下计算：

判断当前时间步的结构弹塑性状态：如果结构为弹性状态，则采用全量法计算当前时间步内的结构状态，并在计算完成后判断该时间步内是否产生界点，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3；

如果结构为塑性状态，则判断当前时间步是否为界点后的第一个时间步或切线刚度是否为0，如果是，则采用全量法计算当前时间步内的结构状态，否则按照当前时刻的切线刚度建立结构刚度矩阵，然后采用增量法计算当前时间步的结构状态；计算完成后判断当前时间步内是否产生界点，如果产生界点则舍弃当前时间步的计算结果，采用增量法精确求解界点时刻和界点处的结构状态，并使用界点处的结构状态代替当前时间步的结构状态，如果没有产生界点或界点处的结构状态计算完成，则进入步骤3；

步骤3，碰撞阶段积分计算：判断当前时间步内是否产生接触碰撞，若未产生，直接执行下一时间步；若产生接触碰撞，则舍弃当前时间步的计算结果，求解精确的碰撞时刻，从该时刻点开始积分步长变为dt₂，并引入碰撞力，根据碰撞点的位置建立碰撞力向量，然后将碰撞力向量分解成刚度矩阵与位移矩阵乘积+阻尼矩阵与速度向量乘积形式，对总体刚度矩阵、阻尼矩阵进行修正；然后采用全量法精确求解该碰撞时刻的结构反应，直到碰撞点分离；

若当前时间步结构进入分离状态，则采用全量法求解精确的分离时刻和该时刻的结构状态；

步骤4，数据处理阶段：弹塑性和碰撞阶段的积分计算完成后，提取结构的位移反应、速度反应或加速度反应并绘制时程曲线；积分计算过程中，对各接触单元的侵彻位移、碰撞力进行存储，绘制接触单元的碰撞力时程或侵彻位移-碰撞力关系曲线；对结构的恢复力进行存储，绘制弹塑性结构的滞回曲线。

2.如权利要求1所述的用于工程结构的地震碰撞反应分析方法，其特征在于，步骤2中使用的全量法和增量法的计算过程如下：

全量动力平衡方程：

增量动力平衡方程：

上式中，M表示质量矩阵，U表示位移矩阵，表示加速度矩阵，表示速度矩阵，C表示阻尼矩阵，K表示割线刚度矩阵，I表示单位阵，表示割线加速度；

在式(1)-(2)中引入对偶变量，将二阶微分方程降为一阶，写成如下形式：

式中，H为一维状态方程的系统矩阵，f为一维状态方程的非齐次项；

对于全量动力平衡方程降阶后的一维状态方程系统矩阵，构造形式如下：

增量动力平衡方程降阶后的一维状态方程系统矩阵与式(4)相同，只需把式(4)中的割线刚度矩阵替换为切线刚度矩阵；

方程(3)具有如下形式的通解：

3.如权利要求1所述的用于工程结构的地震碰撞反应分析方法，其特征在于，步骤2中采用增量法精确求解界点时刻和界点处的结构状态时，弹塑性模型使用折线形模型，假设在t_i～t_i+1时段产生界点，两时刻的结构位移、速度、加速度分别为设界点产生的精确时刻为t_i+τ，界点判定位移为u_p，根据线性加速度假定建立如下方程：

界点判定条件是t_i+τ时刻速度为0，建立如下方程：

方程(6)为一元三次方程，运用数值算法求解[0Δt]内的解，方程(7)为一元二次方程，可直接求解[0Δt]内的解，t_i+τ时刻即为精准界点产生时刻，进而舍弃t_i+1时刻求解结果，以精细积分法求解t_i+τ时刻结果。

4.如权利要求1所述的用于工程结构的地震碰撞反应分析方法，其特征在于，步骤2中采用增量法计算界点时刻和界点处的结构状态时，由于通过界点后的第一个时间步的增量状态未知，因此需要通过全量法迭代求解该时间步，为增量法的启动创造初始条件，具体步骤如下：

第1步，第i积分点的割线刚度和该积分点的结构反应u_i、为已知，根据第i积分点的割线刚度，按照全量精细积分递推公式递推第i+1积分点的结构反应u_i+1、

第2步，根据弹塑性模型滞回规律确定第i+1积分点的割线刚度矩阵K_i+1；

第3步，按照下式计算第i+1积分点的不平衡力：

第4步，根据第2步确定的第i+1积分点的割线刚度K_i+1和第i积分点的结构反应u_i、按照全量精细积分递推公式递推第i+1积分点的结构反应u_i+1、

第5步，根据第4步计算结果，修正第i+1积分点的割线刚度K_i+1；

第6步，执行第3-5步，直至不平衡力小于设定容差。

5.如权利要求1所述的用于工程结构的地震碰撞反应分析方法，其特征在于，步骤3中碰撞时刻的计算方法为：假定t_i时刻第k个接触单元u_k-u_k+1-g_pk<0，即没有接触，而t_i+1时刻u_k-u_k+1-g_pk>0，即已经接触，说明接触开始时刻t_i+τ介于t_i时刻至t_i+1时刻之间，已解出t_i时刻和t_i+1时刻的位移、速度、加速度，根据线性加速度假定，t_i+τ时刻的位移满足：

方程(9)为一元三次方程，利用数值算法求解τ在区间[0Δt]内的解，t_i+τ时刻即为接触开始时刻。