编写线性振动Newmark法程序

以下是一个简单的Matlab程序，用于求解线性振动问题的Newmark法。程序中采用了四节点平面应力单元，可以根据实际情况进行修改。程序中的参数可以根据不同的问题进行调整。

% 定义有限元模型
L = 1; % 梁的长度
h = 0.2; % 梁的高度
b = 0.1; % 梁的宽度
E = 2.1e11; % 杨氏模量
nu = 0.3; % 泊松比
rho = 7800; % 密度
A = b*h; % 截面面积
I = b*h^3/12; % 惯性矩
k = [E*A/L,0,0,-E*A/L,0,0;
     0,E*I/L,0,0,-E*I/L,0;
     0,0,E*I/L,0,0,-E*I/L;
     -E*A/L,0,0,E*A/L,0,0;
     0,-E*I/L,0,0,E*I/L,0;
     0,0,-E*I/L,0,0,E*I/L];

% 定义初始条件和加载条件
dt = 0.01; % 时间步长
t = 0:dt:5; % 时间步
n = length(t); % 时间步数
u0 = zeros(12,1); % 初始位移和速度
v0 = zeros(12,1);
f0 = zeros(12,1); % 初始外力
f0(6) = 1000; % 在第6个节点施加集中单一频率动载荷

gamma = 0.5; % 加速度积分参数
beta = 0.25; % 加速度积分参数
alpha = 0.5*beta/dt^2;
a0 = u0*0; % 初始加速度
d0 = u0*0; % 初始阻尼力

% 初始化输出结果
u = zeros(12,n); % 位移
v = zeros(12,n); % 速度
a = zeros(12,n); % 加速度
f = zeros(12,n); % 外力

% Newmark法求解
for i = 2:n
    f(:,i) = f0; % 更新外力
    a(:,i-1) = a0; % 更新加速度
    d(:,i-1) = 2*gamma*beta/dt*a(:,i-1)+(1-gamma)*d(:,i-1); % 更新阻尼力
    a(:,i) = inv(k+(alpha+d(:,i-1))*eye(12))*(f(:,i)-k*u(:,i-1)-d(:,i-1)*v(:,i-1)-(1-alpha)*a(:,i-1)-alpha*a(:,i-2)); % 求解加速度
    v(:,i) = gamma*dt*a(:,i)+(1-gamma)*v(:,i-1)+dt*d(:,i-1); % 求解速度
    u(:,i) = u(:,i-1)+dt*v(:,i-1)+(0.5-gamma)*dt^2*a(:,i-1); % 求解位移
end

% 输出结果
figure;
plot(t,u(6,:),'b-');
xlabel('time (s)');
ylabel('displacement (m)');
title('Displacement response of cantilever beam under harmonic load');

请注意，此程序仅供参考，实际情况可能需要进行修改和调整。如果您遇到了问题，请参考Matlab官方文档或相关的数值计算书籍，或者向专业人士寻求帮助。