川渝拆除17713551981倾倒数值仿真模拟
1DDA数值仿真模拟
由于混凝土（常态和RCC）破坏机理的复杂性，目前国内外多采用非连续变形分析（Discontinuous Deformation Analysis，简称DDA）来仿真模拟爆破过程，诸如动态破碎、抛掷等。
本节利用DDA主要对三峡工程三期上游RCC围堰仿真模拟其爆破抛掷和倾倒过程。
1DDA原理
非连续变形分析DDA是石根华教授1984年提出的一种平行于有限元和离散元的数值分析方法，能够模拟节理裂隙岩体产生大变形和大位移的数值分析方法。DDA方法以天然存在的不连续面切割岩体形成块体单元，根据系统最小势能原理建立总体平衡方程式，将刚度、质量和荷载矩阵加到联立方程的系数矩阵中去，采用罚函数法强迫块体界面约束求解。由于它能够很好地模拟非连续介质大位移、大变形的静、动力分析等传统有限元方法难以解决的问题，现已日益广泛地应用于大坝、边坡、隧洞的稳定性分析。因为大位移和大变形包含在非连续变形分析中，块体位置、块体形状和块体接触随荷载步或时间步而变化。在块体系统中，块体的运动不允许块体之间受拉和嵌入，对每一种接触有3种模式：张开、滑动和锁定。
DDA将岩体视为非连续块体单元，块体与块体之间用虚拟的弹簧来传递相互的作用力。计算中的每一时步，用如下判据加以判别：
（1）不允许存在法向拉应力，即a，≤0。
（2）剪应力遵循摩尔一库仑（Mohr-Coulomb）准则：
1o，|≤c+ontang（3.3）式中c——节理的粘聚力；
9——内摩擦角。
块体系统运动变形时必须满足两个条件：块体界面间不嵌入和无张拉。该条件的满足通过在各接触位置加上或去掉刚硬弹簧实现。根据经验，弹簧刚度系数P=（10~1000）E，E为块体弹性模量。
1.2计算模型及基本假设
对于围堰1：10模型，采用二维计算模型，计算高
程90~140m（原型）。决定堰体翻转的主动力矩主要为重力矩，为实现围堰的倾倒，需在上游迎水面底部形成
三角形的爆破开口，典型开口断面见图3.3。开口形成后，围堰上部由于失稳而向上游一侧倾倒，其中开口的
深度和角度是决定倾倒成功的关键。
用DDA数值模拟围堰的翻转过程作以下基本假设：
（1）由于混凝土材料爆破后的整体性较好，局部的
爆破开口对整个围堰没有实质性破坏，因此可将开口后
的上部围堰视为刚性。
（2）围堰混凝土为均质材料。
（3）堰体翻转过程中与保留堰体碰撞前后满足动量守恒。
（4）计算中不考虑碰撞后堰体转动以外的其他运动形式。
1.3计算参数
本节中的计算参数包括混凝土的材料参数、节理面参数和计算时步参数。混凝土的材料参数：弹性模量E=25.5GPa，泊松比u=0.17，重度p=2400kg/m3。为简化计算，将
图3.3典型开口断面（单位：m）
水介质的作用简化为静力效应，将材料的重度设为浮重度；结构面强度，摩擦角9=35°，粘聚力c=0.5MPa，抗拉强度2MPa；破坏块体的结构面强度：摩擦角g=25°，粘聚力c
=0；刚硬弹簧系数P=1200GPa；计算时步参数为：总时步为25000步，时间步长△t
=0.00015s。
3.3.1.4爆破荷载
在利用DDA数值模拟倾倒过程时，合理地确定由炸药爆炸产生的载荷对模拟结果的准确性有着非常重要的作用，这包含确定爆破激振力的大小、作用位置和方向、作用时刻和持续时间等方面的内容。对于这样复杂的载荷要输入DDA计算程序中进行模拟计算必定是很复杂的。由于本书研究的主要内容是上部拆除堰体在爆炸载荷作用后的翻转运动过程，因此对其进行合理的简化，将荷载曲线简化为三角形波，加载到峰值压力的升压时间为100us，卸载时间为500us，总计算时间为600us。在密闭容器中，炸药爆炸产生的压力只与爆速和密度有关，而模型试验与原型爆破均在混凝土药室内，炸药重度原型为1100~1200kg/m3，模型试验为900~1000kg/m3，炸药爆速原型要求4500m/s，模型试验为3000~3500m/s。峰值压力po=pD2/8=1125×105~1531.25×106Pa。
假设爆振荷载以冲击荷载方向垂直作用在3个药室的边界面（对于模拟倾倒过程来说，做这样的假设可以满足要求）。荷载作用的次序按药室起爆的时间顺序，1号药室先起爆，2号药室滞后17ms，3号药室滞后1号药室34ms。根据药室尺寸的不同确定3个药室作用荷载的峰值分别为70t、420t、140t。
断裂孔荷载作同样的假设，升压时间100us，卸载时间500us，总计算时间为600us，其峰值压力的确定如表3.2所示。
表3.2断裂孔荷载计算表
10模型倾倒过程模拟结果分析
根据以上所建立的力学模型，采用DDA对围堰1：10模型的倾倒过程进行模拟，倾倒过程如图3.4所示。
计算结果：0.27s断裂面出现裂缝，开始翻转；1.83s拆除堰块与下部残留堰块闭合，发生碰撞；3.06s拆除堰块与地面发生碰撞；3.36s堰块完全倒地。倾倒部分与残留体最近距离为0.8m；顶部矩形块与倾倒部分脱离，距离为0.5m。
对比1：10模型试验结果：第一次模型试验倒地时间分别为第一块3.4s、第二块2.6s，倾倒部分与残留体最近距离为1m，顶部矩形块与倾倒部分脱离，距离为0.75m。
第二次模型试验倒地时间分别为第一块3.0s、第二块3.4s。
图3.4围堰1：10模型倾倒过程模拟
（a）t=375ms；（b）t=750ms；（c）t=1125ms；（d）t=1500ms；（e）t=1875ms；（f）t=2250ms；
（g）t=2625ms；（h）t=3000ms；（i）t=3375ms；（j）t=3750ms
具体结果见表3.3。
值得说明的一点，1：10模型计算参数
是以1：100试验获得的资料进行试算获得，在1：10模型试验前进行预报，预报的结果
与模型试验的结果误差范围为1.96%~
15.0%；从另外一个角度看，也可以说模型实验结果得到DDA方法的验证，更进一步证明了实验结果的可靠性。因此可以采用DDA方法对三峡工程三期上游RCC围堰爆破的倾倒过程进行数值模拟，为爆破设计提供依据。
3.3.3三峡工程三期上游RCC围堰原型倾倒过程数值模拟参照对围模模型的数值模拟过程，对三峡工程三期上游RCC围堰上部拆除堰体的爆破倾倒过程进行数值模拟，采用二维计算模型进行DDA模拟，计算参数采用1：10模型的计算参数，计算模型的模拟高程50~140m。倾倒过程如图3.5所示。
计算结果：0.8s断裂面出现裂缝，开始翻转；5.17s拆除堰块与下部残留堰块闭合，发生碰撞；12.15s拆除堰体与残留堰体碰撞，13.38s堰块触地；15s拆除堰块完全倒地。
需说明的是：堰块在倾倒过程中出现的翻转，是由于堰块倾倒空间高（模拟堰块底高程为50m）所致，当堰前底高程为70m时不会出现堰块倾倒过程中的翻转现象，这在1：100模型试验和模拟计算中得到验证。
综上所述，三峡工程三期上游RCC围堰爆破拆除的倾倒数值仿真模拟结果有以下几个方面。
（1）通过DDA数值模拟，表明在爆破缺口形成的条件下，上部拆除堰体能顺利倾倒。
（2）采用DDA方法进行数值模拟时，除常规的材料参数外，刚硬弹簧系数P和时间步长△t也为影响倾倒过程的两个重要参数，本次计算过程中发现，P取50倍弹模，At取150us时，计算结果与模型试验结果吻合较好。
（3）计算结果表明：倒地时间为3.06s；倾倒部分与残留体最近距离为0.8m；顶部矩形块与倾倒部分脱离，距离为0.5m。对比模型试验结果，倾倒时间误差范围在
1.96%~15.0%。说明用DDA方法对爆破倾倒过程进行模拟是可行的。
（4）参照围堰模型的模拟过程，采用DDA方法对三峡工程三期上游RCC围堰爆破的倾倒过程进行数值模拟，计算结果为：0.8s断裂面出现裂缝，开始翻转；5.17s拆除堰块与下部残留堰块闭合，发生碰撞；12.15s拆除堰体与残留堰体碰撞，13.38s堰块触地；
15s拆除堰块完全倒地。
3.4ANSYS/LS-DYNA堰块碰撞的振动响应仿真模拟
本节所述的振动响应仿真模拟，旨在进行爆破冲击荷载作用下的有限元应力分析。
（1）围堰1：10模型应力分析。针对模型试验中高程130m（原型高程）以上部分与下部出现分离现象进行数值分析，若将顶部矩形块（高程130m以上部分）与下部拆除堰块之间设为一个弱面，计算发现，顶部将出现明显的上抬、向后翻转情况。
图3.5三峡工程三期上游RCC围堰原型爆破倾倒过程数值模拟
（a）t=2.25s；（b）t=4.5s；（c）t=6.75s；（d）t=9.0s；
（e）t=11.25s；（f）t=13.5s；（g）t=15s采用ANSYS有限元进行爆压作用下堰体的应力状态分析。材料参数：混凝土堰块的弹性模量E=25.5GPa，动弹模提高30%取33.2GPa，泊松比u=0.17，重度2400kg/m3，抗压强度17MPa，抗拉强度1.75MPa；底座钢筋混凝土的弹性模量E=30GPa，动弹模39GPa，泊松比u=0.17，重度2500kg/m3；基础为软岩，弹性模量E=1GPa，动弹模
1.3GPa，泊松比u=0.17，重度2200kg/m3。在高程130m处设置一薄弱单元，材料强度分别取原来的90%、80%、60%、50%。将爆破荷载作用转化为作用在断裂面上的冲击荷载，1号、2号、3号药室作用在切口顶部的爆压为
p，=po是）（3.4）式中po——密闭容器爆压，MPa；ro——药室直径，m；R——药室中心至断裂面的距离，m。
计算得出：1号药室p，=6.4MPa；2号药室p，=21MPa；3号药室p，=7.9MPa。
按《水工建筑物抗震设计规格》（DL5073一2000）规定，采用承载能力极限状态式验算构件强度：
Ya式中yo——结构的重要性系数，由结构安全级别确定；g——设计状况系数；ya——结构系数；S（。）——作用效应函数；R（.）—结构抗力函数（混凝土设计强度，动力计算提高30%）。
为方便计算和比较，令xK=尽2，则XK≥1.0即满足要求。校核结果见表3.4。
Yo9ora
表3.41：10模型下游折坡处应力状态校核结果
①1：10模型RCC抗拉强度取低值。
计算结果表明，在爆破荷载作用下，下游折坡处产生拉应力集中，材料强度降低64.5%以上时，最大拉应力超过材料强度极限值，可能产生破坏。
这和实际模型试验的结果相符，在模型混凝土浇筑的过程中，浇筑高程130m以上部位（换算原型）时，将斜模板换成直立模板，混凝土浇筑时间间隔较长，使得两块混凝土之间形成一个薄弱面，在爆破荷载作用下，顶部与下部堰块脱离并上抬。整个试验过程仅一块发生此现象，概率为25%。
（2）三峡工程三期上游RCC围堰原型应力分析。参照1：10模型的有限元分析，采用ANSYS程序计算分析在爆压作用下三峡工程三期上游RCC围堰原型的应力状态。材料参数：碾压混凝土弹模26GPa，泊松比0.167，重度2350kg/m3，抗压强度19.6~
23.5MPa，抗拉强度1.4~2.35MPa；基础岩石弹模取29GPa。动力计算将弹模和强度提高30%。
在高程130m处设置一薄弱单元，材料强度分别取原来的90%、80%、60%、50%。
同样将爆破荷载作用转化为作用在切口顶部的冲击荷载，各药室的爆压为：1号药室p，=
28.9MPa；2号药室p，=97.6MPa；3号药室p，=36.2MPa。采用承载能力极限状态式验算构件强度，校核结果见表3.5。
表3.5三峡工程三期上游RCC围堰原型下游折坡处应力状态校核结果
①原型碾压混凝土强度取低值。
计算结果表明，在爆破荷载作用下，下游折坡处产生拉应力集中，材料强度降低50%以上时，最大拉应力超过材料强度极限值，可能产生破坏。三峡工程三期上游RCC围堰浇筑质量较好，不会出现高程130mm弱面强度低于设计强度50%，所以在实际爆破时，这种现象将不会发生。
（3）1：10模型爆破冲击振动有限元仿真分析。2005年10月24日和2005年11月8日分别进行了1：10模型爆破倾倒试验，在振动测试中均发现振动的最大峰值不是切口爆破和堰体触地后引起的振动，经分析可能是在爆破切口形成后，上部堰体在爆破冲击荷载作用下上抬下座并撞击保留体后引起的振动，冲击振动的频率小于爆破振动的频率，因而后者对需保护建筑物具有更大的潜在威胁。
采用ANSYS/LSDYNA有限元的接触分析功能，模拟堰块下座冲击残留堰体的振动响应。
1）实测振速。2005年11月8日模型试验时，在堰块底座后缘及堰右前侧底板（距离堰块距离11.5m）各布置了振动测点。振动测试成果见表3.6。
表3.6三峡工程三期上游RCC围堰1：10模型拆除试验质点振动测试成果
注根据起爆顺序，B1堰块为先起爆堰块；B2堰块为后起爆堰块，选取两个具有代表性的测点进行比较。
2）上抬高度估算。根据堰块上抬后下座并撞击振动历时，采用自由落体公式估算堰块上抬的高度，如表3.7所示。
采用DDA非连续变形程序模拟计算堰块上抬的高度，算得上部堰块上抬的高度为19.3cm。由于未充分考虑水介质的影响，DDA计算结果和实测的上抬高度存在误差，其值偏大（除表3.7中A2段外）。
表3.7按自由落体估算1：10模型上抬高度
①A2段由于斜向倾倒，误差较大。
3）计算模型及参数。根据试验观测资料及表3.7中的计算，假定上部拆除堰块在爆破荷载作用下，上抬20cm，下座后与下部残留堰块碰撞产生振动。模型计算中，假定上部堰块从20cm处自由落体，以2m/s的速度与残留堰块碰撞，采用ANSYS/LS-DYNA模拟碰撞后的振动响应。
为简化计算，做如下假设：①材料本构关系为线弹性模型；②岩石与混凝土作为均质材料处理；③认为围堰与堰基础结合紧密，不考虑两者的接触关系。
取4m段堰体进行三维数值计算，基础计算范围取上、下游及底部各延伸15m（约3倍的堰高）。计算模型见图3.6。
图3.6围堰1：10计算模型
材料参数：混凝土堰块，弹性模量E=25.5GPa，动弹模提高30%取33.2GPa，泊松比x=0.17，重度24kN/m2；底座混凝土，弹性模量E=30GPa，动弹模39GPa，泊松比u=0.17，重度25kN/m'；基础为软岩，弹性模量E=1GPa，动弹模1.3GPa，泊松比
=0.17，重度22kN/m3。
计算采用4节点固体单元，总共划分1640个单元，2358个节点。计算时采用坐标系为：X轴为垂直河流方向（近似），Y轴为竖直向上，Z轴为顺河流方向。
边界为无反射边界。无反射边界（NonmRefecting Boundary）又称通射边界（Trans miting Boundary）或无反应力界（Silent Boundary），主要应用于无限体或半无限体中，为减小研究对象的尺寸而采用的边界条件。
水的作用考虑为动水压力，转化为附加质量施加。
4）计策钻果分析。假定上部据块以2m/分的速度与下都堰块撞击，采用ANSYSLS DTNA软件计算，对应于振动测试时测点所布置的部位选取相应点的点的节点，节点745位
F根凝士加雌初常高点0金工基址上最后能凝止5品，书点1号检字提高作解概盖离11.5m，速度时程曲线如图3.7~图3.9所示。
图3.7混凝土底座边缘节点745振速时程曲线（a）水平向质点振速V；（b）竖直向质点振速V，
图3.8堰后距离11.5m基础上节点1087振速时程曲线（a）水平向质点振速V；（b）竖直向质点振速V，
图3.9堰前距离11.5m基础上节点1457振速时程曲线（a）水平向质点振速Vz；（b）竖直向质点振速V，
表3.81：10模型质点最大振速结果表
从表3.8可以看出，计算值与实测值存在误差，且计算值偏大，说明采用数值分析预报原型冲击振动是偏于安全的。
综上所述，堰块碰撞后的振动响应仿真模拟分析结果表明：
1）采用ANSYS/LS-DYNA模拟堰块碰撞后的振动响应，计算结果表明混凝土上抬后下座并撞击引起的振动较大，混凝土底座边缘节点的最大振速为15.1cm/s；软岩基础上，堰后距离11.5m处节点的最大振速为0.77cm/s，堰前距离11.5m处节点的最大振速为1.08cm/s。
2）计算结果和模型实验结果对比，计算值偏大，说明采用数值分析预报原型冲击振动是偏于安全的。
（4）三峡工程三期上游RCC围堰原型爆破振动预测。根据前面的计算成果，先采用DDA模拟计算三峡工程三期上游RCC围堰实际爆破上抬的高度，再采用ANSYS/LS DYNA的接触分析功能预报三峡工程三期上游RCC围堰的爆破振动。
1）DDA模拟。
采用DDA模拟计算三峡工程三期上游RCC围堰上抬的高度，计算参数选用1：10模型计算所采用的参数。断裂孔的爆压见表2.1，将爆破荷载作用转化为作用在断裂面上的冲击荷载，1号、2号、3号药室作用在切口顶部的爆压用式（3.4）计算。
计算得各药室的爆压为：1号药室p，=28.9MPa；2号药室p，=97.6MPa；3号药室p，=36.2MPa。
各药室的延迟时间延长√10=3.16倍（按试验获得的时间比尺），荷载作用时间也相应延长。采用DDA模拟，计算得三峡工程三期上游RCC围堰爆破时的上抬高度为6cm。
2）有限元仿真分析。
①节点振速结果分析。根据5.1节的计算成果，分析三峡工程三期上游RCC围堰爆破时上部混凝土块上抬下座冲击引起的振动。基础岩石弹模取29GPa，动弹模提高30%
取37.7GPa。计算选取厂房坝段一典型断面，大坝轴线和围堰轴线距离114m，大坝混凝土弹模28GPa，抗压强度17.5MPa，抗拉强度1.75MPa。
假定上部堰块上抬高度为6cm，并自由落体，以1.10m/s的速度与下部混凝土块撞击（方案1），采用ANSYS/LS-DYNA软件计算。选取节点位于基础上，节点2407位于大坝坝踵处，速度时程曲线如图3.10所示。
方案1计算结果，大坝坝睡处节点2407的竖直向最大振速为0.72cm/s，水平向最大
图3.10方案1大坝坝踵节点2407振速时程曲线
（a）水平向质点振速V；（b）竖直向质点振速V，振速为0.78cm/s。
为安全起见，计算第二种方案，假定上部堰块上抬高度为12cm，并自由落体以1.55m/s的速度与下部混凝土块撞击（方案2），选取同样的节点速度时程曲线如图3.11所示。
图3.11方案2大坝坝踵节点2407振速时程曲线
（a）水平向质点振速V；（b）竖直向质点振速V，方案2按最不利计算，大坝坝踵处节点2407的竖直向最大振速为1.10cm/s，水平向最大振速为1.01cm/s。两方案结果列于表3.9中。
②应力分析。按最不利计算，上部堰块从12cm高处自由落体与下部混凝土块撞击，三峡厂房坝段在振动作用下，特征部位最大动拉应力结果见表3.10。
从表3.10可知：大坝在冲击振动荷载作用下，坝踵处为拉应力集中区，在坝踵处产生的最大动拉应力为0.18MPa。
表3.9三峡工程三期上游RCC围堰原型质点最大振速结果表单位：cm/s
③应力叠加。对模型进行静力分析，静力分项荷载包括坝体自重、静水压力（水位高程135m）、扬压力。实际爆破荷载也会与冲击荷载叠加，爆破荷载峰值速度按2.48cm/s控制。将3种工况应力结果叠加（动力值乘以0.5后再与静力值叠加），按综合应力进行安全评价。“十”为拉应力，“一”为压应力，计算结果见表3.11。
表3.10三峡厂房坝段特征部位
最大动拉应力结果单位：MPa
表3.11特征部位叠加后应力计算结果单位：MPa
计算得出：按最不利情况进行静动叠加后，在坝踵处产生拉应力，但未超过允许值。
（5）总结。
1）通过DDA方法数值模拟，表明在爆破缺口形成的条件下，上部拆除堰体能顺利倾倒。对1：10模型进行模拟，计算的倒地时间3.06s；倾倒部分与残留体最近距离为
0.8m；顶部矩形块与倾倒部分脱离，距离为0.5m。预报倾倒时间值与模型的试验值误差范围在1.96%~15.0%。说明采用DDA方法对三峡工程三期上游RCC围堰爆破的倾倒过程进行数值模拟是可行的。
2）采用DDA方法对三峡工程三期上游RCC围堰原型爆破的倾倒过程进行数值模拟。计算结果：0.8s断裂面出现裂缝，开始翻转；5.17s拆除堰块与下部残留堰块闭合，发生碰撞；12.15s拆除堰体与残留堰体碰撞，13.38s堰块触地，15s拆除堰块完全倒地。
3）有限元分析表明，在爆破荷载作用下，1：10模型下游折坡处产生拉应力集中，材料强度降低64.5%以上时，最大拉应力超过材料强度极限值，可能产生破坏。RCC围堰原型下游折坡处产生拉应力集中，材料强度降低50%以上时，最大拉应力超过材料强度极限值，可能产生破坏。整个试验过程仅一块发生与下部堰块脱离上抬的现象，概率为25%。三峡工程三期上游RCC围堰浇筑质量较好，不会出现高程130m弱面强度低于设计强度50%，所以在实际爆破时，这种现象将不会发生。
4）假定模型的上部堰块上抬高度为20cm，并自由落体，以2m/s的速度与残留堰块碰撞，采用ANSYS/LS-DYNA模拟碰撞后的振动响应。结果表明堰块上抬后下座并撞击引起的振动较大，混凝土底座边缘节点的最大振速为15.1cm/s；软岩基础上，堰后距离11.5m处节点的最大振速为0.77cm/s，堰前距离11.5m处节点的最大振速为1.08cm/s。
计算结果和模型实验结果对比，计算值偏大，说明采用数值分析预报原型冲击振动是偏于安全的。
5）参照DDA模拟分析结果，采用ANSYS/LS-DYNA的接触分析功能，模拟三峡工程三期上游RCC围堰原型爆破的冲击振动，上部堰块上抬高度为6cm，并自由落体，以1.10m/s的速度与下部混凝土块撞击，大坝坝踵处节点2407的竖直向最大振速为0.72cm/s，水平向最大振速为0.78cm/s；若上部堰块上抬高度为12cm，大坝坝踵处节点2407的竖直向最大振速为1.10cm/s，水平向最大振速为1.01cm/s。
6）大坝在冲击振动荷载（上部堰块上抬高度为12cm，并自由落体）作用下，坝踵处为拉应力集中区，在坝踵处产生的最大动拉应力为0.18MPa。将静力、爆破荷载、冲击荷载三种工况按最不利情况进行应力叠加，结果表明在坝踵处产生拉应力，但未超过允许值。
推荐阅读：
