一种土体本构模拟计算方法和系统与流程

本发明属于电网防灾减灾技术领域，具体涉及一种土体本构模拟计算方法和系统。

背景技术：

在传统的结构动力计算分析中，对于建筑于由土或岩石所构成的场地之上的结构物，往往假定场地地基为刚性，即只考虑结构物与地基之间力的相互作用问题，而忽略其变形协调问题。而实际上，由于地基材料的非绝对刚性，甚至相比于结构材料来说其具有强烈的非线性变形性质，在动力荷载作用过程中，结构物与地基之间既存在着力的相互作用，也存在着变形的相互约制，进而引起振动能量的相互传播和交换，使得实际结构物的动力反应与在刚性地基假设下所算得的反应有很大差别。因此把结构和地基基础作为一个整体的开放体系，而不是单独把结构当作一个封闭的！与地基介质之间没有任何能量交换的系统，来研究其在动力荷载作用下的反应，即构成了土－结构动力相互作用问题。

目前，现有外对电气设备的抗震性能开展了较多的研究，并取得了很多重大成果。但大多数的研究都是将地基基础和上部结构分开考虑，把地基作为绝对刚性，上部结构的底部作为固定端进行研究，而在地震过程中电气设备的地基基础不可避免的会发生变形，对上部结构的振动产生影响，这种对上部结构的自振特性、应力及位移的影响还鲜有研究涉及，导致现有方法不能准确计算土体本构模拟。

技术实现要素：

为克服上述现有技术的不足，本发明提出一种土体本构模拟计算方法，包括：

通过对黏性土体进行循环单剪切试验，建立所述土体的本构模型；

获取土体每级剪应变和对应的动剪切模量对所述土动力本构模型进行拟合得到等效线性粘弹性本构模型；

基于土体的弹性边界以及土体与变压器机构接触的实验数据对所述等效线性粘弹性本构模型进行迭代求解得到最优状态变量；

基于最优状态变量确定所述土体的最大动剪切模量。

优选的，其特征在于，所述通过对黏性土体进行循环单剪切试验，建立所述土体的本构模型，包括：

通过对黏性土体进行循环单剪切试验，得到所述土体的动剪切模量—剪应变关系曲线和阻尼比—剪应变曲线；

对所述剪切模量—剪应变关系曲线进行拟合，得到动剪切模量的计算式，并对所述阻尼比—剪应变曲线进行拟合，得到阻尼比的计算式；

以所述动剪切模量和阻尼比的计算式作为土体的本构模型。

优选的，所述动剪切模量的计算式如下：

式中，α为模量曲线衰减参数；γa为剪应变幅值；γref为参考剪应变幅值；gmax为最大动剪切模量；g为当动应变幅值等于γa时的动剪切模量。

所述阻尼比的计算式如下：

式中d为阻尼比，dmasing表示masing滞回准则预计阻尼比，dmin表示小应变阻尼比，b为阻尼比曲线参数。

优选的，所述masing滞回准则预计阻尼比dmasing的计算式如下：

dmasing＝c1dmasing,α＝1.0+c2dmasing,α＝1.0²+c3dmasing,α＝1.0³

式中，c1为一次系数，c2为二次系数，c3为三次系数，dmasing,α＝1.0是α为1.0时的masing滞回准则预计阻尼比；

dmasing,α＝1.0的计算式如下：

优选的，所述获取土体每级剪应变和对应的动剪切模量对所述土动力本构模型进行拟合得到等效线性粘弹性本构模型，包括：

获取土体每级剪应变和对应的动剪切模量对所述土动力本构模型，在工程模拟软件环境中进行拟合，得到等效线性黏弹性模型。

优选的，所述等效线性黏弹性模型的计算式如下：

式中，gmax为最大剪切模量，k为系数，n为材料参数，pa为标准大气压，σ′3为土体单元围压。

优选的，所述基于土体的弹性边界以及土体与变压器机构接触的实验数据对所述等效线性粘弹性本构模型进行迭代求解得到最优状态变量，包括：

设置土体的弹性边界以及土体与变压器结构的接触；

对状态变量进行初始化，并根据所述状态变量、土体与变压器机构接触的实验数据和等效线性粘弹性本构模型计算土体地震反应的最大剪应变，所述状态变量包括：地震前土体单元的围压、与土体剪应变相关的剪切模量比、阻尼比和土体最大剪应变；

对所述最大剪应变的值进行迭代计算，直到收敛，得到最优状态变量。

优选的，所述基于最优状态变量确定所述土体的最大动剪切模量，包括：

根据所述最优状态变量和等效线性粘弹性本构模型计算得到所述土体的最大动剪切模量。

基于同一发明构思，本申请还提供了一种土体本构模拟计算系统，包括：本构模型模块、粘弹性模型模块、迭代计算模块和最大动剪切模量模块；

所述本构模型模块，用于通过对黏性土体进行循环单剪切试验，建立所述土体的本构模型；

所述粘弹性模型模块，用于获取土体每级剪应变和对应的动剪切模量对所述土动力本构模型进行拟合得到等效线性粘弹性本构模型；

所述迭代计算模块，用于基于土体的弹性边界以及土体与变压器机构接触的实验数据对所述等效线性粘弹性本构模型进行迭代求解得到最优状态变量；

所述最大动剪切模量模块，用于基于最优状态变量确定所述土体的最大动剪切模量。

优选的，所述本构模型模块包括：数据曲线单元、关系式单元和本构模型单元；

所述数据曲线单元，用于通过对黏性土体进行循环单剪切试验，得到所述土体的动剪切模量—剪应变关系曲线和阻尼比—剪应变曲线；

所述关系式单元，用于对所述剪切模量—剪应变关系曲线进行拟合，得到动剪切模量的计算式，并对所述阻尼比—剪应变曲线进行拟合，得到阻尼比的计算式；

所述本构模型单元，用于以所述动剪切模量和阻尼比的计算式作为土体的本构模型。

与最接近的现有技术相比，本发明具有的有益效果如下：

本发明提供了一种土体本构模拟计算方法和系统，包括：通过对黏性土体进行循环单剪切试验，建立土体的本构模型；获取土体每级剪应变和对应的动剪切模量对土动力本构模型进行拟合得到等效线性粘弹性本构模型；基于土体的弹性边界以及土体与变压器机构接触的实验数据对等效线性粘弹性本构模型进行迭代求解得到最优状态变量；基于最优状态变量确定土体的最大动剪切模量；本发明考虑了土体与变压器相互作用的影响，计算模型依托试验结果参数取值，计算可信度较高；考虑了黏性土体的影响，同时适用于三维的土体和结构模型，更贴近工程实际。

附图说明

图1为本发明提供的一种土体本构模拟计算方法流程示意图；

图2为本发明涉及的理想的土体动剪应力-剪应变关系曲线示意图；

图3为本发明提供的扰动土的试验结果曲线示意图；

图4为本发明提供的扰动土动剪切模量曲线示意图；

图5为本发明提供的扰动土darendeli模型阻尼比拟合曲线示意图；

图6为本发明提供的一种土体本构模拟计算系统基本结构示意图；

图7为本发明提供的一种土体本构模拟计算系统详细结构示意图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步的详细说明。

本发明提出了一种考虑土-结构相互作用的土体本构模拟计算方法，其中计算方法中土体与结构均为三维结构，与实际工程接近。通过该计算方法，针对特高压变压器安装于黏性土地基上时，可以较为精确的分析在地震载荷作用下，土体的剪应变，该计算方法考虑到了黏性土地基与上部结构的相互作用影响，与工程实际较为一致。

实施例1：

本发明提供的一种土体本构模拟计算方法流程示意图如图1所示，包括：

步骤1：通过对黏性土体进行循环单剪切试验，建立土体的本构模型；

步骤2：获取土体每级剪应变和对应的动剪切模量对土动力本构模型进行拟合得到等效线性粘弹性本构模型；

步骤3：基于土体的弹性边界以及土体与变压器机构接触的实验数据对等效线性粘弹性本构模型进行迭代求解得到最优状态变量；

步骤4：基于最优状态变量确定土体的最大动剪切模量。

步骤1具体包括：

本发明给出了一种考虑土-结构相互作用的土体本构模拟计算方法，首先通过对黏性土进行循环单剪切试验，得到土体的本构关系，通过对试验数据进行分析，得到适合的本构模型。以一组黏性土循环单剪试验为例，理想的土体动剪应力和剪应变关系曲线如图2所示，扰动土的试验结果如图3所示。

其中动剪切模量g和阻尼比d的定义分别可以表示如式(1)和式(2)

g＝(|τd1|+|τd2|)/(|γd1|+|γd2|)(1)

d＝a0/(πat)(2)

式中τd1、τd2分别为正负最大动剪应力；γd1、γd2分别为正负最大动剪应变，a0为滞回圈的面积，表示了在循环一周中土体消耗的能量，at为三角形abc的面积。

本构模型采用土动力本构模型darendeli模型，其动力特性方程表示为式(3)：

其中α为模量曲线衰减参数；γa为剪应变幅值；γref为参考剪应变幅值；gmax为最大动剪切模量；g为当动应变幅值等于γa时的动剪切模量。

根据试验数据记录表，将每级应变对应的动剪切模量汇成动剪切模量—剪应变散点图，并按照式(3)对动剪模量散点进行拟合，如图4所示。

黏性土的阻尼比采用darendeli模型，其动力特性曲线方程见下式(4)-(9)。

dmasing＝c1dmasing,α＝1.0+c2dmasing,α＝1.0²+c3dmasing,α＝1.0³(5)

c1＝-1.1143α²+1.8618α+0.2523(7)

c2＝0.0805α²-0.0710α-0.0095(8)

c3＝-0.0005α²+0.0002α+0.0003(9)

式中d为阻尼比，dmasing表示masing滞回准则预计阻尼比，dmin表示小应变阻尼比，b为阻尼比曲线参数，dmasing,α＝1.0是α为1.0时的masing滞回准则预计阻尼比。

根据试验数据记录表，将每级应变对应的阻尼比汇成动阻尼比—剪应变散点图，并分别按照式(4)和式(5)-式(9)对阻尼比散点进行拟合，如图5所示。

步骤2具体包括：

接下来需利用工程模拟软件环境例如abaqus进行二次开发来创建土的本构模型，使本构模型与试验数据相吻合。结构的抗震反应分析可以分为等效线性和非线性，土的等效线性模型被绝大多数的学者所采用作为土-结构相互作用抗震分析的土的本构模型，并且取得了很好的计算结果。利用abaqus的二次开发功能来创建等效线性黏弹性模型，如式(10)所示，等效线性黏弹性模型可以较好地用来反映实验中土的非线性动力特性。等效线性模型利用迭代的方法来反映土体中非线性的变化，在每次迭代计算中使用前一步剪应变所对应的剪切模量和阻尼比参与计算，然后得到地震反应的等效剪应变，然后又依据现有等效剪应变对应的动剪切模量和阻尼比参与下一次计算，一般迭代3-5次即可取得较好计算结果。

式中：k为系数，n为材料参数，pa为标准大气压，σ′3为土体单元围压。因此一般需要在abaqus中输入k、n、泊松比μ和圆频率w四个材料参数来表示最大动剪切模量。同时还需要用导入土体单元的四个求解状态变量，状态变量statev(1)～statev(4)分别代表地震前土体单元的围压、与土体剪应变相关的剪切模量比g/gmax、阻尼比d和土体最大剪应变γmax。为了与土动力实验数据相对应，即最大动剪切模量gmax与土动力试验中土的最大剪切模量对应，初始阻尼比为土动力实验最小阻尼比dmin。

步骤3具体包括：

然后建立包含土体和变压器的有限元模型，土体模型采用自定义开发的等效线性粘弹性本构模型，通过设置弹性边界以及土体与变压器结构的接触，每次计算完毕后，读取土体的剪应变作为下一个迭代计算的初始数据，如此迭代3-4次后，土体最大剪应变趋于稳定，即视为结果收敛。其中，状态变量statev(1)～statev(4)的初始化包括：statev(1)初始值为土体循环剪切试验的围压；statev(2)取值为1，定义初始的时候栋剪切模量为最大剪切模量，即g/gmax＝1；statev3)初始值为土体循环单剪试验中测得到的阻尼比；statev(4)取值为0，即初始时最大剪应变为0。

最后进行步骤4的计算，基于最优状态变量确定土体的最大动剪切模量。

实施例2：

基于同一发明构思，本发明还提供了一种土体本构模拟计算系统。

该系统基本结构如图6所示，包括：本构模型模块、粘弹性模型模块、迭代计算模块和最大动剪切模量模块；

本构模型模块，用于通过对黏性土体进行循环单剪切试验，建立土体的本构模型；

粘弹性模型模块，用于获取土体每级剪应变和对应的动剪切模量对土动力本构模型进行拟合得到等效线性粘弹性本构模型；

迭代计算模块，用于基于土体的弹性边界以及土体与变压器机构接触的实验数据对等效线性粘弹性本构模型进行迭代求解得到最优状态变量；

最大动剪切模量模块，用于基于最优状态变量确定土体的最大动剪切模量。

土体本构模拟计算系统详细结构如图7所示。

其中，本构模型模块包括：数据曲线单元、关系式单元和本构模型单元；

数据曲线单元，用于通过对黏性土体进行循环单剪切试验，得到土体的动剪切模量—剪应变关系曲线和阻尼比—剪应变曲线；

关系式单元，用于对剪切模量—剪应变关系曲线进行拟合，得到动剪切模量的计算式，并对阻尼比—剪应变曲线进行拟合，得到阻尼比的计算式；

本构模型单元，用于以动剪切模量和阻尼比的计算式作为土体的本构模型。

其中，迭代计算模块，包括：模型设置单元、初始计算单元和迭代单元；

模型设置单元，用于设置土体的弹性边界以及土体与变压器结构的接触；

初始计算单元，用于对状态变量进行初始化，并根据状态变量、土体与变压器机构接触的实验数据和等效线性粘弹性本构模型计算土体地震反应的最大剪应变，状态变量包括：地震前土体单元的围压、与土体剪应变相关的剪切模量比、阻尼比和土体最大剪应变；

迭代单元，用于对最大剪应变的值进行迭代计算，直到收敛，得到最优状态变量。

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

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

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

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

最后应当说明的是:以上实施例仅用于说明本申请的技术方案而非对其保护范围的限制,尽管参照上述实施例对本申请进行了详细的说明,所属领域的普通技术人员应当理解:本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或者等同替换，但这些变更、修改或者等同替换，均在申请待批的权利要求保护范围之内。