非线性方程组的解法
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
非线性方程组的解法
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 几何非线性 物理位置：平衡方程必须按照变形后的几何位置建立 数学描述：物体的变形一般由位移的一阶微分求得， 当变形很大，高阶微分不能忽略 举例：
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 材料非线性 混凝土，钢筋，粘结滑移，时变效应 一般认为，混凝土开裂即存在材料非线性 其他非线性：边界非线性 常见的是接触问题
[K ][∆δ ] = [∆P ]
n n n
[∆δ n ] = [K n ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：……. 缺点：存储两个刚度，存贮量大
非线性方程组的解法
逐步增量法－中点刚度增量法
施加一半荷载，求出相应的位移（中点位移）
[K n−1 ]∆δ n − 1 = 1 [∆Pn ]
非线性方程组的解法
逐步搜索法 可求出极限荷载，无法求解
步骤：1. 加一级荷载 ∆P ，结果发散，退回 2. 加 1 ∆P ，若结果发散，再次退回
2
P −δ
下降段
3. 加
1 ∆P ，……至收敛，得到极限荷载 4
非线性方程组的解法
虚加刚性弹簧法 加虚拟弹簧, 改变原结构形式，新的结构形式 曲线在一定范围内没有下降段
[B] = [M ][R][Q]
其他数值方法
刚体弹簧元
k n [D] = 0
Tresca屈服条件
0 ks
本构矩阵 (2x2)
平面应变下的莫尔-库仑屈服条件
平面应力下的莫尔-库仑屈服条件
其他数值方法
刚体弹簧元
其他数值方法
无网格法
特点 在计算域上，用离散的点拟合场函数，采用基于节点（非单元节点） 的位移插值函数 部分或彻底地消除网格，避免了有限单元法中单元网格的限制 在处理大变形、开裂等问题时，有其优势 采用离散的点拟合场函数的方法 核（权）函数法 移动最小二乘法 径向基函数法 然后建立方程求解
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ， 即时切线刚度 [K 1 ]
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K1 ]−1 [∆P1 ]
[δ 2 ] = [δ1 ] + [∆δ 2 ]
[δ n ] = [δ n−1 ] + [∆δ n ]
节点荷载 [P2 ] ， 即时切线刚度 [K 2 ]
非线性方程组的解法
添加行列法
非线性方程组的解法
强制迭代法 结构总刚是负定的,只要施加负荷载,可求得相应的位移增量
问题:适用范围?
非线性方程组的解法
硬化刚度法 人为增大结构刚度,改善收敛性 一般结构切线刚度可分为
KT = K E − K p
弹性 非线性影响 刚度 刚度
K T = K E − λK p
能量准则
([P ] [δ ]) res
T
两曲线的接近如何描述？
非线性方程组的解法
两个问题： 考虑结构负刚度的一些算法 弹塑性单元应力调整计算
非线性方程组的解法
考虑结构负刚度的一些算法 问题：混凝土结构 P − δ 下降段，由力求位移求解困难 原因：上升段，刚度矩阵正定，线性方程组有解 下降段，刚度矩阵非正定，求解困难 负刚度：非正定的刚度矩阵 几种方法：逐步搜索法，虚加刚性弹簧法，位移加载法， 强制迭代法，硬化刚度法，弧长法
刚体弹簧元 刚体＋弹簧 自由度少，总刚阶数少 可在在极限分析领域运用 思路 弹簧单元刚度矩阵 界面应力应变关系 界面法向和切向的应力应变关系 界面法向和切向刚度
其他数值方法
刚体弹簧元
[K ] = ∫ [B ] [D ][B ]ds
T
k n [D] = 0
0 ks
单元刚度矩阵 (6x6) 本构矩阵 (2x2) 形函数矩阵 (2x6)
网格大小与数目的选择 块体大小，计算精度
其他数值方法
公共面法 块体间接触判断 公共面的定义（如何得到） 假象面 令两块体向薄板靠近 两块体在接触时停止运动 薄板被固定在一个位置，该面即为公共面 块体与公共面间的接触关系
其他数值方法
公共面法 接触类型
其他数值方法
颗粒离散元法 多边形、多面体 圆形、球形 可以以链杆或弹簧连接颗粒 质点－桁架模型
2
2
δ 1 = [δ n −1 ] + ∆δ 1 n− 2 n− 2
K n− 1 确定相应的刚度（中点刚度） 2
以该刚度作为初始刚度求当前位移增量
K n − 1 [∆δ n ] = [∆Pn ] 2
优点: 比两步欧拉折线法精度高
P −δ
问题:在何处加弹簧? 弹簧刚度如何?
非线性方程组的解法
位移控制法 支座位移控制加载,可求支座的反力, 和全部 P − δ 曲线 两个算法:刚度矩阵重新排列法, 添加行列法
非线性方程组的解法
添加行列法
[p1 p2]T为参考矢量, ∆λ 为荷载控制步长系数, [R1 R2]T为不平衡力 问题:外力不按比例加载?
其他数值方法
离散单元法 刚体弹簧元法 无网格法
其他数值方法
离散单元法 Cundall于1971年提出 基本原理:牛顿第二运动定律 建立每个刚体的力平衡方程, 采用动力松弛法求解 特点 岩体或颗粒组合体被模拟成通过边或角相互接触而产 生相互作用 块体之间边界的相互作用可以体现其不连续性和节理 的特征 使用显式积分迭代法,允许大位移、转动和非线性材料
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 割线迭代法
[δ 1 ] = [K 0 ]−1 [P]
[δ 2 ] = [K 1 ]−1 [P ] …… [δ n ] = [K n−1 ]−1 [P]
割线刚度 [K 1 ] 割线刚度 [K 2 ]
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 切线迭代法
有
δFs = K s ∆δ s
t + ∆t
Fs = t F + δFs ≤ c + t Fn tan ϕ
文献中给出了Ks和Kn的建议值
其他数值方法
运动方程 位移和转角、一阶导（速度）、二阶导（加速度） 采用中心差分法，逐步迭代，至达到平衡状态
考虑阻尼（Rayleigh阻尼）
阻尼系数的选择一般依赖经验方法
…
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 等刚度迭代法
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ，
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K 0 ]−1 [∆P1 ]
[δ 2 ] = [δ 1 ] +其他数值方法
无网格法
无网格伽辽金法
由移动最小二乘法得到近似场函数 由变分原理得到基本方程
取
问题:适用范围?
刚度硬化 系数 0 ≤ λ ≤ 1
非线性方程组的解法
弧长法
[K ][∆u ] = ∆λ [P] + [R ]
i i +1 i +1 i
荷载控制 步长系数 同时控制 ∆λ 和 ∆u ,使得
不平衡力
[∆u ]T [∆u ] + ∆λ2 = dS
i +1
在 λ − u 空间中,向量r迭代按一圆弧进行,故称弧长法 ……可以求得 ∆λ
其他数值方法
时步选择 具有质量m和刚度系数k的弹性系统，其运动方程
时间步长需满足
∆t ≤ 2 m / k
块体系统由多个质量不一、刚度不同的块体组成， 理论上取其最小的 ∆tmin，实际计算时取小于 ∆tmin
其他数值方法
接触判断算法 二维相对成熟，三维有待发展 分格接触检索算法 1.将研究区域划分成网格 2.在每个块体内，按照接触算法，检索块体是否接触 3.每次迭代，均需接触检索
4. 区分不同情况,计算应力增量
非线性方程组的解法
弹塑性单元应力调整计算——二维
其他数值方法
其他数值方法
有限单元法的局限 基于连续体力学,相邻边界上位移协调, 对于处理应力或位移间断的问题麻烦. 目前解决办法: 1.重新划分网格 计算工作量大 2.改变单元位移插值函数使单元各向异性 网络依赖性,精度差 3.生死单元 网格太细密
[K n−1 ][∆δ n ] = [∆Pn ]
[∆δ n ] = [K n−1 ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：1.施加第n步荷载增量 [∆Pn ]，利用始点刚度矩阵 [K n−1 ]，求位移增量 [∆δ n ] 2.计算单元应变增量和应力增量，计算总的位移与应力
非线性方程组的解法
弹塑性单元应力调整计算——一维
非线性方程组的解法
弹塑性单元应力调整计算——二维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元高斯点处应力增量和 应力水平 3. 检查各单元高斯点处应力水平: j-1步弹性, j步弹性 j-1步屈服, j步(弹性)卸载 j-1步弹性, j步屈服 j-1步屈服, j步应力继续增长
其他数值方法
质点－桁架模型 弹簧刚度 轴向刚度等效
k1 = 1 Et 2(1 + ν )
ν k2 = Et 2 1 −ν
切向刚度等效
k2 = 1 Et 2(1 + ν )
取 ν = 1/ 3
开裂前，桁架单元的应力－应变关系服从混凝土 单轴应力应变关系
其他数值方法
两个例子
其他数值方法
两个例子
其他数值方法
非线性方程组的解法
弹塑性单元应力调整计算——一维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元应力增量和应力水平 3. 检查各单元应力水平:j-1步弹性, j步弹性 j-1步弹性, j步屈服 j-1步屈服, j步(弹性)卸载 j-1步屈服, j步应力继续增长 4. 区分不同情况,计算应力增量 注意:屈服区分初次屈服和后继屈服
其他数值方法
基本计算过程 1. 分析对象划分成非连续刚体组合 2. 初始运动状态描述 3. 外力作用下，接触形式定义、接触判断、 本构关系 平衡方程求解 4. 确定新的运动状态 5. 运动变化全过程
其他数值方法
力－位移关系
其他数值方法
力－位移关系 法向接触力 Fn = K nδ n 切向接触力 对于库仑摩擦作用 Fs ≤ c + Fn tan ϕ 粘 聚 应 力 摩 擦 角
非线性方程组的解法
逐步增量法
[K ][δ ] = [P]
[K ] 不是常量
代数方程 曲线
[K ]d [δ ] = d [P]
一阶微分方程 欧拉法，龙格-库塔法 分段线性折线
步骤：荷载分成许多小步，每步中，认为结构线性， 刚度矩阵常量 不同的小步中，刚度矩阵可能不同
非线性方程组的解法
逐步增量法－欧拉折线法
…
[δ n ] = [δ n−1 ] + [∆δ n ]
避免了重新计算刚度的麻烦，但迭代次数增加了，收敛速度慢。
非线性方程组的解法
实际计算时，常常同时采用增量法和迭代法
非线性方程组的解法
收敛准则： 不平衡力收敛准则，位移收敛准则
Pres ≤ α P
∆δ k < α δ k
三种范数：各元素绝对值之和，各元素平方根之和，各元素中绝对值最大者
[δ n ] = [δ n−1 ] + [∆δ n ]
[σ n ] = [σ n−1 ] + [∆σ n ]
3.若不是最后荷载，求与本级荷载相应的单元刚度和总体刚度，转入步骤1。 缺点：误差大
非线性方程组的解法
逐步增量法－平均刚度增量法
[K ] = 1 ([K ] + [K ]) 2
n n −1 n
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 几何非线性 物理位置：平衡方程必须按照变形后的几何位置建立 数学描述：物体的变形一般由位移的一阶微分求得， 当变形很大，高阶微分不能忽略 举例：
非线性方程组的解法
主要有两类： 几何非线性，材料非线性 材料非线性 混凝土，钢筋，粘结滑移，时变效应 一般认为，混凝土开裂即存在材料非线性 其他非线性：边界非线性 常见的是接触问题
[K ][∆δ ] = [∆P ]
n n n
[∆δ n ] = [K n ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：……. 缺点：存储两个刚度，存贮量大
非线性方程组的解法
逐步增量法－中点刚度增量法
施加一半荷载，求出相应的位移（中点位移）
[K n−1 ]∆δ n − 1 = 1 [∆Pn ]
非线性方程组的解法
逐步搜索法 可求出极限荷载，无法求解
步骤：1. 加一级荷载 ∆P ，结果发散，退回 2. 加 1 ∆P ，若结果发散，再次退回
2
P −δ
下降段
3. 加
1 ∆P ，……至收敛，得到极限荷载 4
非线性方程组的解法
虚加刚性弹簧法 加虚拟弹簧, 改变原结构形式，新的结构形式 曲线在一定范围内没有下降段
[B] = [M ][R][Q]
其他数值方法
刚体弹簧元
k n [D] = 0
Tresca屈服条件
0 ks
本构矩阵 (2x2)
平面应变下的莫尔-库仑屈服条件
平面应力下的莫尔-库仑屈服条件
其他数值方法
刚体弹簧元
其他数值方法
无网格法
特点 在计算域上，用离散的点拟合场函数，采用基于节点（非单元节点） 的位移插值函数 部分或彻底地消除网格，避免了有限单元法中单元网格的限制 在处理大变形、开裂等问题时，有其优势 采用离散的点拟合场函数的方法 核（权）函数法 移动最小二乘法 径向基函数法 然后建立方程求解
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ， 即时切线刚度 [K 1 ]
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K1 ]−1 [∆P1 ]
[δ 2 ] = [δ1 ] + [∆δ 2 ]
[δ n ] = [δ n−1 ] + [∆δ n ]
节点荷载 [P2 ] ， 即时切线刚度 [K 2 ]
非线性方程组的解法
添加行列法
非线性方程组的解法
强制迭代法 结构总刚是负定的,只要施加负荷载,可求得相应的位移增量
问题:适用范围?
非线性方程组的解法
硬化刚度法 人为增大结构刚度,改善收敛性 一般结构切线刚度可分为
KT = K E − K p
弹性 非线性影响 刚度 刚度
K T = K E − λK p
能量准则
([P ] [δ ]) res
T
两曲线的接近如何描述？
非线性方程组的解法
两个问题： 考虑结构负刚度的一些算法 弹塑性单元应力调整计算
非线性方程组的解法
考虑结构负刚度的一些算法 问题：混凝土结构 P − δ 下降段，由力求位移求解困难 原因：上升段，刚度矩阵正定，线性方程组有解 下降段，刚度矩阵非正定，求解困难 负刚度：非正定的刚度矩阵 几种方法：逐步搜索法，虚加刚性弹簧法，位移加载法， 强制迭代法，硬化刚度法，弧长法
刚体弹簧元 刚体＋弹簧 自由度少，总刚阶数少 可在在极限分析领域运用 思路 弹簧单元刚度矩阵 界面应力应变关系 界面法向和切向的应力应变关系 界面法向和切向刚度
其他数值方法
刚体弹簧元
[K ] = ∫ [B ] [D ][B ]ds
T
k n [D] = 0
0 ks
单元刚度矩阵 (6x6) 本构矩阵 (2x2) 形函数矩阵 (2x6)
网格大小与数目的选择 块体大小，计算精度
其他数值方法
公共面法 块体间接触判断 公共面的定义（如何得到） 假象面 令两块体向薄板靠近 两块体在接触时停止运动 薄板被固定在一个位置，该面即为公共面 块体与公共面间的接触关系
其他数值方法
公共面法 接触类型
其他数值方法
颗粒离散元法 多边形、多面体 圆形、球形 可以以链杆或弹簧连接颗粒 质点－桁架模型
2
2
δ 1 = [δ n −1 ] + ∆δ 1 n− 2 n− 2
K n− 1 确定相应的刚度（中点刚度） 2
以该刚度作为初始刚度求当前位移增量
K n − 1 [∆δ n ] = [∆Pn ] 2
优点: 比两步欧拉折线法精度高
P −δ
问题:在何处加弹簧? 弹簧刚度如何?
非线性方程组的解法
位移控制法 支座位移控制加载,可求支座的反力, 和全部 P − δ 曲线 两个算法:刚度矩阵重新排列法, 添加行列法
非线性方程组的解法
添加行列法
[p1 p2]T为参考矢量, ∆λ 为荷载控制步长系数, [R1 R2]T为不平衡力 问题:外力不按比例加载?
其他数值方法
离散单元法 刚体弹簧元法 无网格法
其他数值方法
离散单元法 Cundall于1971年提出 基本原理:牛顿第二运动定律 建立每个刚体的力平衡方程, 采用动力松弛法求解 特点 岩体或颗粒组合体被模拟成通过边或角相互接触而产 生相互作用 块体之间边界的相互作用可以体现其不连续性和节理 的特征 使用显式积分迭代法,允许大位移、转动和非线性材料
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 割线迭代法
[δ 1 ] = [K 0 ]−1 [P]
[δ 2 ] = [K 1 ]−1 [P ] …… [δ n ] = [K n−1 ]−1 [P]
割线刚度 [K 1 ] 割线刚度 [K 2 ]
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 切线迭代法
有
δFs = K s ∆δ s
t + ∆t
Fs = t F + δFs ≤ c + t Fn tan ϕ
文献中给出了Ks和Kn的建议值
其他数值方法
运动方程 位移和转角、一阶导（速度）、二阶导（加速度） 采用中心差分法，逐步迭代，至达到平衡状态
考虑阻尼（Rayleigh阻尼）
阻尼系数的选择一般依赖经验方法
…
非线性方程组的解法
迭代法－割线迭代法，切线迭代法，等刚度迭代法 等刚度迭代法
[δ 1 ] = [K 0 ]−1 [P]
节点荷载 [P1 ] ，
[∆P1 ] = [P] − [P1 ]
[∆δ 2 ] = [K 0 ]−1 [∆P1 ]
[δ 2 ] = [δ 1 ] +其他数值方法
无网格法
无网格伽辽金法
由移动最小二乘法得到近似场函数 由变分原理得到基本方程
取
问题:适用范围?
刚度硬化 系数 0 ≤ λ ≤ 1
非线性方程组的解法
弧长法
[K ][∆u ] = ∆λ [P] + [R ]
i i +1 i +1 i
荷载控制 步长系数 同时控制 ∆λ 和 ∆u ,使得
不平衡力
[∆u ]T [∆u ] + ∆λ2 = dS
i +1
在 λ − u 空间中,向量r迭代按一圆弧进行,故称弧长法 ……可以求得 ∆λ
其他数值方法
时步选择 具有质量m和刚度系数k的弹性系统，其运动方程
时间步长需满足
∆t ≤ 2 m / k
块体系统由多个质量不一、刚度不同的块体组成， 理论上取其最小的 ∆tmin，实际计算时取小于 ∆tmin
其他数值方法
接触判断算法 二维相对成熟，三维有待发展 分格接触检索算法 1.将研究区域划分成网格 2.在每个块体内，按照接触算法，检索块体是否接触 3.每次迭代，均需接触检索
4. 区分不同情况,计算应力增量
非线性方程组的解法
弹塑性单元应力调整计算——二维
其他数值方法
其他数值方法
有限单元法的局限 基于连续体力学,相邻边界上位移协调, 对于处理应力或位移间断的问题麻烦. 目前解决办法: 1.重新划分网格 计算工作量大 2.改变单元位移插值函数使单元各向异性 网络依赖性,精度差 3.生死单元 网格太细密
[K n−1 ][∆δ n ] = [∆Pn ]
[∆δ n ] = [K n−1 ]−1 [∆Pn ]
[δ n ] = [δ n−1 ] + [∆δ n ]
步骤：1.施加第n步荷载增量 [∆Pn ]，利用始点刚度矩阵 [K n−1 ]，求位移增量 [∆δ n ] 2.计算单元应变增量和应力增量，计算总的位移与应力
非线性方程组的解法
弹塑性单元应力调整计算——一维
非线性方程组的解法
弹塑性单元应力调整计算——二维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元高斯点处应力增量和 应力水平 3. 检查各单元高斯点处应力水平: j-1步弹性, j步弹性 j-1步屈服, j步(弹性)卸载 j-1步弹性, j步屈服 j-1步屈服, j步应力继续增长
其他数值方法
质点－桁架模型 弹簧刚度 轴向刚度等效
k1 = 1 Et 2(1 + ν )
ν k2 = Et 2 1 −ν
切向刚度等效
k2 = 1 Et 2(1 + ν )
取 ν = 1/ 3
开裂前，桁架单元的应力－应变关系服从混凝土 单轴应力应变关系
其他数值方法
两个例子
其他数值方法
两个例子
其他数值方法
非线性方程组的解法
弹塑性单元应力调整计算——一维 1. 已知j-1步迭代结束后的所有量 2. 按弹性材料试算j步各单元应力增量和应力水平 3. 检查各单元应力水平:j-1步弹性, j步弹性 j-1步弹性, j步屈服 j-1步屈服, j步(弹性)卸载 j-1步屈服, j步应力继续增长 4. 区分不同情况,计算应力增量 注意:屈服区分初次屈服和后继屈服
其他数值方法
基本计算过程 1. 分析对象划分成非连续刚体组合 2. 初始运动状态描述 3. 外力作用下，接触形式定义、接触判断、 本构关系 平衡方程求解 4. 确定新的运动状态 5. 运动变化全过程
其他数值方法
力－位移关系
其他数值方法
力－位移关系 法向接触力 Fn = K nδ n 切向接触力 对于库仑摩擦作用 Fs ≤ c + Fn tan ϕ 粘 聚 应 力 摩 擦 角
非线性方程组的解法
逐步增量法
[K ][δ ] = [P]
[K ] 不是常量
代数方程 曲线
[K ]d [δ ] = d [P]
一阶微分方程 欧拉法，龙格-库塔法 分段线性折线
步骤：荷载分成许多小步，每步中，认为结构线性， 刚度矩阵常量 不同的小步中，刚度矩阵可能不同
非线性方程组的解法
逐步增量法－欧拉折线法
…
[δ n ] = [δ n−1 ] + [∆δ n ]
避免了重新计算刚度的麻烦，但迭代次数增加了，收敛速度慢。
非线性方程组的解法
实际计算时，常常同时采用增量法和迭代法
非线性方程组的解法
收敛准则： 不平衡力收敛准则，位移收敛准则
Pres ≤ α P
∆δ k < α δ k
三种范数：各元素绝对值之和，各元素平方根之和，各元素中绝对值最大者
[δ n ] = [δ n−1 ] + [∆δ n ]
[σ n ] = [σ n−1 ] + [∆σ n ]
3.若不是最后荷载，求与本级荷载相应的单元刚度和总体刚度，转入步骤1。 缺点：误差大
非线性方程组的解法
逐步增量法－平均刚度增量法
[K ] = 1 ([K ] + [K ]) 2
n n −1 n