离散元_边界元动力耦合模型在地下结构动力分析中的应用
中国的cae软件
边界元分析(BEM)技术
边界元分析技术是一种数值计算 方法，适用于解决边界积分方程 等问题。
在中国，边界元分析技术的研究 和应用也得到了广泛的发展。中 国在边界元分析技术方面拥有丰 富的理论和实践经验，并将其应 用于多个领域，如结构力学、流 体动力学、电磁场等领域。
中国在边界元分析技术方面也取 得了很多成果，包括开发了多种 边界元分析软件，如Fesafe、 Elmer等，这些软件在中国的工 程应用中发挥了重要作用。
04
CATALOGUE
中国cae软件的应用实例
cae软件在机械制造行业的应用
模具设计优化
01
利用cae软件对模具设计方案进行模拟和优化，提高模具的精
度和使用寿命。
生产线规划与仿真
02
通过cae软件模拟生产线流程，优化生产线的布局和生产效率
。
机械部件疲劳寿命预测
03
利用cae软件对机械部件的疲劳寿命进行预测，提前发现潜在
离散元(DEM)技术
离散元技术是一种数值计算方法，适用于解决离散系统的动力学问题。
在中国，离散元技术的研究和应用也得到了广泛的发展。中国在离散元 技术方面拥有丰富的理论和实践经验，并将其应用于多个领域，如地质
工程、机械工程、土木工程等领域。
中国在离散元技术方面也取得了很多成果，包括开发了多种离散元分析 软件，如UDEC、PFC等，这些软件在中国的工程应用中发挥了重要作 用。
拓展国际市场
加强国际合作与交流，开拓国际市场，提 升中国CAE软件的国际影响力。
加强国际合作与交流
合作研究
鼓励国内企业和研究机构与国际先进企业与研究机构开展合作 研究。
技术交流
积极参加国际技术交流活动，分享中国CAE软件的技术成果和 发展经验。
离散元粘聚力模型
离散元粘聚力模型
离散元粘聚力模型是一种模拟材料离散化行为的力学模型，通常用于分析岩石、土壤等颗粒材料的力学行为。
离散元粘聚力模型基于以下假设：
1.颗粒之间存在粘聚力，这种粘聚力随着颗粒间距的减小而增
大；
2.颗粒之间存在摩擦力，这种摩擦力与颗粒之间的相对运动方向
和速度有关；
3.颗粒之间存在弹性恢复力，这种恢复力与颗粒之间的相对位移
和速度有关。
离散元粘聚力模型通过引入粘聚力和摩擦力等参数，可以模拟颗粒材料的应力-应变关系、强度、破裂等行为。
该模型对于理解颗粒材料的微观结构和宏观行为之间的关系、优化颗粒材料的设计和应用具有重要的意义。
离散元粘聚力模型的优点包括：
1.可以考虑颗粒之间的相互作用和相互影响，更准确地模拟材料
的力学行为；
2.可以模拟材料的破裂和失稳等复杂行为；
3.可以用于分析不同颗粒材料在不同条件下的力学行为。
离散元粘聚力模型的缺点包括：
1.需要大量的参数和数据，且参数的确定较为复杂；
2.对于某些材料的模拟结果可能存在误差；
3.对于动态行为的模拟可能不够准确。
岩土工程中的数值模拟方法及工程应用
岩土工程中的数值模拟方法及工程应用岩土工程是一门研究土体和岩石在水、力和热的作用下行为特性及其在工程实践中应用的学科。
随着计算机技术的不断发展和应用，数值模拟方法已经成为岩土工程中必不可少的研究手段之一。
本文将从有限元方法、离散元方法和边界元方法三个方面探讨岩土工程中常见的数值模拟方法及其工程应用。
一、有限元方法有限元方法是目前最为广泛应用的岩土工程数值模拟方法之一，其主要特点是可以进行非线性和非平衡的分析。
在岩土工程中，有限元方法主要用于模拟岩土体在受力下的变形和破坏过程。
有限元方法的求解过程可以划分为以下三个步骤：1. 离散化——将复杂的物理问题离散化为条形单元进行计算，使得计算变得简单；2. 建立方程——将有限元模型建立为代数方程组，通过求解方程组得到解；3. 处理结果——利用分析结果来展示研究对象的物理特性和行为。
在岩土工程中，有限元法主要用于地下工程和地震工程等方面的研究，比如隧道围岩和坝体安全评价、塑性材料本构模型细化、岩石三轴试验模拟等。
有限元法的应用使得传统规律模型得以精细化，模拟效果更加接近实际情况。
二、离散元方法离散元方法是一种用离散单元来描述物质状态、分析物质运动的力学方法。
离散元方法是一种适用于多体动力学和岩土体力学问题的数值分析方法。
离散元方法的特点是将物体分解成为微小单元进行数值模拟，从而得到宏观上看起来的结果。
在岩土工程中，离散元方法主要用于土体颗粒流、岩体破坏分析、地震工程模拟等方面的研究。
离散元法常用于研究固体、颗粒和流体的耦合问题，如土石流运动规律研究、软黏土土体力学性质研究等。
三、边界元方法边界元方法，也叫边界积分方法，是一种应用在数学物理问题上的计算算法。
该方法不需要离散化处理，只需要在表面上建立边界元网格即可。
在岩土工程中，边界元方法主要用于颗粒间相互作用、地下水流、地震动等方面的研究。
边界元方法的优点是不需要建立离散网格，仅需在边界上建立少量的节点，计算速度较快，且精度较高，由此常用于模拟地下水流动或地震波传播。
离散元_边界元动力耦合模型
2001年1月水 利 学 报SHU IL I XU EBAO 第1期收稿日期:1999212227基金项目:“九五”国家攻关项目.作者简介:金峰(1966-),男,贵州遵义人,教授、博士生导师,主要从事结构动力分析研究.文章编号:055929350(2001)0120023205离散元2边界元动力耦合模型金 峰1,贾伟伟1,王光纶1(11清华大学水利水电工程系,北京 100084)摘 要:本文提出了一种二维变形体离散元与时域边界元的耦合模型,这一模型可以将非连续体的模拟与无限域的模拟统一在一个模型中,可用于在地震波动输入条件下,考虑辐射阻尼的岩体边坡或地下结构等的动力稳定和变形分析,拓宽了离散元动力分析的领域.算例分析表明本耦合分析模型具有较高的精度.关键词:离散元;边界元;耦合中图分类号:O344 文献标识码:A离散元法是一种模拟离散介质的计算方法,自Cundall 在70年代提出以来,在岩石力学、土力学、结构分析等领域的数值模拟中得到广泛应用,是一种新兴的非连续体分析方法[1,2,3].在动力分析中,许多分析均表明辐射阻尼对分析结果影响很大,应该充分重视[4],动力边界元方法由于在其基本解中包含了无限远处的辐射条件,在处理辐射阻尼的影响时十分方便,是分析地下结构动力响应的一种有力的工具[5].动力边界元法可以分为频域边界元法[5]与时域边界元法[6,7],前者主要适用于线性问题.是发展较早、相对成熟的动力边界元方法.时域边界元方法直接在时域内求解,适合解决非线性问题,根据求解的问题具有非线性的特点,本文选择动力时域边界元方法与可变形体离散元耦合,提出了一个二维时域动力边界元———可变形体离散元耦合模型,充分发挥离散元与边界元的优点,将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中,为地下结构和岩质边坡的抗震稳定分析提供了全新的手段.1 离散元2边界元耦合模型图1 离散体系示意111 二维可变形体离散元原理 可变形体离散元法将模拟的计算区域看作是若干可变形块体的组合,见图1.这些块体可以任意平移、旋转,块体之间的相互作用力,用法向和切向弹簧表示,称为接触,接触力的大小由块体的相对位置决定.每个块体又划分为三角形差分网格以模拟变形,每个差分三角形的顶点称为节点,其应变假定为常数,可以由节点的位置确定.因为一旦节点位置确定,便可以得到差分三角形乃至整个块体的变形和应力,从而得到所有的响应历程,因此,所有的计算将围绕节点进行,边界条件也可以通过给定边界节点位移或节点力来实现.离散元计算采用显式步进的方法,首先将计算的过程分为若干等长的时步,静力分析作为动力分—32—析的特例,可以采用临界阻尼以增加收敛速度,此时静力分析的时步不具有真实的物理意义,可看作迭代步.在每一时步内,对所有块体的所有节点分别进行循环计算,每一节点循环计算的主要步骤为:(1)根据上一循环的结果或边界条件,确定本时步初节点的位置和速度等运动量.(2)根据本节点与相邻节点的位置可确定差分网格的应变,根据本构关系求出应力,再积分得到本节点所受弹性力.若为块体边界上的节点,可根据相邻块体的位置确定接触力,从而得到本节点在本时步所受合力的F ′.以上具体公式和算法可参见文献[1,2,3].(3)由牛顿第二定律或动力平衡方程m ¨u +αm u =F ′+m g(1)式中:F ′是接触力F c 、弹性力F e 与刚度阻尼力F d =β(F c +F e )之和;am u 代表质量阻尼力;α、β是Rayleigh 阻尼系数;g 代表重力加速度,最后一项表示重力的影响.由式(1)通过中心差分,可得到本时步(第i 时步)的速度u i = u i -1(1-015αΔt )+F ′m+g Δt (1+015αΔt )(2)进而可以得到本时步末本节点的位置.返回第(1)步,继续下一轮循环计算,最终求得所有时刻所有节点的解.由于任意块体间在运动过程中都有发生各种接触的可能,在计算中需要对所有的块体之间进行接触判断并计算接触力,这是一项十分耗时的工作,也是离散元法的关键技术之一,为减少计算工作量,针对岩体中构造面的特点,已有一些较好的算法,如Cundall 提出的域算法,充分利用生成离散块体时的信息,将接触检索局限于初始顶点构成的域中,特别适合于由构造面切割形成的离散块体系统,由于篇幅所限,在此不能详细介绍,请参阅有关文献[1]和[2].112 二维时域边界元基本原理 二维全平面时域动力边界元方程[6,7,8]c αβu β(S ,t )=∫t 0∫Γu 3αβp β(Q ,τ)-p 3αβu β(Q ,τ)d Γd τ+u I β(3)式中:u 3αβ,p 3αβ为二维全平面时域动力基本解,其具体的表达式可参见文献[6]和[8];u ,p 分别表示位移与面力;u I 表示从无限远处入射的位移场;S ,Q 分别表示源点及场点.α,β=1,2,同时对时间、空间进行离散,可以得到时域边界元方程可写作[6,7,8]:[H ]L {u}0=[G ]0{p}L +{B }L (4)式中:{u}L 、{p}L 分别是L 时刻的位移与面力分量,{B }L =∑L -1l =1[G ]L -l {p}l -[ H ]L -l {u}l +{u I }L (5)式中:[H ]l 、[G ]l 、(l =1,L )共2L 个矩阵均为系数矩阵,他们可由各边界单元的子矩阵集成而得,具体表达式可参见文献[6]和[8].由于{u}l 、{p}l (l =1,2,…,L -1)是1,2,…,L -1时刻的量,在求解L 时刻时{B }L 为已知量,而在{u}L ,{p}L 中有N (总的自由度数)个值已知,另N 个值未知,因此可解此方程组,可求得L 时刻的N 个未知量,依此类推,即可求得全部时程的历程反应.图2 离散元-边界元耦合模型113 耦合模型的实现 为了实现时域边界元和离散元的耦合,必须解决以下问题:(1)离散元模拟的离散体与边界元分析的连续体实现耦合;(2)求解方式的统一是实现耦合的关键;(3)离散元计算时步与边界元计算时步的协调及离散元节点与边界元单元节点的对应可以进一步提高耦合计算的效率.下面,分别介绍解决的方法.(1)首先,将计算区域划分为两个区域,见图2,即主要考虑无限域辐射阻尼和地震输入的连续体区域,称为边界元域,用时域边界元(B EM )模拟,其模拟的重点是无限域,地表水平自由边界在离散一段距离后截断,经过—42—收敛计算证明可以很好地模拟远域的影响[8].另一区域,称为离散元域,用离散元(DEM)模拟,模拟的重点是离散体.在远域与近域之间设立一个界面块(Interface Block),它是一个完整的可变形离散块体,可进一步划分为三角差分网格,本身是连续的,差分节点同时也可以是边界元节点,边界元与离散元的耦合完全在过渡块上进行,这样,就避免了一个边界元节点可能与两个块体相连而造成的非连续问题.(2)在求解方式上,将整个边界元区域看成是一个离散元法的块体,离散元法计算的核心是根据位置、位移、速度等运动量求作用力,为此,改写边界元方程(4),可以得到{F I}L=[R][G]0-1[H]0{u}L-{B}L(6)式中:[R]是将面力{p}L转换成作用力的转换矩阵.这样,只要知道过渡块中与边界元相接的节点的位移,就可以根据式(6)求出边界元区域作用在这些节点的作用力,进而用式(2)计算节点速度、位移,从而完成一次循环,实现了两种模型的耦合.(3)为了保证计算效率和计算稳定性,时域边界元方法的计算时步不能太小,通常为离散元方法计算时步的几十到上百倍.为解决这一难题,我们采用了异步计算的办法.设离散元法的时步为Δτ,则取时域边界元法的时步Δt=K1K2Δτ,即将每个边界元时步划分为K1个小时步,每个小时步由K2个离散元时步组成,并认为一个小时步内边界元域的作用力不变,既每K2个离散元时步才调用一次式(6)更新作用力.在一个边界元时步内,式(6)的各系数[G]L,[H]L,{B}L皆保持不变,{u}L将不断按下式更新:{u}L={u}L-1+{u}k-{u}L-1K1k (k=1,2,…,K1)(7)式中:{u}k是第k个小时步中第一个离散元时步初的节点位移.当一个边界元时步结束后,式(6)中的面力向量也相应按下式更新{p}L=∑K1k=1[R]-1{F I}kK1(6)进而更新{B},再进行下一个边界元时步的计算.为保证计算效率,边界元的节点通常要大大少于耦合边界上离散元的节点,我们采用一个边界单元对应若干离散元节点的办法,一个边界元内相互作用面力按线性分布,再根据力和力矩的平衡,分别计算分配给每个离散元节点的相互作用力{F I}.这样,通过空间和时间的异步,大大提高了计算效率和改善计算稳定性,从而实现了时域边界元方法和离散元方法的耦合.2 模型验证为验证耦合模型的精度,列出以下两个算例,其他算例可见文献[9].211 岩柱算例 如图3所示的一个岩柱,分为离散元块体模拟的4个岩块和一个边界元区域模拟的顶部岩块,最下一个离散岩块固定,在边界元岩块最上的边界作用单位阶跃荷载110MPa,岩块的弹模E=30GPa,密度ρ=2000kg/m3,阻尼比ξ=011,K n=K s=1011N/m,f=1118,c=0.由于在计算条件下,所有构造面保持完整接触,能保证岩柱的连续性,所以可同时采用时域边界元法对整个岩柱按连续体进行分析.图4为岩柱上各点的位移响应及其与时域边界元法的计算结果对比,在整个计算时间内,两者吻合很好.这一算例显示了当离散块体紧密结合成一个连续体时,耦合模型能够与连续体模型分析结果吻合,并可清楚地看到波动在岩柱中的传播过程,波动能够准确地在边界元区域与离散元区域之间和离散块体之间传播,从一个侧面证实了耦合模型及软件的正确性.—52—图3 岩柱示意图4 耦合模型与时域边界元模型的比较图5 所示的半平面上的离散块体212 半平面的离散块受入射波作用 为进一步验证耦合模型在处理波入射动力问题的计算精度,首先计算了如图5所示的半平面上的离散块体在竖直向上的SV 波入射情况下的响应,离散元域除界面块外仅有一个块体,并且假定离散块体之间的摩擦系数与抗拉强度足够大以保证他们之间不会滑动,边界元域共采用了26个单元,单元长度均为10m ,其中6个单元与界面块体接触,在SV 波入射下的位移历程见图6.图6 输入的位移历程图7 两种方法计算A 点水平位移的比较 图7示出了A 点的水平位移响应与完全采用时域边界元方法计算结果进行的比较.两种方法计算结果的对比再次说明耦合模型具有较高的精度,同时说明了耦合模型能够模拟波动从无限远入射并正确输入到离散元域.其他证明耦合模型能够模拟离散块体的开合等大变形行为的算例,因为篇幅所限和耦合模型与普通离散元模型在这些方面并无本质区别,所以不再列出,可以参考文献[9].致 谢 本文得到国家电力公司成都勘测设计研究院肖白云和王仁坤两位总工的支持和帮助,清华大学水利水电工程系张楚汉教授和徐艳杰博士也给予了大量的帮助,在此一并表示感谢.参 考 文 献:[1] Cundall PA ,Hart RD.Development of G eneralized 22D and 32D Distinct Element Programs for Modelling Joint 2ed Rock [R ].ITASCA Consulting Group ,Misc.Paper SL 28521,1985.[2] 鲁军,张楚汉,王光纶,金峰.岩体动静力稳定分析的三维离散单元法[J ].清华大学学报,1996,(10).[3] Zhang Chuhan ,et al.Application of distinct element method in dynamic analysis of high rock slopes and blockystructures.[J ]S oil Dyn.Earthq.Eng.1997,(12).[4] 张楚汉.结构-地基动力相互作用问题[C].结构与介质相互作用理论及其应用.南京:河海大学出版社,1993.[5] Niwa Y ,K obayashi S ,Azuma N.Application of Integral Equation Method to S ome G eomechanical Problems[C].Numerical Method in G eomechanics ,120-131,ASCE ,New Y ork ,1976.—62—[6] Mansur WJ,Brebbia CA.Topics in Boundary Element Research[C].Chap.4,Springer2Verlag World Pub2lishing Company,87-123,1985.[7] 金峰,张楚汉,王光纶.有阻尼的时域边界元方法[J].力学学报,1997,29(15).[8] 任允涛.各向同性与各向异性介质波动问题边界元法及其工程应用[D].北京:清华大学,1996年.[9] 贾伟伟.离散元-边界元动力耦合模型研究及其工程应用[D].北京:清华大学,1999年.Coupling model of distinct element2bound ary elementJ IN Feng1,J IA Wei2wei1,WAN G Guang2lun1(11Tsi nghua U niversity,Beiji ng 100084,Chi na)Abstract:A22D dynamic model coupling the deformation block distinct element method with boundary element method in time domain is established.The model simulates the static and dy2 namic responses of discontinuous rock and the effects of infinite domain simultaneously.It can be used to analyze the static and dynamic stability,deformation of rocky slopes as well as the under2 ground structures,especially when wave propagation input of earthquake and radiation damping must be taken into acciunt.The results of bench mark problems show shows that the precision of this method is high.K ey w ords:distinct element method;boundary element method;coupling model——72。
混凝土结构多物理场耦合分析方法研究
混凝土结构多物理场耦合分析方法研究一、研究背景混凝土结构是现代建筑中常用的结构材料之一，具有强度高、耐久性好、施工方便等优点。
然而，在实际使用中，混凝土结构会受到多种物理场的作用，如荷载、温度、湿度等，这些物理场的作用会相互耦合，影响混凝土结构的安全性和使用寿命。
因此，混凝土结构多物理场耦合分析方法的研究具有重要的理论和实践意义。
二、研究现状目前，混凝土结构的多物理场耦合分析方法主要包括有限元方法、边界元方法、离散元方法等。
其中，有限元方法是最常用的一种方法，它可以将混凝土结构分为有限个小单元进行分析，建立数学模型，求解各个物理场的相互作用，得到混凝土结构的应力、应变等参数。
边界元方法则是将混凝土结构的边界分为有限个小区域进行分析，求解边界上的物理量，然后利用边界条件得到混凝土结构的应力、应变等参数。
离散元方法则是将混凝土结构分为有限个小颗粒进行分析，求解颗粒间的相互作用，得到混凝土结构的应力、应变等参数。
三、研究内容本研究旨在探讨混凝土结构多物理场耦合分析方法的研究，具体研究内容如下：1.建立混凝土结构多物理场耦合分析的数学模型。
根据混凝土结构受到的物理场和相互作用关系，建立相应的数学模型，包括有限元模型、边界元模型、离散元模型等。
2.求解混凝土结构的应力、应变等参数。
利用数学模型，求解混凝土结构在荷载、温度、湿度等物理场作用下的应力、应变等参数，分析混凝土结构的变形、破坏等情况。
3.分析不同物理场的相互作用对混凝土结构的影响。
研究不同物理场的相互作用对混凝土结构的影响，如荷载和温度、湿度和荷载等因素的相互作用，分析不同情况下混凝土结构的稳定性、安全性等参数。
4.优化混凝土结构设计和维护方案。
根据研究结果，对混凝土结构的设计和维护方案进行优化，提高混凝土结构的耐久性和安全性。
四、研究方法本研究采用有限元方法和边界元方法相结合的方法，建立混凝土结构多物理场耦合分析的数学模型，求解混凝土结构的应力、应变等参数，并分析不同物理场的相互作用对混凝土结构的影响。
基于离散元的流固耦合方法在工程中的应用
基于离散元的流固耦合方法在工程中的应用摘要：该文简要介绍了基于离散元的流固耦合方法的基本原理，在工业生产、工程建设等中的应用，阐述了采用基于离散元的流固耦合方法进行数值模拟在处理实际问题时的可行性和优越性，最后对基于离散元法的流固耦合在工程中的更广泛应用做出了总结和展望。
关键词：离散元流固耦合数值模拟应用离散元法（Discrete Element Method，简称DEM）是分析散体的力学行为的数值方法，在对包含离散颗粒的系统进行数值模拟时更加接近实际情况，模拟结果的参考性就越强。
流固耦合现象广泛存在于自然界和人类实践中，对于流固耦合问题，按照固相和液相相互作用的重叠范围，可分为两大类[1]。
第一大类是两相有部分或全部重叠在一起，例如利用气体动力进行固体颗粒的输送（密相气力输送），化工过程中的流化床问题，地层中水、油、气的运动，空气中粉尘的运动等等；第二大类是两相的耦合作用发生在两相交界面处，结构物在空气、水中与液相的相互作用，例如机翼和气流的相互作用，桥梁结构在风荷载作用下的强度分析，船舶与波浪的相互作用分析等等。
对于第二大类的固相介质都是连续的整体，用离散元模拟将不再具有优势，因此该文讨论的基于离散元的流固耦合方法主要针对第一大类问题，即流固两相部分或完全重叠的实际应用。
1 基于离散元的流固耦合方法原理简介在流固系统中，固体是用离散的圆形(或球形）颗粒来模拟的，根据牛顿第二定律的力与位移的关系，采用显式时步循环运算规则，对数值模型进行循环计算，重复应用此定律于颗粒上，颗粒间的接触力通过接触模型得到。
对于流体的描述采用连续性方程和Navier-Stokes方程，假设流体为不可压缩流体且密度不变。
在模型的计算中，流体对颗粒体的作用力以外部体力的方式施加于球体颗粒，它包括两部分，一部分是流体的压力梯度产生的，另一部分是流固之间的相互作用力，该部分的计算采用经验公式求出。
在流体计算模块中采用的是由Patankar和Spalding[2]最先提出的SIMPLE计算方法(Semi-Implicit Method for Pressure Linked Equations)来求解每一个流体单元中流体的压力和速度矢量，SIMPLE被广泛用于不可压缩流体的求解。
离散元法及其应用简介
Fluided bed (Yutaka Tsuji)
Particle Transport in Electrostatic Traveling Wave (Hiroyuki Kawamoto)
陶瓷材料二维BPM模型的建立
生成颗粒
颗粒位置重排
删除“浮点”颗粒
Cluster模型
BPM模型
陶瓷材料DEM模型的校准
模拟残余应力
v＝10m/s、ap＝5μm
v＝5m/s、ap＝15μm
v＝10m/s、ap＝10μm
v＝10m/s、ap＝15μm
v＝10m/s、ap＝15μm 不同切削深度下的残余应力云图
v＝15m/s、ap＝15μm 不同切削速度下的残余应力云图
混凝土泵送
在混凝土工程施工过程中，混凝土的运输和浇筑是一项繁重、关键性的
混凝土泵送
GID建模 区域离散与 网格划分 计算及结果 可视化
建模过程
泵送过程的三维压力分布
典型的泵车弯管破损形式
陶瓷压实成型
重力情况下
-210 -205
微重力情况下
Gravity Microravity
20
Gravity Microgravity
-195
Distribution (%)
0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
Yuanqiang Tan, Hao Zhang, Dongmin Yang et al. Numerical simulation of concrete pumping process and investigation of wear mechanism of the piping wall. Tribology International, 2011 (In Press).
混凝土力学数值仿真研究进展 (修改)
混凝土力学性能数值仿真研究进展吉久茂，焦楚杰（广州大学土木工程学院，广东广州 510006）摘要：回顾了目前研究混凝土数值仿真的常用方法，分别是有限元法，离散元法，边界元法及它们的发展过程。
综述了数值仿真方法的基本原理和应用模型。
侧重叙述了颗粒离散元法和散体细观力学理论及其DEM模型与改进。
介绍了各种算法国内外研究进展及耦合算法的扩展研究情况。
最后指出了一些数值仿真算法的存在问题，并对数值仿真的未来发展即科研实验与数值仿真结合研究提出了希望。
关键词：混凝土数值仿真有限元离散元边界元1．数值仿真的起源与科研意义数值仿真技术也叫计算机模拟。
它是以电子计算机为手段通过数值计算和图像显示的方法，达到对工程问题及物理问题研究的目的。
数值模拟技术最早诞生于20世纪中叶。
由国外知名学者Bruce G.H和PeaceD.W模拟了一维气相不稳定径向和线形流。
1955年，Peace man与Rach ford 研发的交替隐式解法（ADI）取得了数值仿真技术的重大突破。
该算法稳定且计算速度较快，所以在相关领域得到广泛应用。
1959年，Peace man与Douglas第一次进行了两维两相模拟，这标志着现代数值模拟技术的开始。
60年代，Coats K.H和Nielsen R.L首次进行了三维两相模拟。
1968年Breitenbach E.A发表了三维三相模拟解法，这期间的另一项成就就是Peace man 提出了后来通用的Peace man方程，用于二维三维扩散方程的数值解答。
80年代最大的成就是Appleyyard J R和Cheshire I.M发表了嵌套因式分解法。
90年代，Zoltan E.Heinemann提出了PEBI网格，PEBI网格结合了正交网格和角点网格的优点，现在正逐渐成为主流数值模拟网格体系。
关于数值模拟，我们可以理解为利用计算机做实验，特别是近几十年来随着计算机性能的不断提高和力学学科的深入发展，数值模拟方法对问题的分析研究取得重大进展。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2001年2月水 利 学 报SHU IL I XU EBAO 第2期收稿日期:1999212227基金项目:“九・五”国家攻关项目.作者简介:金 峰(1966-),男,贵州遵义人,教授,博士生导师,主要从事结构动力分析.文章编号:055929350(2001)022*******离散元2边界元动力耦合模型在地下结构动力分析中的应用金 峰1,王光纶1,贾伟伟1(11清华大学水利水电工程系,北京 100084)摘 要:利用离散元2边界元耦合模型对溪洛渡工程地下厂房洞室群静、动力响应进行了分析,对离散元模拟地下结构的一些问题进行了探讨,得到了一些结论:离散元模型能够很好地模拟地下洞室群的变形,通过与边界元的耦合,可以模拟辐射阻尼的影响,由于辐射阻尼的影响,地下洞室群的地震响应比地面结构小很多.关键词:离散元;辐射阻尼;地下结构;动力分析中图分类号:TV554 文献标识码:A离散元2边界元动力耦合模型[1]将二维可变形体离散元[2]与二维动力时域边界元[3]两种计算方法耦合起来,充分发挥离散元与边界元的优点,将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中,为地下结构的抗震稳定分析提供了全新的手段.本文利用这一耦合模型对溪洛渡工程地下厂房洞室群进行了静、动力分析,对离散元模拟地下结构的一些问题进行了探讨,得到了一些结论.1 工程简况溪洛渡水电站在长江上游干流金沙江上,是一座以发电为主,兼有拦沙、防洪、漂木、航运等综合效益的巨型水电站.根据国家电力公司成都勘测设计研究院的溪洛渡地下厂房区地质概况报告显示,坝址区位于溪洛渡峡谷中段,山体浑厚,谷坡陡峻,地形完整.经国家地震局烈度评定委员会审定,地震烈度为Ⅷ度.厂房发电系统拟定为库区内地下厂房,左、右岸对称布置,且地下主厂房、主变室和尾水调压室均呈“品”字形分布,由上游至下游依次是地下厂房、主变室和尾水调压室.两岸地下厂房尺寸规模相同,跨度32m ,长414m ,高7417m ,底板高程33515m.厂房区岩体由含斑玄武岩和玄武质角砾熔岩组成,岩体一般新鲜较完整,无大的断层切割,层间错动带一般不发育,产状平缓,挤压紧密,一般干燥无水.层内错动带在局部层内相对较发育,错动带一般宽5～10cm ,挤压紧密,为岩块岩屑型.裂隙以陡倾和缓倾为主,中等倾角裂隙一般不发育.厂房地区地下水活动微弱,厂房区地应力属中等地应力区.2 计算条件和简化假定211 计算断面的选取 限于时间和条件,目前开发的离散元-边界元动力耦合分析软件只能进行二维模拟,因此必须选取适当的二维计算断面.以右岸厂房为例,初步选定通过12号机组中心线的厂房断面,得到一个向两侧无限延伸的断面,如图1所示.然后分别再划定边界元和离散元的离散范围.离散元的离散范围为以洞室群为基准,向下延伸长度大于最大洞室高度,向两侧延伸约两倍最大洞室高度,即图中折线ABCDEA 所包围的区域.折线FABCD G 以下的半无限平面则采用边界元方法来模拟,并假定为各向同性的连续介质,AF 和D G 均为水平线.根据计算经验,各用15个70m 的边界单元离散,并在F 和G 处截断.折线ABCD 用24个边界单元离散,每个单元长度大约为70m ,整个地表FAED G 的边界条件均为自由表面.图1 简化后的计算断面图2 计算采用的离散元网格 计算中对洞室断面作了相应的简化,对于采用离散元模拟的ABCDEA 区域,必须进行块体的划分和材料基本参数的确定.根据成勘院提供的工程地质条件,首先确定选取断面内主要裂隙的参数,根据这些参数可以由计算机自动生成裂隙并以此来对计算区域进行块体划分.岩体中实际存在的裂隙往往是间断性的,存在岩桥,而目前离散元的研究对象必须是被裂隙完全切割的各种形状的块体.因此必须对实际的裂隙进行必要的处理,将较短的裂隙进行延长,使它们相交,将整个岩块比较彻底地切割成各种块体,以便进行计算.为保证如此处理后仍能模拟原有的岩体,将根据裂隙延长的长度对计算的力学参数如摩擦系数f 、粘滞力c 、接触弹性系数K n 、K s 等进行加权修正,以保证等效的抗滑力、弹性力相同.最后形成的离散块体见图2,图2中所示阴影部分表示的块体为耦合模型分析所需的界面块[1].为清晰起见,图2中未示出各离散元块体内部的差分网格,在洞室群附近差分网格的大小均小于4m ,离洞室群稍远处网格的大小均小于10m ,上部山体的网格大小均小于40m ,整个计算模拟的块体数大于500块.整个完成一次静、动力计算在Pentium Ⅲ/450上约需80～100h.212 地应力及其模拟方法 地应力是影响围岩稳定的一个重要因素.溪洛渡地下厂房区存在着较高的初始地应力,在计算模型中应该充分予以考虑.在采用耦合模型进行数值分析中,针对溪洛渡工程的具体情况,可以在竖直方向,通过增大重力加速度,即施加一个大于g 的重力加速度,相当于给岩体施加附加的竖向荷载,有效地模拟竖向地应力,再通过施加水平力来挤紧岩块,来模拟水平方向的应力,这样,通过调整施加的重力加速度及水平力的大小,可以使垂直地应力和水平地应力均达到预期要求的数值.本文主要考虑地下厂房洞室群的稳定性,因此计算中主要考虑洞室群及其附近区域的地应力场.实际的地下厂房洞室群地应力场数据已经由国家电力公司成都勘测设计研究院提供,它是根据实测的地应力加以反演推得的.经过试算,计算模型中在水平方向共施加415×106kN 压力,在竖直方向施加1115g 的重力加速度.图3 简化的施工开挖顺序213 渗流场的模拟 水库蓄水后,水位抬高,会在两岸山体内形成渗流场.地下洞室开挖以后,水会由洞壁渗漏到洞室中,地下水将在洞室周围形成局部汇流区.地下渗流场的存在对洞室的稳定会产生一定的影响,因此计算模型中必须考虑渗流场.耦合模型中主要考虑静水头的影响,即根据各个节点的静水头确定作用在节点上的渗透压力,目前,静水头暂时按渗流面考虑,未能充分模拟洞室本身的排水效应.当然,将来也可按渗流等压面考虑静水头.地震动水压力的影响由于时间所限未予考虑,蓄水前,地下水位很低,水库蓄水后,地下水位上升,洞室才开始受渗透水压力的影响.因此,渗透压力的施加是在开挖完成,且离散元系统计算达到新的稳定以后.214 地下厂房系统开挖顺序 由于耦合模型分析是仿真分析,施工开挖顺序会对分析结果产生较大影响,计算中必须模拟施工开挖的顺序.在施工过程中,溪洛渡地下厂房系统的开挖有一定的先后顺序,须按开挖施工进度表的安排进行.在主厂房开挖的过程中,同时还要进行主变室和调压室的开挖及支护.但实际的开挖施工进度安排比较复杂,且本研究的重点又是动力稳定,故在计算中进行了一定的简化,将洞室群的开挖简化为三步进行,如图3所示,其中Ⅰ、Ⅱ、Ⅲ分别代表第一、二、三步开挖的部分.计算时每一步开挖完成后给予系统一定的时间,使之达到稳定状态,以便进行下一步计算.3 施工期与蓄水期分析311 模拟过程 离散元计算的初始状态假定块体间没有相互作用力.计算开始以后,在自重和地应力(以及可能的地下水渗透压力)作用下,块体之间会产生相互作用力并引起变形,直至达到平衡状态,可称为“固结过程”.因此,采用耦合模型进行计算时,首先要经过一个固结过程,直到系统到达稳定后,才开始施工模拟计算.施工过程分三步开挖来模拟,每步开挖均要在达到稳定后再进行下一步开挖.施工开挖模拟结束后,再施加渗透压力荷载,模拟蓄水过程,为简化起见,蓄水过程假定为一步. 312 变形分析 取一典型a点(见图3),耦合模型计算的该点在施工过程中的水平位移曲线如图4 (a)所示(图中的时间没有实际意义),可以分为三个阶段.ab段为从计算初始状态开始的稳定阶段(固结阶段),即达到施工前初始地应力的状态;bc阶段为开挖阶段,共分三次完成开挖;cd阶段为蓄水阶段.图4(b)为图4(a)的bd段局部放大图,表示了每一步洞室开挖以后,a点发生的回弹变形.由于该点位于洞室侧壁的中部,第一步只开挖洞室的上部,因此第一步a点的回弹位移较小,为013mm;第二步开挖到洞室的中部,这时a点产生了较大的回弹位移,为1014mm;第三步开挖洞室的下部,这一步开挖量较大,a点产生了很大的回弹位移,大小为2915mm.水库蓄水以后,由于上游侧水位高于下游侧,岩体受到渗压的作用,也会产生新的位移,运动趋势为向下游方向.图4 典型的位移曲线其它各点也有类似规律,值得注意的是主变室的位移.它的开挖在第二步一次完成,它的回弹位移也一次完成,规律同主厂房和调压室.但在第三步开挖主厂房和调压室的下部时,受调压室开挖的影响,它会产生向下游的水平位移.取其洞壁上的b点(见图3)进行分析,在第二步开挖中水平回弹位移大小为417mm,由于调压室所影响,它在开挖第三步产生了向下游方向的位移,大小为2010mm,影响是很大的.因此洞室群的开挖要十分注意洞室开挖的相互影响.3个洞室因开挖产生的总的回弹位移大小如图5所示,最大位移发生在调压室的侧壁上,最大回弹的位移值为53mm.另外主厂房下游侧底部也存在向下游的变形,这也是主变室和调压室开挖的影响造成的.4 地震响应分析在开挖与蓄水过程静力分析的基础上,利用离散元2边界元耦合模型分析了地震波波动输入条件下,各洞室的地震响应.地震荷载采用从无限远处波动输入,峰值加速度是采用溪洛度工程的地震危险性分析提供的峰值加速度0118g,地震加速度反应历程由抗震规范规定的标准反应谱反演得到,地震历时24146s,见图6,通过数值积分并进行归零调整得到地震位移历程,波动输入的幅值按该地震波在半平面表面的响应为0118g反演,在SV波竖向输入时,在无限远处输入的地震波动幅值为0109g,此时的地震输入与通常的水平地震输入等价.图5 洞室群开挖回弹总变形图6 设计地震加速度历程411 地表地震响应 首先分析在地震波输入时,山体表面的地震响应.计算结果显示山体地表各点的峰值加速度略有不同,图1中山顶E 点附近的加速度较大,峰值加速度为1194m/s 2,为0118g 的111倍;河谷底部A 点附近的加速度较小,峰值加速度为1136m/s 2,为0118g 的0177倍.这一分布与大量地震实际观测的现象吻合.地表峰值加速度平均值为1163m/s 2,为0118g 的01925倍,略小于0118g ,这与离散块体的阻尼和地震分布有关,总的来说,地表的地震加速度与设计地震相当.图7 耦合模型计算的a点地震加速度图9 耦合模型分析的加速度与设计加速度比值分布图8 常规离散元分析的加速度与设计加速度比值分布412 洞室周边各点的地震响应 图7为耦合模型分析得到的a 点(见图3)地震加速度响应历程,它是由计算得到的位移历程通过数值微分得到的,可以注意到在25s 附近有一较大峰值,这是由于地震时程只有24146s 记录,突然终止造成的虚假响应.其他各点的地震加速度与之相似.图8和图9分别是采用离散元模型与耦合模型计算的各洞室周边各点的加速度峰值与设计地震峰值的比值分布图,其中,离散元的计算网格与图2中耦合模型的离散元网格系统,仅是没有与边界元耦合,界面块体为刚性块,地震输入是在这个刚性块上给定设计地震来施加的,因此它们之间的差别显示了辐射阻尼和地震波动输入的影响.从图中可以分析出以下结论:(1)常规的离散元分析中,各洞室的地震响应与输入地震相近,最大响应的放大倍数为1164,发生在调压室中下部.(2)由于辐射阻尼和地震波动输入,耦合模型计算的各洞室的地震响应与地表响应与常规离散元分析的地震响应相比有非常明显的下降,最大响应仅有设计地震的01357倍,这与实际地震中地下结构地震反应较小的实际观测结果相吻合,说明在溪洛渡地下厂房洞室群这样埋深较大的地下结构,考虑辐射阻尼和地震波动输入能够更好地符合实际情况.(3)耦合模型计算中,由于地震波从下向上入射,洞室底部是迎波面,加速度响应较大,而顶拱的加速度响应较小.残余位移的分布呈上大下小的分布,最大残余位移发生在调压室顶部,最大残余位移值为10mm.地震残余位移较小,对洞室群的围岩稳定威胁较小.413 地震波不同入射方向的影响 考虑到地震波的输入方向受到多种因素的影响,难以准确确定,从工程安全的角度出发,对地震波的不同入射方向进行了分析,给出了SV波与竖直方向夹角为30°和-30°入射的结果.可以分析出以下结论:(1)不同入射角工况条件下,洞室周边的地震响应略有不同.(2)在各种入射角条件下,一般迎波面的响应比背波面大,在地震从下向上的入射条件下,洞室底部的加速度响应较大,而顶拱的加速度响应较小.(3)残余位移的分布仍呈上大下小的分布.各种输入条件下的残余位移略有不同,其中30°输入条件下,不仅加速度响应较大,残余位移也大,最大值为29mm,也发生在调压室顶部,对地下洞室群的整体围岩稳定威胁不大.5 主要结论采用离散元-边界元动力耦合模型对溪洛渡工程地下厂房洞室群的静、动力响应进行了分析,得出以下主要结论:(1)开挖过程中,由于应力回弹作用,地下厂房会发生一定的变形.变形的最大值发生在调压室侧壁的中段,最大值为53mm.(2)开挖过程中,三个洞室之间会产生相互的影响,因此在实际施工时要合理安排开挖顺序.(3)利用上述模型分析溪洛渡工程地下厂房洞室群,得出了地下洞室群地震响应明显小于地表响应的结论,与地震实际观测的一般结论吻合.(4)不同的地震入射方向,对计算结果有一定的影响,在目前的条件下,建议在工程分析中应该计算可能的不同入射方向,选取最危险的情况进行设计.(5)溪洛渡工程地下厂房洞室群的围岩完整性较好,具有较高的抗震稳定能力,在设计地震作用下发生围岩整体失稳破坏的可能性不大.致谢 本文的研究工作得到了国家电力公司成都勘测设计研究院肖白云和王仁坤两位高工的支持和帮助,清华大学水利水电工程系张楚汉教授和徐艳杰博士也给予了大量的帮助,在此一并表示感谢.参 考 文 献:[1] 金峰,王光纶,贾伟伟.离散元2边界元动力耦合模型[J].水利学报,2001,(1).[2] Cundall P A.The Measurment and Analysis of Acceleration in Rock Slopes[D].Ph. D.Thesis,Universityof London,Imperial College Science and Technology,1971.[3] Mansur W J,Brebbia C A.Topics in Boundary Element Research(Ed.by Brebbia, C. A.)[M].Springer2Verlag World Publishing Company,1985.Application of distinct element2bound ary element coupling modelin underground structure dynamic analysisJ IN Feng1,WAN G Guang2lun1,J IA Wei2wei1(11Tsi nghua U niversity,Beiji ng 100084,Chi na)Abstract:Employing the coupling model of distinct element2boundary element,the static anddynamic analysis of underground power plant of Xiluodu Project is conducted.Some problems indistinct element simulation of underground structure are discussed.The conclusion that the dis2tinct element method can simulate the deformation of underground structures is drawn.The ef2fects of radiation damping can also be modeled,if distinct elements can be coupled with boundaryelements.The seismic responses of underground structures will be much smaller than that ofstructures on the ground surface,due to effects of radiation damping.K ey w ords:distinct element;boundary element;coupling model;underground structure;dy2namic analysis。