连续体结构强非线性仿真——离散实体单元法

第一章 绪论
　　1.1 背景及意义
　　在工程科学研究领域，对于复杂的力学问题和物理问题，采用解析方法对其数学模型进行计算时，由于数学方程的非线性和求解域的复杂性很难得到精确的结果。为了克服解析方法的局限性，长期以来，科研人员创建和发展了一种新的解决方法——数值计算方法。特别是近几十年来，随着电子科学与计算机技术的飞速发展和广泛使用，各种新型数值计算方法相继提出，为解决结构的大变形、强材料非线性、冲击和断裂等复杂力学问题提供了新的解决途径[1]。目前，在土木工程领域，科学问题的研究主要采用的方法包括科学试验、理论分析和数值计算，三种研究方法相辅相成推进土木工程学科的发展。由于结构试验的费用较高并且一些力学性能参数很难测得等限制条件[2]，数值计算方法已经成为理论实践和工程设计中不可缺少的重要手段。作为科学问题的重要研究工具，数值分析在土木工程学科中具有重要的地位。
　　在结构的寿命服役期内，可能会遭受爆炸、火灾、强震和恐怖袭击等突发灾害事件，这将引起结构的局部破坏和整体倒塌问题[3]。这些问题通常伴随着结构的几何非线性、材料非线性、断裂和冲击等复杂力学问题，在极端荷载下，结构将从连续体转变为非连续体，原有的内力平衡被打破，结构在变形过程中寻求新的内力传递路径，直至结构出现新的平衡状态。对于这类强非线性问题的全过程数值仿真研究多年来一直是研究的热点和难点。
　　在计算力学领域，数值计算方法可分为两大类：连续介质力学方法和非连续介质力学方法。连续介质力学方法是将分析的系统简化为数学意义上的连续体并基于变分原理得到问题的解答，其中有限元方法(Finite Element Method，FEM)因理论更加完善、拥有大量成熟的商业软件，近几十年在结构工程、岩土工程、流体力学、电磁学等领域已得到了广泛应用。FEM 的基本特点为：(1) 求解域是连续的，必须满足位移连续条件和变形协调方程；(2) 需要划分网格并且结构的力学响应对网格存在极强的依赖性；(3) 需要求解非线性方程组和进行刚度矩阵运算，并且伴随结果不收敛问题，非线性问题计算效率偏低；(4) *关键的一点是FEM难以处理不连续和断裂等复杂力学问题。FEM 的理论基础是变分原理和连续体力学，由于网格畸变严重、单元消失和位移场不连续等困难，如图1.1(a)～(c) 所示，采用该方法较难计算结构的大变形、断裂和强材料非线性等力学行为，需要引入其他方法或技术对其进行改进和修正。如Wang 等[4] 重新定义了离散的弱梯度算子，基于此算子建立了二维泊松方程，用于解决二维几何大变形问题。Kopacz 等[5] 根据耗散定理修正了流体方程，用于计算流体中粒子之间的接触碰撞问题。Cloud 等[6] 研究了FEM 中非线性浅壳的平衡问题，Desai 等[7] 建立了薄层单元用于模拟非连续截面的断裂问题。当结构涉及几何非线性和材料非线性问题时，计算过程中需要求解非线性方程组与迭代计算，不仅耗时，而且很难控制数值计算的收敛性与稳定性。另外，如果此时结构再经历断裂等非连续问题，基于连续体力学的数值方法仍然要求离散结构体系满足连续条件，将造成连续介质方法的计算更加复杂，求解更加困难。
　　图1.1 FEM 中存在的主要问题
　　1.2 数值计算方法现状
　　在土木工程中，结构在荷载作用下产生大变形从而进入非线性工作阶段，直至发生破坏的过程中，通常涉及结构的几何非线性、材料非线性和断裂等复杂行为。在计算分析中进行结构的大变形、强非线性和断裂问题的全过程仿真，从而实现结构从连续体向非连续体的转换，这对现有的数值计算方法是一项巨大的挑战。为了克服 FEM 的局限性，各国学者对此做出了很多研究工作，开始尝试采用非连续介质计算方法解决连续体向非连续体发展的问题，提出了各种各样的数值模型和数值方法，主要包括有限质点法、无网格法、非连续变形分析法、数值流形法和离散单元法等。这些数值方法彼此取长补短，相互借鉴又相互融合，为数值计算领域注入新的血液，推动了数值方法的发展，使得我们掌握的分析方法不断拓展与深入。下面将阐述近年来数值计算方法的发展与现状。
　　1.2.1 扩展有限元法
　　扩展有限元法 (Extended Finite Element Method，XFEM) 是目前处理材料断裂问题的杰出代表之一。XFEM 由美国学者 Belytschko[8] 和 Black[9] 于1999年提出，首次应用XFEM 模拟了线弹性裂纹的扩展。基于传统FEM，当采用XFEM 模拟裂纹扩展时，计算对象的连续区域仍然应用传统FEM，在包括不连续边界的较狭窄的计算区域内，修正传统有限元的位移近似函数，增加了对不连续边界的描述函数。在满足单位分解的前提下，在位移近似函数中增加反应裂纹间断特性的函数项，称为富集函数 (Enrichment Funciton)[10]。同时，采用水平集方法(Level Set Method，LSM) 或快速推进法 (Fast Marching Method，FMM)建立单元间断界面[11]，使裂纹间断特性的描述独立于有限元网格。与传统FEM相比，*根本的区别是XFEM 克服了裂纹尖端高应力和变形集中区进行高密度网格划分带来的困难，以及当裂纹扩展时也不需要对网格重新划分[12]。
　　自XFEM 提出以来，后经不同学者的发展和改进，目前XFEM 已经广泛应用于各种断裂问题的研究中。Daux 等[13] 建立了反映分叉裂纹单元间断性质的联结函数，对分叉裂纹和多裂纹交叉进行了XFEM 模拟。Rethore 等[14] 提出了一种基于拉格朗日守恒格式的用于计算二维裂纹扩展的动态应力强度因子的技术。陈亚宾[15] 采用XFEM 对素混凝土中裂纹开裂进行了系统地研究。茹忠亮等[16] 建立了预留裂纹的钢筋混凝土梁的三维扩展有限元模型，应用软化模型对钢筋混凝土梁的复合断裂过程进行了模拟分析。Song 等[17] 比较了XFEM、元素删除法(Element Deletion Method) 和互元裂纹法 (Interelement Crack Method)处理裂纹扩展与裂纹分支的优势与不同。Fires 等[18] 对XFEM 的单位分解函数、 近似函数和水平集方法等基本思想和公式进行了描述，并介绍了XFEM 的程序实现和应用展望。
　　1.2.2 凝聚有限元法
　　另外一种连续体计算领域模拟裂纹扩展的数值模型为凝聚裂纹模型(Cohesive Zones Model，CFM)，也称为凝聚有限元法(Cohesive Finite Element Method，CFEM)。凝聚模型的概念首先由苏联学者 Barenblatt[19] 于1962 年提出，用于解决裂纹尖端的应力在理论上无穷大的问题。该模型假设在材料裂纹尖端存在凝聚力损伤区，当损伤区内应力满足凝聚力开裂准则时，单元沿着其边界发生分离，从而模拟裂纹的扩展[20]。后经过Hillerborg[21]、Xu[22]、Needleman[23]、Turon[24]、Ferté[25] 和张志春 [26] 等学者对CFM 的发展和修正，将CFM 与FEM 相结合，采用CFEM 对二维和三维的线弹性和弹塑性裂缝的扩展与裂缝的分叉进行了模拟分析。近年来，CFEM 用于微观尺度和多尺度断裂模型的研究中，He[27] 和Guo[28]分别基于CFEM 建立了确定镁合金断裂韧性的微观方法，以及开发了嵌入式原子超弹性本构模型，通过原子信息确定中尺度和宏观尺度材料的力学行为。对于裂纹扩展问题的研究，XFEM 和CFEM 都取得了较好的研究成果，但是无论是XFEM 还是CFEM 都摆脱不了计算网格的限制。当材料中出现多条裂纹时，由于裂纹复杂的模拟机制和高密度网格等因素，这两类改进的FEM 很难进行多条裂纹的扩展模拟。为了克服FEM 这类网格类数值模拟方法解决裂纹扩展问题的局限性，近年来许多学者将研究重点转移到质点类、粒子类和无网格方法上，并做出了许多重要的研究工作。
　　1.2.3 有限质点法
　　美国普渡大学Ting 教授[29] 基于向量力学和数值计算结合的概念提出了直接用离散的 “点值” 和质点运动定律来描述结构行为的构想，在结构物理模型的基 础上，引用广义向量式力学作为运动和变形的准则，发展了向量式结构与固体力学[30]。2006年，罗尧治教授[31] 基于向量式结构与固体力学的基本概念，提出了面向结构工程的数值分析方法——有限质点法(Finite Particle Method，FPM)[32]，针对结构的非线性问题和不连续行为进行分析。该方法在空间上，将结构描述为一群质点的集合，质点间采用单元连接，用牛顿第二定律取代连续体的偏微分方程来描述质点的运动[33]。在时间上，将运动历程划分为一系列途径单元[34]。在内力计算上，采用虚拟的逆向运动获得单元的纯变形，从而采用材料力学公式计算内力。在运动公式求解上，采用中央差分的显示积分法[35]。
　　喻莹等[36] 推导了FPM 计算杆系结构的动力响应、几何非线性和材料非线性问题的方程，进行了双层网壳结构在周期荷载作用下破坏全过程的仿真研究[37]。杨超等[38] 研究了FPM 计算平面固体几何大变形问题。张鹏飞等[39] 在FPM 中建立了四面体实体单元，用于三维固体弹塑性问题的分析。屈曲、褶皱等失效行为是结构受力过程中经历了几何非线性和材料非线性的一类特殊现象，杨超[40]和罗尧治[41] 根据张力场理论，采用FPM 对薄壳结构的褶皱形态进行了研究，分析了各质点的受力状态，重点关注了褶皱区域的处理。王震等[42] 推导了三角形薄壳单元FPM 的内力方程，对薄壳结构的屈曲、后屈曲、破碎和碰撞等问题在FPM中的处理方式进行了讨论。
　　FPM 采用多质点的力学模型用于数值计算，近年来主要用于解决空间结构的复杂力学问题，采用统一的计算框架对结构的力学行为进行分析[43]，计算中不区分线性和非线性问题，对于强非线性问题不需要转换计算模块[44]。该方法不存在微分方程的假设，克服了函数描述的困难[45]。与FEM相比，FPM在处理大变形、断裂、碰撞等问题方面具有较大的优势[46]。
　　1.2.4 无网格法
　　无网格法(Meshless Method) 的研究起始于20 世纪70 年代，其核心内容为形函数的建立理论，这也是无网格法与FEM 的根本区别之一。无网格法采用离散的点分解求解区域，通过离散的点构造近似函数，可以彻底或部分消除FEM计算网格的限制[47]，因此称为无网格法。该方法不需要划分分析对象的计算网格，在保证计算精度的前提下，减小了复杂力学行为的计算难度。但是，由于无网格法中近似函数普遍比较复杂并且一般不具有插值属性，因此无网格法的计算量较大并且本质边界条件的施加比较困难[48]。发展至今，无网格法相继已经提出了十多种，但是各种无网格法的本质区别为采用何种加权余量法和试探函数进行微分方程求解[49]，比如光滑质点流体动力学法(Smooth Particle Hydrodynamics，SPH)、重构核粒子法(Reproducing Knernel Particle Mehtod，PKPM)、无单元伽辽金法(Element Free Galerkin Method，EFGM) 和小波粒子法 (Wavelet Particle Method，WPM)[50]等。1977年Lucy[51] 和 Gingold[52] 基于空间场函数和核函数的概念，首次提出了SPH，并采用该方法研究了天体星系爆炸问题。为了解决SPH 的计算精度和计算稳定性方面的问题，Johnson[53] 提出了标准平滑算法，Vignjevic[54]、Swegle[55]、Dyka[56] 和 Chen[57] 提出了 SPH 不稳定的因素和相应的稳定性优化方案。目前，该方法已经成功应用于高速冲击、爆炸、裂纹扩展等问题研究中。Belyts

《连续体结构强非线性仿真——离散实体单元法》系统地总结了作者近年来关于离散实体单元法的研究成果。《连续体结构强非线性仿真——离散实体单元法》共六章，第一章介绍了数值计算方法发展历程；第二章介绍了三维离散实体单元法的物理模型和基本原理、运动方程的推导、接触本构方程的建立、阻尼和计算时步的确定；第三章介绍了离散实体单元法的应变能密度的计算、弹簧刚度系数与材料弹性常数关系式的确定、几何大变形分析；第四章介绍了离散实体单元法理想弹塑性和双线性等向强化弹塑性本构模型的建立、流动准则、加卸载准则和接触力增量计算流程的建立、弹塑性大变形分析；第五章介绍了离散实体单元法边界球元的分类、各类边界球元弹簧刚度系数的确定及边界效应分析；第六章介绍了离散实体单元法裂纹扩展准则、断裂软化模型的建立及动态断裂分析。

目录
前言
第一章　绪论　1
1.1　背景及意义　1
1.2　数值计算方法现状　2
1.2.1　扩展有限元法　2
1.2.2　凝聚有限元法　3
1.2.3　有限质点法　4
1.2.4　无网格法　4
1.2.5　非连续变形分析法　5
1.2.6　数值流形法　6
1.2.7　离散单元法　7
1.3　离散实体单元法的优势　9
1.4　本书主要内容　9
第二章　三维离散实体单元法的基本原理与推导　12
2.1　离散单元法　12
2.1.1　基本思想与假设　12
2.1.2　基本特点　12
2.1.3　动态松弛法　13
2.1.4　计算循环　14
2.1.5　离散单元法的不足　14
2.2　三维离散实体单元法　15
2.2.1　计算特点与物理模型　15
2.2.2　运动方程的建立与推导　18
2.2.3　接触本构方程的建立与推导　19
2.2.4　弹簧的本构关系　24
2.2.5　阻尼的确定　24
2.2.6　计算时步的确定　27
2.3　离散实体单元法计算软件的开发　29
2.3.1　编程语言和编程平台　29
2.3.2　程序数据结构　30
2.3.3　计算流程　31
2.3.4　接触判断算法　32
2.4　离散实体单元法与有限单元法区别与联系　35
2.5　本章小结　37
第三章　离散实体单元法的几何大变形分析　38
3.1　引言　38
3.2　离散实体单元法分析思路　39
3.3　弹簧刚度系数与材料弹性常数关系式的确定　41
3.3.1　球元间的弹簧类型　41
3.3.2　局部坐标下球元间相对位移的计算　42
3.3.3　离散实体单元法的应变能密度计算　43
3.3.4　弹簧系统表示的应变能密度计算　46
3.3.5　基于能量守恒原理弹簧刚度系数的确定　49
3.4　算例分析与验证　51
3.4.1　三维弹性块体结构的均布荷载受压分析　51
3.4.2　悬臂梁受端弯矩的大变形分析　57
3.4.3　刚架受对边集中力作用的大变形分析　60
3.4.4　Williams双杆体系　64
3.4.5　平面门式刚架大变形　67
3.4.6　空间刚架　71
3.5　本章小结　74
第四章　离散实体单元法的弹塑性分析　75
4.1　引言　75
4.2　离散实体单元法分析思路　77
4.3　屈服方程的建立　78
4.4　离散实体单元法理想弹塑性本构模型　80
4.4.1　流动准则　80
4.4.2　弹塑性接触本构方程　82
4.4.3　接触力增量计算流程　83
4.4.4　加卸载准则　84
4.5　离散实体单元法双线性等向强化弹塑性本构模型　85
4.5.1　弹塑性接触本构方程　85
4.5.2　接触力增量计算流程与加卸载准则　87
4.6　弹塑性大变形算例分析与验证　89
4.6.1　方钢管受轴向动力荷载的弹塑性大变形分析　89
4.6.2　薄板在轴压荷载下的动力屈曲计算　93
4.6.3　单调加载下冷弯薄壁型钢剪力墙的抗剪性能分析　98
4.6.4　Williams双杆体系　103
4.6.5 六角星形穹顶结构　104
4.6.6　空间直角刚架　107
4.7　本章小结　109
第五章　离散实体单元法的边界效应分析　111
5.1　引言　111
5.2　离散实体单元法计算模型的边界球元　114
5.3　边界面球元弹簧刚度与弹性常数关系式的推导　115
5.3.1　边界面球元受力状态分析　115
5.3.2　边界面球元的应变能密度计算　115
5.3.3　基于能量守恒原理的边界面球元弹簧刚度系数的确定　118
5.4　边界棱球元弹簧刚度与弹性常数关系式的推导　120
5.4.1　边界棱球元的受力状态分析　120
5.4.2　边界棱球元的应变能密度与弹簧刚度系数的推导　121
5.5　边界角球元弹簧刚度与弹性常数关系式的推导　123
5.5.1　边界角球元的受力状态分析　123
5.5.2　边界角球元的应变能密度与弹簧刚度系数的推导　123
5.6　算例分析与验证　125
5.6.1　薄板的弹塑性弯曲分析　125
5.6.2　开裂薄板的屈曲分析　132
5.7　本章小结　137
第六章　离散实体单元法的断裂分析　138
6.1　引言　138
6.2　离散单元法分析思路.139
6.3　裂纹的类型　141
6.4　基于连续介质力学的裂纹扩展准则　142
6.4.1　*大应力准则　142
6.4.2　应变能密度因子准则　143
6.4.3　应变能释放率准则　145
6.5　离散单元法的裂纹扩展准则　146
6.6　离散实体单元法的断裂模型　147
6.6.1　双线性软化模型　147
6.6.2　三线性软化模型　150
6.7　构件断裂分析与验证.152
6.7.1　双悬臂梁裂缝开展试验模拟　152
6.7.2　不对称集中荷载作用下圆管的裂纹动态扩展　157
6.7.3　矩形梁受拉扭作用断裂分析　160
6.8　考虑杆件断裂的单层网格结构倒塌数值模拟研究　163
6.8.1　模型建立及程序验证　164
6.8.2　单层网格结构倒塌性能影响因素分析　170
6.9　本章小结　178
参考文献　179