基于ABAQUS的土体地基液化研究方法

技术领域

本发明涉及一种基于ABAQUS的土体地基液化研究方法，属于水利工程技术 领域。

背景技术

土体液化引起的建筑物破坏是地震灾害的主要形式之一。在以往的地震灾 难中，如1964年的美国阿拉斯加地震、1976年的唐山地震、2008年的汶川地 震等，都出现了大面积的地层液化现象，并引起了相应建筑物的极大破坏。因 此，地震作用下土体地基的液化判别方法一直是抗震减灾设计中的重要问题之 一。

目前，国内相关设计规范对土体液化的分析方法主要是以经验分析法为主， 基本框架及数据来源都是上世纪发生的几次地震的经验总结，因而在理论和应 用上都存着一些不足之处，如液化影响折减系数的“分档”取值跳跃性较大、 倾斜场地的液化判别未见到相应标准等，对工程结构的安全性以及工程造价等 方面产生了一定的影响。随着液化理论和计算方法的不断发展，结合数值方法 对土体液化问题进行深入研究，已成为提高土体地基液化分析的合理性和可靠 性的一个重要课题。

ABAQUS是目前国际上功能强大的有限元软件之一，其在前处理、计算精度 和复杂工程问题模拟等方面的优势已得到广泛认可，但其在土体材料的液化模 拟方面还存在一些不足。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种基于ABAQUS的土体地 基液化研究方法，解决现有技术中液化场地设计中土层液化影响折减系数的分 档取值范围跳跃性较大的技术问题，并给出了倾斜场地的液化判别标准。

为达到上述目的，本发明是所采用的技术方案是：1.基于ABAQUS的土体地 基液化研究方法，包括如下步骤：

步骤一：在ABAQUS软件中二次开发Byrne液化模型：根据土体在地震过程 中的液化变形性质，结合ABAQUS软件中二次开发平台，将Byrne液化模型嵌入 ABAQUS的计算分析中，实现ABAQUS对土体地基液化过程的数值模拟；

步骤二：水平场地液化影响折减系数取值精细化处理：分析液化判别标准 中规范法与数值模拟法间各参数的相互关系，通过逐渐改变抗液化指数I_N值， 对水平场地地基模型液化过程的数值模拟，根据液化土层各参数变化规律，给 出水平场地液化影响折减系数的精细化处理计算公式；

步骤三：倾斜场地的液化判别研究：在步骤二的基础上，对倾斜场地地基 模型进行模拟，逐渐改变地基倾斜角度，根据液化土层的液化程度随倾斜角度 的变化情况，给出不同倾角倾斜场地的抗液化指数与超孔隙水压比之间的对应 关系，并根据此对应关系对倾斜场地液化指数进行修正，达到对倾斜场地进行 液化判别的目的。

步骤一的具体步骤如下：

步骤101)：假定单元的超孔隙水压力增量Δu与体积应变增量Δε_vd相关，其 公式的简要表达式如下：

Δu＝E_rΔε_vd(1)

E_r为饱和土体的回弹模量，其表达式如下所示：

$E_r=\frac{{({{\sigma }_{v0}}^{\prime }-u)}^{1-m}}{{\mathrm{mk}}_2{\left({{\sigma }_{v0}}^{\prime }\right)}^{n-m}}---\left(2\right)$

式中，σ_v0′为初始铅直向有效应力；u为超孔隙水压力,i为计算迭代次 数；m、n、k₂为土体相关的经验常数；

Δε_vd为体积应变增量，其表达式如下所示：

${\mathrm{\Delta ϵ}}_{vd}={\mathrm{\gamma C}}_1\exp (-C_2\frac{{ϵ}_{vd}}{\gamma })---\left(3\right)$

式中，ε_vd为土单元累计体积应变，i为计算迭代次数；γ为第N次 应力循环的剪应变；C₁、C₂为材料参数，可由标准贯入试验锤击数(N₁)₆₀确定， 定义为C₂＝0.4/C₁，其中，(N₁)₆₀＝N/σ_v′，N为液化土层实测标准 贯入锤击数，σ_v′为上覆土层有效压力；

步骤102)：在每一次应力循环后对土体材料的弹性模量进行折减，达到对 其体剪切模量修正的目的，计算公式如下所示：

$E_n=E_0{\left(\frac{{{\sigma }_{v0}}^{\prime }-u}{{{\sigma }_{v0}}^{\prime }}\right)}^x---\left(4\right)$

式中，E_n为第n周应力循环弹性模量；E₀为地震初始弹性模量；x为与土性有关 参数；当土体完全液化时，取此时的应力循环弹性模量为地震初始弹性模量的 1‰；

步骤103)：根据土体地震液化情况对模型的静、动力分析采用分开计算：

首先，在动力分析之前，对模型进行静力的流固耦合计算，并进行地应力 平衡，此时得到的有效应力场、渗流场为震前的初始状态，将其结果导出作为 动力计算的初始条件；

然后，对模型进行动力分析，基于ABAQUS二次开发平台，编写VUSDFLD子 程序，在计算过程中通过场变量传递，读入静力部分的结果并赋予模型，并根 据液化模型对土体地基的地震液化情况进行动力分析，给出地震过程中振动孔 压上升模式和有效应力的变化过程，进而对其液化情况进行分析。

步骤二的具体步骤如下：

步骤201)：以某一水平场地地基模型为研究对象，根据模型网格图确定水 平主向、水平次向以及竖向的深度；

步骤202)：将土层分为三层，上层和底层分别为非液化土层，中层为可液 化土层；

步骤203)：确定水平场地地基模型参数，包括液化判别的标准贯入锤击数 基准值N₀、地下水位埋深d_w、黏粒含量百分率ρ_c；ABAQUS时程分析中阻尼按瑞 利阻尼施加，模型各土层材料的物理参数如表1所示：

表1地基模型各土层材料物理参数

步骤204)：设计水平场地地基模型的地震分组为第一组，设计地震基本加 速度、地震反应谱特征周期、持续时间；

步骤205)：考虑三向地震波同时作用，根据地震影响系数曲线，人工合成 三条地震波加速度时程曲线，根据设计规范，确定加速度峰值比例；

步骤206)：根据液化土层的抗液化指数，采用分档的方法对液化影响减系 数进行取值，土层的抗液化指数计算公式如下所示：

I_N＝N/N_cr(5)

式中，N为液化土层实测标准贯入锤击数，N_cr为液化土层液化判别标准贯入锤 击数临界值，其公式如下所示：

$N_{cr}=N_0\beta [ln(0.6d_s+1.5)-0.1d_w]\sqrt{3/{\rho }_c}---\left(6\right)$

式中，N₀为液化判别的标准贯入锤击数基准值；d_s为饱和土标准贯入点深 度；d_w为地下水位；ρ_c为黏粒含量百分率；β为调整系数；

根据式(5)计算液化土层实测标准贯入锤击数N，将N代入式(3)计算 超孔隙水压力增量Δu，从而得到为超孔隙水压力u，根据如下公式计算液化土 层的超孔隙水压比r_u：

$r_u=\frac{u}{{{\sigma }_{v0}}^{\prime }}\times 100\%---\left(7\right)$

式中，σ_v0′为初始铅直向有效应力；

将式(7)代入式(4)，可得液化影响折减系数α与超孔隙水压比r_u之间的 关系为

α＝(1-r_u)^x(8)；

步骤207)：取地基的抗液化指数I_N的变化范围为0.1～1之间，根据建立的 计算流程，模拟场地模型的地震液化情况，记录可液化土层中三个监测点在液 化过程中超孔隙水压比r_u的变化情况，定义三点的加权平均值为液化土层的超孔 隙水压比值，则可得到液化土层的超孔隙水压比r_u与抗液化指数I_N、液化影响折 减系数α与抗液化指数I_N的对应关系曲线，由关系曲线得出：液化影响折减系 数α是随着抗液化指数I_N连续变化，表明采用数值方法对规范法液化影响折减 系数的取值进行精细化处理是可行的；

步骤208)：应用Matlab工具箱对曲线数据进行拟合，去除α与I_N关系曲线 中超孔隙水压比r_u＝1的I_N所对应的点，得到曲线拟合关系式，如下所示：

α＝-0.61I_N²+1.67I_N-0.4(9)

根据α与I_N关系曲线和式(9)可对地基模型的液化判别标准进行精细化处 理。

步骤三的具体步骤如下：

步骤301)：数值模拟条件与步骤二相同，将地基左侧抬高，以达到倾斜目 的；

步骤302)：对不同倾角的液化土层进行计算，记录不同倾角下液化土层的 超孔隙水压比r_u与抗液化指数I_N对应关系，并绘出结果关系曲线；

步骤303)：取水平场地模型计算结果为基数，对比倾斜场地中液化土层超 孔隙水压比r_u(β)随倾角β的变化情况；

步骤304)：应用Matlab工具箱对r_u(β)与β进行关系曲线拟合，得到倾斜场 地中液化土层超孔隙水压比r_u(β)关于水平场地液化土层超孔隙水压比r_u(0)和倾 角β的函数关系式为

$r_u\left(\beta \right)=\{\begin{array}{cc}{r}_u\left(0\right)& 0\le I_N\le 0.2\\ -0.027-0.0016\beta +1.019r_u\left(0\right)& 0\le I_N\le 1.0\end{array}---\left(11\right)$

步骤305)：液化土层的超孔隙水压比r_u(β)和场地的液化指数I_lE都是对场地 模型液化程度的表征，认为两者在数值上成正比的关系，则可得倾斜场地液化 指数为：

$I_{lE}\left(\beta \right)=\{\begin{array}{cc}{I}_{lE}\left(0\right)& 0\le I_N\le 0.2\\ -0.027-0.0016\beta +1.019r_u\left(0\right)& 0\le I_N\le 0.2\end{array}---\left(12\right)$

式中，β为倾斜场地的倾斜角度，I_lE(β)为倾角为β的倾斜场地液化指数，I_lE(0)为 水平场地液化指数。

与现有技术相比，本发明所达到的有益效果是：

(1)本发明将土体液化分析和ABAQUS分析平台进行了有机结合，实现了地 基场地的地震液化过程在ABAQUS中的数值模拟，具有一定的科研和实用价值；

(2)本发明通过对水平场地的液化过程分析，给出了规范法中液化影响系 数分档取值问题的精细化处理方案，使设计方案更趋合理化。

(3)本发明通过对不同倾斜角度液化场地的液化过程分析，给出了不同倾 角倾斜场地的抗液化指数与超孔隙水压比之间的对应关系；进而指出现行规范 中的液化判别方法应用于倾斜场地时需要进行修正，并给出了一种修正方法， 以达到对倾斜场地进行液化判别的目的。这可以对“规范法”不适用于倾斜场 地的液化判别这一不足之处做出补充，提高了工程的安全性。

(4)地基液化判别和抗液化设计一直是土木工程中的重要问题。本发明的 研究思路及研究方法可以为实际工程的抗液化设计提供较好的参考和借鉴作 用；对涉及到岸坡液化问题的相关工程，应用本发明的研究内容，可以使得抗 液化设计方案更趋合理。

附图说明

图1是ABAQUS的土体液化计算流程图。

图2是液化场地模型网格图。

图3是人工生成地震波X方向加速度时程曲线。

图4是人工生成地震波Y方向加速度时程曲线。

图5是人工生成地震波Z方向加速度时程曲线。

图6是水平场地超孔隙水压比与抗液化指数关系曲线。

图7是液化影响折减系数与抗液化指数关系曲线。

图8是倾斜场地超孔隙水压比与抗液化指数关系曲线。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明 本发明的技术方案，而不能以此来限制本发明的保护范围。

基于ABAQUS的土体地基液化研究方法，包括如下步骤：

步骤一：在ABAQUS软件中二次开发Byrne液化模型：根据土体在地震过程 中的液化变形性质，结合ABAQUS软件中二次开发平台，将Byrne液化模型嵌入 ABAQUS的计算分析中，实现ABAQUS对土体地基液化过程的数值模拟。Byrne液 化模型是在Finn液化模型的基础上进行简化得到的，其都是在Mohr-Coulomb 模型的基础上增加了振动孔压的上升模式。具体步骤如下：

步骤101)：假定单元的超孔隙水压力增量Δu与体积应变增量Δε_vd相关，其 公式的简要表达式如下：

Δu＝E_rΔε_vd(1)

E_r为饱和土体的回弹模量，其表达式如下所示：

$E_r=\frac{{({{\sigma }_{v0}}^{\prime }-u)}^{1-m}}{{\mathrm{mk}}_2{\left({{\sigma }_{v0}}^{\prime }\right)}^{n-m}}---\left(2\right)$

式中，σ_v0′为初始铅直向有效应力；u为超孔隙水压力，i为计算 迭代次数；m、n、k₂为土体相关的经验常数；

Δε_vd为体积应变增量，其表达式如下所示：

${\mathrm{\Delta ϵ}}_{vd}={\mathrm{\gamma C}}_1\exp (-C_2\frac{{ϵ}_{vd}}{\gamma })---\left(3\right)$

式中，ε_vd为土单元累计体积应变，i为计算迭代次数；γ为第N次 应力循环的剪应变；C₁、C₂为材料参数，可由标准贯入试验锤击数(N₁)₆₀确定， 定义为C₂＝0.4/C₁，其中，(N₁)₆₀＝N/σ_v′，N为实测标准贯入锤击 数，σ_v′为上覆土层有效压力。

步骤102)：在计算过程中，土体的应力-应变关系遵循Mohr-Coulomb本构 模型，为考虑在应力循环过程中土体有效应力变化对其力学性能的影响，在每 一次应力循环后对土体材料的弹性模量进行折减，达到对其体剪切模量修正的 目的，计算公式如下所示：

$E_n=E_0{\left(\frac{{{\sigma }_{v0}}^{\prime }-u}{{{\sigma }_{v0}}^{\prime }}\right)}^x---\left(4\right)$

式中，E_n为第n周应力循环弹性模量；E₀为地震初始弹性模量；x为与土性有关 参数，对于砂土材料可取x＝0.5。当土体完全液化时，取此时的应力循环弹性模 量为地震初始弹性模量的1‰。

步骤103)：根据土体地震液化情况对模型的静、动力分析采用分开计算：

首先，在动力分析之前，对模型进行静力的流固耦合计算，并进行地应力 平衡，此时得到的有效应力场、渗流场为震前的初始状态，将其结果导出作为 动力计算的初始条件；

然后，对模型进行动力分析，基于ABAQUS二次开发平台，编写VUSDFLD子 程序，在计算过程中通过场变量传递，读入静力部分的结果并赋予模型，并根 据液化模型对土体地基的地震液化情况进行动力分析，给出地震过程中振动孔 压上升模式和有效应力的变化过程，进而对其液化情况进行分析。

通过上述方法即可实现土体地基液化过程在ABAQUS中的数值模拟，计算流 程如图1所示。

步骤二：水平场地液化影响折减系数取值精细化处理：分析液化判别标准 中规范法与数值模拟法间各参数的相互关系，通过逐渐改变抗液化指数I_N值， 对水平场地地基模型液化过程的数值模拟，根据液化土层各参数变化规律，给 出水平场地液化影响折减系数的精细化处理计算公式。

下面分别对规范法液化影响折减系数取值和数值模拟法液化影响折减系数 取值进行分析。

(1)“规范法”液化影响折减系数取值

目前，液化场地工程设计时，一般采用“规范法”，根据液化土层的抗液化 指数，采用分档的方法对液化影响减系数进行取值，土层的抗液化指数计算公 式如下所示：

I_N＝N/N_cr(5)

式中，N为液化土层实测标准贯入锤击数，N_cr为液化土层液化判别标准贯入锤 击数临界值，其公式如下所示：

$N_{cr}=N_0\beta [ln(0.6d_s+1.5)-0.1d_w]\sqrt{3/{\rho }_c}---\left(6\right)$

式中，N₀为液化判别的标准贯入锤击数基准值；d_s为饱和土标准贯入点深度(m)； d_w为地下水位(m)；ρ_c为黏粒含量百分率；β为调整系数，与设计地震相关。

分档的取值方法对工程的设计、施工存在着较大的影响，因为在实际工程 中某些参数的微小变化，就可能会导致取值结果出现在不同的档位，从而造成 设计结果出现质的改变，对工程安全和工程成本产生重大影响。因此，对“规 范法”中液化影响折减系数的取值方法进行精细化处理，是液化土体抗液化设 计的一个重要方向。

(2)数值模拟法液化影响折减系数取值

根据有效应力原理，在数值模拟方法中，为了判别所研究土体是否液化， 一般直接根据计算得到的土体有效应力是否为“0”进行判别；也有研究者根据 超孔隙水压力与初始有效应力之比r_u来判断土体是否液化，称为超孔隙水压比 法，r_u的表达式如下所示：

$r_u=\frac{u}{{{\sigma }_{\nu 0}}^{\prime }}\times 100\%---\left(7\right)$

式中，σ_v0′为初始铅直向有效应力；

将式(7)代入式(4)，可得液化影响折减系数α与超孔隙水压比r_u之间的 关系为

α＝(1-r_u)^x(8)

对比数值模拟方法和“规范法”，从式(3)、(5)、(8)中不难发现，在材 料各项力学参数确定的前提下，场地的超孔隙水压比r_u、抗液化指数I_N、液化影 响折减系数α是一一对应的。

因此，可结合“规范法”和数值模拟方法，通过研究某一液化层的抗液化 指数I_N随超孔隙水压比r_u的变化规律，探讨对常规设计中液化影响折减系数α的 取值进行精细化处理的可能性。下面结合具体实施例对本步骤作具体说明，步 骤如下：

步骤201)：以某一水平场地地基模型为研究对象，模型网格图如图2所示， 根据模型网格图确定水平主向、水平次向以及竖向的深度：假设竖向(Z方向) 深度为20m，水平主向(X方向)、水平次向(Y方向)均为40m。

步骤202)：将土层分为三层，上层和底层分别为非液化土层，中层为可液 化土层：假设可液化土层厚度为4m，其中心位置距顶面距离为4m。

步骤203)：确定水平场地地基模型参数，包括液化判别的标准贯入锤击数 基准值N₀、地下水位埋深d_w、黏粒含量百分率ρ_c：水平场地地基模型参数根据 “大唐南京电厂工程场地地震安全性评价工作报告”取值，液化判别标准贯入 锤击数基准值为8击，地下水位埋深统一按0.5m考虑，粘粒含量取3％。

ABAQUS时程分析中阻尼按瑞利阻尼施加，模型各土层材料的物理参数如表 1所示：

表1地基模型各土层材料物理参数

步骤204)：设计水平场地地基模型的地震分组为第一组，设计地震基本加 速度、地震反应谱特征周期、持续时间：地震基本加速度取0.15g，地震反应谱 特征周期取0.5s，持续时间为20s。

步骤205)：考虑三向地震波同时作用，根据地震影响系数曲线，人工合成 三条地震波加速度时程曲线，根据设计规范，确定加速度峰值比例：如图3至 图5所示，加速度峰值比例为：X向：Y向：Z向＝1.00：0.85：0.65。

步骤206)：根据液化土层的抗液化指数，采用分档的方法对液化影响减系 数进行取值，具体为：根据式(5)计算液化土层实测标准贯入锤击数N，将N 代入式(3)计算超孔隙水压力增量Δu，从而得到为超孔隙水压力u的值，代入 式(7)计算液化土层的超孔隙水压比r_u的值，将液化土层的超孔隙水压比r_u的 值代入式(4)，可得液化影响折减系数α与超孔隙水压比r_u之间的关系，即 α＝(1-r_u)^x；

步骤207)：取地基的抗液化指数I_N的变化范围为0.1～1之间，根据建立的 计算流程，模拟场地模型的地震液化情况，记录可液化土层中三个监测点(图2 中a、b、c三点)在液化过程中超孔隙水压比r_u的变化情况，定义三点的加权平 均值为液化土层的超孔隙水压比值，则可得到：液化土层的超孔隙水压比r_u与抗 液化指数I_N的对应关系曲线，如图6所示；以及液化影响折减系数α与抗液化 指数I_N的对应关系曲线，如图7所示。由图7给出的关系曲线得出：液化影响 折减系数α是随着抗液化指数I_N连续变化，表明采用数值方法对规范法液化影 响折减系数的取值进行精细化处理是可行的。

步骤208)：应用Matlab工具箱对曲线数据进行拟合，如图7中虚线点A所 示，在I_N处于0.2～0.3之间时地基场地的r_u值已达到1。因此，在拟合过程中 去除α与I_N关系曲线中超孔隙水压比r_u＝1的I_N所对应的点，得到曲线拟合关系 式，如下所示：

α＝-0.61I_N²+1.67I_N-0.4(9)

根据图7和式(9)可对地基模型的液化判别标准进行精细化处理，得到实 例中标准贯入点深度为4m的水平地基，其土层液化影响系数取值标准，如表2 所示：

表2土层液化影响折减系数

目前，“规范法”中对液化场地的液化等级判别主要通过液化指数的方法，

对含有n层液化土层的地基，液化指数计算公式如下所示：

$I_{lE}=\underset{i}{\overset{n}{\Sigma }}(1-\frac{N_i}{N_{cri}})d_i{W}_i---\left(10\right)$

式中，I_lE为该地基的液化指数；d_i为第i层液化土层厚度(m)；W_i为第i层液化 土层单位土层深度的层位影响权函数值(m^-1)；N_i为第i层液化土层实测标准贯 入锤击数；N_cri为第i层液化土层液化判别标准贯入锤击数临界值。

但是，在现有规范中，尚未发现对于非平整场地或倾斜土层的液化判别标 准，因此进一步对倾斜场地的液化判别进行研究。

步骤三：倾斜场地的液化判别研究：在步骤二的基础上，对倾斜场地地基 模型进行模拟，逐渐改变地基倾斜角度，根据液化土层的液化程度随倾斜角度 的变化情况，给出不同倾角倾斜场地的抗液化指数与超孔隙水压比之间的对应 关系，并根据此对应关系对倾斜场地液化指数进行修正，达到对倾斜场地进行 液化判别的目的。本发明分别取地基倾斜角度为：2°、5°、8°、10°、15°进行判别 研究。具体步骤如下：

步骤301)：数值模拟条件与步骤二相同，将地基左侧抬高，以达到倾斜目 的；地基右侧地下水位埋深按0.5m考虑；

步骤302)：取I_N＝0.5时为例，对不同倾角的液化土层进行计算，记录不同 倾角下液化土层的超孔隙水压比r_u与抗液化指数I_N对应关系，如表3所示，并绘 出结果关系曲线，如图8所示；

表3不同倾斜角度和抗液化指数条件下的地基超孔隙水压比

从表3和图8可以看出，倾斜场地的液化判别标准与水平场地相比存在一 定的差异，“规范法”中的液化判别标准在倾斜场地处应进行修正。

步骤303)：取水平场地模型计算结果为基数，对比倾斜场地中液化土层超 孔隙水压比r_u(β)随倾角β的变化情况。如表4所示。

步骤304)：应用Matlab工具箱对r_u(β)与β进行关系曲线拟合，得到倾斜场 地中液化土层超孔隙水压比r_u(β)关于水平场地液化土层超孔隙水压比r_u(0)和倾 角β的函数关系式为

$r_u\left(\beta \right)=\left(\begin{array}{cc}{r}_u\left(0\right)& 0\le I_N\le 0.2\\ -0.027-0.0016\beta +1.019r_u\left(0\right)& 0\le I_N\le 1.0\end{array}\right)---\left(11\right)$

步骤305)：液化土层的超孔隙水压比r_u(β)和场地的液化指数I_lE都是对场地 模型液化程度的表征，认为两者在数值上成正比的关系，则可得倾斜场地液化 指数为：

$I_{lE}\left(\beta \right)=\left(\begin{array}{cc}{I}_{lE}\left(0\right)& 0\le I_N\le 0.2\\ -0.027-0.0016\beta +1.019r_u\left(0\right)& 0\le I_N\le 0.2\end{array}\right)---\left(12\right)$

式中，β为倾斜场地的倾斜角度，I_lE(β)为倾角为β的倾斜场地液化指数，I_lE(0)为 水平场地液化指数。

在实际工程中，应用本发明方法，可以得到不同倾角倾斜场地的抗液化指 数与超孔隙水压比之间的对应关系；进而通过对“规范”中水平场地液化判别 方法进行修正、达到对倾斜场地进行液化判别的目的，对“规范法”不适用于 倾斜场地的液化判别这一不足之处做出补充，提高了工程的安全性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通 技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变 形，这些改进和变形也应视为本发明的保护范围。