CN115270515B - 基于边界面模型的土体状态数据预测方法、装置和设备 - Google Patents

Abstract

本发明公开了一种基于边界面模型的土体状态数据预测方法、装置和设备，该方法首先对每个增量步进行子增量步划分，基于当前子增量步的步长和初始土体状态数据进行八次应力试算，得到八组土体增量数据；然后对八组土体增量数据进行加权求和，得到当前子增量步结束时的结束土体状态数据；之后基于结束时土体状态数据和八组土体增量数据计算局部截断误差；若局部截断误差不小于预设阈值，则缩小当前子增量步的步长，并返回基于当前子增量步的步长和初始土体状态数据进行八次应力试算的步骤，重新计算。本发明提供的技术方案，进一步提高了各个子增量步的截断误差控制精度，从而进一步提高了土体状态数据的计算准确度。

Description

基于边界面模型的土体状态数据预测方法、装置和设备

技术领域

本发明涉及岩土力学领域，具体涉及一种基于边界面模型的土体状态数据预测方法、装置和设备。

背景技术

海洋资源的有效开发与利用，已经成为当今全球共同关注的焦点问题。由于海洋环境的特殊性，需要在海上安装或建造各类构筑物才能保障在海洋环境中开展的生产活动，此类构筑物不仅要承受结构自重，还要承受风、浪、流等导致的循环荷载。在国际范围内，想要有效地分析构筑物循环承载特性，需要建立分析构筑物循环承载的数值模型，若要建立此类数值模型，则需要先建立一个可以反映土体循环动力特性的本构模型，然后提出与之匹配的应力积分流程。

目前，专利文件CN113779775A提出的一种边界面塑性模型(通俗地讲，边界面塑性模型或称为边界面模型，反映了土体应力与应变之间的关系)，在模拟土体循环动力特性的领域已有一定的应用。然而该边界面塑性模型本质上属于增量型本构模型，增量可以近似理解为在增量步(一段微小时间步长)内的变化量，例如给定一个应变增量Δε，根据本构关系可以计算出相应的应力增量Δσ，而实际应用时需要应力应变的全量关系，那么就需要基于边界面塑性模型的增量关系计算全量关系。由于增量型本构关系的数学表达式本质上属于常微分方程组，因此在绝大多数情况下无法通过直接积分得到精确的应力应变全量关系的解析表达式。求解常微分方程的方法通常可分为隐式算法和显式算法两大类，与之对应的求解本构方程组的应力积分算法也可分为隐式积分算法和显式积分算法。

与隐式积分算法相比，显式积分算法没有迭代求解多元非线性方程组的繁琐过程，其计算效率更高，但由于局部截断误差的存在，计算过程产生的累积误差将会更大，导致显式算法的精度相对较低。但是，在应用显示积分算法时，只要控制应变增量步的大小在一定范围内，则能够提高该算法的收敛性，并提高其计算精度。又因为该算法的效率较高，因此显式积分算法不失为一种较为理想的计算土体状态数据的办法。

基于此，专利文件CN114492094A公开了一种土体状态数据的显式积分计算方法，该方法将一个增量步再次划分为多个子增量步，从第一个子增量步开始，获取该子增量步初始的土体状态数据，然后基于上述边界面塑性模型进行两次应力试算，并根据两次应力试算的结果来判断第一个子增量步的长度划分是否合理，是否使积分收敛，如果不合理则重新划分第一个子增量步的长度，如果合理则计算第一个子增量步结束时的土体状态数据，然后将第一个子增量步结束时的土体状态数据作为下一个子增量步的初始土体状态数据，计算下一个子增量步，直至一个增量步的长度被多个子增量步划分完为止，最后一个子增量步结束时的土体状态数据即为当前增量步结束时的土体状态数据。通过上述步骤，从而实现了将子增量步的大小控制在一定范围内，结合本构关系计算各个子增量步结束时的土体状态数据，并基于各个子增量步结束时的土体状态数据计算各个增量步结束时的土体状态数据，最终完成迭代。后续利用各个增量步的结束时的土体状态数据即可进行应力应变的全量关系计算。

在此基础上，只要使各个子增量步的截断误差控制更加精确，计算土体状态数据的准确度会更高。因此，如何提高各个子增量步的截断误差控制精度，从而进一步提高计算土体状态数据的准确度，是值得研究的问题。

发明内容

有鉴于此，本发明实施方式提供了一种基于边界面模型的土体状态数据预测方法、装置和设备，提高了各个子增量步的截断误差控制精度，从而进一步提高了计算土体状态数据的准确度。

根据第一方面，本发明实施例提供了一种基于边界面模型的土体状态数据预测方法，所述方法包括：对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据；基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，所述土体增量数据中包括应力增量和内变量增量；计算所述初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，所述目标土体增量数据是基于所述八组土体增量数据进行加权求和的结果；基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，并判断所述局部截断误差是否小于预设阈值；若所述局部截断误差不小于所述预设阈值，则缩小所述当前子增量步的步长，并返回所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算的步骤，重新计算。

可选地，所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，包括：基于当前子增量步的步长计算当前子增量步的应变增量；若当前是第一次应力试算，则基于所述初始土体状态数据中的应力值、内变量值和所述应变增量计算得到应力增量，并基于所述初始土体状态数据和硬化准则计算得到内变量增量；若当前不是第一次应力试算，则将所述初始土体状态数据中应力值和内变量值分别与融合应力增量和融合内变量增量对应相加，并利用相加后的应力值、内变量值和所述应变增量计算本次试算的应力增量，所述融合应力增量和所述融合内变量增量分别是前面每次应力试算得到的应力增量和内变量增量按照影响比例缩小后的加权融合值；将所述初始土体状态数据中的应变值与影响应变增量相加，得到调整后的土体状态数据，并利用所述调整后的土体状态数据和硬化准则计算得到本次试算的内变量增量，所述影响应变增量是按照所述影响比例缩小后的应变增量；其中，所述影响比例从第二次试算开始逐次递增，直至第八次试算递增为1。

可选地，计算所述融合应力增量和所述融合内变量增量的步骤包括：将前面每次应力试算得到的应力增量和内变量增量按照影响比例缩小；将每次应力试算缩小后的应力增量加权融合得到所述融合应力增量，并将每次应力试算缩小后的内变量增量加权融合得到融合内变量增量；其中，各次试算对应的融合权重随试算次数呈递增顺序。

可选地，所述影响比例从第二次试算开始逐次递增，直至第八次试算递增为1，包括：从第二次试算至第八次试算，对应的影响比例为

1、1。

可选地，所述基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，包括：通过下式计算所述局部截断误差，

式中，R^(k)表示第k个子增量步的局部截断误差，

表示第n个增量步中第k个子增量步的结束土体状态数据中的应力值，/>

表示第n个增量步中第k个子增量步的结束土体状态数据中的内变量值，Δσ₁～Δσ₈分别是试算得到的八组土体增量数据中的应力增量、Δκ₁～Δκ₈分别是八组土体增量数据中的内变量增量，E_σ和E_κ分别是计算中间量。

可选地，基于所述八组土体增量数据加权求和得到所述目标土体增量数据的步骤包括：通过下一子增量步的预测应力值对当前子增量步初始土体状态数据中的应变值和应力值进行泰勒展开，并对比同类项系数确定八组土体增量数据的目标权重集合，所述目标权重集合中同种增量的权重值的和为1；利用所述目标权重集合与所述八组土体增量数据进行加权求和，得到所述目标土体增量数据。

可选地，所述方法还包括：若所述局部截断误差小于所述预设阈值，则将所述结束土体状态数据作为下一子增量步的初始土体状态数据进行迭代计算，直至当前增量步的步长被多个子增量步划分完毕，并将最后一个子增量步的结束土体状态数据作为当前增量步的结束土体状态数据进行保存。

根据第二方面，本发明实施例提供了一种基于边界面模型的土体状态数据预测装置，所述装置包括：初始化模块，用于对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据；试算模块，用于基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，所述土体增量数据中包括应力增量和内变量增量；土体状态数据预测模块，用于计算所述初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，所述目标土体增量数据是基于所述八组土体增量数据进行加权求和的结果；误差分析模块，用于基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，并判断所述局部截断误差是否小于预设阈值；步长修正模块，用于若所述局部截断误差不小于所述预设阈值，则缩小所述当前子增量步的步长，并返回所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算的步骤，重新计算。

根据第三方面，本发明实施例提供了一种基于边界面模型的土体状态数据预测设备，包括：存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器中存储有计算机指令，所述处理器通过执行所述计算机指令，从而执行第一方面，或者第一方面任意一种可选实施方式中所述的方法。

根据第四方面，本发明实施例提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于使所述计算机从而执行第一方面，或者第一方面任意一种可选实施方式中所述的方法。

本申请提供的技术方案，具有如下优点：

本申请提供的技术方案，针对当前子增量步，基于初始土体状态数据进行八次应力试算得到八个应力增量和八个内变量增量的增量数据，然后对八组数据进行加权得到目标土体增量数据，相比两次应力试算取平均的方法，增量数据准确度更高。然后将初始土体状态数据与目标土体增量数据相加，得到当前子增量步结束时的结束土体状态数据。然后基于结束土体状态数据和八组土体增量数据计算局部截断误差，使得局部截断误差更小，更容易满足误差判别条件。如果局部截断误差小于预设阈值则保留当前子增量步的结束土体状态数据，如果局部截断误差大于预设阈值则缩小当前子增量步的步长，重新计算当前子增量步的结束土体状态数据。从而通过几轮试算，进一步提高了当前子增量步的截断误差控制精度。然后按照上述方案对每一子增量步进行迭代，得到每一个子增量步结束时的结束土体状态数据，进而可以递推出更加准确的当前增量步结束时的土体状态数据。同理按照上述计算步骤，提高了各个增量步预测的土体状态数据的准确度，进而提高土体应力应变全量关系的计算准确度。

此外，由于边界面模型中硬化规律通常与应变相关，硬化规律矩阵受应变增量大小影响，因此在迭代试算过程中，本发明实施例逐步增加应变增量对硬化规律矩阵的影响，使得每次试算的增量数据更加准确，进而提高每个子增量步计算结束土体状态数据的准确率，在进行局部截断误差判别时，更容易收敛，从而最终进一步提高了创建应力应变全量关系的准确度。

附图说明

通过参考附图可以更加清楚地理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：

图1示出了现有技术中基于两次应力试算预测土体状态数据的几何示意图；

图2示出了本发明一个实施方式中一种基于边界面模型的土体状态数据预测方法的步骤示意图；

图3示出了本发明一个实施方式中一种基于边界面模型的土体状态数据预测方法的流程示意图；

图4示出了本发明一个实施方式中一种基于边界面模型的土体状态数据预测装置的结构示意图；

图5示出了本发明一个实施方式中一种基于边界面模型的土体状态数据预测设备的结构示意图。

具体实施方式

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

为了更好的解释本发明实施方式相比现有技术的改进点，下面首先对现有技术基于两次试算调整子增量步步长进行土体状态数据计算的过程和相关术语进行简单介绍：

将欧拉积分格式与多步法相结合，进行根据误差自动划分子增量步的显式积分算法。该算法将单一增量步划分为多个子增量步，并根据误差自动调整子增量步长，该算法的几何意义如图1所示。假设第k个子增量步初始时的土体状态数据包括应力σ_(k)、应变ε_(k)，则此时弹塑性刚度矩阵可根据应力σ_(k)和应变ε_(k)计算得到，以D₁ ^ep标记。图1中D₁ ^ep即为应力应变曲线在状态点(σ_(k)，ε_(k))处的斜率，若此处以切线代替曲线，则通过子增量步的应变增量步长Δε_(k)进行第一次应力试算，得到第一个应力增量Δσ₁，第一次试算后的土体状态数据更新为σ_(k)+Δσ₁和ε_(k+1)，基于边界面塑性模型计算更新后土体状态数据的刚度矩阵为D₂ ^ep，然后以D₂ ^ep取代第一次试算的刚度矩阵D₁ ^ep(图中将D₂ ^ep位置的切线下移到初始位置)，进行第二次应力试算，得到第二次试算的应力增量Δσ₂。此时以两次试算的应力增量的均值对应力状态数据进行更新，即σ_(k+1)＝σ₍k₎+1/2(Δσ₁+Δσ₂)，内变量的预测过程类似，在此不再赘述。上述过程即为基于两次应力试算通过上一子增量步的土体状态数据预测下一子增量步土体状态数据的过程，其几何意义为在每个子增量步以折线代替积分曲线，从而利用多段折线逼近真实的全量应力应变关系曲线。由于折线法存在局部截断误差，需限制子增量步长Δε_(k)的大小才能使折线代替积分曲线的误差足够小。因此，显式积分算法需要根据两次试算的应力增量差值Δσ₁-Δσ₂的大小，判断子增量步长是否满足误差要求，若不满足要求，则根据应力增量差值缩小子增量步长然后重新计算，若满足要求，则保存本次预测的数据并进入下一子增量步的计算。本申请的技术方案，则是在上述技术方案的基础上，进行了改进。

请参阅图2和图3，在一个实施方式中，本发明实施例提供的一种基于边界面模型的土体状态数据预测方法，具体包括以下步骤：

步骤S101：对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据。

步骤S102：基于当前子增量步的步长和初始土体状态数据进行八次应力试算，得到八组土体增量数据，土体增量数据中包括应力增量和内变量增量。

步骤S103：计算初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，目标土体增量数据是基于八组土体增量数据进行加权求和的结果。

步骤S104：基于结束土体状态数据和八组土体增量数据计算局部截断误差，并判断局部截断误差是否小于预设阈值。

步骤S105：若局部截断误差不小于预设阈值，则缩小当前子增量步的步长，并返回基于当前子增量步的步长和初始土体状态数据进行八次应力试算的步骤。

具体地，在本实施例中，针对第k个子增量步，基于初始土体状态数据进行第一次应力试算，然后利用试算的第一应力增量和第一内变量增量更新初始土体状态数据，得到用于第二次试算的土体状态数据，然后利用第二次试算的土体状态数据进行第二次应力试算，得到第二应力增量和第二内变量增量，再次更新土体状态数据，直至进行到第八次应力试算，最终得到八个应力增量和八个内变量增量的增量数据。之后对八个应力增量进行加权融合，并对八个内变量增量进行加权融合，得到目标土体增量数据。相比两次应力试算取平均的方法，增量数据准确度更高。然后将初始土体状态数据与目标土体增量数据相加，预测得到当前子增量步结束时的结束土体状态数据。之后开始进行误差判断步骤：基于结束土体状态数据和八组土体增量数据计算局部截断误差，使得局部截断误差更小，更容易满足误差判别条件，如果局部截断误差小于预设阈值，则保留当前子增量步的结束土体状态数据，进行下一个子增量步计算。如果局部截断误差大于预设阈值则缩小调整当前子增量步的步长，重新计算当前子增量步的结束土体状态数据。从而通过几轮试算，进一步提高了当前子增量步的步长精度。然后按照上述方案对每一子增量步进行迭代，得到每一个子增量步结束时的结束土体状态数据，进而可以递推出当前增量步更加准确的结束土体状态数据。同理按照上述计算步骤，提高了各个增量步预测的土体状态数据的准确度，进而提高后续创建应力应变全量关系的准确度。

具体地，在一实施例中，上述步骤S102，具体包括如下步骤：

步骤一：基于当前子增量步的步长计算当前子增量步的应变增量。

步骤二：若当前是第一次应力试算，则基于初始土体状态数据中的应力值、内变量值和应变增量计算得到应力增量，并基于初始土体状态数据和硬化准则计算得到内变量增量。

步骤三：若当前不是第一次应力试算，则将初始土体状态数据中应力值和内变量值分别与融合应力增量和融合内变量增量对应相加，并利用相加后的应力值、内变量值和应变增量计算本次试算的应力增量，融合应力增量和融合内变量增量分别是前面每次应力试算得到的应力增量和内变量增量按照影响比例缩小后的加权融合值。

步骤四：将初始土体状态数据中的应变值与影响应变增量相加，得到调整后的土体状态数据，并利用调整后的土体状态数据和硬化准则计算得到本次试算的内变量增量，影响应变增量是按照影响比例缩小后的应变增量。其中，影响比例从第二次试算开始逐次递增，直至第八次试算递增为1。

具体地，考虑到边界面模型中硬化规律通常与应变相关，硬化规律矩阵受应变增量大小影响，本实施例在迭代试算过程中，还提出了逐步增加应变增量对硬化规律矩阵的影响的方案，从而进一步提高了八次试算的准确度，减小局部截断误差，提高积分收敛度，从而提高最终预测的土体状态数据的准确度。具体地方案原理解释为：

在本实施例中，提出了参数“影响比例”，该“影响比例”用于逐步增加应变增量对硬化规律矩阵的影响。如果本次试算是第一次应力试算，由于本次试算距离上一子增量步最近，此时硬化规律矩阵Δλ还不考虑应变增量的影响，第一次试算与现有技术相同，均是先基于初始土体状态数据中的应力值、内变量值和应变增量计算得到应力增量Δσ₁，基于初始土体状态数据和硬化准则计算得到内变量增量Δκ₁，计算公式如下：

第一次试算的应力增量Δσ₁可表示为

式中，σ^(k) _n,κ^(k) _n分别是第n个增量步的第k个子增量步的初始应力数据和初始内变量数据，

是第一次应力试算需要的刚度矩阵，基于边界面塑性模型得到，同时根据硬化规律矩阵Δλ(即硬化准则)计算，得到内变量增量Δκ₁的表达式如下：

如果本次试算不是第一次试算，那么内变量增量与应力增量的计算就需要逐步引入应变增量的影响。为了表达的更加清楚，本实施例从公式的角度对上述步骤三进行说明，假设本次试算是第三次试算，则将应变增量Δε^(k)逐步引入硬化准则计算公式如下：

式中，

分别是第二次试算结束时，作为第三次试算初始状态的应力值和内变量值，“影响比例*Δε^(k)”即为影响应变增量，影响比例随着试算次数的增加而递增，例如第三次试算的影响比例大于第二次试算的影响比例，直至第八次试算，影响比例的值递增为1。

相应地，每次试算使用的应力值与内变量值中也增加前一次试算获得增量，以保证应变增量与应力增量、内变量增量同步增加。因此，

的计算方法如下式：

式中，“权重*影响比例*Δσ₁+权重*影响比例*Δσ₂”即为第三试算时的融合应力增量，“权重*影响比例*Δκ₁+权重*影响比例*Δκ₂”即为第三次试算时的融合内变量增量。

第二次试算至第八次试算与上述第三次试算的原理相同，在此不再赘述。

此外，在本实施例中，各次试算对应的融合权重随试算次数呈递增顺序。例如：第三次试算公式为

其中，权重2＞权重1，以使每次试算更新土体状态数据时，更加考虑最近一次试算的影响，提高土体状态数据更新的准确率。

具体地，在本实施例中，从第二次试算至第八次试算，对应的影响比例依次设定为

1、1，从而在硬化准则计算中使应变增量的影响更均匀地施加。

基于上述步骤，为了更加清楚的描述，通过初始土体状态数据进行八次应力试算，进而更新第k个子增量步结束时的土体状态数据的具体计算流程如下：

(1)在第n增量步的开始时，土体状态数据由上一增量步结束时的应力数据、应变数据和内变量数据表示，其中第n个增量步的第k个子增量步的土体状态数据则由第n个增量步的上一个子增量步的土体状态数据表示，记为

此处上、下标括号中的数字分别代表子增量步数与增量步数；

(2)给定第n子增量步的应变增量Δε_(n)，定义T(0≤T≤1)为子增量步计算的“伪时间”，该伪时间只用于表征将增量步划分为子增量步的长度，并不指代真实时间，例如在第一个子增量步开始前，T＝0，△T⁽¹⁾＝1，每将增量步划分出一个子增量步，T均会累加，直至增加到1为止。

(3)计算子增量步的应变增量步长，对于第k个子增量步，应变增量步长Δε^(k)＝ΔT^(k)Δε_(n)，即初始状态△T⁽¹⁾＝1，Δε^(k)＝Δε_(n)。这样做的目的是为了提高计算效率，第一次计算时子增量步的步长即完整的增量步的长度，如果计算不收敛再将子增量步划小；

(4)基于边界面塑性模型计算第一次应力试算所需的刚度矩阵

然后进行第一次应力试算，得到第一应力增量表示为

同时根据硬化规律矩阵Δλ进行计算，得到第一内变量增量Δκ₁的表达式如下：

(5)从第二次试算开始，逐步增加应变增量对硬化规律矩阵的影响，此次计入影响比例为1/6，相应地，应力值与内变量值中也增加前一次试算获得增量的1/6，以保证应变增量与应力增量、内变量增量同步增加。因此，以应力值

与内变量值/>

计算第二次试算所需的刚度矩阵/>

并进行第二应力试算，则第二次试算的第二应力增量与第二内变量增量可由以下公式进行计算；

(6)第三次试算计入影响比例2/9，相应地，应力值与内变量值中也增加前两次试算获得增量的2/9，并以1:2的加权比例进行计算，以保证应变增量与应力增量、内变量增量同步增加，以应力值

与内变量值/>

计算第三次试算所需的刚度矩阵/>

并进行第三应力试算，则第三次试算的第三应力增量与第三内变量增量可由以下公式进行计算；

(7)第四次试算计入影响比例3/7，相应地，应力值与内变量值中也增加前三次试算获得增量的3/7，并以一定的的加权比例进行计算，以保证应变增量与应力增量、内变量增量同步增加，以应力值

与内变量值/>

计算第四次试算所需的刚度矩阵/>

并进行第四次应力试算，则第四次试算的第四应力增量和第四内变量增量可由以下公式进行计算；

(8)第五次试算计入影响比例2/3，以应力值

与内变量值/>

计算第五次试算所需的刚度矩阵

并进行第五次应力试算，第五应力增量与第五内变量增量计算公式如下；

(9)第六次试算计入影响比例3/4，同理，以应力值

与内变量值

计算第六次试算所需的刚度矩阵/>

并进行第六次应力试算，则第六次试算的第六应力增量与第六内变量增量可由以下公式进行计算；

(10)第七次试算计入影响比例1，相应地，应力值与内变量值中也增加前六次试算获得增量的1倍，并以一定的加权比例进行计算，以保证应变增量与应力增量、内变量增量同步增加，以应力值

与内变量值/>

计算第七次试算所需的刚度矩阵/>

并进行第七试算，则第七次试算的第七应力增量与第七内变量增量可由以下公式进行计算；

(11)第八次试算同样计入影响比例1，相应地，应力值与内变量值中也增加前七次试算获得增量的1倍，并以一定的加权比例进行计算，以保证应变增量与应力增量、内变量增量同步增加，以应力值

与内变量值/>

计算第八次试算所需的刚度矩阵/>

并进行第八试算，则第八次试算的应力增量Δσ₈与内变量增量Δκ₈可由以下公式进行计算；

具体地，在一实施例中，上述步骤S104，具体包括如下步骤：

步骤五：通过下式计算局部截断误差，

式中，R^(k)是第k个子增量步的局部截断误差，

是第n个增量步中第k个子增量步的结束土体状态数据中的应力值，/>

是第n个增量步中第k个子增量步的结束土体状态数据中的内变量值，Δσ₁～Δσ₈分别是八组土体增量数据中的应力增量、Δκ₁～Δκ₈分别是八组土体增量数据中的内变量增量。

具体地，本实施例计算局部截断误差时，与各个应力增量和内变量增量相乘的融合系数基于预测的下一子增量步的预测应力值对当前子增量步初始土体状态数据中的应变值和应力值进行泰勒展开，并对比同类项系数得到。从而使局部截断误差中，各个增量的加权结果更加贴近真实的增量情况，进一步提高了局部截断误差计算的准确率。

具体地，在一实施例中，当基于上述八次试算得到八个应力增量和八个内变量增量后，先通过步骤五计算得到局部截断误差，然后进行误差判别流程。若局部截断误差小于容许误差STOL(预设阈值)，则说明当前子增量步的计算收敛，可进入下一步骤，若局部截断误差不小于容许误差STOL，则当前子增量步的计算不收敛，需按下式缩小当前子增量步长△T^(k)，即将当前增量步划分的更加精细，同时返回第一次应力试算，重新计算该子增量步。

ΔT^(k)＝max{ωΔT^(k)，ΔT_min}

式中，ω是缩小比例，R^(k)是局部截断误差，ΔT_min是最小子增量步长。

若局部截断误差小于容许误差STOL，则当前子增量步计算收敛，不需要重新计算，使T＝T+△T，按下式更新当前子增量步结束时的土体状态数据，包括应力状态数据与内变量状态数据，然后计算下一子增量步。

在本实施例中，上式中的

和

即为目标土体增量数据。本实施例更新土体状态数据时只采用了前七次的应力增量和内变量增量，第八次增量数据只用于局部截断误差的判别，是否将第八次增量数据也添加到上式中可基于实际情况灵活选择，本实施例并不以此为限。

此外，在本实施例中，上式中加权运算的权重值，同样通过下一子增量步的预测应力值对当前子增量步初始土体状态数据中的应变值和应力值进行泰勒展开，并对比同类项系数确定，上述权重值组成权重集合，同种增量的权重值的和约等于1，即各个应力增量对应的权重值求和约为1，内变量增量的权重值同理。

在每个子增量步计算完成时，如果T＜1，则继续计算下一子增量步，并且根据当前子增量步的局部截断误差增大下一子增量步长△T^(k+1)，以提高计算效率，子增量步长的放大比例ω按下式计算

为使下一子增量步结束后T≤1，△T^(k+1)还需满足

T^(k+1)＝min{ωT^(k)，1-T^(k)}

然后进行下一子增量步的计算。直至当前增量步的步长被多个子增量步划分完毕。以最后一个子增量步的结束土体状态数据作为当前增量步的结束土体状态数据保存，此时第n增量步结束。

通过上述步骤，本申请提供的技术方案，针对当前子增量步，基于初始土体状态数据进行八次应力试算得到八个应力增量和八个内变量增量的增量数据，然后对八组数据进行加权得到目标土体增量数据，相比两次应力试算取平均的方法，增量数据准确度更高。然后将初始土体状态数据与目标土体增量数据相加，得到当前子增量步结束时的结束土体状态数据。然后基于结束土体状态数据和八组土体增量数据计算局部截断误差，使得局部截断误差更小，更容易满足误差判别条件，如果局部截断误差小于预设阈值则保留当前子增量步的结束土体状态数据，如果局部截断误差大于预设阈值则缩小当前子增量步的步长，重新计算当前子增量步的结束土体状态数据。从而通过几轮试算，进一步提高了当前子增量步的步长精度。然后按照上述方案对每一子增量步进行迭代，得到每一个子增量步结束时的结束土体状态数据，进而可以递推出当前增量步更加准确的结束土体状态数据。同理按照上述计算步骤，提高了各个增量步预测的土体状态数据的准确度，进而提高应力应变全量关系的计算准确度。

此外，由于边界面模型中硬化规律通常与应变相关，硬化规律矩阵受应变增量大小影响，因此在迭代试算过程中，本发明实施例逐步增加应变增量对硬化规律矩阵的影响，使得每次试算的增量数据更加准确，进而提高每个子增量步计算结束土体状态数据的准确率，在进行局部截断误差判别时，更容易收敛，从而最终进一步提高了创建应力应变全量关系的准确度。

如图4所示，本实施例还提供了一种基于边界面模型的土体状态数据预测装置，装置包括：

初始化模块101，用于对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据。详细内容参见上述方法实施例中步骤S101的相关描述，在此不再进行赘述。

试算模块102，用于基于当前子增量步的步长和初始土体状态数据进行八次应力试算，得到八组土体增量数据，土体增量数据中包括应力增量和内变量增量。详细内容参见上述方法实施例中步骤S102的相关描述，在此不再进行赘述。

土体状态数据预测模块103，用于计算初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，目标土体增量数据是基于八组土体增量数据进行加权求和的结果。详细内容参见上述方法实施例中步骤S103的相关描述，在此不再进行赘述。

误差分析模块104，用于基于结束土体状态数据和八组土体增量数据计算局部截断误差，并判断局部截断误差是否小于预设阈值。详细内容参见上述方法实施例中步骤S104的相关描述，在此不再进行赘述。

步长修正模块105，用于若局部截断误差不小于预设阈值，则缩小当前子增量步的步长，并返回基于当前子增量步的步长和初始土体状态数据进行八次应力试算的步骤。详细内容参见上述方法实施例中步骤S105的相关描述，在此不再进行赘述。

本发明实施例提供的一种基于边界面模型的土体状态数据预测装置，用于执行上述实施例提供的一种基于边界面模型的土体状态数据预测方法，其实现方式与原理相同，详细内容参见上述方法实施例的相关描述，不再赘述。

通过上述各个组成部分的协同合作，本申请提供的技术方案，针对当前子增量步，基于初始土体状态数据进行八次应力试算得到八个应力增量和八个内变量增量的增量数据，然后对八组数据进行加权得到目标土体增量数据，相比两次应力试算取平均的方法，增量数据准确度更高。然后将初始土体状态数据与目标土体增量数据相加，得到当前子增量步结束时的结束土体状态数据。然后基于结束土体状态数据和八组土体增量数据计算局部截断误差，使得局部截断误差更小，更容易满足误差判别条件，如果局部截断误差小于预设阈值则保留当前子增量步的结束土体状态数据，如果局部截断误差大于预设阈值则缩小当前子增量步的步长，重新计算当前子增量步的结束土体状态数据。从而通过几轮试算，进一步提高了当前子增量步的步长精度。然后按照上述方案对每一子增量步进行迭代，得到每一个子增量步结束时的结束土体状态数据，进而可以递推出当前增量步更加准确的结束土体状态数据。同理按照上述计算步骤，提高了各个增量步预测的土体状态数据的准确度，进而提高后续创建应力应变全量关系的准确度。

此外，由于边界面模型中硬化规律通常与应变相关，硬化规律矩阵受应变增量大小影响，因此在迭代试算过程中，本发明实施例逐步增加应变增量对硬化规律矩阵的影响，使得每次试算的增量数据更加准确，进而提高每个子增量步计算结束土体状态数据的准确率，在进行局部截断误差判别时，更容易收敛，从而最终进一步提高了创建应力应变全量关系的准确度。

图5示出了本发明实施例的一种基于边界面模型的土体状态数据预测设备，该设备包括处理器901和存储器902，可以通过总线或者其他方式连接，图5中以通过总线连接为例。

处理器901可以为中央处理器(Central Processing Unit，CPU)。处理器901还可以为其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等芯片，或者上述各类芯片的组合。

存储器902作为一种非暂态计算机可读存储介质，可用于存储非暂态软件程序、非暂态计算机可执行程序以及模块，如上述方法实施例中的方法所对应的程序指令/模块。处理器901通过运行存储在存储器902中的非暂态软件程序、指令以及模块，从而执行处理器的各种功能应用以及数据处理，即实现上述方法实施例中的方法。

存储器902可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储处理器901所创建的数据等。此外，存储器902可以包括高速随机存取存储器，还可以包括非暂态存储器，例如至少一个磁盘存储器件、闪存器件、或其他非暂态固态存储器件。在一些实施例中，存储器902可选包括相对于处理器901远程设置的存储器，这些远程存储器可以通过网络连接至处理器901。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。

一个或者多个模块存储在存储器902中，当被处理器901执行时，执行上述方法实施例中的方法。

上述土体状态数据预测设备具体细节可以对应参阅上述方法实施例中对应的相关描述和效果进行理解，此处不再赘述。

本领域技术人员可以理解，实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，实现的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)、随机存储记忆体(Random Access Memory，RAM)、快闪存储器(Flash Memory)、硬盘(Hard Disk Drive，缩写：HDD)或固态硬盘(Solid-State Drive，SSD)等；存储介质还可以包括上述种类的存储器的组合。

虽然结合附图描述了本发明的实施例，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下作出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

Claims (9)

1.一种基于边界面模型的土体状态数据预测方法，其特征在于，所述方法包括：

对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据；

基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，所述土体增量数据中包括应力增量和内变量增量；

计算所述初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，所述目标土体增量数据是基于所述八组土体增量数据进行加权求和的结果；

基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，并判断所述局部截断误差是否小于预设阈值；

若所述局部截断误差不小于所述预设阈值，则缩小所述当前子增量步的步长，并返回所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算的步骤，重新计算；

所述基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，包括：

通过下式计算所述局部截断误差，

式中，R^(k)表示第k个子增量步的局部截断误差，

表示第n个增量步中第k个子增量步的结束土体状态数据中的应力值，/>

表示第n个增量步中第k个子增量步的结束土体状态数据中的内变量值，Δσ₁～Δσ₈分别是八组土体增量数据中的应力增量、Δκ₁～Δκ₈分别是八组土体增量数据中的内变量增量，E_σ和E_κ分别是计算中间量。

2.根据权利要求1所述的方法，其特征在于，所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，包括：

基于当前子增量步的步长计算当前子增量步的应变增量；

若当前是第一次应力试算，则基于所述初始土体状态数据中的应力值、内变量值和所述应变增量计算得到应力增量，并基于所述初始土体状态数据与硬化准则计算得到内变量增量；

若当前不是第一次应力试算，则将所述初始土体状态数据中应力值和内变量值分别与融合应力增量和融合内变量增量对应相加，并利用相加后的应力值、内变量值和所述应变增量计算本次试算的应力增量，所述融合应力增量和所述融合内变量增量分别是前面每次应力试算得到的应力增量和内变量增量按照影响比例缩小后的融合值；

将所述初始土体状态数据中的应变值与影响应变增量相加，得到调整后的土体状态数据，并利用所述调整后的土体状态数据和硬化准则计算得到本次试算的内变量增量，所述影响应变增量是按照所述影响比例缩小后的应变增量；

其中，所述影响比例从第二次试算开始逐次递增，直至第八次试算递增为1。

3.根据权利要求2所述的方法，其特征在于，计算所述融合应力增量和所述融合内变量增量的步骤包括：

将前面每次应力试算得到的应力增量和内变量增量按照影响比例缩小；

将每次应力试算缩小后的应力增量加权融合得到所述融合应力增量，并将每次应力试算缩小后的内变量增量加权融合得到融合内变量增量；

其中，各次试算对应的融合权重随试算次数呈递增顺序。

4.根据权利要求2所述的方法，其特征在于，所述影响比例从第二次试算开始逐次递增，直至第八次试算递增为1，包括：

从第二次试算至第八次试算，对应的影响比例为

1、1。

5.根据权利要求1所述的方法，其特征在于，基于所述八组土体增量数据加权求和得到所述目标土体增量数据的步骤包括：

通过下一子增量步的预测应力值对当前子增量步初始土体状态数据中的应变值和应力值进行泰勒展开，并对比同类项系数确定八组土体增量数据的目标权重集合，所述目标权重集合中同种增量的权重值的和为1；

利用所述目标权重集合与所述八组土体增量数据进行加权求和，得到所述目标土体增量数据。

6.根据权利要求1所述的方法，其特征在于，所述方法还包括：

若所述局部截断误差小于所述预设阈值，则将所述结束土体状态数据作为下一子增量步的初始土体状态数据进行迭代计算，直至当前增量步的步长被多个子增量步划分完毕，并将最后一个子增量步的结束土体状态数据作为当前增量步的结束土体状态数据进行保存。

7.一种基于边界面模型的土体状态数据预测装置，其特征在于，所述装置包括：

初始化模块，用于对当前增量步进行子增量步划分，并获取当前子增量步的初始土体状态数据；

试算模块，用于基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算，得到八组土体增量数据，所述土体增量数据中包括应力增量和内变量增量；

土体状态数据预测模块，用于计算所述初始土体状态数据与目标土体增量数据的和，得到当前子增量步结束时的结束土体状态数据，所述目标土体增量数据是基于所述八组土体增量数据进行加权求和的结果；

误差分析模块，用于基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，并判断所述局部截断误差是否小于预设阈值；

步长修正模块，用于若所述局部截断误差不小于所述预设阈值，则缩小所述当前子增量步的步长，并返回所述基于当前子增量步的步长和所述初始土体状态数据进行八次应力试算的步骤，重新计算；

所述基于所述结束土体状态数据和所述八组土体增量数据计算局部截断误差，包括：

通过下式计算所述局部截断误差，

式中，R^(k)表示第k个子增量步的局部截断误差，

表示第n个增量步中第k个子增量步的结束土体状态数据中的应力值，/>

表示第n个增量步中第k个子增量步的结束土体状态数据中的内变量值，Δσ₁～Δσ₈分别是八组土体增量数据中的应力增量、Δκ₁～Δκ₈分别是八组土体增量数据中的内变量增量，E_σ和E_κ分别是计算中间量。

8.一种基于边界面模型的土体状态数据预测设备，其特征在于，包括：

存储器和处理器，所述存储器和所述处理器之间互相通信连接，所述存储器中存储有计算机指令，所述处理器通过执行所述计算机指令，从而执行如权利要求1-6任一项所述的方法。

9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于使所述计算机执行如权利要求1-6任一项所述的方法。

Images