ldda动接触力的迭代算法

ldda动接触力的迭代算法

在许多实际工程问题中，有必要模拟接触面的开口运动。例如:将拱坝和地基作为整体分析求解坝肩抗震稳定时,在坝肩附近,存在各种构造面,这些构造面相互切割而形成可能滑动块体,可能滑动块体的边界面必为接触缝面。另外,为了减少温度应力的影响,拱坝的浇筑均是分块进行的为了较为合理地考虑坝、地基(包括可能滑动块体不连续面的影响)的相互作用,接触缝隙的数值模拟是极为重要的。接触缝隙在地震作用下的开合属于接触非线性问题,这是区别于材料非线性和几何非线性的第三类非线性问题,即边界条件非线性。对这一问题的研究始于100年前Hertz对二次曲面挤压问题的研究。起初,采用的是解析方法,近30年来,随着计算机技术的发展,接触问题的数值解法受到了众多研究者的关注,发展了各种求解接触问题的数值计算方法,包括:(1)以Lagrange乘子为未知量的各类Lagrange乘子或修正的Lagrange乘子方法;(2)罚函数方法;(3)线性补偿方法;(4)接触单元方法;(5)动接触力方法等。清华大学的刘晶波教授结合时域显式有限元与人工透射边界,初步应用研究了动接触力模型LagrangeDiscontinousDeformationAnalysis(LDDA)是在接触界面上采用拉格朗日乘子,而接触判断使用DiscontinousDeformationAnalysis(DDA)规则,在块体内部进行有限元细分的DDA方法。众所周知,不连续变形分析(DDA)是由石根华博士提出的,主要是用来求解岩石介质中不连续体的变形和运动。DDA理论体系较完整,常用于计算边坡的滑动,岩石块体的倾覆,隧洞洞室的崩塌,其主要缺点是:在接触缝面应用了压缩弹簧以模拟接触刚度,而弹簧常数不易确定;其次,因使用块体单元,块体元内部的应力是常应力,应力精度不高。LDDA是对DDA的改进,接触缝面上使用了拉格朗日乘子,其物理意义代表接缝法向和切向接触力,另外,块体内部使用常规有限元离散,应力精度较高。接缝张开时,拉格朗日乘子为0,接缝闭合时,真实的接触力即为拉格朗日乘子的值,因此,LDDA方法对接缝的处理更真实、更自然。LDDA与动接触力模型的不同之处是:构成接缝的两侧不要求是节点对,可以对坝体和基础使用不同的有限元网格,因此,坝体单元可据需要进行加密;LDDA可以计算大滑移,LDDA的这种特性可用来重现溃坝的过程。因无节点对的限制,LDDA的接触是利用主从节点来实现的,主从节点在计算过程中不是固定不变的,它们随接缝表面的相对运动而变化,计算时要不断地搜索主从接触点,并对其进行接触判断;LDDA对全部接触力进行平衡迭代,因此,所求出的动接触力精度较高。LDDA方法LDDA分析中至关重要的是接触力的迭代求解,它是LDDA能否成功解决问题的关键。本文将提出一个新的动接触力的迭代求解算法,该算法的松驰因子是依据块体接触的物理特性而非按照数学理论选用的。算例表明,该算法具有收敛速度快,收敛精度高的特点。虽然本文的算例是针对拱坝进行的,但从理论和方法上,本文算法对缝面接触问题的求解是普遍适用的。因此,本文为接触问题的计算提供了一个有效的方法和工具。1质量矩阵与阻尼矩阵的显式积分格式利用虚功原理式中,M、C、K分别是系统的质量矩阵、阻尼矩阵和刚度矩阵;有限元方程式(1a),可利用中心差分法离散,中心差分格式为:式(2)和式(3)中,上标t表示时间,Δt为时间间隔。将式(2)和式(3)代入式(1a)可得质量矩阵与阻尼矩阵均为对角矩阵时的显式积分格式。当阻尼矩阵为非对角矩阵时,李小军等建议了一种非对角阻尼矩阵的显式积分格式将式(4)代入式(1a)可得:从式(5)可见,只要质量矩阵为对角阵,则式(5)是解耦的,已知当前时刻的位移、速度和拉格朗日乘子向量,则可求得下一时步的总体位移向量。t+Δt时刻的速度和加速度可由下列递推公式求出:若令:则式(5)可写为:将式(10)代入式(1b)得:已知当前时步的位移和速度,由式(11)可求当前时步拉格朗日乘子(即接触力)。对接触状态的判断应遵循下列准则:(1)脱离状态:(2)粘着状态:(3)滑动状态:其中,Λ2ldda的接触刚度计算方法在上述LDDA方法求解多变形体的接触问题中,至关重要的一步是获得接触面上的接触力,但接触力与接触状态有关,因此,一般需要迭代求解。已有不少研究者提出了各种求解方法,包括实验误差法(1)假设初始状态为接触状态,由式(11)求接触力,由式(12)对接触力进行修正得Λ(2)由式(10)求位移u(3)求接触间隙:u(4)由下式迭代求解接触力:(5)停止迭代的条件是u上述迭代计算的第(4)步,Λ式(14)中,m这一接触刚度的计算可推广至LDDA,设i为主节点,i′为从节点,则:式(15)中,m3有缝坝施工及求解某拱坝为抛物线双曲拱坝,大坝最低建基面高程925m,坝顶高程1135m,最大坝高210m。顶拱拱冠处厚10m,拱冠梁底厚55m,厚高比0.262。坝顶中心线弦长523.82m,弦高比2.494。水库正常蓄水位1130m;水库运行低水位(死水位)1120m,水体容重1.0t/m坝体横缝布置及编号见图1,横缝的抗拉强度为0,抗剪强度:摩擦系数0.8,粘滞力0。本考题计算:在坝体自重,正常蓄水位水荷载,坝前淤砂压力,温度荷载(温降)作用下,拱坝的应力、位移和横缝的开合。计算时步0.001s,瑞利阻尼系数α=0.908、β=0.000959。本考题计算不模拟基础,在坝与地基交接面上采用固定约束。共有块体单元264个,节点数1108,15条横缝,总自由度1716。本考题虽然没有考虑地震作用,但显然是接触非线性问题,拱坝横缝将在水压力等静力荷载作用下,渐渐闭合,最终达到静力平衡而稳定。现实中的拱坝,在蓄水前已对横缝进行封拱灌浆,要模拟这一施工过程显然太复杂。本考题将作两种简化,与此有两种计算方案,其一是拱坝横缝和坝体自重、水砂压力、温度载荷同时作用,导致横缝渐渐闭合,最终拱坝成为整体的受力结构;其二是先求无缝坝在上述静力荷载作用下的静力反应,然后将这一静力反应作为初始条件作用于有缝拱坝上,这时的横缝也将渐渐闭合至真实状态。上述求解过程中,横缝间的接触力也经历一个不断调整至真实值的过程。关键的是横缝间接触力的求解要应用上述迭代求解算法。在不使用改进的Uzawa迭代法的情况下,作者应用不同的积分方案图2是两方案下拱冠梁顶拱上游面位移时程;图3是两方案下拱冠梁顶拱法向接触力时程。从上述时程图可见,大约10s后,两方案的位移和接触力均趋于稳定,但方案2的收敛速度比方案1要快得多。可见,对于拱坝这类结构,先求整体坝的静力反应,然后将其作为初始条件,求解分缝坝的响应是一个比较好的求解方法。进一步分析可以看出,有缝坝在静力荷载作用下,除1缝、2缝和14缝、15缝顶拱上游侧有0.2mm左右的张开外,其余各处,横缝都是闭合的。虽然本算例只是简单地将坝基作为刚性处理,但也能基本反应有缝拱坝的受力特点。即,在静荷载作用下,拱坝横缝基本上是闭合的。进一步分析还可看出,整体坝和分缝坝在静力荷载作用下,坝体顺河向位移等值线图基本一致。最大顺河向位移发生在顶拱拱冠下游,其值为4.165cm(见图4)。同时,在静力荷载作用下,整体坝和分缝坝坝体最大主应力等值线图大致相同,整体坝和分缝坝的最小主应力等值线图也大致相同。从图5可见,最大主应力最大值为7.18MPa,发生在上游面坝踵附近。从有缝坝与整体坝的位移、应力和接触力这三方面的比较可见,LDDA基本反应了拱坝横缝在大坝蓄水后闭合,拱坝拱圈作为整体承受水荷载作用的这一事实。从而也验证了LDDA和本文提出的关于接触力迭代算法的正确性。上述算例,使用本文方法(方案2)在作者所用的PC机上(HPd530,奔IVCPU2.80GHz,1GB内存)进行计算花时114s,而使用动接触力方法计算时,达到相同计算精度,在相同的微机上需时920s,这表明本文算法有较高的计算效率。其主要原因是动接触力算法采用了静动