Fractal derivative Nishihara rheology model and its application
-
摘要: 软土具有典型的非线性流变特性,传统黏弹塑性元件模型难以高效描述其三阶段流变行为。采用分形阻尼器取代传统阻尼器,考虑软土的损伤演化,提出了改进的能完整描述包含衰减流变、稳定流变及加速流变在内的三阶段黏弹塑性流变特性的分形导数西元流变模型,并给出了解析解。模型具有精度高和计算简捷的特点,与试验结果吻合良好,证明了其合理性。基于ABAQUS平台,编制了分形导数西元黏弹塑性流变模型的UMAT子程序,并通过软土三轴流变试验和路堤工程实例验证了子程序的有效性。结果表明,所建立的模型及开发的子程序能准确地反映软土的完整三阶段非线性流变特性。Abstract: The soft soil exhibits typical nonlinear rheological characteristics, while the traditional element models cannot accurately describe them. By replacing the conventional damper with a fractal damper and considering the damaging effects of soils, an improved nonlinear fractal derivative Nishihara viscoelastic-plastic rheology model which can fully describe the decay stage, steady stage and accelerated stage of rheology is proposed, and the analytical solution is given. The proposed model is characterized by high precision and simple calculation, and is in good agreement with the test ones, which proves its rationality. Based on the ABAQUS platform, the corresponding UMAT subroutine of the fractal derivative Nishihara viscoelastic-plastic rheology model is compiled, and the validity of the subroutine is verified by the triaxial creep tests and an embankment case. The comparison results show that the established model and developed subroutine can accurately reflect the complete three-stage viscoelastic-plastic rheological characteristics of soft soils.
-
Keywords:
- fractal derivative /
- rheological property of soft soil /
- viscoelastic-plasticity /
- ABAQUS /
- UMAT
-
0. 引言
随着交通强国战略的稳健推进,中国高铁、跨海大桥、过江通道及隧道工程建设需求增长迅猛,许多临江滨海和穿山越岭工程不可避免地建于软弱岩土材料之上。实际工程中岩土材料的时效变形乃至失稳破坏都与其流变特性密切相关[1-2],研究岩土材料的流变行为,建立合理的流变模型并将其拓展到数值分析软件中去解决实际问题具有十分重要的理论意义和工程应用价值。
目前用于描述岩土材料流变特性的模型主要有经验模型[3]、理论模型[4]和元件模型[5],元件模型因参数较少且物理意义明确而具有更强的普适性。国内外学者通过传统线性元件的多样化组合对岩土材料的流变特性开展了大量的研究,但线性元件组合模型在描述岩土材料非线性流变特性方面存在着一定的缺陷。为了消除传统元件模型的弊端,研究学者通过将传统元件与经验模型或新型非线性元件相结合等多种方法实现了非线性流变行为的模拟[1, 6],但所建立模型参数较多,工程应用不便;另一些学者通过引入分数阶微积分理论来简化流变模型[7-15],实现了用较少的参数达到较高拟合精度的显著成效,但同时,分数阶导数流变模型的卷积积分算法对计算机内存要求高,存在计算效率低的缺点;近年来,作为一种局部算子的非卷积积分分形导数理论被引入到岩土材料的流变模型求解过程中[16-19],研究发现,分形导数流变模型在提高流变变形计算精度的同时能有效地节省计算成本,值得进一步推广应用。
黏弹塑性元件模型一般由黏弹性元件与塑性元件串联而成,可以实现加速阶段流变变形的表征,目前应用最为广泛的是西元模型[5]。但传统塑性元件属于线性元件,难以准确描述岩土材料在加速流变阶段的非线性流变特性[17]。近年来,随着非线性元件的建立以及损伤理论的发展,国内外学者在岩土工程材料非线性加速流变阶段的研究越来越完善。Lu等[10]建立了能表征土体黏弹塑性流变特性的3D分数阶流变模型。薛东杰等[11]在西元模型的基础之上建立了考虑温度和体积应力的非线性分数阶流变损伤模型。Wu等[12]将传统元件和分数阶阻尼器串联,得到了可以描述盐岩完整流变阶段的非线性分数阶黏弹塑性流变模型。张胜利等[13]将分数阶Merchant模型与黏塑性体损伤模型串联,构建了考虑温度的盐岩分数阶黏弹塑性流变损伤模型。孙增春等[14]基于分数阶微积分理论建立了砂-粉混合料的塑性本构模型。Yin等[19]提出了考虑流变损伤效应的分形导数黏弹塑性流变模型。张亮亮等[20]引入非线性元件和流变损伤,描述了岩体的非线性黏弹塑性变形特性。骆亚生等[21]通过引入新型触发式黏塑性元件,建立了可以描述重塑黄土的蠕变、静力及动力特性的黏弹塑性本构模型。上述研究使得岩土材料的完整三阶段非线性黏弹塑性流变模型的建立日趋完善。
由于土体应力-应变关系的非线性、荷载作用的复杂性及边界条件的多样性,用解析解求解实际工程的流变问题比较困难,而数值模拟则可以实现各种复杂工程问题分析。其中,ABAQUS作为一种大型通用有限元分析软件,凭借其强大的非线性求解能力及前、后处理功能被广泛应用于岩土工程领域。ABAQUS本身带有几种常用的土体本构模型,也可以通过用户自定义UMAT程序编写更多材料的本构模型,因地制宜地解决多种复杂土工分析难题[22-25]。在岩土工程领域,尽管分形导数理论具备通过数学方法简化复杂求解过程的能力,目前的研究中仍很少涉及。而且现有的黏弹塑性模型大多围绕岩石开展,对于土体的适用性仍待考究。此外,如何将改进的新型黏弹塑性流变模型通过二次开发集成于数值模拟软件中去解决复杂工程问题仍有待完善。
为了精确地描述软土材料的完整三阶段非线性黏弹塑性流变变形规律,本文采用分形导数阻尼器来改进传统西元模型,引入损伤因子描述软土加速流变的演化过程,提出一个分形导数黏弹塑性流变模型,并给出模型的解析解。进一步结合ABAQUS二次开发模块,通过FORTRAN编制了分形导数西元黏弹塑性流变模型的UMAT子程序,对临江软土的三轴流变试验及路基工程进行分析,以验证模型及UMAT程序的准确性。
1. 分形导数西元流变模型
1.1 传统流变元件力学特性
岩土材料典型流变曲线如图 1所示。以屈服应力σs为界限,当σ<σs时,流变曲线呈现只包含衰减和稳态流变的Ⅰ型衰减型;当σ⩾σs时,流变曲线呈现包含加速流变的Ⅱ型非衰减型,二者均表现出显著的非线性特性。传统流变模型的线性Hooke弹簧元件及Newton阻尼器元件的本构方程可以统一表述为σ(t)∼dαε(t)/dtβ,其中α,β表示整数型导数阶,这些线性元件均无法准确描述流变曲线的非线性特性。分数阶微分型元件Abel阻尼器的引入实现了岩土材料的非线性流变变形的描述,其本构方程与传统元件表达相同,区别在于导数阶0⩽α=β⩽1。当α=β=0时,Abel阻尼器退化为Hooke弹簧;当α=β=1时,Abel阻尼器则进化为Newton阻尼器,可以很好地表述从理想固体到理想流体之间材料的非线性渐变过程[7]。但分数阶微分型本构关系的存在导致求解过程涉及到卷积积分,为了规避上述问题,本文将分形导数理论引入到流变模型的构建过程中。
1.2 分形导数理论与分形阻尼器
Chen[16]指出分形导数是一种无卷积积分的局部算子,并给出时间分形导数公式:
df(t)dtβ=lim (1) Cai等[17]基于上述分形导数理论,提出了一种新的分形阻尼器元件,其本构方程为
\sigma (t) = \eta \frac{{{{\text{d}}^\alpha }\varepsilon (t)}}{{{\text{d}}{t^\beta }}}{\text{ (}}\alpha {\text{ = 1, }}0 \leqslant \beta \leqslant 1) 。 (2) 式中:η为分形阻尼器的黏滞系数; \alpha ,β为分形导数阶。
图 2为应力保持100 kPa不变的情况下,不同导数阶β及黏滞系数η的分形阻尼器与Abel阻尼器的流变曲线(分别以 \eta = 40{\text{ kPa}} \cdot {{\text{s}}^\beta } 和 \beta = 0.3 为例)。从图 2中可以看出,分形阻尼器具有同Abel阻尼器相似的描述材料非线性行为的能力,可以很好地表征岩土材料的黏弹性变形和衰减规律。此外,对比二者本构方程还可以发现,Abel阻尼器导数阶 0 \lt \alpha = \beta {\text{ < 1}} 均为分数阶,而分形阻尼器导数阶 \alpha = 1 为整数阶, 0 \leqslant \beta \leqslant {\text{1}} 为分数阶,这一改变使得分形阻尼器能有效规避卷积积分的计算过程,提高计算效率。
1.3 土体损伤演化
Wang等[18]引入损伤变量来描述岩土材料的非线性加速流变特性,取得了理想的效果。考虑损伤效应的分形阻尼器黏度系数η(t)可以表示为
\eta (t) = \eta (1 - {{\text{e}}^{ - \lambda {t^\beta }}})。 (3) 式中,\lambda 为分形阻尼器的损伤因子。
图 3显示了 \sigma = 50{\text{ kPa}} , \eta = 20{\text{ MPa}} \cdot {{\text{s}}^\beta } 时考虑损伤演化的分形阻尼器的流变曲线。从图 3中可以看出,考虑损伤效应的分形阻尼器具有准确表征非线性材料不同程度损伤导致的加速流变的良好能力。
1.4 一维分形导数西元流变模型
传统西元模型流变方程为
\varepsilon (t) = \frac{\sigma }{{{E_0}}} + \frac{\sigma }{{{E_1}}}\left( {1 - {{\text{e}}^{ - \frac{{{E_1}}}{{{\eta _1}}}t}}} \right) + \frac{{\left\langle {\sigma - {\sigma _{\text{s}}}} \right\rangle }}{{{\eta _2}}}t 。 (4) 式中:E0,E1分别为Hooke弹簧及黏弹性元件的弹性模量;η1、η2分别为黏弹性及黏塑性元件的黏滞系数; {\sigma _{\text{s}}} 为屈服应力。
传统西元模型中的Bingham体属于线性元件,无法准确描述岩土体非线性加速流变的特征。Zhou等[7]将西元模型中的Newton阻尼器换成Abel阻尼器,建立了分数阶西元流变模型方程:
\epsilon (t)=\frac{\sigma }{{E}_{0}}+\frac{\sigma }{{\eta }_{1}}{\displaystyle \sum _{n=1}^{\infty }\frac{1}{\Gamma (1+\beta n)}}{\left(-\frac{{E}_{1}}{{\eta }_{1}}{t}^{\beta }\right)}^{n}+\frac{\langle \sigma -{\sigma }_{s}\rangle }{{\eta }_{2}}\frac{{t}^{\beta }}{\Gamma (1+\beta )}。 (5) 式中: \Gamma (*)为伽马函数。
分数阶西元流变模型能够弥补传统模型无法精确描述加速流变的不足,得到了广泛的推广与应用,但它存在卷积积分计算过程复杂的缺陷。
鉴于此,本文采用分形阻尼器代替分数阶西元模型中的Abel阻尼器,且引入岩土材料的损伤效应,构建改进的分形导数西元流变模型,如图 4所示,旨在保证模拟精度的前提下提高计算效率。
根据Merchant体和Bingham体的串联法则,分形导数西元模型的总应力-应变关系为
\left.\begin{array}{l} \varepsilon=\varepsilon^{\mathrm{e}}+\varepsilon^{\mathrm{ve}}+\varepsilon^{\mathrm{vp}}, \\ \sigma=\sigma^{\mathrm{e}}=\sigma^{\mathrm{ve}}=\sigma^{\mathrm{vp}} 。 \end{array}\right\} (6) 式中: \sigma 为总应力;ε为总应变;上标e、ve、vp分别代表Hooke弹簧、黏弹性和黏塑性元件。
分形导数西元模型黏弹性应力-应变关系式为
\left.\begin{array}{l} \sigma^{\mathrm{e}}=E_0 \varepsilon^{\mathrm{e}}, \\ \varepsilon^{\mathrm{ve}}=\varepsilon_{\mathrm{H}}=\varepsilon_{\mathrm{f}}, \\ \sigma^{\mathrm{ve}}=E_1 \varepsilon^{\mathrm{ve}}+\eta_1 \frac{\partial \varepsilon^{\mathrm{ve}}}{\partial t^\beta} 。 \end{array}\right\} (7) 式中:ε为应变;下标H、f分别代表Merchant体中的Hooke弹簧及分形阻尼器;E1和η1分别为二者对应的弹性模量和黏滞系数。
分形导数西元模型黏塑性应力-应变关系式为
\left.\begin{array}{ll} \varepsilon^{\mathrm{vp}}=0 & \left(\sigma<\sigma_{\mathrm{s}}\right), \\ \sigma^{\mathrm{vp}}=\sigma-\sigma_{\mathrm{s}}=\eta_2(t) \frac{\partial \varepsilon^{\mathrm{vp}}}{\partial t^\beta} & \left(\sigma \geqslant \sigma_{\mathrm{s}}\right) \circ \end{array}\right\} (8) 式中: {\eta _2}(t) = {\eta _2}{{\text{e}}^{ - \lambda {t^\beta }{\text{ }}}} 为Bingham体中考虑材料随时间损伤演化的分形阻尼器的黏滞系数。
结合式(6)~(8),进行Laplace变换整合,即可得到分形导数西元流变模型本构方程:
\varepsilon(t)=\left[\frac{1}{E_0}+\frac{1}{E_1}\left(1-\mathrm{e}^{-\frac{E_1 t^\beta}{\eta_1}}\right)\right] \sigma+\frac{\left\langle\sigma-\sigma_{\mathrm{s}}\right\rangle}{\lambda \eta_2}\left(\mathrm{e}^{\lambda t^\beta}-1\right)。 (9) 1.5 三维分形导数西元流变模型
分形导数西元流变模型在三维应力状态下的总应力-应变关系可以表达为
\left.\begin{array}{l} \varepsilon_{i j}=\varepsilon_{i j}^{\mathrm{e}}+\varepsilon_{i j}^{\mathrm{ve}}+\varepsilon_{i j}^{\mathrm{vp}}, \\ \sigma_{i j}=\sigma_{i j}^{\mathrm{e}}=\sigma_{i j}^{\mathrm{ve}}=\sigma_{i j}^{\mathrm{vp}} \quad \circ \end{array}\right\} (10) 式中: {\sigma _{ij}} 为总应力; {\varepsilon _{ij}} 为总应变。
三维应力及应变张量关系可表达为
\sigma_{i j}=s_{i j}+\delta_{i j} \sigma_{\mathrm{m}}, \quad \varepsilon_{i j}=e_{i j}+\delta_{i j} \varepsilon_{\mathrm{m}} \text{,} (11a) \sigma_{\mathrm{m}}=\frac{1}{3}\left(\sigma_1+\sigma_2+\sigma_3\right)=\frac{1}{3} I_1 ; \quad \varepsilon_{\mathrm{m}}=\frac{1}{3}\left(\varepsilon_1+\varepsilon_2+\varepsilon_3\right)=\frac{1}{3} I_1^{\prime}, (11b) s_{i j}=\sigma_{i j}-\delta_{i j} \sigma_{\mathrm{m}}=\sigma_{i j}-\frac{1}{3} I_1 \delta_{i j} ; e_{i j}=\varepsilon_{i j}-\delta_{i j} \varepsilon_{\mathrm{m}}=\varepsilon_{i j}-\frac{1}{3} I_1^{\prime} \delta_{i j}。 (11c) 式中: {\sigma _{\text{m}}} , {\varepsilon _{\text{m}}} 分别为平均应力和应变;I1,Iʹ1分别为第一应力张量和应变张量不变量;sij和eij分别为偏应力和偏应变张量; {\sigma _{ij}} 和 {\varepsilon _{ij}} 分别为应力和应变张量;δij为Kronecker函数。
在三维状态下,由广义Hooke定律可知,弹性体应力与应变之间存在以下转换关系:
\left.\begin{array}{l} e_{i j}=\frac{s_{i j}}{2 G} , \\ \varepsilon_{k k}=\frac{\sigma_{k k}}{3 K} 。 \end{array}\right\} (12) 式中: {\sigma _{kk}} , {\varepsilon _{kk}} 分别代表应力和应变张量第一不变。
故弹性部分应变可表示为
\varepsilon _{ij}^e = \frac{{{s_{ij}}}}{{2G}} + \frac{{{\sigma _m}{\delta _{ij}}}}{{3K}} 。 (13) 弹性模量E、剪切模量G、体积模量K与泊松比 \nu 之间满足
\left.\begin{array}{l} K=\frac{E}{3(1-2 v)}, \\ G=\frac{E}{2(1+v)} 。 \end{array}\right\}。 (14) 分形导数Bingham体黏塑性部分的三维本构关系可以表达为
\frac{\partial \varepsilon^{\mathrm{vp}}}{\partial t^\beta}=\frac{1}{\eta_2(t)}\left\langle\mathit{\Phi}\left(F-F_0\right)\right\rangle 。 (15) 式中:F为屈服函数; {F_0} 为屈服函数初始值;Ф函数取幂指数为1的幂函数形式。
当F-F0 ≥ 0时,由相关联流动法则可以得到
\varepsilon^{\mathrm{vp}}=\frac{\left(\mathrm{e}^{\lambda t^\beta}-1\right)}{\lambda \eta_2}\left\langle\mathit{\Phi}\left(F-F_0\right)\right\rangle 。 (16) 根据Drucker-Prager屈服准则,土体屈服函数可表达为
\left.\begin{array}{l} F=\sqrt{J_2}+a \sigma_{\mathrm{m}}, \\ a=\frac{2 \sin \varphi}{\sqrt{3}(3-\sin \varphi)}, \\ F_0=\frac{6 \cos \varphi}{\sqrt{3}(3-\sin \varphi)} c_{0^{\circ}} \end{array}\right\} (17) 式中:J2为偏应力第二不变量;c0和φ分别为土体的有效黏聚力和内摩擦角。
假设土体弹性应变仅由应力球张量产生,流变变形由应力偏张量产生,则三维应力状态下分形西元模型本构方程可表达为
\varepsilon_{i j}(t)=\frac{s_{i j}}{2 G_0}+\frac{s_{i j}}{2 G_1}\left(1-\mathrm{e}^{-\frac{G_1}{\eta_1^{\prime}} t^\beta}\right)+\frac{\sigma_{\mathrm{m}} \delta_{i j}}{3 K}+\frac{\left(\mathrm{e}^{\lambda t^\beta}-1\right)}{2 \lambda \eta_2^{\prime}}\left\langle\mathit{\Phi}\left(F-F_0\right)\right\rangle 。 (18) 式中,G0,G1分别为Merchant体第一个和第二个弹簧的剪切模量; {\eta '_1} , {\eta '_2} 分别为对应的三维状态黏滞系数。
一维与三维状态黏滞系数转化关系为
\left.\begin{array}{l} \eta_1^{\prime}=\frac{\eta_1}{2(1+v)}, \\ \eta_2^{\prime}=\frac{\eta_2}{2(1+v)} 。 \end{array}\right\} (19) 对于等围压三轴流变试验,围压 {\sigma _3} 已知,且 {\sigma _2} = {\sigma _3} ,轴向应力 {\sigma _1} 为预设常数,可以得到:
{\sigma }_{\text{m}}=\frac{{\sigma }_{1}+2{\sigma }_{3}}{3}\text{;}{s}_{11}={\sigma }_{1}-{\sigma }_{\text{m}}\text{;}\sqrt{{J}_{2}}=\frac{{\sigma }_{1}-{\sigma }_{3}}{\sqrt{3}} 。 (20) \begin{gathered} \varepsilon_{11}(t)=\frac{\sigma_1-\sigma_3}{3 G_0}+\frac{\sigma_1-\sigma_3}{3 G_1}\left(1-\mathrm{e}^{-\frac{G_1}{\eta_1^{\prime}} t^\beta}\right)+\frac{\sigma_1+2 \sigma_3}{9 K}+ \\ \frac{\left\langle\sigma_1-\sigma_3-\sigma_s\right\rangle}{3 \lambda \eta_2^{\prime}}\left(\mathrm{e}^{\lambda t^\beta}-1\right) 。 \end{gathered} (21) 把式(18)~(20)代入式(21)可以得到
\begin{aligned} & \varepsilon_{11}(t)=\frac{2\left(\sigma_1-\sigma_3\right)(1+v)}{3 E_0}+\frac{2\left(\sigma_1-\sigma_3\right)(1+v)}{3 E_1}\left(1-\mathrm{e}^{-\frac{E_1}{\eta_1} t^\beta}\right)+ \\ & \quad \frac{\left(\sigma_1+2 \sigma_3\right)(1-2 v)}{3 E_0}+\frac{\left\langle\sigma_1-\sigma_3-\sigma_{\mathrm{s}}\right\rangle}{3 \lambda \eta_2}\left(\mathrm{e}^{\lambda t^\beta}-1\right) 。 \end{aligned} (22) 2. 分形导数西元流变模型增量形式
分形导数元件的应用使得元件模型只需要较少的参数即可实现岩土材料的复杂非线性特性研究,但在ABAQUS软件模型库中还未得到开发,为实现其在ABAQUS中的应用,利用FORTRAN语言编制了分形导数西元流变模型对应的UMAT子程序,将其嵌入到ABAQUS中,以期实现岩土材料的完整三阶段流变过程分析。
2.1 有限元增量形式推导
时间 {t_n} 对应的第n步分形导数西元流变模型的总应力-应变增量关系为
\Delta \varepsilon_n=\Delta \varepsilon_n^{\mathrm{e}}+\Delta \varepsilon_n^{\mathrm{v}}+\Delta \varepsilon_n^{\mathrm{vp}} ; \quad \Delta \sigma_n=\Delta \sigma_n^{\mathrm{e}}=\Delta \sigma_n^{\mathrm{v}} 。 (23) 矩阵A为
\boldsymbol{A}=\left[\begin{array}{cccccc} 1 & -v & -v & 0 & 0 & 0 \\ -v & 1 & -v & 0 & 0 & 0 \\ -v & -v & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2(1+v) & 0 & 0 \\ 0 & 0 & 0 & 0 & 2(1+v) & 0 \\ 0 & 0 & 0 & 0 & 0 & 2(1+v) \end{array}\right] 。 (24) 对于弹簧元件和分形导数黏弹性体,可以得到
\Delta \varepsilon^{\mathrm{e}}=\frac{\boldsymbol{A}}{E} \Delta \sigma ; \Delta \varepsilon_n^{\mathrm{ve}}=\Delta t_n\left[(1-\theta) \dot{\varepsilon}_n^{\mathrm{ve}}+\theta \dot{\varepsilon}_{n+1}^{\mathrm{ve}}\right] \text{,} (25a) \Delta \varepsilon_n^{\mathrm{ve}}=B_n \Delta t_n \dot{\varepsilon}_n^{\mathrm{ve}}+\boldsymbol{C}_\boldsymbol{n} \Delta \sigma_n (25b) \begin{gathered} \dot{\varepsilon}_n^{\mathrm{ve}}=\frac{\partial \varepsilon_n^{\mathrm{ve}}}{\partial t}=\beta \frac{\sigma}{\eta_1} t^{\beta-1} ; \quad B_n=1-\frac{E_1 \Theta \Delta t_n \beta t_n^{\beta-1}}{\eta_1} ; \\ C_n=\frac{\Theta \boldsymbol{A} \Delta t_n \beta t_n^{\beta-1}}{\eta_1} 。 \end{gathered} (25c) 对于分形导数黏塑性体,可以得到
\Delta \varepsilon_n^{\mathrm{vp}}=\Delta t_n \dot{\varepsilon}_n^{\mathrm{vp}}+\boldsymbol{D}_\boldsymbol{n} \Delta \sigma_n , (26a) \dot{\varepsilon}_n^{\mathrm{vp}}=\frac{\partial \varepsilon_n^{\mathrm{vp}}}{\partial t}=\beta t_n^{\beta-1} \frac{\sigma}{\eta_2} \mathrm{e}^{t^\beta} ; \quad \boldsymbol{D}_\boldsymbol{n}=\frac{\partial \dot{\varepsilon}^{\mathrm{vp}}}{\partial \sigma}=\frac{\Theta A \Delta t_n \beta t_n^{\beta-1}}{\eta_2} \mathrm{e}^{\lambda t^\beta} 。 (26b) 结合式(23)~(26)可以得到分形导数西元流变模型的三维增量形式如下:
\Delta \varepsilon_n=\frac{\boldsymbol{A} \Delta \sigma_n}{E_0}+B_n \Delta t_n \dot{\varepsilon}^{\mathrm{ve}}+\boldsymbol{C}_\boldsymbol{n} \Delta \sigma_n \quad\left(F-F_0<0\right) \text{,} (27a) \begin{gathered} \Delta \varepsilon_n=\frac{\boldsymbol{A} \Delta \sigma_n}{E_0}+B_n \Delta t_n \dot{\varepsilon}^{\mathrm{ve}}+\boldsymbol{C}_\boldsymbol{n} \Delta \sigma_n+\Delta t_n \dot{\varepsilon}^{\mathrm{vp}}+\boldsymbol{D}_\boldsymbol{n} \Delta \sigma_n \\ \left(F-F_0 \geqslant 0\right) 。 \end{gathered} (27b) 2.2 应力更新算法
有限元计算中雅可比矩阵D的定义为 \boldsymbol{D} = \frac{{\partial \Delta \sigma }}{{\partial \Delta \varepsilon }} 。第n步应力与应变增量之间存在以下关系式:
\Delta \sigma_n=E \boldsymbol{A}^{-1} \Delta \varepsilon_n^{\mathrm{e}}=\boldsymbol{D}^{\boldsymbol{\mathrm{e}}} \Delta \varepsilon_n^{\mathrm{e}} 。 (28) 式中: {\boldsymbol{D}^\boldsymbol{\mathrm{e}}} 为弹性刚度矩阵。
根据式(25),(28)可以得到
\Delta \sigma_n\left[\left(\boldsymbol{D}^{\boldsymbol{\mathrm{e}}}\right)^{-1}+C_n\right]=\Delta \varepsilon_n-B_n \Delta t_n \dot{\varepsilon}_n^{\mathrm{ve}} \text{,} (29a) \boldsymbol{D}^{\boldsymbol{\mathrm{ev}}}=\left(\frac{\boldsymbol{A}}{E_0}+\boldsymbol{C}_\boldsymbol{n}\right)^{-1} ; \Delta \sigma_{\mathrm{n}}=\boldsymbol{D}^{\boldsymbol{\mathrm{ev}}}\left(\Delta \varepsilon_{t n}-B_{\mathrm{n}} \Delta t_{\mathrm{n}} \dot{\varepsilon}_{\mathrm{tn}}^{\mathrm{ve}}\right) 。 (29b) 有限元读取上一步应力 {\sigma _n} 与应变量 {\varepsilon _n} , 给出一个应变增量 \Delta {\varepsilon _{n{\text{ + 1}}}} ,先假设应变增量 \Delta {\varepsilon _{n{\text{ + 1}}}} 与黏弹性试应变增量 \Delta {\varepsilon _{tn{\text{ + 1}}}} 相等,按照黏弹性计算公式(30)计算求得黏弹性矩阵 {\boldsymbol{D}D^\boldsymbol{\mathrm{ev}}} 、试应力 {\sigma _{tn{\text{ + 1}}}} 及其增量 \Delta {\sigma _{tn{\text{ + 1}}}} :
\Delta \sigma_{t n+1}=\boldsymbol{D}^{\boldsymbol{\mathrm{ev}}}\left(\Delta \varepsilon_{t n+1}-B_n \Delta t_n \dot{\varepsilon}_{t n+1}^{\mathrm{ve}}\right) ; \quad \sigma_{t n+1}=\sigma_n+\Delta \sigma_{t n+1} 。 (30) 判断材料单元是否达到屈服状态,如果不发生屈服,即F-F0 < 0,此时试应力即为实际应力,本步计算步结束。若F-F0≥0,根据式(26)计算损伤演化产生的塑性应变增量 \Delta \varepsilon _{}^{{\text{vp}}} ,按照下式更新黏弹塑性矩阵 {\boldsymbol{D}^\boldsymbol{\mathrm{evp}}} :
\Delta \sigma_n\left[\left(\boldsymbol{D}^{\boldsymbol{\mathrm{e}}}\right)^{-1}+\boldsymbol{C}_{\boldsymbol{n}}+\boldsymbol{D}_{\boldsymbol{n}}\right]=\Delta \varepsilon_n-B_n \Delta t_n \dot{\varepsilon}_n^{\mathrm{ve}}-\Delta t_n \dot{\varepsilon}_n^{\mathrm{vp}} \text{,} (31a) \boldsymbol{D}^{\boldsymbol{\text {evp }}}=\left(\frac{A}{E_0}+\boldsymbol{C}_\boldsymbol{n}+\boldsymbol{D}_\boldsymbol{n}\right)^{-1} 。 (31b) 然后根据下式计算真实应力值:
\Delta \sigma_{n+1}=\boldsymbol{D}^{\boldsymbol{\mathrm{evp}}}\left(\Delta \varepsilon_{t n+1}-B_n \Delta t_n \dot{\varepsilon}_{t n+1}^{\mathrm{ve}}-\Delta t_n \dot{\varepsilon}_{n+1}^{\mathrm{vp}}\right) \text{,} (32a) \sigma_{n+1}=\Delta \sigma_n+\Delta \sigma_{n+1} 。 (32b) 将上述分行导数西元流变模型的增量形式通过FORTRAN编制成UMAT子程序,计算流程图如图 5所示。
3. 模型参数识别
3.1 三轴流变试验
针对江苏张靖皋长江大桥锚碇基础软土开展三轴流变试验[19],试验用土的基本参数如表 1所示。试样直径为39.1 mm,高度为80 mm。先分别将试样在100,200,400 kPa围压 {\sigma _3} 下进行排水固结;然后,每个围压下分别开展五级轴向偏压q的排水剪切试验,剪切速率为0.01%/min;最后,保持应力不变观察流变特性。试验得到各级围压及应力水平下的流变曲线,如图 6所示。
表 1 江苏软土基本参数Table 1. Properties of Jiangsu soft clay土体 γ /(kN·m-3) e0 w/% wL/% Ιp/% 软土 18.6 0.9 35.4 38 15 从图 6可以看出,不同围压下软土的五级流变曲线呈类似的变化规律:同等围压下轴向应变随偏应力增大而逐渐增加;每级围压低偏应力条件下(前四级)流变曲线呈现衰减型,高偏应力条件下(第五级)轴向应变短时间内迅速增大,试样很快发生破坏,流变曲线呈现非衰减型;同等偏应力下轴向应变随围压增大呈逐渐降低的趋势。
3.2 三轴流变试验模型参数识别
为了验证本文所提分形导数西元流变模型对软土流变特性的表征效果,通过MATLAB软件对三轴流变试验结果进行拟合,得到模型参数,并与传统Burgers、广义Kelvin、西元及分数阶西元模型进行对比。五种模型在200 kPa围压条件下的五级偏应力拟合参数如表 2,3所示,拟合效果如图 7所示。
表 2 Burgers与广义Kelvin模型拟合参数Table 2. Parameters of Burgers and generalized Kelvin model模型 q/
kPaE0/
MPaE1/
MPaη1/
(GPa·sβ)β1 η2/
(GPa·sβ)E2/
MPaβ2 Burgers 117 19.9 34.9 77500 1 780 — 1 220 9.5 76.5 33000 1510 330 6.4 45.3 52800 1040 459 5.5 45.0 38700 2220 550 — — — — 广义Kelvin 117 19.9 45.3 510 1 10400 71.8 1 220 9.5 33.1 18300 850 95.5 330 6.4 60.7 17700 763 54.4 459 5.5 69.5 1300 10300 37.6 550 — — — — — 表 3 西元及分形导数西元模型拟合参数Table 3. Parameters of Nishihara and fractal Nishihara models模
型q/
kPaE0/
MPaE1/
MPaη1/
(GPa·sβ)β1 η2/
(GPa·sβ)β2 λ 西
元117 19.9 30 1000 1 — 1 — 220 9.5 36 5500 330 6.4 33 2100 459 5.5 29 2300 550 4.0 37 6.9 40 分数阶西元 117 19.9 5.60 0.06 0.001 — — — 220 9.5 11.50 0.01 0.002 330 6.4 12.30 0.01 0.002 459 5.5 11.90 0.01 0.003 550 4.0 0.51 0.02 0.17 0.20 0.17 分形西元 117 19.9 4 0.40 0.22 — — — 220 9.5 3.8 0.70 330 6.4 1.0 0.55 459 5.5 0.9 0.47 550 4.0 38 0.15 3 0.55 0.06 从表 2,3可以看出,同等应力状态下5种模型的E0相同,作为流变曲线瞬时弹性变形的影响参数,它既不受流变模型类型的影响又与流变时间无关;同等围压条件下,E0随偏应力增大呈逐渐减小的趋势。从表 3还可以观察到,当分数阶和分形导数西元流变模型的导数阶为1时,即退化为西元模型。
Burgers及广义Kelvin模型没有塑性元件,只能描述衰减型流变。西元模型在应力 \sigma 未达到屈服应力 {\sigma _{\text{s}}} 前,退化为Merchant模型,可用于描述衰减型流变;当应力 \sigma 达到屈服应力 {\sigma _{\text{s}}} 后,塑性元件开始发挥作用,可用于描述包含加速流变在内的非衰减型流变。
从图 7(a)~(d)可以看出,Burgers、广义Kelvin及西元模型能一定程度上描述低应力(前四级)状态下的衰减型流变特性,但对黏弹性变形的初始阶段拟合效果较差,且稳态阶段斜率存在偏差,这种偏差会直接导致预测误差随着时间的增长而增大,不容忽视。而分数阶及分形导数西元流变模型则可以高精度地描述软土前四级状态下的衰减流变特性。从图 7(e)可以看出,西元模型无法准确描述高应力(第五级)状态下的非衰减型流变特性,尤其是加速流变阶段。而分数阶和分形导数西元流变模型因为导数阶的存在能较好地描述非衰减型流变特性。此外,分形西元模型则在分数阶西元模型的基础上进一步改进,更精确地呈现了土体在高应力状态下的非衰减型黏弹塑性流变特性。即,本文提出的分形导数西元流变模型具备准确表征软土完整三阶段流变特性的优良能力且计算更加简洁。
图 8给出了200 kPa围压第五级偏应力条件下分形导数西元流变模型中的两个导数阶对流变曲线的影响。从图 8中可以看出,第一个和第二个导数阶分别影响黏弹性阶段和黏塑性阶段的流变特性。通过对试验数据拟合发现,给定任意一个导数阶初始值时,不同应力状态下的流变曲线拟合得到的导数阶β1、β2分别介于0.20~0.24和0.53~0.56,当指定β1=0.22,β2= 0.55时,通过对不同围压及应力状态流变试验数据拟合发现精度R2 > 0.9,因此,江苏软土的分形导数西元流变模型导数阶可取为β1=0.22,β2=0.55。本文所提分形导数西元流变模型参数相对容易确定,计算简单,应用便捷,具有巨大的推广潜力和应用价值。
4. 模型数值验证
4.1 三轴流变试验验证
采用开发的分形导数西元流变模型UMAT子程序对围压200 kPa条件下的三轴流变试验[19]进行数值模拟,分别模拟围压施加、偏压施加及流变分析过程。UMAT子程序将材料参数E0,E1,η1,η2,β1,β2和λ分别存储于PROPS(1)~PROPS(7)中,按表 3进行赋值,泊松比 \nu 取0.3,并将土体黏弹性应变率、黏塑性应变率、黏弹性应变、黏塑性应变及总应变5个状态变量分别存储于STATEV(1)~STATEV(5)中。有限元模型及流变分析位移云图如图 9所示,分析结果与试验数据对比如图 10所示。
从图 10可以看出,本文编制的分形导数西元流变模型UMAT子程序不仅能很好地计算低应力状态下的衰减型流变变形,还能比较准确地计算高应力状态下的非衰减型流变变形。较传统流变模型具有更优的拟合效果,且能够实现三阶段流变全过程的完整描述,值得推广应用。
4.2 路堤算例数值验证
为了验证分形导数西元流变模型在实际工程中的适用性,对天津某道路工程软土路堤[22]进行了数值模拟。路堤宽30 m,软土厚度为16 m,填土高度5.5 m,软土参数如表 4所示。为了减少计算时间,采用了深16 m、宽45 m的半模型。将路堤荷载简化为从0到110 kPa线性增加的均布荷载,施加宽度为15 m。模型采用平面应力单元,右侧边界限制水平位移,底部边界限制竖向位移。土层简化为粉质黏土层,软土流变参数由文献[22]中所给出的围压100 kPa偏应力60 kPa条件下的三轴流变试验曲线拟合得到,见表 5。数值模拟过程分三步:先对地基土建立初始应力条件;然后将路堤载荷施加在软土地基上,计算其随荷载施加产生的沉降;最终维持荷载不变,计算其随时间变化产生的流变变形。图 11为地基土体数值模拟结束时(400 d)的沉降云图,图 12为地基土体顶部路堤中心线以下沉降的数值模拟和监测结果的对比效果。
表 4 土体参数Table 4. Properties of soils土体 深度/m c/kPa γ /(kN·m-3) ΙL Ιp 淤泥 1 12.5 15.9 1.09 32.7 粉质黏土 15 15.1 17.9 1.13 22.3 表 5 土体流变参数Table 5. Rheological parameters of soilsσ3/kPa q/kPa E0/MPa E1/MPa η1/ (GPa·sβ) β R2 100 60 2.48 0.08 12.2 0.56 0.967 从图 12可以看出,本文所提的分形导数西元流变模型数值模拟分析结果与实测值较为吻合,说明编制的分形导数西元流变模型UMAT子程序合理可靠,适用于分析软土的长期流变特性。研究结果可为ABAQUS中描述完整流变过程的流变模型的二次开发及应用提供参考。
5. 结论
本文基于分形导数理论建立了分形导数西元黏弹塑性流变模型,给出了模型的解析解,并将其推广至三维。以ABAQUS软件为平台,编制了相应的UMAT子程序,通过软土的三轴流变试验及路堤工程实例对模型及其开发程序进行了验证,得到以下4点结论。
(1)采用分形阻尼器代替传统模型中的Newton阻尼器和Abel阻尼器,引入损伤因子,提出了精度高且计算简便的改进型分形导数西元黏弹塑性流变模型,并给出了一维及三维解析解。
(2)通过三轴流变试验分别对传统元件模型和分形导数西元黏弹塑性模型参数进行识别,得出了分形导数西元黏弹塑性流变模型的导数阶与应力状态无关的结论。江苏软土的分形导数西元黏弹塑性流变模型导数阶可分别取为0.22和0.55。
(3)通过三轴试验与传统模型及所提模型拟合结果对比发现,与传统模型相比,分形导数西元黏弹塑性模型不仅能准确地描述低应力条件下的黏弹性衰减型流变特性,还能很好地呈现高应力条件下完整三阶段黏弹塑性非衰减型流变全过程。
(4)推导出了分形导数西元黏弹塑性流变模型解析解的增量型式,编制了相应的UMAT子程序,开展了三轴流变试验及路堤工程数值模拟。数值模拟结果与测试结果吻合度较好,进一步验证了模型及开发程序的合理性,可为ABAQUS中流变模型的二次开发及应用提供参考。
利益冲突声明/Conflict of Interests:所有作者声明不存在利益冲突。All authors disclose no relevant conflict of interest.作者贡献/Authors' Contributions:尹倩、竺明星参与试验设计、理论推导及论文撰写;龚维明、戴国亮参与论文修改及提供资金资助;尹倩、王博臣参与数值模拟编程及分析;尹倩、阮静、刘强、张凡参与工程案例分析及论文修改。所有作者均阅读并同意最终稿件的提交。YIN Qian and ZHU Mingxing were involved in the experimental design, theoretical derivation and writing; GONG Weiming and DAI Guoliang were involved in the paper revision and financial support; YIN Qian and WANG Bochen were involved in the numerical simulation programming and analysis; YIN Qian, RUAN Jing, LIU Qiang and ZHANG Fan were involved in the engineering case study and manuscript revision. All authors have read and agreed to the final manuscript submission. -
表 1 江苏软土基本参数
Table 1 Properties of Jiangsu soft clay
土体 γ /(kN·m-3) e0 w/% wL/% Ιp/% 软土 18.6 0.9 35.4 38 15 表 2 Burgers与广义Kelvin模型拟合参数
Table 2 Parameters of Burgers and generalized Kelvin model
模型 q/
kPaE0/
MPaE1/
MPaη1/
(GPa·sβ)β1 η2/
(GPa·sβ)E2/
MPaβ2 Burgers 117 19.9 34.9 77500 1 780 — 1 220 9.5 76.5 33000 1510 330 6.4 45.3 52800 1040 459 5.5 45.0 38700 2220 550 — — — — 广义Kelvin 117 19.9 45.3 510 1 10400 71.8 1 220 9.5 33.1 18300 850 95.5 330 6.4 60.7 17700 763 54.4 459 5.5 69.5 1300 10300 37.6 550 — — — — — 表 3 西元及分形导数西元模型拟合参数
Table 3 Parameters of Nishihara and fractal Nishihara models
模
型q/
kPaE0/
MPaE1/
MPaη1/
(GPa·sβ)β1 η2/
(GPa·sβ)β2 λ 西
元117 19.9 30 1000 1 — 1 — 220 9.5 36 5500 330 6.4 33 2100 459 5.5 29 2300 550 4.0 37 6.9 40 分数阶西元 117 19.9 5.60 0.06 0.001 — — — 220 9.5 11.50 0.01 0.002 330 6.4 12.30 0.01 0.002 459 5.5 11.90 0.01 0.003 550 4.0 0.51 0.02 0.17 0.20 0.17 分形西元 117 19.9 4 0.40 0.22 — — — 220 9.5 3.8 0.70 330 6.4 1.0 0.55 459 5.5 0.9 0.47 550 4.0 38 0.15 3 0.55 0.06 表 4 土体参数
Table 4 Properties of soils
土体 深度/m c/kPa γ /(kN·m-3) ΙL Ιp 淤泥 1 12.5 15.9 1.09 32.7 粉质黏土 15 15.1 17.9 1.13 22.3 表 5 土体流变参数
Table 5 Rheological parameters of soils
σ3/kPa q/kPa E0/MPa E1/MPa η1/ (GPa·sβ) β R2 100 60 2.48 0.08 12.2 0.56 0.967 -
[1] 孙钧. 岩土材料流变及其工程应用[M]. 北京: 中国建筑工业出版社, 1999. SUN Jun. Rheology of Geotechnical Materials and its Engineering Application[M]. Beijing: China Architecture & Building Press, 1999. (in Chinese)
[2] HOU R B, ZHANG K, TAO J, et al. A nonlinear creep damage coupled model for rock considering the effect of initial damage[J]. Rock Mechanics and Rock Engineering, 2019, 52(5): 1275-1285. doi: 10.1007/s00603-018-1626-7
[3] MESRI G, FEBRES-CORDERO E, SHIELDS D R, et al. Shear stress-strain-time behaviour of clays[J]. Géotechnique, 1981, 31(4): 537-552. doi: 10.1680/geot.1981.31.4.537
[4] LEROUEIL S, KABBAJ M, TAVENAS F, et al. Stress–strain–strain rate relation for the compressibility of sensitive natural clays[J]. Géotechnique, 1985, 35(2): 159-180. doi: 10.1680/geot.1985.35.2.159
[5] NISHIHARA M. Rheological properties of rocks I and II[J]. Doshisha Engineering Review, 1958, 8, 85-115.
[6] 谷任国, 邹育, 房营光, 等. 引入非线性瞬时弹性模量的软土流变模型[J]. 岩土力学, 2018, 39(1): 237-241. GU Renguo, ZOU Yu, FANG Yingguang, et al. Rheological model of soft soils using nonlinear instantaneous elastic modulus[J]. Rock and Soil Mechanics, 2018, 39(1): 237-241. (in Chinese)
[7] ZHOU H W, WANG C P, HAN B B, et al. A creep constitutive model for salt rock based on fractional derivatives[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(1): 116-121. doi: 10.1016/j.ijrmms.2010.11.004
[8] WANG L Y, ZHOU F X. Analysis of elastic-viscoplastic creep model based on variable-order differential operator[J]. Applied Mathematical Modelling, 2020, 81: 37-49. doi: 10.1016/j.apm.2019.12.007
[9] ZHOU F X, WANG L Y, LIU H B. A fractional elasto-viscoplastic model for describing creep behavior of soft soil[J]. Acta Geotechnica, 2021, 16(1): 67-76. doi: 10.1007/s11440-020-01008-5
[10] LU D C, LIANG J Y, DU X L, et al. Fractional elastoplastic constitutive model for soils based on a novel 3D fractional plastic flow rule[J]. Computers and Geotechnics, 2019, 105: 277-290. doi: 10.1016/j.compgeo.2018.10.004
[11] 薛东杰, 路乐乐, 易海洋, 等. 考虑温度和体积应力的分数阶蠕变损伤Burgers模型[J]. 岩石力学与工程学报, 2021, 40(2): 315-329. XUE Dongjie, LU Lele, YI Haiyang, et al. A fractional Burgers model for uniaxial and triaxial creep of damaged salt-rock considering temperature and volume-stress[J]. Chinese Journal of Rock Mechanics and Engineering, 2021, 40(2): 315-329. (in Chinese)
[12] WU F, LIU J, ZOU Q L, et al. A triaxial creep model for salt rocks based on variable-order fractional derivative[J]. Mechanics of Time-Dependent Materials, 2021, 25(1): 101-118. doi: 10.1007/s11043-020-09470-0
[13] 张胜利, 梁卫国, 肖宁, 等. 考虑温度的盐岩分数阶黏弹塑性蠕变损伤模型[J]. 岩石力学与工程学报, 2022, 41(增刊2): 3198-3209. ZHANG Shengli, LIANG Weiguo, XIAO Ning, et al. Fractional viscoelastic-plastic creep damage model of salt rock considering temperature[J]. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(S2): 3198-3209. (in Chinese)
[14] 孙增春, 刘汉龙, 肖杨. 砂-粉混合料的分数阶塑性本构模型[J]. 岩土工程学报, 2024, 46(8): 1596-1604. doi: 10.11779/CJGE20230666 SUN Zengchun, LIU Hanlong, XIAO Yang. Fractional-order plasticity model for sand-silt mixtures[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(8): 1596-1604. (in Chinese) doi: 10.11779/CJGE20230666
[15] XU X B, CUI Z D. Investigation of a fractional derivative creep model of clay and its numerical implementation[J]. Computers and Geotechnics, 2020, 119: 103387.
[16] CHEN W. Time–space fabric underlying anomalous diffusion[J]. Chaos, Solitons & Fractals, 2006, 28(4): 923-929.
[17] CAI W, CHEN W, XU W X. Characterizing the creep of viscoelastic materials by fractal derivative models[J]. International Journal of Non-Linear Mechanics, 2016, 87: 58-63.
[18] WANG R, ZHUO Z, ZHOU H W, et al. A fractal derivative constitutive model for three stages in granite creep[J]. Results in Physics, 2017, 7: 2632-2638.
[19] YIN Q, ZHAO Y, GONG W M, et al. A fractal order creep-damage constitutive model of silty clay[J]. Acta Geotechnica, 2023, 18(8): 3997-4016.
[20] 张亮亮, 王晓健. 岩石黏弹塑性损伤蠕变模型研究[J]. 岩土工程学报, 2020, 42(6): 1085-1092. doi: 10.11779/CJGE202006012 ZHANG Liangliang, WANG Xiaojian. Viscoelastic-plastic damage creep model for rock[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(6): 1085-1092. (in Chinese) doi: 10.11779/CJGE202006012
[21] 骆亚生, 赵程斌, 孙哲, 等. 基于塑性元件微元化的重塑黄土黏弹塑性本构模型[J]. 岩土工程学报, 2024, 46(3): 624-631. doi: 10.11779/CJGE20221452 LUO Yasheng, ZHAO Chengbin, SUN Zhe, et al. Visco-elastoplastic constitutive model for remolded loess based on infinitesimalization of plastic elements[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(3): 624-631. (in Chinese) doi: 10.11779/CJGE20221452
[22] 袁宇, 刘润, 邱长林, 等. 与应力相关的软土蠕变本构模型的建立和应用[J]. 天津大学学报(自然科学与工程技术版), 2018, 51(7): 711-719. YUAN Yu, LIU Run, QIU Changlin, et al. Establishment and application of creep constitutive model related to stress level of soft soil[J]. Journal of Tianjin University (Science and Technology), 2018, 51(7): 711-719. (in Chinese)
[23] AI Z Y, ZHAO Y Z, DAI Y C, et al. Three-dimensional elastic–viscoplastic consolidation behavior of transversely isotropic saturated soils[J]. Acta Geotechnica, 2022, 17(9): 3959-3976.
[24] 朱文波, 戴国亮, 王博臣, 等. 吸力式沉箱基础底部土体卸荷蠕变及其长期抗拔承载特性研究[J]. 岩土力学, 2022, 43(3): 669-678. ZHU Wenbo, DAI Guoliang, WANG Bochen, et al. Unloading creep of soft clay and long-term uplift bearing characteristics of suction caisson foundation[J]. Rock and Soil Mechanics, 2022, 43(3): 669-678. (in Chinese)
[25] 元捷衡, 王元战, 王轩, 等. 碱渣土蠕变特性试验及数值模拟研究[J]. 土木工程学报, 2024, 57(5): 107-116. YUAN Jieheng, WANG Yuanzhan, WANG Xuan, et al. Experimental and numerical simulation studies on creep behavior of soda residue soil[J]. China Civil Engineering Journal, 2024, 57(5): 107-116. (in Chinese)
-
期刊类型引用(2)
1. 刘超凡,杜强,曲立强,刘旭梅. 基于显微镜观测的不同钙源下微生物固化砂强度提升机理的宏细观研究. 内蒙古农业大学学报(自然科学版). 2025(02): 63-71 . 百度学术
2. 刘小军,朱芮,王东亚,孙跃. 冻融循环下MICP固化风积砂的性能劣化试验研究. 中国安全生产科学技术. 2025(04): 71-78 . 百度学术
其他类型引用(4)
-
其他相关附件