Duhamel积分的基本原理是将任意载荷视为一系列脉冲激励的叠加,系统的响应可以表示成对这些脉冲激励的响应函数的叠加。在单自由度系统中,Duhamel积分的形式为:
其中, 是系统在时间 的位移响应, 是外载荷函数, 是质量, 是无阻尼自然频率, 为有阻尼时结构的频率, 是阻尼比.
当外荷载 不能用解析函数表达时,则式(1)需按等时间间隔 进行离散,用数值方法来计算:
如图1所示,在 积分时,每一个小矩形的高为 。注意这不是斐波那契数列那样的简单相加,因为函数中的 随时在变化。因此,需要定义一个高阶函数,比如
def f1(n):
def f2(x):
return x * n
return f2
times3 = f1(3)
print(times3(10)) # 输出:30, 10的3倍
式(2)的高阶函数
# 杜哈梅积分函数,荷载为简谐荷载sin(t)
def func_duhame(zeta,omg,t):
def func1(tau):
y = np.sin(tau) * np.exp( -zeta*omg*(t-tau) ) * np.sin(omg*(t- tau ) )
return 1/ omg * y
return func1
计算图1中 时刻的动位移,可以用 这两个高阶函数。
f = func_duhame(zeta,omg,t) # 返回已经代入t的式(1)
l = np.arange(0.0, t, delta_tau)
q1 = map(f,l) # 计算每个小矩形的高
q2 =map(lambda x: x*delta_tau, q1) # 计算每个小矩形的面积
res = reduce(lambda x,y: x+y,q2) # 每个小矩形的面积累加得到积分结果
完整代码
import matplotlib.pyplot as plt
from functools import reduce
import numpy as np
# 杜哈梅积分函数
deffunc_duhame(zeta,omg,t):
deffunc1(tau):
y = np.sin(tau) * np.exp( -zeta*omg*(t-tau) ) * np.sin(omg*(t- tau ) )
return1/ omg * y
return func1
if __name__ == "__main__":
plt.style.use("ggplot")
plt.rcParams['font.sans-serif'] = ['MiSans'] # 这个字体能够支持中英文字符
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
omg = 1.0 # 自振频率
zeta = 0.05 # 固有阻尼比
dt = 0.005 # 时间间隔
t0 = np.arange(0.0, 40, dt)
step = len(t0)
dpms = np.zeros( step )
for i,t inenumerate(t0):
f = func_duhame(zeta,omg,t)
if i==0:
l = [0, t0[0]]
else:
l = t0[0:i+1]
res = reduce(lambda x,y: x+y, map(lambda x: x*dt, map(f,l)) )
dpms[i] = res
fig, ax1 = plt.subplots(nrows=1, ncols=1, sharex=True, )
ax1.plot(t0, dpms, linestyle='--', color='#1661ab',label='数值法')
ax1.set( title= '位移时程')
ax1.grid()
ax1.legend(frameon=True, # 启用边框
framealpha=0.8, # 设置透明度
fancybox=False, # 边框圆角与否
shadow=True, # 是否启用阴影
facecolor='lightgrey', # 背景颜色
edgecolor='dimgrey', # 边框颜色
loc= 'upper left')
fig.savefig("32.png",dpi=300)
plt.show()