第一篇：弹簧单元
序：我要写一期python和Abaqus与有限元的文章，从弹簧单元、杆单元一直到实体单元，通过简单的实例用python编程，Abaqus验证结果。
例：弹簧串联受外力作用，具体数值如下图所示，求：a）总刚；b）节点2与节点3的位移；c）节点1的反力；d）弹簧内力。
一、有限元法求解
步骤1：离散化
单元 |
节点i |
节点j |
1 |
1 |
2 |
2 |
2 |
3 |
步骤2：写单刚
步骤3：写总刚
步骤4：边界条件
本例中，u1=0，F2=0，F3=1000N，代入上述方程
步骤5：求方程，解u2和u3
利用上述方程不难解出u2=10m，u3=15m，具体不再赘述。
步骤6：后处理，求节点1反力F1与弹簧内力f1、f2
取出相应的方程可求得F1=-1000N，f1=1000N（拉），f2=1000N（拉）。
二、python求解
import numpy as np
# 离散化
ele_node = np.array([[1, 2], [2, 3]])
node_coord = np.array([0, 1, 2])
# 计算单刚
dimen = 1 # 1d
k1 = 100
k2 = 200
def fun(K):
Ke = np.array([[K, -K], [-K, K]])
return Ke
# 总自由度数
ndof = dimen * len(node_coord)
# 外力向量
F = np.array([0, 0, 1000])
K = np.array([k1, k2])
# 初始化总刚
K_t = np.zeros((ndof, ndof))
# 计算总刚
for e in range(0, len(ele_node)):
n1 = ele_node[e][0] - 1
n2 = ele_node[e][1] - 1
K_t[np.ix_([n1, n2], [n1, n2])] += fun(K[e])
print(K_t)
# 解u2和u3
k=K_t[1:3, 1:3]
f=F[1:3]
k=np.mat(k)
f=np.mat(f)
u=k.I*f.T
print(u)
# 解支反力F
u=np.mat(u)
U=np.r_[[[0]],u]
U=np.mat(U)
F=K_t*U
print(F)
# 解内力
u1=U[0:2]
u2=U[1:3]
f1=fun(k1)*u1
f2=fun(k2)*u2
print(f1)
print(f2)
计算结果如下图
三、Abaqus求解
步骤1：离散化
步骤2：创建弹簧单元、施加边界条件
步骤3：求解，后处理
结果：
1.整体刚度矩阵如下图：
2.节点2位移u2=10m，节点3位移u3=15m，如下：
3.节点1的支反力F1=-1000N，如下：
4.弹簧单元内力均为1000N，如下：
结论：
不难发现，手工计算、python计算、Abaqus计算结果完全一致。
我看后台有朋友回复我做设计的都开始搞有限元了，哈哈！其实这几年我一直搞专用车力学工作，当时起了个专用车设计名字打算分享一些与专用车设计相关的好文章。我发现写文章比日常工作费力的多，不知下一篇是不是又要2年以后了。
