最近做結構動力學的編程作業。

大部分人都是用matlab做的，看了下網上的相關代碼也都是matlab。人生苦短，我用python，這次作業我選擇用python來寫（其實是我也不會寫matlab），在這裡給大家分享一下。

這次先開個頭，講一下用於結構時程分析的一種直接積分方法：Newmark-beta法，在Python中的實現。

在本文中就不再闡述Newmark-beta法的具體原理和推導了，具體可見《結構動力學》。

Newmark-beta的求解過程如下：

（1）基本數據準備和初始條件計算：

1）選擇時間步長 、參數 和 ，並計算積分常數：

# 積分步長一般取0.02s
t=0.02
# 選擇參數β和γ（下面的取法可保證該積分無條件穩定），並計算積分常數
gama = 0.5
beta = 0.25
p = [] # 積分常數
p.append(1 / beta / t ** 2)
p.append(gama / beta / t)
p.append(1 / beta / t)
p.append(1 / 2 / beta - 1)
p.append(gama / beta - 1)
p.append(t / 2 * (gama / beta - 2))
p.append(t * (1 - gama))
p.append(gama * t)


2）確定運動的初始值｛u0｝、｛v0｝和｛a0｝（初始位移、速度、加速度）：

# 存儲各時間點的位移、速度以及加速度
u = []
v = []
a = []
# 給u、v、a添加運動初始值，為簡單起見，在本文中默認計算結構為葫蘆串模型，n為質點數（層數）
u.append(np.zeros(n))
v.append(np.zeros(n))
a.append(np.zeros(n))
# 確定每個質點的初始運動狀態
for i in range(n):
u[0] = u0
v[0] = v0
a[0] = a0


（2）形成剛度矩陣K、質量矩陣M、阻尼矩陣C

（3）形成等效剛度矩陣：

# 形成等效剛度矩陣K_
K_ = K + p[0] * M + p[1] * C


（4）計算 時刻的等效荷載：

# 計算ti+1時刻的等效荷載
P_ = P[:, i] + np.dot(M, (p[0] * u[-1] + p[2] * v[-1] + p[3] * a[-1])) + np.dot(C, (
p[1] * u[-1] + p[4] * v[-1] + p[5] * a[-1]))


（5）計算 時刻的位移：

# 計算ti+1時刻的位移
u.append(np.dot(np.linalg.inv(K_), P_))


（6）計算 時刻的加速度和速度：

# 計算ti+1時刻的加速度和速度
a.append(p[0] * (u[-1] - u[-2]) - p[2] * v[-1] - p[3] * a[-1])
v.append(v[-1] + p[6] * a[-2] + p[7] * a[-1])


循環（4）-（6）計算步驟，可以得到線彈性體系在任意時刻的動力反應，對於非線性問題，則應循環（2）-（6）計算步驟完成計算。

將以上六步歸結為一個函數，代碼如下：

"""
函數說明:Newmark-beta積分方法

Parameters:
t - 時間步長
K - 剛度矩陣
M - 質量矩陣
C - 阻尼矩陣
P - 荷載向量
u0 - 初始位移
v0 - 初始速度
a0 - 初始加速度
Returns:
a, v, u - 荷載時長內的結構響應,分別為加速度、速度和位移

Modify:
2018-07-29
"""
def newmark_beta(t, K, M, C, P, u0, v0, a0):
# 確定質點數量
n = len(K)
# 確定總時間步
nt = len(P[0])
# 選擇參數β和γ，並計算積分常數
gama = 0.5
beta = 0.2
p = [] # 積分常數
p.append(1 / beta / t ** 2)
p.append(gama / beta / t)
p.append(1 / beta / t)
p.append(1 / 2 / beta - 1)
p.append(gama / beta - 1)
p.append(t / 2 * (gama / beta - 2))
p.append(t * (1 - gama))
p.append(gama * t)
# 存儲各時間點的位移、速度以及加速度
u = []
v = []
a = []
# 給u、v、a添加運動初始值
u.append(np.zeros(n))
v.append(np.zeros(n))
a.append(np.zeros(n))
# 確定每個質點的初始運動狀態
for i in range(n):
u[0] = u0
v[0] = v0
a[0] = a0
# 形成等效剛度矩陣K_
K_ = K + p[0] * M + p[1] * C
# 遞推求解
for i in range(nt):
# 計算ti+1時刻的等效荷載
P_ = P[:, i] + np.dot(M, (p[0] * u[-1] + p[2] * v[-1] + p[3] * a[-1])) + np.dot(C, (
p[1] * u[-1] + p[4] * v[-1] + p[5] * a[-1]))
# 計算ti+1時刻的位移
u.append(np.dot(np.linalg.inv(K_), P_))
# 計算ti+1時刻的加速度和速度
a.append(p[0] * (u[-1] - u[-2]) - p[2] * v[-1] - p[3] * a[-1])
v.append(v[-1] + p[6] * a[-2] + p[7] * a[-1])
return a, v, u

推薦閱讀：