序：我要寫一期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年以后了。