0 引言

石峁遗址，位于陕西省榆林市神木市高家堡镇石峁村的秃尾河北侧山峁上，地处陕北黄土高原北部边缘，黄河一级支流秃尾河从其西南侧穿过^[1]。石峁遗址距今约四千年左右，是探寻中华文明起源的窗口，为中国文明起源形成的多元性和发展过程提供了全新的研究资料。如图1所示，石峁遗址主要包括皇城台、门址、角台、城墙，从位置上来看，主要包括外城、内城和皇城台^[2]。其中外城为不规则弧形城墙，与内城东南部墙体共同构造形成了相对独立的区域；内城以“皇城台”为中心依山势而建，形成一个不规则椭圆形空间；“皇城台”则是一处四面包砌护坡的台城，是当时石峁城址的中心^[3]。

在几千年的岁月变迁中石峁城墙遭受了大量诸如裂缝、倾斜、坍塌、泛盐、风化剥蚀等病害，特别是外城东门墩台墙体已出现严重倾斜现象，迫切地需要对其进行加固保护。

目前针对石峁遗址研究，主要集中在考古^[1-2]、历史^[3]、出土文物及城建技术^[4-7]等方面，对石峁城墙稳定性研究还未完全展开。而对其他城墙如西安明城墙的稳定性研究中，大多数学者通过建立有限元模型进行稳定性^[8-9]及动力分析^[10]，以上研究认为砌体是一种连续性的材料，这与实际不相符合，且城墙本质为带有包砖的夯土墙体，这与石峁城墙的组成结构有较大差异。与石峁城墙研究相近的，是对石砌挡墙研究，主要也是利用有限元对其稳定性进行了分析^[11-13]。有限元将石砌挡墙看着均值材料，无法考虑灰缝强度及其开裂过程。离散单元法则可以模拟不同块体的相互作用，块体可以发生滑动、张开，甚至脱离，更接近墙体内部的实际情况，能真实地表达砌体结构中的裂缝扩展过程^[14]。

综上，本研究以外城东门墩台城墙为研究对象，验证了离散单元法建模的有效性，基于离散单元法对城墙的稳定性进行评估，研究了风化对稳定性的影响，在此基础上提出了加固方法。

1 外城东门

石峁遗址外城东门位于外城东北部，是遗址东北部地区地势最高处^[4]。外城东门视野开阔，位置险要，是一座规模宏大、结构复杂、设施齐备，在同时期城门筑造中技术最为先进的城门遗迹，被誉为“华夏第一门”，同时也是在中国城建史中瓮城建造的开端，对了解我国早期城建史有着重要意义^[5-6]。

外城东门主要由南北两墩台组成，如图2所示。南、北墩台只剩下内侧夯土和外城墙。北墩台内侧夯土平面尺寸约为16 m×14 m，最高处约为6.7 m，外侧城墙厚2.2~4.1 m。南墩台内侧夯土平面尺寸约17 m×6.8 m、外侧城墙厚约4.2 m，最高处约5.6 m。由于墩台内部为夯土，外侧城墙不仅承担自身重力，更主要承担夯土的侧向土压力，起到了重力式挡土墙作用。

外城东门典型病害类型为墩台墙体向外倾斜，如图3所示。经现场查勘，目前南北墩台均出现了向外倾斜的情况，墩台墙体倾斜数据如表1所示，南北墩台墙体偏角达到5°~9°，墙体顶部位移15 cm以上，倾倒危险较大，故需尽快对其进行保护加固。

2 计算方法

石峁城墙是由毛石和黄土泥浆组成的砌体结构。由于砌体结构中的砌块和灰缝的力学性能差异巨大，砌体力学响应呈现出高度非线性，砌体结构计算分析一直是一项具有挑战性的任务。

2.1 砌体结构建模方法

目前，根据砌体结构建模的精细程度可分为连续模型和分离模型两类^[15]。砌体结构计算模型如图4所示。连续模型将灰缝和砌块看作均质的等效单元，采用均值等效的砌体本构模型。分离模型将灰缝和砌块分开建模，包括简化模型和复杂模型。简化模型将砌块作为连续单元，而将灰缝看作接触单元。复杂模型是将砌块和灰缝均作为连续单元，砌块与灰缝间设置接触单元^[16]。

连续模型适用于砌体整体结构性能分析，建模和计算成本较低，但无法得到细致的结果。复杂模型适用于小型砌块试样的开裂破坏行为分析，建模过程繁琐，计算成本高。简化模型在建模和计算成本方面介于前两者之间，可用较低的计算成本模拟砌体结构的整体受力性能和灰缝的破坏特征^[17]。

鉴于石峁城墙组成砌块的随机性和复杂性，为了兼顾计算的准确性和高效性，采用离散单元法建立简化模型进行研究分析。

2.2 离散单元法计算原理

离散单元法（DEM）是1971年由Cundall提出的，主要是为了解决非连续介质大变形问题^[18]。采用离散单元法模拟砌体可以将砌体单元作为离散块体，将灰缝看成接触，所以可以直接考虑砌体的截面形态^[19]。在数值分析过程中，可以捕捉计算模型的大变形、灰缝错位和块体分离，因此可以用于分析砌体结构中裂缝扩展和倒塌机制^[20]。

在接触面处，块体通过多组点接触相互连接。点接触的主要优点是其通用性和简单性，能够处理块之间的边和边及点和边等各种类型的几何相互作用。由此将法向力或剪切力从一个块传递到另一个块。并且可以考虑块体大的运动，包括在循环加载下分离和重新闭合的情况^[21]。

接触面的力学模型如图5所示，接触面的法向应力和剪切应力的计算公式如下：

式中：σ、τ为接触面法向和剪切应力；k_n为法向刚度；k_s为剪切刚度；u_n为法向应变；u_s为剪切应变。

当接触面法向应力大于法向粘结强度f_t时，块体发生分离，同时法向应力为零。切向接触采用Mohr-Coulomb模型，当剪切应力大于剪切强度时，块体发生相互滑动。计算公式如下：

式中：φ为接触面的内摩擦角；c为接触面粘聚力。

2.3 稳定性计算方法

稳定性计算主要采用强度折减法^[22-23]。强度折减法是将岩土体的抗剪强度指标（c和tanφ）同时除以折减系数F_s，从而得到一组新的强度指标，然后进行数值分析，反复折减计算至岩土体达到临界破坏状态，此时的折减系数即为边坡稳定性系数^[24]。和传统方法相比较，强度折减法可以分析复杂的边坡形式，不需要假定滑移面的形式，并可以计算岩土体和结构的相互作用^[25]。适合用于石峁城墙的稳定性计算。

根据GB 50330—2013《建筑边坡工程技术规范》规定：边坡稳定性状态分为稳定、基本稳定、欠稳定和不稳定四种状态。考虑石峁城墙的永久性及其文物价值，稳定性状态划分如表2所示。

3 方法验证

通过室内模型试验和数值模拟试验相互对比分析，验证离散元建模方法的有效性。

3.1 室内模型试验

东门墙体宽度约为2~4 m，高度约为3~6 m。依据取样尺寸及室内场地试验条件，将试验模型比例确定为1∶4.5，试验模型如图6所示，模型高900 mm，顶长800 mm，底长1 000 mm，顶厚300 mm，底厚550 mm。石墙由人工分层砌筑。砌筑的石块长70~110 mm，宽15~25 mm，高10~22 mm。砌筑灰缝8~12 mm。砌筑墙体的石块与灰缝所用黄土均取自石峁遗址当地，经试验测定，灰缝材料参数为：内摩擦角29°，粘聚力48 kPa。

整个模型位于木架平台上，平台一侧铰接固定，另一侧通过千斤顶加载，使得平台一侧升高。石墙随平台升高发生倾斜，直至墙体倒塌。试验采用摄像机进行实时拍摄，记录倒塌过程，同时在墙体底部及侧面均放置尺子记录墙体倒塌时平台的倾斜角度。

试验过程如图7所示，木板与地面倾斜角度为40°时，墙体发生失稳破坏。破坏现象为墙体整体向前倾倒，上半部灰缝出现裂纹，触地后全部破坏；下半部墙体小部分留存在木板上，其余大部分倾倒。

3.2 数值模拟

根据室内试验模型尺寸，使用UDEC数值模拟软件建立的模型如图8所示，模型共186个块体，单元9 819个，接触3 146个。接触采用摩尔库伦模型，材料参数与室内试验参数相同。模拟时，逐渐增大支撑底板的水平夹角α，计算研究不同倾斜角条件下墙体的失稳破坏情况。

模拟结果表明，当底板水平夹角小于42°时，墙体均处于静止状态；当等于42°时，墙体开始转动，最终发生倒塌破坏。倒塌过程如图9所示，墙体三分之一处出现一道裂缝将墙体分为上、下两部分，下部墙体完整，上部墙体内有大量的裂缝；最终墙体完全倒塌，下部墙体依旧完整，而上部墙体的不同块体相互分离，墙体完全破碎。

3.3 结果分析

通过比较室内模型试验和数值模拟试验结果可知，数值模拟的墙体倒塌角度和试验角度相差2°，数值模拟墙体的倒塌过程、裂缝出现及其分布特征均和室内试验相近。说明采用离散单元法建立的简化模型可以用于石峁遗址东门墙体研究。

4 东门模型计算

4.1 计算模型

采用UDEC数值模拟软件进行计算分析。以图2中断面A-A为研究对象，将其简化为二维平面应变问题，城墙的岩石块体组成是通过对现场典型断面拍照，然后采用AutoCAD描绘而成。城墙概化模型高3.5 m，宽2.5 m，岩石块体共401块，岩石块体间的灰缝接触共7 960个，城墙断面如图10所示。从图中可以看出，城墙两侧是两片叶墙，叶墙的岩石块体扁平，块体高约10~20 cm，宽度约为30~50 cm，灰缝大致呈水平状。两片叶墙中间的岩石块体形状及大小较为随机，灰缝呈不规则分布。

城墙计算模型和网格划分如图11所示，模型墙体高3.5 m，宽度为2.5 m，墙后土体宽度4 m，墙体底部的地基厚度为1.5 m。墙体网格尺寸为0.08 m，土体网格尺寸为0.15 m，单元网格共261 841个。考虑到后期将在墙后土体顶部设置步道供游人行走参观，将人群荷载和步道自重荷载进行荷载组合，同时考虑文物特殊性，最终取墙后土体竖向外荷载为15 kPa。模型底部边界约束竖向位移，模型左侧和右侧约束水平位移，其余为自由边界。

4.2 材料参数

材料参数主要通过室内试验确定。岩石和墙后土体均采用摩尔库伦本构模型，通过室内岩石直剪和单轴抗压试验确定岩石力学参数，通过三轴试验确定墙后土体力学参数，试验结果如表3所示。

土体与砌体接触采用摩尔库伦模型，接触面内摩擦角为20°，剪切刚度和法向刚度为1 GPa^[26]。灰缝接触采用摩尔库伦模型，通过块体剪切试验确定灰缝的粘聚力、内摩擦角和剪切刚度，通过压缩试验确定灰缝的法向刚度，通过拉伸试验确定灰缝的法向强度。试验如图12所示，岩石块体尺寸为300 mm×200 mm×100 mm，灰缝厚度为20 mm。灰缝采用石峁黄土和水搅拌制作而成，根据当地石砌经验，需将黄土加水搅拌至可塑状态。经测定石峁黄土液限为25.9%，塑限为16.8%，结合对当地匠人的走访调查，本试验水和黄土质量比取0.19∶1。试件制作方法为：先放置好底层条石，在其上铺一层泥浆，手工整平后放置上层条石。砌筑完成后，在自然条件下养护28 d，待灰缝自然风干后进行试验。考虑到石峁城墙的实际高度，竖向压力取50 kPa、100 kPa和150 kPa。由于灰缝泥浆抗拉强度很低，本次拉伸试验没有测得有效数据，因此灰缝法向强度取零。试验结果如表4所示。

4.3 计算工况

石峁城墙历经约4 000年，石块和灰缝材料均发生了不同程度风化，特别是灰缝，几乎所有灰缝表面已经粉化脱落，形成了深度为2~6 cm的灰缝沟，如图13所示。说明墙体材料力学性能因为风化已经发生退化，其强度和刚度均产生了一定程度的降低。

由于室内试验的岩石取自现场，其测定的参数代表了现在即风化后的岩石力学性能，而灰缝为现制的黄土泥浆，其测定值代表了城墙刚建造时的灰缝力学性能。显然，室内试验确定的灰缝参数无法直接用于现在城墙的稳定性分析。由于文物的特殊性，现场无法进行原位试验。因此，为了研究灰缝风化对城墙稳定性的影响并确定灰缝的风化程度，本研究提出通过风化折减系数f来表征不同的风化程度。假定风化后灰缝刚度和强度通过乘以f进行折减，主要研究如下三个方面内容：

1）研究灰缝刚度变化对城墙稳定性影响。保持灰缝强度不变，对灰缝刚度进行折减，然后计算城墙稳定性。刚度折减计算公式如下：

式中：k_sf、k_nf为风化折减后的灰缝刚度。

2）研究灰缝强度变化对城墙稳定性影响。保持灰缝刚度不变，对灰缝强度进行折减，然后计算城墙稳定性。强度折减计算公式如下：

式中：c_f、φ_f为风化折减后的粘聚力和内摩擦角。

3）根据模拟和实际城墙变形破坏特征，通过反演分析法得出现在灰缝力学参数，以便为城墙加固设计提供依据。

5 稳定性分析

5.1 灰缝刚度对稳定性影响

为了研究灰缝刚度变化对城墙的影响，保持灰缝强度不变仅对灰缝刚度进行折减。不同灰缝刚度折减后城墙稳定性如表5所示。从表中可以看出，随着刚度减小城墙稳定性系数也同时减小，但由于城墙的稳定性系数变化很小，稳定性状态并没有发生变化。可以认为，城墙灰缝刚度减小对城墙稳定性没有影响。

不同灰缝刚度城墙最大剪应力分布如图14所示。从图中可以看出，城墙应力呈不均匀分布，在临空侧墙脚应力集中，这主要是由于水平土压力导致城墙偏心引起的。城墙最大集中应力随着灰缝刚度降低逐渐增大，但其变化幅度较小，分别为0.95 MPa、0.99 MPa、1.06 MPa和1.14 MPa。城墙所有岩石块体均未发生破坏。

不同灰缝刚度折减后模型水平位移如图15所示。从图中可以看出，城墙水平位移均为上面大下面小的倾覆位移模式。随着灰缝刚度减少城墙的水平位移逐渐增大，城墙顶部水平位移分别为6.51 mm、6.9 mm、7.3 mm和8.4 mm。

通过上述分析可以得出，灰缝刚度降低增大了城墙的侧向位移，但对墙体的稳定性没有影响。

5.2 灰缝强度对稳定性影响

为了研究灰缝强度变化对城墙的影响，保持灰缝刚度不变仅对灰缝强度进行折减。不同灰缝强度稳定性如表6所示。从表中可以看出，随着折减系数减小，稳定性系数也不断减少。城墙灰缝强度降低，导致墙体从稳定状态逐渐过渡为不稳定状态。灰缝强度是影响城墙稳定的主要原因。

不同灰缝强度城墙最大剪应力分布如图16所示。从图中可以看出，城墙应力分布都有应力集中现象，没有风化的城墙应力集中位于临空侧的墙脚，随着灰缝强度降低应力集中位置逐渐向城墙中间转移。这主要和城墙滑移变形有关，最大剪切应力主要位于发生滑移变形附近。最大剪应力随着灰缝强度降低而增大，分别为0.95 MPa、2.71 MPa、3.08 MPa和5.36 MPa。在不同情况下，城墙所有岩石块体均未发生破坏。

不同灰缝强度模型水平方向位移如图17所示。从图中可以看出，随着灰缝强度降低，墙体向临空面的位移增大，当灰缝强度降低至0.50时，城墙本体出现裂缝。

风化折减系数为0.5和0.25时模型速度向量如图18所示。从图中可以看出，由于灰缝强度降低，城墙无法阻止墙后土体向墙体方向运动，在土体内形成倾斜的滑移面，下滑的楔形土体给墙体提供了侧向土压力。值得注意的是，经典的主动土压力滑移面通过墙脚，而石峁城墙的墙后土体内的滑移面高于墙脚，且高度随着风化折减系数减小而增高。这是因为城墙含有大量灰缝，不能将城墙看为刚性挡墙，且灰缝强度越低，城墙整体刚度越小。

风化导致墙体的变形及破坏类型如图19所示。从图中可以看出，城墙变形破坏呈倾覆-滑移模式。如图19a所示，当灰缝强度为原强度的0．5倍时，虽然城墙处于基本稳定状态，但城墙底部出现裂缝，城墙上部出现整体向外移动。如图19b所示，当灰缝强度为原强度的0．25倍时，城墙底部出现倾斜裂缝，导致上部墙体向外整体移动，城墙上部外叶墙内侧出现竖向裂缝，导致城墙的外叶墙和城墙主体发生分离，如果进一步发展，城墙将会倒塌。石峁东门多处出现城墙外倾现象，主要和灰缝强度降低有关。

5.3 东门城墙灰缝风化程度

目前，东门城墙存在向外变形但并未形成明显的竖向裂缝，也没发生整体失稳破坏。数值模拟结果得出城墙稳定性系数为1.0~1.05，根据表2可以判断城墙处于欠稳定状态。通过反演分析法，使用UDEC模拟不同的风化折减系数，得到不同折减系数下城墙的变形及破坏状态，然后与城墙现状进行对比，反复试算至状态相近，得出当前城墙的风化折减系数为0.45~0.48，灰缝的强度参数：内摩擦角为14°，粘聚力为21 kPa，并以此作为后续加固设计依据。

6 加固方法

通过对城墙稳定性分析可知，在水平土压力作用下，随着城墙风化，灰缝力学性能退化，城墙稳定性逐渐降低并向临空一侧运动。因此，必须采取措施对已经发生侧向变形的城墙进行加固，以防止其发生失稳倒塌。可采取两种方法进行加固：提高灰缝力学性能的灰缝注浆法和阻止墙体进一步变形的外侧支撑法。

6.1 灰缝注浆法

灰缝注浆法是目前砌体古建筑加固较为常用的方法^[27-28]。针对石峁城墙，注浆材料应采用矿物聚合物灌浆材料。加固包括四个步骤：1）钻孔。沿城墙侧面水平灰缝和顶部灰缝钻圆形孔，孔径20 mm，间距0.5 m；2）安装注浆花管。将花管插入孔内，长度与钻孔深度相同；3）注浆。将浆液通过注射器注入灰缝，为保证浆液饱满需加压反复注射。施工时应从下向上逐层注射，侧面完成后再进行顶部注射；4）注浆口封堵。浆液完成后，用黄土泥浆将孔口封闭。当注浆完成后，灰缝强度提高（粘聚力c=300 kPa、内摩擦角φ=30°），稳定性系数为5.7，模型整体水平位移如图20所示，最大水平位移为1.8 mm。

6.2 外侧支撑法

外侧支撑主要是由立柱、斜撑、水平联系板组成的空间结构，可以提供较大的水平支撑力阻止城墙向临空一侧移动。为了和环境协调，支撑材料可选用木材，施工步骤包括：1）构件制作。构件全部在加工厂制作，统一运至现场；2）现场安装。按照水平联系板、立柱和支撑顺序进行分段安装。设置支撑后城墙稳定性系数为6.4，模型整体水平位移如图21所示，最大水平位移为1.6 mm。

6.3 加固方法比选

从加固机理和数值模拟效果来看，灰缝注浆和侧向支撑两种方法均能减小侧向变形并提高城墙稳定性，但两者各有优缺点。

灰缝注浆法优点是对城墙外貌无干扰，缺点是施工繁琐，工期较长。特别是由于城墙内部的灰缝方向杂乱，注浆孔直接到达城墙内部非常困难，加固效果很难保证。侧向支撑法优点是体系受力明确，施工简单，效果可控，且后期拆卸方便，是可逆的加固方法。缺点是对城墙的外貌有干扰。

同时由于陕北神木地区当地暴雨天气频发，这对灰缝注浆法的效果影响较大。因此，综合考虑后，决定采用侧向支撑法。城墙灰缝风化后灰缝强度下降，会导致墙体从稳定状态逐渐过渡为不稳定状态，从而增大对木支撑压力，为保证支护结构安全，取稳定性系数为2.0进行计算，通过计算得出不同灰缝强度下墙体达到稳定所需支撑力大小，如表7所示。同时结合实际情况，最终选择2 000 mm×200 mm×125 mm方木和直径180 mm圆木建造木支撑对墙体进行支护，如图22所示。支撑中木材的抗压强度为14 MPa，抗拉强度为10 MPa。通过分析计算，木支撑压应力为2.5 MPa，支撑间木板拉应力为6.8 MPa，均满足当前支护需求。

在支护加固后，采用人工方式对东门已支护加固后墙体进行监测，位移监测仪器采用科力达KTS-442R10C型全站仪。外城东门支护加固后墙体位移监测布设监测点共计6个测点，数据统计周期为十二个月，各监测点水平位移变化范围0~2.1 mm，与数值模拟结果相近。从图23中可以看出，城墙在支护加固后，墙体水平位移增长逐渐变缓，城墙处于稳定状态。

综合数值模拟结果和监测数据可以得出结论：城墙在支护加固后侧向变形趋于稳定，取得良好的效果。

7 结论

本研究首先通过室内模型试验和数值模拟试验验证离散元建模方法的有效性。然后利用离散单元法对石峁遗址外城东门挡墙类城墙进行了变形及稳定性分析，研究结论如下：

1）由于岩石块体强度较高，石峁城墙的变形及稳定性主要取决于灰缝性能。风化导致灰缝刚度和强度降低，灰缝刚度降低增加了城墙的侧向变形但不影响其稳定性，而灰缝强度降低是城墙变形失稳的主要原因。

2）在土压力作用下城墙变形及破坏呈倾覆滑移模式。灰缝强度降低为原来一半时，城墙虽然处于基本稳定状态，但城墙底部出现了倾斜裂缝，随着灰缝强度进一步降低，城墙外叶墙内侧出现竖向裂缝，最终将会导致城墙倒塌。

3）根据城墙变形现状和数值模拟结果，得出城墙处于欠稳定状态，利用反演分析法确定城墙的风化折减系数为0.45~0.48，灰缝的粘聚力为21 kPa，内摩擦角为14°，为加固设计提供了依据。

4）石峁城墙加固可以采用灰缝灌浆和侧向支撑两种方法。通过比选两种方法优缺点，同时考虑工程实际情况，最终采用了侧向支撑进行加固，取得了良好的效果。

参考文献