离散单元法（精选九篇）
离散单元法 篇1
1 沥青混合料三维模型
1.1 级配球单元
沥青混合料中集料为不规则的多面体,由于计算机运算能力及多面体的接触判断分析还比较困难[2],而且研究表明材料的宏观特性诸如变形、应力等可以通过调整颗粒的密实度与界面弹簧的变形与强度特性,由刚性颗粒间的滑移、分离以及颗粒转动来获得较好的近似,无需精确模拟颗粒自身的变形[3]。下面以AC16为例说明集料与级配球单元的转化关系。
假设19.0,…,2.36,1.18号筛对应i=1,…,6,7,则单元的体积为:
由表1知集料的质量剩余率为P1=0,P2=5,P3=11,P4=14,P5=22,P6=14,P7=9.5,第i号筛剩余集料单元的个数可以表示为:
若生成1 000个球单元,其单元个数如表1所示。
1.2 模型建立
为了给球单元设定生成边界并且在模型模拟实验中给球单元施加约束力,本文创建了墙单元——圆柱壁面和正方形壁面,在墙单元内部随机生成球单元,根据沥青混合料胶浆理论及离散单元法中的本构关系,在球单元间设置平行粘结模型模拟沥青胶浆的作用。球单元间发生相对运动时,弹簧间可以承受力和力矩。
通过施加一定的速度给上底墙单元采用控制空隙率的方法成型模型,模型成型后设置墙体单元的速度为零,经过一定步长的循环消除球单元间的不平衡力。成型后模型如图1所示。
2 模型抗压强度实验
根据T0713-2000公路工程沥青及沥青混合料实验规程中沥青混合料单轴压缩试验(圆柱体法)的试验方法[4],通过给上底墙单元一定的速度压缩模型模拟试验中的轴向力,此时要反复调整上下墙单元的速度及墙单元的刚度;当墙单元刚度过大时,球单元不易穿透墙单元,当墙单元速度过小时,计算模型的效率过低。反复调整墙单元的刚度及速度,当模型的空隙率为0.30时[5],球单元从墙体中穿出,认为模型承受的荷载达到最大的轴向应力,此时模型已破坏。监测模型的抗压强度如图2所示。
根据模型抗压强度实验中求得的最大轴向应力,将其分为10级荷载,分别取0.1σ,…,0.7σ七级荷载作为模拟实验的荷载。给上底墙单元施加一定的速度并通过一定的迭代循环给模型加载到第一级荷载,并以相同的速度循环卸载回零,循环一定的次数,监控加载第一级荷载时模型的轴向应变(ΔL1),依次进行第1,…,7级荷载的加卸载过程。
依据公式:
其中,σ5,ΔL5分别为第五级荷载下的轴向应力及相应的应变,由实验室测得沥青混合料的回弹模量一般在195 MPa～1 917.4 MPa之间[6],本文由离散元模型模拟测得模型回弹模量Es=1 625 MPa,故本文建立的模型比较符合实际。
3 结语
1)通过沥青混合料与模型球单元间转化可以控制模型中不同球单元的大小及体积,可以较好的对不同级配的沥青混合料进行模拟,而且模型较接近实际情况[2]。
2)由于模型开始破坏时只是个别球单元从墙单元穿出,破坏不易监控,又因为模型中球单元只是模拟粗集料,故可以采用控制模型空隙率作为模型成型及模型破坏的判断标准。
参考文献
[1]W.G.B.Z.You.Discrete Element Modeling to Predict the Mod-ulus of Asphalt Concrete Mixtures[J].Journal of Materials inCivil Engineering,2004,16(2):140-145.
[2]周长红.沥青混合料非连续力学计算模型的研究[D].大连:大连理工大学,2008.
[3]张楚汉,金风.岩土和混凝土离散—接触—断裂分析[M].北京:清华大学出版社,2008.
[4]JIG E20-2011,公路工程沥青及沥青混合料实验规程[S].
[5]刘玉.沥青混合料细观结构模型的离散单元数值模拟[D].郑州:河南工业大学,2005.
离散数学一单元 篇2
一单元测试题
1.将下列命题翻译成符号逻辑形式
（1）银行利率一降低，股价随之上扬。
（2）尽管银行利率降低，股价却没有上升。
（3）占据空间的、有质量而且不断变化的对象称为物质
（4）如果一个整数能背6整出，那么它就能被2或3整除。如果一个整数能被3整
除，那么它的各位数字之和也能被3整除。
2.判断下面各语句是否是命题，如果是命题，说出它的真值。
（1）可导的实函数都是连续函数。
（2）凡是都有例外。
（3）白天比夜晚时间长
（4）两个三角形全等当且仅当它们的对应角相等。
3.简述命题的定义。
4.简述原子命题的定义。
5.下列公式中，（）不是永真式。（单选，写清楚每个属于什么公式）
A.(P∧Q)→QB.P→(P∨Q)
C．（P→Q）↔（~Q→~P）D.(~P∨Q)∧(~(~P∨~Q))
5.下列语句，是命题的有（）（多选）
1）美国的首都是纽约。2）你喜欢日本吗？3）我们一定要解放台湾！
4）所有实数都是整数。3）如果3>2，那么有人不死。
6.构造公式的真值表，判断哪些是永真式，矛盾式，和可满足式
（1）(P→(Q→R))↔((P∧Q)→R)
（2）(P∧(P∧Q))↔~P
（3）~（P∨Q）→R
7.如果P∨QQ∨R，能否判断PR？如果P∧QR∧Q,能否判断PR？如果~P~R能否判断PR。
8．判断下面等式是否是等价式：P→（Q∨R）（P→Q）∨（P→R）
9．求下列两式的对偶式
（1）（P∧~Q）∨（R∧T）∨F
（2）~（P∨~（Q∨R））∧（R∧~Q）
10．分别利用真值表法和等价变换法求下列公式的主合取范式及主析取范式。
（1）P→（R∧（Q→P））
（2）（P→（Q∧R））∧（~P→（~Q∧~R））
11．证明（P→Q）∧（Q→R）P→R
12．证明R→S是{P→（Q→S），~R∨P，Q}的逻辑结果（使用直接法，CP规则法，和反证法）
13．求公式（P→（R∨P））∧（Q ↔P）的主合取范式和主析取范式。
离散单元法 篇3
关键词：弹簧元法 离散元法 弹簧刚度 边坡工程
Abstract：Aiming at the continuous-discontinuous failure process of rock and soil materials in slope engineering，a novel deformable block discrete element method which combined spring element method（SEM） and discrete element method（DEM）together is presented. Compared with the accustomed element in traditional finite element method（FEM），the element in SEM is described as a spring system that contained some orthogonal generalized springs.This generalized springs are defined in 3D space， which means that each spring can has two or three spring stiffness.How to determine the generalized spring stiffness for continuous material is the difficult and most important in SEM.With the triangle element as an example，the basic theory of SEM is introduced in detail.Assuming the relationship between the generalized spring deformation and the element strain，the generalized spring stiffness can be obtained directly by comparing the elastic strain energy of the element and the elastic potential energy of the spring system.The Poisson and shear stiffness coefficients were defined as system parameters to describe the relationship between different generalized springs.The SEM can consider the Poisson effect accurately for any Poisson’s ratio material；and the result using SEM are the same with using traditional FEM.This method does not need to know the expression of the element-stiffness-matrix. It can be used in 4-node element；and the stiffness expressions of springs are given clearly. With the SEM used to compute the block deformation and the contact-spring used to calculate the interaction between blocks，the combined SEM/DEM program can be used to simulate the failure process of rock and soil material from continuous to discontinuous.The SEM and DEM are combined in the motion equation of each node in each element.The contact-spring in DEM satisfied specific strength criterion.When the contact-spring force exceeded its limit，the material became discontinuous from continuous. The combined SEM/DEM program is implemented easily in the continuum-based discrete element method（CDEM）program.The simulation of a homogeneous soil slope under gravity shows that the SEM is performed as good as FEM when using line elastic constitutive and reasonable when using Mohr-Coulomb strength criterion. The simulation of a bedrock and overburden layer slope shows that the combined program is suitable to simulate the slope failure process.
Key Words：Spring element method（SEM）；Discrete element method（DEM）；Spring stiffness；Slope engineering
离散单元法 篇4
李海峰[4]、Shin[5]等研究者建立了熔融气化炉移动床区域的一维数学模型,数学模型主要包括气-固-液的质量、能量守恒方程,研究得到移动床区域炉气-炉料的温度分布,以及炉气的组分分布。王凤[6,7]在研究中将炉料简化为多孔介质,在没有考虑温度场的情况下,模拟熔融气化炉内的气流分布。王成善[8,9]借鉴高炉数学模型的成功经验,建立了稳态下移动床区域热过程的二维数学模型,数学模型主要包括气固液各相的质量、动量及能量守恒方程。通过对熔融气化炉中移动床区域的数值模拟,得到气固相的速度场及温度场; 并研究不同布料方式、移动床区域的高度以及死料柱对各组分的温度、速度、浓度分布规律的影响。
从目前的研究现状分析可知,研究者们主要基于连续介质模型,对移动床区域的热过程以及炉气的流动过程进行研究。然而,在连续介质模型中,将离散的炉料颗粒处理为拟流体时存在一些问题[10,11]: 颗粒相本构方程难以确定、无法进行微观层次上的研究。由于拟流体模型在研究炉料流动中存在不足之处,在高炉炉料流动的研究中,离散单元法( DEM) 广泛的被采用。Zhou研究了模型高炉内炉料颗粒的运动[12],炉料流动的DEM模拟结果与实验结果吻合的较好; 为了减少计算时间,加快了炉内颗粒的流动速度; 研究发现随着炉料消耗速率的增加,死料柱的几何尺寸变小。由于磨损导致焦炭颗粒变小,从而影响高炉的透气性; 因此为了弄清焦炭的磨损机理,Taihei[13]基于EDM的模拟高炉内炉料的应力分布; 研究颗粒之间的应力分布能有效的决定高炉内焦炭应该达到的强度; 同时在考虑焦炭强度的情况下,能采用合理的加料方式。Allert基于DEM-CFD模型研究了非球形颗粒和气体流量对小尺度高炉模型中炉料流动的影响[14,15]。与炉料采用球形颗粒相比,采用非圆球形颗粒时,死料柱几何尺寸稍微的变大。
从以上分析可知,基于DEM模型对熔融气化炉移动床区域炉料流动过程的研究较少。在本文中,对移动床区域炉料流动的进行DEM模拟,经过研究得到炉料的流动规律,以及炉内炉料颗粒之间的法向力分布。
1 炉料颗粒流动的DEM模型
1. 1 颗粒运动方程
离散单元法( DEM) 的基本思想是把不连续散体分离为相互独立的单元,根据单元之间的相互作用和牛顿运动定律,计算每一个时间步长内所有单元的受力和位移,通过对每个单元的微观运动进行跟踪计算,可得到整个研究对象的宏观及微观运动规律。颗粒的运动主要包括线性运动和转动,由下面的方程描述。
式中m、v、I、#、r分别为颗粒的质量( kg) 、速度( m/s) 、转动惯量( kg·m2) 、角速度( rad /s) 、半径( m) ;ni为颗粒i的质心指向接触点的法向单位向量; Fgi、Fcij分别为颗粒所受到的体积力( N) ,颗粒ij之间的接触力( N) 。
1. 2 颗粒间接触力
在DEM中,颗粒间的作用力用弹簧-阻尼器-摩擦器表示[16]。在该模型中,弹簧模拟颗粒的形变,阻尼器模拟碰撞过程中的阻尼影响,摩擦器模拟颗粒间的滑动摩擦力。采用Hertz-Mindlin( no slip) 模型计算颗粒之间的接触力[17],接触力主要包括法向力、切向力以及摩擦力。
1. 3 摩擦系数
炉料颗粒的摩擦系数包括静摩擦系数和滚动摩擦系数,通过大量的文献调研,炉料颗粒之间的静摩擦系数的取值约为0. 5[18—21]。DEM模型采用球形颗粒模拟移动床区域内不规则颗粒的运动,并不能有效的模拟颗粒形状因素对炉料运动行为产生的影响。在DEM模型中,通过增加滚动摩擦系数来模拟形状因素对炉料运动的影响,然而滚动摩擦系数对移动床区域内炉料运动过程的影响还不清楚,本节将对此进行研究。
通过炉料自然堆积过程的DEM模拟,得到不同滚动摩擦系数下熔融气化炉炉料的自然堆积角度,如图2 和图3 所示,在静摩擦系数值为0. 5 时,随着滚动摩擦系数的增加,炉料的自然堆积高度和角度增加。当滚动摩擦系数在0. 01 ~ 0. 2 之间时,随着滚动摩擦系数的增加,堆积角度近似线性的增加,当滚动摩擦系数从0. 2 增加到0. 4 时,堆积角度从45°增加到47°,堆积角度的变化较小。参考高炉炉料的自然堆积角度在40° ~ 45°之间[22],选取滚动摩擦系数0. 15 作为炉料运动DEM模型的参数。
得到DEM模拟的参数如表1 所示。
2 计算区域
图4 为熔融气化炉移动床区域流料流动的计算区域,实际物理模型呈圆台形状。由于对称性以及为了减少计算量,取移动床区域的一半作为计算区域。炉料颗粒在移动床表面连续的加入,风口回旋区处炉料颗粒不断的消耗,从而形成一个缓慢下降的移动床区域。
3 计算结果及分析
3. 1 炉料流动规律及DEM模型的验证
得到炉料流动的DEM模拟结果,如图5 所示。为了得到料层下降过程的变化规律,在移动床表面加上一层物性相同颜色不同的颗粒作为示踪颗粒层,得到不同时刻移动床区域中料层的形状见图5( a) 。为了更好的分析炉料的运动,提取示踪颗粒层的位置数据放在同一张图中,即炉料流动的时间线,如图5( b) 所示。
图6 为炉料流动的实验研究结果。从移动床区域中炉料流动的DEM模拟结果与实验研究结果可知,示踪颗粒层的流动形态和时间线呈现出相似的变化规律。随着料层下降,示踪颗粒层逐渐发生变形,靠近壁面的示踪颗粒层下降速度较慢; 风口回旋区以上炉料下降速度较大,风口回旋区和壁面对示踪颗粒层的影响逐渐变得更加明显; 示踪颗粒层逐渐呈现出“V”形状。从以上分析可知,DEM模型的计算结果与实验中炉料的运动规律基本上相同,从而验证了DEM模型计算结果的正确性与合理性,DEM模型能够准确的模拟移动床区域炉料的运动过程。
在DEM模型中采用实际炉料颗粒的粒径及物性,物理模型尺寸为原型的1 /10; 实验研究中的实验模型为原型的1 /20。DEM模拟的物理模型以及实验模型的几何形状与原型移动床区域的几何形状相似。由于DEM模拟的计算量大,在计算中为了减少计算量,加快了炉料的下降速度[12]; 因此当示踪颗粒层下降到相同位置时,DEM模拟中的时间远远小于实验中所用时间。
从图5 可知,根据炉料流动速度,将移动床区域分为四个不同的子区域。区域①( 径向方向r < 350mm,高度方向h > 300 mm) ,即稳定流动区域: 料层的水平下降,向下流动速度呈现相对均匀的速度,料层之间并没有出现相互掺混的情况。在区域②( r >350 mm,h > 300 mm) ,即壁面附近区域: 由于壁面摩擦力和阻力的影响,料层的流动速度逐渐减小,并且料层颗粒出现明显的变形和相互掺混。在区域③( r < 350 mm,h < 300 mm) ,即死料柱区域; t = 500 s和t = 550 s两层示踪颗粒层几乎重合在一起,说明示踪颗粒层的向下流动速度非常小,下降到该区域时几乎停止下降运动,形成一个弧形示踪颗粒层即死料柱的轮廓。死料柱呈现为半椭圆形,高度约为210 mm,宽度约为390 mm。在区域④( r > 350 mm,h < 300 mm) ,即活跃区域: 死料柱的存在以及熔融气化炉半径的减小,使炉料的流动区域变窄,导致该区域炉料流动速度的加快。活跃区的速度明显高于其它区域。从t = 400 ~ 550 s的时间线可以看出,示踪颗粒层呈现出“V”字形状,“V”字形料层的最下端靠近风口回旋区。在活跃区域内料层的向下流动速度较快,壁面和死料柱的影响使靠近壁面以及死料柱之上的料层向下流动速度减慢,因此呈现出“V”字型的料层形状。
3. 2 移动床区域内颗粒之间法向力分布
炉料流动的实验研究,并不能获得移动床区域内炉料颗粒之间的法向力分布,而炉料颗粒之间的法向力是影响颗粒破碎的因素之一。因此通过炉料流动的DEM模拟,得到移动床区域内炉料的应力分布,如图7 所示。在死料柱中及死料柱以上区域中,颗粒之间的法向力较大。在稳定流动区域中和活跃区中,颗粒之间的法向力较小。随着靠近壁面附近区域,颗粒之间的法向力逐渐增加。
图8 为移动床区域内颗粒之间的平均法向应力分布。图8( a) 为径向方向r =0 ~50 mm及r =400 ~450 mm区域内,在高度方向上颗粒之间法向力分布。在r = 0 ~ 50 mm区域内,随着高度的降低,颗粒之间法向力逐渐增加; 在高度h = 600 ~ 700 mm之间,炉料颗粒之间的法向力较小约为120 N; 在h =210 ~ 700 mm范围内,随着炉料靠近死料柱,炉料颗粒之间的法向力增加到700 N; 在h = 0 ~ 210 mm之间,处于死料柱区域内,颗粒之间的法向力增加约为800 N。在r = 400 ~ 450 mm区域内,在h = 600 ~700 mm范围之间,炉料颗粒之间的法向力约为120N; h = 400 ~ 600 mm范围靠近壁面,由于壁面的阻力的影响,炉料流动速度减小,颗粒之间法向力增加约为200 N; 在h = 0 ~ 400 mm之间,随着炉料进入活跃区,炉料颗粒的快速运动使炉料之间的接触力较弱,颗粒之间颗粒之间法向力降低为150 N。
在死料柱区域内,由于死料柱内的颗粒几乎静止不动,该区域炉料颗粒持续承受的重力最大,导致颗粒之间法向力较大; 该区域的炉料颗粒起着支撑整个移动床的作用,因此该区域的炉料颗粒应该达到足够的强度; 在实际生产中加入的块煤转变为半焦,其强度并不大而容易劣化成小颗粒,导致死料柱区域的孔隙率非常小; 气流经过该区域的阻力非常大导致进入该区域的炉气量就会非常少,造成气流分布的均匀性较差; 同时高温炉气不能顺畅的流出移动床区域,炉气聚集在风口回旋区及附近区域,使该区域炉气温度升高,导致风口烧损等一系列问题,影响生产的稳定。从图5( a) 中示踪颗粒层运动轨迹可知,中心区域的料层最终下降到死料柱及其附近区域,因此在中心区域多加入一些高强度焦炭,焦炭最终会下降到死料柱及其附近区域,从而提高该区域炉料的强度和该区域的空隙率; 这将会提高整个移动床区域气流分布的均匀性,有利于生产的稳定。
图8( b) 为在高度方向h = 0 ~ 50 mm、h = 350 ~400 mm及h = 650 ~ 700 mm区域内,在径向方向上颗粒之间法向力分布。在h = 650 ~ 700 mm区域内,即处于移动床上部区域,炉料颗粒之间的法向力分布保持不变且较小约为120 N。随着炉料的下降,在区域h = 350 ~ 400 mm中,径向方向r = 0 ~200 mm范围处于死料柱上部附近区域,该范围内颗粒之间法向力较大,约为500 N; 在r = 200 ~ 400 mm范围内,随着远离死料柱,颗粒之间法向力逐渐减小为200 N; 在r = 400 ~ 580 mm范围内,随着靠近壁面,颗粒之间的法向力增加到350 N。在区域h =0 ~ 50 mm中,即处于移动床底部,r = 0 ~ 300 mm范围处于死料柱内部,颗粒之间法向力较大,为600 ~850 N; r = 300 ~ 450 mm范围内,随着靠近活跃区,颗粒之间法向力急剧减小为160 N。
从以上分析可知,死料柱区域内颗粒之间法向力最大,为700 ~ 850 N; 其次为壁面附近区域,为300 ~ 400 N; 稳定流动区和活跃区内颗粒之间法向力最小,约为150 N; 在稳定流动区中,随着靠近死料柱,颗粒之间法向力逐渐增加。
3. 3 风口回旋区炉料消耗速率对应力分布的影响
图9 为不同风口回旋区炉料消耗速率时炉内颗粒之间法向力分布规律,随着风口回旋区炉料消耗速率的增加,炉内颗粒之间法向力分布规律并没有发生明显的变化。不同风口回旋区炉料消耗速率时都呈现出: 死料柱区域法向力最大,其次为壁面附近区域,稳定流动区和活跃区内颗粒之间法向力最小;在稳定流动区中,随着靠近死料柱,颗粒之间法向力逐渐增加。因此在任何条件下,都应该在中心加入高强度焦炭颗粒,以保证死料柱区域炉料的强度和提高该区域的空隙率,从而有利于整个移动床区域气流分布的均匀性。
4 结论
对熔融气化炉移动床区域炉料流动过程进行了DEM模拟,模拟了得到移动床区域炉料的运动和颗粒之间法向力的分布规律。得到以下主要结论:
( 1) 建立移动床区域炉料流动的DEM模型,通过炉料自然堆积过程的DEM模拟,确定了颗粒之间滚动摩擦系数为0. 15。DEM模拟和实验研究的炉料运动规律基本上相同,验证了模型计算结果的正确性与合理性。
( 2) 根据炉料流动速度,移动床区域可分为四个子区域: 稳定流动区域、壁面附近区域、死料柱区域以及活跃区域。在移动床区域的中心区域形成半椭圆形状的死料柱,高度约为210 mm,宽度约为390 mm。
( 3) r = 0 ~ 50 mm区域内,随着高度的降低并靠近死料柱区域,颗粒之间法向力逐渐从120 N增加到700 N,在死料柱区域内颗粒之间法向力约为800 N。死料柱区域内颗粒之间法向力最大,为700 ~850 N; 其次为壁面附近区域,为300 ~ 400 N; 稳定流动区和活跃区内颗粒之间法向力最小,约为150 N在稳定流动区中,随着靠近死料柱,颗粒之间法向力逐渐增加。风口回旋区炉料消耗速率对炉内颗粒之间法向力分布影响较小。
在死料柱中及其以上的区域,料层法向应力最大,因此在中心区域多加入一些高强度焦炭,焦炭最终会下降到死料柱及其附近区域,从而提高该区域炉料的强度和该区域的空隙率,有利于整个移动床区域气流分布的均匀性。
摘要：建立熔融还原炼铁(coal reduction extreme,COREX)工艺熔融气化炉移动床区域炉料流动离散单元法(discrete element method,DEM)模型。通过炉料自然堆积过程的DEM模拟,确定了固体颗粒之间的滚动摩擦系数。炉料流动的DEM模拟结果与试验研究结果吻合的较好。研究了移动床区域炉料的流动规律、死料柱形状及尺寸以及颗粒间法向应力的分布规律,同时分析了风口回旋区炉料消耗速率对应力分布的影响。研究结果表明:熔融气化炉移动床区域可分为四个不同的子区域,死料柱区域法向力最大,其次为壁面附近区域,稳定流动区和活跃区内颗粒之间法向力最小。
离散元法模拟砂浆流变性研究 篇5
关键词：砂浆,剪切屈服强度,离散元法
1 离散元模型
离散元方法的基本假定是模拟材料由离散的颗粒组成, 颗粒本身是刚性的, 颗粒间相互作用为一种动态的平衡, 即内部力一直处于平衡状态。每个模拟单元的运动遵循牛顿第二定律和转矩方程, 因此颗粒的位移和旋转量可根据下列控制方程获得:
其中, Fi为颗粒间接触力矢量;Fin为法向接触矢量;Fis为切向接触矢量;m为颗粒质量;为颗粒平移加速度矢量;gi为重力加速度;Mi为作用于颗粒上的合力矩;I为转动惯量;6) ωi为角加速度。
2 模型细观参数与宏观参数的关系
2.1 流变参数测量方法
Chu于1997年提出一种简易测量水泥浆、砂浆和新拌混凝土流变参数的方法———拖拽式粘度计。其工作原理:将已知半径的钢球置于宾汉姆流体中并施加外力使其以速度v匀速上升, 此时剪切速率γ和剪切应力τ可通过式 (4) , 式 (5) 计算:
通过改变速度v获得剪切速率γ及剪切应力τ关系曲线, 得到对应的宾汉姆方程:
其中, γ为剪切速率;v为球体的速率;r为球体半径;τ为剪切应力;F为宾汉姆介质对球体的作用力;τ0为剪切屈服值;η为粘度系数。
2.2 接触模型细观参数的确定
2.2.1 粒径和密度
新拌砂浆可看作由水泥浆和细骨料拌和而成, 细骨料均匀分布于水泥浆中, 并且细骨料的最大粒径为5 mm。因此将新拌砂浆简化为单一粒径6 mm的球形单元。
为了使简化后的砂浆单元集合质量与实际砂浆的质量相同, 需要对砂浆单元的密度进行修正, 采用等效砂浆单元密度:
其中, ρ'为等效砂浆单元密度;ρ为实际砂浆密度;ε为砂浆模型单元的孔隙率。单一粒径球体的堆积孔隙率受到其空间排列结构的影响, 本文中砂浆模拟单元孔隙率为单元集合体在自重作用下的孔隙率。经过多次的数值模拟, 孔隙率取值为0.4。
2.2.2 泊松比和剪切模量
砂浆在宏观尺度上表现为一种均匀连续的材料, 可以看作是一种不可压缩流体。因此在离散元模拟中, 新拌砂浆模拟单元的泊松比ν取值为0.5。对于砂浆的杨氏模量或者是剪切模量的取值现阶段还没有成熟的理论。本文剪切模量取值参考黄绵松[1]的研究, 考虑时间步长与模拟单元剪切模量的反比关系以及模拟单元在自重作用下的变形量, 通过大量的模拟实验最终剪切模量G取值为106Pa。
2.2.3 Parallel-Bond接触模型参数
2.3 模型计算及结果分析
测试砂浆剪切屈服值模型如图1所示, 模拟单元在自重作用下均匀分布于边长为10 cm的容器中。由于颗粒堆积存在边界条件, 必须考虑粮仓效应, 即当颗粒体堆积高度达到一定值后, 容器底部所受压强不随颗粒的增加而增加, 与液体在容器中所表现的性质完全不同。为了减小粮仓效应影响使颗粒体内部压强与流体静力学压强一致, 建立模型时必须考虑颗粒堆积高度与容器底部尺寸的比值。图2为颗粒堆积模型内部力链网络分布图, 黑色力链随堆积深度增加变粗, 同等深度处力链分布均匀, 因此粮仓效应可不予考虑。
通过控制单因素变量获得宏观剪切屈服值和细观接触模型参数之间的定量关系结果如图3~图6所示。图3显示新拌砂浆的剪切屈服值与模拟单元密度成线性关系。从离散元方法解释砂浆的剪切屈服值, 其本质是颗粒间相互作用的屈服, 而每个模拟单元都遵循牛顿第二定律式 (1) , 即颗粒间接触力与模拟单元密度呈线性关系。因此剪切屈服值随模拟单元密度线性增加符合内在机理。
图4, 图5显示剪切屈服强度与平行粘结切向强度成线性关系, 而与摩擦系数成分段线性关系。这是因为砂浆模型初始为受压状态, 如图2所示只存在黑色压力链不存在红色拉力链, 平行粘结强度对模型结构无影响, 因此成线性关系, 然而摩擦系数的改变会对模型颗粒的空间分布及受力状态产生影响。摩擦力与摩擦系数成正比, 当摩擦系数较小时, 颗粒在自重状态下既可克服摩擦, 砂浆堆积模型结构较密实, 摩擦系数对剪切屈服值的影响如图5中A区域所示。随着摩擦系数增大, 摩擦力达到一定值时, 部分颗粒无法在自重作用下克服摩擦, 结构较疏松, 摩擦对剪切屈服强度影响如图5中B区域所示。
本文砂浆模型引入法向强度来模拟净浆的粘结作用, 其对剪切屈服值的影响如图6所示。显然法向强度对屈服值的影响远小于切向强度。新拌砂浆剪切屈服强度主要由摩擦力和平行粘结强度两部分组成, 其中摩擦力与颗粒堆积体内部应力有关。颗粒体内部压强与流体静力学压强一致, 随砂浆深度线性增加, 因此剪切屈服值公式中摩擦力部分引入砂浆深度参数H来表征其随砂浆深度线性变化的特征。砂浆宏观参数 (屈服值) 与细观参数 (质量、平行粘结强度和摩擦系数) 之间的关系公式:
其中, τ0为剪切屈服值;H为砂浆深度;ρ'为砂浆模型单元密度;为切向强度;f为摩擦系数。
3 对比试验
本文通过砂浆流动度试验来验证砂浆离散元模型及其剪切屈服值与细观接触模型参数之间定量关系的准确性。试验选取两种细度模数的人工砂, 改变水灰比配置不同稠度的砂浆, 并使用NXS-11B型旋转粘度仪分别测得各组砂浆的剪切屈服值, 同时对各组砂浆进行流动度试验, 如图7a) 所示, 各组试验数据见表1。之后, 根据粘度计测得的剪切屈服值及式 (8) 、式 (9) 选择适当砂浆模型细观参数, 建立砂浆离散元模型计算获得相对应流动度, 如图7b) 所示, 最后对模型计算与试验结果进行比较。由于砂浆模型内部剪切屈服值不一致, 通过试算选取模型中心, 即砂浆深度为3 cm处屈服值来表征整个模型剪切屈服值。由图8可知物理试验及模型计算测得的流动度存在一致性。
4 结语
本文运用PFC3D软件建立新拌砂浆离散模型, 并通过模拟拖拽式粘度计试验获得细观接触模型参数与宏观流变参数 (剪切屈服值) 的定量关系, 其关系式表明砂浆的剪切屈服强度主要由摩擦和平行粘结切向强度两部分组成, 这与砂浆流变行为的本质 (即剪切变形) 是相一致的。
拖拽式粘度计模拟试验显示剪切屈服值随砂浆深度线性增加, 其原因在于组成剪切屈服值的摩擦力与砂浆离散元模型内部应力有关。砂浆颗粒体内部压强与流体静力学压强一致, 即压强随深度线性增大, 从而导致摩擦力与砂浆深度成正比。
新拌砂浆流动度试验与模型计算结果的一致性验证了砂浆离散元模型及细观接触模型参数与剪切屈服强度定量关系式的正确性。
参考文献
离散单元法 篇6
随着电力市场化的进行,现今电力系统运行的一个重要任务就是在保证全网电压质量的前提下,使系统运行更经济。而无功功率的合理分布可以有效地减少网损,提高电网的运行经济性。传统的无功优化算法如线性规划、非线性规划、动态规划和混合整数规划都不适于处理离散变量,且由于局部收敛常常达不到全局最优点。近年来,一些人工智能方法如遗传算法、Tabu搜索方法、混沌优化方法等[1,2,3,4]在电力系统无功优化中的研究与应用也取得了一定的进展。
针对离散粒子群与非线性内点法方法的特点,本文提出了一种新的混合策略来求解无功优化问题。该策略主要分为两部分:首先不考虑无功优化中的离散约束,采用内点法求解原问题的初始解;然后将原无功优化问题分解为连续优化和离散优化两个子问题,分别采用内点法和离散粒子群算法交替求解,使两者的优化结果互为基础、相互利用,从而保证了本文混合策略的整体寻优效率。通过IEEE 30和IEEE118节点系统的仿真计算,表明本文混合策略在处理离散无功优化方面具有明显的优势。
1 无功优化的数学模型
电力系统无功优化(ORPF)的数学模型可表示为:
式中:f为系统有功网损最小的目标函数;g为系统潮流约束;Z=[X,UC,UD]为系统变量,其中,X为系统状态变量(负荷节点电压幅值和发电机注入无功功率);UC为连续控制变量(发电机节点电压);UD为离散控制变量(无功补偿装置的无功补偿容量和可调变压器分接头);Zmin和Zmax为系统变量的运行限制约束。由于电力系统无功优化大多是在满足支路潮流不等式约束的情况下进行的,优化后的支路潮流变化不会很大。如果初始条件比较恶劣,考虑到支路潮流不等式约束的复杂性,通常可以把不等式约束作为一个惩罚因子,放入目标函数中。本文为了简单明了,不考虑带有惩罚因子的目标函数。
2 非线性内点法[5,6,7]
非线性内点法直接求解连续非线性规划问题的主要优点是计算时间对问题的规模不敏感,所具有的多项式时间复杂性在计算大规模非线性问题时很有优势。设非线性规划问题
对式(2)所描述的问题,可先引入松弛变量了l和u,将目标函数f(x)改造为障碍函数,得到优化问题B:
其中:扰动因子μ>0。当l或u靠近边界时,以上函数趋于无穷大,因此满足以上障碍目标函数的极小解只能在可行域内找到。这样通过目标函数的变换把含有不等式限制的优化问题变成只含等式限制的优化问题,因此可以直接用拉各朗日乘子法求解。
优化问题的拉各朗日函数
该问题极小值存在的必要条件是拉各朗日函数(4)对所有变量及乘子的偏导数为0。
用牛顿法解上述KKT条件,得到线性化的修正方程组为:
求解方程组(5)得到迭代修正量,不断的迭代更新,从而得到近似最优解。
定义对偶间隙:
Fiacco和Mc Cornic证明在一定的条件下,如果x*是优化问题A的最优解,那么当Gap→0时,产生的序列{xi}能收敛至x*。
3 粒子群优化算法简述
3.1 基本粒子群算法
粒子群优化算法最初是由Kennedy和Eberhart博士[8]于1995年受人工生命研究结果启发,在模拟鸟群觅食过程中的迁徙和群集行为时提出的一种基于群体智能的演化计算技术。假设在N维空间中有m个粒子组成一个群,待优化问题的变量数决定了解空间的维数。任何一个粒子被看成是在N维空间里的点Xi=(Xi1,Xi2,,Xin),第i个粒子迭代到目前为止最好的位置(即具有最佳适应度)被称为个体最好粒子,记为pBesti,而在全部粒子迭代到当前为止具有最佳适应度的粒子则被称为全局最好粒子标记为gBesti,粒子位置变化的速度记为Vi=(Vi1,Vi2,,Vin)。各粒子根据式(7)更新自己的速度和位置。
式中:k表示第k次迭代;c1,c2为学习因子;r1,r2是[0,1]间的随机数;σ为控制速度权重的约束因子;wk为惯性权重。进化过程中按式(2)线性减少[7]:
式中:itermax为最大进化代数;wmin和wmax为wk取值的范围。使得PSO在开始时探索较大的区域,较快地定位最优解的大致位置,随着ω逐渐减小,粒子速度减慢,开始精细的局部搜索(这里ω类似于模拟退火中的温度参数)。该方法加快了收敛速度,提高了PSO算法的性能。
3.2 离散粒子群优化算法
Clerc提出了一种针对离散变量优化问题的离散粒子群优化(DPSO)算法[9]。DPSO算法中,粒子也是通过追踪两个极值粒子进行迭代寻优的。然而,离散变量优化问题的搜索空间和可行解均为离散值构成的集合,针对DPSO算法,由于篇幅问题,本文仅给出其基本公式详细请参考文献[10],如下:
DPSO算法基本公式如下:
上述定义的运算算子可以满足DPSO算法公式的需要。由于DPSO算法在结构上与PSO算法具有统一性,可知DPSO算法与PSO算法具有一致的迭代寻优原理。
4 混合型PSO算法在无功优化中的应用
MPSO求解无功优化步骤:
针对内点法和离散粒子群算法的特点,本文提出了一种有效的混合策略:先忽略离散约束,将ORPF转化为一个连续非线性规划问题,采用内点法直接求解;在此初始解的基础上,根据变量的性质,将ORPF分解为连续优化和离散优化2个子问题,分别采用内点法和离散粒子群法交替求解,直到获得满意的优化结果。在连续优化子问题中,以离散优化结果为基础,假定电容器、电抗器的无功补偿容量和有载调压变压器分接头不变,控制变量为发电机端电压,该子问题相当于离散控制变量为常数,相应不含其上、下限和离散约束。在离散优化子问题中,以连续优化结果为基础,假定发电机端电压不变,控制变量为电容器、电抗器的无功补偿容量和有载调压变压器分接头,该子问题又相当于连续控制变量为常数,相应不含其上、下限约束。
新型混合策略的具体步骤可以表述如下:
MPSO算法求解无功优化问题,可按以下步骤进行:
a)输入原始数据,获取系统节点信息和支路信息,获取控制变量的个数及各自的取值范围,获取粒子群的群体规模m等参数。
b)当迭代次数k=0时,在可行域内随机生成m个粒子的位置xi,每个粒子有Nt+Nc维。其中,变压器变比和电容器容量排在前Nt+Nc维,组成的集合生成初始的离散解值集,同时在一定范围内生成每个粒子相对应的初始速度vi。
c)将上一步结果作为内点法的初始条件调整发电机的出力,进一步对全网无功进行优化。
d)计算各粒子的适应度值,将离散的整数值换算成实际的数值,然后一起代入潮流计算公式(2)和目标函数式(1)计算获得F(xi)值。
e)令所有粒子适应度值F(xi)中的最小值为Fmin,Fmin
f)终止条件判断——迭代终止条件根据具体问题一般选为最大迭代次数或粒子群迄今为止搜索到的最优位置满足预定最小适应阈值。判断是否达到终止条件,若是,则输出计算结果,否则k=k+1,转步骤d。
5 系统算例
本文采用IEEE30节点系统和IEEE118节点系统作为试验系统。为便于比较分析,我们采用三种优化方法:(1)本文所介绍的混合策略算法;(2)把离散变量进行连续化处理并用非线性原对偶内点法求解的算法;(3)在方法二所求得的优化结果基础上,将离散变量的优化值就近归整,再进行潮流计算。运用Matlab7.0进行编程计算的结果如下:
1)IEEE30节点系统
该节点系统包括6台发电机、4台有载调压变压器、4个无功补偿点,节点和支路数据参见文献[12]。优化结果如表1、表2。
由表2可知,连续优化法计算得到的网损为0.045 8,是理想值,比本文算法要小,且各变量也在约束范围之内,但是在实际调度运行中是无法操作的;简单归整法使网损变为0.045 9,但是导致某些节电电压越限,如:27节点的电压计算值为1.100 7,超出了节点电压的上限,这是不可行的;本文算法的网损为0.050 5,与连续优化法接近,而且离散变量的值能较好地满足离散分级变化的要求,同时各变量均在约束范围之内。因此,本文所提出的方法在解决电力系统离散无功优化问题方面还是比较有效的。
2)IEEE118节点系统
该节点系统包括36台发电机、10个无功补偿节点和8台可调变压器,表3和表4列出了运用各种算法优化后的离散变量值和网络损耗值。限于篇幅,连续变量值略写。
由表3、表4知,对IEEE118节点系统运用方法一和方法二优化后,网络损耗值接近,且所有的节点电压、发电机和无功补偿装置的无功出力均在约束范围内,但是方法一满足离散分级变化的要求,更加符合实际情况。方法三求得的网损最小,但有12个节点的电压和6个发电机的无功出力超出上下限范围,显然是不可行解。可见,本文提出的无功优化混合策略较传统的无功优化算法能更合理有效地处理混合整数无功优化问题。
6 结论
在优化计算中离散粒子群算法易于处理离散变量,而内点法在求解大规模非线性连续优化问题时,具有收敛速度快的突出优点。本文提出的混合策略充分利用了这两种算法的优势,采用离散粒子群算法处理离散变量,采用内点法处理连续变量,使二者互为基础,相互利用,从而大幅度提高了混合策略的寻优效率。IEEE30和IEEE118节点系统计算结果表明,本文混合策略实用性强,性能稳定,能有效地求解无功优化这一类大规模混合整数非线性规划问题。
摘要：基于离散粒子群优化算法与内点法,提出了一种新颖的混合策略来求解电力系统无功优化问题:不考虑无功优化中的离散约束,采用内点法求解得到初始解;根据优化变量的不同性质将无功优化问题分解为离散优化和连续优化2个子问题,并采用离散粒子群优化算法和内点法交替求解,使两者的优化结果互为基础、相互利用,从而保证了混合策略的整体寻优效率。以IEEE30和IEEE118节点作为试验系统,与常规的离散优化算法做比较,验证了该算法的正确性和有效性。
关键词：电力系统,无功优化,混合整数规划,内点法,离散粒子群算法
参考文献
[1]张勇军,任震,李邦峰.电力系统无功优化调度研究综述[J].电网技术,2005,29(2):50-56.ZHANG Yong-jun,REN Zhen,LI Bang-feng.Survey on Optimal Reactive Power Dispatch of Power Systems[J].Power System Technology,2005,29(2):50-56.
[2]马晋韬,Lai L L,杨以涵.遗传算法在电力系统无功优化中的应用[J].中国电机工程学报,1995,15(5):347-353.MA Jin-tao,Lai L L,YANG Yi-han.Application of Genetic Algorithm in Reactive Power Optimization[J].Proceedings of the CSEE,1995,15(5):347-353.
[3]王洪章,熊信艮,吴耀武.基于改进Tabu搜索算法的电力系统无功优[J].电网技术,2002,26(1):15-18.WANG Hong-zhang,XIONG Xin-yin,WU Yao-wu.Power System Reactive Power optimization Based on Modified Tabu Search Algorithm[J].Power System Technology,2002,26(1):15-18.
[4]赵涛,雄信银,吴耀武.基于混沌优化算法的电力系统无功优化[J].继电器,2003,31(3):20-25.ZHAO Tao,XIONG Xin-yin,WU Yao-wu.Reactive Power Optimization of Power System Based on Chaos Optimization Algorithm[J].Relay,2003,31(3):20-25.
[5]程莹,刘明波.含离散控制变量的大规模电力系统无功优化[J].中国电机工程学报,2002,22(5):54-60.CHENG Ying,LIU Ming-bo.Reactive-Power Optimiza-tion of Large-scale Power Systems with Discrete Control Variables[J].Proceedings of the CSEE,2002,22(5):54-60.
[6]韦化,李滨,杭乃善,等.大规模水-火电力系统最优潮流的现代内点算法实现[J].中国电机工程学报,2003,23(6):13-18.WEI Hua,LI Bin,HANG Nai-shan,et al.An Implementation of Interior Point Algorithm for Large-scale Hydro-thermal Optimal Power Flow Problems[J].Proceedings of the CSEE,2003,23(6):13-18.
[7]LIU Ming-bo,Tao S K.An Extended Nonlinear Primal-dual Interior Point Algorithm for Reactive Power Optimization of Large-scale Power Systems with Discrete Control Variables[J].IEEE Trans on Power Systems,2002,17(4):982-991.
[8]Kennedy J,Eberhart R.Particle Swarm Optimization[A].in:IEEE International Conference on Neural Network[C].Perth(Australia):1995.
[9]Maurice Clerc.Discrete Particle Swarm Optimization[EB/OL].http://www.Mauriceclerc.net.2001.
离散单元法 篇7
PID控制器问世至今已有近70 年历史, 它以其结构简单、稳定性好、工作可靠、调整方便而成为工业控制的主要技术之一。当被控对象的结构和参数不能完全掌握, 或得不到精确的数学模型时, 控制理论的其它技术难以采用时, 系统控制器的结构和参数必须依据经验和现场调试来确定, 这时应用PID控制技术最为方便。 即当我们不完全了解一个系统和被控对象﹐或不能通过有效的测量手段来获得系统参数时, 最适合用PID控制技术。
随着数字化的发展, 对PID控制器的各种离散方法的层出不穷, 各种离散方法都有本身的优点和局限性, 本文就其采用双线性变换法离散存在的问题进行分析, 从而得出PID控制器不能采用双线性变换法离散的原因所在。
1 PID控制器及双线性变换法
具有比例、 积分和微分控制规律的控制称为PID控制器。 调节PID的参数, 可实现在系统稳定的前提下, 兼顾系统的带载能力和抗扰能力, 同时, 在PID控制器中引入积分项, 系统增加了一个零极点, 使之成为一阶或一阶以上的系统, 这样系统阶跃响应的稳态误差就为零。 PID控制器的传递函数为:
式中, Kp为比例系数;Ti为积分时间常数;Td为微分时间常数。
双线性变换法也叫做梯形积分法, 或称作突斯汀 (Tustin) 变换法。用近似, 也就是用梯形面积之和近似代替数值积分。其离散公式为:
映射关系为:
其特点是:如果D (s) 稳定, 则相应的D (z) 也稳定。 所得D (z) 的频率响应在低频段与D (s) 频率响应相近, 而在高频段相对于D (s) 的频率响应有严重畸变[1]。
2 数字PID控制器的设计
现选取PID控制器参数:Kp=1, Ti=1, Td=0.1, 则:
被控对象为:
采用闭环控制, 其采样控制系统如图1 所示。
在连续域控制系统中, 系统的开环传递函数为:
得到系统的闭环传递函数为:
当输入信号为单位阶跃信号r (t) =1 (t) 时, 闭环系统的输出响应如图2 所示:
该系统是一个I型系统, 所以当输入为阶跃信号时, 系统的稳态误差为0。 由于PID控制器中的积分环节I为系统提供一个位于原点的开环极点, 提高系统型别保证系统稳定性。系统的动态性能良好, 超调量 σ%=0, 调节时间ts=3.4s, 上升时间tr=2.9s。
以下将采用双线性变换法的连续域—离散化方法, 设计数字PID控制器。 并通过在闭环系统的输出曲线, 分析PID控制器采用双线性法进行离散存在的问题。 设 ωc为连续系统的截止频率, 根据式
可求出 ωc=0.953, 根据可知, 采样周期T取0.18s~0.5s为宜。 为了便于研究采样系统中统一取采样周期T=0.5s。
根据双线性变换法的理论, 将PID离散化后的控制器为:
被控对象为:
采样控制系统的单位阶跃响应如图3 所示:
可见采用双线性变换后, 系统在很长时间内输出量很小, 系统输出经一段时间后上下扰动加剧, 并随着时间的延长, 这种扰动的幅值不断扩大, 严重破坏系统的性能, 最后使系统崩溃。
3 PID控制器采用双线性变换法离散存在的问题分析
PID控制器实际上是由比例环节P, 积分环节I, 微分环节D叠加构成。 那么究竟是其中哪一个环节使得输出产生如此剧烈的振荡, 显然, 比例环节不可能使系统随时间发散, 它只会改变系统增益。 因此, 分别只采用单个积分或微分环节进行离散化后观察系统性能, 以此找出产生振荡的原因。
当只采用积分环节时, 。系统的单位阶跃响应曲线如图4所示:
当只采用微分环节时, 。系统的单位阶跃响应曲线如5所示:
从图4 及5 可以看出, 当只采用积分环节时系统的稳态误差为0, 超调量 σ%=31, 调节时间ts=11.6s, 上升时间tr=1.3s。 积分环节只改变了系统的动态性能, 使系统动态性能变差, 但没有改变系统的稳态性能。 而微分环节的控制器却使整个闭环系统的稳态性能, 动态性能变得极差。 系统的输出经过一段时间后扰动幅值不断加剧, 最终会使整个系统完全崩溃。
由双线性变换离散公式:
可知, S平面的左半平面[Re (s) <0]映射到Z平面, 其关系如下:
因为T>0, 上面的不等式可以简化为:
令z=σ+jω, 则上式为:
即 σ2+ω2<1。 这相应于Z平面单位圆内部。 因此, 双线性变换将S平面上整个左半平面映射到Z平面上以原点为圆心的单位圆内部 (这是Z平面上的稳定区) 。 同理可知双线性变换将S右半平面映射到单位圆外, 将S平面上的虚轴映射到单位圆上 (映射关系如图6 所示) [2,3]。 由于微分环节经双线性变换离散后, 有一个z=-1 的极点和一个z=1零点。 z=-1 的极点对应S平面上虚轴无穷远处, z=1 的零点对应S平面上原点。 z=-1 的极点使系统的单位阶跃响应为等幅振荡曲线, 同时z=1 零点使系统的单位阶跃响应的超调量增大。 所以, 只采用微分环节, 系统的单位阶跃响应在很长时间内输出量很小, 但系统输出经一段时间后上下扰动加剧, 并随着时间的延长, 这种扰动的幅值不断扩大。
4 结论
因此, 从上面仿真及原理性分析可知, 采用双线性变换法离散设计的数字PID控制器中真正使系统不稳定的根源在于微分环节。所以, 可以得出一般性的结论, 即含微分环节的控制器不能采用双线性法离散。
参考文献
[1]高金源.计算机控制系统——原理、设计与实现[M].北京:北京航天航空大学出版社, 2001.
[2]席爱民.计算机控制系统[M].北京:高等教育出版社, 2004.
离散单元法 篇8
边坡是指由于天然地质或工程地质的作用形成的具有一定倾斜度的地质体,按成因共分为两种:自然边坡和人工边坡[1]。例如,意大利1963年发生的瓦依昂水库库岸滑坡,造成了严重的损失,水库也因此失效[2]。由此可见,对边坡失稳进行有效的计算刻不容缓。
目前,边坡稳定分析的方法较多,主要有定性分析法、定量分析法、非确定分析法。本文将利用有限元和离散元对同一边坡的稳定性进行计算分析,并基于计算结果对两种计算方法进行对比。
2 有限元法在边坡算例中的应用
2.1 有限元强度折减法公式
其中,C,φ均为岩土体参数;F为折减系数;C″,φ″均为计算新参数。
2.2 边坡算例
选取国内某矿,边坡尺寸如图1所示。该边坡围岩材料属性见表1。
由于边坡为线性工程,假设沿边坡方向为无限延伸,因此计算模型取二维平面为断面就可满足要求。计算模型见图2。
当折减系数F=1.4时,在坡脚开始出现塑性区(见图3,图4);随着F继续增大,X方向最大位移继续增大,塑性区继续扩展(见图5,图6);当F=3.0时,计算不能收敛,水平位移达到最大值1 016 mm(见图7),塑性区经坡脚贯通至坡顶(见图8),而这也是边坡即将破坏的重要标志,所以有限元计算得到的此边坡的安全系数为2.8。
3 离散元法介绍及边坡算例
离散单元法的突出优势是能够方便的处理非连续介质力学问题。在岩土计算力学方面,由于离散元法能更真实的表达节理岩体的几何特点,便于处理非线性变形和破坏都集中在节理面上的岩体破坏问题,因此被广泛应用于模拟边坡、滑坡等力学行为。离散元法还可以在颗粒体模型基础上通过随机生成方法建立具有复杂几何结构的模型,并通过单元间多种连接方式来体现土壤等多相介质间的不同物理关系,从而更有效的模拟分离等非连续现象,成为分析和处理岩土工程问题的有效方法[3]。
本文对图1中的同一个边坡模型采用M-C模型,经过8个阶段的迭代,得到安全系数为2.52。边坡破坏时的位移矢量图见图9。
4 有限元与离散元计算结果比较
从上文的有限元和离散元计算结果可以看出,有限元计算得到的安全系数2.8略大于离散元计算得到的安全系数2.52。这是由于采用了不同的材料模型以及有限元和离散元计算原理的不同所致。
4.1 计算原理的比较
有限元法是将连续的目标体划分成大量微小单元,但依然保持原有性质。
离散单元法是以牛顿第二定理为理论依据。将目标体看作为刚体,并按照整个目标体的节理裂隙互相镶嵌排列,在空间每个岩块有自己的位置并处于平衡状态。
4.2 在计算有若干节理的岩质边坡方面的比较
运用有限元塑性极限分析方法研究节理岩体边坡的稳定性,含成组节理岩体的力学行为受控于节理面的方位和强度参数,从宏观上将节理岩体视为均质各向异性体,在有限元塑性极限分析中引入节理岩体各向异性的屈服条件,建立了节理岩体边坡有限元塑性极限分析的上、下限法数学模型。而对于节理组比较多,节理比较密集或者比较复杂的时候有限元方法则无能为力,并且其计算结果很难收敛。而用离散元则可以根据节理的参数,计算得到安全系数2.47,小于上文用离散元计算得到的没有节理的完整岩质边坡的安全系数2.52。得到的位移矢量图如图10所示,将其与图9比较可以发现,其破坏形式无节理的情况相似,说明此边坡因为节理的强度比较高,并没有在节理部位发生相对滑移等失稳破坏,但是由于节理的存在,其安全系数有所降低。
5 结语
1)有限元法在处理连续性问题时明显优于离散元;2)离散元在处理有节理的边坡等非连续性问题时,比有限元具有更大的优越性。尤其是对节理组数多并且复杂的情况,离散元法计算过程相对简单,计算结果相对更加精确;3)通过有节理边坡和无节理边坡离散元计算结果的比较可以发现,该节理的存在对边坡的破坏形式产生了一定的影响,但是由于节理刚度比较高,没有产生意外的失稳破坏;4)在用强度折减法计算边坡的安全系数时,计算至不收敛所得到的F值并非真正的安全系数,变形和位移也并非完全符合实际情况,安全系数要在一定基础上予以折减,减小的程度将视情况而定;5)在分析岩质边坡时,由于经常要考虑到实际中的地形地貌和节理的影响,以及考虑到材料的非连续性,离散元在这方面的优势比有限元法要大得多。
参考文献
[1]刘佑荣,唐辉明.岩体力学[M].武汉:中国地质大学出版社,1999:173-179.
[2]张咸恭,王思敬,张卓元,等.中国工程地质学[M].北京:科学出版社,2000:186-215.
[3]卢廷浩.岩土数值分析[M].北京:中国水利水电出版社,2008:60-78.
离散单元法 篇9
铁路电力线路为电气化铁道的行车信号设备供电。它的电气特点[1,2,3,4,5]:所带负荷为一级, 负荷分为车站负荷 (集中且容量大) 和区间负荷 (分散且容量小) ;供电臂40~60 km;是架空线和电缆的混合连接, 线路参数不均匀、不对称。
电气化铁道馈线发生故障后产生的向线路两端传播的暂态行波中包含着故障距离信息, 采集馈线上发生故障后的暂态行波, 通过信号处理手段对行波信号分析, 能够实现馈线故障的精确定位。在现有的三种主要线路故障测距方法 (阻抗法、故障分析法、行波发) 中, 行波法因测距精度高、受系统运行方式影响小等优点而倍受欢迎。行波测距法中年来人们对此进行大量研究, 提出了小波变换、导数法、波形匹配法、相关法等, 其中以小波变换法最为普遍。文献[7]提出了基于小波变换法的暂态行波信号的奇异性检测理论, 利用故障电流行波线模分量的模极大值性质的测距方案。文献[8]采用小波变换法对行波信号进行分析, 根据现代双端行波故障测距原理, 提出了一种与波速无关的测距新方案。文献[9]提出了基于小波包提取算法和相关分析的双端行波测距方法。在使用小波变换法时, 小波基和尺度的选择是一个需要认真研究的问题。
在上述识别行波波头的方法中, 导数法对噪声比较敏感, 但该方法简单快捷, 容易实现, 不失为一种比较实用的方法。导数法在离散信号处理中反应信号变化的变化率会受到信号自身幅度变化的影响[10], 对此, 本文在离散求导的应用基础上, 结合求模极大值法, 并增加了判定的约束函数, 可以有效地消除信号自身幅度值对信号变化率的影响, 进一步提高了故障测距的准确度。
1 算法理论分析
1.1 离散求导的突变点检测
输电线路发生故障后, 将会产生沿线路向两端传输的电流和电压行波, 在行波的波头起始处会表现出一突变值, 这种突变反应在数学上就是一很高的斜率, 因而可以用导数来表示。
将导数原理应用于行波波头检测, 其基本思路是通过在检测点测到行波的一阶导数是否满足设定的阈值及相邻模极大值点之间的约束条件。在离散信号处理中, 导数是以差分形式来表示的, 以一阶导数为例, 有
式中:K1为预先设定的阈值;Y (i) 、Y (i-1) 为采样值;Ts为采样周期。
K1的取值受采样频率的影响, 采样频率较大的情况下, K1的取值也偏大。本文使用的采样频率为100 MHz, Ts=10 ns。原始信号数据Y为一维数组, 求一介导数后的结果得到一维数组D, 取D数组的元素, 构造两个一维数组D1和D2, 对其进行运算处理求其极大值点序号。运算表达式如式 (2) 。
其中,
再对式 (2) 的结果求得极大值点序号后, 对极大值点处导数值进行取绝对值运算, 得到一个行向量W。
1.2 阈值约束查找波头位置点
取得电压行波信号的离散求导的模极大值点之后, 得到一个电压信号突变点数据W;为了提高导数法的可靠性, 同时也使得其具有一定的抗干扰能力, 对W向量进行循环运算, 并增加约束判定函数来查找突变的模极大值点。判定函数表达式如式 (3) 。
其中:W (j+1) , W (j-1) , W (j) 表示相邻的模极大值点斜率;K2为限阈值;L为W向量的长度;N为初始行波信号之前的噪声长度。
采集到的数据在检测到行波的起始波之前会有一段噪声, 因此设定阈值时要减去N长度个点。同时N的选择会影响离散求导模极大值法识别波头点的位置, 从而影响测距精度。如果N值取得过大, 则会漏掉一些模极大值的点, 反之取值过小的话, 会误将噪声信号中的极大值点判定为波头点的位置, 根据对采集到的故障数据分析经验, 一般取N的长度在200个点左右。
2 现场行波数据的采集
为验证该算法的有效性及其在实际故障下的波头识别精度, 在天津南仓变电站将实验设备现场挂网采集数据。限于篇幅, 这里仅简单介绍下实验设备及系统的动作过程。故障发生后, 馈线动作保护, 断路器跳闸, 此时, 负荷开关合闸, 对实验设备中的高压电容充电到设定电压后分段, 控制电路触发使得高压脉冲开关合闸, 对高压电容器充电后, 高压电容器再对故障线路进行放电, 送出高压脉冲, 馈线下方的传感器采集故障电压波形, 传输给高速采集卡, 录入存储故障波形数据, 使用本文的方法识别波头。操作步骤流程如图1所示。
图2所示为选取的两个实测217#馈线上倒闸时, 行波电压数据的波形。
3 本文算法在波头识别中的应用
由现场实验所得行波数据, 将离散求导模极大值的方法应用于图2所示的倒闸时馈线行波数据, 进行行波波头的识别, 标记出行波波头的位置点, 将得出的结果用于测距, 然后与线路实际距离进行比较, 验证其准确性。应用本文的算法进行处理后, 所得数据波形结果见图3所示。
图3中圆圈标记的位置点即波头位置, 从图中可以看出在标记的起始波的波头 (第一个圆圈处) 后面, 会有多次波的反射。表1中给出了本自动算法处理行波数据后所得结果。
由表1中的数据可知第一个反射波波头和起始波的波头之间的时间间隔Δt, 取加拿大的B.C.Hydro行波定位系统中, 线模分量的行波波速v=2.95×108为本线路上的行波传播速度。
由式 (4) 可计算出线路距离。已知天津南昌站上海方向距离为9.26 km。上文所述选取的两个现场实测数据的分析结果记录到表2中。
按照类似的方法, 对石家庄南新城变电站216号馈线停送电时检测到的行波波形进行分析, 定位其初始行波与线路末端反射行波起始点, 并根据天津南仓变电站数据分析时使用的波速计算故障距离, 结果记录到表3中。
已知石家庄216线馈线长度为22.19 km, 由表3所示, 可得行波测距系统定位误差最大为 (22278-22190) m=88 m。
由表2和表3可知自动测试算法所得测距结果误差在100 m以内, 因此对于行波波头的识别有效, 并能够较准确的实现测距。
4 结论
铁路变配电所实测数据结果表明, 将基于离散求导阈值约束法应用于复杂干扰情况下的行波故障测距中:1) 方法原理简单, 数据计算量小, 具有灵活性。2) 能够有效查找出行波信号经故障点反射回来的幅值较弱的此类波头位置点。
更进一步深入研究, 提高精确度, 并结合其他行波测距算法, 及现场实验所设计的硬件设备, 做成一套完整的故障测距装置, 实现超高速继电保护和精确故障测距, 能够较好地解决电气化铁道馈线保护和快速、精确定位问题。为行波保护提供了很好的理论和实践依据。
参考文献
[1]葛耀中.新型继电保护与故障测距原理与技术[M].西安:西安交通大学出版社, 1996:4-8, 219-220, 266.GE Yao-zhong.New protective relay and principle and technology of fault location[M].Xi'an:Xi'an Jiaotong University Press, 1996:4-8, 219-220, 266.
[2]李群湛, 连级三, 高仕斌.高速铁路电气化工程[M].成都:西南交通大学出版社, 2005.LI Qun-zhan, LIAN Ji-san, GAO Shi-bin.High speed railway electrification project[M].Chengdu:Southwest Jiaotong University Press, 2005.
[3]薛永端, 徐丙垠, 陈平, 等.铁路自闭/贯通线路故障行波技术[J].电力设备, 2004, 5 (5) :41-44.XUE Yong-duan, XU Bing-yin, CHEN Ping, et al.Railway self closing/through transmission line faulttravelling wave technology[J].Electric Power Equipment, 2004, 5 (5) :41-44.
[4]蔡玉梅, 何正友, 王志兵, 等.行波法在10kV铁路自闭/贯通线故障测距中的应用[J].电网技术, 2005, 29 (1) :15-19.CAI Yu-mei, HE Zheng-you, WANG Zhi-bing, et al.Travelling wave method used in railway self closing/through10kV transmission line fault location[J].Power System Technology, 2005, 29 (1) :15-19.
[5]刘明光, 程辉海, 潘光瑞, 等.贯通线故障自动处理系统研究[J].电气化铁道, 2005 (1) :8-10.LIU Ming-guang, CHENG Hui-hai, PAN Guang-rui, et al.Automatic processing system research of power supply through line fault[J].Electrified Railway, 2005 (1) :8-10.
[6]李强, 王银乐.高压输电线路的故障测距方法[J].电力系统保护与控制, 2009, 37 (23) :193-195.LI Qiang, WANG Yin-le.Fault location methods for high voltage power transmission lines[J].Power System Protection and Control, 2009, 37 (23) :193-195.
[7]马丹丹, 王晓茹.基于小波模极大值的单端行波故障测距[J].电力系统保护与控制, 2009, 37 (3) :56-58.MA Dan-dan, WANG Xiao-ru.Single terminal methods of traveling wave fault location based on wavelet modulus maxima[J].Power System Protection and Control, 2009, 37 (3) :56-58.
[8]尹晓光, 宋琳琳, 尤志.与波速无关的输电线路双端行波故障测距研究[J].电力系统保护与控制, 2011, 39 (1) :36-38.YIN Xiao-guang, SONG Lin-lin, YOU Zhi.Study of locating for transmission line double terminal traveling waves unrelated to wave speed[J].Power System Protection and Control, 2011, 39 (1) :36-38.
[9]周湶, 卢毅, 廖瑞金, 等.基于小波包提取算法和相关分析的电缆双端行波测距[J].电力系统保护与控制, 2012, 40 (1) :1-3.ZHOU Quan, LU Yi, LIAO Rui-jin, et al.Double terminal traveling wave fault location for cable based on the wavelet packet extraction algorithm and correlation analysis[J].Power System Protection and Control, 2012, 40 (1) :1-3.