接触面损伤演化过程的数值模型

拟。局部坐标 系 中，节 理 单 元 的 应 力－节 点 位 移 关 系为
｛σ｝＝
烅烄烆στｎｓ烍烌烎＝
ｈ１［λ］烅烄烆ΔΔυｕ烍烌烎＝
１ ｈ
［λ］［Ｍ
］｛δ｝ｅ
（１）
第１期
唐 世 斌 ，等 ：接 触 面 损 伤 演 化 过 程 的 数 值 模 型
１４７
式中 γｓ 和εｎ 分 别 为 接 触 单 元 的 剪 应 变 和 法 向 应 变，τｓ 和σｎ 分别为切 向 剪 应 力 和 法 向 正 应 力，｛δ｝ｅ 是单元节点 位 移，［λ］为 刚 度 系 数 矩 阵，由 下 式 确 定：
（大连理工大学 岩石破裂与失稳研究所，大连 １１６０２４）
摘 要：结合接触面在细观尺度上的非均匀性及其损伤演化规律，提出了一种可以 考 虑 接 触 面 损 伤 演 化 过 程 的 数 值模型。数值模拟得到的剪切应力－剪切 位 移 关 系 表 明：随 着 剪 切 位 移 的 增 加，接 触 面 上 产 生 的 损 伤 单 元 导 致 曲 线斜率逐渐降低；当接触面被完全剪断后，模型中的剪 切 应 力 出 现 了 一 定 程 度 的 下 降，此 后 的 剪 应 力 在 摩 擦 力 的 作用下基本保持不变。模型所能承受的最大剪应力随着垂直应力的增加而增大。整个损伤过程揭示了非均匀性 的存在是非线性变形的根本，而高的剪应力是产生剪切损伤的根源。损伤的产生 使 得 应 力 在 损 伤 点 得 以 释 放 ，并 转移到其前端未损伤的接触单元上。这种损伤－应力调整的循环 规 律 贯 穿 整 个 变 形 过 程 ，直 至 模 型 中 的 接 触 面 被 完全剪断。
应的损伤变量 Ｄｓ 表示为
烄０，
Δｕ ＜ Δｕτ０
Ｄｓ（Δｕ）＝ 烅烆１－λｈｓτ０Δｓｕｕ， Δｕ ≥ Δｕτ０
（１０）
式中 τｓｕ 表示剪切损伤残余强度，Δｕτ０ 是剪切弹 性
极限对应的水平相对位移。
如 果 单 元 产 生 剪 切 损 伤 ，则 把 单 元 标 记 为 滑 移
［ ］ ［Ｍ］＝
１ ２
－Ｚ１ ０ －Ｚ２ ０ Ｚ２ ０ Ｚ１ ０ ０ －Ｚ１ ０ －Ｚ２ ０ Ｚ２ ０ Ｚ１
（３）
式中 Ｚ１ ＝１－２ｌｘ′， Ｚ２ ＝１＋２ｌｘ′。 根 据 虚 功 原 理 ，由 上 述 方 程 转 换 而 成 的 有 限 元
方程为
∫ ［Ｆ］ｅ＝
ｔ ｈ
ｌ／２
［Ｍ］Ｔ［λ］［Ｍ］ｄｘ′｛δ｝ｅ＝［Ｋ］ｅ｛δ｝ｅ （４）
收 稿 日 期 ：２００９－０２－２２；修 改 稿 收 到 日 期 ：２０１０－０６－１７． 基 金 项 目 ：国 家 重 点 基 础 研 究 发 展 计 划 （２００７ＣＢ２０９４０４）；
国家自然科学基金（５１００４０２０，５０８０４００６，５０７７８０４６， ５０８０９０１２）资助项目 ． 作者简介：唐世斌＊ （１９８０－），男，博士后 （Ｅ－ｍａｉｌ：ｔａｎｇ＿ｓｈｉｂｉｎ＠１２６．ｃｏｍ）．
（ｂ）单 元 法 向 处 于 压 缩 状 态 研究表明 ：当 ［１０－１２］ 接 触 面 处 于 压 缩 状 态 时 ，其
法向和切向刚度系数都随着法向应力的变化而变
１４８
计算力学学报
第 ２８ 卷
化。为了反映刚度系 数 随 应 力 的 变 化 关 系，本 数 值
模型采用如下形式 ： ［１０］
关 键 词 ：损 伤 演 化 ；接 触 单 元 ；破 坏 过 程 ；数 值 模 拟 中 图 分 类 号 ：ＴＵ４４３ 文 献 标 志 码 ：Ａ
１ 引 言
岩体结构、混 凝 土 浇 注 层 之 间 的 粘 结 带、钢 筋 与 混 凝 土 之 间 的 过 渡 层 、建 筑 结 构 与 地 基 之 间 的 接 触面以及土与柱桩之间的粘结层等都对结构的稳 定性起着重要影响。这些接触面在荷载作用下可 能会出现局部脱 开、滑 动、错 位 及 张 闭 等 非 连 续 变 化现象，产生不同 于 接 触 面 两 边 材 料 的 力 学 响 应， 按照实体结构建立的理论方法不能很好地分析接 触面的力学行为。考虑到岩体中存在的节理特性， Ｇｏｏｄｍａｎ［１］于 １９６８ 年 提 出 了 线 性 的 无 厚 度 节 理 单 元 ，从 而 奠 定 了 接 触 面 数 值 计 算 的 研 究 基 础 。 为 了 防止 Ｇｏｏｄｍａｎ单元 的 嵌 入 行 为，李 守 德 等 ［２］ 采 用 如 下 方 法 进 行 处 理 ，即 单 元 受 压 时 将 法 向 刚 度 取 一 大值，受 拉 时 则 取 一 小 值。Ｚｉｅｎｋｉｅｗｉｃｚ［３］于 １９７０ 年 提 出 了 二 维 等 参 节 理 单 元 模 式 ，并 将 其 推 广 到 三 维问题 的 研 究 中。Ｄｅｓａｉ［４］提 出 了 一 种 薄 层 单 元， 在 单 元 中 采 用 了 防 止 单 元 嵌 入 的 控 制 方 法 ，在 一 定 程度上克 服 了 Ｇｏｏｄｍａｎ 单 元 容 易 发 生 嵌 入 的 缺 点 ，但 在 计 算 中 涉 及 到 的 单 元 厚 度 难 以 确 定 以 及 没 有反 映 法 向 和 切 向 的 耦 合 影 响 ，是Ｄｅｓａｉ单 元 的 不
上述成果为接触面力学行为的研究提供了理 论依据。但这些理论模型是建立在接触面为均匀 材 料 假 设 基 础 上 的 ，在 反 映 其 本 质 力 学 特 性 方 面 有 所 欠 缺 。 本 文 拟 在 前 人 研 究 结 果 的 基 础 上 ，结 合 接 触 面 的 非 均 匀 性 特 性 ，从 损 伤 力 学 的 角 度 建 立 模 拟 接触面的破裂过 程 以 及 拉 开、滑 动 的 数 值 模 型，以 期为更深入研究接触面的力学特性提供参考依据。
的非均匀性，本文 假 设 单 元 材 料 性 质 服 从 Ｗｅｉｂｕｌｌ
分布，按照蒙特 卡 罗 （Ｍｏｎｔｅ Ｃａｒｌｏ）方 法 生 成 各 单
元的参数值。假设模型中有 Ｎ 个接触单元，则可按
照 Ｗｅｉｂｕｌｌ分布函数生成 Ｎ 个 单 元 的 切 向 和 法 向
刚度系数以及相应的强度参数。
２．２ 接触面单元损伤演化以及滑动模拟 按 照 损 伤 程 度 ，接 触 面 的 运 动 状 态 可 以 分 为 两
（１－Ｄ）σ～
＝
烄（１－Ｄｓ（Δｕ））λｓ０Δｕ／ｈ 烅烆（１－Ｄｎ（Δｖ））λｎ０Δｖ／ｈ
（６）
式中 Ｄｓ（Δｕ）是 切 向 刚 度 损 伤 变 量，Ｄｎ（Δｖ）是 法 向刚度损伤变量，λｓ０ 和λｎ０ 分别为切向和 法 向 的 初 始刚度系数。
已有研究表 明：虽 然 岩 石、混 凝 土 等 非 均 匀 材
第２８卷 第１期 ２０１１ 年 ２ 月
计算力学学报 Ｃｈｉｎｅｓｅ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｍｅｃｈａｎｉｃｓ
Ｖｏｌ．２８，Ｎｏ．１ Ｆｅｂｒｕａｒｙ ２０１１
文 章 编 号 ：１００７－４７０８（２０１１）０１－０１４６－０６
接触面损伤演化过程的数值模型
唐世斌＊ ， 唐春安， 梁正召
［ ］ 烄λｎ
＝λｎ０
１－
σｎ ΔＶλｎ０ ＋σｎ
－２
（ ） 烅
λｓ 烆
＝λｓ０
１－σｎＲｔａｆτｎｓδ
２
（９）
式中 ΔＶｕ 为法向最大压缩量，Ｒｆ 为破坏比，δ为内
摩擦角。
（ｃ）单 元 处 于 剪 切 状 态
当接触面单元的剪 应 力τｓ 达 到 由 莫 尔 －库 仑 准 则确定的剪切破坏阀值，即｜τｓ｜≥τｓｕ ＝Ｃ－ｆσｎ 时 ，单 元 产 生 剪 切 损 伤 ，从 而 导 致 切 向 刚 度 弱 化 。相
Ｃ－ｆσｎ 时，接触面产生剪切损伤。无论拉伸损伤还 是剪切损伤，单 元 在 损 伤 后 仍 具 有 一 定 的 残 余 强

度 。在 本 文 中 ，接 触 单 元 的 应 力 以 拉 升 为 正 ，压 缩 为
负 ；位 移 差 以 拉 长 为 正 ，缩 短 为 负 。以 下 将 按 照 各 种
应力状态分别阐述接触单元的损伤原理。
－ｌ／２
式中 ［Ｋ］ｅ 是接 触 单 元 在 局 部 坐 标 系 下 的 刚 度 矩
阵 ，转 换 为 整 体 坐 标 系 的 刚 度 矩 阵 为
［Ｋ］Ｇ ＝ ［Ｔ］－１［Ｋ］ｅ［Ｔ］
（５）
式中 ［Ｋ］Ｇ 是全局坐标系下的刚度矩阵，［Ｔ］为与
节 理 倾 角θ 相 关 的 坐 标 转 换 矩 阵 。
为 了 模 拟 节 理 法 向 和 切 向 刚 度 系 数 ，以 及 强 度
（ａ）单 元 法 向 处 于 拉 伸 状 态
当σｎ ＜σｎｕ 时，法向应力与法向位移差之间 保
持线性关系，而当σｎ ≥σｎｕ 时，单元产生拉伸损 伤， 从而导致法向和切向刚度同时弱化。
法向刚度损伤变量 Ｄｎ 可描述为
烄０ Ｄｎ（Δｖ）＝ 烅１－λｈｎσ０Δｎｒｖ
Δｖ ＜ Δｖｎ０ Δｖｎ０ ≤ Δｖ ≤ Δｖｎｕ
则 需 判 断 单 元 是 否 产 生 了 嵌 入 和 滑 动 ，并 进 行 相 应 的 刚 度 调 整 。具 体 处 理 方 法 如 下 。
２．２．１ 未开裂接触单元损伤演化过程 由 Ｌｅｍａｉｔｒｅ损 伤 原 理 可 知，未 开 裂 接 触 单 元
的损伤本构关系为
｛σ｝＝烅烄烆στｎｓ烍烌烎＝
２ 接触面损伤演化模型
２．１ 接 触 面 单 元 本 构 关 系 的 有 限 元 模 型 利 用 有 限 元 对 接 触 面 进 行 数 值 模 拟 时 ，通 常 采
用有厚度和无厚度的接触单元。有厚度接触单元 能 克 服 无 厚 度 单 元 的 嵌 入 问 题 ，且 能 很 好 地 分 析 接
触 面 厚 度 对 工 程 结 构 的 影 响 。 因 此 ，本 文 选 用 有 厚 度、带强度 的 接 触 单 元 进 行 接 触 面 力 学 特 性 的 模
［ ］ ［λ］＝ λｓｓ λｓｎ λｎｓ λｎｎ
（２）
式中 λｓｓ 为切向刚度系 数，λｎｎ 为 法 向 刚 度 系 数，λｓｎ
和λｎｓ 为考虑切向和 法 向 耦 合 效 应 的 刚 度 系 数。由
于目前很难测定 耦 合 系 数 项 的 数 值，一 般 假 定λｓｎ
和λｎｓ 为零。
系 数 矩 阵 ［Ｍ ］由 下 式 确 定 ：
料 在 宏 观 上 表 现 为 非 线 性 ，但 是 在 细 观 层 次 上 却 可
以用具有残余强度的弹脆性本构关系描述［９］。本 数
值 模 型 中 认 为 接 触 面 的 损 伤 也 具 有 同 样 的 性 质 ：当
接触面的法向应 力 达 到 其 拉 伸 应 力 阀 值，即σｎ ≥
［σ］＝σｔ 时，单 元 产 生 拉 伸 损 伤；当 剪 应 力 达 到 摩 尔 －库伦准则所规定的剪切强度，即｜τｓ｜≥ ［τ］＝
理。
拉 伸 损 伤 状 态 下 ，认 为 切 向 刚 度 具 有 与 法 向 刚
度 相 同 的 损 伤 量 ，即
Ｄｓ（Δｕ）＝ Ｄｎ（Δｖ）
（８）
由式（８）可以看出，若 Δｖｎｕ ＝ Δｖｎ０，则 为 理 想 的弹脆性模型；若σｎｒ ＝σｎｕ ，则为理想弹塑性模型。
因 此，参数 Δｖｎｕ 和σｎｕ 使模型在理想弹塑性与弹脆 性之间具有更好的灵活性。
足 。 ［５］ 卢廷浩 等 ［６］ 结合接触面单剪实验，提 出 了 接 触 面 薄 层 单 元 耦 合 本 构 模 型 ，从 而 使 计 算 结 果 更 符 合实际。此后，胡黎明 根 ［７］ 据试验结果建立了 描 述 接触面应变软化和剪胀现象的损伤力学模型。任 懿 等 ［８］ 通过时域展开技术与一个简单的等效 假 定， 得到了递推格式的粘弹性单向节理岩体的等效本 构关系。
（７）
烆１
Δｖ ＞ Δｖｎｕ
式中 σｎｒ 表示拉伸损伤残余强度，Δｖｎ０ 是拉伸 弹 性
极限对应 的 垂 直 相 对 位 移，Δｖｎｕ 是 最 大 相 对 垂 直 位移，即当单元的垂直位移达到 Δｖｎｕ 时，单 元 完 全 失去承载能力，单元 真 正 被 拉 开，并 在 数 值 模 型 中
将此单元标记 为 分 离 状 态，便 于 后 续 接 触 判 断 处
个 阶 段 ，即 破 坏 前 的 损 伤 积 累 阶 段 和 破 坏 后 的 拉 开 或滑移阶段。因此，本 文 将 接 触 单 元 分 为 未 开 裂 和 已开裂两种单元：对 于 未 开 裂 接 触 单 元，在 荷 载 作
用过程中将按照带拉伸的库伦摩尔准则进行破坏
判断并做相应的损 伤 处 理；对 于 已 开 裂 接 触 单 元，