本发明属于土体动力多维本构，涉及一种贯通式土体等效线性化多维本构二次开发方法及装置。
背景技术：
1、土体动力本构用来描述在动载荷作用下土体的力学性能的变化特性，是进行土结构相互作用体系地震响应分析的重要基础。目前，等效线性化方法是实际工程中考虑土体非线性的最常用的方法，现有的等效线性化软件如shake91、deepsoil等均是基于一维分层土模型来进行计算，其计算结果无法满足工程场地的三维有限元分析结果精度。然而，随着城市综合体、建筑群和基础设施群建设持续聚集，高效准确的工程场地三维有限元分析方法和实现方式亟需发展。
2、土体等效线性化本构模型都主要用于描述水平单向地震动作用下土体发生剪切变形时的动力特性，通过多次迭代后得到等效剪切模量与等效阻尼比，进而应用于工程场地的有限元分析。随着深覆土地下结构、大跨度空间结构、异型空间结构的发展及抗震分析研究的深入，水平与竖向地震动共同作用的多维地震动场景成为研究的新方向和新难点。由于竖向地震动作用下土体发生压缩变形，此时仍将多维地震动作用下的土体考虑为各向同性材料显然是不合适的，会导致数值模拟结果产生较大误差。此外，现有工程动力分析软件尚无法实现面向三维有限分析的土体等效线性化本构模型。
3、专利cn112632674a公开了一种土体本构模拟计算方法和系统，包括：通过对黏性土体进行循环单剪切试验，建立土体的本构模型；获取土体每级剪应变和对应的动剪切模量对土动力本构模型进行拟合得到等效线性粘弹性本构模型；基于土体的弹性边界以及土体与变压器机构接触的实验数据对等效线性粘弹性本构模型进行迭代求解得到最优状态变量；基于最优状态变量确定土体的最大动剪切模量。该专利虽然基于土体试验建立了等效线性化本构模型，但是其只考虑了土体的水平剪切变形导致的动剪切模量与阻尼比的衰减，无法反映多维地震动作用下土体真实的动剪切模量和阻尼比的变化情况，进而不能实现多维地震动作用下三维土体的力学行为仿真。
技术实现思路
1、本发明的目的就是为了克服上述现有土体等效线性化方法技术存在的三维模型下、多维地震动作用下分析精度不足等至少一种缺陷而提供一种贯通式土体等效线性化多维本构二次开发方法及装置，本发明拓展了等效线性化土体本构在多维动力响应分析中的应用场景，同时实现了前处理、后处理工作的集成化，有效减少人工操作断点，提高模型分析效率，为寻求分析高精度和分析低时间成本之间的平衡提供了有效方案。
2、本发明的目的可以通过以下技术方案来实现：
3、本发明的技术方案之一在于，提供一种贯通式土体等效线性化多维本构二次开发方法，该方法包括以下步骤：
4、s1、土体动力特性曲线的拟合与使用：
5、采集动剪切模量-剪应变、阻尼比-剪应变数据点并通过python读取，通过python中的非线性回归进行动剪切模量-剪应变、阻尼比-剪应变的曲线拟合，求得拟合参数取值，并对拟合曲线进行优度分析；
6、s2、基于python的abaqus模型输入参数前处理：
7、通过abaqus获取土体模型的单元数量与装配件的命名抬头并输入python程序中，并在python程序中赋值初始的动剪切模量比、阻尼比与剪应变，生成作为参数输入的文件；
8、s3、基于fortune的多维地震动等效线性化本构：
9、通过abaqus读取步骤s2中参数输入的文件，读入动剪切模量比、阻尼比与剪应变，同时对土体参数赋值初始剪切模量，以及圆频率、泊松比与土体类型，通过运行fortune程序计算刚度矩阵、abaqus分析步之间的应力与应变关系，通过运行python脚本获取并输出最大剪应变与土体类型识别符；
10、s4、基于python的等效线性化本构的迭代与收敛判断：
11、通过python将最大剪应变与土体类型识别符整合至一个数据框架内，计算等效剪应变范数，判断是否需要进行下一次的迭代，迭代时重新生成作为参数输入的文件，并重复进行步骤s3和s4的迭代计算过程，直至迭代完成，生成用于地震响应计算的参数输入的文件。
12、动剪切模量-剪应变、阻尼比-剪应变拟合曲线中的拟合参数取值取决于实际土体的动力特性。通常通过动三轴试验或共振柱试验来获取相关数据点，而后对数据点进行非线性拟合。以往采用origin、matlab、excel等软件进行拟合工作，这些软件缺乏与abaqus软件之间的原生接口，使得目前由试验数据点拟合为动剪切模量-剪应变、阻尼比-剪应变曲线，再到运用于abaqus二次开发过程中具有较多操作断点，难以通过运行脚本实现贯通式的数据处理工作。
13、abaqus软件是当下世界上最先进的大型通用有限元软件之一，拥有成熟的非线性动力相互作用问题的处理体系。同时，考虑到使用者的特定分析需求，abaqus软件提供了大量用户子程序接口，用户可以开发出abaqus软件自身并不具备的材料、单元、载荷等，以及干预abaqus软件的内核计算过程，从而实现高效的自定义建模计算。通过umat材料用户子程序接口，可以实现三维模型在多维地震动作用下的土体动力特性演变。
14、进一步地，步骤s1中动剪切模量-剪应变、阻尼比-剪应变数据点通过土体材性试验采集，土体材性试验包括动三轴试验或共振柱试验。
15、进一步地，步骤s1中动剪切模量-剪应变、阻尼比-剪应变拟合表达式如下：
16、
17、
18、其中，gd为动剪切模量，g0为初始剪切模量，γd为剪应变，a、b为拟合参数，d为阻尼比，dmin为基本阻尼比，d0、β为拟合参数。
19、进一步地，步骤s1中对拟合曲线进行优度分析，优度分析选自简化卡方检验、残差平方和、r2检验中的一种或多种。
20、进一步地，步骤s3中将土体假设为横观各向同性材料，刚度矩阵如下：
21、
22、其中，ex为水平弹性模量，vxy为水平面上的泊松比，ez为竖向弹性模量，vxz为竖向面上的泊松比，gx为水平剪切模量，gz为竖向剪切模量。
23、进一步地，步骤s3中fortune程序利用kelvin粘弹性模型来模拟土体的非线性行为，三维情况下的应力-应变关系张量表达式如下：
24、
25、其中，σii为正应力张量，σij为切应力张量，εii为正应变张量，εij为切应变张量，为正应变率张量，为切应变率张量，k为体积模量，nk为体积粘滞系数，g为剪切模量，ng为剪切粘滞系数，εv为体积应变，为体积应变率。
26、进一步地，步骤s4中等效剪应变与等效剪应变范数计算方法如下：
27、
28、
29、其中，为第i次迭代等效剪应变，为第i次迭代剪应变，||γeq||i为第i次迭代等效剪应变范数，i为迭代次数，n为土体单元个数。
30、进一步地，步骤s4中通过将迭代前后等效剪应变范数的相对误差与给定容差相对比，判断是否需要进行下一次的迭代，误差对比计算方法如下：
31、
32、其中，||γeq||i为第i次迭代等效剪应变范数，||γeq||i+1为第i+1次迭代等效剪应变范数，i为迭代次数。
33、作为优选的技术方案，所述给定容差为3-7％。
34、进一步地，步骤s4中若需要进行再一次的迭代计算，则基于步骤s1中的动剪切模量-剪应变、阻尼比-剪应变的拟合曲线，计算地震动作用后各个土体单元的动剪切模量与阻尼比，生成作为参数输入的文件，重复进行步骤s3和s4的迭代计算过程；
35、若无需进行再一次的迭代计算，则基于步骤s1中的动剪切模量-剪应变、阻尼比-剪应变的拟合曲线，计算地震动作用后各个土体单元的动剪切模量与阻尼比，生成用于地震响应计算的参数输入的文件。
36、作为优选的技术方案，所述方法包括以下步骤：
37、s1、土体动力特性曲线的拟合与使用：
38、s1-1、输入土体材性试验或文献查阅所得到的动剪切模量-剪应变、阻尼比-剪应变数据点，作为拟合曲线的初始条件；
39、s1-2、通过python中的非线性回归进行动剪切模量-剪应变、阻尼比-剪应变的曲线拟合，求得拟合参数取值，并对拟合曲线进行优度分析；
40、s2、基于python的abaqus模型输入参数前处理：
41、s2-1、通过abaqus查询并记录土体模型的单元数量；
42、s2-2、通过abaqus查询并记录装配件的命名抬头；
43、s2-3、将单元数量与命名抬头作为变量输入python程序中，并在python程序中赋值初始的动剪切模量比、阻尼比与剪应变，取值为(1，0，0)，生成作为参数输入的文件；
44、s3、基于fortune的多维地震动等效线性化本构：
45、s3-1、通过abaqus读取步骤s2中参数输入的文件，读入动剪切模量比、阻尼比与剪应变，同时对土体参数赋值初始剪切模量，以及圆频率、泊松比与土体类型；
46、s3-2、将土体假设为横观各向同性材料，计算刚度矩阵；
47、s3-3、通过运行fortune程序实现abaqus分析步之间的应力与应变关系计算；
48、s3-4、通过运行python脚本获取并输出最大剪应变与土体类型识别符；
49、s4、基于python的等效线性化本构的迭代与收敛判断：
50、s4-1、通过python将最大剪应变与土体类型识别符整合至一个数据框架内；
51、s4-2、基于读入的最大剪应变计算等效剪应变，进一步计算等效剪应变范数；
52、s4-3、通过将迭代前后等效剪应变范数的相对误差与给定容差相对比，判断是否需要进行下一次的迭代；
53、s4-4、若需要进行再一次的迭代计算，则输出继续迭代提示，并基于步骤s1中的动剪切模量-剪应变、阻尼比-剪应变的拟合曲线，计算地震动作用后各个土体单元的动剪切模量与阻尼比，生成作为参数输入的文件，重复进行步骤s3和s4的迭代计算过程；
54、s4-5、若无需进行再一次的迭代计算，则输出无需迭代提示，并基于步骤s1中的动剪切模量-剪应变、阻尼比-剪应变的拟合曲线，计算地震动作用后各个土体单元的动剪切模量与阻尼比，生成用于地震响应计算的参数输入的文件。
55、本发明的技术方案之一在于，提供一种贯通式土体等效线性化多维本构二次开发装置，该装置实现所述的方法，所述装置包括曲线拟合模块、参数前处理模块、二次开发模块、迭代收敛模块和输入输出模块，所述曲线拟合模块对土体动力特性曲线进行拟合，所述参数前处理模块对土体模型的单元数量与装配件的命名抬头进行参数前处理，所述二次开发模块对土体单元的最大剪应变进行计算，所述迭代计算模块对土体单元的等效剪应变范数进行计算并判断是否收敛需要迭代，所述输入输出模块对采集的数据，以及所述方法本身产生的数据进行输入输出。
56、本发明的技术方案之一在于，提供一种贯通式土体等效线性化多维本构二次开发设备，包括存储器、处理器，以及存储于存储器上的计算机程序，所述处理器执行所述计算机程序以实现所述的方法的步骤。
57、本发明的技术方案之一在于，提供一种计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现所述的方法的步骤。
58、本发明的技术方案之一在于，提供一种计算机程序产品，包括计算机程序，该计算机程序被处理器执行时实现所述的方法的步骤。
59、与现有技术相比，本发明具有以下有益效果：
60、(1)本发明联合fortune、python语言和abaqus软件平台，搭建了满足三维模型动力响应分析的土体等效线性化动力本构力学矩阵和实现途径，填补了现有土体等效线性化建模在三维动力响应分析的空白，克服了土体等效线性化方法在三维动力响应分析过程中跨软件平台调用、计算精度和计算时间复杂度等方面的不足；基于fortune语言，通过abaqus软件的umat子程序接口，实现了其三维模型构建和有限元模拟，满足了多维地震动作用下各计算单元等效剪应变的计算与输出；
61、(2)本发明在python软件中实现动三轴试验、共振柱试验等相关土体材性试验数据点的曲线参数拟合及优度分析，进一步地将该功能集成至基于python的等效线性化迭代计算之中，提高了方法的连贯性与效率；
62、(3)本发明利用python软件编写脚本，实现了提交abaqus计算、abaqus后处理、等效线性化迭代计算与收敛判定、生成迭代后的输入参数文件、再提交abaqus计算的完整循环体系，实现了迭代计算与有限元分析的贯通式流程；经过程序封装后可降低人员操作与学习成本，增强了程序的易用性；
63、(4)本发明的土体多维本构是基于等效线性化方法建立的，通过将土体考虑为横观各向同性材料，实现了水平向与竖向的等效线性化计算，克服了现有等效线性化软件无法考虑竖向地震动作用下等效线性化分析的缺陷，能够反映多维地震动作用下的土体等效线性化性能变化，为合理描述土体的动力特性提供了参考，同时对指导工程实际问题具有重要的科学价值。
1.一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，该方法包括以下步骤：
2.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s1中动剪切模量-剪应变、阻尼比-剪应变数据点通过土体材性试验采集，土体材性试验包括动三轴试验或共振柱试验。
3.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s1中动剪切模量-剪应变、阻尼比-剪应变拟合表达式如下：
4.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s1中对拟合曲线进行优度分析，优度分析选自简化卡方检验、残差平方和、r2检验中的一种或多种。
5.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s3中将土体假设为横观各向同性材料，刚度矩阵如下：
6.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s3中fortune程序利用kelvin粘弹性模型来模拟土体的非线性行为，三维情况下的应力-应变关系张量表达式如下：
7.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s4中等效剪应变与等效剪应变范数计算方法如下：
8.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s4中通过将迭代前后等效剪应变范数的相对误差与给定容差相对比，判断是否需要进行下一次的迭代，误差对比计算方法如下：
9.根据权利要求1所述的一种贯通式土体等效线性化多维本构二次开发方法，其特征在于，步骤s4中若需要进行再一次的迭代计算，则基于步骤s1中的动剪切模量-剪应变、阻尼比-剪应变的拟合曲线，计算地震动作用后各个土体单元的动剪切模量与阻尼比，生成作为参数输入的文件，重复进行步骤s3和s4的迭代计算过程；
10.一种贯通式土体等效线性化多维本构二次开发装置，其特征在于，该装置实现如权利要求1至9中任一项所述的方法，所述装置包括曲线拟合模块、参数前处理模块、二次开发模块、迭代收敛模块和输入输出模块，所述曲线拟合模块对土体动力特性曲线进行拟合，所述参数前处理模块对土体模型的单元数量与装配件的命名抬头进行参数前处理，所述二次开发模块对土体单元的最大剪应变进行计算，所述迭代计算模块对土体单元的等效剪应变范数进行计算并判断是否收敛需要迭代，所述输入输出模块对采集的数据，以及所述方法本身产生的数据进行输入输出。
