CN110852004A - 力学渗流耦合的非饱和土数值模拟方法及系统 - Google Patents

力学渗流耦合的非饱和土数值模拟方法及系统 Download PDF

Info

Publication number
CN110852004A
CN110852004A CN201910781119.XA CN201910781119A CN110852004A CN 110852004 A CN110852004 A CN 110852004A CN 201910781119 A CN201910781119 A CN 201910781119A CN 110852004 A CN110852004 A CN 110852004A
Authority
CN
China
Prior art keywords
calculation
unsaturated soil
time step
stress
finite element
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201910781119.XA
Other languages
English (en)
Inventor
万柯
冯文凯
廖小芳
王晓兰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujian Geological Engineering Survey Institute
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN201910781119.XA priority Critical patent/CN110852004A/zh
Publication of CN110852004A publication Critical patent/CN110852004A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供了一种力学渗流耦合的非饱和土数值模拟方法及系统，属于非饱和土数值模拟领域。该方法在合理表征非饱和土吸力相关的弹塑性本构行为的巴塞罗那基本模型(BBM)的基础上，发展基于非饱和的Cosserat‑Biot连续体的水力‑力学有限元数值模拟方法。通过采用更加合理地表征非饱和土吸力相关的弹塑性行为的本构模型，及更准确反映土体渗流‑变形耦合行为的连续体模型，提高非饱和土的力学渗流耦合数值模拟的准确性。

Description

力学渗流耦合的非饱和土数值模拟方法及系统
技术领域
本申请涉及非饱和土数值模拟领域，具体而言，涉及一种力学渗流耦合的非饱和土数值模拟方法及系统。
背景技术
非饱和土的变形是由其内部应力和孔隙流体的耦合作用引起的。为了准确合理地描述非饱和土的水力-力学本构行为，其本构模型中采用了由应力、孔隙流体参量共同定义的有效应力状态变量。国际上非饱和土研究领域的知名专家Bishop、Fredlund、Coleman等人在其研究成果中强调了非饱和土的本构模型中引入孔隙水渗流效应相关的吸力、有效应力参量的重要性。在此基础上，Alonso和Gens提出了针对非饱和土的著名的Alonso-Gens本构模型，或称为巴塞罗那基本模型(BBM)，并在此本构模型中将净应力和吸力定义为基本应力状态变量。此后BBM开始广泛地应用于非饱和土的研究和工程领域，此模型也在广泛地应用中得到了更加深入地修正和发展。
在二十世纪初，Cosserat兄弟首先提出引入了“微旋转”自由度的Cosserat连续体，而后Eringen总结并发展了Cosserat连续体具体的理论框架。Cosserat连续体具有在一定程度上考虑土体的细观微结构特征的优势，在经典连续体平移自由度的基础上引入了“微旋转”自由度，并且引入了具有微结构“特征长度”意义的内尺度参数。并且Cosserat连续体对有限元模型引入了正则化机制，具有消除对应变局部化现象数值模拟产生的病态网格依赖性的优势。
非饱和土的数值模型研究是在饱和土数值模型理论的基础上开展的。基于饱和土分析的Biot理论框架，英国的Zienkiewicz率先发展了非饱和土的数值模型，然而在他们的模型中由于应用了被动空气压力假定，而忽略了非饱和土孔隙气体压力的变化和孔隙气体的流动，导致了此模型的误差和应用上的局限性。并且很多国内学者如李锡夔、武文华、吴礼舟、徐炎兵等也进一步发展了固、液、气三相耦合分析的非饱和土有限元数值模拟，使非饱和土的数值模拟在国内的土力学研究和工程领域得到了广泛的应用和长足的进步。
然而，当前已有的非饱和土数值模拟研究对合理地表征土体细观尺度的土颗粒旋转这一细观特征，以及准确地表征表征吸力相关的弹塑性本构行为，并由此开展的非饱和土力学-渗流耦合效应引起的非饱和土性质研究是有所欠缺的，导致非饱和土力学渗流耦合数值模拟的准确性低。
发明内容
有鉴于此，本发明提供了一种力学渗流耦合的非饱和土数值模拟方法及系统，旨在提高非饱和土的力学渗流耦合数值模拟的准确性。
第一方面，本申请提供一种力学渗流耦合的非饱和土数值模拟方法，包括如下步骤：
S100：根据非饱和土的地质模型建立有限元数值模型；
S200：根据非饱和土物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；
S300：根据所述有限元数值模型、参数特征确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；
S400：以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
结合第一方面，在第一方面的第一种可能的实现方式中，所述以均匀时间步长递增时间的方式开展有限元循环计算，获取分析的计算步骤包括
S410：施加当前时间步的力学、孔隙流体的边界荷载增量。
S420：对有限单元计算应变增量、更新单元内的孔隙流体压力、流速。
S430：在单元内实施BBM一致性算法的计算，更新单元应力，计算单元一致性切线模量。
S440：计算单元节点的广义内力以及节点的广义不平衡力。
S450：判断节点的广义不平衡力是否满足平衡要求，如满足平衡要求则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算。
S460：计算单元刚度阵，由单元刚度阵累加获得总体刚度阵，由总体刚度阵和节点的广义不平衡力计算更新节点广义位移。
S470：判断节点的广义位移是否满足收敛条件，如满足则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算。
S480：跳出迭代循环至下一时间步时相应获得当前时间下模型的位移、应力、应变、孔隙流体压力、流速。
S490：判断是否达到总的计算时间步，未达到则开展新的时间步计算，达到则跳出有限元计算并获取相应的分析结果。
结合第一方面的第一种可能的实现方式，在第一方面的第二种可能的实现方式中，BBM一致性算法的计算步骤为：
基于Cosserat连续体理论，对非饱和土的固体骨架的变形计算，其运动微分方程可表示为：
Figure BDA0002176624330000031
式1中ui和ωi分别为固体骨架的位移和微转角，σjk和μjk分别为传统应力和偶应力，bi和gi分别为体力和体力偶，ρ和Ic分别是密度和微惯性矩，lc为内尺度参数，δik和eijk分别为克罗内克尔张量和置换张量，表示梯度算子。
考虑非饱和土中孔隙水和孔隙气体在渗流过程中满足的质量守恒方程和动量守恒方程，孔隙水和孔隙气体满足的控制方程为
Figure BDA0002176624330000042
Figure BDA0002176624330000043
其中,
Figure BDA0002176624330000044
Figure BDA0002176624330000045
Figure BDA0002176624330000046
式2～6中下标w和a分别表示孔隙水和气，Sw、Sa是饱和度，ρw、ρa表示密度，
Figure BDA0002176624330000047
表示达西速度，n为孔隙率。
针对单元内的应力更新计算，考虑在一个时间增量步上的塑性变形，当前应力表示为：
Figure BDA0002176624330000048
式7中Δσe为当前时间步应力增量的弹性预测值，
Figure BDA0002176624330000049
为应力作用引起的塑性应变增量，Δεs为吸力引起的应变增量，De为线弹性本构模量矩阵。
应力增量的弹性预测值表示为
Δσe＝DeΔε -式8
式8中Δε为当前时间步总的应变增量；
应用关联塑性流动准则，塑性应变增量
Figure BDA0002176624330000051
则表示为
Figure BDA0002176624330000052
式9中Δλ为塑性乘子，s为偏应力向量。
基于牛顿-拉弗森方法迭代计算当前时间步的塑性乘子Δλ，考虑在从第v步迭代到第v+1步时，屈服面方程应得到满足，即
由此，当前时间步的塑性乘子表示为
Figure BDA0002176624330000054
Δλv+1＝Δλv+δ(Δλ)v -式12
式10～12中，下标v和v+1分别代表第v步和第v+1步的变量，屈服面偏导数F,λ表示为
Figure BDA0002176624330000055
式13中F,p为加载湿陷屈服面对等效塑性的偏导数。
当前时间步应力向量的计算式为
σv+1＝σv-(2G*s+b2m)vδ(Δλ)v -式14
当前等效塑性应变的更新表示为
Figure BDA0002176624330000056
其中，修正的剪切模量和体积模量分别表示为
Figure BDA0002176624330000061
Figure BDA0002176624330000062
其中，相关参数表示为
Figure BDA0002176624330000064
Figure BDA0002176624330000066
推导得到的应力对应变和吸力偏导量的一致性弹塑性切线模量矩阵
Figure BDA0002176624330000067
和
Figure BDA0002176624330000068
表示如下：
Figure BDA0002176624330000069
Figure BDA00021766243300000610
Figure BDA00021766243300000611
式22～24中，
Figure BDA00021766243300000612
Figure BDA00021766243300000613
结合第一方面的第二种可能的实现方式，在第一方面的第三种可能的实现方式中，有限元数值模型包括有限元节点坐标、网格类型、网格节点编号信息、积分点信息、模型的荷载及约束。
第二方面，本申请提供一种力学-渗流耦合的非饱和土数值模拟系统，包括
建模模块，所述建模模块用于根据非饱和土的地质模型建立有限元数值模型；
参数获取模块，所述参数获取模块用于根据非饱和土物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；
计算模型，所述计算模块用于根据所述有限元数值模型、参数特征确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；
分析模块，所述分析模块用于以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
第三方面，本申请提供一种电子设备，包括存储器和处理器，所述存储器中存储有可在处理器上运行的计算机程序，所述处理器执行所述计算及程序时实现前述的力学渗流耦合的非饱和土数值模拟方法的步骤。
第四方面，本申请提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现前述的的力学渗流耦合的非饱和土数值模拟方法的步骤。
本发明的有益效果是：本发明提供的力学渗流耦合的非饱和土数值模拟方法，该方法在合理表征非饱和土吸力相关的弹塑性本构行为的巴塞罗那基本模型(BBM)的基础上，发展基于非饱和的Cosserat-Biot连续体的水力-力学有限元数值模拟方法。通过采用更加合理地表征非饱和土吸力相关的弹塑性行为的本构模型，及更准确反映土体渗流-变形耦合行为的连续体模型，提高非饱和土的力学渗流耦合数值模拟的准确性。
附图说明
为了更清楚地说明本申请实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1是本申请实施例提供的力学渗流耦合的非饱和土数值模拟方法流程图。
图2是本申请实施例提供的非饱和土力学渗流耦合的有限元方法流程图。
图3是本申请实施例提供的非饱和土三轴压缩试验模型图
图4是本申请实施例提供的三轴压缩试验应力应变曲线图。
图5是本申请实施例提供的压胀破坏模式的土柱的数值模拟等效塑性应变与真实试验破坏形态对比图。
图6是本申请实施例提供的剪切破坏模式的土柱的数值模拟等效塑性应变与真实试验破坏形态对比图。
图7是本申请实施例提供的剪切破坏模式下土柱的孔隙水及孔隙气体压力分布图。
图8是本申请实施例提供的剪切破坏模式下土柱的孔隙水及孔隙气体流线分布图。
图9是本申请实施例提供的不同单元密度的土柱等效塑性应变分布图。
具体实施方式
为使本申请实施方式的目的、技术方案和优点更加清楚，下面将结合本申请实施方式中的附图，对本申请实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本申请一部分实施方式，而不是全部的实施方式。基于本申请中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本申请保护的范围。
因此，以下对在附图中提供的本申请的实施方式的详细描述并非旨在限制要求保护的本申请的范围，而是仅仅表示本申请的选定实施方式。基于本申请中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本申请保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
实施例
基于非饱和土实验的结果分析，Alonso和Gens提出了针对非饱和土的著名的巴塞罗那基本模型(Barcelona Basic Model，简称BBM)，并在此本构模型中将净应力和吸力定义为基本应力状态变量。
请参阅图1和图2，本申请提供一种力学渗流耦合的非饱和土数值模拟方法，包括如下步骤：
S100：根据非饱和土的地质模型建立有限元数值模型。其中，有限元数值模型包括有限元节点坐标、网格类型、网格节点编号信息、积分点信息、模型的荷载及约束。
S200：根据非饱和土物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数。
S300：根据所述有限元数值模型、参数特征确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义。
S400：以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
在一些具体的实施方案中，以均匀时间步长递增时间的方式开展有限元循环计算，获取分析的计算步骤包括
S410：施加当前时间步的力学、孔隙流体的边界荷载增量。
S420：对有限单元计算应变增量、更新单元内的孔隙流体压力、流速。
S430：在单元内实施BBM一致性算法的计算，更新单元应力，计算单元一致性切线模量。
其中，BBM一致性算法的计算步骤为：
基于Cosserat连续体理论，对非饱和土的固体骨架的变形计算，其运动微分方程可表示为：
Figure BDA0002176624330000101
式1中ui和ωi分别为固体骨架的位移和微转角，σjk和μjk分别为传统应力和偶应力，bi和gi分别为体力和体力偶，ρ和Ic分别是密度和微惯性矩，lc为内尺度参数，δik和eijk分别为克罗内克尔张量和置换张量，表示梯度算子；
考虑非饱和土中孔隙水和孔隙气体在渗流过程中满足的质量守恒方程和动量守恒方程，孔隙水和孔隙气体满足的控制方程为
Figure BDA0002176624330000103
其中,
Figure BDA0002176624330000105
Figure BDA0002176624330000106
Figure BDA0002176624330000107
式2～6中下标w和a分别表示孔隙水和气，Sw、Sa是饱和度，ρw、ρa表示密度，
Figure BDA0002176624330000111
表示达西速度，n为孔隙率。
针对单元内的应力更新计算，考虑在一个时间增量步上的塑性变形，当前应力表示为：
Figure BDA0002176624330000112
式7中Δσe为当前时间步应力增量的弹性预测值，
Figure BDA0002176624330000113
为应力作用引起的塑性应变增量，Δεs为吸力引起的应变增量，De为线弹性本构模量矩阵。
应力增量的弹性预测值表示为
Δσe＝DeΔε -式8
式8中Δε为当前时间步总的应变增量；
应用关联塑性流动准则，塑性应变增量
Figure BDA0002176624330000114
则表示为
Figure BDA0002176624330000115
式9中Δλ为塑性乘子，s为偏应力向量；
基于牛顿-拉弗森方法迭代计算当前时间步的塑性乘子Δλ，考虑在从第v步迭代到第v+1步时，屈服面方程应得到满足，即
Figure BDA0002176624330000116
由此，当前时间步的塑性乘子表示为
Δλv+1＝Δλv+δ(Δλ)v -式12
式10～12中，下标v和v+1分别代表第v步和第v+1步的变量，屈服面偏导数F,λ表示为
Figure BDA0002176624330000121
式13中F,p为加载湿陷屈服面对等效塑性的偏导数；
当前时间步应力向量的计算式为
σv+1＝σv-(2G*s+b2m)vδ(Δλ)v -式14
当前等效塑性应变的更新表示为
Figure BDA0002176624330000122
其中，修正的剪切模量和体积模量分别表示为
Figure BDA0002176624330000123
Figure BDA0002176624330000124
其中，相关参数表示为
Figure BDA0002176624330000125
Figure BDA0002176624330000126
Figure BDA0002176624330000127
Figure BDA0002176624330000128
推导得到的应力对应变和吸力偏导量的一致性弹塑性切线模量矩阵
Figure BDA0002176624330000129
和
Figure BDA00021766243300001210
表示如下：
Figure BDA00021766243300001211
Figure BDA00021766243300001212
Figure BDA0002176624330000131
式22～24中，
Figure BDA0002176624330000132
Figure BDA0002176624330000133
S440：计算单元节点的广义内力以及节点的广义不平衡力。
S450：判断节点的广义不平衡力是否满足平衡要求，如满足平衡要求则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算。
S460：计算单元刚度阵，由单元刚度阵累加获得总体刚度阵，由总体刚度阵和节点的广义不平衡力计算更新节点广义位移。
S470：判断节点的广义位移是否满足收敛条件，如满足则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算。
S480：跳出迭代循环至下一时间步时相应获得当前时间下模型的位移、应力、应变、孔隙流体压力、流速。
S490：判断是否达到总的计算时间步，未达到则开展新的时间步计算，达到则跳出有限元计算并获取相应的分析结果。
本发明针对非饱和土吸力相关的弹塑性力学行为和孔隙流体渗流过程的耦合分析，采用了巴塞罗那基本模型(BBM)，发展基于固、液、气三相的Cosserat-Biot连续体模型的有限元方法及其数值过程。针对有限元计算，提出Cosserat-Biot连续体的巴塞罗那基本模型的一致性算法：应力更新的返回映射算法和一致性弹塑性切向模量矩阵的计算。研究成果为全面和合理地分析非饱和土的水力-力学耦合行为提供了参考，通过采用更加合理地表征非饱和土吸力相关的弹塑性行为的本构模型，及更准确反映土体渗流-变形耦合行为的连续体模型，提高非饱和土的力学渗流耦合数值模拟分析的精度和准确性。
在一种具体的实施方案中，为了验证本发明的方法对非饱和土的力学-渗流耦合响应模拟分析的有效性和准确性，考虑一个非饱和土的三轴压缩试验的模拟问题，在平面应变条件下的考虑尺寸为200mm×100mm的矩形土柱，在土柱的上下边界由刚性板施加位移控制的轴向压力σ2＝100kPa作用，对土柱的左右边界施加固定的围压，请参阅图3。
将土样的初始时刻孔隙水饱和度(含水率)设定为均匀分布的，初始含水率分别设定为5.9％、13.3％和17.6％三种情况。将土样的上下边界设置为排水(吸水)边界，而左右边界设定为不排水边界。图4给出了含水率为5.9％、13.3％和17.6％的土样三轴压缩试验的应力应变曲线，以及本文有限元方法模拟的曲线。图4中有限元方法所得的曲线和三轴压缩试验应力-应变曲线的对比，说明有限元计算结果相对于试验结果误差较小。
图5和图6分别给出了压胀破坏模式和剪切破坏模式下有限元计算所得的等效塑性应变分布图表示的破坏现象与试验所得的破坏现象的对比，说明本文的有限元计算能够有效地模拟土体的压胀破坏和剪切破坏现象。图5和图6均说明了本申请提供的有限元计算方法的有效性和精确性。
本发明的数值模拟方法舍弃了被动空气压力假定，能够模拟孔隙气体压力的变化，并且能够模拟由土体变形引起的超孔隙水压力的变化。图7和图8给出的非饱和土柱剪切破坏下的孔隙水、气体的压力分布图及流线图均说明了此点。
相较于经典连续体模型，Cosserat连续体模型考虑了土体细观尺度下微结构的转动效应，引入了正则化机制。本发明所发展的非饱和土Cosserat-Biot连续体模型的有限元方法能够有效地消除剪切带宽度依赖于网格尺寸的病态网格依赖性。图9分别给出了采用15*15、30*30、40*40三种网格密度计算得到的正方形土柱的等效塑性应变分布对比。如图9所示采用本申请提供的方法模拟不同单元密度情况下非饱和土的应变局部化现象，剪切带宽度基本不变，说明该方法引入正则化机制的有效性。
本申请提供一种力学渗流耦合的非饱和土数值模拟系统，包括
建模模块，所述建模模块用于根据非饱和土的地质模型建立有限元数值模型；
参数获取模块，所述参数获取模块用于根据非饱和土物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；
计算模型，所述计算模块用于根据所述有限元数值模型、参数特征确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；
分析模块，所述分析模块用于以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
本申请提供一种电子设备，包括存储器和处理器，所述存储器中存储有可在处理器上运行的计算机程序，所述处理器执行所述计算及程序时实现前述的力学-渗流耦合的非饱和土数值模拟方法的步骤。
本申请提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现前述的力学-渗流耦合的非饱和土数值模拟方法的步骤。
在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device，PLD)(例如现场可编程门阵列(Field Programmable GateArray，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(Altera Hardware DescriptionLanguage)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(RubyHardware Description Language)等，目前最普遍使用的是VHDL(Very-High-SpeedIntegrated Circuit Hardware Description Language)与Verilog2。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。
控制器可以按任何适当的方式实现，例如，控制器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application Specific Integrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式，控制器的例子包括但不限于以下微控制器：ARC 625D、Atmel AT91SAM、Microchip PIC18F26K20以及Silicone Labs C8051F320，存储器控制器还可以被实现为存储器的控制逻辑的一部分。
本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
上述实施例阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本申请各个实施例或者实施例的某些部分所述的方法。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本申请的优选实施方式而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (7)

1.一种力学渗流耦合的非饱和土数值模拟方法，其特征在于，包括如下步骤：
S100：根据非饱和土的地质模型建立有限元数值模型；
S200：根据非饱和土的物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；
S300：根据有限元数值模型、参数特征、动静态模拟要求确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；
S400：以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
2.根据权利要求1所述的力学渗流耦合的非饱和土数值模拟方法，其特征在于，以均匀时间步长递增时间的方式开展有限元循环计算，获取分析的计算步骤包括
S410：施加当前时间步的力学、孔隙流体的边界荷载增量；
S420：对有限单元计算应变增量、更新单元内的孔隙流体压力、流速；
S430：在单元内实施BBM一致性算法的计算，更新单元应力，计算单元一致性切线模量；
S440：计算单元节点的广义内力以及节点的广义不平衡力；
S450：判断节点的广义不平衡力是否满足平衡要求，如满足平衡要求则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算；
S460：计算单元刚度阵，由单元刚度阵累加获得总体刚度阵，由总体刚度阵和节点的广义不平衡力计算更新节点广义位移；
S470：判断节点的广义位移是否满足收敛条件，如满足则跳出迭代循环至下一时间步，反之不满足则在当前时间步继续迭代计算；
S480：跳出迭代循环至下一时间步时相应获得当前时间下模型的位移、应力、应变、孔隙流体压力、流速；
S490：判断是否达到总的计算时间步，未达到则开展新的时间步计算，达到则跳出有限元计算并获取相应的分析结果。
3.根据权利要求2所述的力学渗流耦合的非饱和土数值模拟方法，其特征在于，BBM一致性算法的计算步骤为：
基于Cosserat连续体理论，对非饱和土的固体骨架的变形计算，其运动微分方程可表示为：
Figure FDA0002176624320000021
式1中ui和ωi分别为固体骨架的位移和微转角，σjk和μjk分别为传统应力和偶应力，bi和gi分别为体力和体力偶，ρ和Ic分别是密度和微惯性矩，lc为内尺度参数，δik和eijk分别为克罗内克尔张量和置换张量，
Figure FDA0002176624320000022
表示梯度算子；
考虑非饱和土中孔隙水和孔隙气体在渗流过程中满足的质量守恒方程和动量守恒方程，孔隙水和孔隙气体满足的控制方程为
Figure FDA0002176624320000024
其中,
Figure FDA0002176624320000026
Figure FDA0002176624320000027
式2～6中，下标w和a分别表示孔隙水和气，Sw、Sa是饱和度，ρw、ρa表示密度，表示达西速度，n为孔隙率；
针对单元内的应力更新计算，考虑在一个时间增量步上的塑性变形，当前应力表示为：
Figure FDA0002176624320000032
式7中，Δσe为当前时间步应力增量的弹性预测值，
Figure FDA0002176624320000033
为应力作用引起的塑性应变增量，Δεs为吸力引起的应变增量，De为线弹性本构模量矩阵；
应力增量的弹性预测值表示为
Δσe＝DeΔε -式8
式8中，Δε为当前时间步总的应变增量；
应用关联塑性流动准则，塑性应变增量
Figure FDA0002176624320000034
则表示为
Figure FDA0002176624320000035
式9中Δλ为塑性乘子，s为偏应力向量；
基于牛顿-拉弗森方法迭代计算当前时间步的塑性乘子Δλ，考虑在从第v步迭代到第v+1步时，屈服面方程应得到满足，即
Figure FDA0002176624320000036
由此，当前时间步的塑性乘子表示为
Figure FDA0002176624320000037
Δλv+1＝Δλv+δ(Δλ)v -式12
式10～12中，下标v和v+1分别代表第v步和第v+1步的变量，屈服面偏导数F,λ表示为
Figure FDA0002176624320000041
式13中，F,p为加载湿陷屈服面对等效塑性的偏导数；
当前时间步应力向量的计算式为
σv+1＝σv-(2G*s+b2m)|vδ(Δλ)v -式14
当前等效塑性应变的更新表示为
Figure FDA0002176624320000042
其中，修正的剪切模量和体积模量分别表示为
Figure FDA0002176624320000043
Figure FDA0002176624320000044
其中，相关参数表示为
Figure FDA0002176624320000045
Figure FDA0002176624320000046
推导得到的应力对应变和吸力偏导量的一致性弹塑性切线模量矩阵
Figure FDA0002176624320000049
和
Figure FDA00021766243200000410
表示如下：
Figure FDA00021766243200000411
Figure FDA00021766243200000412
Figure FDA0002176624320000051
式22～24中，
Figure FDA0002176624320000052
Figure FDA0002176624320000053
4.根据权利要求3所述的力学渗流耦合的非饱和土数值模拟方法，其特征在于，有限元数值模型包括有限元节点坐标、网格类型、网格节点编号信息、积分点信息、模型的荷载及约束。
5.一种力学渗流耦合的非饱和土数值模拟系统，其特征在于，包括
建模模块，所述建模模块用于根据非饱和土的地质模型建立有限元数值模型；
参数获取模块，所述参数获取模块用于根据非饱和土物理力学参数、土工经验参数获得非饱和土BBM模型所需的岩石力学和水力学参数；
计算模型，所述计算模块用于根据有限元数值模型、参数特征确定时间步长、荷载增量、计算步数、计算模式、计算结果输出定义；
分析模块，所述分析模块用于以均匀时间步长递增时间的方式开展有限元循环计算，获取分析结果。
6.一种电子设备，其特征在于，包括存储器和处理器，所述存储器中存储有可在处理器上运行的计算机程序，所述处理器执行所述计算及程序时实现上述权利要求1-4任一项所述的力学渗流耦合的非饱和土数值模拟方法的步骤。
7.一种计算机可读存储介质，其特征在于，其上存储有计算机程序，计算机程序被处理器执行时实现权利要求1-4任一项所述的力学渗流耦合的非饱和土数值模拟方法的步骤。
CN201910781119.XA 2019-08-23 2019-08-23 力学渗流耦合的非饱和土数值模拟方法及系统 Pending CN110852004A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910781119.XA CN110852004A (zh) 2019-08-23 2019-08-23 力学渗流耦合的非饱和土数值模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910781119.XA CN110852004A (zh) 2019-08-23 2019-08-23 力学渗流耦合的非饱和土数值模拟方法及系统

Publications (1)

Publication Number Publication Date
CN110852004A true CN110852004A (zh) 2020-02-28

Family

ID=69595542

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910781119.XA Pending CN110852004A (zh) 2019-08-23 2019-08-23 力学渗流耦合的非饱和土数值模拟方法及系统

Country Status (1)

Country Link
CN (1) CN110852004A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949065A (zh) * 2021-03-04 2021-06-11 长江水利委员会长江科学院 模拟层状岩体力学行为的双尺度方法、装置、存储介质及设备
CN113916746A (zh) * 2021-10-14 2022-01-11 深圳大学 一种基于元胞自动机的非饱和地层入渗分析方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014002977A1 (ja) * 2012-06-25 2014-01-03 国立大学法人名古屋大学 空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラム
WO2015124633A1 (en) * 2014-02-19 2015-08-27 Repsol, S.A. Method implemented in a computer for the numerical simulation of a porous medium
CN107958113A (zh) * 2017-11-23 2018-04-24 国家电网公司 一种非饱和膨胀土地基上杆塔基础稳定性数值分析方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014002977A1 (ja) * 2012-06-25 2014-01-03 国立大学法人名古屋大学 空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラム
WO2015124633A1 (en) * 2014-02-19 2015-08-27 Repsol, S.A. Method implemented in a computer for the numerical simulation of a porous medium
CN107958113A (zh) * 2017-11-23 2018-04-24 国家电网公司 一种非饱和膨胀土地基上杆塔基础稳定性数值分析方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
S.P. XIAO等: "A bridging domain method for coupling continua with molecular dynamics", COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING, vol. 193, no. 17, pages 1645 - 1669 *
万柯 等: "Biot-Cosserat 连续体-离散颗粒集合体模型的 非饱和土连接尺度方法", 应用力学学报, vol. 30, no. 3, pages 297 - 303 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949065A (zh) * 2021-03-04 2021-06-11 长江水利委员会长江科学院 模拟层状岩体力学行为的双尺度方法、装置、存储介质及设备
CN113916746A (zh) * 2021-10-14 2022-01-11 深圳大学 一种基于元胞自动机的非饱和地层入渗分析方法
CN113916746B (zh) * 2021-10-14 2023-08-04 深圳大学 一种基于元胞自动机的非饱和地层入渗分析方法

Similar Documents

Publication Publication Date Title
JP5911077B2 (ja) 空気と水と土骨格の連成計算装置および連成計算方法並びに連成計算プログラム
Farhat et al. Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: Momentum and energy conservation, optimal discretization and application to aeroelasticity
Farhat et al. Robust and provably second‐order explicit–explicit and implicit–explicit staggered time‐integrators for highly non‐linear compressible fluid–structure interaction problems
Jeremić et al. Numerical simulation of fully saturated porous materials
Sidarta et al. Constitutive modeling of geomaterials from non-uniform material tests
JP4441693B2 (ja) 水と土骨格の連成計算装置および水と土骨格の連成計算方法
Chandra et al. A robust composite time integration scheme for snap-through problems
Wüchner et al. A framework for stabilized partitioned analysis of thin membrane–wind interaction
Bobby et al. Data-driven performance-based topology optimization of uncertain wind-excited tall buildings
Bojanowski Numerical modeling of large deformations in soil structure interaction problems using FE, EFG, SPH, and MM-ALE formulations
Brun et al. External coupling software based on macro-and micro-time scales for explicit/implicit multi-time-step co-computations in structural dynamics
CN110852004A (zh) 力学渗流耦合的非饱和土数值模拟方法及系统
Jeremić et al. Template elastic‐plastic computations in geomechanics
Fu et al. Wind resistant size optimization of geometrically nonlinear lattice structures using a modified optimality criterion method
Barcelos et al. A Schur–Newton–Krylov solver for steady-state aeroelastic analysis and design sensitivity analysis
Gu et al. Integration of peridynamic theory and OpenSees for solving problems in civil engineering
He et al. New speedup algorithms for nonlinear dynamic time history analysis of supertall building structures under strong earthquakes
Wollny et al. A hierarchical sequential ALE poromechanics model for tire‐soil‐water interaction on fluid‐infiltrated roads
Beran et al. A reduced order cyclic method for computation of limit cycles
Fritzen et al. Space–time model order reduction for nonlinear viscoelastic systems subjected to long-term loading
Kostic et al. An efficient beam-column element for inelastic 3D frame analysis
Rathi et al. Sequential stochastic response surface method using moving least squares-based sparse grid scheme for efficient reliability analysis
Zhuge et al. Research on the multi-shear spring model for circular-section steel piers
Yaghoubi et al. A parallel solution method for structural dynamic response analysis
Allen Grid adaptation for unsteady flow computations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20220221

Address after: 610000 1 East three road, Erxian bridge, Chengdu, Sichuan

Applicant after: Chengdu University of Technology

Applicant after: Fujian Geological Engineering Survey Institute

Address before: Three road 610059 Sichuan city of Chengdu province Chenghua District Erxian Qiaodong No. 1

Applicant before: Chengdu University of Technology