abaqus手册_专栏 | 几何非线性的刚度矩阵求解—有限元理论基础及Abaqus内部实现方式研究系列22...

傷城~ 2023-01-01 08:55 63阅读 0赞

**作者介绍**

![9bdf28813524fe81d7af710d1815e93c.png][]

**snowwave02**

博士，高级工程师

snowwave02团队：设计仿真领域的软件开发团队，由软件、机械、物理等专业人员组成，10年以上CAD/CAE软件开发经验，精通Abaqus二次开发，承接过多个航天、航空、船舶、机械等行业大型设计仿真类项目，具有丰富的实战经验。 文末附作者团队《基础理论->Abaqus操作->matlab编程》 **免费**视频教程，讲解线性的应变在商软或者自主软件中在单元中实现。 本文由作者原创发布于技术邻平台，转载需作者授权。

**概述**

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过 (1)   基础理论 (2)   商软操作 (3)   自编程序 三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。 有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

![58422e67cb50865c3984d16a64a7aec9.png][]

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

**几何非线性的刚度矩阵求解**

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

![12c4114b9895151034deafdfd0545a5b.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][] **1.1 几何非线性的刚度矩阵推导理论** 在前面17章：几何非线性的物理含义中，我们提到如果是非线性系统，应变能W随t的变化就是个非线性过程。每个时刻点可以求出一个斜率，这个斜率最终会形成当前时刻点的刚度矩阵。

![8fad7383fba0a69e257d648803264c8a.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][] 求导后得到的刚度K：

![160a67b61fbce425b024116aa1ef461d.png][]

也就是刚度矩阵将分为两块： (1)   上式的前面一部分称为材料刚度阵，依然是以前的BDB形式，只不过B换成了当前时刻的应变位移矩阵 (2)   后面新增项一般称为几何刚度阵，在Abaqus中称为初始应力矩阵(initial stress stiffness)。

##  ##

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

##  ##

## **1.2.1 理想的求解方式** ##

##  ##

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

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

![703daae8da08eaef9d000092b0538be3.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][] ![c089c2efe18186988ee0caf126dc2b2d.png][] 在增量迭代法，此时时刻t和t+Δt分别表示为当前增量步的开始和结束时间。此时W求导表示为增量的形式：

![0dd03845bea0e20bfd1295b510070397.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][] 上面的应力应变理论上将都应该取当前增量步的结束时间的结果。先将应力增量转换为应变和本构关系，得到：

![68e3a7519afaf8d853aa9f18a6ed981e.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][]

### **1.2.2 解题困难点** ###

![42f65c1d7b397d1a77e0d38e2e4b07c1.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][]

### **1.2.3 困难点的解决方法** ###

解决上述困难点的基本原则：对一个非线性问题来说，在增量步中的斜率K即使计算不精确，只要不是偏差太大，依然收敛，而且满足精度范围。K只影响迭代的次数，不影响结果的精度。 关于这个原则，具体的解释和算例验证可以看我们下面的视频(复制至浏览器打开)： https://www.jishulink.com/college/video/c13034  Abaqus用户子程序UMat详解与开发工具(未完，待续)：2.4 整体理论02-Jacobian和应力矩阵的含义 当然这个K的近似也不能偏差太大，要不然迭代次数太多或者根本无法收敛。本着这个基本原则，两个困难点的解决方法分别如下：

![8cc57b691512729191f2ea394de85cb0.png][]

###  ###

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

###  ###

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

###  ###

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

![14c7e8423301eb5388a6e2b54fe29008.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][] 其中，Abaqus的UEL的U和dU表示的是当前增量步最后时刻的全量和增量的缘故。刚度矩阵只与增量步的量有关，迭代步的量只是对增量值的修正，不直接反应到刚度矩阵中，固UEL输入参数没有迭代相关的任何量。Abaqus无论是内部单元还是外部单元都是通过输入增量步最后时刻的位移全量和增量来求应变增量。具体的UEL说明可参看下方视频(复制至浏览器打开)： https://www.jishulink.com/college/video/c14948 深入浅出有限元：基础理论->Abaqus操作->matlab编程  (2)困难点2：应力增量的求法通过材料函数实现，内部材料函数接口和Abaqus的另一个UMAT接口是一致的，UMAT接口如下：

![4889227030b9f09d442eea0c26b53a68.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][]

其中的上述四个关键参数中，输入的STRESS是增量步开始时刻的应力全量，STRAN和DSTRAN分别是增量步开始时刻和结束时刻的应变全量和应变增量，需要通过这三个量来求出结束时刻的应力增量，

![c089c2efe18186988ee0caf126dc2b2d.png][] ，累加加到STRESS上，从而输出到程序中，而求解的过程猜测无论是Abaqus内部材料还是UMAT都采用了开始时刻的材料本构关系。譬如下方为塑性材料的UMAT的应力增量求解的matlab代码。

![f5075cb7902f980bf49b1fba2a314c1c.png][]

![c089c2efe18186988ee0caf126dc2b2d.png][]

![4e176f28b2c5392167cfc2ed70980068.png][]

##  ##

## **1.3 算例验证** ##

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

![c089c2efe18186988ee0caf126dc2b2d.png][]

![38519b6c0c754554fa32a07d35291533.png][]

两者结果相差0.08%，基本一致，从而部分验证了本文对Abaqus几何非线性下刚度矩阵实现方式的猜测，不能说完全验证，毕竟几何非线性除了上面提到的两个困难点，想要和商软完全一致，涉及的问题很多，但至少在如果上面两个困难点的解决方法和Abaqus不一致，那么得到的结果也不会一致。该算例的视频解说及iSolver操作演示如下(复制至浏览器打开)： http://www.jishulink.com/college/video/c12884   **视频解说及操作演示：** **20理论系列文章17-几何非线性的物理含义**

![01cfeee2f679d67dc05ae0758c453a31.png][]

**详见文末课程介绍。**

**优质课程**

**《自主结构有限元求解器开发框架iSolver及有限元理论介绍**》

![c1a782e9a9c436275aa55c1b45170af9.png][]

该系列视频通过算例实际操作介绍自主结构有限元求解器开发框架iSolver，并采用Abaqus快速校核结果。iSolver为一个完全自主的通用有限元程序开发框架，结果和Abaqus相当，可以快速集成客户的自研有限元算法和分析流程，也可授权作为客户软件中的一个模块共同销售，帮助客户实现自主CAE软件的商业化包装和推广。

软件具体介绍和下载传送门(复制至浏览器打开)： http://www.jishulink.com/content/post/337351 **第一部分(0~9)：单独功能介绍** 1：单元篇 2：材料篇 3：载荷和边界篇 4：分析篇 **第二部分(10~19)：案例介绍** 10：集成篇-介绍iSolver集成到其它软件 11：分析案例篇 **第三部分(20~29)：理论介绍** 20：有限元理论基础及Abaqus内部实现方式研究系列文章 **《系列免费课程》**

![dc3d97843a72309adf571d2ff2b6b462.png][]

\-END-

##  ##

#  #

##  ##

`` **![605e16fe02b77690f08d24470cec1b9b.gif][]** **点击阅读原文，查看iSolver介绍！**

[9bdf28813524fe81d7af710d1815e93c.png]: /images/20221120/5982ccbdbee140a09cee4d12714633e8.png
[58422e67cb50865c3984d16a64a7aec9.png]: /images/20221120/eb71908cf5a2401f8a76f88bcec7b27b.png
[12c4114b9895151034deafdfd0545a5b.png]: /images/20221120/01af44e01e0746ae9b1ebf474e3a08b5.png
[c089c2efe18186988ee0caf126dc2b2d.png]: /images/20221120/e0f5b4c100c54412bfddf3b4cb164348.png
[8fad7383fba0a69e257d648803264c8a.png]: /images/20221120/860612a553a040e5a8a4f3a85fcd2216.png
[160a67b61fbce425b024116aa1ef461d.png]: /images/20221120/2531ae9001964234a79c008c45f8368b.png
[703daae8da08eaef9d000092b0538be3.png]: https://img-blog.csdnimg.cn/img_convert/703daae8da08eaef9d000092b0538be3.png
[0dd03845bea0e20bfd1295b510070397.png]: https://img-blog.csdnimg.cn/img_convert/0dd03845bea0e20bfd1295b510070397.png
[68e3a7519afaf8d853aa9f18a6ed981e.png]: https://img-blog.csdnimg.cn/img_convert/68e3a7519afaf8d853aa9f18a6ed981e.png
[42f65c1d7b397d1a77e0d38e2e4b07c1.png]: https://img-blog.csdnimg.cn/img_convert/42f65c1d7b397d1a77e0d38e2e4b07c1.png
[8cc57b691512729191f2ea394de85cb0.png]: https://img-blog.csdnimg.cn/img_convert/8cc57b691512729191f2ea394de85cb0.png
[14c7e8423301eb5388a6e2b54fe29008.png]: https://img-blog.csdnimg.cn/img_convert/14c7e8423301eb5388a6e2b54fe29008.png
[4889227030b9f09d442eea0c26b53a68.png]: https://img-blog.csdnimg.cn/img_convert/4889227030b9f09d442eea0c26b53a68.png
[f5075cb7902f980bf49b1fba2a314c1c.png]: https://img-blog.csdnimg.cn/img_convert/f5075cb7902f980bf49b1fba2a314c1c.png
[4e176f28b2c5392167cfc2ed70980068.png]: https://img-blog.csdnimg.cn/img_convert/4e176f28b2c5392167cfc2ed70980068.png
[38519b6c0c754554fa32a07d35291533.png]: https://img-blog.csdnimg.cn/img_convert/38519b6c0c754554fa32a07d35291533.png
[01cfeee2f679d67dc05ae0758c453a31.png]: https://img-blog.csdnimg.cn/img_convert/01cfeee2f679d67dc05ae0758c453a31.png
[c1a782e9a9c436275aa55c1b45170af9.png]: https://img-blog.csdnimg.cn/img_convert/c1a782e9a9c436275aa55c1b45170af9.png
[dc3d97843a72309adf571d2ff2b6b462.png]: https://img-blog.csdnimg.cn/img_convert/dc3d97843a72309adf571d2ff2b6b462.png
[605e16fe02b77690f08d24470cec1b9b.gif]: https://img-blog.csdnimg.cn/img_convert/605e16fe02b77690f08d24470cec1b9b.gif