结构动力响应求解的数值方法

本篇文章主要介绍结构动力学中结构响应求解的数值方法：杜哈梅积分、中心差分法和NewMark-β法，并利用python实现。对于matlab版本，将python代码直接扔给AI，应该不难得到，故此处不再给出。

• 地震波时程文件，单位为重力加速度g。

数值算法

杜哈梅积分

中心差分法

NewMark-β法

示例

#### README ####
# 杜哈梅积分、中心差分法和Newmark-β法的python实现。
# 采用国际单位制（N、t、mm、s...）

mode= input("请输入选择的数值方法：")

import numpy as np
from scipy.linalg import eig
import matplotlib.pyplot as plt

# 结构参数
# 质量矩阵
M = np.diag([120,120,120,120,120,120,120,110])

# 刚度矩阵
k = np.array([300,260,260,260,200,200,200,200])*1e3
K = np.zeros((len(k),len(k)))
for i in range(len(k)-1):
    K[i,i] = k[i]+k[i+1]
    K[i+1,i] = -k[i+1]
    K[i,i+1] = -k[i+1]

K[len(k)-1,len(k)-1] = k[len(k) - 1]

# 阻尼矩阵
a0 = 0.6296
a1 = 0.003075
C = a0*M + a1*K

# 地震动时程曲线
time_step = 0.005
file_path = 'ElCentro.txt'
ground_acc = np.loadtxt(file_path)
ground_acc = ground_acc*9.8*1e3
time = np.arange(len(ground_acc))*time_step
x0 = 0
v0 = 0


if mode == "Duhamel":
    #### 采用杜哈梅积分 ####
    # 由于杜哈梅积分只适用于单自由度，这里先利用振型分解法将位移解耦，计算每个振型的响应时程，再叠加。
    # 因为python计算的特征值并没有按大小排序，而matlab是排序了的。
    eigenvalues, eigenvectors = eig(K, M)
    sorted_indices = np.argsort(np.real(eigenvalues))
    sorted_eigenvalues = eigenvalues[sorted_indices]
    sorted_eigenvectors = eigenvectors[:, sorted_indices]

    # 初始化响应数组
    num_modes = len(sorted_eigenvalues)
    response = np.zeros(len(ground_acc))

    # 逐个振型进行计算并叠加响应
    for mode_idx in range(num_modes):
        mode = np.real(sorted_eigenvectors[:, mode_idx])  
        mode = mode / mode[-1]  # 规范化处理：顶层振型位移为1
        
        # 计算该振型的广义质量、刚度和阻尼矩阵
        M1 = mode @ M @ mode
        K1 = mode @ K @ mode
        C1 = mode @ C @ mode

        # 计算该振型的固有频率、阻尼比和阻尼频率
        w1 = np.sqrt(K1 / M1)
        s1 = C1 / (2 * M1 * w1)
        wd1 = np.sqrt(1 - s1**2) * w1

        # 杜哈梅积分计算响应
        x_mode = np.zeros(len(ground_acc))
        for i in range(len(ground_acc)):
            temp = 0
            for j in range(i + 1):
                p1 = -mode @ M @ np.ones(len(k)) * ground_acc[j]
                temp += np.exp(-s1 * w1 * (i - j) * time_step) * np.sin(wd1 * (i - j) * time_step) * p1 *time_step/ (M1 * wd1)
            x_mode[i] = temp

        # 将该振型的响应加到总响应中。能这么做的原因是因为所有振型中顶层的位移都是1，故直接将振型响应时程直接相加即可。
        # 根据这一原理，我们可以直接求出任一层的相对位移时程曲线。
        response += x_mode

    # 绘制顶层位移随时间变化曲线
    plt.figure()
    plt.plot(np.concatenate((time, [len(k) * time_step])), np.concatenate((np.array([0]), response)))
    plt.xlabel(r"$\it{t} (s)$")
    plt.ylabel(r"$\it{x} (mm)$")
    plt.savefig("顶层位移计算结果.png")
    plt.show()
elif mode == "CentralDifference":
    #### 中心差分法 ####
    C_acc = np.zeros((len(k), len(time)))  
    C_v = np.zeros((len(k), len(time)))   
    C_x = np.zeros((len(k), len(time)))    

    # 初始位移、速度
    C_x[:, 0] = x0  
    C_v[:, 0] = v0  

    C_acc[:, 0] = np.linalg.inv(M)@(-M @ np.ones(len(k)) * ground_acc[0] -C @ C_v[:, 0] - K@C_x[:, 0])
    C_x_minus1 = C_x[:, 0] - time_step * C_v[:, 0] + C_acc[:, 0] * time_step**2 / 2

    C_K = M / (time_step**2) + C / (2 * time_step)
    C_p_a = M / (time_step**2) - C / (2 * time_step)
    C_p_b = K - 2 * M / (time_step**2)

    inv_C_K = np.linalg.inv(C_K)

    for i in range(len(time) - 1):
        if i == 0:
            pi = -M @ np.ones(len(k)) * ground_acc[i]  
            C_pi = pi - C_p_a @ C_x_minus1 - C_p_b @ C_x[:, i]
        else:
            pi = -M @ np.ones(len(k)) * ground_acc[i]  
            C_pi = pi - C_p_a @ C_x[:, i - 1] - C_p_b @ C_x[:, i]
        
        C_x[:, i + 1] = inv_C_K @ C_pi  

        if i > 0:
            C_v[:, i] = (C_x[:, i + 1] - C_x[:, i - 1]) / (2 * time_step)  
            C_acc[:, i] = (C_x[:, i + 1] + C_x[:, i - 1] - 2 * C_x[:, i]) / (time_step**2) 

    plt.figure(2)
    plt.plot(time, C_x[len(k)-1,:])  
    plt.xlabel(r"$\it{t} (s)$")
    plt.ylabel(r"$\it{x} (mm)$")
    plt.show()
elif mode == "Newmark":
    #### Newmark-β法 ####
    N_acc = np.zeros((len(k), len(time)))  
    N_v = np.zeros((len(k), len(time)))   
    N_x = np.zeros((len(k), len(time)))    

    N_x[:, 0] = x0  
    N_v[:, 0] = v0  
    
    # 选用无条件稳定的参数
    beta = 0.25
    gamma = 0.5

    N_md = M + gamma * C * time_step + beta * K * time_step **2

    N_acc[:, 0] = np.linalg.inv(M) @(-M @ np.ones(len(k)) * ground_acc[0] -C @ N_v[:, 0] - K@N_x[:, 0])
    inv_N_md = np.linalg.inv(N_md)

    for i in range(len(time)-1):
        temp_v = N_v[:,i] + (1-gamma)*N_acc[:,i]*time_step
        temp_x = N_x[:,i] + N_v[:,i]*time_step +(0.5 - beta)*N_acc[:,i]*time_step**2
        N_acc[:,i+1] = inv_N_md @ (-M @ np.ones(len(k)) * ground_acc[i+1] -C @ temp_v - K@temp_x)
        N_v[:,i+1] = temp_v + gamma* N_acc[:,i+1]*time_step
        N_x[:,i+1] = temp_x + beta * N_acc[:,i+1] *time_step**2

    plt.figure(3)
    plt.plot(time, N_x[len(k)-1,:])  # 绘制顶层位移随时间变化
    plt.xlabel(r"$\it{t} (s)$")
    plt.ylabel(r"$\it{x} (mm)$")
    plt.show()

由于杜哈梅积分利python实现起来非常耗时，所以此处也给出其matlab版本。

clear all
clc
% 结构参数
M = diag([120, 120, 120, 120, 120, 120, 120, 110]);  % 质量矩阵
k = [300, 260, 260, 260, 200, 200, 200, 200] * 10^3;  % 刚度矩阵
K = zeros(length(k), length(k));  % 初始化刚度矩阵

% 刚度矩阵
for i = 1:length(k)-1
    K(i, i) = k(i) + k(i+1);
    K(i+1, i) = -k(i+1);
    K(i, i+1) = -k(i+1);
end
K(end, end) = k(end);  

% 阻尼矩阵
a0 = 0.6296;
a1 = 0.003075;
C = a0 * M + a1 * K;  

% 地震动时程曲线
time_step = 0.005;
ground_acc = load('ElCentro.txt');  % 读取地震加速度数据
ground_acc = ground_acc * 9.8 * 10^3;  % 转换为单位为 mm/s^2 的地震加速度
time = (0:length(ground_acc)-1) * time_step;  

% 计算特征值和特征向量
[eigenvectors, eigenvalues] = eig(K, M);

% 初始化总响应
num_modes = length(eigenvalues);
response = zeros(length(ground_acc), 1);
response_test = zeros(length(ground_acc),length(k));

% 逐个振型进行计算并叠加响应
for mode_idx = 1:num_modes
    mode = real(eigenvectors(:, mode_idx));  % 当前振型
    mode = mode / mode(end);  % 规范化处理：顶层振型位移为1

    % 计算该振型的广义质量、刚度和阻尼矩阵
    M1 = mode' * M * mode;
    K1 = mode' * K * mode;
    C1 = mode' * C * mode;

    % 计算该振型的固有频率、阻尼比和阻尼频率
    w1 = sqrt(K1 / M1);
    s1 = C1 / (2 * M1 * w1);
    wd1 = sqrt(1 - s1^2) * w1;

    % 计算该振型的响应
    x_mode = zeros(length(ground_acc), 1);  % 地震响应时程
    for i = 1:length(ground_acc)
        temp = 0;
        for j = 1:i
            p1 = -mode' * M * ones(length(k), 1) * ground_acc(j);
            temp = temp + exp(-s1 * w1 * (i - j) * time_step) * sin(wd1 * (i - j) * time_step) * p1 *time_step/ (M1 * wd1);
        end
        x_mode(i) = temp;
    end

    % 将该振型的响应加到总响应中
    response = response + x_mode;
end

% 绘制顶层位移随时间变化曲线
figure;
plot([time, time(end) + time_step], [0; response]);
xlabel('\it{t}\rm{(s)}');
ylabel('\it{x}\rm{(mm)}');
saveas(gcf, 'Duhamel_Results_matlab.png');

结果

三种方法得到的结果如下图所示。（以顶层相对位移为例）

杜哈梅积分结果(matlab):

中心差分法结果：

Newmark-β法结果：