EQX 快速上手指南（一）线性单自由度体系减震分析

2020-03-04 本文已影响0人 RODS动力有限元

EQX 快速上手

EQX 简介

EQX 是基于Matlab的多自由度体系地震响应求解工具箱。可以很方便地对线性/非线性结构进行地震作用下的动力响应分析。
使用 EQX 请关注微信公众号 “EQ减震” 后回复 “eqx” 获得下载地址。

简单算例

模型

单自由度惯容减震结构体系。运动方程为
$\mathbf{M \ddot{u}}(t) + \mathbf{C \dot{u}}(t) + \mathbf{K u}(t) = -\mathbf{M_p E \ a_g}(t)$
当不设置惯容系统时
$\begin{aligned} \mathbf{M} &= \begin{bmatrix} m \end{bmatrix}\\ \mathbf{C} & = \begin{bmatrix} c \end{bmatrix}\\ \mathbf{K} & = \begin{bmatrix} k \end{bmatrix}\\ \mathbf{M_p} & = \begin{bmatrix} m \end{bmatrix}\\ \mathbf{E} & = \begin{pmatrix} 1 \end{pmatrix}\\ \end{aligned}$
当设置惯容系统（SPIS-II）时
$\begin{aligned} \mathbf{M} &= \begin{bmatrix} m & 0\\ 0 & m_{in} \end{bmatrix}\\ \mathbf{C} & = \begin{bmatrix} c &0\\ 0 & c_d \end{bmatrix}\\ \mathbf{K} & = \begin{bmatrix} k + k_s & -k_s \\ -k_s & k_s \end{bmatrix}\\ \mathbf{M_p} & = \begin{bmatrix} m&0\\ 0&0 \end{bmatrix}\\ \mathbf{E} & = \begin{pmatrix} 1\\ 1 \end{pmatrix}\\ \end{aligned}$

EQX求解

EQX 线性响应求解函数

x = eqss( M,K,C,Mp,ag,dt ); % 状态空间精确积分法
x = eqrk4( M,K,C,Mp,ag,dt ); % 4阶龙格-库塔法
x = eqnmk( M,K,C,Mp,ag,dt ); % Newmark-β法
...

使用 EQX 求解后可得到一个数值矩阵
$\mathbf{x} = \begin{pmatrix} \mathbf{u}^T & \mathbf{\dot u} ^T \end{pmatrix}$

提取第 $i$ 自由度的位移和速度响应时程

u_i = x(:,i);
v_i = x(:,N+i); % N 为总自由度数

EQX求解示例

读取地震动数据

ag = load('EQ-S-1.txt');
dt = 0.005;
n = length(ag);
t = linspace(0,dt*(n-1),n);

设置结构参数

zeta = 0.02;
m = 1.0;
k = m*100;
c = 2.0*zeta*sqrt(m*k);
omg = sqrt(k/m);
T = 2.0*pi/omg;

求解不设置惯容系统时的响应

M0 = m;
K0 = k;
C0 = c;
Mp0 = m;
x0 = eqss( M0,K0,C0,Mp0,ag,dt );

求解设置惯容系统后的响应

% 惯容系统参数设置
mu = 0.2;
kp = mu/(1.0-mu);
xi = mu/2.0*sqrt(kp/2.0);

m_in = mu*m;
k_s = kp*k;
c_d = 2.0*xi*sqrt(m*k);

% 建立矩阵并求解
M = [m 0;
     0 m_in];
C = [c 0;
     0 c_d];
K = [k+k_s -k_s;
     -k_s   k_s];
Mp = [m 0;
      0 0];
x = eqss( M,K,C,Mp,ag,dt );

结果可视化

位移响应

u0 = x0(:,1); % Deformation of original structure
u = x(:,1); % Deformation of damped structure

h = figure('Name','位移响应');
h.Units = 'normalized';
h.Position = [0.05,0.2,0.6,0.4];
plot(t,u0,'DisplayName','原结构');
hold on;
plot(t,u,'DisplayName','惯容减震结构');
hold on;
xlabel('Time');
ylabel('Deformation');
grid on;
hold on;
legend;

sdof_disp.png

滞回曲线

u_in = x(:,2); % Deformation of inerter
f = k_s*(u - u_in); % Force of inerter system
h1 = figure('Name','惯容系统滞回曲线');
h1.Units = 'normalized';
h1.Position = [0.65,0.2,0.3,0.4];
plot(u,f);
xlabel('变形');
ylabel('力');
grid on;

sdof_hy.png