(1)
北京大学学报(自然科学版) 第62卷 第3期 2026年5月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 62, No. 3 (May 2026)
doi: 10.13209/j.0479-8023.2025.034
国家自然科学基金(42174133)资助
收稿日期: 2025–04–09;
修回日期:2025–05–01
摘要 针对最大有效力矩准则的作用范围缺乏定量解释的问题, 基于弹塑性梁模型开展研究。首先, 根据均匀弹塑性梁在纯弯曲状态下的正应力分布解析解, 推导出准则的单元边长近似等于梁的高度。然后, 采用COMSOL 有限元软件, 模拟双层和三层弹塑性梁的弯曲屈服过程, 结果表明上述理论关系对叠层梁(代表沉积地层)同样有效。最后, 通过分析梁上近似膝折带塑性变形区域的几何比例关系, 结合野外露头测量膝折带宽度, 定量地估算受该准则控制的地层厚度, 实现最大有效力矩准则中单元边长的定量化。研究结果可为相关地质构造的应力–应变场恢复和变形过程分析提供理论依据。
关键词 最大有效力矩准则; 弹塑性梁; 膝折带; 有限元方法; COMSOL
地壳中岩石因受到挤压而发生变形时, 上部脆性域内主要发生受库伦–摩尔准则控制的脆性破 裂[1], 破裂带的共轭角为锐角, 且角平分线指向最大主压应力(s1)方向。脆性域之下的韧性域内则局部发育塑性变形[1–2], 形成诸如膝折带[3]、伸展褶劈理[4]和低角度剪切带[5]等变形局部化的间隔性共轭剪切构造现象, 但它们指向最大主压应力(s1)方向的共轭角为稳定的 110°±10°钝角。为了解释这些与库伦–摩尔准则不相容的现象, 郑亚东等[5–10]提出控制变形局部化的最大有效力矩准则(maximum effec-tive moment criterion, MEMC)。
最大有效力矩准则认为, 由于地壳和岩石在现实中是非均匀的, 共轭韧性剪切带是局部塑性变形的产物, 变形带之间的区域不发生变形, 而是在力矩作用下发生旋转[10]。这种只旋转不变形的区域在其对应的观测尺度上都可以视为一个边长为 L的有限小方块单元, 因此作用在单元周缘产生变形带的切向力能够产生力矩效应, 并在变形带与 s1 的夹角a 等于 54.7°时取得最大值。由于最大有效力矩的方向夹角 54.7°及其倍数 109.4°为一不变量[10], 因此在实际应用中可以很方便地用于确定应力–应变场的方向, 其约 110°钝角的平分线指向最大主压应力方向的有效性得到大量野外观察[11–23]及实验结果[24–28]的证实。然而, 最大有效力矩准则在估算应力场大小方面的应用很少, 这是因为虽然该准则指出发生变形局部化时的差应力(s1–s3)等于材料的屈服强度 ss, 但材料的屈服强度随变形时的环境要素(温度、围压和应变速率)而改变[10,29]。显然, 对于某个具体的变形环境, 方块单元的边长 L决定了温度、围压和应变速率的变化范围, 是确定方块单元屈服强度(差应力大小)的关键参数。然而, 目前对单元边长 L的定量研究比较匮乏, 妨碍了对该准则控制的岩石变形过程进行定量分析。
弹塑性梁在弯曲过程中会产生约 110°钝角的塑性变形局部化区域[24], 符合最大有效力矩准则。因此, 本文以力学结构相对简单的弹塑性梁为基本模型框架, 基于理论推导和有限元数值模拟, 分析梁内由最大有效力矩准则控制的局部化塑性变形几何要素, 探究单元边长 L的定量关系以及应用于野外测量的方法。
韧性域内形成的膝折带、伸展褶劈理和低角度剪切带(图 1)等间隔性共轭剪切构造现象的共同特征是局部发育的塑性变形所夹钝角稳定为 110°± 10°, 角平分线指向最大主压应力(s1)方向。为解释这一现象, 郑亚东等[5–10]根据膝折带和共轭伸展褶劈理等都具有物质线旋转的基本特征, 分析构造带内任意有限小方块边界的应力状态(图 2), 获得使方块单元旋转的有效力矩 Meff [6]:
(1)
式中, H为力偶臂, L为单元边长, s和t分别为作用在单元边界上的正应力和剪应力, s1 为最大主(压)应力, s3 为最小主(压)应力。
最大有效力矩准则控制的是局部塑性变形, 因此单元周缘进入塑性变形状态的条件是差应力(s1–s3)达到材料的屈服强度ss。同时, 对于具体的共轭韧性剪切带, 方块单元的边长 L也是不变量, 因此当变形带与s1 的夹角a等于 54.7°时, Meff 取得最大值, 导致沿最大力矩方向能够有效地形成共轭韧性剪切带, 故称为最大有效力矩准则。在此基础上, 童亨茂等[30]认为式(1)是s1 方向与先存面理或层理平行时的一种特殊情况, 提出先存面理与s1 存在夹角𝛽时的修正公式:
(2)
然而, 最大有效力矩准则并没有对其定义的“有限小方块”的尺度进行约束, 即方块单元的边长L缺少定量关系, 因此无法确定准则控制的地层范围。
前人研究观察到钢管在热加工(例如 700℃)的弯曲变形过程中, 可能产生代表材料稳定性遭受破坏的 1%～4%的不均匀塑性变形[24], 符合最大有效力矩准则的特征: 变形局部化和标志性的约 110°钝角[25](图 3)。可以将钢管抽象为梁, 地层也往往被简化为平板状的梁。由此得到启发, 鉴于最大有效力矩准则在定义上具有普适性, 可以利用力学结构相对简单的梁模型对其开展深入的研究。同时, 最大有效力矩准则控制的地质构造是应力集中条件下产生的塑性变形带, 因此可以认为梁的材料性质总体上表现为弹塑性[31], 在应力达到弹塑性材料的屈服强度后才发生局部化的塑性变形。
图 4(a)用一个长度为 a, 高度为 b, 厚度为 d的细长梁(a/b>6)代表水平地层。当梁在两端弯矩 M作用下发生的变形为纯弯曲时, 具有以下特性[32]: 1)沿 x方向(梁轴线), 梁上的剪力为 0, 且弯矩保持恒定值 M; 2)梁的横截面(平行于 y方向且垂直于 xy平面)为平截面, 即原来垂直于梁中线的平面变形后, 仍然保持为平面并垂直于弯曲后的中线; 3)平截面上只有垂直于截面(即 x方向)的正应力 sx, 且平行于截面(即 y方向)的正应力 sy 和剪应力 txy 始终为 0。
(a)理想膝褶带的几何特征示意图(据文献[3]修改); (b)共轭伸展褶劈理(extensional crenulation cleavage, ECC)的角度关系(据文献[4]修改); (c)伸展褶劈理(点虚线)向上发育进入脆性域并最终形成低角度剪切带(拆离断层)的图解(据文献[5]修改)。MF: 糜棱岩前缘(mylonitic front), 为糜棱岩剪切组构发育的上界
图1 具有 110°钝角的间隔性共轭剪切构造
Fig. 1 Spaced conjugate shear structures with a 110° obtuse angle
图2 最大有效力矩图示(据文献[6]修改)
Fig. 2 Diagram of the maximum effective moment (modified after Ref. [6])
图3 数值模拟计算得到的局部化等效塑性应变随钢管弯曲程度的变化(据文献[25]修改)
Fig. 3 Numerical simulation results showing the variation of localized equivalent plastic strain with the steel tube bending (modified after Ref. [25])
如果图 4(a)中发生纯弯曲的梁由均匀的各向同性的理想弹塑性材料构成, 则当弯矩 M小于弹性极限弯矩时, 横截面上仅有的正应力 sx 的分布特征为从中间向两侧线性增大, 并在梁的上下边缘处达到最大值(图 4(c))。当弯矩 M等于弹性极限弯矩时, 梁边缘处的 sx 达到弹塑性材料的屈服极限, 由弹性变形转为理想塑性变形, 可以继续变形, 但不再能承担应力的增加(即维持 sx =ss 的状态不变), 这时抵抗进一步弯矩的能力由内部尚属弹性的材料来承担。此后, 进一步加大弯矩 M会导致 sx 从梁上下边缘向内部逐渐超过材料的屈服极限而由表及里地相继进入塑性状态, 此时梁横截面上的 sx 分布就会变成图 4(d)所示状态。
(a)理想弹塑性梁只在受弯矩 M 作用时发生纯弯曲, 箭头所示方向的 M 会导致梁向下弯, 因此梁的上半部分为压缩变形, 下半部分为拉伸变形; (b)梁上边缘中间位置((a)中虚线框)的塑性变形三角形微元示意图; (c)未发生屈服时的正应力(sx)分布; (d)从梁边缘向内部逐渐发生屈服时的正应力分布
图4 纯弯曲理想弹塑性梁的正应力示意图
Fig. 4 Schematic diagram of normal stress for an ideal elastoplastic beam under pure bending
由于纯弯曲梁上各处的弯矩相等, 因此上述塑性变形区域的出现和扩展在梁上各处是同时发生的。为不失一般性, 取塑性变形刚出现的梁中间上边缘的微元(图 4(a)中虚线框)进行分析。在放大微元的图 4(b)中, 令最大有效力矩准则控制的局部塑性变形带的长度为Ds, 则塑性小三角形的高Dy= Ds×sina。由于此时Dy 是一个较小的量, 所以其中心与矩心 O(在 sx 为 0 的中线上)的距离为 b/2, 则Dy上的正应力(图 4(d))绕 O的总力矩为
(3)
同时, 由最大有效力矩准则的计算公式(式(1))可知, 塑性小三角形斜边Ds上的剪切作用产生的总力矩应为
(4)
由于该塑性小三角形处于平衡状态, 而三角形底边(梁的上边缘)上的力矩为 0(临空面上的剪应力为 0), 因此 Ms 与 Me 应该数值相等但旋转方向相反。将式(3)和(4)联立, 得到单元边长 L的表达式为
(5)
由于最大有效力矩准则的 a固定为 54.7°, 且(s1–s3)就是材料塑性变形时的屈服强度 ss, 所以进一步得到 L的关系式:
(6)
即对于一个理想弹塑性材料的纯弯曲梁, 当因其上下边缘的应力达到屈服强度而形成约 55°的局部塑性变形带时, 最大有效力矩准则控制的块体单元边长 L的几何尺度等于梁的高度。
为验证上述单层弹塑性梁的理论分析结果对实际沉积地层(可视为叠层梁结构)的有效性, 本研究采用 COMSOL 有限元软件, 对双层和三层梁模型开展二维(平面应力)数值模拟实验。
如图 5(a)所示, 设置三层梁(弹塑性–纯弹性–弹塑性)模型的长、高和厚度(垂直纸面)分别为 a, b和d。同时, 设置 m和 n两个参数, 分别表示梁内弹塑性–纯弹性以及纯弹性–弹塑性两个层界面与梁下边缘的距离。显然, 当 n=0 时, 此三层模型自然转为一个双层模型(弹塑性–纯弹性)。通过逐渐增加施加在梁两端的弯矩 M, 使弹塑性梁向下弯曲的程度逐渐加大, 其中的纯弹性层只发生弹性形变, 可以模拟构造变形中不发生塑性变形的硬层(能干层), 弹塑性层则代表会发生塑性变形的软层(非能干层), 在应力达到屈服强度的部位由弹性变形转为塑性变形。
模型中的纯弹性和弹塑性材料并未设置成实际岩石, 而是参照构造物理模拟的思路(例如, 用黄油代替膏岩、岩盐和泥岩等, 用黏土代替碳酸盐岩和碎屑岩等[28]), 将杨氏模量和屈服强度按照一定的比例缩小, 加速数值模拟过程。其中, 纯弹性层较为简单, 设置杨氏模量 E=50MPa, 泊松比 ν=0.3; 弹塑性层则是在纯弹性层基础上增加屈服强度 ss (该参数的设置对有限元计算的收敛性影响很大, 需要特殊考虑)。由于本研究是通过梁两端的弯矩 M来控制梁的变形, 为避免边界效应的影响, 需要使得塑性变形发生在距离梁两端足够远的中部区域。因此, 设置弹塑性层的屈服强度从两端的 10MPa 向中间线性地递减, 在(5/12a, 7/12a)的中间区域形成一个屈服强度很低(0.1MPa)的易变形区(图 5(b))。同时, 前面进行理论推导时采用的理想弹塑性模型所产生的塑性区域会无限变形, 很难控制 M的增量幅度, 容易导致数值计算结果不收敛。因此, 本文采用线性强化弹塑性模型, 通过将材料的线性强化率E′设置为比杨氏模量 E小两个数量级(即 E′=E/100) (图 5(c))来近似地模拟理性弹塑性材料。
(a)梁模型的几何设置; (b)弹塑性层的屈服强度设置; (c)弹塑性层采用的线性强化塑性模型中的应力–应变关系
图5 COMSOL 中梁模型的参数设置
Fig. 5 Parameter settings for the beam model in COMSOL
根据上述几何模型和材料模型, 在 COMSOL有限元软件中创建二维平面应力条件下的双层和三层梁模型, 完成网格剖分和计算。经过多次尝试, 综合考虑收敛性、计算精度和耗时, 本研究的网格采用自由四边形单元, 具体参数列于表 1。
通过改变梁的长度(a)、整体高度(b)、弹塑性–弹性层界面高度(m)和弹性–弹塑性层界面高度(n), 可以计算得到一系列不同几何特征的梁随弯矩 M变化的形变位移以及多种应力–应变分布数据集。其中, 有效塑性应变(eep)可以直观地度量塑性变形的程度。
表1 网格剖分所使用的自由四边形单元参数
Table 1 Parameters of free quadrilateral elements used for meshing
参数名称数值 最大单元边长0.006 m 最小单元边长2.4×10−5 m 最大单元增长率1.1 曲率因子0.2 狭窄区域分辨率1
图 6(a)展示一个双层梁在 M=280N·m 时的塑性变形情况, 其中梁的长度 a=1.2m, 整体高度 b=0.1m, 弹塑性层和弹性层各占梁的一半高度(即弹塑性–弹性层界面高度 m=0.05m, 弹性–弹塑性层界面高度 n=0m)。可以看到, 梁上出现塑性变形的主要为弹塑性层的中间易变形区域(即图 5(b)中 x=0.5～0.7m 范围内), 且左右对称。需要注意的是, 虽然eep从梁的上边缘往内部整体上逐渐减小, 但在中间易变形区域的两端, 材料属性的“突变”造成应力集中, 导致该处的eep明显局部增大。
由于模型的对称性, 选取图 6(a)左侧的塑性应变局部化区域作为代表开展具体分析。通过对该区域进行放大(图 7), 可以从叠加的白色塑性应变等值线观察到类似膝折带的局部化结构, 其变形方向与梁边缘的夹角接近 55°。我们注意到, 计算得到的最大主压应力方向也平行于梁边缘(图 8(a)), 因此可以确认该变形局部化结构受到最大有效力矩准则的控制。同时, 该区域的 x方向正应力sx 几乎为一固定值(图 8(b)), 近似地满足理想塑性变形时应力不再增加的特征, 验证了将线性强化率 E′设置为E/100 的有效性。进一步地, 按照上述纯弯梁理论推导的方法, 在近似膝折带结构的边缘构建一个高为Dy, 斜边为Ds 的三角形区域(图 7 和 8)来计算最大有效力矩的单元边长 L。由式(3)和(4)可知, Ds 和Dy 上的力矩分别为
(7)
(8)
令
, 即可得到
(9)
式中, s1, s3和sx均为模拟结果, b'为Dy的矩心距。
(a)双层梁模型的等效塑性应变(eep)分布(左右对称); (b)三层梁模型的等效塑性应变分布(除左右对称(故只画出左半部分)外, 还上下对称)
图6 梁模型计算得到的等效塑性应变分布
Fig. 6 Equivalent plastic strain distribution computed from the beam models
q为膝折相关形迹之间的弦长(膝折带初始宽度), 红色三角形为计算 L 值的塑性变形区域, 白色等效塑性应变等值线的间距为 0.18×10−3
图7 中间易变形区域左侧的塑性应变局部化区域放大图
Fig. 7 Enlarged view of the localized plastic strain zone on the left side of the central weak deformation region
需要注意的是, 与理想弹塑性单层梁不同, Dy在叠层梁上的矩心(正应力sx =0 的点)往往不在梁的轴线上, 因此需要先找到数值计算结果中正应力为零的点的坐标, 再计算其与Dy的距离b'。
(a)最大主压应力方向(蓝色箭头)平行于梁边缘; (b)正应力 sx 的分布特征, 红色等值线的间距为 0.01MPa
图8 最大主压应力方向与正应力大小的特征
Fig. 8 Characteristics of the maximum principal compressive stress direction and the normal stress magnitude
固定 a=1.2m, b=0.1m, 改变 m值可获得一系列层厚度不同的双层梁模型来进行数值模拟, 并计算相应的 L值。图 9 显示, 随弯矩 M和层厚度的变化, L在其平均值(约 0.101m)附近很小的范围内波动, 可以认为等于双层梁的高度 0.1m。
固定 a=1.2 m, b=0.1 m; 图例中的白色部分为弹性层, 下同
图9 双层梁的 L计算值随弯矩和弹性层厚度的变化
Fig. 9 Variation of calculated L values of the double-layer beam with bending moment and the elastic layer thickness
图 6(b)展示一个三层梁模型的模拟结果, 梁的长度和高度与图 6(a)中双层梁相同, 顶部和底部弹塑性层厚度均为 0.03m, 中部弹性层厚度为 0.04m, 弯矩 M=260N·m。可见, 模型顶部和底部弹塑性层的等效塑性应变相对于梁的中线上下对称, 并且具有与双层模型类似的分布特征, 即顶层和底层都出现近似膝折带的局部变形区域。因此, 采取与双层梁模型相同的方法来分析和计算 L值。图 10 展示在保持三层梁模型中的弹性层厚度(m−n=0.04m)不变时, L值随 M值和弹性层高度(改变顶部和底部弹塑性层的厚度)的变化情况, 可见 L值的变化范围很小, 平均值约为 0.091m, 明显稍小于双层梁模型的模拟结果, 体现中间“硬层”的影响。进一步的模拟分析表明, 在一定的范围内, L值随中间弹性层厚度增大(不超过梁高度的一半)而大致呈线性减小趋势(图 11)。
固定 a=1.2 m, b=0.1 m; 保持中间弹性层的厚度不变(m−n=0.04 m)
图10 三层梁的 L计算值随弹性层位置的变化
Fig. 10 Variation of calculated L values of the triple-layer beam with the intermediate elastic layer position
图11 三层梁的 L 平均值随中间弹性层厚度的变化
Fig. 11 Variation of the averaged L values of the triple-layer beam with the intermediate elastic layer thickness
综合上述双层梁和三层梁的数值模拟结果, 可知最大有效力矩准则控制的单元边长 L与梁的长度a无关, 只与梁的高度 b相关, 双层梁的 L值几乎等于 b, 三层梁的 L值则会因为中间硬层的影响而略小于 b。因此, 可以推论, 在相当大的变形尺度范围内, 当一套沉积地层可以近似为细长梁时, 其实际地层厚度(梁的高度)在几何尺度上即为最大有效力矩准则控制的单元边长 L。
叠层梁的模拟结果表明, 最大有效力矩准则控制的单元边长 L并不是发生塑性变形的“软层”厚度, 而是包括“硬层”的整条梁的高度。那么, 在野外如何通过对“软层”中类似膝折带的变形局部化区域的观测来较准确地估计实际梁结构中包含的地层及其厚度呢？
对于实际的地质变形, 最可能直接测量到的是两个构造形迹之间的距离。如果可以将该变形地层简化为一个弯曲的梁, 则构造形迹之间的距离为弯曲梁的一段弦长 q(图 12)。特别地, 对于最大有效力矩准则控制的变形局部化构造(例如膝折带), 由于变形带所夹钝角为约 110°的不变量, 与材料无关, 因此 q对应的圆心角也是固定的, 因此 q与梁中性面的曲率半径 R具有某种固定的比例关系, 可记为 q=kR。
图12 纯弯曲梁的中性面曲率半径(R)与某段变形之间的距离(弦长)(q)示意图
Fig. 12 Schematic diagramof the neutral surface radius of curvature (R) and the chord length (q) of a deformed segment in a pure bending beam
当图 12 为纯弯梁时, 其中性面的曲率半径 R与正应力sx有如下关系[32]:
(10)
由于塑性变形首先发生在梁的边缘(y=b/2), 为一自由面, 其上的剪应力为 0, 因此有 sx =s1=ss, 从而式(10)变为
(11)
这说明在材料性质(杨氏模量和屈服强度)确定的情况下, 曲率半径 R与梁的高度 b成正比。
由于 q=kR与材料无关, 因此可以通过数值模拟结果确定 k值。如图 7 所示, 形迹 1 代表近似膝折带变形一侧的边界, 形迹 2 代表材料分界处(应力集中)的另一侧边界, 因此可以将形迹 1 和形迹 2 分别与梁上边缘的交点之间的距离 q视为两个变形构造之间的弦长。在双层梁模型中, 保持其他参数不变, 通过更改高度 b, 可以模拟得到一系列等效塑性应变, 并通过测量得到相应的 q值, 且由式(11)计算得到对应的 R值。图 13 显示, q与 R 有非常好的线性关系, 与理论推导结果一致, 并得到 k=3.94×10−4。
如果在野外可以测量受最大有效力矩准则控制的构造变形的横向宽度 q, 就可以利用 q=kR计算得到对应梁结构的曲率半径 R, 再在式(11)中代入实际岩石的杨氏模量和屈服强度, 最终计算得到梁结构的高度 b, 即最大有效力矩准则控制的单元边长L。需要指出的是, 由于本研究主要模拟梁刚发生塑性变形时的情况, 梁的弯曲程度其实很小(变形程度小), 因此本文中的 k值适用于膝折带刚开始形成时的最短宽度; 在实际构造中, 随着膝折带从平缓变为陡立, 有一个宽度加大的过程(图 14(a)～(c)), 即野外实际测量的膝折带宽度往往大于初始变形时的最短宽度。下面以新疆巴楚地区大阪塔格山南段褶皱的野外剖面相关数据[28]为例进行说明。
图13 不同高度双层梁的塑性变形局部化宽度(q)与该区域对应的曲率半径(R)的关系
Fig. 13 Relationship between the localized plastic deformation width (q) and the corresponding radius of curvature (R) in the double-layer beam with varying heights
巴楚地区的地层大致可以分为 3 层, 中间层是发育膝折褶皱的奥陶系、志留系和上新统的碳酸盐岩、砂岩、泥岩和碎屑岩等, 厚约 1500～3500ｍ; 上覆新近系～第四系未成岩的松散黄色细砂, 厚约几百至 2000m; 下伏寒武系膏岩、岩盐和泥岩等能干性较弱的岩石, 厚二百至近千米, 成为广泛分布且稳定的主要滑脱层[28]。可见, 由于上下层的岩石强度相对于中间层小很多, 可将发育膝折褶皱的中间地层视为多层的梁结构。闫淑玉等[28]对巴楚隆起西北端的大阪塔格山南段褶皱进行测绘, 并根据箱状褶皱冲起后形成的构造样式(图 14(e)), 认为该褶皱是膝折型褶皱的剥蚀残留(图 14(f))。
(a)水平 s1 压缩作用下, 初始宽度为 q 的膝折带开始形成, b= 35.3°, 符合最大有效力矩准则; (b)膝折带边界固定, 但带内的地层发生旋转, 当 a=b时, 带内地层与边界垂直, 此过程中膝折带宽度保持为 q 不变; (c)在保持a=b的条件下, 膝折带边界发生迁移, 原来水平的地层被纳入膝折带后发生旋转, 使得膝折带宽度增大; (d)边界迁移结束后(比如抵达新的材料边界), 膝折带内地层恢复旋转, 直到 a=2b时地层陡立(约70°); (e)后续构造缩短变形不再通过膝折带的边界扩展和地层旋转来吸收, 需要通过断层来调节, 形成箱型脱顶褶皱(box-shaped lift-off fold); (f)新疆巴楚地区大阪塔格山南段膝折褶皱的野外露头解析, 虚线为恢复的被剥蚀部分
图14 新疆巴楚地区大阪塔格山南段膝折褶皱的解析示意图(据文献[28]修改)
Fig. 14 Analytical diagram of the kink fold in the southern section of Dabantage Mountain, Bachu area, Xinjiang (modified after Ref. [28])
如图 14(f)所示, 野外露头上膝折带的宽度 q'约为 330m。由于该膝折褶皱的核部地层陡立, 对应膝折带发育的第三阶段(图 14(d))。对比图 14(a)与(d), 可假设膝折带刚发育时的最短宽度 q (图 14(a))大致为野外观察值 q'的一半, 即 165m 左右, 则可计算得到梁结构的曲率半径:
考虑到自然界中碳酸盐岩的屈服强度变化范围较大, 简单起见, 本文用膝折褶皱中的砂岩来估算梁结构的材料参数。砂岩的杨氏模量 E可取值为 30GPa, 因其弹性极限基本上等于抗压极限, 往往近似地将其抗压强度(通常可取值为 75MPa[32])作为无围压时的屈服强度。但是, 有围压时的岩石屈服强度或弹性极限会明显提高(即韧性增强), 例如, 实验测得围压达到 75MPa 时(对应约 3km 的埋深), 抗压强度最大可以提升到约两倍[32]。膝折褶皱中的砂岩显然是在有围压的条件下发生变形, 因此不妨将梁结构的屈服强度定为最大值 150MPa。将 R= 4.2×104 m, ss=150×106 Pa 和 E=3×1010Pa 代入式(11), 计算得到梁结构的高度 b≈4200m。根据穿过巴楚隆起西段的东北向地震深度剖面资料[28], 在大阪塔格山南段区域, 较完整的中间地层厚度可达 3000～ 4000m, 可见 4200m 的梁高度估算值是合理的。如果能在实验室测到当地碳酸盐岩的屈服强度, 由于碳酸盐岩相对于砂岩偏软, 因此可以推测估算的 b值会更符合中间地层的实际厚度。
本文基于理论分析并结合有限元数值模拟, 利用弹塑性梁模型, 对最大有效力矩准则控制的单元边长 L进行定量研究, 得到如下结论。
1)对于均匀的理想弹塑性梁, 可由纯弯曲条件下的正应力分布解析解, 推导得到单元边长 L约等于梁的高度 b。
2)对于双层和三层弹塑性梁, 则通过有限元数值方法模拟弯曲过程, 得到符合最大有效力矩准则角度关系的近似膝折带的塑性变形区域, 并分析得到的单元边长 L仍然对应梁的高度。这一几何定量关系对野外实际地层也是有效的。
3)对于类似叠层梁的实际地层, 通过在野外露头测量膝折带宽度, 可以估计受最大有效力矩准则控制的地层总厚度。
参考文献
[1] Ramsay J G, Carreras J, Cobbold P R, et al. Shear zone geometry: a review. Journal of Structural Geology, 1980, 2(1/2): 83–99
[2] 侯泉林, 程南南, 石梦岩, 等. 不同构造层次岩石变形准则的融合与发展. 岩石学报, 2018, 34(6): 1792–1800
[3] 张波, 李生福, 张进江, 等. 膝褶、膝褶带、共轭膝褶带——一种可能的新型油气构造样式. 天然气工业, 2010, 30(2): 32–39, 71, 136–137
[4] 郑亚东. 共轭伸展褶劈理夹角的定量解析. 地学前缘, 1999, 6(4): 391–395
[5] Zheng Y, Wang T, Ma M, et al. Maximum effective moment criterion and the origin of low-angle normal faults. Journal of Structural Geology, 2004, 26(2): 271–285
[6] 郑亚东, 王涛, 王新社. 新世纪构造地质学与力学的新理论——最大有效力矩准则. 自然科学进展, 2005, 15(2): 16–22
[7] 郑亚东, 王涛, 王新社. 最大有效力矩准则的理论与实践. 北京大学学报(自然科学版), 2007, 43(2): 145–156
[8] 郑亚东, 王涛, 王新社. 最大有效力矩准则及相关地质构造. 地学前缘, 2007, 14(4): 49–60
[9] 郑亚东, 张进江, 王涛. 实践中发展的最大有效力矩准则. 地质力学学报, 2009, 15(3): 209–217
[10] 郑亚东, 张进江, 张波. 新世纪构造地质学两大支柱理论: 最大有效力矩准则与变位形分解. 地质力学学报, 2022, 28(3): 319–337
[11] Park R G. Shear-zone deformation and bulk strain in granite-greenstone terrain of the western Superior Pro-vince, Canada. Precambrian Research, 1981, 14(1): 31–47
[12] Spencer J E. Role of tectonic denudation in warping and uplift of low-angle normal faults. Geology, 1984, 12(2): 95–98
[13] Sibson R H, Robert F, Poulsen K H. High-angle reverse faults, fluid-pressure cycling, and mesothermal gold-quartz deposits. Geology, 1988, 16(6): 551–555
[14] Malavieille J. Late orogenic extension in mountain belts: insights from the Basin and Range and the late Paleozoic Variscan belt. Tectonics, 1993, 12(5): 1115–1130
[15] John B E, Foster D A. Structural and thermal con-straints on the initiation angle of detachment faulting in the southern Basin and Range: the Chemehuevi Mountains case study. Geological Society of America Bulletin, 1993, 105(8): 1091–1108
[16] Ujiie K, Maltman A J, Sánchez-Gómez M. Origin of deformation bands in argillaceous sediments at the toe of the Nankai accretionary prism, southwest Japan. Journal of Structural Geology, 2004, 26(2): 221–231
[17] Camerlo R H, Benson E F. Geometric and seismic interpretation of the Perdido fold belt: northwestern deep-water Gulf of Mexico. AAPG Bulletin, 2006, 90 (3): 363–386
[18] 莫午零, 郑亚东, 张文涛, 等. 柴达木盆地油泉子储油构造分析. 石油与天然气地质, 2007, 28(3): 324–328
[19] Zhang H, Cheng S, Wang G, et al. Tectonic inversion in sediment-hosted copper deposits: the Luangu area, West Congo Basin, Republic of the Congo. Minerals (Basel), 2024, 14(11): 1061
[20] 陈鹏, 施炜. 南秦岭造山带韧性剪切系中–晚侏罗世运动学分析与力学机制探讨. 地球科学进展, 2015, 30(1): 69–77
[21] Pei S, Zhang H, Su J, et al. Ductile gap between the Wenchuan and Lushan earthquakes revealed from the two-dimensional Pg seismic tomography. Scientific Reports, 2014, 4(1): 6489
[22] 张沛全, 李冰溯. 最大有效力矩准则约束下的广西罗阳山地震、地质与地貌效应. 地质力学学报, 2012, 18(1): 79–90
[23] 张波, 张仲培, 张进江, 等. 盆地腹地隆起带内高角度变形带构造样式新解. 北京大学学报(自然科学版), 2012, 48(1): 79–91
[24] Hallai J F, Kyriakides S. On the effect of Lüders bands on the bending of steel tubes. Part I: Experiments. In-ternational Journal of Solids and Structures, 2011, 48 (24): 3275–3284
[25] Hallai J F, Kyriakides S. On the effect of Lüders bands on the bending of steel tubes. Part II: Analysis. Inter-national Journal of Solids and Structures, 2011, 48(24): 3285–3298
[26] Vernooij M G C, Kunze K, den Brok B. ‘Brittle’ shear zones in experimentally deformed quartz single cry-stals. Journal of Structural Geology, 2006, 28(7): 1292–1306
[27] 申发龙, 张进江. 单轴压缩剪应变局部化带分布的数值模拟及其对最大有效力矩准则的验证. 北京大学学报(自然科学版), 2013, 49(5): 873–879
[28] 闫淑玉, 张进江, 张波, 等. 新疆巴楚地区共轭膝折带的物理模拟研究. 大地构造与成矿学, 2011, 35 (1): 24–31
[29] 郑亚东. 结构面力学性质的定量鉴定. 地质力学学报, 2005, 11(3): 197–203
[30] 童亨茂, 王明阳, 郝化武, 等. 最大有效力矩准则的理论拓展. 地质力学学报, 2011, 17(4): 312–321
[31] 郑亚东, 张进江, 王涛, 等. 控制变形局部化的最大有效力矩准则//2016 中国地球科学联合学术年会. 北京, 2016: 2114–2116
[32] 王仁, 丁中一, 殷有泉. 固体力学基础. 北京: 地质出版社, 1979
Quantitative Study on Unit Length of Maximum Effective Moment Criterion Based on Elastoplastic Beam Model
Abstract Aiming at the lack of quantitative interpretation of the operational scale of the maximum effective moment criterion, this study is carried out based on an elastoplastic beam model. First, using the analytical solution for the normal stress distribution in a homogeneous elastoplastic beam under pure bending, the unit length of the criterion is derived to be approximately equal to the beam height. Subsequently, the bending and yielding process of double- and triple-layer elastoplastic beams are simulated using COMSOL multiphysics software. The results confirm that the above theoretical relationship is also valid for laminated beams (representing sedimentary strata). Finally, by analyzing the geometric scaling relationship of the kink band‑like plastic deformation zones in beams and com-bining with the measurement of kink band width in field outcrops, the thickness of strata governed by this cri- terion is quantitatively estimated, thereby achieving the quantification of the unit length in the maximum effective moment criterion. The results provide a theoretical basis for reconstructing stress‑strain fields and analyzing defor-mation processes in related geological structures.
Key words maximum effective moment criterion; elastoplastic beam; kink band; finite element method; COMSOL