close all; clear all;
orient tall
data_fn=load('newhall_360.dat');
g=9.81; % gravity
dt=0.02; % time increment
t_end=30; %end time (sec) for plotting
[m,n]=size(data_fn);
k=1;
for i=1:m
for j=1:n
Ag(k)=data_fn(i,j)/100; %change units from cm/s^2 to m/s^2
Time(k)=dt*(k-1);
k=k+1;
end
end
N=t_end/dt+1;
plot(Time(1:N),Ag(1:N))
xlabel('Time (sec)')
ylabel('Acceleration (m/sec^2)')
title('1994 Northridge Newhall Record (Fault Normal Direction)')
%% Problem 1
strain = [9.8, 74, 124.2, 179.6];
Po = [2.02, 7.26, 10.93, 17.73];
uo = [0.327, 2.457, 4, 5.9];
wd = [0.763, 17.25, 38, 76.03];
k2 = wd./(pi*uo.^2);
k1 = sqrt((Po./uo).^2 - k2.^2);
beff = wd./(2*pi.*k1.*uo.^2);
iso_matrix = zeros(length(strain), 3);
iso_matrix(:,1) = k1';
iso_matrix(:,2) = k2';
iso_matrix(:,3) = beff';
% for i = 1:size(iso_matrix,2)
% f = iso_matrix(:,i);
% plot(strain', f)
% hold on
% end
% xlabel('Strain (%)')
% ylabel('K1, K2, B_{eff}')
% legend('k1', 'k2', 'B_{eff}')
%% Problem 1b
%w1 and w2
ug2dot = Ag; %convert to in/s^2
neq = length(Ag);
np = 2 ^ (ceil(log(neq)/log(2)));
ug2dot(neq+1:np) = zeros(1, np-neq);
ntt = length(ug2dot);
delta_w = (2*pi)/((ntt - 1)*dt);
w1 = 0:delta_w:ntt/2 * delta_w;
w2 = -(ntt/2 - 1)*delta_w : delta_w: -delta_w;
%compute fourier transform for Ag (ground acceleration)
Fourier_ground_motion = fft(ug2dot);
abs_fourier = abs(Fourier_ground_motion);
%range of strain values
new_strain = 0:0.1:300;
%finding omega
index_w1 = find(abs_fourier(1:size(w1,2)) == max(abs_fourier(1:size(w1,2))));
dominate_frequency = w1(index_w1);
nb = 10; %number of bearings
for i = 1:size(iso_matrix,2)
interpol_value(:,i) = spline(strain', iso_matrix(:,i), new_strain');
end
list = ["K1 (\Omega)", "K2 (\Omega)", "B_{eff}"];
title_list = ["Equivalent storage stiffness", "loss stiffness", "Performance Index"];
figure
for i = 1:size(iso_matrix,2)
subplot(3,1, i)
plot(new_strain', interpol_value(:,i))
xlabel("Strain")
ylabel(list(i))
title(title_list(i))
end
saveas(gcf, "stiffness.png")
%natural frequency per strain
%natural frequency
mass = 125e3; %units in kips.s^2/in (1 N = 2.248*10^-4, m/s^2 = 39.37 in/s^2)
%damping ratio per strain
% zeta_strain = w_o .* interpol_value(:,3)/dominate_frequency;
w = [w1 , w2];
initial_strain = 9.8;
%initial values
rubber_height = 5.25*0.0254; %unit is in inches
max_iteration = 1000;
tol = 1e-6;
for i = 1:max_iteration
k_inter = (interp1(new_strain, interpol_value(:,1), initial_strain)*1000*0.454*9.81)/0.0254;
w_inter = sqrt((nb * k_inter)/mass);
zeta_inter = interp1(new_strain, interpol_value(:,3), initial_strain);
zeta_strain = w_inter * zeta_inter/dominate_frequency;
%frequency ratio
beta = w ./w_inter;
H = -1./ (w_inter^2 * ((1 - beta.^2) + (2 * zeta_strain .* beta * 1i)));
%compute fourier transform for Ag (ground acceleration)
[a_fourier, v_fourier, u_fourier] = response_frequency_domain(w_inter, zeta_strain, ug2dot, dt);
strain_update = (max(abs(u_fourier))/rubber_height)*100;
%compute fourier transform for Ag (ground acceleration)
Fourier_ground_motion1 = fft(ug2dot + a_fourier);
abs_fourier1 = abs(Fourier_ground_motion1);
%figure
%plot(w1, abs_fourier(1:N))
%finding omega
index_w2 = find(abs_fourier1(1:size(w1,2)) == max(abs_fourier1(1:size(w1,2))));
dominate_frequency = w1(index_w2);
if abs(strain_update - initial_strain) <= tol
figure
plot(Time(1:N), u_fourier(1:N))
xlabel("Time (sec)")
ylabel("Relative Displacement (m)")
title("Time History of the displacement response")
saveas(gcf, "time_history.png")
figure
plot(w1, abs_fourier(1:size(w1,2)))
xlabel("Frequency")
ylabel("Acceleration")
title("Forcing Frequency at the converged shear strain")
saveas(gcf, "forcing_frequency.png")
break
end
initial_strain = strain_update;
end
disp("The converged shear strain (%)")
disp(strain_update)
disp("The equivalent storage stiffness (K1) at the converged shear strain (N/m)")
disp((interp1(new_strain, interpol_value(:,1), strain_update)*1000*0.454*9.81)/0.0254)
disp("The loss stiffness (K2) at the converged shear strain (N/m)")
disp((interp1(new_strain, interpol_value(:,2), strain_update)*1000*0.454*9.81)/0.0254)
disp("Natural Frequency at the converged shear strain (rad/s):")
disp(w_inter)
disp("Effective Damping ratio at the converged shear strain (%):")
disp(zeta_inter * 100)
disp("Damping ratio at the converged shear strain (%):")
disp(zeta_strain * 100)
%% Problem 2
nbs = 15;
total_rubber_thickness = 0.15; %convert the rubber thickness to inches
rubber_diameter = 0.5; %convert the rubber diameter to inches
weight = 3*(150 * 47.88) * (25 * 2) *(25 * 4) * (0.3048)^2;
Pressure = weight/ (nbs * (pi * rubber_diameter^2/4));
% % matrix of pressure, Geff and damping
prop_matrix = [0 , 500 , 1000, 1500; 96, 102.4, 96, 89.6; 11.5, 12, 12.5,14]; %extracted from the plots
%
Geff = interp1(prop_matrix(1,:), prop_matrix(2,:), 0.000145038*Pressure); %Geff interpolated for the pressure calculated
%
building_damping = interp1(prop_matrix(1,:), prop_matrix(3,:), 0.000145038*Pressure); %Geff interpolated for the pressure calculated
%
keff = 6894.76*(Geff * (nbs * (pi * rubber_diameter^2/4)))/total_rubber_thickness;
T_period = 2*pi* sqrt((weight/9.81)/keff);
disp("Isolation period of the building")
disp(T_period)
w_o2 = sqrt(keff/(weight/9.81));
%compute fourier transform for Ag (ground acceleration)
Fourier_ground_motion2 = fft(ug2dot);
abs_fourier2 = abs(Fourier_ground_motion2);
%finding omega
index_w3 = find(abs_fourier2(1:size(w1,2)) == max(abs_fourier2(1:size(w1,2))));
dominate_frequency2 = w1(index_w3);
zeta_strain = 0.01*building_damping*w_o2/dominate_frequency2;
for i = 1:max_iteration
%frequency ratio
%compute fourier transform for Ag (ground acceleration)
[a_fourier1, v_fourier1, u_fourier1] = response_frequency_domain(w_o2, zeta_strain, ug2dot, dt);
%compute fourier transform for Ag (ground acceleration)
Fourier_ground_motion2 = fft(ug2dot + a_fourier1);
abs_fourier3 = abs(Fourier_ground_motion2);
%figure
%plot(w1, abs_fourier(1:N))
%finding omega
index_w4 = find(abs_fourier3(1:size(w1,2)) == max(abs_fourier3(1:size(w1,2))));
dominate_frequency3 = w1(index_w4);
if abs(dominate_frequency3 - dominate_frequency2) <= tol
disp("maximum relative displacement of the building in meters")
max_dis = max(abs(u_fourier1));
disp(max_dis)
break
end
dominate_frequency2 = dominate_frequency3;
zeta_strain = 0.01*building_damping*w_o2/dominate_frequency2;
end
没有合适的资源?快使用搜索试试~ 我知道了~
减轻地震灾害matlab代码.zip
共47个文件
m:14个
jpg:11个
png:9个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 86 浏览量
2024-05-18
17:37:27
上传
评论
收藏 1.06MB ZIP 举报
温馨提示
1.版本:matlab2014/2019a/2021a 2.附赠案例数据可直接运行matlab程序。 3.代码特点:参数化编程、参数可方便更改、代码编程思路清晰、注释明细。 4.适用对象:计算机,电子信息工程、数学等专业的大学生课程设计、期末大作业和毕业设计。
资源推荐
资源详情
资源评论
收起资源包目录
减轻地震灾害matlab代码.zip (47个子文件)
减轻地震灾害matlab代码
Earthquake_Hazard_Mitigation-main
newhall_4
earthquake_loop.png 33KB
motion1DOF.m 218B
newhall_fn_plot.m 3KB
earthquake_loop.fig 87KB
harmonic_loop.fig 112KB
acc_dis.fig 89KB
equivalent_linear.png 46KB
ddis_loop.fig 51KB
newhall_360.dat 30KB
acc_dis.png 43KB
motionDOF.m 175B
harmonic_loop.png 25KB
ddis_loop.png 30KB
motion2DOF.asv 326B
duhamel_dis.png 40KB
motion2DOF.m 326B
duhamel.m 547B
newhall_3
forcing_frequency.png 32KB
newhall_fn_plot.m 7KB
stiffness.png 32KB
response_frequency_domain.m 825B
newhall_360.dat 30KB
time_history.png 37KB
newhall_fn_plot.asv 7KB
newhall-1
newhall_fn_plot.m 4KB
NewmarkIntegrator.m 621B
response_frequency_domain.m 788B
disp_three.jpg 37KB
acc_three.jpg 48KB
time_domain4096.jpg 49KB
frequency_domain4096.jpg 48KB
frequency_domain16348.jpg 49KB
newhall_360.dat 30KB
duhamel.m 543B
duhamel.asv 543B
newhall_2
plotResponse.asv 1KB
newhall_fn_plot.m 1KB
response_2.jpg 59KB
fixed.jpg 33KB
response_frequency_domain.m 788B
fixed_2b.jpg 31KB
response1.jpg 52KB
plotResponse.m 1KB
newhall_360.dat 30KB
relative2b.jpg 38KB
newhall_fn_plot.asv 2KB
relative_disp.jpg 37KB
共 47 条
- 1
资源评论
Matlab科研辅导帮
- 粉丝: 2w+
- 资源: 7668
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功