基于超对偶数微分的无约束应力更新算法

路德春; 石安毓; 周鑫; 杜修力

doi:10.11779/CJGE20240138

基于超对偶数微分的无约束应力更新算法 English Version

1.
北京工业大学岩土与地下工程研究所, 北京 100124
2.
清华大学土木工程系, 北京 100084

基金项目:

国家自然科学基金项目 52025084

国家重点研发计划课题项目 2023YFC3009301

中国科协青年人才托举工程项目 2023QNRC001

中国博士后科学基金面上项目 2022M721884

详细信息

作者简介:
路德春(1977—)，男，博士，教授，主要从事岩土与城市地下工程等方面的教学和科研工作。E-mail: dechun@bjut.edu.cn

中图分类号: TU452
计量
- 文章访问数: 353
- HTML全文浏览量: 22
- PDF下载量: 91
出版历程
- 收稿日期: 2024-02-18
- 网络出版日期: 2024-07-25
- 刊出日期: 2025-05-31

An unconstrained stress update algorithm based on hyper-dual step derivative approximation

1.
Institute of Geotechnical and Underground Engineering, Beijing University of Technology, Beijing 100124, China
2.
Department of Civil Engineering, Tsinghua University, Beijing 100084, China

摘要

摘要: 加卸载判断和繁琐的解析求导运算一直是制约先进弹塑性模型数值应用的瓶颈问题。研究提出一种基于超对偶数微分方法的无约束应力隐式更新算法，有效解决了上述计算难点。针对加卸载判断问题，新算法利用光滑函数代替弹塑性本构方程组中的Karush-Kuhn-Tucker条件，将受不等式约束的非线性应力积分方程组问题，转化为无约束的最小化问题，计算时无需加卸载判断。针对导数计算问题，新算法利用超对偶数微分方法代替解析求导，获得光滑函数的1阶导数以及塑性势函数的1阶和2阶导数，用于构造非线性计算的迭代公式，以保证局部应力更新迭代和全局平衡迭代的二次收敛速度。数值算例表明，相较于其它数值微分方法，超对偶数微分方法不受截断误差和减法消去误差影响，计算结果等同于解析求导。最后，基于所提算法编写了光滑莫尔库仑塑性模型的UMAT子程序，并通过3个典型边值问题的数值分析，验证了算法的有效性和收敛性速度。
- 应力更新算法 /
- 塑性模型 /
- 超对偶数数值微分 /
- 一致性切线刚度矩阵 /
- 有限元法
Abstract: The loading/unloading judgment and analytical derivative operations have been the bottlenecks restricting numerical application of elastoplastic models. An unconstrained implicit stress update algorithm is proposed based on the hyper-dual step derivative approximation, which solves the above calculation difficulties. For the problem of loading/unloading judgment, in the new algorithm, the nonlinear stress integral equations with inequality constraints are transformed into an unconstrained minimization problem by using the smooth function to replace the Karush-Kuhn-Tucker conditions. Thus, there is no need for loading/unloading judgement during the calculation. To solve the problem of derivative evaluation, the algorithm uses the hyper-dual step derivative approximation instead of the analytical derivative to obtain the 1st derivative of the smooth function and the 1st and 2nd derivatives of the plastic potential function, which are used to construct iterative formulas for nonlinear calculation, ensuring the quadratic convergence speed of local stress update iterations and global equilibrium iterations. Numerical examples demonstrate that, compared with other numerical differentiation methods, the hyper-dual step derivative approximation is free from truncation errors and subtraction cancellation errors, and its computational results are almost equivalent to the analytical derivations. Finally, based on the proposed algorithm, a UMAT subroutine of smooth Mohr-Coulomb plasticity model is programmed. The effectiveness and convergence speed are verified through numerical analyses of three typical boundary value problems.
- stress update algorithm /
- plastic model /
- hyper-dual number numerical differentiation /
- consistent tangent stiffness matrix /
- finite element method

HTML全文

0. 引言

石油是国家的经济命脉,持续稳定的石油供给是经济与社会可持续发展的重要保障条件^[1]。由于具有安全性能高、成本低、环保等优点,地下水封石油洞库已成为中国战略石油储备的一种重要方式^[2]。地下水封石油洞库的场址应选择在区域稳定性好,岩体质量好,且有稳定地下水的区域^[3]。地下水封石油洞库储备油品品种主要是进口原油,库址宜选择在滨海地区,且距离石油码头或炼油厂距离较近位置处^[4]。滨海地区的地下水封石油洞库建设应该重视海水入侵问题,避免不必要的损失。然而,国内对地下水封石油洞库的海水入侵研究还不充分。

对海水入侵的研究始于19世纪,大部分学者认为海水入侵的本质是滨海地区人为超量开采地下水引起地下水位大幅度下降,海水与淡水之间的动态平衡被破坏,导致咸淡水界面向陆地方向移动^[5]。海水中存在的Cl^-等会对地下工程的结构设施造成一定程度的腐蚀,从而能影响地下工程的长期稳定性和使用寿命。根据海水入侵的机理,为防治海水入侵,可通过控制沿海地区地下水开采,沿海岸带布置人工回灌补给井等针对海水入侵较为有效的防治对策^[6]。

国内外许多学者针对地下水封石油洞库与海水入侵的问题进行了一些研究。Lim等^[7]以韩国地下水封石油洞库为背景,用水文地球化学和同位素指标综合分析滨海裂隙含水层洞库周围海水入侵,得出海水入侵程度的不同是由于研究地点的水文地质异质性造成不同的流入途径。Park等^[8]通过对韩国沿海地区地下水封石油洞库的研究,认为在洞库的建造与运营过程中会产生不同程度的海水入侵现象,且可以采用Cl^-浓度表征海水入侵程度。张彬等^[9]、Li等^[10]研究表明科学设置水幕系统可以有效地降低海水入侵风险。

地下水封石油洞库必须建在一定地下水位线以下以保证洞库周围的地下水压力大于洞内储存介质的压力。为研究地下水封石油洞库的水封可靠性及渗流情况,王者超等^[11-13]依托黄岛地下水封石油洞库工程,分别获得了在不同水幕压力、洞室间液面差、洞室间距和裂隙倾角条件下的地下水头分布和流体渗流情况。在岩体渗流特性中,洞库涌水量规模与空间分布特征是一个至关重要方面^[14]。

本文以某地下水封石油洞库建设项目为背景,在调研了地质勘察资料、现场试验结果的基础上,模拟分析了无水幕系统、水平水幕系统、水平和竖直双水幕系统条件下对Cl^-浓度场、水封性和涌水量影响,为在滨海地区建造地下水封油库面临海水入侵问题的研究提供参考依据。

1. 基本原理

1.1 控制方程

海水入侵的对流-弥散方程由两个部分组成,一部分是地下水流运动方程,另一部分是溶质运移方程。在数值模拟中,采用多物理场耦合方法对两个方程进行耦合。

Huyakorn等^[15]基于质量守恒定律和达西定律,推导出地下水流运动控制方程,表达式如下：

$\begin{matrix} \frac{\partial}{\partial x_{i}} [K_{i j} (\frac{\partial H}{\partial x_{j}} + \frac{ρ_{s} - ρ_{0}}{ρ_{0}} \frac{\partial Z}{\partial x_{j}})] \\ = S_{f} \frac{\partial H}{\partial t} + θ \frac{\partial ρ}{\partial C} \frac{\partial C}{\partial t} - \frac{ρ}{ρ_{0}} q_{s} \end{matrix}$ 。

(1)

式中 $K_{i j}$ 为渗透系数张量（m/s）, $i, j$ =1,2,3； $H$ 为淡水水头（m）； $x_{i}, x_{j}$ 为笛卡尔坐标, $i, j$ =1,2,3； $ρ_{s}$ 为海水密度（kg/m³）； $ρ_{0}$ 为淡水密度（kg/m³）； $ρ$ 为地下水（淡水和海水的合液体）密度（kg/m³）； $Z$ 为计算点高程（m）； $S_{f}$ 为单位储水系数（1/m）； $t$ 为模拟时间（s）； $q_{s}$ 为源（汇）单位体积流量（1/s）； $θ$ 为多孔介质的有效孔隙率； $C$ 为溶质质量浓度（kg/m³）。

溶质质量在多孔介质中通过地下水流动(对流)、分子扩散和机械弥散进行运输。地下水中溶质质量的运移可以用以下对流-弥散方程,表达式如下^[16]：

$\frac{\partial C}{\partial t} + \nabla \cdot (C v) - \nabla \cdot (D \nabla C) = 0$ 。

(2)

式中 $v$ 为流速（m/s）； $D$ 为水动力弥散系数（m²/s）。式（1）,（2）中溶质浓度和淡水水头为待求解函数。

1.2 海水入侵控制原理

对于海水入侵来说,含水层中物质浓度的变化主要是由两种因素引起的：①由于弥散作用,包括分子扩散和机械弥散；②由于液体平均整体运动而引起的物质通量,称为对流^[17]。对地下水封洞库建设而言,初始状态下如图1（a）所示,一方面海水密度大于淡水密度,两者密度差而引起海水向陆地的渗流,促使海水向陆地含水层的入侵,同时海水浓度大于淡水浓度,浓度差也会促使海水向淡水弥散；另一方面,由于陆地含水层的地下淡水保持较高的水头,即淡水位高于海水位,淡水向大海的渗流可以阻止海水向陆地的入侵,两者最终达到一个动态平衡。洞室开挖无水幕系统条件下,开挖洞室会造成地下水位下降,当淡水地下水位降落到一定位置时,淡水水力坡度的减小造成淡水流向海水的对流作用减弱,而淡水与海水之间的弥散作用增强,从而造成海水入侵,如图1（b）所示。有水幕系统条件下,通过水幕系统提高淡水地下水位,增强淡水流向海水的对流作用,减弱淡水与海水之间的弥散作用,从而抑制海水入侵,如图1（c）所示。具体建设中应根据洞库具体条件确定合适的水幕系统,以防止海水入侵。

图 1 洞库海水入侵控制原理示意图

Figure 1. Schematic graph of control principle of seawater intrusion in cavern

下载: 全尺寸图片幻灯片

1.3 海水入侵控制标准

洞室内的海水入侵标准选用《岩土工程勘察规范：GB50021—2001》（2009版）有关地下水对钢筋混凝土中钢筋及钢结构的腐蚀性作为依据^[18]。由于Cl^-的相对原子含量35.5,将规范中mg/L除以相对原子含量35.5转换为mol/m³,如表1。为确保地下水封洞库的使用寿命,洞室内Cl^-最大浓度必须与初始浓度相同,即完全阻碍海水入侵的影响。

表 1 钢筋腐蚀等级转换表

Table 1. Conversion table of steel corrosion grade

腐蚀等级	Cl^-/(mol·m^-3)
微腐蚀	<3
弱腐蚀	3~14.3
中度腐蚀	14.3~143
强腐蚀	>143

下载: 导出CSV

| 显示表格

2. 工程背景

2.1 工程概况

此地下水封石油洞库由相互独立的5条储油洞室组成。储油洞室的顶面标高确定为-62 m,底标高-92 m,洞室埋深约为92～212 m,总长度为2704 m。储油洞室断面为22 m×30 m直墙三心圆拱结构,洞室间距为30.8 m（1.4倍洞跨）。水幕巷道的底面标高确定为-37 m,顶标高-30 m。施工巷道横断面形式采用直墙圆拱结构,宽度为10 m,高度为10 m,施工巷道沿着库区南侧及西侧周边布置。施工巷道入口地面标高为30 m,末端地面标高为-92 m（储油洞室底层）。

2.2 工程地质条件

库区地表大致呈NE20°—30°方向展布,按照丘陵的分布,可分为3个丘陵区,3个丘陵区的高度自北向南递减,它们的最高高程分别为152.6,110.9,81.0 m。角闪辉长岩体是本区最大岩体,面积约2.25 km²,与围岩呈侵入接触关系。岩石呈灰黑色、绿黑色,风化后呈浅灰色。库址区无断裂发育,对库区水文地质影响的断层主要为F1、F8断层。F1带宽30 m,走向NE30°,倾向SE,倾角76°,位于库区东部。断层F8出露于选区南部,呈NW向延伸。

2.3 水文地质条件

根据地质勘查报告,库区除南侧与陆地相连外,西侧、北侧、东侧均为海域,库区周边没有水系发育。根据水文地质测绘得到的水井以及泉水水位,结合钻探得到的孔内水位,综合周边海域极值潮水位,稳定水位值为-3.8 m。地下水接受大气降水的垂直补给,主要以大气降水为补给源。地下水排泄以蒸发和侧向径流排泄方式为主。

根据气象站资料统计,地区的多年平均降水量为977.5～1316.6 mm。海域属强潮海区,其潮汐型态系数为0.26,根据潮汐分类,拟建工程海域的潮汐性质属于正规半日潮。潮汐情况见表2。潮汐波动曲线函数关系为^[10]

$H_{t} = (80 sin \frac{2 π t}{15} + 155) \sin （ 4 π t ） - 20 \sin \frac{2 π t}{365} + 4.73$ 。

(3)

表 2 海区潮汐表

Table 2. Tidal tables in sea area

海区	潮差/m			潮水位/m		最大表层流速/(m·s^-1)
海区	最大	最小	平均	高潮	低潮	涨潮	退潮
海域附近	6.50	3.01	4.73	4.42	-3.80	0.84	0.89

下载: 导出CSV

| 显示表格

3. 海水入侵的影响分析

3.1 模拟方法

图2为库区几何模型示意图,洞室由南到北编号为^#1～^#5,模型长3200 m,宽3000 m,最高点标高156 m,底部标高-250 m。图3为模型与平均海平面的相对位置图,平均海平面位置在0 m。图4（a）为水平水幕系统布置示意图,水平水幕系统底面标高确定为-37 m,顶标高-30 m,四周超出主洞室10 m,孔径0.1 m,间距为10 m。图4（b）为水平和竖直双水幕系统布置示意图,竖直水幕孔上方超出主洞室10 m,下方超出主洞室20 m,孔径0.1 m,间距为10 m,布置在洞室最外围（四周）的水幕巷道中,四周距离洞室50 m。

图 2 库区几何模型示意图

Figure 2. Diagram of geometric model for reservoir area

下载: 全尺寸图片幻灯片

图 3 平均海平面的相对位置

Figure 3. Relative position of mean sea level

下载: 全尺寸图片幻灯片

图 4 水幕系统布置示意图

Figure 4. Schematic diagram of layout of water curtain system

下载: 全尺寸图片幻灯片

在有限元软件中通过多物理场耦合实现渗流场和溶质运移场的耦合,渗流场满足达西定律,溶质运移场主要是溶质对流–弥散作用,即满足式（1）,（2）。库址区岩体以花岗岩为主,洞室段获取钻孔内角闪辉长岩渗透系数平均值为7.5×10^-8 m/s。断层附近钻孔处围岩的渗透性则相对较强,渗透系数平均值可达4.5×10^-7 m/s。扩散系数与岩石（土壤）的孔隙度相关性很强,岩石的扩散系数变化范围很大,查阅相关资料及经验^[19-20],断层扩散系数取1.5×10^-10 m²/s。模型计算参数：花岗岩密度为2600 kg/m³,渗透系数为7.8×10^-8 m/s,扩散系数为1×10^-11 m²/s,储水系数为3.3×10^-6 m,淡水密度为1000 kg/m³,海水密度为1025kg/m³,孔隙率为1.5%。

关于初始条件的设置,在未开挖洞室运行到稳定状态得到的压强和浓度模拟结果作为本次模拟研究结果的初始条件。海水入侵模型分达西渗流场和溶质运移场两个物理场,初始及边界条件也有两套。渗流场中初始地下水压力按自重计算。溶质运移中海水边界根据地下水监测,浓度设置为该地区海水中的浓度值,浓度值大小为500 mol/m³,模型内地下水浓度设置为3 mol/m³。

对于地下水渗流的边界条件,模型的东面、西面、北面三面环海,压力设置为P₁= $ρ_{s}$ g(H_t +z₀-z),z₀为地面高程（m）。模型顶面为降雨入渗边界,由现场气象站统计年降水量,取降水量平均值为1000 mm。考虑到当地的蒸发量,并且查阅资料存在降雨入渗补给系数0.15,综合考虑取值降雨入渗系数为1×10^-9 m/s。模型的底部设置为无流动边界。模型的南面为内陆,根据研究区域的初始地下水深度,压力值为P₂=ρ₀g（z-3.8）。主洞室设置为零压力边界。水平水幕孔压力设置为0.4 MPa,竖直水幕孔压力P₃= $- ρ_{0}$ gh,其中,h为竖直水幕孔高度。对于溶质运移场的边界条件,主洞室开挖后为流出边界,即Cl^-可以随地下水渗流入洞内；模型顶面大气降水中Cl^-浓度为3 mol/m³；模型底部设为无通量边界。水幕系统注入溶液浓度均为3 mol/m³,P₁,P₂,P₃单位均为Pa。

为研究水幕系统对海水入侵过程的影响,首先建立无水幕系统模型,然后建立水平水幕系统模型进行分析,最后在水平水幕的基础上增设竖直水幕系统进行分析,计算工况见表3。

表 3 模拟工况

Table 3. Working coditions of DEM simulations

工况	水平水幕孔间距/m	水平水幕孔压力/MPa	竖直水幕孔间距/m	竖直水幕孔孔底高程超底板深度/m
1	无水幕
2	10	0.4	无竖直水幕
3	10	0.4	10	10

下载: 导出CSV

| 显示表格

3.2 模拟结果

（1）海水入侵控制效果

由图5（a）,（b）,（c）3个工况下各洞室50 a Cl^-浓度云图看出：由于受地质因素、降水因素等诸多条件的影响,Cl^-在地下水的运移过程中不是均匀扩散的,造成洞室东端比西端受海水入侵速度更快,并且^#1～^#5洞室入侵程度依次减小。由图5（a）无水幕系统条件下各洞室Cl^-浓度看出,5个洞室均受海水入侵影响严重。由图5（b）增设水平水幕系统条件下各洞室Cl^-浓度看出,水平水幕系统对于抑制海水入侵具有一定效果,但^#1～^#3洞室东端仍能受到比较严重的海水入侵影响。由图5（c）水平和竖直双水幕系统条件下各洞室Cl^-浓度看出,竖直水幕系统对海水入侵有明显的抑制作用。

图 5 3个工况下各洞室50 a Cl^-浓度云图

Figure 5. Cloud maps of 50-year Cl^- concentration of cavern under three working conditions

下载: 全尺寸图片幻灯片

由图5可知,^#1洞室在5个洞室中受海水影响最严重,取^#1洞室为代表性洞室进行分析。图6为3个工况下,^#1洞室在25,50,75,100 a 4个时间段内Cl^-浓度变化曲线图。由图6（a）,（b）中沿^#1洞室轴线不同时间Cl^-浓度变化曲线图看出,洞室内Cl^-浓度变化趋势基本相同。从入侵时间来看,由于洞室开挖造成地下水位的下降,引起淡水水力坡度减小,淡水流向海水的对流作用减弱,海水与淡水之间的弥散作用逐渐增强,从而引起海水不断入侵,使得洞室在开挖0～100 a之间Cl^-浓度逐渐增加。如图6（c）在设置水平和竖直双水幕系统后,淡水在水力梯度作用下向海水的对流作用增强,海水入侵阻力增大,海水与淡水之间的弥散作用减弱,使得洞室在0～100 a之间达到完全抑制海水入侵的效果。

图 6 3个工况下沿^#1洞室轴线不同时间Cl^-浓度变化曲线图

Figure 6. Variation curve of Cl^- concentration at different time along axis of cavern No. 1 under three working conditions

下载: 全尺寸图片幻灯片

从入侵程度来看（以洞室50 a使用期为节点）,如图6（a）无水幕系统条件下沿^#1洞室轴线Cl^-浓度变化曲线图可知,洞室内最大Cl^-浓度达到约28 mol/m³,洞室沿长轴方向距离约为300 m处于中度腐蚀程度；如图6（b）可知,洞室内最大Cl^-浓度达到约19 mol/m³,洞室沿长轴方向距离约为130 m处于中度腐蚀标准；如图6（c）,洞室内最大Cl^-浓度与初始浓度相同,入侵距离也没有增加,说明在水平和竖直双水幕系统设置下,可以达到完全阻碍海水入侵的作用。

（2）水封性评价

抑制地下水封洞库的海水入侵,应在保证洞室水封性的前提下,即要求水封洞库必须长期处于地下水位以下一定深度。

由图7（a）可知,在不设置水幕系统的条件下,地下水位线会降落到洞室顶部高程以下,洞库上方最低地下水位会下降约至-75 m,距离洞室顶端大约-13 m,根据水封原则,无法保证洞库的水封性,故仅靠天然地下水补给无法达到对储油洞室的水封作用。为确保洞室水封性,首先需要设置水平水幕系统,由图7（b）可知,设置水平水幕系统后,所有洞室上方的地下水位都不会降落到顶拱以下,即洞库上方没有形成与大气连通的非饱和区,洞库上方保留了一定厚度的水盖层。除最外两侧洞室地下水位略有下降,中间洞室地下水位上升至0 m处,可以满足水封性要求。由图7（c）可知,增设竖直水平水幕系统后对洞室地下水位的影响范围主要为洞室两侧区域,两侧区域地下水位提高到0 m,洞室上方基本恢复到初始地下水位,同样满足洞室的水封要求。

图 7 垂直主洞室轴线方向地下水位图（m）

Figure 7. Groundwater level map in the direction perpendicular to the axis of the main cavern

下载: 全尺寸图片幻灯片

（3）洞库涌水量统计

图8为3个工况下各个洞室涌水量直方图,无水幕系统条件下,由于^#1,^#5洞室只有一侧受其他洞室影响,而其他洞室涌水量受左右两侧洞室影响,所以^#1,^#5洞室涌水量比中间洞室大。又^#1洞室受到断层F1、F8影响,断层使涌水断面面积增大,故^#1洞室涌水量最大。

图 8 3个工况下各个洞室涌水量直方图

Figure 8. Histogram of water inflow in each cavern under three working conditions

下载: 全尺寸图片幻灯片

对比3个工况可知,无水幕系统条件下,洞库总涌水量约868.88 m³/d；增设水平水幕系统会大大增加洞室涌水量,在水平水幕系统的作用下,洞室稳定总涌水量共计3045.48 m³/d。在水平和竖直双水幕系统共同作用下,洞室稳定总涌水量共计3776.14 m³/d,竖直水幕系统对洞室总涌水量的贡献不大。

6. 结论

本文以某地下水封油库工程为例,建立了研究区地下水封洞库三维数值模型,分析了无水幕系统、水平水幕系统、水平和竖直双水幕系统条件下Cl浓度场、水封性、涌水量规律,得出下3点结论。

（1）为确保地下水封石油洞库工程不受海水入侵的影响,只设置水平水幕系统达不到完全抑制海水入侵洞室,必须在水平水幕系统的基础上增设竖直水幕系统,才能达到完全抑制海水入侵洞室的效果。

（2）水平水幕系统条件下,除最外两侧洞室地下水位略有下降,中间洞室地下水位上升至0 m处；增设竖直水幕后,会提高洞室两侧区域的地下水位,地下水位可以基本恢复到初始地下水位,两种条件下均能够保证洞室的水封性,但库区无水幕系统条件下,地下水位将不满足水封要求。

（3）在无水幕系统、水平水幕系统、水平和竖直双水幕系统3种工况下,洞室稳定总涌水量分别约为868.88,3045.48,3776.14 m³/d,并且竖直水幕对洞室总涌水量的贡献要远小于水平水幕。

图 1 莫尔库仑屈服函数

Figure 1. Mohr-Coulomb yield function

下载: 全尺寸图片幻灯片

不同 ${\beta _2}$ 值下的平滑曲线

Figure 2. Smooth curves with different values of β₂

下载: 全尺寸图片幻灯片

图 3 算法流程图

Figure 3. Flow chart of algorithm

下载: 全尺寸图片幻灯片

图 4 误差随扰动值h的变化规律

Figure 4. Change of errors with perturbation value h

下载: 全尺寸图片幻灯片

图 5 悬臂式基坑开挖模拟

Figure 5. Simulation of excavation of cantilever foundation pit

下载: 全尺寸图片幻灯片

图 6 墙体位移和内力分布对比

Figure 6. Comparison of wall displacement and internal force distribution

下载: 全尺寸图片幻灯片

图 7 拉剪荷载下的三维圆柱体

Figure 7. 3D cylinder under tensile-shear loads

下载: 全尺寸图片幻灯片

图 8 解析和数值一致性切线刚度矩阵的收敛性对比

Figure 8. Comparison of convergence behavior of analytical and numerical consistent tangent stiffness matrices

下载: 全尺寸图片幻灯片

图 9 条形基础承载力模拟

Figure 9. Simulation of bearing capacity of strip foundation

下载: 全尺寸图片幻灯片

图 10 位移-反力曲线对比

Figure 10. Displacement-reaction force curves

下载: 全尺寸图片幻灯片

表 1 不同方法的微分公式

Table 1 Differentiation formulas for different methods

数值微分方法	${{{\text{d}}f} \mathord{\left/ {\vphantom {{{\text{d}}f} {{\text{d}}x}}} \right. } {{\text{d}}x}} = \left( \cdot \right)$	截断误差	${{{{\text{d}}^2}f} \mathord{\left/ {\vphantom {{{{\text{d}}^2}f} {{\text{d}}{x^2}}}} \right. } {{\text{d}}{x^2}}} = \left( \cdot \right)$	截断误差
FDM	$\frac{{f(x + h) - f(x)}}{h}$	$O(h)$	$\frac{{f(x + 2h) - 2f(x + h) + f(x)}}{{{h^2}}}$	$O(h)$
CDM	$\frac{{f(x + h) - f(x - h)}}{{2h}}$	$O({h^2})$	$\frac{{f(x + h) - 2f(x) + f(x - h)}}{{{h^2}}}$	$O({h^2})$
CSDA	$\frac{{{\Re _i}\left[ {f\left( {x + ih} \right)} \right]}}{h}$	$O({h^2})$	$\frac{2}{{{h^2}}}\left[ {f\left( x \right) - {\Re _i}\left[ {f\left( {x + ih} \right)} \right]} \right]$	$O({h^2})$
HDSDA	$\frac{{{\Im _{{\varepsilon _1}}}\left[ {f\left( {x + {h_{\text{1}}}{\varepsilon _1} + {h_{\text{2}}}{\varepsilon _2}} \right)} \right]}}{{{h_{\text{1}}}}}$	无	$\frac{{{\Im _{{\varepsilon _1}{\varepsilon _2}}}\left[ {f\left( {x + {h_{\text{1}}}{\varepsilon _1} + {h_{\text{2}}}{\varepsilon _2}} \right)} \right]}}{{{h_{\text{1}}}{h_{\text{2}}}}}$	无

下载: 导出CSV

表 2 用于边值问题数值分析的模型参数

Table 2 Model parameters for numerical analysis of boundary value problems

边值问题	$E$ /kPa	$\nu$	${c_0}/{\text{kPa}}$	$\varphi$ /(°)	$\psi$ /(°)	${H_{\text{p}}}$ /kPa
基坑开挖	66000	0.2	0.1	30	0.1	1×10^-8
3D圆柱	30000	0.3	50	12	0.1	1×10^-8
条形基础	30000	0.3	0.1	30	0.1	1×10^-8

下载: 导出CSV

参考文献(24)

[1]	BORJA R I, LEE S R. Cam-clay plasticity, part 1: implicit integration of elasto-plastic constitutive relations[J]. Computer Methods in Applied Mechanics and Engineering, 1990, 78(1): 49-72. doi: 10.1016/0045-7825(90)90152-C
[2]	BORJA R I. Cam-clay plasticity, part Ⅱ: implicit integration of constitutive equation based on a nonlinear elastic stress predictor[J]. Computer Methods in Applied Mechanics and Engineering, 1991, 88(2): 225-240. doi: 10.1016/0045-7825(91)90256-6
[3]	贾善坡, 陈卫忠, 杨建平, 等. 基于修正Mohr-Coulomb准则的弹塑性本构模型及其数值实施[J]. 岩土力学, 2010, 31(7): 2051-2058. JIA Shanpo, CHEN Weizhong, YANG Jianping, et al. An elastoplastic constitutive model based on modified Mohr-Coulomb criterion and its numerical implementation[J]. Rock and Soil Mechanics, 2010, 31(7): 2051-2058. (in Chinese)
[4]	ZHAO J D, SHENG D C, ROUAINIA M, et al. Explicit stress integration of complex soil models[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2005, 29(12): 1209-1229. doi: 10.1002/nag.456
[5]	HONG H K, LIU C S. Internal symmetry in the constitutive model of perfect elastoplasticity[J]. International Journal of Non-Linear Mechanics, 2000, 35(3): 447-466. doi: 10.1016/S0020-7462(99)00030-X
[6]	REZAIEE-PAJAND M, NASIRAI C. On the integration schemes for Drucker-Prager's elastoplastic models based on exponential maps[J]. International Journal for Numerical Methods in Engineering, 2008, 74(5): 799-826. doi: 10.1002/nme.2178
[7]	SLOAN S W, ABBO A J, SHENG D C. Refined explicit integration of elastoplastic models with automatic error control[J]. Engineering Computations, 2001, 18(1/2): 121-194. doi: 10.1108/02644400110365842
[8]	SIMO J C, TAYLOR R L. A return mapping algorithm for plane stress elastoplasticity[J]. International Journal for Numerical Methods in Engineering, 1986, 22(3): 649-670. doi: 10.1002/nme.1620220310
[9]	MIRA P, TONNI L, PASTOR M, et al. A generalized midpoint algorithm for the integration of a generalized plasticity model for sands[J]. International Journal for Numerical Methods in Engineering, 2009, 77(9): 1201-1223. doi: 10.1002/nme.2445
[10]	SCALET G, AURICCHIO F. Computational methods for elastoplasticity: an overview of conventional and less-conventional approaches[J]. Archives of Computational Methods in Engineering, 2018, 25(3): 545-589. doi: 10.1007/s11831-016-9208-x
[11]	SIMO J C, HUGHES T J R. Computational Inelasticity[M]. New York: Springer Verlag, 1998.
[12]	范庆来, 栾茂田, 杨庆. 修正剑桥模型的隐式积分算法在ABAQUS中的数值实施[J]. 岩土力学, 2008, 29(1): 269-273. FAN Qinglai, LUAN Maotian, YANG Qing. Numerical implementation of implicit integration algorithm for modified Cam-clay model in ABAQUS[J]. Rock and Soil Mechanics, 2008, 29(1): 269-273. (in Chinese)
[13]	ARMERO F, PÉREZ-FOGUET A. On the formulation of closest-point projection algorithms in elastoplasticity: Part I The variational structure[J]. International Journal for Numerical Methods in Engineering, 2002, 53(2): 297-329. doi: 10.1002/nme.278
[14]	STARMAN B, HALILOVIČ M, VRH M, et al. Consistent tangent operator for cutting-plane algorithm of elasto-plasticity[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 272: 214-232. doi: 10.1016/j.cma.2013.12.012
[15]	陈洲泉, 陈湘生, 庞小朝, 等. 角点型非共轴模型的半隐式应力积分算法及应用[J]. 岩土工程学报, 2023, 45(3): 521-529. doi: 10.11779/CJGE20220050 CHEN Zhouquan, CHEN Xiangsheng, PANG Xiaochao, et al. Semi-implicit integration algorithm for non-coaxial model based on vertex theory and its application[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(3): 521-529. (in Chinese) doi: 10.11779/CJGE20220050
[16]	KRABBENHOFT K, LYAMIN A V. Computational cam clay plasticity using second-order cone programming[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 209/210/211/212: 239-249.
[17]	ZHOU X, LU D C, SU C C, et al. An unconstrained stress updating algorithm with the line search method for elastoplastic soil models[J]. Computers and Geotechnics, 2022, 143: 104592. doi: 10.1016/j.compgeo.2021.104592
[18]	PEREZ F A. Numerical differentiation for local and global tangent operators in computational plasticity[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 189(1): 277-296. doi: 10.1016/S0045-7825(99)00296-0
[19]	CHOI H, YOON J W. Stress integration-based on finite difference method and its application for anisotropic plasticity and distortional hardening under associated and non-associated flow rules[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 345: 123-160. doi: 10.1016/j.cma.2018.10.031
[20]	SU C C, LU D C, ZHOU X, et al. An implicit stress update algorithm for the plastic nonlocal damage model of concrete[J]. Computer Methods in Applied Mechanics and Engineering, 2023, 414: 116189. doi: 10.1016/j.cma.2023.116189
[21]	ZHANG N, LI X, WANG D. Smoothed classic yield function for C2 continuities in tensile cutoff, compressive cap, and deviatoric sections[J]. International Journal of Geomechanics, 2021, 21(3): 4021005. doi: 10.1061/(ASCE)GM.1943-5622.0001910
[22]	MENETREY P, WILLAM K J. Triaxial failure criterion for concrete and its generalization[J]. ACI Structural Journal, 1995, 92(3): 311-318.
[23]	FIKE J A. Multi-Objective Optimization Using Hyper-Dual Numbers[D]. Palo Alto: Stanford University, 2013.
[24]	ZHOU X, SHI A Y, LU D C, et al. A return mapping algorithm based on the hyper dual step derivative approximation for elastoplastic models[J]. Computer Methods in Applied Mechanics and Engineering. 2023, 417: 116418. doi: 10.1016/j.cma.2023.116418

施引文献(7)

期刊类型引用(6)

1.	王敬奎，乔丽苹，王菲，王者超，李崴. 海岛地下水封石油洞库海水入侵防控. 山东大学学报(工学版). 2025(02): 134-142 . 百度学术
2.	乔丽苹，王菲，王者超，李成. 地下水封石油洞库竖直水幕设计参数及其影响分析. 岩土工程学报. 2024(07): 1525-1533 . 本站查看
3.	魏松源，李涵硕，彭振华，王者超，李崴. 某地下水封洞库竖直水幕试验分析与水幕系统有效性. 隧道与地下工程灾害防治. 2024(04): 81-89 . 百度学术
4.	曲建军，曹彪，宋大钊，杨连枝，何生全. 近海地下石油储备库海水入侵风险研究. 水文地质工程地质. 2023(06): 184-192 . 百度学术
5.	蒋中明，钟兵，万发. 水封石油洞库污染物运移规律研究. 岩土工程学报. 2023(12): 2529-2536 . 本站查看
6.	王秀英. 海滨地区水文环境变化及衍生灾害研究现状分析. 中国石油大学胜利学院学报. 2021(02): 24-28 . 百度学术