有限差分法matlab程序_专栏 | 显式求解Step By Step—有限元理论基础及Abaqus内部实现方式研究系列24...

作者介绍

snowwave02

博士，高级工程师

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

概述

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

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

显式求解 Step By Step

动力学问题是将力的方程和运动学方程耦合在一起的理论，在实际问题中，运动过程无时不刻的与力相关，因此数学上两个方程联立求解是最理想的。但数值计算中却难以实现，也没必要联立求解，只需先求一个方程，然后再求另一个方程，只要时间步长足够短，那么精度依然可以保证。有限元亦是如此解耦，有限元两个时刻点之间会认为力的状态不变，那么可以分两步求解。(1)   在每个时刻点只求力学方程，得到运行学的初始条件。(2)   在两个时刻点之间由于力的状态不变，那么可以只按运动学来求的运动学量。本文首先研究商用有限元中最常用的显式求解算法中心差分法的理论，并给出了一个弹簧显式动力学分析的Step by Step例子，通过这个例子猜测了Abaqus中采用中心差分法求解显式动力学问题的过程，然后在自编有限元程序iSolver采用同样的算法，验证iSolver的结果和Abaqus完全一致，从而证明Abaqus的内部算法和我们设想的一致。具体验证过程也可以参考我们的演示录像(复制至浏览器打开)：https://www.jishulink.com/college/video/c128844分析篇.3-弹簧显示动力学分析

1.1 中心差分法的理论

中心差分法的标准理论可查看相应的论文，由于和Abaqus的中心差分法比较接近，所以在此不累述。

1.2 Abaqus中心差分法的理论

注：本节公式均摘自《Abaqus Theory Manual 2.4.5 Explicit dynamic analysis》1.2.1 差分公式

取i-0.5时刻的速度和时刻的加速度，则有下式，

式中，

表示速度，

表示加速度，

表示时间步大小。带入i时刻的位移，则有，

其中，i时刻的加速度可根据牛顿定律计算得出，即，                                                 

式中，M表示集中质量阵，F表示外力，I表示内力。针对初始化、处理某些约束条件以及在后处理过程中，Abaqus对速度有特殊的处理，如下式，                             

1.2.2 初始化

在初始时，即(t=0)，除非用户自定义，一般情况下速度和加速度都设为0。假定有下式，                                         

带入式(1)，则有                                         

将i=-1带入式(4)中，也可得到式(6)，故前后是一致的。初始化流程如下，

1.2.3 整体流程

以某一时刻t为例，整体流程如下：

• 其中，内力和时间步是同时求的，这个我们会在后续的VUEL文章中讲到。需要特别注意初始化时内力和时间步对应时刻点和整体流程的差别，因此我们调整了求时间步和内力的顺序以和Abaqus一致。

1.3 Abaqus的实现验证

1.3.1 模型例子：弹簧的显式动力学分析

Example 1：弹簧的显式动力学分析(模型详见附件Job-ExplicitSpring-AddF-NoBC.inp)创建一根线弹簧，在弹簧两侧各加两个点质量，无约束，右端X方向拉力。参数如下：尺寸：X方向长度L=1；质量：各向同性，大小为100；刚度：弹簧刚度为1000；力：分析时间内恒定大小，为800；时间设置：总时间为1，时间增量固定为0.2。

1.3.2 增量步零(0s)

1.3.2.1 加速度对比

1.3.2.2 速度对比

1.3.2.3 位移对比

1.3.3 增量步一(0.2s)

1.3.3.1 加速度对比

1.3.3.2 速度对比

1.3.3.3 位移对比

1.3.4 增量步二、三、四

与第一个增量步类似，不再累述。

1.3.5 增量步五(1.0s)

1.3.5.1 加速度对比

1.3.5.2 速度对比

1.3.5.3 位移对比

1.3.6 验证结果

可以发现iSolver和Abaqus完全一致。

总结

本文概要性地介绍了Abaqus中心差分法的理论以及算法实现的整体流程，并通过简单弹簧显式动力学分析算例与Abaqus计算结果进行对比，验证了算法和整体流程的正确性。后续文章，我们会逐步深入显式动力学的一些细节，敬请期待。如果有任何其它疑问或者项目合作意向，也欢迎联系我们：snowwave02 From www.jishulink.comemail: snowwave02@qq.com

优质课程

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

该系列视频通过算例实际操作介绍自主结构有限元求解器开发框架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内部实现方式研究系列文章《系列免费课程》

-END-

点击阅读原文，查看iSolver介绍！