中心差分法求相对位移时程曲线

源码依附在文章最后，读者可跳过讲解部分可直接下载。

中心差分法原理

从上述资料中可以看出，以线弹性体系为例，有：

$$\begin{equation} \label{eq1} X _ { i + 1 } = ( \frac { M } { \Delta t ^ { 2 } } + \frac { C } { 2 \Delta t } ) ^ { - 1 } \left[ - M \{1\} \ddot { X } _ { g , i } - ( \frac { M } { \Delta t ^ { 2 } } - \frac { C } { 2 \Delta t } ) X _ { i - 1 } - ( K - \frac { 2 M } { \Delta t ^ { 2 } } ) X _ { i } \right] \end{equation}$$

只要我们知道了第i时刻、i-1时刻结构体系的位移(严格来说叫相对位移，下同)和第i时刻地震动加速度$\ddot { X } _ { g , i }$，就可以求出第i+1时刻结构的位移。上述资料假设了第0时刻的位移、速度为0，现不妨更为普遍的设其为$X_{0}$和$\dot{X}_{0}$,则有：

$$\begin{equation} \label{eq2} {\dot X}_0=\frac{X_1-X_{-1}}{2\Delta t} \end{equation}$$ $$\begin{equation} \label{eq3} {\ddot X}_0=\frac{X_1-2X_0+X_{-1}}{\left(\Delta t\right)^2} \end{equation}$$

联立即得：

$$\begin{equation} \label{eq4} X_{-1}=X_0-\frac12{\dot X}_0\Delta t+\frac12{\ddot X}_0\left(\Delta t\right)^2 \end{equation}$$

式中${\ddot X}_0$可根据第0时刻的多自由度结构体系运动微分方程求得:

$$\begin{equation} \label{eq5} M{\ddot X}_0+C{\dot X}_0+KX_{0}=-M\left\{1\right\}{\ddot x}_{g,0} \end{equation}$$

由此，启动条件便确定了，再根据(1)式递推即可。

多自由度结构体系的中心差分法的MATLAB源码如下：

function [x, v, a] = MDOF_CentralDifference(m, k, c, P, x0, v0, dt, n)
    x = zeros(size(m,1), n);  % 位移
    v = zeros(size(m,1), n);  % 速度
    a = zeros(size(m,1), n);  % 加速度

    x(:, 1) = x0;  % 设置初始位移
    v(:, 1) = v0;  % 设置初始速度
    a(:, 1) = m\(P(:, 1) - c * v(:, 1) - k * x(:, 1));  % 计算初始加速度

    x_1 = x(:, 1) - dt .* v(:, 1) + dt^2 / 2 .* a(:, 1);  % 计算 X_(-1)
    P_eq = zeros(size(m));  % 初始化等效力

    k_eq = m ./ dt^2 + c ./ (2 * dt);  % 有效弹簧系数
    b1 = m ./ dt^2 - c ./ (2 * dt);  % 一阶阻尼项系数
    b2 = k - 2 .* (m ./ dt^2);  % 二阶质量项系数

    P_eq(:, 1) = P(:, 1) - b1 * x_1 - b2 * x(:, 1);  % 计算第一步的等效力
    x(:, 2) = k_eq\P_eq(:, 1) ;  % 计算第二步的位移

    for i = 2:n-1
        P_eq(:, i) = P(:, i) - b1 * x(:, i-1) - b2 * x(:, i);  % 计算等效力
        x(:, i+1) = k_eq\P_eq(:, i);  % 计算位移
        v(:, i) = (x(:, i+1) - x(:, i-1)) ./ (2 * dt);  % 计算速度
        a(:, i) = (x(:, i+1) - 2 * x(:, i) + x(:, i-1)) ./ dt^2;  % 计算加速度
    end

    P_eq(:, n) = P(:, n) - b1 * x(:, n-1) - b2 * x(:, n);
    xx = k_eq\P_eq(:, n);
    v(:, n) = (xx - x(:, n-1)) ./ (2 * dt);  % 计算最后一步的速度
    a(:, n) = (xx - 2 * x(:, n) + x(:, n-1)) ./ dt^2;  % 计算最后一步的加速度
end

示例

已知各阶振型及周期如下图所示，

试利用中心差分法求其在ELCentro波作用下的位移时程曲线。

---

结构基本信息

采用国际单位制：t\N\mm\s;

质量矩阵

易知结构质量矩阵为对角矩阵，有:

M = diag([120, 120, 120, 120, 120, 120, 120, 110]);

刚度矩阵

采用层剪切模型，易知n个自由度糖葫芦串模型的刚度矩阵表达式如下所示：

$$\begin{equation} K=\left[\begin{array}{cccccc} k_1+k_2 & -k_2 & & & & \\ -k_2 & k_2+k_3 & -k_3 & & & \\ & -k_3 & \cdots & & & \\ & & & k_{n-2}+k_{n-1} & -k_{n-1} & \\ & & & -k_{n-1} & k_{n-1}+k_n & -k_n \\ & & & & -k_n & k_n \end{array}\right] \end{equation}$$

故生成刚度矩阵如下：

k = [300, 260, 260, 260, 200, 200, 200, 200]*1000 ; %各层层间刚度
K = zeros(8,8);
for i = 1:7       %层剪切模型的刚度矩阵
    K(i,i) = k(i) + k(i+1);
    K(i+1,i) = -k(i+1);
    K(i,i+1) = -k(i+1);
end
K(8,8) = k(8);

阻尼矩阵

采用瑞利阻尼，即$C=\alpha_0M+\alpha_1K$,其中两个系数$\alpha_0$和$\alpha_1$由下式确定：

$$\begin{equation} \alpha _ { 0 } = \frac { 2 ( \xi _ { j } \omega _ { i } - \xi _ { i } \omega _ { j } ) \omega _ { i } \omega _ { j } } { \omega _ { i } ^ { 2 } - \omega _ { j } ^ { 2 } } \end{equation}$$ $$\begin{equation} \alpha _ { 1 } = \frac { 2 ( \xi _ { i } \omega _ { i } - \xi _ { i } \omega _ { i } ) } { \omega _ { i } ^ { 2 } - \omega _ { j } ^ { 2 } } \end{equation}$$

其中，$\omega _ { i }$，$\xi _ { i }$表示第i振型的频率和阻尼比。
本题采用第1、2振型的频率和阻尼比，求解得到$\alpha_0 = 0.6296$和$\alpha_1 = 0.003075$.

C = 0.6296*M + 0.003075*K; 
xi = 0.05;  %各阶振型的阻尼比

读取地震动加速度时程信息

值得注意的是，ElCentro.txt文件中，加速度单位为g，也就是$9.8m/s^2$，注意转化。

a=textread('ElCentro.txt');  % 先将地震动加速度记录存放在a中，用于下面的处理
[i,j]=size(a);                   % 地震数组维数
RecordLength=i*j;               % 记录步数
A=zeros(RecordLength,2);        % 定义处理后的地震动加速度时程数组，有两列，第一列存放时间，第二列存放对应地震动加速度值
A(1:RecordLength,1)= 0:0.01:(RecordLength-1)*0.01;  % 在A中输出时间变量,式中的0.01为记录间隔
for n1 = 1:1:i                                           % 用双重for循环将a中的记录按顺序放入A中
    for n2=1:1:j
        A((n1-1)*j+n2,2)=a(n1,n2)*9.8*1000;     % g(9.8 m/s2) 转化为 mm/s2
    end
end

定义地震动等效动荷载矩阵

即：$P=-M\left\{1\right\}{\ddot x}_g$.

onevector = ones(size(M,2),1);
P = zeros(size(M,1),RecordLength);
for i = 1:RecordLength
    P(:,i) = -(M*onevector).*A(i,2);
end

调用中心差分法函数进行求解

x0=zeros(8,1);    % 相对运动的初位移(数值积分用)
v0=zeros(8,1);    % 相对运动的初速度(数值积分用)
dt = 0.01;  %数值步长(=地震动记录间隔)
[x,v,a] = MDOF_CentralDifference(M,K,C,P,x0,v0,dt,RecordLength);

绘制第1层的位移时程曲线

figure
layer_num = 1;  
plot(A(:,1),x(layer_num,:).')       %绘制第 layer_num 层的相对位移时程曲线
title(['Time-history curve of relative displacement for Layer ', num2str(layer_num)])   % 图形标题
xlabel('time (s)')                     % 图形x坐标
ylabel('relative displacement (mm)')  % 图形y坐标

上方layer_num是层号，可随意调整，输出该层位移的时程曲线。

---

文件下载

点击下载：ElCentro.txt、MDOF_CentralDifference.m、MDOF_RelativeDisp.m.
将上述三个文件放在同一个文件夹下即可在MATLAB上运行，没安装MATLAB的同学可以试试MATLAB online.