專業方向系列-00-Python與有限元初探

既生喻何生亮發表於2018-10-22
Python

案例1

給出4個彈簧的勁度係數，離散後，求其總的剛度矩陣。
程式碼：

import numpy as np
k1, k2, k3, k4 = 500, 250, 2000, 1000
ki = np.array([k1, k2, k3, k4])
p = 1000
K = np.matrix(np.zeros((4, 4)))
eindex = [K[0:2, 0:2], K[1:3, 1:3], K[0:3:2, 0:3:2], K[2:, 2:]]
for i in range(4):
    ke = ki[i] * np.array([[1, -1], [-1, 1]])
    eindex[i] += np.matrix(ke)   
print(K)

結果如下：
[[ 2500. -500. -2000. 0.]
[ -500. 750. -250. 0.]
[-2000. -250. 3250. -1000.]
[ 0. 0. -1000. 1000.]]

案例2

兩根杆的例項。
程式碼分析：

import numpy as np
# 系統的初始引數如下：
A1, A2, L1, L2, E1, E2 = 2400, 600, 300, 400, 7e4, 2e5
# 計算每根杆的剛度
k1 = (E1*A1/L1)*np.array([[1, -1], [-1, 1]])
print(`k1`, k1)
k2 = (E2*A2/L2)*np.array([[1, -1], [-1, 1]])
print(`k2`, k2)
ki = [k1, k2]
# 初始化總的剛度矩陣、力以及位移
K, F, d = np.zeros((3, 3)), np.zeros((3, 1)), np.zeros((3, 1))
print(K)
print(F)
print(d)
# 組裝總的剛度矩陣
eindex = [K[0:2, 0:2], K[1:, 1:]]
for i in range(2):
    eindex[i] += ki[i]
print(K)
# 給出力向量
F[1] = 200000
print(F)
# 求解過程
KK = K[1, 1]
FF = F[1]
dd = FF/KK
d[1] = dd
print(d)

計算結果：

  1. 每根杆的剛度矩陣
    k1 [[ 560000. -560000.]
    [-560000. 560000.]]
    k2 [[ 300000. -300000.]
    [-300000. 300000.]]
  2. 總的剛度矩陣
    [[ 560000. -560000. 0.]
    [-560000. 860000. -300000.]
    [ 0. -300000. 300000.]]
  3. 力向量
    [[ 0.]
    [200000.]
    [ 0.]]
  4. 求解的位移向量
    [[0. ]
    [0.23255814]
    [0. ]]

相關文章