kelvin弹性梁的matlab

伤痕累累

最近用matlab编写kelvin弹性梁仿真程序。里面总是存在小bug，计算到一定时间段就会发生跳跃。不知道有没有做过这方面的，求指导！

凌绝顶

什么是kelvin弹性梁？LZ能否普及一下知识。

伤痕累累

回复 2 # 凌绝顶 的帖子

我不知道怎么样上传图片。你可以百度一下。kelvin粘弹性梁。

伤痕累累

回复 2 # 凌绝顶 的帖子

粘弹性梁

欧阳中华

.
   可以展开介绍一下粘弹性梁和一般梁的区别...

伤痕累累

我给我的差分程序弄上来给大家看看。主要是不会上传文件个图片。
function nextsolv3
%%差分的转化没有问题，但是为什么会出这样的结果呢，欢迎高手给予指导！！！
%我的邮箱：caitylnbug@126.com
%论文里给的时间和空间的步长分别是 dt=1e-6,h=1e-3,做出了倍周期和混沌，运行的时间是700s
clear all
clc
%kelvin弹性梁数学模型，由Hamilton原理推导，经无量纲后方程如下（D表示微分）：
%D^2v(x,t)/Dt^2+2*r*D^2v(x,t)/DtDx+(r^2-1)*D^2v(x,t)/Dx^2+...
%vf^2**D^4v(x,t)/Dx^4+a*D^5v(x,t)/Dx^4*Dt=3/2*k1^2*(Dv(x,t)/Dx)^2*D^2v(x,t)/Dx^2
%初始条件v(0,t)=v(1,t)=0,D^2v(0,t)/Dx^2=D^2v(1,t)/Dx^2=0
%边界条件v(x,0)=K*x*(1-x);Dv(x,0)/Dt=0;K=0.001
%参数下面有给出
reltol=1e-7;
dt=0.01;h=0.1;%时间步长dt,空间步长h
T=0.2;L=1;
x=0:h:L;
y=0:dt:T;%时间点
n=length(x)%空间离散结点个数
m=length(y)%时间离散节点个数
a=0.001;k1=2000;r0=3;w=3.5;r1=0.4;vf=0.8;
D=0.001;
v=zeros(n+4,m+3);%初始化
for j=3:n+2
    for i=3:m+1
        v(n+2,i)=0;%右边界
        v(3,i)=0;%左边界
        r(i)=r0+r1*sin(w*y(i-1));%轴向梁运动速度
        dr(i)=r1*w*cos(w*y(i-1));%轴向速度的导数
        v(j,2)=D*x(j-2)*(1-x(j-2));%下边界
        v(n+3,i)=2*v(n+2,i)-v(n+1,i);%边界有限差分的处理
        v(2,i)=2*v(3,i)-v(4,i);%边界有限差分的处理
        c0=2/dt^2+2*(r(i)^2-1)/h^2-6*vf^2/h^4;c1=1/dt^2;c2=2*r(i)/(4*h*dt);c3=dr(i)/(2*h);c4=(r(i)^2-1)/h^2;c5=vf^2/h^4;c6=a/(2*dt*h^4);c7=1.5*k1^2/(4*h^4);
        v(j,i)=1/c0*(c1*(v(j,i+1)+v(j,i-1))+c2*(v(j+1,i+1)-v(j-1,i+1)-v(j+1,i-1)+v(j-1,i-1))+c3*(v(j+1,i)-v(j-1,i))+c4*(v(j+1,i)+v(j-1,i))+c5*(v(j-2,i)-4*v(j-1,i)-4*v(j+1,i)+v(j+2,i))+...
            c6*(v(j+2,i+1)-4*v(j+1,i+1)+6*v(j,i+1)-4*v(j-1,i+1)+v(j-2,i+1)-v(j+2,i-1)+4*v(j+1,i-1)-6*v(j,i-1)+4*v(j-1,i-1)-v(j-2,i-1))-c7*(v(j+1,i)-v(j-1,i))^2*(v(j+1,i)-2*v(j,i)+v(j-1,i)));
    end
end
V=v(3:n+2,2:m+1)
plot(y,V((n+1)/2,:),'b-')%v(0.5,t)点的振幅
figure
plot(y,V)end
V=v(3:n+2,2:m+1)
plot(y,V((n+1)/2,:),'b-')%v(0.5,t)点的振幅
figure
plot(y,V)

凌绝顶

回复 6 # 伤痕累累 的帖子

我根据你上面提供的代码将你的模型敲进Mathtype里，然后截图上传如下：

不知道我输入的对不对。我有个问题，就是你的程序中为什么会出现r0，r1和w呢？另外，请将你自己遇到的问题讲清楚一点。

凌绝顶

检查一下你的数值方法是不是有奇异点

伤痕累累

回复 7 # 凌绝顶 的帖子

首先要谢谢你。
r=r0+r1*sin（w*t），这是变化的轴向速度。
dr是r的一阶导数。
现在就是里面的程序问题。经过无量转化之后，matlab里，取不了原论文里那么大的步长，否则
v=zeros(n+4,m+3);%初始化     会显示超出矩阵的最大维数（当然我也借用过矩阵维数令v=sparse（1e-8,1e-8）,还是不行）。我后来想，作者应该不是用matlab做的程序，应该是VC或者FROTRAN做的。

伤痕累累

回复 8 # 凌绝顶 的帖子

奇异点在计算的过程中应该怎么检查呢？

凌绝顶

回复 9 # 伤痕累累 的帖子

应该不至于超出矩阵的最大维数的限制吧。你按照文献中离散得到的矩阵有多大啊？是不是出现了Index exceeds matrix dimensions的提示？？这并不是超过了matlab允许的最大的矩阵维数的意思，你应该检查一下你的矩阵操作过程是不是有问题。

凌绝顶

回复 10 # 伤痕累累 的帖子

就是看看计算过程中会不会出现在数学上没有意义的点，比如分母为零啊等。

凌绝顶

动力学方程中vf、a、k1以及初始条件中的k，它们有什么物理意义？

伤痕累累

回复 13 # 凌绝顶 的帖子

给你的邮箱留下我发给你看看。站内信

凌绝顶

回复 14 # 伤痕累累 的帖子

我的邮箱见站内信