《有限元仿真实践原理》1D网格划分 (1)
本文摘要（由AI生成）：
本文主要探讨了刚度矩阵在结构分析和有限元分析中的核心概念、重要性和应用。刚度矩阵描述了结构在受力时的位移响应，是结构分析的关键属性。文章介绍了刚度矩阵的推导方法，包括直接法、变分法和加权残值法，并强调其在预测结构行为中的关键作用。通过Rod单元和Quad4单元的例子，文章详细说明了刚度矩阵的推导和应用。此外，还介绍了有限元软件利用刚度矩阵求解结构方程的过程，并提供了相关练习以加深对刚度矩阵和有限元分析的理解。本文为工程师和研究人员提供了有效进行结构分析和设计的指导。
刚度矩阵推导
什么是刚度以及为什么我们在FEA中需要它？
刚度‘K’定义为力/长度（单位N/mm）。物理意义 – 刚度等于产生单位位移所需要的力。刚度取决于几何形状以及材料属性。
考虑3个几何尺寸完全相同的二力杆 – 铸铁、低碳钢和铝。如果我们测量产生1mm位移所需要的力，铸铁需要的力最大，然后依次是低碳钢和铝，即KCI> KMS > KAl。
现在考虑3个相同材料不同截面的二力杆。同样，产生单位位移所需要的力是不同的。所以，刚度不仅依赖于材料，也依赖于几何形状。
刚度矩阵的重要性 - 对于结构分析，刚度是一个非常重要的属性。
线性静态分析的方程是[F] = [K] [D]。力通常是已知的，位移是未知的，而刚度是单元的特有属性。这就意味着如果我们用公式表达一个给定形状的刚度矩阵，比如线、四边形或者四面体，那我们就可以通过网格划分来表达任何几何形状并使用方程F = K D求解。
公式表达刚度矩阵的方法：
直接法
变分法
加权残值法
直接法很容易理解，但是很难用电脑程序表达。而变分法和加权残值法很难理解，但是从编程的角度来说很简单。这就是为什么所有的软件要么使用变分法，要么使用加权残值法。
直接法推导二力杆单元的刚度矩阵
直接法推导刚度矩阵的方法：
对于一个给定的单元，假设有n个自由度（比如，一个quad4单元的所有自由度 = 4*6 = 24）。
步骤1）假设第一个自由度 ≠ 0，并且其它所有自由度 = 0。生成方程1。
步骤2）假设第二个自由度 ≠ 0，并且其它所有自由度 = 0。生成方程2。
……
步骤n）假设第n个自由度 ≠ 0，并且其它所有自由度 = 0。生成方程n。
步骤n+1）所有方程相加，1 + 2 + 3 + 4 ………..+ n。
步骤n+1会给我们一个非常广义的刚度矩阵方程。
Rod单元
刚度矩阵的属性
刚度矩阵的阶数与总自由度有关。
一个奇异刚度矩阵表示结构未被约束并有刚体运动。
刚度矩阵的各列是一个节点力的平衡集 合，节点力产生对应自由度的位移。
一个对称的刚度矩阵表示力与位移成正比。
矩阵的对角线项全是正值表示一个力指向左边就不能产生向右的位移。只有结构在不稳定的情况下，对角线项才会是零或者负值。
rod单元只支持拉伸或压缩，不支持剪力、弯矩和扭矩。在上述方程中，刚度矩阵的阶数是2x2，未知量的数量是2。
未知量的数量 = 自由度的数量 – 被“单点约束”约束了的自由度的数量（对于固定的节点，自由度记为0）
通常与模型中的总自由度相比，约束的数量是可以忽略不计的，所以未知量的数量近似等于总自由度。
刚度矩阵的阶数 = 总自由度 x 总自由度
Quad4单元
自由度/节点 = 6
总自由度 = 6X4 = 24,
刚度矩阵阶数 = 24 x 24
|F|24x1 = |K|24x24 |δ|24x1
这就意味着求解一个quad4单元的问题，软件内部将会求解24个方程。
Tetra10 单元
自由度/节点 = 3
总自由度 = 3*10 = 30
刚度矩阵阶数 = 30 x 30
|F|30x1 = |K|30x30 |δ|30x1
对于一个给定的有限元模型，软件需要求解多少个方程呢？
假设模型中只有薄壳单元，共20,000节点（6 自由度/节点）。
总自由度 = 20000*6 =120,000
刚度矩阵阶数 =120,000x120,000
有限元软件会在内部求解的方程数 = 120,000
固体力学中变分法有限元与加权残值法有限元的等效表格
问题：在一个1D的rod单元的一端（自由端）施加一个集中拉力，另一端固定。
对应的微分方程和边界条件为：
AE( d2u/dx2) = 0 u |x =0 =0
AE ( du/dx) |x=L = P
练习1：对Rod单元施加一个拉力
E= 2.1 * e5 N/mm2, u = 0.30（钢）
解析结果：
应力 = F/A = 100,000/(3.14*25*25) = 50.95 N/mm2
位移 = FL / AE = 100,000*500 / (3.14*25*25*2.1E05) = 0.12 mm
怎么使用软件求解：下面是所有商业软件要求的基本步骤。但是对于不同的软件，步骤的顺序和输入文件的格式可能会有所不同。
步骤1 梁截面
定义圆形截面（直径50 mm）
步骤2 划分网格
创建节点
手工输入节点坐标（0，0，0），（500，0，0）
点击上述2个节点创建Rod单元（指定圆形截面）
创建材料并赋予给这个单元
步骤3 边界条件
指定分析类型（线性静态分析）
沿X方向施加集中力 = 100,000 N
施加约束 = 约束x位移
步骤4 求解
创建求解集/工况
求解
步骤5 后处理
显示位移和应力结果
播放变形动画
比较解析结果和软件求解结果
软件内部是怎么求解这个问题的？
uj = 0.12 mm ε = uj / L = 0.12/500 =0.00024638
σ = E *ε = 2.068 E05 * 0.00024638 = 50.9N/mm2
练习2：Rod单元不能承载剪力
所有步骤都与练习1一样，除了把力施加在y轴方向。比较有限元结果和手工计算结果。有限元结果是位移很大，应力为零。这是因为rod单元没有能够承载垂直力的自由度（即保持系统平衡）。
为了承载垂直力（ΣFy = 0 且 ΣMz = 0），我们需要一个支持Y轴移动和Z轴转动的单元。
Beam和Bar单元支持这些自由度，如果你使用一个beam单元或者bar单元重做一次这个练习，你将会得到正确的结果。
实际情况中，几乎不会只使用一个单元做分析。通常是很多单元组合使用。
当模型中包含了很多单元时，软件会把各个单元的刚度矩阵组合到一起。组合的刚度矩阵的阶数等于模型的总自由度。
