上一篇我们提到,在土力学领域的有限元计算中,计算机每一步的计算本质都是在求解下面的方程:

即每一个步长的刚度矩阵乘以节点位移向量等于荷载向量,在一个步长内,我们假定刚度矩阵不变,用优化的牛顿法(Modified Newton-Raphson Method)迭代可以一步步求出收敛的数值解。每一步的刚度矩阵都是在运算开始时确定的。

优化的牛顿拉夫森法运算示意图

那么问题来了,这个每个步长的整体刚度矩阵

要怎么算?

整体刚度矩阵

其实是很多单元刚度矩阵(Element Stiffness Matrix)的合成,合成的规则在此省略一万字,相信大家都有过做题的经历(画大矩阵的痛苦回忆)。每个单元刚度矩阵,写在母单元里面长这样:

其中的核心矩阵

取决于材料的本构模型。譬如我们最熟悉的就是理想线弹性固体材料(Linear Elastic Solid Element),该矩阵满足如下关系,
,其中,
是基于线弹性本构模型的刚度矩阵。在平面应变问题里面该式写开来就是:

E是弹性模量,μ自然就是泊松比啦。那其他种类单元的D矩阵长啥样?比如膜单元的

,其实是这样的:
,这个只有二阶是因为膜单元上面每个节点只有两个自由度……其他形式譬如梁单元(Beam Element)的单元刚度矩阵可以在书上查阅(是五阶的,一根梁上各点有五个自由度哦)。这些虽然是不同的单元类型,因为都是线弹性的,所以只要确定杨氏模量和泊松比,那刚度矩阵就基本确定了。

那么非线性材料怎么办?比如土本来就是一种弹塑性材料,其刚度一定是非线性的。我们以大家耳熟能详的剑桥模型(Cam Clay Model)为例,其刚度矩阵满足如下关系:

其中K是体积模量,

是偏应力和平均主应力的比值,
是平均主应力和偏应力的比值,

等……等一下,为什么画风突然复杂了起来?刚刚说好的线弹性材料D矩阵不变来着的,现在换了个非线性的剑桥模型,刚度矩阵一下子跟p和q相关了,换句话说我们需要知道当前应力状态才能确定当前的刚度矩阵。那就是说在有限元计算时,每算一个步长,刚度矩阵就要重新确定一次。而且进一步来讲,即便在一个步长内,随着p和q的变化,我们的刚度矩阵理论上来讲也是在不断变化的!

哦,那当前应力状态怎么确定呢?我们仔细观察了一下上式,发现其本质好像是个微分方程,所以我们唯一能做的就是在已知边界条件的情况下,求解微分方程

解这个微分方程的目标就是为了在这个步长给定了初始刚度矩阵的前提下,通过逐步推导,确定下一个步长的初始刚度矩阵

一个普通的微分方程

,胡乱举个例子,比如
这种形式,对于计算机来讲,不光函数式写不出来,其精确解更是相当难算的,只能借助于边界条件去靠近。这里先介绍一下一种大一学生秒懂的笨办法,欧拉法(Euler Scheme),它是一种很不精确的计算微分方程的方法。

比如函数y在

区间上满足
,并且有
这个初始的条件(我们管这个叫边界条件),目标是计算
,那么我们必须先要确定在此区间内一系列等间距的的函数点值。即
, 假设这些点间隔都为h,从左端
开始,欧拉法就是简单粗暴地通过计算
就能够把
估算出来啦!放进上述算式那就是
。其中h当然是
了,是一个很小的数。既然我们知道了
,那么接下去
等等在区间内一系列的点都可以依样画葫芦地确定了,即每一个点都可以通过它前面的点推算出来。最后就像多米诺骨牌一样,我们借助最后一块骨牌
终于可以确定
了,从而终于可以确定下一个初始刚度矩阵了(气喘吁吁)。
用第一块骨牌(本步长的初始刚度矩阵)来推出最后一块骨牌(下一个步长的刚度矩阵)

说到欧拉法,我总能想起很多人小时候玩过的一个游戏,叫做“传令兵”。在教室里第一个人说一句命令,然后转述给下一个人,一直到最后一个同学公布命令,看看和第一位同学有多少差距。由于每次转述都能产生误差,并且该误差能以同阶形式进行积累(而非高一阶),因而最后结果的精确度可见一斑。正所谓一步错步步错。用数学语言来说,就是“其局部截断误差为二阶,因而欧拉法只有一阶精度”。

对于学有限元的小伙伴来说,欧拉法也并不是一无是处,毕竟相比其他的算法,它占用计算机内存很小……几秒钟就算完了(先不提算得对不对,你就说快不快吧)。

但是考虑到其精度感人,所以欧拉法的改进方法(Modified Euler Scheme)就呼之欲出啦。改进的欧拉法其实就是在欧拉法基础上多算了一步。比如我们借助左端已知的

想要算第一个点
,那么我们先用欧拉法算一下,得到一个不那么精确的函数值,我们暂且叫它试验值
,然后把这个值代入下式可以算出更为精确的
:

在上述例子中,

随后以此类推,得到更精确的

最后推算出

用这种运算方法得到的结果具有二阶精度,比欧拉法要高一阶,因而在有限元运算中十分受欢迎。至于为什么通过这多一步的操作,我们就可以得到更高一阶精度,这点可以请大家在课后查阅论证。

因而,在土力学有限元应用里,改进欧拉法里面的那个间隔“h”其实就是一个小小的应变增量,把一个步长,分解成许多个小应变增量,运用该算法(Substepping algorithm)就可以计算该步长内的刚度矩阵变化了。

python 拟牛顿法 求非线性方程_有限元简单科普之——改进的欧拉法相关推荐

  1. python 拟牛顿法 求非线性方程_C语言实现迭代法求非线性方程的根

    迭代法求非线性方程的根 迭代法是一种逐次逼近法.它是求解代数方程,超越方程及方程组的一种基本方法,但存在收敛性及收敛快慢的问题. 为了用迭代法求非线性方程f(x) = 0的近似根: 1.首先需要将此方 ...

  2. python 拟牛顿法 求非线性方程_9-非线性优化

    本章内容: 介绍了无约束和有约束两种类型的非线性优化 一.无约束优化 无约束优化的一般形式为: 如果要求最大值,则 1.1 fminunc() 介绍:是MATLAB求解无约束优化的主要函数,算法有:信 ...

  3. python怎么求指数_求指数 python

    softmax用于多分类过程中最后一层,将多个神经元的输出,映射到(0, 1)区间内,可以看成概率来理解,从而来进行多分类! softmax函数如下: 更形象的如下图表示: softmax 直白来说就 ...

  4. python打造excel神器_超简单:用Python让Excel飞起来

    前言 如何获取学习资源 章Python快速上手 1.1为什么要学习用Python控制Excel 1.2Python编程环境的搭建 1.2.1安装Python官方的编程环境IDLE 1.2.2安装与配置 ...

  5. python编程求导数_面向对象编程 —— java实现函数求导

    首先声明一点,本文主要介绍的是面向对象(OO)的思想,顺便谈下函数式编程,而不是教你如何准确地.科学地用java求出函数在一点的导数. 一.引子 defd(f) :defcalc(x) : dx= 0 ...

  6. python 表达式求值_简单算术表达式求值

    本文主要探讨简单的数学算术表达式求值算法的原理和实现. 1. 约束 本文只是探讨简单的算术表达式的求值算法,为了将主要精力放在算法思想的探讨和实现上,避免陷入对其他不是直接相关的细节的过多思考,所以提 ...

  7. python编程求导数_用python怎么计算导数最简单?

    谢邀,请恕我微积分学得不扎实,我记得常数的一阶导数均为0. 如果列表中传入的为含变量x的式子,代码可能如下. from sympy import Symbol, diff x = Symbol('x' ...

  8. python迭代法求解非线性方程_荐【数学知识】非线性方程求解的二分法以及牛顿迭代法...

    [数学知识]非线性方程求解的二分法以及牛顿迭代法 本博客不谈及理论推导,只提供代码实现,关于理论推导,大家可以查看其它博客文章. 导入包 import sys import math import s ...

  9. python列表求平均值_长篇文讲解:Python要求O(n)复杂度求无序列表中第K的大元素实例...

    本文内容主要介绍了Python要求O(n)复杂度求无序列表中第K的大元素实例,具有很好的参考价值,希望对大家有所帮助.一起跟随小编过来看看吧! 昨天面试上来就是一个算法,平时基本的算法还行,结果变个法 ...

最新文章

  1. 启动Genymotion时报错Failed to initialize backend EGL display
  2. 每日一皮:听说学琵琶的都很文弱...
  3. Java中“==”的使用,以及“==”和equal的比较
  4. wenbao与windows命令
  5. 单网段DHCP服务器的架设
  6. 为什么截屏不能分享微信_为什么腾讯可以在移动端QQ做到闪照,而在Windows桌面端做不到?...
  7. openai-gpt_您可以使用OpenAI GPT-3语言模型做什么?
  8. 前端学习(2932):vue中的v-show
  9. adb无法连接安卓手机
  10. 数码管显示实验一 编写程序让8只数码管同时显示零
  11. EularProject 39:给周长推断构成直角三角形个数
  12. .net core sorteddictionary 排序_#键盘排序——为什么我们的键盘字母不是按照ABCD的顺序排列?...
  13. 传智播客Java 关键字,标识符,注释
  14. 菜鸟学习oracle一看就会
  15. php 漏洞扫描,10个最佳PHP代码安全扫描程序来查找漏洞
  16. 智能家居语音控制系统项目毕业答辩
  17. mysql 添加表字段并添加数据,MySQL为表的所有字段添加数据
  18. MindAR初体验——一款js实现的AR库
  19. Linux自建RustDesk中继服务器
  20. STM32机器人控制开发教程No.3 使用遥控控制电机/舵机(基于HAL库)

热门文章

  1. 解决方案:SpringBoot分布式项目跨域
  2. IO多路复用中select、poll、epoll之间的区别
  3. nginx设置http强制跳转https
  4. dubbo服务降级与限流
  5. PYTHON——多线程:Thread类与线程函数
  6. Struts2 自定义验证器
  7. [译]CSV 注入:被人低估的巨大风险
  8. 阿里云释放数据能力 开启大数据元年
  9. Acoustic Echo Cancellation (AEC) 回音消除技术探索
  10. 测测你的显示器灰阶显示