newmark法计算结构响应

clc
clear
%load dizhengbo4.txt

%function [accveldsp]=transresp1（kkccmmfdbcdofntdtq0dq0）

%ds=xlsread（‘ds.xlsx‘）;
% wave1=dizhengbo4;n    
% wave=wave1‘;
kk=xlsread（‘k_b.xlsx‘）;
mm=xlsread（‘m.xlsx‘）;

w = xlsread（‘rate.xlsx‘1）; %频率
kesi  = [5e-25e-2];
c_alpha = 2*w（1）*w（2）*（kesi（2）*w（1）-kesi（1）*w（2））/（w（1）^2-w（2）^2）;
c_belta = 2*（kesi（1）*w（1）-kesi（2）*w（2））/（w（1）^2-w（2）^2）;
cc = c_alpha*kk+ c_belta*mm;  % 阻尼阵：采用瑞利结构阻尼

nt=20;
dt=0.01;
time=0:dt:nt*dt;

dof=30;
q0=zeros（dof1）;
dq0=zeros（dof1）;
bcdof=zeros（dof1）;



fd=zeros（dofnt）;
q0=zeros（dof1）;
for i=1:1
    q0（281）=1;
end


% q0=xlsread（‘q0.xlsx‘）;
% dq0=xlsread（‘dq0.xlsx‘）;
        
[sdofn2]=size（kk）;

dsp=zeros（sdofnt）;
acc=zeros（sdofnt）;
vel=zeros（sdofnt）;


dsp（:1）=q0;
vel（:1）=dq0;

alpha=0.5;
beta=0.5;

acc（:1）=inv（mm）*（fd