采用弹塑性损伤力学的砖石建筑遗产结构性能非线性有限元分析

王浩宇，淳庆，张承文，华一唯

(东南大学 建筑学院，江苏 南京 210096)

砖石砌体建筑遗产是我国物质文化遗产的重要组成部分，记载了我国历史的变迁和文化的传承.砖石砌体建筑的材料大多由硅酸盐和碳酸盐组成，容易遭到外界因素的影响，从而产生表面破损和结构性的损伤，因此，需要运用强现代化技术，对砖石砌体建筑进行有效的保护和修缮，避免砖石砌体建筑在受到外界影响时发生损坏.

文献[1-2]研究砖石砌体结构建筑遗产表面防护修复工艺.彭程[3]研究表面做色涂料.然而，砖石砌体古建筑的保护和修缮不仅仅是对其表面、外观的修复和防护，更需要对其主体结构的承载能力进行复核与验算，以确保结构的安全性.古建筑存在较多的病害，包括外墙的裂缝、破损，内部受力构件的老化、开裂等，为了较准确地模拟古建筑当前的受力情况、计算其当前承载能力，需要采用适用于砖石砌体的损伤本构模型.张亚东等[4]基于Hoek-Brown准则的岩体损伤本构模型认为岩石材料承载能力分为弹性和损伤两个部分，但忽略了岩石结构抗压时的塑性行为.鄢俊彪等[5]从能量角度分析岩石的变形破坏特性，建立基于三剪能量屈服准则的应变软化岩石的统计损伤本构模型.黄榜彪等[6]开发一种含参数的砌体本构模型.Liu等[7]开发基于能量耗散的损伤本构模型，描述岩石在循环荷载作用下的力学行为.Roca等[8]对砌体结构的正交各向异性进行研究，并通过试验模拟了该模型对不同砌体结构的适用性.Fu等[9]建立砌体结构的正交各向异性塑性损伤本构模型.PEL等[10]研究正交各向异性损伤模型在砌体结构中的运用.

国内外关于砖石砌体材料的弹塑性损伤本构模型鲜有研究.本文在Drucker-Prager屈服准则的基础上开发了砖石砌体的弹塑性损伤本构模型，模拟砖石砌体结构损伤的发展.

1 中华门城堡概况

中华门城堡是目前中国保存最完好、规模最大的堡垒瓮城，建于公元1386年(洪武十九年)，是在南唐旧城正南门基础上扩建而成，内向有三道瓮城，从南到北分别为南闸楼、中闸楼、北闸楼，地面上三道瓮城拱券门与前正门一线贯通.中华门城堡主体历经长期的侵蚀、风化和人为的损伤，加之年久失修，各面墙体存在多处裂缝、砖块缺失与孔洞等损伤.现存的中华门东西长116 m，南北深128 m，占地14 848 m2，由主城台、瓮中城27个藏兵洞、东西马道及一条登城曲道组成.主城台是南京13座外城门中唯一一座两层结构的城台，总高19.42 m.其中：第一层为城墙，北侧东西宽约58.6 m，南北深约52.4 m，高约9.5 m，由中央主城门道、两侧对称分布的6个藏兵洞组成；第二层为城楼下部的基座，设有7个藏兵洞，北侧较底层向内收缩，长约3.5 m，留出一条东西向通道；基座顶部东西宽约65.2 m，南北长46.2 m，城台突出城墙34.6 m，高约9.9 m.中华门城堡现状图，如图1所示.

(a)南立面 (b)北立面

中华门城堡使用的砌体材料主要有两种：条石与城砖.条石主要用于下部城门墩主体，城砖主要用于上部墙体、台顶铺砖、女墙、雉碟等部位.中华门主城台的南侧全部为条石砌体，北侧二层女墙下部也均为条石砌体，北侧二层为城砖砌体.城门拱券厚约1 m，内部采用碎砖三合土回填.

2 砖石砌体结构弹塑性损伤本构模型

对于年久失修的砖石砌体结构，本构模型如果只采用弹性，或者弹塑性分析，承载力、刚度较原结构偏高，在对这些文物建筑进行保护修缮前的结构受力分析时，需要考虑材料性能的退化.在Drucker-Prager塑性屈服准则的基础上加入损伤变量，模拟结构受拉和受压损伤的出现与发展.

2.1 连续力学损伤理论

在外部荷载作用下，砌体材料内部微裂缝及微孔洞扩展、汇合、破坏的过程也是材料力学性质逐渐劣化的过程.基于塑性力学的模型无法描述内部微裂缝发展造成损伤的影响，如刚度退化等劣化现象.根据连续损伤力学理论[11]，此过程可以用损伤变量d描述.

2.2 应力-应变关系

假定损伤是各向同性的，根据应变等价原理，Cauchy应力σ为

(1)

由于材料具有明显的塑性效应，且满足小变形假定，总应变ε为

ε=εe+εp.

(2)

式(2)中：εe为弹性应变；εp为塑性应变，包含由微孔洞、微裂纹引起的各种不可逆损伤变形.

有效应力表示作用在微缺陷之间的净面积上的应力，根据应变等价原理，有效应力为

(3)

式(3)中：D为弹性本构矩阵.因此，Cauchy应力为

σ=(1-d)Dε.

(4)

2.3 塑性模型

(5)

式(5)中：Fp是塑性势函数；κ是强化函数，用于表示屈服变量与塑性变形的关系，从而确定屈服面的扩展规律.

屈服准则采用Drucker-Prager塑性屈服准则，即

(6)

式(6)中：α是材料内摩擦；I1和J2为应力不变量.

考虑κ与总塑性变形程度相关，引入等效塑性应变增量(dεp)，即

(7)

由于屈服准则与静水压力无关，即偏塑性应变增量与塑性应变增量相同，强化参数为

(8)

假设等效塑性应变率符合相关联流动法则，则塑性应变率为

(9)

式(9)中：dλ是塑性乘子.

2.4 损伤模型

为了区分材料在受拉和受压损伤演化过程的不同，根据文献[12-13]的建议，对应力采取谱分解方式，即

(10)

式(10)中：〈〉为MaCaulay括号；pi为有效应力主方向的单位列矢量.

(11)

(12)

式(12)中：f0和ft分别为受压初始屈服强度和抗拉强度；E为材料的弹性模量.

损伤演化方程[12-13]为

(13)

式(13)中：A-，B-为受压损伤演化方程的模型参数，由单轴受压损伤曲线标定；A+为受拉损伤演化方程的模型参数，为了避免计算结果的网格敏感性，A+定义[14]为

(14)

3 有限元计算分析

3.1 有限元模型

基于Fortran语言编译节2的弹塑性损伤本构模型，将文件嵌入ANSYS的usermat子程序中，对中华门主城台进行受力性能的有限元模拟分析.建模时，几何尺寸根据现场测绘获得，假设不考虑原有结构缺陷，材料参数根据现场无损检测并参考文献[15]，按偏保守的原则取值，城砖的弹性模量为2 250 MPa，容重为20 kN·m-3，极限抗压强度为0.70 MPa，极限抗拉强度为0.10 MPa，内摩擦角为30°；条石的弹性模量为3 750 MPa，容重为25 kN·m-3，极限抗压强度为1.3 MPa，极限抗拉强度为0.11 MPa，内摩擦角为30°；碎砖三合土弹性模量为20 MPa，容重为19 kN·m-3.碎砖三合土采用solid 45单元，共计100 439个，条石与城砖都采用solid 185单元，共计167 019个.选取两个工况进行静力分析：工况一考虑恒荷载的作用；工况二考虑恒荷载和活荷载的组合作用.其中，恒荷载仅考虑结构自质量，活荷载为3.5 kN·m-3的均布活荷载.有限元模型与网格的划分，如图2所示.

图2 有限元模型与网格的划分

3.2 损伤计算结果与验证

拱券损伤云图、二层城砖砌体墙及侧墙损伤云图、条石砌体墙损伤云图，分别如图3～5所示.由图3～5可知：工况一，二的受拉和受压损伤分布情况几乎相同，工况二损伤部位发展得更多.拱券受拉损伤分布在南侧拱券与城砖砌体墙的连接处、东西两侧拱券与城砖砌体墙的连接处，工况一的最大值为0.952，工况二的最大值为0.991，均出现在二层中部两拱券与砖墙交接处，这些部位受拉裂缝较多，材料承载能力下降显著.拱券受压损伤分布较少，只在东西两侧拱券底部内侧和二层两侧拱券底部内侧，工况一的最大值为0.007，工况二的最大值为0.033，说明拱券受压状况较好，因受压而出现的结构损伤发展较少.

(a)工况一受拉损伤 (b)工况二受拉损伤

(a)工况一拉损伤 (b)工况二受拉损伤

二层城砖砌体墙及侧墙受拉损伤发生在两侧和南侧墙的两侧顶部，工况一的最大值为0.998，工况二的最大值为0.999.工况一情况下，结构受压并未达到损伤阈值，损伤并未开会演化，工况二情况下，受压损伤只在北侧墙底部存在一些，最大值为0.012.

(a)工况一体受拉损伤 (b)工况二受拉损伤

条石砌体墙受拉损伤分布在东西两侧及南侧墙体上部，工况一的最大值为0.946，工况二的最大值为0.980，最大值出现在南侧墙体与东侧墙体交接处，说明该处材料受拉损伤严重.受压损伤主要分布在南侧墙体底部，工况一的最大值为0.003，工况二的最大值为0.011，说明南侧墙体底部受压损伤出现处于演化的初始阶段，但并没有成为结构受损的主要部位.

南侧条石砌体墙底部局部受压损伤图与残损现状，如图6所示.图6(a)红框中受压损伤较大处对应右侧中华门城堡南侧墙体下部出现的损伤现状，部分条石砌块与砂浆被压碎.南侧条石砌体墙受拉损伤图与上侧裂缝，如图7所示.图7(a)红框中受拉损伤较大处对应右侧中华门城堡南侧墙体上部出现的裂缝.东侧墙体受拉裂缝与受拉损伤图，如图8所示.图8(a)红框中受拉损伤较大处对应右侧中华门城堡东侧墙体上部出现的裂缝.

(a)局部受压损伤图 (b)残损现状

(a)受拉损伤图 (b)上侧裂缝

(a)受拉损伤图 (b)受拉裂缝

综上所述，采用基于Drucker-Prager弹塑性损伤本构模型的中华门城堡损伤分析计算结果与中华门现状较为吻合，这说明模型能够较好地模拟砖石砌体结构的力学行为.因此，采用文中开发的材料本构模型进行中华门城堡的有限元静力分析，更符合结构现状.受拉和受压损伤分布能够很好地反应结构整体的损伤状况，以便在修缮加固时对症下药.

3.3 静力计算结果

土体沉降图，如图9所示.墙体和拱券竖向位移图，如图10所示.由图10知：工况一土体最大位移处为东侧土体的南面顶部(26.38 mm)，墙体最大位移处为东侧墙体顶部(5.83 mm)；工况二土体最大位移处为东侧土体的南面顶部(29.27 mm)，相较工况一增大10.96%，墙体最大位移处为东侧侧墙体顶部(5.84 mm)，相较工况一增大0.17%.由于活荷载直接作用在土体表面，并且土的刚度较低所以土体位移增大较多，而马道侧墙体条石的刚度较高且位移主要由土的侧推造成，所以加上活荷载后并没有显著增大.

(a)工况一 (b)工况二

(a)工况一 (b)工况二

拱券第一、第三主应力云图，如图11所示.第一主应力云图的拉应力的分布情况主要在二层藏兵洞南面拱券与隔墙体连接处，工况一的最大值为0.101 MPa，工况二较工况一增大了0.28%，另外，北侧拱券与隔墙体连接处、两侧拱券与隔墙体连接处也存在一些拉应力分布；第三主应力云图的压应力的分布情况主要在拱券上部和拱券底部内侧区域，工况一的最大值为0.755 MPa，工况二较工况一增大了6.97%.由于二层南侧墙体与南侧条石砌体墙之间存在较多填土，二层拱券部分受力区域较靠南面.工况二拉应力的最大值超过材料极限值.因此，二层拱券、两侧拱券与隔墙体连接处易发生受拉破坏.

(a)工况一第一主应力 (b)工况二第一主应力

二层墙体第一、第三主应力云图，如图12所示.第一主应力云图的拉应力的分布情况主要在南侧墙体顶部及与东西两侧底部，工况一的最大值为0.098 MPa，工况二较工况一增大2.39%.第三主应力云图的压应力的分布情况主要集中在侧墙体及拱券落脚区域，工况一的最大值为0.771 MPa，工况二相较工况一增大5.12%.工况二拉应力的最大值超过材料极限值，东西两侧墙体底部易发生受拉破坏.

(a)工况一第一主应力 (b)工况二第一主应力

条石砌体墙第一、三主应力云图，如图13所示.第一主应力云图的拉应力的分布情况主要在东西两侧的条石砌体墙外侧，南侧墙体拉应力主要集中在两侧与顶部，工况一的最大值为0.112 MPa，工况二相较工况一增大了0.57%；第三主应力云图的压应力的分布情况主要集中在南立面墙体与东西两侧墙体的底部，工况一的最大值为1.542 MPa，工况二相较工况一增大了7.41%.工况二最大压应力与拉应力均超过材料极限值，南侧条石砌体墙底部易发生受压破坏，东侧马道侧墙体顶部易发生受拉破坏.

(a)工况一第一主应力云图 (b)工况二第一主应力云图

综上所述，工况一与工况二的情况下，结构的受压状况良好，受压应力较大处为两侧拱券和二层两侧拱券的拱顶与内侧拱脚，二层城砖砌体墙南侧墙体底部及拱券落脚，南侧条石砌体墙底部等，整体结构并未出现较多受压损伤，局部受压损伤处于损伤初始演化阶段；结构整体受拉应力较大处为拱券与隔墙连接处，二层城砖砌体墙南侧顶部及东西两侧底部，南面条石砌体墙顶部及两侧，这些地方也是受拉损伤值较大的地方，并且多处处于损伤演化末段，材料承载能力下降显著，这些受拉损伤的宏观表现形式为结构出现的裂缝.

4 结论

中华门城堡是南京明城墙的重要组成部分，极具代表性的堡垒瓮城，对该类型建筑遗产的结构进行非线性有限元分析，能更准确弄清结构的受力情况及残损病害的成因，从而为修缮保护提供科学依据.

1)基于Drucker-Prager塑性屈服准则，引入受拉和受压损伤变量，开发了考虑材料性能劣化的适用于砖石砌体结构建筑遗产的材料本构模型，并通过ANSYS的usermat子程序，实现软件的二次开发，可直接用于砖石砌体建筑遗产受力性能的有限元分析.

2)对工况一，二的中华门城堡进行静力分析，分析表明，工况二相较城堡最大拉应力提高0.28%～2.39%，最大压应力提高5.12%～7.41%.在工况一下，中华门城堡除南侧条石砌体墙底部受压应力超过材料极限以外，其余结构受压状况良好，并且受压损伤处于演化的初始阶段；结构多处受拉应力超过材料极限，并且二层城砖砌体墙两侧、南侧条石砌体墙上侧角部等部位受拉损伤处于演化末段，材料承载能力下降显著.

3)对于长期服役、年久失修的砖石砌体建筑遗产，采用提出的弹塑性损伤本构模型进行受力性能的有限元模拟分析，能够通过受拉和受压损伤变量反映结构整体的损伤状况.根据中华门当前荷载作用下受拉和受压损伤的情况可知，最大受拉损伤区域主要在二层城砖砌体墙两侧及条石砌体墙南面上侧角部，建议对这些部位进行结构安全监测或采取加固措施，防止裂缝的继续扩展.最大受压损伤区域主要在南侧砌体墙与东西两侧墙交接处底部，建议对此处进行结构安全监测并采取加固措施，以满足受压承载力的要求.