有限元学习必须知道：隐式与显式有限元算法

2017-02-27 by:CAE仿真在线来源:互联网

导语:对于有限元的学习,最最关键的其实不是本构方程,也不是屈服准则,而在于对求解算法的理解。本人根据长期学习经验,在美国作访问学者学习经历,纯属兴趣,在这里做一点肤浅的总结,完全原创。

隐式与显式有限元最大的区别在于是否迭代,是否所有的物理量在同一时刻获得。采用隐式迭代求解平衡方程(位移、速度和加速度)、而不管是否用隐式与显式的方法(前向或者后向欧拉求解方法)求解本构方程(应力和应变)叫做隐式有限元;用显式时间积分的方法求解叫做显式有限元。

首先,对于本构方程的求解,通常分为前向和后向欧拉算法。对于后向欧拉算法求解弹塑性问题,所有的物理量(包括等效塑性应变增量、N+1迭代步的应变和应力以及相关依赖于solution的状态变量)均是同时求解获得,因为涉及到多个物理量,而通常情况下他们是相互依赖、相互成为函数,所以必须通过牛顿迭代同时求解几个方程(如采用试应力方程、屈服函数径向返回算法(对于各向异性,也叫回映算法,最近点的投射算法)联合求解等效塑性应变增量)。对于前向欧拉,直接由N时刻的应力和应变求出N+1时刻的应力和应变,无需迭代。

其次,对于平衡方程的求解,通常分为隐式和显式有限元算法。对于隐式有限元算法,由应力平衡方程+边界条件变分之后获得的刚度方程KU=F,隐式求解必须引入雅可比矩阵(二次收敛、只影响计算速率、不影响数值精度;K又称为雅可比),其是实时更新的,是N+1时刻的应力、应变以及状态变量(如损伤内变量)的函数,隐式求解是很robust的,确保了计算精度,但是不足之处在于计算非常expensive,每次迭代都要计算K的逆矩阵,也容易产生数值收敛性问题,目前解决的方法有弧长法、粘性阻尼法等,个人认为粘性阻尼法效果最好。

对于显示算法,采用时间积分,用t+1时刻的积分点应力、应变,获得t+1时刻的节点位移,无需迭代求解,也不需要雅可比矩阵(应力对应变偏导数);如果硬是要有,连续雅可比,基于本构模型而不是刚度方程推导近似的连续雅可比。对于显示算法,单元高斯积分点应力、应变的求解可用前向或者后向欧拉方法,然后通过时间积分求取节点位移。本质上,平衡方程中位移的迭代求解与本构方程中的应力、应变求解没有关联,这点很容易造成误解,很多时候将前、后欧拉算法视为显式和隐式的区别,大大错误。通常应用较广的显示算法纽马克法、威尔逊-sita法,其中改变纽马克法中的两个参数,可以实现隐式与显式求解,其中alpha=0.5和beta=0是中心差分法(二阶精度)。目前一个大的误区认为只有显示算法可以求解动力学问题,隐式只能求解准静态问题(如低速冲击),alpha=0.5和beta=0.25就是隐式,所有的物理量在t+1时刻同时求解,通常ABAQUS软件中所说的隐式动力学求解采用了斯坦福大学Hilber、HUGHES院士(现在德克萨斯大学奥斯丁分校)和加州大学伯克利分校Taylor院士提出的无条件稳定隐式差分算法,可以求解低速动力学问题,缺点是不适合含阻尼的求解、计算效率不高;alpha=0.5和beta=0时的纽马克法更适合求解动力学问题,主要原因在于比隐式求解计算效率更高,不足之处在于其是条件稳定,时间增量过大位移解容易震荡,根本原因是差分算法的条件稳定导致的,时间增量必须非常小(其值越大,一方面不稳定、另一方面计算误差也更大),其依赖于波速、弹性模量和最小单元网格尺寸,这是显式算法计算最耗时的地方。相对于隐式算法,显式算法的功能更强大,通常计算依赖于率的变形和应力,也可以求解稳态问题,如alpha=0.5和beta=0时,对于刚度方程中引入阻尼矩阵后,叫做动态松弛法,可以解决静力学问题的一些稳态问题(如重力、预应力引起的初始应力)。此外,一些准静态的剪切自锁问题,本质上有解,但是用牛顿法求解失效,中心差分引入质量矩阵后,可以获得正常的解。需要注意的是,时间积分算法通常采用Lumped集中对角质量矩阵而不是一致质量矩阵,以提高计算效率。总体来说,由于计算效率的问题,隐式时间积分算法ABAQUS-Standard特别适合于低速冲击问题;对于高速冲击问题,由于存在不连续非线性接触的动响应过程,隐式算法解决不好,使用显式时间算法ABAQUS-Explicit更好。此外,对于瞬态和稳态热传导问题,半离散的抛物线方程,中心差分法可以较好获得温度分布。

对于依赖于率的粘塑性问题(对于本质上的粘性材料),与弹塑性材料的根本区别在于,一般来说是一致性条件不满足(排除弹塑性材料在高温下的软化问题,对于这种问题,屈服条件也可以满足),即屈服条件不满足,N+1时刻的物理量不用回映到N+1时刻的屈服面上,粘塑性模型成为过应力模型,显示和隐式算法都可以求解。对于依赖于率的本构模型,其可解决模拟高速冲击、爆炸、弹道射击问题时存在的动态应变局部化问题(对于动态问题,平衡方程丧失双曲线特性;对于静态问题,平衡方程失去椭圆性),解决网格尺寸效应,其实质上是引入了适当的阻尼迟滞效应。需要注意的是,对于大变形(又称为有限变形)问题,Cauchy应力率和速度梯度(包括客观和对称的扭曲张量率D、不客观和反对称的spin旋转张量W两个部分)均是不客观的,为解释刚体旋转(如纯剪切变形就包含刚体旋转),在共旋坐标系下面求解真实应力和应变,应力和应变积分求解的时候应首先求解客观性的Jaumann应力率(相对于真实应力,空间坐标系),相对于Second Piola-Kirchhoff应力是Truesdell率(材料坐标系)。ABAQUS软件对于大变形问题已经做了旋转。

对于一些耦合场问题,由于计算量非常大,同时要求解太多物理量,如热-流-固耦合,要求解位移、压力、温度,采用纯隐式算法或显式算法基本不太获得收敛或准确的解,这时候可采用混合的隐式与显式有限元格式mixed implicit-explicit partitoning方法,将刚度矩阵和阻尼矩阵分成两个部分,在同一区域采取不同算法,提高计算效率和精度、稳定性和收敛性。

开放分享：优质有限元技术文章,助你自学成才