反应谱求解

反应谱求解反应谱定义参考书籍线性加速度法线性加速度法求解程序响应谱的求解示例

反应谱定义

反应谱：是在给定的地震加速度作用期间，单质点系统的最大反应随系统自振周期变化的曲线。

更直观的定义：一组具有相同阻尼、不同自振周期的单质点系统，在某一地震时程作用下的最大反应，为该地震的反应谱，分为加速度反应谱、速度反应谱与位移反应谱。

参考书籍

接触反应谱主要源于力学分析工作中的谱分析，主要是抗震分析。

有关于反应谱和地震相关的学习主要参考书籍《地震动的谱分析入门》，作者为日本的大崎顺彦，写的非常好。

这本书里面给出了地震动加速度激励下线性单质点弹簧阻尼系统的响应的求解方法和Fortran代码，用的是线性加速度法，是一种直接积分法。

线性加速度法

$x$ $\omega$ $h$ $\ddot{y}(t)$ ，则系统的运动方程可表示为：

\ddot{x}+2h\omega\dot{x}+{\omega}^2x=-\ddot{y}(t)\qquad(a)

$\ddot{y}(t)$ $\Delta t$ $t_1$ $t_2$ $t_i$ $\ddot{y}_1$ $\ddot{y}_2$ $\ddot{y}_i$ $\ddot{y}_i$ $\ddot{y}_{i+1}$ $t_i$ $t_i$ $t_{i+1}$ $\tau$ $\Delta\ddot{y}=\ddot{y}_{i+1}-\ddot{y}_i$ $\ddot{y}(t)$ 表示为：

\ddot{y}(t)=\frac{\Delta\ddot{y}}{\Delta t}+\ddot{y}_t

在该区间内公式1成为：

\ddot{x}(\tau)+2h\omega\dot{x}(\tau)+\omega^2x(\tau)=-\frac{\Delta\ddot{y}}{\Delta t}\tau+\ddot{y}\qquad 0\le\tau\le\Delta t

$t_i$ 的初始时刻初始条件为

\tau=0,x=x_i,\dot{x}=\dot{x}_i

$\tau=\Delta t$ $t_{i+1}=t_i+\Delta t$ 时的相对位移及相对速度为：

x_{i+1}=A_{11}x_i+A_{12}\dot{x}_i+B_{11}\ddot{y}_i+B_{12}\ddot{y}_{i+1}\qquad\quad\\ \dot x_{i+1}=A_{21}x_i+A_{22}\dot{x}_i+B_{21}\ddot{y}_i+B_{22}\ddot{y}_{i+1}\qquad(b)

其中：

A_{11}=e^{-h\omega \Delta t}(cos\omega_d\Delta t+\frac{h\omega}{\omega_d}sin\omega_d\Delta t)

A_{12}=e^{-h\omega \Delta t}\frac{1}{\omega_d}sin\omega_d\Delta t

A_{21}=-e^{-h\omega \Delta t}\frac{\omega^2}{\omega_d}sin\omega_d\Delta t

A_{22}=e^{-h\omega \Delta t}(cos\omega_d\Delta t-\frac{h\omega}{\omega_d}sin\omega_d\Delta t)

B_{11}=e^{-h\omega \Delta t}[(\frac{1}{\omega^2}+\frac{2h}{\omega^3\Delta t})cos\omega_d\Delta t+(\frac{h}{\omega\omega_d}-\frac{1-2h^2}{\omega^2\omega_d\Delta t})sin\omega_d\Delta t]-\frac{2h}{\omega^3\Delta t}

B_{12}=e^{-h\omega \Delta t}[-\frac{2h}{\omega^3\Delta t}cos\omega_d\Delta t+\frac{1-2h^2}{\omega^2\omega_d\Delta t}sin\omega_d\Delta t]-\frac 1{\omega^2}+\frac{2h}{\omega^3\Delta t}

B_{21}=e^{-h\omega \Delta t}[-\frac{1}{\omega^2\Delta t}cos\omega_d\Delta t-(\frac{h}{\omega\omega_d\Delta t}+\frac{1}{\omega_d})sin\omega_d\Delta t]+\frac{1}{\omega^2\Delta t}

B_{22}=e^{-h\omega\Delta t}[\frac{1}{\omega^2\Delta t}cos\omega_d\Delta t+\frac{h}{\omega\omega_d\Delta t}sin\omega_d\Delta t]-\frac{1}{\omega^2\Delta t}

$x_{i+1}$ 后，由（a）式可求得绝对加速度反应：

(\ddot x+\ddot y)_{i+1|}=-(2h\omega\dot x_{i+1}+\omega^2x_{i+1})\qquad(c)

当给定t=0时，反应的初值为：

x_1=0\\ \dot{x}_1=-\ddot{y}_1\Delta t\\ (\ddot x+\ddot y)=2h\omega\ddot y_1\Delta t

时，按（b）、（c）式就可逐次计算得到反应值（我在实际计算时全取0）。

另外，在计算（b）式各系数时，定义以下变量以便编制程序：

E=e^{-h\omega\Delta t}\\ SS=-h\omega sin\omega\Delta t-\omega_dcos\omega_d\Delta t\\ CC=-h\omega cos\omega_d\Delta t+\omega_dsin\omega_d\Delta t\\

S1=(E\cdot SS+\omega_d)/\omega^2\\ C1=(E\cdot CC+h\omega)/\omega^2

S2=(E\Delta t\cdot SS+h\omega S1+\omega_d\cdot C1)/\omega^2\\ C2=(E\Delta t\cdot CC+h\omega C1-\omega_d)\cdot S1)/\omega^2

S3=\Delta t\cdot S1-S2\\ C3=\Delta t\cdot C1-C2

线性加速度法求解程序


47
1
function [acc,vel,dis]=sdof_response(h,f,dt,ddy)
2
%sdof_response 求解单质点阻尼系统在地震加速度激励下的响应；
3
% h:阻尼比
4
% f：系统固有频率（Hz）
5
% dt：地震加速度时程采样时间间隔
6
% ddy：地震加速度时程
7
8
n_sample=length(ddy);
9
acc=zeros(1,n_sample);
10
vel=zeros(1,n_sample);
11
dis=zeros(1,n_sample);
12
13
w=2*pi*f;
14
w2=w*w;
15
hw=h*w;
16
wd=w*sqrt(1-h*h);
17
wdt=wd*dt;
18
E=exp(-hw*dt);
19
cwdt=cos(wdt);
20
swdt=sin(wdt);
21
a11=E*(cwdt+hw*swdt/wd);
22
a12=E*swdt/wd;
23
a21=-E*w2*swdt/wd;
24
a22=E*(cwdt-hw*swdt/wd);
25
26
ss=-hw*swdt-wd*cwdt;
27
cc=-hw*cwdt+wd*swdt;
28
s1=(E*ss+wd)/w2;
29
c1=(E*cc+hw)/w2;
30
s2=(E*dt*ss+hw*s1+wd*c1)/w2;
31
c2=(E*dt*cc+hw*c1-wd*s1)/w2;
32
s3=dt*s1-s2;
33
c3=dt*c1-c2;
34
b11=-s2/wdt;
35
b12=-s3/wdt;
36
b21=(hw*s2-wd*c2)/wdt;
37
b22=(hw*s3-wd*c3)/wdt;
38
39
acc(1)=0;
40
vel(1)=0;
41
dis(1)=0;
42
43
for i=2:n_sample
44
    dis(i)=a11*dis(i-1)+a12*vel(i-1)+b11*ddy(i-1)+b12*ddy(i);
45
    vel(i)=a21*dis(i-1)+a22*vel(i-1)+b21*ddy(i-1)+b22*ddy(i);
46
    acc(i)=-2*hw*vel(i)-w2*dis(i);
47
end

响应谱的求解

通过上述方法可以求得任意单个固有频率的系统在地震激励下的响应（含位移、速度、加速度），根据反应谱定义，输入一系列不同的待求解频率，分别取其反应的最大值，即可求得反应谱。

以加速度反应谱为例，可以通过以下代码求出：


12
1
function  rsacc = tts(h,freqs,acc,dt)
2
%TTS 求解给定的地震加速度时程的加速度反应谱
3
%h 系统阻尼比
4
%freqs 待求解的频率点
5
%acc 地震加速度时程
6
%dt 地震时程的采样时间间隔
7
n=length(freqs);
8
rsacc=zeros(1,n);
9
 for i=1:n
10
     rsacc(i)=max(abs(sdof_response(h,freqs(i),dt,acc)));
11
 end 
12
end

示例

以下以一段随机数作为地震波，采样频率100Hz，然后求其响应谱。


17
1
fs=100; %采样频率
2
T=60;   %时长
3
ddy=rand(fs*T,1)-0.5; %时程
4
h=0.02; %阻尼比
5
dt=1/fs;
6
freqs=logspace(log10(0.1),log10(100),100);
7
accResponseSpectrum=tts(h,freqs,ddy,dt);
8
subplot(2,1,2);
9
plot(freqs,accResponseSpectrum);
10
xlabel('频率(Hz)');
11
ylabel('加速度反应谱值');
12
grid on
13
subplot(2,1,1);
14
plot(0:dt:T-dt,ddy);
15
xlabel('时间(s)');
16
ylabel('加速度');
17
grid on;