newmark
在结构动力分析的时候刚度矩阵,阻尼矩阵,质量矩阵以及荷载矩阵已经建立好了,用NEWMARK做时间历程的分析.怎么能引入边界条件,比如简支梁,两端的位移等于0,是不是将刚度矩阵相应的主元置1,相应的行列元素为0就可以了,其他的矩阵怎么变化 .简单的方法是将刚度矩阵对应约束自由度的对角线元素的值置个比较大的数.. .. 也可以将M、K、C相应的行列划去 刚性支撑没有划线的好，还是划线吧 边界条件应该在建立刚度矩阵的时候就确定。 离散的时候就是根据边界条件来取形函数的阿？所以已经考虑边界条件了，你说的是初始条件吗？/ 一般认为，边界条件是指支承结构体系的支座，它直接影响结构体系的振动特性。
而初始条件是指激振开始时的状态，是起始时间的条件。
所以，边界条件和初始条件不是一个概念。 下面的是我做的结构动力学的作业，用matlab编写的。我没有遇到楼主说的边界条件问题。是否说的是开始计算时候的初始条件？
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%Newmark β method
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%结构参数
m=*1.0e+3; %质量
k=*1.0e+5; %刚度
cn=length(m); %层数
I=diag(ones(cn)); %单位向量
ksi1=0.05; %对应于第一振型的阻尼比
ksi2=0.07; %对应于第二振型的阻尼比
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%地震参数
dzb=dzhbo; %地震波
xs=200/max(abs(dzb)); %地震地面加速度放大系数xs=A'max/Amax
ag=dzb*0.01*xs; %地面运动加速度
dt=0.02; %步长
ndzh=400; %地震时步的步数
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%m,k
M=diag(m); %生成广义质量阵
K=zeros(cn); %生成整体刚度阵
for i=1:cn-1
K(i,i)=k(i)+k(i+1); %主对角元素
K(i,i+1)=-k(i+1); %次对角元素
K(i+1,i)=-k(i+1) ; %次对角元素
end
K(cn,cn)=k(cn); %最后元素
=eig(K,M); %返回矩阵K,M的广义特征值和广义特征向量
d=diag(sqrt(D)); %求解自振圆频率
for i=1:cn
=min(d); %对w 进行排序
xgd(:,i)=Y(:,j);
d(j)=max(d)+1;
end
w=dl; %频率向量
T=2*pi./w; %周期向量
Y=xgd; %振型矩阵
alfa1=2*w(1)*w(2)*(ksi1*w(2)-ksi2*w(1))/(w(2)^2-w(1)^2);%瑞利阻尼系数alfa1
alfa2=2*(ksi2*w(2)-ksi1*w(1))/(w(2)^2-w(1)^2); %瑞利阻尼系数alfa2
C=alfa1*M+alfa2*K; %阻尼矩阵
for j=1:cn
Y(:,j)=Y(:,j)/Y(cn,j);
gama(j)=(Y(:,j))'*M*I/((Y(:,j))'*M*Y(:,j)); %j振型的振型参与系数
end
for i=0:ndzh
P(:,i+1)=-M*I*ag(i+1); %地震作用
end
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%γβ参数的选取和系数的计算
gama=1/2;beta=1/4; %gama和beta选取
c0=1/(beta*dt^2);
c1=gama/(beta*dt);
c2=1/(beta*dt);
c3=1/(2*beta)-1;
c4=gama/beta-1;
c5=dt*(gama/beta-2)/2;
c6=dt*(1-gama);
c7=gama*dt;
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%地震作用计算
pk=K+c0*M+c1*C; %形成有效刚度矩阵
pk=inv(pk); %有效刚度矩阵逆阵
u0=zeros(cn,1); %初始位移
du0=zeros(cn,1); %初始速度
ddu0=zeros(cn,1); %初始加速度
for i=0:ndzh
ff=P(:,i+1)+M*(c0*u0+c2*du0+c3*ddu0)+C*(c1*u0+c4*du0+c5*ddu0);
u1=pk*ff;
ddu1=c0*(u1-u0)-c2*du0-c3*ddu0;
du1=du0+c6*ddu0+c7*ddu1;
u0=u1;
ddu0=ddu1;
du0=du1;
A1(:,i+1)=u1;
dA1(:,i+1)=du1;
ddA1(:,i+1)=ddu1;
end
A1; %速度阵
dA1; %位移阵
ddA1; %加速度阵
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%求层间剪力
F=M*ddA1+P; %地震作用力
V=flipud(cumsum(flipud(F))); %结构层间剪力
%%%%%%%%%%%%%%%%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&%%%%%%%%%%%%%%%%%%%%%%%%
%绘图
t=0:dt:dt*ndzh;
for i=1:cn
subplot(cn,1,i);
plot(t,A1(i,:));
grid on
titlen=['第',int2str(i),'层位移'];
title(titlen);
xlabel('时间t');
ylabel('位移');
end
figure
for i=1:cn
subplot(cn,1,i);
plot(t,dA1(i,:));
grid on
titlen=['第',int2str(i),'层速度'];
title(titlen);
xlabel('时间t');
ylabel('速度a');
end
figure
for i=1:cn
subplot(cn,1,i);
plot(t,ddA1(i,:));
grid on
titlen=['第',int2str(i),'层加速度'];
title(titlen);
xlabel('时间t');
ylabel('加速度a');
end
figure
for i=1:cn
subplot(cn,1,i);
plot(t,V(i,:));
grid on
titlen=['第',int2str(i),'层间剪力'];
title(titlen);
xlabel('时间t');
ylabel('层间剪力N');
end 楼主,你的dzhbo是如何输入 刚性支撑下还是划线法比较好
页:
[1]