有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解

(原创,转载请注明出处)

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图1有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图2==概述==

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图3本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式,通过

(1)   基础理论

(2)   商软操作

(3)   自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来,同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟,商用有限元软件同样也是采用这些成熟的有限元理论,只是在实际应用过程中,商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题,且各家软件的修正方法都不一样,每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式,但都只是提一下用什么方法修正,很多没有具体的实现公式。商用软件对外就是一个黑盒子,除了开发人员,使用人员只能在黑盒子外猜测内部实现方式。

0.png

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法,另一方面我们自己编程实现结构有限元求解器,通过自研求解器和商软的结果比较来验证我们的猜测,如同管中窥豹一般来研究的修正方法,从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元,所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式,同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨,同时由于水平有限可能有许多的理论错误,欢迎交流讨论,也期待有更多的合作机会。

自主结构有限元求解器iSolver介绍视频:

1.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图6

                                             http://www.jishulink.com/college/video/c12884

==第22篇:几何非线性的刚度矩阵求解==

几何非线性在界面上是很容易设置的,但商软内部的处理相当复杂,我们从最基本的刚度矩阵的求解出发,看看在几何非线性设置后,刚度矩阵具体是怎么实现的。本文首先介绍几何非线性下的刚度矩阵的理论推导和计算机求解方法,说明理想的求解方式的困难点和猜测Abaqus内部的解决方法。最后利用一个简单的算例通过对比iSolver和Abaqus的结果,部分验证我们对Abaqus几何非线性的刚度矩阵的实现方式的猜测。

2.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图8

1.1 几何非线性的刚度矩阵推导理论

在前面17章:几何非线性的物理含义中,我们提到如果是非线性系统,应变能W随t的变化就是个非线性过程。每个时刻点可以求出一个斜率,这个斜率最终会形成当前时刻点的刚度矩阵。

QQ图片20200604193539.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图10

求导后得到的刚度K:

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图11

1.png

也就是刚度矩阵将分为两块:

(1)   上式的前面一部分称为材料刚度阵,依然是以前的BDB形式,只不过B换成了当前时刻的应变位移矩阵

(2)   后面新增项一般称为几何刚度阵,在Abaqus中称为初始应力矩阵(initial stress stiffness)。

1.2 几何非线性的刚度矩阵计算机求解

1.2.1 理想的求解方式

理论上受力曲线是一条光滑曲线,计算机没法求解曲线上每个时刻点的结果,只能求解部分有限间隔点的结果。非线性问题不是一条直线,所以需要多次迭代才能实现。

非线性问题求解有多种方法主要分为以下几类:增量法、迭代法、增量迭代法和弧长法。

 

2.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图14有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图15

在增量迭代法,此时时刻t和t+Δt分别表示为当前增量步的开始和结束时间。此时W求导表示为增量的形式:

3.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图17

上面的应力应变理论上将都应该取当前增量步的结束时间的结果。先将应力增量转换为应变和本构关系,得到:

4.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图19

1.2.2 解题困难点

2.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图21

1.2.3 困难点的解决方法

解决上述困难点的基本原则:对一个非线性问题来说,在增量步中的斜率K即使计算不精确,只要不是偏差太大,依然收敛,而且满足精度范围。K只影响迭代的次数,不影响结果的精度。关于这个原则,具体的解释和算例验证可以看我们下面的视频:

https://www.jishulink.com/college/video/c13034 Abaqus用户子程序UMat详解与开发工具(未完,待续):2.4 整体理论02-Jacobian和应力矩阵的含义

当然这个K的近似也不能偏差太大,要不然迭代次数太多或者根本无法收敛。本着这个基本原则,两个困难点的解决方法分别如下:

3.png

1.2.4 Abaqus中的内部实现方式猜测

当然,我们不知道Abaqus的内部计算方式,但可以通过两个外部接口可以猜测Abaqus的内部实现方式。

(1)困难点1:应变增量与增量步结束时刻的位移增量的关系,Abaqus可以通过单元函数实现,内部单元函数接口和Abaqus的子程序UEL的接口是一致的,如下表示:

6.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图24

其中,Abaqus的UEL的U和dU表示的是当前增量步最后时刻的全量和增量的缘故。刚度矩阵只与增量步的量有关,迭代步的量只是对增量值的修正,不直接反应到刚度矩阵中,固UEL输入参数没有迭代相关的任何量。Abaqus无论是内部单元还是外部单元都是通过输入增量步最后时刻的位移全量和增量来求应变增量。具体的UEL说明可参看下方视频:

https://www.jishulink.com/college/video/c14948 深入浅出有限元:基础理论->Abaqus操作->matlab编程

 (2)困难点2:应力增量的求法通过材料函数实现,内部材料函数接口和Abaqus的另一个UMAT接口是一致的,UMAT接口如下:

7.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图26

其中的上述四个关键参数中,输入的STRESS是增量步开始时刻的应力全量,STRAN和DSTRAN分别是增量步开始时刻和结束时刻的应变全量和应变增量,需要通过这三个量来求出结束时刻的应力增量有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图27,累加到STRESS上,从而输出到程序中,而求解的过程猜测无论是Abaqus内部材料还是UMAT都采用了开始时刻的材料本构关系。譬如下方为塑性材料的UMAT的应力增量求解的matlab代码。

8.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图29

2.png

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图31有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图321.3 算例验证

从UEL和UMAT两个接口只能猜测Abaqus内部单元几何非线性的刚度矩阵的实现方式,为了证明我们的猜想,我们创建了一个几何非线性的悬臂梁的简单算例,采用iSolver和Abaqus分别计算,而iSolver的两个困难点按我们前面所述实现,最终的位移结果如下:

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图33

2.png

两者结果相差0.08%,基本一致,从而部分验证了本文对Abaqus几何非线性下刚度矩阵实现方式的猜测,不能说完全验证,毕竟几何非线性除了上面提到的两个困难点,想要和商软完全一致,涉及的问题很多,但至少在如果上面两个困难点的解决方法和Abaqus不一致,那么得到的结果也不会一致。该算例的视频解说及iSolver操作演示如下:

http://www.jishulink.com/college/video/c12884 20理论系列文章17-几何非线性的物理含义。

1.4 联系方式

如果有任何其它疑问或者项目合作意向,也欢迎联系我们:

SnowWave02 From www.jishulink.com

email: snowwave02@qq.com

有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图351.5 以往的系列文章

以往的系列文章:

1.5.1 ========第一阶段========

第一篇:S4壳单元刚度矩阵研究。介绍Abaqus的S4刚度矩阵在普通厚壳理论上的修正。

http://www.jishulink.com/content/post/338859

第二篇:S4壳单元质量矩阵研究。介绍Abaqus的S4和Nastran的Quad4单元的质量矩阵。

http://www.jishulink.com/content/post/343905

第三篇:S4壳单元的剪切自锁和沙漏控制。介绍Abaqus的S4单元如何来消除剪切自锁以及S4R如何来抑制沙漏的。

http://www.jishulink.com/content/post/350865

第四篇:非线性问题的求解。介绍Abaqus在非线性分析中采用的数值计算的求解方法。

http://www.jishulink.com/content/post/360565

第五篇:单元正确性验证。介绍有限元单元正确性的验证方法,通过多个实例比较自研结构求解器程序iSolver与Abaqus的分析结果,从而说明整个正确性验证的过程和iSolver结果的正确性。

https://www.jishulink.com/content/post/373743

第六篇:General梁单元的刚度矩阵。介绍梁单元的基础理论和Abaqus中General梁单元的刚度矩阵的修正方式,采用这些修正方式可以得到和Abaqus梁单元完全一致的刚度矩阵。

https://www.jishulink.com/content/post/403932

第七篇:C3D8六面体单元的刚度矩阵。介绍六面体单元的基础理论和Abaqus中C3D8R六面体单元的刚度矩阵的修正方式,采用这些修正方式可以得到和Abaqus六面体单元完全一致的刚度矩阵。

https://www.jishulink.com/content/post/430177

第八篇:UMAT用户子程序开发步骤。介绍基于Fortran和Matlab两种方式的Abaqus的UMAT的开发步骤,对比发现开发步骤基本相同,同时采用Matlab更加高效和灵活。

https://www.jishulink.com/content/post/432848

第九篇:有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图36有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图37有限元理论基础及Abaqus内部实现方式研究系列22: 几何非线性的刚度矩阵求解的图38编写线性UMAT Step By Step。介绍基于Matlab线性零基础,从零开始Step by Step的UMAT的编写和调试方法,帮助初学者UMAT入门。

http://www.jishulink.com/content/post/440874

第十篇:耦合约束(Coupling constraints)的研究。介绍Abaqus中耦合约束的原理,并使用两个简单算例加以验证。

https://www.jishulink.com/content/post/531029

1.5.2 ========第二阶段========

第十一篇:自主CAE开发实战经验第一阶段总结。介绍了iSolver开发以来的阶段性总结,从整体角度上介绍一下自主CAE的一些实战经验,包括开发时间预估、框架设计、编程语言选择、测试、未来发展方向等。

http://www.jishulink.com/content/post/532475

第十二篇:几何梁单元的刚度矩阵。研究了Abaqus中几何梁的B31单元的刚度矩阵的求解方式,以L梁为例,介绍General梁用到的面积、惯性矩、扭转常数等参数在几何梁中是如何通过几何形状求得的,根据这些参数,可以得到和Abaqus完全一致的刚度矩阵,从而对只有几何梁组成的任意模型一般都能得到Abaqus完全一致的分析结果,并用一个简单的算例验证了该想法。

http://www.jishulink.com/content/post/534362

第十三篇:显式和隐式的区别。介绍了显式和隐式的特点,并给出一个数学算例,分别利用前向欧拉和后向欧拉求解,以求直观表现显式和隐式在求解过程中的差异,以及增量步长对求解结果的影响。

http://www.jishulink.com/content/post/537154

第十四篇:壳的应力方向。简单介绍了一下数学上张量和Abaqus中壳的应力方向,并说明Abaqus这么选取的意义,最后通过自编程序iSolver来验证壳的应力方向的正确性。

https://www.jishulink.com/content/post/1189260

第十五篇:壳的剪切应力。介绍了壳单元中实际的和板壳近似理论中的剪切应力,也简单猜测了一下Abaqus的内部实现流程,最后通过一个算例来验算Abaqus中的真实的剪切应力。

https://www.jishulink.com/content/post/1191641

第十六篇:Part、Instance与Assembly。介绍了Part、Instance与Assembly三者之间的关系,分析了Instance的网格形成原理,并猜测Abaqus的内部组装实现流程,随后针对某手机整机多part算例,通过自编程序iSolver的结果比对验证我们的猜想。

https://www.jishulink.com/content/post/1195061

第十七篇:几何非线性的物理含义。介绍了几何非线性的简单的物理含义,并通过几何非线性的悬臂梁Abaqus和iSolver的小应变情况的结果,从直观上理解几何非线性和线性的差异。

https://www.jishulink.com/content/post/1198459

第十八篇:几何非线性的应变。首先从位移、变形和应变的区别说起,然后通过一维的简单例子具体介绍了几何非线性下的应变的度量方式,并给出了工程应变、 真实应变、Green应变三者一维情况下在数学上的表达方式。

https://www.jishulink.com/content/post/1201375

第十九篇:Abaqus几何非线性的设置和后台。首先介绍了几何非线性一般的分类,然后详细说明了Abaqus中几何非线性的设置方式和常用单元的分类,最后以一个壳单元的简单算例为对象,可以发现应变理论、Abaqus和iSolver三者在线性、小应变几何非线性和大应变几何非线性三种情况下都完全一致,从而验证Abaqus几何非线性后台采用的应变和我们的预想一致。

http://www.jishulink.com/content/post/1203064

第二十篇:UEL用户子程序开发步骤。本文首先简单的讨论了UEL的一般含义,并详细的介绍了基于Fortran和Matlab两种方式的UEL的开发步骤,对比发现开发步骤基本相同,但Matlab更加高效和灵活。

https://www.jishulink.com/content/post/1204261

1.5.3 ========第三阶段========

第二十一篇:自主CAE开发实战经验第二阶段总结。从实战角度介绍自主CAE在推广和工程化应用的过程中的体会,同时说明一个CAE平台最重要的两个特点:可扩展和易维护。

https://www.jishulink.com/content/post/1204970

(21条)
默认 最新
老鼠,我想问一下,在一个增量步内只有一次迭代,UEL却调用了两次。如果说最后一次调用是为了确定下一个增量步的初始刚度矩阵,但是abaqus在新的增量步开始时是还会调用一次UEL,并且增量步内第一次调用的UEL中DU和上衣一增量步结束时的一致,这就导致新增量步初始时的U和上一增量步结束时的U不同,且刚好差一个DU。这对于大变形材料,比如橡胶会导致计算的增量步初始刚度矩阵出现误差,我想问一下这是什么原因造成的呢,还是UEL本来就是这么调用的,麻烦老师了。
评论 4 点赞
回复
好的我试一下
评论 点赞
回复
首先解释一下ABAQUS除了第一个初始增量步以外,对于一个线性问题,每个增量步依然需要调用两次uel呢?首先你在监视器中看到的1是不平衡迭代的次数,不平衡迭代1次就是迭代了2次。ABAQUS帮助文档中明确说了它只有在迭代2次后才会进入是线性方程还是非线性方程的判断,也就是说最开始它是先把所有方程都当非线性方程求解的,对于非线性方程收敛的判断有两个,一个是rhs<容差,一个是Du<容差(注意这里的Du指迭代步前后的Du,uel中的Du是本次增量步内所有迭代步的累加值),也就是说对于一个即使是线性方程的问题,第一次调用uel即使达到了rhs收敛,但Du一定不收敛,所以会迭代第二次,也就是你所说的两次迭代Du相同,这说明Du收敛了。 说完这些再来说初始增量步,如果你足够细心你会发现初始增量步也是调用两次UEL,但初始增量步因为不能继承上一个增量步的刚度矩阵K,因此初始增量步的第一次调用uel是只计算amatrx不计算du的,第二次调用才使用第一次调用得到的amatrx计算Du,那么问题来了,按照Du判断收敛,这依然不满足标准,为什么ABAQUS不进行第三次调用呢?因为调用第二次之后ABAQUS会判断是否为线性方程,如果是则直接判断收敛。因此求解线性方程,初始增量步也只调用两次uel。
评论 1 点赞
查看全部4条回复 >
老师 我不是力学专业的 看了您的课后很是受用 我想问编写屈服准则之前我应该看力学方面的哪些书籍 像刚度矩阵 这类的知识我都不知道在哪入手 您能给些建议吗
评论 1 点赞
回复
如果是初学者可以看logan的那本有限元入门的书
评论 点赞

查看更多评论 >

点赞 50 评论 30 收藏 7
关注