基于扩展有限元法的坝-基系统地震响应及破坏分析方法

基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法
技术领域
1.本发明涉及结构动力学数值模拟技术领域，具体涉及一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法。


背景技术：

2.诸如大坝等大型工程建筑物的抗震安全性能是设计、施工、运维等各方最为关注的要素之一，建筑物在地震作用下的响应分析及破坏分析是不可或缺的重要环节。建筑物的地震动态响应及破坏分析涉及多个学科分支，如材料力学、结构动力学、断裂损伤力学等，而数值方法也是当前不可或缺的分析工具，常用的数值方法有有限单元法、离散元法、扩展有限元法等。
3.目前诸多学者已经开展了关于坝
‑
基系统的地震动态响应分析与破坏分析的研究，包括室内试验和相关的数值模型研究，得出了许多有价值的结论。但是这些研究通常是将响应分析和破坏分析分别考虑，缺乏从动载荷引起结构响应至破坏的全过程性分析，并且几乎很少涉及考虑坝
‑
基系统材料的流固两相动态响应分析，其中主要的难点在于当前两相介质动力理论研究的复杂性，以及结构动态断裂问题分析及动态裂纹追踪模拟的局限性。


技术实现要素：

4.本发明的目的是为了充分发挥扩展有限元法处理断裂问题时不连续界面独立于网格的特性，而提出的一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法。
5.基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法，它包括以下步骤：
6.步骤一：建立坝
‑
基系统实体模型，确定各区域的材料信息并赋予相应的材料属性，采用标准有限元法的网格生成算法划分单元网格；
7.步骤二：引入位移逼近模式和压力逼近模式，基于初始裂纹信息更新单元类型信息和节点附加自由度信息；
8.步骤三：定义满足位移逼近模式和压力逼近模式的试函数，构造控制方程弱形式，离散化处理得到基于扩展有限元的控制方程组；
9.步骤四：构造饱和多孔介质粘弹性人工边界，完成人工边界地震动载荷输入；
10.步骤五：进行时域积分计算，得到每个载荷步下坝
‑
基系统状态；
11.步骤六：对时域内每个载荷步进行坝
‑
基系统地震响应及破坏分析，追踪坝体裂纹演化过程。
12.在步骤二中，具体包括以下步骤：
13.1：引入间断函数和裂尖附加函数构造结点位移逼近模式，引入水压附加加强函数构造结点压力逼近模式；
14.2：基于初始裂纹信息，确定所有单元类型和节点类型，其中单元类型包含常规单元、贯穿单元、裂尖单元三种，结点类型包含常规结点、间断加强结点、裂尖奇异加强结点三
种。
15.在步骤1中，构造的结点位移逼近模式和结点压力逼近模式如下：
[0016][0017][0018]
其中，x代表域内的任意点，c代表域内所有结点集合，c
cut
为所有被裂纹贯穿单元的结点集合，c
tip
为裂尖影响范围内结点集合；u
i
、a
i
、b
i
分别为结点i的位移常规自由度、间断加强自由度、裂尖奇异加强自由度；p
i
、c
i
分别为结点i的常规孔隙压力和加强的孔隙压力部分；n
i
(x)为点x的形函数；h(x)为广义的heaviside函数；φ(x)是与局部坐标相关的矩阵形式的加强函数；β
t
是坐标转换矩阵；ψ(x)为结点水压附加加强函数。
[0019]
上述函数h(x)、φ(x)、ψ(x)的形式表达具体如下：
[0020][0021][0022][0023]
其中，(r,θ)是点x的局部极坐标；κ为材料相关的常数；l(x)为点x的水平集函数值，x
γ
为裂纹γ
c
上的点，表示计算域内任意点x到γ
c
的最短距离，即为满足的那个点，n是γ
c
的外法向向量，sign()是符号函数，当ξ大于0时sign(ξ)值为1，否则sign(ξ)为
‑
1。
[0024]
在步骤三中，具体包括如下步骤：
[0025]
1：将坝体和地基视为多孔介质，在拉格朗日框架下，考虑饱和情况，忽略流相加速度，推导建立坝
‑
基系统控制方程，具体形式如下：
[0026][0027][0028]
其中，为梯度算子；σ为柯西应力张量；f
b
为平均重度，f
b
＝ρg，g为重力加速度；χ为阻尼系数；u为位移向量；ρ为平均密度，ρ＝(1
‑
n)ρ
s
+nρ
f
，ρ
s
和ρ
f
分别为固相和流相密度，n是孔隙率；k
f
为渗透系数；f
fb
为流相的重度，f
fb
＝ρ
f
g；p是平均压力，α是p的影响系数，α＝1
‑
k
t
/k
s
，k
t
、k
s
分别为多孔介质体积弹性模量和固相体积弹性模量；q为材料相关常数，
[0029]
2：定义满足位移逼近模式和压力逼近模式的试函数u
h
和p
h
，将相应的试函数变分
形式δu
h
和δp
h
分别乘以控制方程，考虑边界约束条件，在域内积分展开得到相应的弱形式；
[0030]
3：将位移逼近模式和压力逼近模式带入控制方程的弱形式，进行离散化处理，得到一组基于扩展有限元法的控制方程组，具体形式如下：
[0031][0032][0033]
其中，u为结点的位移自由度列阵；p为结点孔隙水压自由度列阵；m为质量矩阵，m
rs
＝∫
ω
ρ(n
r
)
t
n
s
dω；c为阻尼矩阵，c
rs
＝∫
ω
χ(n
r
)
t
n
s
dω；k为刚度矩阵，k
rs
＝∫
ω
(b
r
)
t
db
s
dω；o为耦合项矩阵，o
rp
＝∫
ω
(b
r
)
t
αmn
p
dω；w为与渗流作用相关的惯性项矩阵，w
ps
＝∫
ω
(b
p
)
t
k
f
ρ
f
n
s
dω；h
pp
＝∫
ω
(b
p
)
t
k
f
b
p
dω为渗流相关矩阵；s为压缩性相关矩阵，s
pp
＝∫
ω
(n
p
)
t
n
p
/qdω；f
ext
为外载荷列阵，f
int
为内载荷列阵，上下标r、s、表征位移自由度类型，即标准的、间断加强的、裂尖奇异加强的三种类型；上下标p、ζ表征压强自由度类型，即标准的和加强的两种类型；n
r
、n
s
、为位移形函数，n
p
、是压力场的形函数；b是相应形函数的偏导矩阵；符号表示函数的跳跃，即函数在裂纹面上侧的值减去下侧的值；是裂纹面的单位法向；为裂纹中流体的平均流速。
[0034]
在步骤四中，具体包括如下步骤：
[0035]
1：在坝底各方向上延伸约两倍的坝高距离作为参与计算的近域地基范围，在所述的近域地基范围的边界上建立人工边界；
[0036]
2：考虑波动以柱面波形式传播，采用柱坐标系(r,θ,z)，得到两相体多孔介质粘弹性人工边界条件，具体边界条件如下：
[0037]
(1)构建人工边界的法向边界条件，考虑平面内p波的波动问题，假设p波从扰动源位置向四周传播，忽略流相加速度，得到平面内法向应力(σ
r
)边界条件和流量(q
r
)边界条件，具体形式为：其中，u
r
对应柱面坐标轴下在r方向的位移，c
p
为p波波速。
[0038]
(2)构建人工边界的平面内切向边界条件，考虑平面内剪切波的波动情况，即扰动源在θ方向振动，剪切方向水压为0，总应力与有效应力相同，得到人工边界在平面内切向应力(τ
rθ
)边界条件，具体形式为：其中，u
θ
对应柱面坐标轴下在θ方向的位移，其中g为剪切模量，c
s
为剪切波波速。
[0039]
(3)构建人工边界的出平面切向边界条件，考虑出平面剪切波的波动情况，即扰动源在z方向振动，得到人工边界在出平面切向应力(τ
zr
)边界条件，具体形式为：其中，u
z
对应柱面坐标轴下在z方向的位移。
[0040]
3.依据所述的人工边界条件，在人工边界结点上施加的等效结点力f
b
，所述等效
结点力的具体形式如下：
[0041][0042]
其中，u
b
为人工边界上结点的位移场；n为边界外法线方向向量；σ
b
为人工边界结点处的自由场应力张量；k
b
和c
b
分别为人工边界的弹簧刚度和阻尼系数，对应人工边界法向应力条件时取k
b
＝2μ/r
b
、c
b
＝ρc
p
，其中r
b
为结点到扰动源的距离；对应人工边界平面内剪切向应力条件时取k
b
＝3g/(2r
b
)、c
b
＝ρc
s
；对应人工边界出平面剪切向应力条件时取k
b
＝g/2r
b
、c
b
＝ρc
s
；a
b
为结点的影响面积；
[0043]
4.输入地震波，得到等效结点力f
b
在各个方向上的具体形式，完成人工边界地震动输入。
[0044]
在步骤五中，采用广义newmark法进行时域积分计算。
[0045]
所述广义newmark法的具体格式如下：
[0046][0047][0048][0049]
其中，参数β、γ、η是newmark参数，一般在0到1之间取值；t时刻的量u
t
、及p
t
、为已知，t+δt时刻对应的量u
t+δt
、及p
t+δt
、待求解。
[0050]
在步骤五中，将广义newmark法格式代入步骤三中的控制方程组，得到t+δt时刻的非线性方程组，采用迭代算法进行求解，即可获得当前载荷步下坝
‑
基系统全部结点的位移和水压。
[0051]
在步骤五中，具体包括如下步骤：
[0052]
1.根据步骤五中所获得的当前载荷步下全部结点的位移和水压，得到坝
‑
基系统的应力分布情况；
[0053]
2.基于断裂准则，判断当前载荷步下坝
‑
基系统是否发生裂纹萌生和裂纹扩展，将更新后的坝
‑
基系统状态作为下一个载荷步的初始状态，步骤五中的步骤2的具体步骤如下：
[0054]
(1)遍历坝
‑
基系统模型的所有单元，分析每个单元的应力状态；
[0055]
(2)采用最大主拉应力准则作为断裂准则，判断所有单元是否有出现新的断裂破坏现象；
[0056]
(3)如果坝
‑
基系统模型的所有单元都不满足步骤(2)中所述的断裂准则，则保留所有单元的类型信息和结点加强信息，进入下一个载荷步；
[0057]
(4)如果坝
‑
基系统模型中有单元满足步骤(2)中所述的断裂准则，更新潜在的裂纹路径，直至当前载荷步下所有单元都不满足步骤(2)中所述的断裂准则，更新所有单元的类型信息和结点加强信息，进入下一个载荷步；
[0058]
(5)将当前载荷步下更新后的坝
‑
基系统状态作为下一个载荷步的初始状态，进入下一个载荷步，更新步骤三的步骤3中所述控制方程组中相关矩阵、列阵，循环所述步骤五和步骤六，直到加载完成。
[0059]
与现有技术相比，本发明具有如下技术效果：
[0060]
1)发明考虑了坝
‑
基系统流固耦合作用效应，建立了两相体多孔介质扩展有限元模型，构造了位移逼近模式和压强逼近模式，可以揭示不连续界面对整体位移响应、孔压分布以及流体流动等方面的影响；
[0061]
2)本发明给出了两相体多孔介质粘弹性人工边界的构造过程，直观地阐述了坝
‑
基系统的地震动载荷输入方法；
[0062]
3)本发明充分发挥了扩展有限元法处理断裂问题时不连续界面独立于网格的特性，实现了坝
‑
基系统的地震响应及破坏全过程分析，计算效率高，在工程建筑物抗震安全性能分析方面具有很强的实用性。
附图说明
[0063]
下面结合附图和实施例对本发明作进一步说明：
[0064]
图1是本发明实施例提供的一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法的流程图；
[0065]
图2是本发明实施例提供的一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法的原理图；
[0066]
图3是本发明实施例提供的实体模型结构示意图；
[0067]
图4是本发明实施例提供的实体模型网格划分示意图；
[0068]
图5是本发明实施例提供的水平方向地震波加速度时程曲线；
[0069]
图6是本发明实施例提供的竖直方向地震波加速度时程曲线；
[0070]
图7是本发明实施例提供的计算结果的水平方向坝顶位移响应时程；
[0071]
图8是本发明实施例提供的计算结果的竖直方向坝顶位移响应时程；
[0072]
图9是本发明实施例提供的计算结果的坝顶水平方向加速度反应时程；
[0073]
图10是本发明实施例提供的计算结果的坝顶竖直方向加速度反应时程；
[0074]
图11是本发明实施例提供的计算结果的坝踵处孔压时程曲线；
[0075]
图12是本发明实施例提供的计算结果的坝趾处孔压时程曲线；
[0076]
图13是本发明实施例提供的计算结果的坝踵处的最大主应力地震响应时程曲线；
[0077]
图14是本发明实施例提供的计算结果的最大主应力云图(t＝2.63s)；
[0078]
图15是本发明实施例提供的计算结果的坝体裂纹扩展示意图。
具体实施方式
[0079]
如图1所示，一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法，包括如下步骤：
[0080]
(1)建立坝
‑
基系统实体模型，确定各区域的材料信息并赋予相应的材料属性，采用标准有限元法的网格生成算法划分单元网格；
[0081]
(2)引入位移逼近模式和压力逼近模式，基于初始裂纹信息更新单元类型信息和
节点附加自由度信息；
[0082]
(3)定义满足位移逼近模式和压力逼近模式的试函数，构造控制方程弱形式，离散化处理得到基于扩展有限元的控制方程组；
[0083]
(4)构造饱和多孔介质粘弹性人工边界，完成人工边界地震动载荷输入；
[0084]
(5)采用广义newmark法进行时域积分计算，求解得到每个载荷步下坝
‑
基系统状态；
[0085]
(6)对时域内每个载荷步进行坝
‑
基系统地震响应及破坏分析，追踪坝体裂纹演化过程。
[0086]
进一步，上述步骤(2)中，引入间断函数和裂尖附加函数构造结点位移逼近模式，引入水压附加加强函数构造结点压力逼近模式；基于初始裂纹信息，确定所有单元类型和节点类型，其中单元类型包含常规单元、贯穿单元、裂尖单元三种，结点类型包含常规结点、间断加强结点、裂尖奇异加强结点三种。其中，构造的结点位移逼近模式和结点压力逼近模式如下：
[0087][0088][0089]
其中，x代表域内的任意点，c代表域内所有结点集合，c
cut
为所有被裂纹贯穿单元的结点集合，c
tip
为裂尖影响范围内结点集合；u
i
、a
i
、b
i
分别为结点i的位移常规自由度、间断加强自由度、裂尖奇异加强自由度；p
i
、c
i
分别为结点i的常规孔隙压力和加强的孔隙压力部分；n
i
(x)为点x的形函数；h(x)为广义的heaviside函数；φ(x)是与局部坐标相关的矩阵形式的加强函数；β
t
是坐标转换矩阵；ψ(x)为结点水压附加加强函数。其中，函数h(x)、φ(x)、ψ(x)的形式表达具体如下：
[0090][0091][0092][0093]
其中，(r,θ)是点x的局部极坐标；κ为材料相关的常数；l(x)为点x的水平集函数值，x
γ
为裂纹γ
c
上的点，表示计算域内任意点x到γ
c
的最短距离，即为满足的那个点，n是γ
c
的外法向向量，sign()是符号函数，当ξ大于0时sign(ξ)值为1，否则sign(ξ)为
‑
1。
[0094]
进一步，上述步骤(3)的具体步骤如下：
[0095]
1)基于biot理论，将坝体和地基视为多孔介质，在拉格朗日框架下，考虑饱和情
况，忽略流相加速度，由动量守恒原理和质量守恒原理，推到建立坝
‑
基系统控制方程，具体形式如下：
[0096][0097][0098]
其中，为梯度算子；σ为柯西应力张量；f
b
为平均重度，f
b
＝ρg，g为重力加速度；χ为阻尼系数；u为位移向量；ρ为平均密度，ρ＝(1
‑
n)ρ
s
+nρ
f
，ρ
s
和ρ
f
分别为固相和流相密度，n是孔隙率；k
f
为渗透系数；f
fb
为流相的重度，f
fb
＝ρ
f
g；p是平均压力，α是p的影响系数，α＝1
‑
k
t
/k
s
，k
t
、k
s
分别为多孔介质体积弹性模量和固相体积弹性模量；q为材料相关常数，
[0099]
2)定义满足位移逼近模式和压力逼近模式的试函数u
h
和p
h
，将相应的试函数变分形式δu
h
和δp
h
分别乘以控制方程，考虑边界约束条件，在域内积分展开得到相应的弱形式；
[0100]
3)将位移逼近模式和压力逼近模式带入控制方程的弱形式，进行离散化处理，得到一组基于扩展有限元法的控制方程组，具体形式如下：
[0101][0102][0103]
其中，u为结点的位移自由度列阵；p为结点孔隙水压自由度列阵；m为质量矩阵，m
rs
＝∫
ω
ρ(n
r
)
t
n
s
dω；c为阻尼矩阵，c
rs
＝∫
ω
χ(n
r
)
t
n
s
dω；k为刚度矩阵，k
rs
＝∫
ω
(b
r
)
t
db
s
dω；o为耦合项矩阵，o
rp
＝∫
ω
(b
r
)
t
αmn
p
dω；w为与渗流作用相关的惯性项矩阵，w
ps
＝∫
ω
(b
p
)
t
k
f
ρ
f
n
s
dω；h
pp
＝∫
ω
(b
p
)
t
k
f
b
p
dω为渗流相关矩阵；s为压缩性相关矩阵，s
pp
＝∫
ω
(n
p
)
t
n
p
/qdω；f
ext
为外载荷列阵，f
int
为内载荷列阵，上下标r、s、表征位移自由度类型，即标准的、间断加强的、裂尖奇异加强的三种类型；上下标p、ζ表征压强自由度类型，即标准的和加强的两种类型；n
r
、n
s
、为位移形函数，n
p
、是压力场的形函数；b是相应形函数的偏导矩阵；符号表示函数的跳跃，即函数在裂纹面上侧的值减去下侧的值；是裂纹面的单位法向；为裂纹中流体的平均流速。
[0104]
进一步，上述步骤(4)的具体步骤如下：
[0105]
1)在坝底各方向上延伸约两倍的坝高距离作为参与计算的近域地基范围，在所述的近域地基范围的边界上建立人工边界；
[0106]
考虑波动以柱面波形式传播，得到两相体多孔介质粘弹性人工边界条件，具体边界条件如下：
①
考虑平面内p波的波动问题，假设p波从扰动源位置向四周传播，忽略流相加速度，得到平面内法向应力(σ
r
)边界条件和流量(q
r
)边界条件，具体形式为：其中，u
r
对应柱面坐标轴下在r方向的位移，
c
p
为p波波速；
②
考虑平面内剪切波的波动情况，即扰动源在θ方向振动，剪切方向水压为0，总应力与有效应力相同，得到人工边界在平面内切向应力(τ
rθ
)边界条件，具体形式为：其中，u
θ
对应柱面坐标轴下在θ方向的位移，c
s
为剪切波波速，g为剪切模量；
③
考虑出平面剪切波的波动情况，即扰动源在z方向振动，得到人工边界在出平面切向应力(τ
zr
)边界条件，具体形式为：其中，u
z
对应柱面坐标轴下在z方向的位移。
[0107]
2)依据所述的人工边界条件，在人工边界结点上施加的等效结点力f
b
，所述等效结点力的具体形式如下：
[0108][0109]
其中，u
b
为人工边界上结点的位移场；n为边界外法线方向向量；σ
b
为人工边界结点处的自由场应力张量；k
b
和c
b
分别为人工边界的弹簧刚度和阻尼系数，对应人工边界法向应力条件时取k
b
＝2μ/r
b
、c
b
＝ρc
p
，其中r
b
为结点到扰动源的距离；对应人工边界平面内剪切向应力条件时取k
b
＝3g/(2r
b
)、c
b
＝ρc
s
，对应人工边界出平面剪切向应力条件时取k
b
＝g/2r
b
、c
b
＝ρc
s
；a
b
为结点的影响面积。
[0110]
3)输入地震波，得到等效结点力f
b
在各个方向上的具体形式，完成人工边界地震动输入。
[0111]
进一步，在步骤(5)中，引入广义newmark法格式，代入步骤c3中所述的控制方程组，得到t+δt时刻的非线性方程组，采用迭代算法进行求解，获得当前载荷步下坝
‑
基系统全部结点的位移和水压。其中，采用广义newmark法的具体格式如下：
[0112][0113][0114][0115]
其中，参数β、γ、η是newmark参数，一般在0到1之间取值；t时刻量u
t
、、及p
t
、已知，t+δt时刻对应量u
t+δt
、及p
t+δt
、待求解。
[0116]
进一步，根据步骤(5)中所获得的当前载荷步下全部结点的位移和水压，得到坝
‑
基系统的应力分布情况；基于断裂准则，判断当前载荷步下坝
‑
基系统是否发生裂纹萌生和裂纹扩展，将更新后的坝
‑
基系统状态作为下一个载荷步的初始状态，具体步骤如下：
[0117]
1)遍历坝
‑
基系统模型的所有单元，分析每个单元的应力状态；
[0118]
2)采用最大主拉应力准则作为断裂准则，判断所有单元是否有出现新的断裂破坏现象；
[0119]
3)如果坝
‑
基系统模型的所有单元都不满足中所述断裂准则，则保留所有单元的类型信息和结点加强信息，进入下一个载荷步；
[0120]
4)如果坝
‑
基系统模型中有单元满足所述断裂准则，则更新潜在的裂纹路径，直至当前载荷步下所有单元都不满足所述断裂准则，更新所有单元的类型信息和结点加强信息，进入下一个载荷步；
[0121]
5)将当前载荷步下更新后的坝
‑
基系统状态作为下一个载荷步的初始状态，进入下一个载荷步，更新步骤(3)中控制方程组中相关矩阵、列阵，循环上述步骤(5)和步骤(6)，直到加载完成。
[0122]
下面结合具体实施例对本发明的技术方案作进一步说明。
[0123]
本发明的一种基于扩展有限元法的坝
‑
基系统地震响应及破坏分析方法，建模、求解及分析流程如图2。
[0124]
实施例：
[0125]
本实施例以某混凝土重力坝受地震载荷作用为例，利用本发明的分析方法进行地震响应及破坏分析。大坝坝高103m，大坝的顶部宽度为14.8m，坝底宽为70m，坝段厚16m，上游坝面为垂直坡面，下游坝面由两段坡度分别为8.20:1和1.21:1的坡面构成，两段坡面在高程为66.5m处相交，在坝底各方向上延伸约两倍的坝高距离作为计算的近域地基范围，大坝计算模型如图3，网格划分如图4；坝体材料的弹性模量为31.03gpa，泊松比为0.15，孔隙率为0.19，固相密度为2463kg/m3，固相压缩模量为36.0gpa，液相密度为1000kg/m3，流相压缩模量为3gpa，抗拉强度为2.9mpa；坝基材料的弹性模量为20.00gpa，泊松比为0.2，孔隙率为0.19，固相密度为2700kg/m3，其余材料参数与坝基材料取值相同；输入的地震波如图5、图6所示，其中地震波水平向的峰值加速度为4.26m/s2，竖直向峰值加速度为3.05m/s2；时域积分采用newmark法，积分参数β取0.5，γ取0.75，η取0.75，时间步长δt取0.01s，总计算时间为10s。
[0126]
如图7～图9所示，分别为坝顶的位移地震响应时程曲线和加速度响应时程曲线，坝顶的最大水平位移为0.4380m，最大竖直向位移为0.5875m，水平向和竖直向的加速度峰值分别为18.61m/s2和9.03m/s2；坝踵和坝趾处的孔压时程曲线如图11、图12所示，坝踵处孔压的峰值较坝趾处更高，4.20s时坝踵处孔压达到峰值1.25mpa，4.37s时坝趾处孔压达到峰值0.87mpa；图13给出了坝踵的最大主应力地震响应时程曲线，最大主应力峰值为3.24mpa；同理，可给出坝
‑
基系统任意部位的位移响应、加速度响应、孔压时程、应力时程等。在此基础上对坝
‑
基系统进行地震破坏分析，追踪坝体裂纹发展过程，图14给出了2.63s时坝体的最大主应力分布云图，此时坝踵附近区域受压，坝趾和下游坝面变坡度处受拉应力，下游坝面变坡度处最先满足断裂准则，裂纹开始扩展；图15为最终得到的坝体裂纹扩展示意图，裂纹产生的区域及传播路径与振动台模型试验及相关文献的结果比较一致，验证了本方法的有效性。