首页/文章/ 详情

杜哈梅积分的数值方法

6月前浏览721

 

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.040, 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()

 


来源:数值分析与有限元编程
UM
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2025-05-10
最近编辑:6月前
太白金星
本科 慢慢来
获赞 12粉丝 28文章 367课程 0
点赞
收藏
作者推荐

挠度曲线微分方程的应用

力学概念|推导挠度曲线的微分方程有限元计算只能得到梁单元两端的节点位移和内力,单元任意截面的位移和内力需要通过插值得到。梁的挠曲线近似微分方程为弯矩、剪力和载荷集度之间的微分关系式(1)对求导并结合(2)(3)得由式(5)可以看出,如果在梁的某一区间上载荷集度为0,则挠曲线方程为关于的三次多项式;如果载荷集度为非零常数,则挠曲线方程为关于的四次多项式;如果载荷集度为关于的一次多项式,则挠曲线方程为关于的五次多项式。依此类推。而构造梁单元只使用了三次Hermite多项式,因此有限元计算的结果是近似值。为了使得求解过程一般化,将区间映射至标准区间,映射关系为其中,记由链式求导法则则结合(1)(4)(5)有q为0的情况梁在区间不受分布力作用,即。将区间映射到区间,如图1所示。在区间,挠度是关于的三次多项式,在区间,挠度是关于的三次多项式。记区间两个节点挠度和转角分别为和,和,注意到在节点1和2处的值分别是和,利用挠度和转角间的微分关系,构造挠度的Hermite插值多项式为其中令,便可由(8)得到中点的挠度,弯矩和剪力。q为常数的情况在区间,挠度是关于的四次多项式,在区间,挠度是关于的四次多项式。记区间两个节点挠度和转角分别为和,和,以及中间点的挠度,注意到在节点1和2处的值分别是和,利用挠度和转角间的微分关系,构造挠度的Hermite插值多项式为其中q为线性分布的情况在区间,挠度是关于的五次多项式,在区间,挠度是关于的五次多项式。记区间两个节点挠度和转角分别为和,和,以及中间点的挠度和转角,注意到在节点1、2和3处的值分别是和,利用挠度和转角间的微分关系,构造挠度的Hermite插值多项式为其中举例以长度为1,端部集中力为1的悬臂梁为例,节点位移应该是,将梁划分为25个小段,画出弯矩图,剪力图和变形图。Python代码importmatplotlib.pyplotaspltimportnumpyasnpdefmain():plt.style.use("ggplot")#绘图风格plt.rcParams['font.sans-serif']=['MiSans']#这个字体能够支持中英文字符plt.rcParams['axes.unicode_minus']=False#正常显示负号EI=1l=1d=np.array([0,0,-0.333,-0.5])#节点位移xi=np.linspace(-1,1,25)#梁划分为25个小段len=np.size(xi)#插值基函数N1=0.25*(1-xi)**2*(2+xi)N2=0.25*(1-xi)**2*(1+xi)N3=0.25*(1+xi)**2*(2-xi)N4=0.25*(1+xi)**2*(xi-1)#插值基函数的二阶导数N1_2=0.25*6*xiN2_2=0.25*(-2+6*xi)N3_2=-0.25*6*xiN4_2=0.25*(2+6*xi)#插值基函数的三阶导数N1_3=0.25*6*np.ones(len)N2_3=0.25*6*np.ones(len)N3_3=-0.25*6*np.ones(len)N4_3=0.25*6*np.ones(len)v=N1*d[0]+0.5*l*N2*d[1]+N3*d[2]+0.5*l*N4*d[3]v_2=N1_2*d[0]+0.5*l*N2_2*d[1]+N3_2*d[2]+0.5*l*N4_2*d[3]v_3=N1_3*d[0]+0.5*l*N2_3*d[1]+N3_3*d[2]+0.5*l*N4_3*d[3]m=-(2/l)**2*EI*v_2fs=(2/l)**3*EI*v_3x1=np.array(xi[0])x2=np.array(xi[-1])xx=np.hstack((x1,xi,x2))m1=np.array([0])mm=np.hstack((m1,m,m1))fs1=np.array([0])fs2=np.hstack((fs1,fs,fs1))fig,(ax0,ax1,ax2)=plt.subplots(nrows=3,ncols=1,sharex=True,figsize=(12,6))ax0.plot([xi[0],xi[-1]],[0,0])ax0.plot(xx,mm)ax0.set(title='弯矩图')ax0.grid()ax1.plot([xi[0],xi[-1]],[0,0])ax1.plot(xx,fs2)ax1.set(title='剪力图')ax1.grid()ax2.plot([xi[0],xi[-1]],[0,0])ax2.plot(xi,v)ax2.set(title='变形图')ax2.grid()fig.savefig("test.png")plt.show()if__name__=="__main__":main()以长度为1,均布荷载为8的简直梁为例,节点位移应该是,再将跨中线位移-0.104组成插值位移。将梁划分为25个小段,画出弯矩图,剪力图和变形图。importmatplotlib.pyplotaspltimportnumpyasnpdefmain():plt.style.use("ggplot")#绘图风格plt.rcParams['font.sans-serif']=['MiSans']#这个字体能够支持中英文字符plt.rcParams['axes.unicode_minus']=False#正常显示负号EI=1l=1d=np.array([0,-0.333,0,0.333,-0.104])#插值位移xi=np.linspace(-1,1,25)#梁划分为25个小段#插值基函数N1=-0.25*(xi-1)**2*(2*xi**2+3*xi)N2=-0.25*(xi-1)**2*(xi**2+xi)N3=-0.25*(1+xi)**2*(2*xi**2-3*xi)N4=0.25*(1+xi)**2*(xi**2-xi)N5=(xi-1)**2*(1+xi)**2#插值基函数的二阶导数N1_2=-0.25*(24*xi**2-6*xi-8)N2_2=-0.25*(12*xi**2-6*xi-2)N3_2=-0.25*(24*xi**2+6*xi-8)N4_2=0.25*(12*xi**2+6*xi-2)N5_2=12*xi**2-4#插值基函数的三阶导数N1_3=-0.25*(48*xi-6)N2_3=-0.25*(24*xi-6)N3_3=-0.25*(48*xi+6)N4_3=0.25*(24*xi+6)N5_3=24*xiv=N1*d[0]+0.5*l*N2*d[1]+N3*d[2]+0.5*l*N4*d[3]+N5*d[4]v_2=N1_2*d[0]+0.5*l*N2_2*d[1]+N3_2*d[2]+0.5*l*N4_2*d[3]+N5_2*d[4]v_3=N1_3*d[0]+0.5*l*N2_3*d[1]+N3_3*d[2]+0.5*l*N4_3*d[3]+N5_3*d[4]m=-(2/l)**2*EI*v_2fs=(2/l)**3*EI*v_3print(m)x1=np.array(xi[0])x2=np.array(xi[-1])xx=np.hstack((x1,xi,x2))m1=np.array([0])mm=np.hstack((m1,m,m1))fs1=np.array([0])fs2=np.hstack((fs1,fs,fs1))fig,(ax0,ax1,ax2)=plt.subplots(nrows=3,ncols=1,sharex=True,figsize=(12,6))ax0.plot([xi[0],xi[-1]],[0,0])ax0.plot(xx,mm)ax0.set(title='弯矩图')ax0.grid()ax1.plot([xi[0],xi[-1]],[0,0])ax1.plot(xx,fs2)ax1.set(title='剪力图')ax1.grid()ax2.plot([xi[0],xi[-1]],[0,0])ax2.plot(xi,v)ax2.set(title='变形图')ax2.grid()fig.savefig("test1.png",dpi=600)plt.show()if__name__=="__main__":main()来源:数值分析与有限元编程

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈