结构力学基础概念：结构的动力分析：结构动力学的数值方法_2024-08-03

本文链接：https://blog.csdn.net/chenjj4003/article/details/148411985

结构力学基础概念：结构的动力分析：结构动力学的数值方法

结构动力学基础

基本概念

结构动力学的定义

结构动力学是研究结构在动态载荷作用下的响应和行为的学科。它关注结构的振动特性，包括频率、振幅、相位等，以及这些特性如何影响结构的安全性和性能。结构动力学不仅适用于建筑结构，也广泛应用于机械、航空、桥梁、船舶等工程领域。

动力学方程的建立

动力学方程是描述结构动力行为的基本数学模型。对于一个线性系统，动力学方程通常表示为：

$M\ddot{u} + C\dot{u} + Ku = F(t)$

其中：

$M$ 是质量矩阵，反映结构各部分的质量分布。
$C$ 是阻尼矩阵，表示能量耗散的效应。
$K$ 是刚度矩阵，反映结构的弹性性质。
$u$ 是位移向量， $\dot{u}$ 和 $\ddot{u}$ 分别是速度和加速度向量。
$F (t)$ 是随时间变化的外力向量。

振动理论

自由振动分析

自由振动分析研究的是当结构受到初始扰动后，在没有外力作用下的振动行为。对于无阻尼系统，自由振动方程简化为：

$M\ddot{u} + Ku = 0$

该方程的解通常为简谐振动，即：

$A\cos(\omega t + \phi)$

其中 $A$ 是振幅， $\omega$ 是角频率， $\phi$ 是初相位。通过求解特征值问题，可以得到系统的固有频率和振型。

示例：求解单自由度系统的固有频率和振型

假设一个单自由度系统，质量 $m = 1$ kg，刚度 $k = 100$ N/m。

import numpy as np

# 定义质量矩阵和刚度矩阵
M = np.array([[1]])
K = np.array([[100]])

# 求解特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.inv(M) @ K)

# 计算固有频率
omega = np.sqrt(eigenvalues)
f = omega / (2 * np.pi)

# 输出结果
print("固有频率：", f)
print("振型：", eigenvectors)

受迫振动分析

受迫振动分析考虑结构在持续外力作用下的响应。外力可以是周期性的，如风力、波浪力，也可以是非周期性的，如地震力。受迫振动方程为：

$M\ddot{u} + C\dot{u} + Ku = F(t)$

示例：求解单自由度系统在正弦外力作用下的响应

假设一个单自由度系统，质量 $m = 1$ kg，刚度 $k = 100$ N/m，阻尼 $c = 10$ Ns/m，受到 $F(t)=10\sin(10t)$ 的外力作用。

import numpy as np
import matplotlib.pyplot as plt

# 定义参数
m = 1
k = 100
c = 10
F0 = 10
omega_d = 10

# 定义时间向量
t = np.linspace(0, 10, 1000)

# 定义外力向量
F = F0 * np.sin(omega_d * t)

# 定义速度和加速度的初始条件
u0 = 0
v0 = 0

# 使用欧拉法求解微分方程
u = np.zeros_like(t)
v = np.zeros_like(t)
u[0] = u0
v[0] = v0
dt = t[1] - t[0]
for i in range(1, len(t)):
    a = (F[i] - c*v[i-1] - k*u[i-1]) / m
    v[i] = v[i-1] + a*dt
    u[i] = u[i-1] + v[i]*dt

# 绘制位移响应图
plt.plot(t, u)
plt.xlabel('时间 (s)')
plt.ylabel('位移 (m)')
plt.title('受迫振动响应')
plt.show()

振动的阻尼效应

阻尼效应描述了能量在振动过程中的耗散。阻尼可以是结构内部的材料阻尼，也可以是外部环境的空气阻尼或液体阻尼。阻尼的存在会减小振动的振幅，延长振动的衰减时间。

示例：比较不同阻尼比下单自由度系统的自由振动响应

假设一个单自由度系统，质量 $m = 1$ kg，刚度 $k = 100$ N/m，阻尼比分别为 $0.1$ 和 $0.5$ 。

import numpy as np
import matplotlib.pyplot as plt

# 定义参数
m = 1
k = 100
c1 = 0.1 * 2 * np.sqrt(m*k)  # 阻尼比0.1
c2 = 0.5 * 2 * np.sqrt(m*k)  # 阻尼比0.5
u0 = 1
v0 = 0

# 定义时间向量
t = np.linspace(0, 10, 1000)

# 使用欧拉法求解微分方程
u1 = np.zeros_like(t)
v1 = np.zeros_like(t)
u1[0] = u0
v1[0] = v0
dt = t[1] - t[0]
for i in range(1, len(t)):
    a1 = (-c1*v1[i-1] - k*u1[i-1]) / m
    v1[i] = v1[i-1] + a1*dt
    u1[i] = u1[i-1] + v1[i]*dt

u2 = np.zeros_like(t)
v2 = np.zeros_like(t)
u2[0] = u0
v2[0] = v0
for i in range(1, len(t)):
    a2 = (-c2*v2[i-1] - k*u2[i-1]) / m
    v2[i] = v2[i-1] + a2*dt
    u2[i] = u2[i-1] + v2[i]*dt

# 绘制位移响应图
plt.plot(t, u1, label='阻尼比0.1')
plt.plot(t, u2, label='阻尼比0.5')
plt.xlabel('时间 (s)')
plt.ylabel('位移 (m)')
plt.title('自由振动响应')
plt.legend()
plt.show()

通过上述示例，我们可以观察到阻尼比对振动响应的影响，阻尼比越大，振动衰减越快。

数值方法在结构动力学中的应用

有限元法

有限元法的基本原理

有限元法（Finite Element Method, FEM）是一种广泛应用于结构动力学分析的数值方法。它将复杂的结构分解为多个简单的小单元，每个单元的力学行为可以用简单的数学模型描述。通过在每个单元上应用局部的微分方程，然后将这些方程耦合起来，形成整个结构的全局方程，从而求解结构的动力响应。

原理概述

离散化：将连续体结构离散为有限个单元，每个单元用节点表示。
形函数：定义单元内部位移与节点位移之间的关系。
刚度矩阵：基于单元的形函数和材料属性，计算单元的刚度矩阵。
全局方程：将所有单元的刚度矩阵组合成全局刚度矩阵，形成结构的动力学方程。

动力学问题的有限元分析

在结构动力学中，有限元法可以用来求解结构在动态载荷下的响应，包括位移、速度、加速度和应力等。动力学问题的有限元分析通常涉及以下步骤：

建立模型：定义结构的几何形状、材料属性和边界条件。
离散化：将结构离散为有限个单元。
求解动力学方程：基于有限元法，求解结构的动力学方程。
后处理：分析和可视化求解结果。

时间积分方法

显式时间积分

显式时间积分是一种直接求解动力学方程的方法，它不需要求解线性方程组，因此计算效率较高。显式方法基于当前时间步的位移和速度来预测下一时间步的状态，适用于瞬态动力学问题。

方法描述

中心差分法：是一种常用的显式时间积分方法，它使用中心差分来近似加速度。
稳定性条件：显式方法的稳定性受到时间步长的限制，时间步长必须足够小以保证数值稳定性。

隐式时间积分

隐式时间积分方法在求解动力学方程时，需要在每个时间步求解一个线性方程组，因此计算成本较高，但其稳定性条件较宽松，适用于长时间的动力学分析。

方法描述

Newmark方法：是一种常用的隐式时间积分方法，通过调整参数β和γ来控制数值精度和稳定性。
稳定性分析：隐式方法的稳定性主要由方程的条件数决定，通常不受时间步长的限制。

模态分析

模态叠加法

模态叠加法是结构动力学分析中的一种重要方法，它基于结构的模态（固有频率和振型）来求解结构的动力响应。通过将结构的动力响应表示为模态振型的线性组合，可以简化计算过程。

方法步骤

模态提取：求解结构的固有频率和振型。
模态响应：计算每个模态在动态载荷下的响应。
响应叠加：将所有模态的响应叠加起来，得到结构的总响应。

模态综合技术

模态综合技术是一种高效的模态分析方法，它通过选择关键模态来近似结构的总响应，从而减少计算量。

技术应用

子空间迭代法：用于提取结构的关键模态。
模态综合：将关键模态的响应叠加，近似结构的总响应。

模态分析在动力学中的应用

模态分析在结构动力学中有着广泛的应用，包括振动控制、结构优化和故障诊断等。通过模态分析，可以了解结构的动态特性，为设计和分析提供重要信息。

非线性动力学分析

非线性动力学的基本概念

非线性动力学分析考虑了结构的非线性行为，包括材料非线性、几何非线性和接触非线性等。非线性动力学问题的求解通常比线性问题复杂，需要更高级的数值方法。

非线性特性

材料非线性：材料的应力-应变关系不是线性的。
几何非线性：结构的变形影响其刚度。
接触非线性：结构部件之间的接触行为是非线性的。

非线性动力学的数值模拟

非线性动力学的数值模拟通常采用增量迭代法，逐步求解结构在动态载荷下的非线性响应。

模拟步骤

增量加载：将动态载荷分解为多个小的增量。
迭代求解：在每个增量下，通过迭代求解非线性方程，直到满足收敛条件。
结果更新：更新结构的状态，包括位移、速度和加速度等。

非线性动力学分析的挑战与解决方案

非线性动力学分析面临的主要挑战包括收敛性问题、计算成本和非线性行为的复杂性。为了解决这些挑战，可以采用以下策略：

改进的迭代算法：如Newton-Raphson法，提高非线性方程的收敛速度。
自适应时间步长：根据结构响应的动态特性自动调整时间步长，减少计算量。
并行计算：利用多核处理器或分布式计算资源，加速非线性动力学分析的计算过程。

示例：有限元法求解简支梁的动力响应

假设我们有一根简支梁，长度为10米，截面为矩形，宽度为0.1米，高度为0.2米，材料为钢，弹性模量为200GPa，密度为7850kg/m³。梁受到一个正弦波动态载荷的作用，载荷频率为1Hz，幅值为10kN。我们使用有限元法求解梁的动力响应。

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt

# 材料和几何参数
E = 200e9  # 弹性模量，单位：Pa
rho = 7850  # 密度，单位：kg/m³
b = 0.1  # 宽度，单位：m
h = 0.2  # 高度，单位：m
L = 10  # 长度，单位：m
n_elements = 10  # 单元数量
n_nodes = n_elements + 1  # 节点数量
element_length = L / n_elements  # 单元长度

# 动态载荷参数
f = 1  # 载荷频率，单位：Hz
A = 10e3  # 载荷幅值，单位：N
t_max = 10  # 分析时间，单位：s
dt = 1e-3  # 时间步长，单位：s
t = np.arange(0, t_max, dt)

# 动态载荷
F = A * np.sin(2 * np.pi * f * t)

# 刚度矩阵和质量矩阵
K = np.zeros((n_nodes, n_nodes))
M = np.zeros((n_nodes, n_nodes))

# 计算刚度矩阵和质量矩阵
for i in range(n_elements):
    # 单元刚度矩阵
    k = (E * b * h / element_length) * np.array([[1, -1], [-1, 1]])
    # 单元质量矩阵
    m = (rho * b * h * element_length / 6) * np.array([[2, 1], [1, 2]])
    # 更新全局矩阵
    K[i:i+2, i:i+2] += k
    M[i:i+2, i:i+2] += m

# 边界条件
K[0, :] = 0
K[-1, :] = 0
K[:, 0] = 0
K[:, -1] = 0
K[0, 0] = 1
K[-1, -1] = 1

# 时间积分方法：Newmark方法
beta = 0.25
gamma = 0.5
alpha_m = beta / dt**2
alpha_k = beta / dt + gamma / 2

# 初始化
u = np.zeros(n_nodes)
v = np.zeros(n_nodes)
a = np.zeros(n_nodes)

# 时间积分
for i in range(len(t) - 1):
    # 预测位移和速度
    u_pred = u + v * dt + 0.5 * a * dt**2
    v_pred = v + (0.5 - gamma) * a * dt
    # 计算加速度
    a_pred = spsolve(K + alpha_m * M + alpha_k * K, -M @ (u_pred - u) / dt**2 - K @ (v_pred - v) / dt - F[i])
    # 更新状态
    u = u_pred + (1 - 2 * beta) * dt**2 * a_pred
    v = v_pred + (1 - gamma) * dt * a_pred
    a = a_pred

# 可视化结果
plt.figure()
plt.plot(t, u[n_nodes//2], label='位移')
plt.plot(t, v[n_nodes//2], label='速度')
plt.plot(t, a[n_nodes//2], label='加速度')
plt.legend()
plt.show()

示例描述

上述代码示例展示了如何使用有限元法和Newmark时间积分方法求解简支梁的动力响应。首先，定义了梁的几何和材料参数，以及动态载荷的参数。然后，计算了梁的刚度矩阵和质量矩阵，应用了边界条件。接着，使用Newmark方法进行时间积分，逐步求解梁在动态载荷下的位移、速度和加速度响应。最后，使用matplotlib库可视化了梁中点的位移、速度和加速度响应。

通过这个例子，我们可以看到有限元法和时间积分方法在结构动力学分析中的应用，以及如何通过编程实现这些方法。