弹簧单元基础篇

例：弹簧串联受外力作用，具体数值如下图所示，求：a）总刚；b）节点2与节点3的位移；c）节点1的反力；d）弹簧内力。

一、有限元法求解

步骤1：离散化

步骤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计算结果完全一致。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删