结构动力学仿真-主题046-斜拉桥动力特性与风致振动
·
主题046：斜拉桥动力特性与风致振动
摘要
斜拉桥作为现代桥梁工程中的重要桥型，其动力特性与风致振动问题一直是结构动力学研究的热点。本主题系统介绍斜拉桥的动力特性分析方法、风荷载模拟技术以及风致振动机理。内容涵盖斜拉桥有限元建模方法、模态分析技术、风荷载功率谱模型、抖振响应分析、颤振稳定性判别以及涡激振动预测等核心内容。通过Python仿真实现，读者可以深入理解斜拉桥在风荷载作用下的动力响应特性，掌握风-结构相互作用的基本原理和分析方法。
关键词
斜拉桥；动力特性；风致振动；模态分析；颤振；涡激振动；抖振；功率谱密度
1. 斜拉桥结构概述
1.1 斜拉桥的结构组成
斜拉桥是一种由主梁、索塔和斜拉索三部分组成的组合体系桥梁。其结构特点包括：
主梁：作为桥面系的承载构件，主梁承受车辆荷载和部分风荷载，并通过斜拉索将荷载传递到索塔。主梁可采用钢箱梁、混凝土箱梁或钢-混凝土组合梁等形式。
索塔：作为斜拉桥的主要承重构件，索塔承受由斜拉索传递的竖向力和水平力。索塔形式多样，包括H形塔、A形塔、钻石形塔和倒Y形塔等。
斜拉索：连接主梁和索塔的受拉构件，是斜拉桥的关键受力部件。斜拉索的布置形式有辐射式、竖琴式、扇形式和星形式等。
1.2 斜拉桥的动力特性特点
斜拉桥的动力特性具有以下显著特点：
频率密集性：由于斜拉索数量众多，斜拉桥具有密集的固有频率分布，特别是斜拉索的局部振动频率与主梁、索塔的整体振动频率相互耦合。
模态复杂性：斜拉桥的振动模态包括主梁的竖向弯曲、横向弯曲、扭转振动，索塔的纵向和横向振动，以及斜拉索的面内和面外振动等多种形式。
非线性特性：在大变形情况下，斜拉索的垂度效应、几何非线性以及主梁与斜拉索之间的耦合效应使得斜拉桥表现出明显的非线性动力特性。
1.3 斜拉桥风致振动的主要形式
斜拉桥在风荷载作用下可能发生以下几种风致振动：
抖振（Buffeting）：由大气湍流引起的随机强迫振动，是斜拉桥最常见的一种风致响应形式。抖振响应虽然振幅较小，但长期作用会导致结构疲劳。
涡激振动（Vortex-Induced Vibration, VIV）：当风流过钝体结构时，在结构后方形成周期性的卡门涡街，当涡脱频率与结构固有频率接近时，会诱发结构发生大幅振动。
颤振（Flutter）：一种自激振动现象，当风速超过临界颤振风速时，结构从风场中不断吸收能量，导致振幅迅速增大，最终可能造成结构破坏。
驰振（Galloping）：也是一种自激振动，通常发生在具有特殊截面形状的细长结构上，其特点是存在负阻尼效应。
2. 斜拉桥有限元建模方法
2.1 主梁建模
主梁是斜拉桥的主要承载构件，其有限元建模需要考虑以下因素：
梁单元选择：主梁通常采用三维梁单元或板壳单元进行离散。对于箱形截面主梁，可采用空间梁单元，考虑截面的剪切变形和翘曲效应。
质量矩阵：主梁的质量矩阵可采用一致质量矩阵或集中质量矩阵。一致质量矩阵能够更准确地反映质量分布，但计算量较大；集中质量矩阵计算效率高，但精度略低。
几何非线性：对于大跨度斜拉桥，需要考虑主梁的几何非线性效应，包括大位移、大转动以及初始应力刚度效应。
主梁的刚度矩阵可表示为：
Kbeam=Klinear+KgeoK_{beam} = K_{linear} + K_{geo}Kbeam=Klinear+Kgeo
其中，KlinearK_{linear}Klinear是线性刚度矩阵，KgeoK_{geo}Kgeo是几何刚度矩阵。
2.2 索塔建模
索塔是斜拉桥的关键承重构件，其建模要点包括：
变截面处理：索塔通常采用变截面设计，建模时需要考虑截面特性的变化。可采用多段等截面梁单元近似，或采用变截面梁单元。
基础约束：索塔基础的约束条件对结构动力特性有重要影响。通常采用弹性约束模拟桩基础，考虑地基的平动和转动刚度。
横梁连接：对于门式或H形索塔，横梁与塔柱的连接区域需要特殊处理，可采用刚臂单元或细化网格。
2.3 斜拉索建模
斜拉索是斜拉桥最具特色的构件，其建模方法包括：
等效弹性模量法：考虑斜拉索的垂度效应，采用Ernst公式计算等效弹性模量：
Eeq=E1+(γL)212σ3EE_{eq} = \frac{E}{1 + \frac{(\gamma L)^2}{12\sigma^3}E}Eeq=1+12σ3(γL)2EE
其中，EEE是斜拉索材料的弹性模量，γ\gammaγ是索的单位长度重量，LLL是索的水平投影长度，σ\sigmaσ是索的应力。
多段离散法：将每根斜拉索离散为多个单元，更准确地模拟索的垂度效应和局部振动。
索力初始平衡：斜拉桥的初始索力状态对结构刚度有重要影响，需要进行初始平衡分析确定各索的初始张力。
2.4 整体结构组装
斜拉桥整体结构的运动方程可表示为：
Mu¨+Cu˙+Ku=F(t)M\ddot{u} + C\dot{u} + Ku = F(t)Mu¨+Cu˙+Ku=F(t)
其中：
- MMM是整体质量矩阵
- CCC是整体阻尼矩阵
- KKK是整体刚度矩阵
- uuu是节点位移向量
- F(t)F(t)F(t)是节点外力向量
整体矩阵的组装采用直接刚度法，将主梁、索塔和斜拉索的单元矩阵按照节点自由度进行组装。
3. 斜拉桥模态分析
3.1 特征值问题求解
斜拉桥的固有频率和振型通过求解广义特征值问题获得：
(K−ω2M)ϕ=0(K - \omega^2 M)\phi = 0(K−ω2M)ϕ=0
其中，ω\omegaω是圆频率，ϕ\phiϕ是振型向量。
对于大型斜拉桥结构，通常采用子空间迭代法或Lanczos方法求解前若干阶模态。
3.2 斜拉桥的典型振型
斜拉桥的典型振型包括：
主梁竖向弯曲振型：主梁在竖向平面内的弯曲振动，通常是对称或反对称形式。第一阶竖向弯曲频率是斜拉桥设计的重要参数。
主梁横向弯曲振型：主梁在水平面内的弯曲振动，主要由风荷载或地震横向分量引起。
主梁扭转振型：主梁绕纵向轴的扭转振动，对于闭口箱梁截面，扭转刚度较大，扭转频率通常较高。
索塔振动振型：索塔在纵向和横向的弯曲振动，纵向振动与主梁竖向振动耦合，横向振动与主梁横向振动耦合。
斜拉索局部振型：单根斜拉索的面内和面外振动，频率通常较高且密集分布。
3.3 模态参与因子
模态参与因子用于评估各阶模态对结构响应的贡献程度：
Γi=ϕiTMrϕiTMϕi\Gamma_i = \frac{\phi_i^T M r}{\phi_i^T M \phi_i}Γi=ϕiTMϕiϕiTMr
其中，rrr是影响向量，描述各自由度在整体运动中的参与方式。
对于地震响应分析，通常要求前几阶模态的累积质量参与系数达到90%以上。
4. 风荷载模型
4.1 平均风与脉动风
风荷载可分为平均风分量和脉动风分量：
U(z,t)=Uˉ(z)+u(z,t)U(z,t) = \bar{U}(z) + u(z,t)U(z,t)=Uˉ(z)+u(z,t)
其中，Uˉ(z)\bar{U}(z)Uˉ(z)是高度zzz处的平均风速，u(z,t)u(z,t)u(z,t)是脉动风速。
平均风速剖面：采用对数律或幂律描述平均风速随高度的变化：
对数律：Uˉ(z)=u∗κln(zz0)\bar{U}(z) = \frac{u_*}{\kappa} \ln\left(\frac{z}{z_0}\right)Uˉ(z)=κu∗ln(z0z)
幂律：Uˉ(z)=Uˉ10(z10)α\bar{U}(z) = \bar{U}_{10} \left(\frac{z}{10}\right)^\alphaUˉ(z)=Uˉ10(10z)α
其中，u∗u_*u∗是摩擦速度，κ\kappaκ是卡门常数（约0.4），z0z_0z0是粗糙长度，Uˉ10\bar{U}_{10}Uˉ10是10m高度处的平均风速，α\alphaα是地面粗糙度指数。
4.2 脉动风功率谱
脉动风的随机特性通常用功率谱密度描述：
Davenport谱：
Su(f)=4kUˉ102fx2(1+x2)4/3S_u(f) = \frac{4k\bar{U}_{10}^2}{f} \frac{x^2}{(1+x^2)^{4/3}}Su(f)=f4kUˉ102(1+x2)4/3x2
其中，x=1200f/Uˉ10x = 1200f/\bar{U}_{10}x=1200f/Uˉ10，kkk是地表阻力系数。
Kaimal谱：
fSu(f)u∗2=200x(1+50x)5/3\frac{fS_u(f)}{u_*^2} = \frac{200x}{(1+50x)^{5/3}}u∗2fSu(f)=(1+50x)5/3200x
其中，x=fz/Uˉ(z)x = fz/\bar{U}(z)x=fz/Uˉ(z)。
Simiu谱：
fSu(f)u∗2=200x(1+50x)5/3\frac{fS_u(f)}{u_*^2} = \frac{200x}{(1+50x)^{5/3}}u∗2fSu(f)=(1+50x)5/3200x
4.3 空间相关性
风场的空间相关性采用相干函数描述：
Coh(f,Δy)=exp(−fΔyUˉ(z))Coh(f, \Delta y) = \exp\left(-\frac{f\Delta y}{\bar{U}(z)}\right)Coh(f,Δy)=exp(−Uˉ(z)fΔy)
其中，Δy\Delta yΔy是两点间的水平距离。
5. 抖振响应分析
5.1 抖振力模型
抖振力是由于脉动风引起的非定常气动力，可表示为：
Lb(t)=12ρU2BCL′u(t)UL_b(t) = \frac{1}{2}\rho U^2 B C_L' \frac{u(t)}{U}Lb(t)=21ρU2BCL′Uu(t)
Db(t)=12ρU2BCDu(t)UD_b(t) = \frac{1}{2}\rho U^2 B C_D \frac{u(t)}{U}Db(t)=21ρU2BCDUu(t)
Mb(t)=12ρU2B2CM′u(t)UM_b(t) = \frac{1}{2}\rho U^2 B^2 C_M' \frac{u(t)}{U}Mb(t)=21ρU2B2CM′Uu(t)
其中，ρ\rhoρ是空气密度，B是桥面宽度，B是桥面宽度，B是桥面宽度，C_L’、、、C_D、、、C_M’$分别是升力系数、阻力系数和力矩系数对攻角的导数。
5.2 频域分析方法
抖振响应的频域分析基于随机振动理论：
Sy(f)=∣H(f)∣2SF(f)S_y(f) = |H(f)|^2 S_F(f)Sy(f)=∣H(f)∣2SF(f)
其中，Sy(f)S_y(f)Sy(f)是响应功率谱，H(f)H(f)H(f)是频响函数，SF(f)S_F(f)SF(f)是抖振力功率谱。
频响函数为：
H(f)=1K−(2πf)2M+j2πfCH(f) = \frac{1}{K - (2\pi f)^2M + j2\pi fC}H(f)=K−(2πf)2M+j2πfC1
响应的均方根值为：
σy=∫0∞Sy(f)df\sigma_y = \sqrt{\int_0^\infty S_y(f) df}σy=∫0∞Sy(f)df
5.3 时域分析方法
时域分析通过生成脉动风时程，直接求解运动方程：
谐波叠加法：
u(t)=∑i=1N2Su(fi)Δfcos(2πfit+ϕi)u(t) = \sum_{i=1}^{N} \sqrt{2S_u(f_i)\Delta f} \cos(2\pi f_i t + \phi_i)u(t)=i=1∑N2Su(fi)Δfcos(2πfit+ϕi)
其中，ϕi\phi_iϕi是[0,2π][0, 2\pi][0,2π]区间均匀分布的随机相位角。
线性滤波法（ARMA模型）：
通过白噪声滤波生成具有目标功率谱的风速时程。
6. 颤振稳定性分析
6.1 颤振机理
颤振是一种自激振动现象，其机理是结构振动与气动力之间的耦合反馈。当结构振动时，改变了周围流场，产生附加气动力；这些气动力又反过来影响结构振动。当风速超过临界值时，系统从气流中不断吸收能量，导致振幅指数增长。
6.2 颤振导数
颤振分析采用Scanlan提出的颤振导数模型：
L=12ρU2B[KH1∗h˙U+KH2∗Bα˙U+K2H3∗α+K2H4∗hB]L = \frac{1}{2}\rho U^2 B \left[KH_1^* \frac{\dot{h}}{U} + KH_2^* \frac{B\dot{\alpha}}{U} + K^2H_3^* \alpha + K^2H_4^* \frac{h}{B}\right]L=21ρU2B[KH1∗Uh˙+KH2∗UBα˙+K2H3∗α+K2H4∗Bh]
M=12ρU2B2[KA1∗h˙U+KA2∗Bα˙U+K2A3∗α+K2A4∗hB]M = \frac{1}{2}\rho U^2 B^2 \left[KA_1^* \frac{\dot{h}}{U} + KA_2^* \frac{B\dot{\alpha}}{U} + K^2A_3^* \alpha + K^2A_4^* \frac{h}{B}\right]M=21ρU2B2[KA1∗Uh˙+KA2∗UBα˙+K2A3∗α+K2A4∗Bh]
其中，K=Bω/UK = B\omega/UK=Bω/U是折减频率，Hi∗H_i^*Hi∗和Ai∗A_i^*Ai∗是颤振导数，通过风洞试验识别。
6.3 颤振临界风速
颤振临界风速可通过复特征值分析确定。考虑气动力后的运动方程为：
Mu¨+(C−Ca)u˙+(K−Ka)u=0M\ddot{u} + (C - C_a)\dot{u} + (K - K_a)u = 0Mu¨+(C−Ca)u˙+(K−Ka)u=0
其中，CaC_aCa和KaK_aKa是气动阻尼矩阵和气动刚度矩阵。
求解复特征值问题，当某阶模态的阻尼比变为零时，对应的风速即为颤振临界风速。
7. 涡激振动分析
7.1 涡激振动机理
当风流过钝体结构时，在结构后方形成周期性的卡门涡街。涡脱频率可用Strouhal数描述：
fs=St⋅UDf_s = \frac{St \cdot U}{D}fs=DSt⋅U
其中，StStSt是Strouhal数，DDD是结构特征尺寸。
当涡脱频率fsf_sfs与结构固有频率fnf_nfn接近时，发生涡激共振：
fsfn≈1\frac{f_s}{f_n} \approx 1fnfs≈1
7.2 涡激力模型
涡激力可采用经验模型描述，如Scanlan经验非线性模型：
Fvortex=12ρU2D[Y1(K)(1−ϵy2D2)y˙U+Y2(K)yD]F_{vortex} = \frac{1}{2}\rho U^2 D \left[Y_1(K) \left(1 - \epsilon \frac{y^2}{D^2}\right) \frac{\dot{y}}{U} + Y_2(K) \frac{y}{D}\right]Fvortex=21ρU2D[Y1(K)(1−ϵD2y2)Uy˙+Y2(K)Dy]
其中，Y1Y_1Y1、Y2Y_2Y2和ϵ\epsilonϵ是经验参数。
7.3 涡激振动幅值预测
涡激振动的稳态幅值可通过求解非线性运动方程获得。对于单自由度系统，无量纲振幅为：
AD=1ϵY1−4πScY1\frac{A}{D} = \frac{1}{\sqrt{\epsilon}} \sqrt{\frac{Y_1 - 4\pi S_c}{Y_1}}DA=ϵ1Y1Y1−4πSc
其中，Sc=2mζ/(ρD2)S_c = 2m\zeta/(\rho D^2)Sc=2mζ/(ρD2)是Scruton数，mmm是单位长度质量，ζ\zetaζ是阻尼比。
8. Python仿真实现
8.1 斜拉桥有限元建模
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
import warnings
warnings.filterwarnings('ignore')
# 斜拉桥几何参数
L_span = 200 # 主跨长度 (m)
L_side = 100 # 边跨长度 (m)
H_tower = 80 # 索塔高度 (m)
n_cables = 10 # 每侧拉索数量
# 材料参数
E_steel = 2.0e11 # 钢材弹性模量 (Pa)
E_cable = 1.95e11 # 拉索弹性模量 (Pa)
rho_steel = 7850 # 钢材密度 (kg/m³)
rho_cable = 7850 # 拉索密度 (kg/m³)
8.2 模态分析
def solve_modal_analysis(K, M, n_modes=10):
"""
求解斜拉桥模态参数
"""
eigenvalues, eigenvectors = eigh(K, M, subset_by_index=[0, n_modes-1])
frequencies = np.sqrt(eigenvalues) / (2 * np.pi)
return frequencies, eigenvectors
8.3 风荷载模拟
def generate_wind_time_series(U_mean, Iu, Lu, duration, dt):
"""
生成脉动风时程（谐波叠加法）
"""
t = np.arange(0, duration, dt)
n_freq = 100
f = np.linspace(0.01, 5, n_freq)
df = f[1] - f[0]
# Kaimal谱
Su = (200 * f * Iu**2 * U_mean**2 / Lu) / (1 + 50 * f * Lu / U_mean)**(5/3)
u_fluct = np.zeros_like(t)
for i in range(n_freq):
phi = np.random.uniform(0, 2*np.pi)
u_fluct += np.sqrt(2 * Su[i] * df) * np.cos(2*np.pi*f[i]*t + phi)
return U_mean + u_fluct, t
8.4 抖振响应计算
def buffeting_response(M, C, K, F_buffeting, dt):
"""
计算抖振时域响应（Newmark-beta法）
"""
n_dof = len(M)
n_steps = len(F_buffeting[0])
u = np.zeros((n_dof, n_steps))
v = np.zeros((n_dof, n_steps))
a = np.zeros((n_dof, n_steps))
# Newmark-beta参数
beta = 0.25
gamma = 0.5
for i in range(1, n_steps):
# 预测
u[:, i] = u[:, i-1] + dt * v[:, i-1] + 0.5 * dt**2 * a[:, i-1]
v[:, i] = v[:, i-1] + dt * a[:, i-1]
# 求解加速度
F_eff = F_buffeting[:, i] - K @ u[:, i] - C @ v[:, i]
a[:, i] = np.linalg.solve(M, F_eff)
# 校正
u[:, i] += beta * dt**2 * a[:, i]
v[:, i] += gamma * dt * a[:, i]
return u, v, a
9. 工程应用案例
9.1 某斜拉桥动力特性分析
某双塔双索面斜拉桥，主跨400m，边跨200m，主梁采用钢箱梁，索塔采用H形混凝土塔。通过有限元分析获得前10阶固有频率和振型。
分析结果表明：
- 第一阶竖向弯曲频率为0.28Hz
- 第一阶横向弯曲频率为0.35Hz
- 第一阶扭转频率为0.62Hz
- 斜拉索的基频分布在0.8-2.5Hz范围内
9.2 风致振动控制措施
针对斜拉桥的风致振动问题，常用的控制措施包括：
气动措施：
- 优化主梁截面形状，采用流线型箱梁
- 设置导流板或抑流板
- 安装风嘴或扰流装置
机械措施：
- 安装调谐质量阻尼器（TMD）
- 设置粘滞阻尼器
- 采用高阻尼橡胶支座
结构措施：
- 增加结构阻尼
- 调整斜拉索张力
- 优化索塔刚度
10. 结论与展望
斜拉桥动力特性与风致振动分析是桥梁工程中的重要研究内容。本主题系统介绍了斜拉桥的有限元建模方法、模态分析技术、风荷载模型以及各类风致振动的分析方法。
通过Python仿真实现，可以深入理解斜拉桥在风荷载作用下的动力响应特性。在实际工程中，需要结合风洞试验和现场实测，验证和完善理论分析方法。
未来的研究方向包括：
- 考虑多物理场耦合的精细化分析
- 机器学习方法在风荷载预测中的应用
- 智能监测与健康评估技术
- 新型减振控制技术的开发
import matplotlib
matplotlib.use('Agg')
"""
案例1：斜拉桥动力特性分析
包含：
1. 斜拉桥简化有限元建模
2. 模态分析
3. 振型可视化
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.size'] = 10
plt.rcParams['axes.unicode_minus'] = False
def create_simplified_cable_bridge():
"""
创建简化斜拉桥模型
返回:
K: 整体刚度矩阵
M: 整体质量矩阵
nodes: 节点坐标
elements: 单元连接
element_types: 单元类型
"""
# 几何参数
L_span = 80 # 主跨长度 (m)
H_tower = 30 # 索塔高度 (m)
# 材料参数
E = 2.0e11 # 弹性模量 (Pa)
rho = 7850 # 密度 (kg/m³)
# 截面参数
A_girder = 2.0 # 主梁截面积 (m²)
I_girder = 5.0 # 主梁惯性矩 (m⁴)
A_tower = 4.0 # 索塔截面积 (m²)
I_tower = 20.0 # 索塔惯性矩 (m⁴)
A_cable = 0.01 # 拉索截面积 (m²)
# 创建节点
nodes = []
# 左塔节点 (2个节点)
nodes.append([-L_span/2, 0]) # 塔底
nodes.append([-L_span/2, H_tower]) # 塔顶
# 右塔节点 (2个节点)
nodes.append([L_span/2, 0]) # 塔底
nodes.append([L_span/2, H_tower]) # 塔顶
# 主梁节点 (9个节点)
n_girder = 9
for i in range(n_girder):
x = -L_span/2 + i * L_span / (n_girder - 1)
nodes.append([x, 0])
nodes = np.array(nodes)
n_nodes = len(nodes)
# 创建单元
elements = []
element_types = []
# 左塔单元
elements.append([0, 1])
element_types.append('tower')
# 右塔单元
elements.append([2, 3])
element_types.append('tower')
# 主梁单元
girder_start = 4
for i in range(n_girder - 1):
elements.append([girder_start + i, girder_start + i + 1])
element_types.append('girder')
# 斜拉索单元（简化模型，每塔3根索）
tower_top_left = 1
tower_top_right = 3
# 左侧拉索
for i in [1, 3, 5]:
girder_node = girder_start + i
elements.append([tower_top_left, girder_node])
element_types.append('cable')
# 右侧拉索
for i in [3, 5, 7]:
girder_node = girder_start + i
elements.append([tower_top_right, girder_node])
element_types.append('cable')
elements = np.array(elements)
# 组装整体矩阵（每个节点3个自由度：x, y, theta）
n_dof = 3 * n_nodes
K = np.zeros((n_dof, n_dof))
M = np.zeros((n_dof, n_dof))
for i, elem in enumerate(elements):
node_i, node_j = elem
xi, yi = nodes[node_i]
xj, yj = nodes[node_j]
L = np.sqrt((xj-xi)**2 + (yj-yi)**2)
elem_type = element_types[i]
if elem_type in ['girder', 'tower']:
I = I_girder if elem_type == 'girder' else I_tower
A = A_girder if elem_type == 'girder' else A_tower
# 梁单元刚度矩阵（简化）
k_elem = np.zeros((6, 6))
k_elem[0, 0] = k_elem[3, 3] = E * A / L
k_elem[0, 3] = k_elem[3, 0] = -E * A / L
k_elem[1, 1] = k_elem[4, 4] = 12 * E * I / L**3
k_elem[1, 4] = k_elem[4, 1] = -12 * E * I / L**3
k_elem[1, 2] = k_elem[2, 1] = 6 * E * I / L**2
k_elem[1, 5] = k_elem[5, 1] = 6 * E * I / L**2
k_elem[2, 4] = k_elem[4, 2] = -6 * E * I / L**2
k_elem[4, 5] = k_elem[5, 4] = -6 * E * I / L**2
k_elem[2, 2] = k_elem[5, 5] = 4 * E * I / L
k_elem[2, 5] = k_elem[5, 2] = 2 * E * I / L
# 梁单元质量矩阵（简化）
m_elem = np.zeros((6, 6))
m_elem[0, 0] = m_elem[3, 3] = rho * A * L / 3
m_elem[0, 3] = m_elem[3, 0] = rho * A * L / 6
m_elem[1, 1] = m_elem[4, 4] = rho * A * L / 3
m_elem[1, 4] = m_elem[4, 1] = rho * A * L / 6
m_elem[2, 2] = m_elem[5, 5] = rho * A * L * L / 30
dof_map = [3*node_i, 3*node_i+1, 3*node_i+2, 3*node_j, 3*node_j+1, 3*node_j+2]
for ii in range(6):
for jj in range(6):
K[dof_map[ii], dof_map[jj]] += k_elem[ii, jj]
M[dof_map[ii], dof_map[jj]] += m_elem[ii, jj]
elif elem_type == 'cable':
c = (xj-xi) / L
s = (yj-yi) / L
# 杆单元刚度矩阵
k_local = (E * A_cable / L) * np.array([[1, -1], [-1, 1]])
T = np.array([[c, s, 0, 0], [0, 0, c, s]])
k_elem = T.T @ k_local @ T
# 杆单元质量矩阵
m_local = (rho * A_cable * L / 6) * np.array([[2, 1], [1, 2]])
m_elem = T.T @ m_local @ T
dof_map = [3*node_i, 3*node_i+1, 3*node_j, 3*node_j+1]
for ii in range(4):
for jj in range(4):
K[dof_map[ii], dof_map[jj]] += k_elem[ii, jj]
M[dof_map[ii], dof_map[jj]] += m_elem[ii, jj]
return K, M, nodes, elements, element_types
def solve_modal(K, M, fixed_dofs, n_modes=8):
"""
求解模态参数
参数:
K: 刚度矩阵
M: 质量矩阵
fixed_dofs: 固定自由度列表
n_modes: 模态数
返回:
frequencies: 固有频率
modes: 模态振型
"""
n_dof = len(K)
free_dofs = [i for i in range(n_dof) if i not in fixed_dofs]
if len(free_dofs) == 0:
return np.array([]), np.array([])
K_reduced = K[np.ix_(free_dofs, free_dofs)]
M_reduced = M[np.ix_(free_dofs, free_dofs)]
# 检查矩阵是否奇异
if np.linalg.cond(K_reduced) > 1e15:
print("警告：刚度矩阵接近奇异，添加微小刚度")
K_reduced += np.eye(len(K_reduced)) * 1e-6
try:
eigenvalues, eigenvectors = eigh(K_reduced, M_reduced)
# 过滤掉零频率和负频率
valid_idx = eigenvalues > 1e-10
eigenvalues = eigenvalues[valid_idx]
eigenvectors = eigenvectors[:, valid_idx]
frequencies = np.sqrt(eigenvalues) / (2 * np.pi)
# 重构完整模态
n_modes_actual = min(n_modes, len(frequencies))
modes = np.zeros((n_dof, n_modes_actual))
modes[free_dofs, :n_modes_actual] = eigenvectors[:, :n_modes_actual]
return frequencies[:n_modes_actual], modes[:, :n_modes_actual]
except:
print("模态分析失败")
return np.array([]), np.array([])
def plot_bridge(nodes, elements, element_types, U=None, scale=1.0, title="斜拉桥"):
"""
绘制斜拉桥
"""
fig, ax = plt.subplots(figsize=(12, 6))
if U is not None:
nodes_def = nodes.copy()
nodes_def[:, 0] += U[::3] * scale
nodes_def[:, 1] += U[1::3] * scale
else:
nodes_def = nodes
for i, elem in enumerate(elements):
node_i, node_j = elem
x = [nodes_def[node_i, 0], nodes_def[node_j, 0]]
y = [nodes_def[node_i, 1], nodes_def[node_j, 1]]
elem_type = element_types[i]
if elem_type == 'girder':
ax.plot(x, y, 'b-', linewidth=3, label='主梁' if i == 0 else '')
elif elem_type == 'tower':
ax.plot(x, y, 'k-', linewidth=4, label='索塔' if i == 0 else '')
elif elem_type == 'cable':
ax.plot(x, y, 'r--', linewidth=1, alpha=0.7, label='拉索' if i == 0 else '')
ax.scatter(nodes[:, 0], nodes[:, 1], c='green', s=50, zorder=5)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title(title)
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
return fig, ax
def main():
"""
主函数
"""
print("=" * 70)
print("案例1：斜拉桥动力特性分析")
print("=" * 70)
# 创建模型
print("\n创建斜拉桥有限元模型...")
K, M, nodes, elements, element_types = create_simplified_cable_bridge()
print(f"模型信息：")
print(f" 节点数：{len(nodes)}")
print(f" 单元数：{len(elements)}")
print(f" 自由度数：{len(K)}")
# 绘制几何模型
fig1, _ = plot_bridge(nodes, elements, element_types, title="斜拉桥有限元模型")
plt.savefig('case1_bridge_geometry.png', dpi=150, bbox_inches='tight')
print("\n桥梁几何图形已保存")
plt.close()
# 边界条件：固定两个塔底
fixed_dofs = [0, 1, 2, 6, 7, 8] # 节点0和2的三个自由度
# 模态分析
print("\n" + "=" * 70)
print("模态分析")
print("=" * 70)
frequencies, modes = solve_modal(K, M, fixed_dofs, n_modes=8)
if len(frequencies) > 0:
print(f"\n前{len(frequencies)}阶固有频率：")
for i in range(len(frequencies)):
period = 1 / frequencies[i]
print(f" 第{i+1}阶：f = {frequencies[i]:.4f} Hz，T = {period:.4f} s")
# 绘制频率分布
fig2, ax2 = plt.subplots(figsize=(10, 6))
mode_numbers = np.arange(1, len(frequencies) + 1)
ax2.bar(mode_numbers, frequencies, color='steelblue', edgecolor='black')
ax2.set_xlabel('模态阶数')
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('斜拉桥固有频率分布')
ax2.grid(True, alpha=0.3, axis='y')
for i, f in enumerate(frequencies):
ax2.text(i + 1, f + 0.01, f'{f:.3f}', ha='center', va='bottom', fontsize=9)
plt.tight_layout()
plt.savefig('case1_frequency_distribution.png', dpi=150, bbox_inches='tight')
print("\n频率分布图形已保存")
plt.close()
# 绘制模态振型
n_plot = min(6, len(frequencies))
n_rows = (n_plot + 1) // 2
fig3, axes = plt.subplots(n_rows, 2, figsize=(14, 4*n_rows))
if n_rows == 1:
axes = axes.reshape(1, -1)
for i in range(n_plot):
ax = axes[i // 2, i % 2]
mode = modes[:, i]
scale = 10
nodes_def = nodes.copy()
nodes_def[:, 0] += mode[::3] * scale
nodes_def[:, 1] += mode[1::3] * scale
for j, elem in enumerate(elements):
node_i, node_j = elem
x_orig = [nodes[node_i, 0], nodes[node_j, 0]]
y_orig = [nodes[node_i, 1], nodes[node_j, 1]]
x_def = [nodes_def[node_i, 0], nodes_def[node_j, 0]]
y_def = [nodes_def[node_i, 1], nodes_def[node_j, 1]]
ax.plot(x_orig, y_orig, 'k--', alpha=0.3, linewidth=1)
elem_type = element_types[j]
if elem_type == 'girder':
ax.plot(x_def, y_def, 'b-', linewidth=2)
elif elem_type == 'tower':
ax.plot(x_def, y_def, 'k-', linewidth=3)
elif elem_type == 'cable':
ax.plot(x_def, y_def, 'r-', linewidth=1, alpha=0.7)
ax.scatter(nodes[:, 0], nodes[:, 1], c='gray', s=30, alpha=0.5, zorder=5)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title(f'第{i+1}阶模态 (f = {frequencies[i]:.3f} Hz)')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
for i in range(n_plot, n_rows * 2):
axes[i // 2, i % 2].axis('off')
plt.tight_layout()
plt.savefig('case1_mode_shapes.png', dpi=150, bbox_inches='tight')
print("模态振型图形已保存")
plt.close()
else:
print("模态分析未获得有效结果")
print("\n" + "=" * 70)
print("案例1运行完成！")
print("=" * 70)
if __name__ == "__main__":
main()
# ============================================
# 生成GIF动画
# ============================================
def create_animation():
"""创建仿真结果动画"""
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10)
ax.set_ylim(-2, 2)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Response', fontsize=12)
ax.set_title('Dynamic Response Animation', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
line, = ax.plot([], [], 'b-', linewidth=2)
def init():
line.set_data([], [])
return line,
def update(frame):
x = np.linspace(0, 10, 100)
y = np.sin(x - frame * 0.2) * np.exp(-frame * 0.01)
line.set_data(x, y)
return line,
anim = FuncAnimation(fig, update, init_func=init, frames=50, interval=100, blit=True)
output_dir = os.path.dirname(os.path.abspath(__file__))
anim.save(f'{output_dir}/simulation_animation.gif', writer='pillow', fps=10)
print(f"动画已保存到: {output_dir}/simulation_animation.gif")
plt.close()
if __name__ == '__main__':
create_animation()
import matplotlib
matplotlib.use('Agg')
"""
案例2：风荷载模拟与风致振动
包含：
1. 脉动风场生成（谐波叠加法）
2. 风谱分析（Davenport谱、Kaimal谱）
3. 抖振响应分析
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.size'] = 10
plt.rcParams['axes.unicode_minus'] = False
def davenport_spectrum(f, U10, z=10):
"""
Davenport风速谱
参数:
f: 频率 (Hz)
U10: 10m高度处平均风速 (m/s)
z: 高度 (m)
返回:
S: 功率谱密度 (m²/s)
"""
# 粗糙度长度（开阔地带）
z0 = 0.03
# 摩擦速度
kappa = 0.4 # 卡门常数
u_star = U10 * kappa / np.log(z / z0)
# 归一化频率
x = 1200 * f / U10
# Davenport谱
S = 4 * u_star**2 * x**2 / (f * (1 + x**2)**(4/3))
return S
def kaimal_spectrum(f, U, z):
"""
Kaimal风速谱
参数:
f: 频率 (Hz)
U: 平均风速 (m/s)
z: 高度 (m)
返回:
S: 功率谱密度 (m²/s)
"""
# 归一化频率
f_star = f * z / U
# Kaimal谱
S = 200 * f * U**2 / z / (1 + 50 * f_star)**(5/3)
return S
def generate_wind_time_series(U_mean, z, duration, dt, spectrum_type='kaimal'):
"""
生成脉动风时程（谐波叠加法）
参数:
U_mean: 平均风速 (m/s)
z: 高度 (m)
duration: 模拟时长 (s)
dt: 时间步长 (s)
spectrum_type: 谱类型 ('davenport' 或 'kaimal')
返回:
t: 时间数组
U: 风速时程
u_fluct: 脉动风速时程
"""
t = np.arange(0, duration, dt)
n_points = len(t)
# 频率范围
f_max = 5 # 最大频率 (Hz)
f_min = 0.01 # 最小频率 (Hz)
n_freq = 500 # 频率点数
f = np.linspace(f_min, f_max, n_freq)
df = f[1] - f[0]
# 计算功率谱
if spectrum_type == 'davenport':
S = davenport_spectrum(f, U_mean, z)
else:
S = kaimal_spectrum(f, U_mean, z)
# 谐波叠加
u_fluct = np.zeros(n_points)
for i in range(n_freq):
phi = np.random.uniform(0, 2 * np.pi)
amplitude = np.sqrt(2 * S[i] * df)
u_fluct += amplitude * np.cos(2 * np.pi * f[i] * t + phi)
U = U_mean + u_fluct
return t, U, u_fluct
def coherence_function(f, delta_y, U_mean, C_y=10):
"""
横向相干函数
参数:
f: 频率 (Hz)
delta_y: 横向距离 (m)
U_mean: 平均风速 (m/s)
C_y: 衰减系数
返回:
coh: 相干系数
"""
coh = np.exp(-C_y * f * delta_y / U_mean)
return coh
def generate_spatial_wind_field(U_mean, z, y_positions, duration, dt):
"""
生成空间相关风场
参数:
U_mean: 平均风速 (m/s)
z: 高度 (m)
y_positions: 横向位置数组 (m)
duration: 模拟时长 (s)
dt: 时间步长 (s)
返回:
t: 时间数组
U_field: 风速场 (n_positions x n_time)
"""
n_positions = len(y_positions)
t = np.arange(0, duration, dt)
n_points = len(t)
# 频率范围
f_max = 5
f_min = 0.01
n_freq = 500
f = np.linspace(f_min, f_max, n_freq)
df = f[1] - f[0]
# 计算功率谱
S = kaimal_spectrum(f, U_mean, z)
# 生成空间相关风场
U_field = np.zeros((n_positions, n_points))
for i in range(n_freq):
# 生成随机相位
phi = np.random.uniform(0, 2 * np.pi, n_positions)
# 计算相干矩阵
C = np.zeros((n_positions, n_positions))
for j in range(n_positions):
for k in range(n_positions):
delta_y = abs(y_positions[j] - y_positions[k])
C[j, k] = coherence_function(f[i], delta_y, U_mean)
# Cholesky分解
try:
L = np.linalg.cholesky(C + np.eye(n_positions) * 1e-10)
except:
L = np.eye(n_positions)
# 生成相关随机变量
xi = np.cos(phi)
eta = L @ xi
amplitude = np.sqrt(2 * S[i] * df)
for j in range(n_positions):
U_field[j, :] += amplitude * eta[j] * np.cos(2 * np.pi * f[i] * t)
U_field += U_mean
return t, U_field
def buffeting_force(U, B, D, rho_air=1.225):
"""
计算抖振力（准定常理论）
参数:
U: 风速 (m/s)
B: 桥面宽度 (m)
D: 桥面高度 (m)
rho_air: 空气密度 (kg/m³)
返回:
F_L: 升力 (N/m)
F_D: 阻力 (N/m)
F_M: 力矩 (N·m/m)
"""
# 气动系数（简化）
C_L = 0.1 # 升力系数
C_D = 0.6 # 阻力系数
C_M = 0.02 # 力矩系数
q = 0.5 * rho_air * U**2 # 动压
F_L = q * B * C_L
F_D = q * D * C_D
F_M = q * B**2 * C_M
return F_L, F_D, F_M
def sdof_buffeting_response(m, c, k, F, dt):
"""
计算单自由度抖振响应（Newmark-beta法）
参数:
m: 质量
c: 阻尼
k: 刚度
F: 外力时程
dt: 时间步长
返回:
u: 位移响应
v: 速度响应
a: 加速度响应
"""
n_steps = len(F)
u = np.zeros(n_steps)
v = np.zeros(n_steps)
a = np.zeros(n_steps)
# Newmark-beta参数
beta = 0.25
gamma = 0.5
# 等效刚度
k_eff = k + gamma / (beta * dt) * c + 1 / (beta * dt**2) * m
for i in range(1, n_steps):
# 预测
u_pred = u[i-1] + dt * v[i-1] + 0.5 * dt**2 * a[i-1]
v_pred = v[i-1] + dt * a[i-1]
# 等效载荷
F_eff = F[i] + (1 / (beta * dt**2) * m + gamma / (beta * dt) * c) * u[i-1] + \
(1 / (beta * dt) * m + (gamma/beta - 1) * c) * v[i-1] + \
(0.5/beta - 1) * m * a[i-1]
# 求解
u[i] = F_eff / k_eff
a[i] = (u[i] - u[i-1]) / (beta * dt**2) - v[i-1] / (beta * dt) - (0.5/beta - 1) * a[i-1]
v[i] = v[i-1] + dt * ((1 - gamma) * a[i-1] + gamma * a[i])
return u, v, a
def main():
"""
主函数
"""
print("=" * 70)
print("案例2：风荷载模拟与风致振动")
print("=" * 70)
# 参数设置
U_mean = 20 # 平均风速 (m/s)
z = 50 # 高度 (m)
duration = 60 # 模拟时长 (s)
dt = 0.01 # 时间步长 (s)
print(f"\n风场参数：")
print(f" 平均风速：{U_mean} m/s")
print(f" 高度：{z} m")
print(f" 模拟时长：{duration} s")
# 1. 生成脉动风时程
print("\n" + "=" * 70)
print("1. 脉动风时程生成")
print("=" * 70)
t, U, u_fluct = generate_wind_time_series(U_mean, z, duration, dt, 'kaimal')
# 计算统计特性
U_std = np.std(u_fluct)
turbulence_intensity = U_std / U_mean * 100
print(f"\n脉动风统计特性：")
print(f" 平均风速：{np.mean(U):.3f} m/s")
print(f" 风速标准差：{U_std:.3f} m/s")
print(f" 湍流强度：{turbulence_intensity:.2f}%")
print(f" 最大风速：{np.max(U):.3f} m/s")
print(f" 最小风速：{np.min(U):.3f} m/s")
# 绘制风速时程
fig1, axes1 = plt.subplots(2, 1, figsize=(12, 8))
axes1[0].plot(t, U, 'b-', linewidth=0.8)
axes1[0].axhline(y=U_mean, color='r', linestyle='--', label='平均风速')
axes1[0].set_xlabel('时间 (s)')
axes1[0].set_ylabel('风速 (m/s)')
axes1[0].set_title('风速时程')
axes1[0].legend()
axes1[0].grid(True, alpha=0.3)
axes1[1].plot(t, u_fluct, 'g-', linewidth=0.8)
axes1[1].axhline(y=0, color='r', linestyle='--')
axes1[1].set_xlabel('时间 (s)')
axes1[1].set_ylabel('脉动风速 (m/s)')
axes1[1].set_title('脉动风速时程')
axes1[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case2_wind_time_series.png', dpi=150, bbox_inches='tight')
print("\n风速时程图形已保存")
plt.close()
# 2. 风谱分析
print("\n" + "=" * 70)
print("2. 风谱分析")
print("=" * 70)
# 频率范围
f = np.linspace(0.01, 5, 500)
# 计算理论谱
S_davenport = davenport_spectrum(f, U_mean, z)
S_kaimal = kaimal_spectrum(f, U_mean, z)
# 从时程估计谱（周期图法）
n_fft = 2**12
u_fft = np.fft.fft(u_fluct[:n_fft])
psd = np.abs(u_fft)**2 / n_fft * dt
freqs = np.fft.fftfreq(n_fft, dt)
# 只取正频率
positive_idx = freqs > 0
freqs_pos = freqs[positive_idx]
psd_pos = psd[positive_idx] * 2
# 绘制风谱
fig2, ax2 = plt.subplots(figsize=(10, 6))
ax2.loglog(f, S_davenport, 'b-', linewidth=2, label='Davenport谱')
ax2.loglog(f, S_kaimal, 'r-', linewidth=2, label='Kaimal谱')
ax2.loglog(freqs_pos[:len(freqs_pos)//2], psd_pos[:len(psd_pos)//2],
'g--', alpha=0.7, label='估计谱（周期图）')
ax2.set_xlabel('频率 (Hz)')
ax2.set_ylabel('功率谱密度 (m²/s)')
ax2.set_title('风速功率谱密度')
ax2.legend()
ax2.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('case2_wind_spectrum.png', dpi=150, bbox_inches='tight')
print("风谱图形已保存")
plt.close()
# 3. 空间风场
print("\n" + "=" * 70)
print("3. 空间相关风场")
print("=" * 70)
y_positions = np.linspace(-20, 20, 5) # 横向位置
t_spatial, U_field = generate_spatial_wind_field(U_mean, z, y_positions, duration, dt)
print(f"\n空间风场信息：")
print(f" 横向位置数：{len(y_positions)}")
print(f" 横向位置：{y_positions} m")
# 绘制空间风场
fig3, ax3 = plt.subplots(figsize=(12, 6))
for i, y in enumerate(y_positions):
ax3.plot(t_spatial[:1000], U_field[i, :1000],
label=f'y = {y:.0f} m', alpha=0.8)
ax3.set_xlabel('时间 (s)')
ax3.set_ylabel('风速 (m/s)')
ax3.set_title('空间相关风速时程（前10秒）')
ax3.legend()
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case2_spatial_wind_field.png', dpi=150, bbox_inches='tight')
print("空间风场图形已保存")
plt.close()
# 4. 抖振响应分析
print("\n" + "=" * 70)
print("4. 抖振响应分析")
print("=" * 70)
# 桥梁参数
B = 20 # 桥面宽度 (m)
D = 3 # 桥面高度 (m)
m = 20000 # 单位长度质量 (kg/m)
L = 100 # 跨度 (m)
# 结构动力参数
f_n = 0.5 # 固有频率 (Hz)
omega_n = 2 * np.pi * f_n
k = m * omega_n**2
xi = 0.02 # 阻尼比
c = 2 * xi * m * omega_n
print(f"\n结构参数：")
print(f" 桥面宽度：{B} m")
print(f" 桥面高度：{D} m")
print(f" 单位长度质量：{m} kg/m")
print(f" 固有频率：{f_n} Hz")
print(f" 阻尼比：{xi*100}%")
# 计算抖振力
F_L, F_D, F_M = buffeting_force(U, B, D)
# 简化为竖向力（升力）
F_buffeting = F_L * L # 总力
# 计算响应
u, v, a = sdof_buffeting_response(m*L, c*L, k*L, F_buffeting, dt)
# 响应统计
u_max = np.max(np.abs(u))
u_rms = np.std(u)
print(f"\n抖振响应统计：")
print(f" 最大位移：{u_max:.4f} m")
print(f" 位移均方根：{u_rms:.4f} m")
print(f" 最大速度：{np.max(np.abs(v)):.4f} m/s")
print(f" 最大加速度：{np.max(np.abs(a)):.4f} m/s²")
# 绘制响应
fig4, axes4 = plt.subplots(3, 1, figsize=(12, 10))
axes4[0].plot(t, u * 1000, 'b-', linewidth=0.8) # 转换为mm
axes4[0].set_ylabel('位移 (mm)')
axes4[0].set_title('竖向位移响应')
axes4[0].grid(True, alpha=0.3)
axes4[1].plot(t, v, 'g-', linewidth=0.8)
axes4[1].set_ylabel('速度 (m/s)')
axes4[1].set_title('竖向速度响应')
axes4[1].grid(True, alpha=0.3)
axes4[2].plot(t, a, 'r-', linewidth=0.8)
axes4[2].set_xlabel('时间 (s)')
axes4[2].set_ylabel('加速度 (m/s²)')
axes4[2].set_title('竖向加速度响应')
axes4[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case2_buffeting_response.png', dpi=150, bbox_inches='tight')
print("抖振响应图形已保存")
plt.close()
# 5. 响应谱分析
print("\n" + "=" * 70)
print("5. 响应谱分析")
print("=" * 70)
# 计算位移响应谱
n_fft = 2**12
u_fft = np.fft.fft(u[:n_fft])
S_u = np.abs(u_fft)**2 / n_fft * dt
freqs = np.fft.fftfreq(n_fft, dt)
positive_idx = freqs > 0
freqs_pos = freqs[positive_idx]
S_u_pos = S_u[positive_idx] * 2
# 绘制响应谱
fig5, ax5 = plt.subplots(figsize=(10, 6))
ax5.semilogy(freqs_pos[:len(freqs_pos)//2], S_u_pos[:len(S_u_pos)//2] * 1e6,
'b-', linewidth=1.5)
ax5.axvline(x=f_n, color='r', linestyle='--', label=f'固有频率 {f_n} Hz')
ax5.set_xlabel('频率 (Hz)')
ax5.set_ylabel('位移功率谱密度 (mm²/Hz)')
ax5.set_title('位移响应功率谱密度')
ax5.legend()
ax5.grid(True, alpha=0.3)
ax5.set_xlim([0, 2])
plt.tight_layout()
plt.savefig('case2_response_spectrum.png', dpi=150, bbox_inches='tight')
print("响应谱图形已保存")
plt.close()
print("\n" + "=" * 70)
print("案例2运行完成！")
print("=" * 70)
if __name__ == "__main__":
main()
# ============================================
# 生成GIF动画
# ============================================
def create_animation():
"""创建仿真结果动画"""
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10)
ax.set_ylim(-2, 2)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Response', fontsize=12)
ax.set_title('Dynamic Response Animation', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
line, = ax.plot([], [], 'b-', linewidth=2)
def init():
line.set_data([], [])
return line,
def update(frame):
x = np.linspace(0, 10, 100)
y = np.sin(x - frame * 0.2) * np.exp(-frame * 0.01)
line.set_data(x, y)
return line,
anim = FuncAnimation(fig, update, init_func=init, frames=50, interval=100, blit=True)
output_dir = os.path.dirname(os.path.abspath(__file__))
anim.save(f'{output_dir}/simulation_animation.gif', writer='pillow', fps=10)
print(f"动画已保存到: {output_dir}/simulation_animation.gif")
plt.close()
if __name__ == '__main__':
create_animation()
import matplotlib
matplotlib.use('Agg')
"""
案例3：颤振与涡激振动分析
包含：
1. 颤振稳定性分析
2. 涡激振动分析
3. 临界风速计算
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.size'] = 10
plt.rcParams['axes.unicode_minus'] = False
def scanlan_flutter_derivatives(U, B, omega):
"""
Scanlan颤振导数模型（简化）
参数:
U: 风速 (m/s)
B: 桥面宽度 (m)
omega: 圆频率 (rad/s)
返回:
H1, H2, H3, A1, A2, A3: 颤振导数
"""
# 无量纲风速
K = omega * B / U
# 简化颤振导数（典型值）
H1 = -0.5 / (1 + 0.5*K**2) # 竖向阻尼
H2 = 0.1 * K / (1 + K**2) # 竖向刚度
H3 = 0.05 / (1 + K**2) # 竖向耦合
A1 = -0.3 / (1 + 0.3*K**2) # 扭转阻尼
A2 = 0.2 * K / (1 + K**2) # 扭转刚度
A3 = 0.1 / (1 + K**2) # 扭转耦合
return H1, H2, H3, A1, A2, A3
def flutter_analysis(m, I_alpha, xi_h, xi_alpha, omega_h, omega_alpha, B, rho, U_range):
"""
颤振分析
参数:
m: 单位长度质量 (kg/m)
I_alpha: 单位长度质量惯性矩 (kg·m²/m)
xi_h: 竖向阻尼比
xi_alpha: 扭转阻尼比
omega_h: 竖向圆频率 (rad/s)
omega_alpha: 扭转圆频率 (rad/s)
B: 桥面宽度 (m)
rho: 空气密度 (kg/m³)
U_range: 风速范围 (m/s)
返回:
results: 分析结果字典
"""
n_U = len(U_range)
# 存储结果
frequencies = np.zeros((n_U, 2))
damping_ratios = np.zeros((n_U, 2))
for i, U in enumerate(U_range):
# 计算颤振导数
H1, H2, H3, A1, A2, A3 = scanlan_flutter_derivatives(U, B, omega_h)
# 气动阻尼和刚度
c_h_aero = 0.5 * rho * U * B * H1
c_alpha_aero = 0.5 * rho * U * B**3 * A1
k_h_aero = 0.5 * rho * U**2 * H2
k_alpha_aero = 0.5 * rho * U**2 * B**2 * A2
# 总阻尼和刚度
c_h = 2 * xi_h * m * omega_h + c_h_aero
c_alpha = 2 * xi_alpha * I_alpha * omega_alpha + c_alpha_aero
k_h = m * omega_h**2 + k_h_aero
k_alpha = I_alpha * omega_alpha**2 + k_alpha_aero
# 计算等效频率和阻尼比
omega_h_eff = np.sqrt(k_h / m)
omega_alpha_eff = np.sqrt(k_alpha / I_alpha)
xi_h_eff = c_h / (2 * m * omega_h_eff)
xi_alpha_eff = c_alpha / (2 * I_alpha * omega_alpha_eff)
frequencies[i, 0] = omega_h_eff / (2 * np.pi)
frequencies[i, 1] = omega_alpha_eff / (2 * np.pi)
damping_ratios[i, 0] = xi_h_eff
damping_ratios[i, 1] = xi_alpha_eff
return {
'U': U_range,
'frequencies': frequencies,
'damping_ratios': damping_ratios
}
def find_critical_flutter_speed(m, I_alpha, xi_h, xi_alpha, omega_h, omega_alpha, B, rho, U_max=100):
"""
计算临界颤振风速
参数:
m, I_alpha, xi_h, xi_alpha, omega_h, omega_alpha, B, rho: 结构参数
U_max: 最大搜索风速
返回:
U_cr: 临界颤振风速
"""
U_range = np.linspace(1, U_max, 200)
for U in U_range:
H1, H2, H3, A1, A2, A3 = scanlan_flutter_derivatives(U, B, omega_h)
# 气动阻尼
c_h_aero = 0.5 * rho * U * B * H1
c_alpha_aero = 0.5 * rho * U * B**3 * A1
# 总阻尼
c_h = 2 * xi_h * m * omega_h + c_h_aero
c_alpha = 2 * xi_alpha * I_alpha * omega_alpha + c_alpha_aero
# 检查是否失稳
if c_h < 0 or c_alpha < 0:
return U
return U_max
def vortex_shedding_frequency(U, D, St=0.2):
"""
计算涡激脱落频率
参数:
U: 风速 (m/s)
D: 特征尺寸 (m)
St: 斯特劳哈尔数
返回:
f_v: 涡激脱落频率 (Hz)
"""
f_v = St * U / D
return f_v
def viv_amplitude(U, D, f_n, xi, St=0.2, C_L0=0.5, rho_air=1.225):
"""
计算涡激振动振幅（简化模型）
参数:
U: 风速 (m/s)
D: 特征尺寸 (m)
f_n: 结构固有频率 (Hz)
xi: 结构阻尼比
St: 斯特劳哈尔数
C_L0: 升力系数幅值
rho_air: 空气密度
返回:
A: 振幅 (m)
"""
# 涡激脱落频率
f_v = vortex_shedging_frequency(U, D, St)
# 频率比
ratio = f_v / f_n
# 锁定区域
if 0.9 <= ratio <= 1.1:
# 锁定区域内的振幅
m_star = 1000 # 无量纲质量
A_star = C_L0 / (4 * np.pi * St**2 * xi * m_star)
A = A_star * D
else:
# 锁定区域外的振幅（很小）
A = 0.01 * D * np.exp(-5 * abs(ratio - 1))
return A
def viv_response(U, D, m, c, k, St=0.2, C_L0=0.5, rho_air=1.225):
"""
计算涡激振动响应
参数:
U: 风速 (m/s)
D: 特征尺寸 (m)
m: 质量
c: 阻尼
k: 刚度
St: 斯特劳哈尔数
C_L0: 升力系数
rho_air: 空气密度
返回:
A: 振幅
f_v: 涡激频率
"""
omega_n = np.sqrt(k / m)
f_n = omega_n / (2 * np.pi)
xi = c / (2 * m * omega_n)
f_v = vortex_shedding_frequency(U, D, St)
# 升力幅值
q = 0.5 * rho_air * U**2
F_L0 = q * D * C_L0
# 频率比
ratio = f_v / f_n
# 响应幅值（考虑锁定效应）
if 0.9 <= ratio <= 1.1:
# 锁定区域 - 大幅振动
H = 1 / (2 * xi) # 放大因子
A = F_L0 / k * H
else:
# 非锁定区域 - 小幅振动
H = 1 / np.sqrt((1 - ratio**2)**2 + (2 * xi * ratio)**2)
A = F_L0 / k * H * 0.1 # 衰减因子
return A, f_v
def simulate_viv(U, D, m, c, k, duration=50, dt=0.01, St=0.2, C_L0=0.5):
"""
模拟涡激振动时程
参数:
U: 风速 (m/s)
D: 特征尺寸 (m)
m: 质量
c: 阻尼
k: 刚度
duration: 模拟时长
dt: 时间步长
St: 斯特劳哈尔数
C_L0: 升力系数
返回:
t: 时间数组
u: 位移响应
"""
t = np.arange(0, duration, dt)
n_steps = len(t)
omega_n = np.sqrt(k / m)
f_n = omega_n / (2 * np.pi)
f_v = vortex_shedding_frequency(U, D, St)
# 涡激力
q = 0.5 * 1.225 * U**2
F_L0 = q * D * C_L0
F_viv = F_L0 * np.sin(2 * np.pi * f_v * t)
# 锁定效应调制
ratio = f_v / f_n
if 0.9 <= ratio <= 1.1:
lock_in = 1.0
else:
lock_in = 0.1 * np.exp(-5 * abs(ratio - 1))
F_viv *= lock_in
# Newmark-beta求解
u = np.zeros(n_steps)
v = np.zeros(n_steps)
a = np.zeros(n_steps)
beta = 0.25
gamma = 0.5
k_eff = k + gamma / (beta * dt) * c + 1 / (beta * dt**2) * m
for i in range(1, n_steps):
F_eff = F_viv[i] + (1 / (beta * dt**2) * m + gamma / (beta * dt) * c) * u[i-1] + \
(1 / (beta * dt) * m + (gamma/beta - 1) * c) * v[i-1] + \
(0.5/beta - 1) * m * a[i-1]
u[i] = F_eff / k_eff
a[i] = (u[i] - u[i-1]) / (beta * dt**2) - v[i-1] / (beta * dt) - (0.5/beta - 1) * a[i-1]
v[i] = v[i-1] + dt * ((1 - gamma) * a[i-1] + gamma * a[i])
return t, u
def main():
"""
主函数
"""
print("=" * 70)
print("案例3：颤振与涡激振动分析")
print("=" * 70)
# 结构参数
m = 20000 # 单位长度质量 (kg/m)
I_alpha = 500000 # 单位长度质量惯性矩 (kg·m²/m)
B = 20 # 桥面宽度 (m)
D = 3 # 桥面高度 (m)
f_h = 0.3 # 竖向频率 (Hz)
f_alpha = 0.6 # 扭转频率 (Hz)
omega_h = 2 * np.pi * f_h
omega_alpha = 2 * np.pi * f_alpha
xi_h = 0.01 # 竖向阻尼比
xi_alpha = 0.01 # 扭转阻尼比
rho = 1.225 # 空气密度 (kg/m³)
print(f"\n结构参数：")
print(f" 单位长度质量：{m} kg/m")
print(f" 质量惯性矩：{I_alpha} kg·m²/m")
print(f" 桥面宽度：{B} m")
print(f" 桥面高度：{D} m")
print(f" 竖向频率：{f_h} Hz")
print(f" 扭转频率：{f_alpha} Hz")
print(f" 竖向阻尼比：{xi_h*100}%")
print(f" 扭转阻尼比：{xi_alpha*100}%")
# 1. 颤振分析
print("\n" + "=" * 70)
print("1. 颤振稳定性分析")
print("=" * 70)
U_range = np.linspace(5, 80, 100)
results = flutter_analysis(m, I_alpha, xi_h, xi_alpha, omega_h, omega_alpha, B, rho, U_range)
# 绘制颤振频率曲线
fig1, axes1 = plt.subplots(2, 1, figsize=(10, 8))
axes1[0].plot(results['U'], results['frequencies'][:, 0], 'b-', linewidth=2, label='竖向')
axes1[0].plot(results['U'], results['frequencies'][:, 1], 'r-', linewidth=2, label='扭转')
axes1[0].set_xlabel('风速 (m/s)')
axes1[0].set_ylabel('频率 (Hz)')
axes1[0].set_title('颤振频率曲线')
axes1[0].legend()
axes1[0].grid(True, alpha=0.3)
axes1[1].plot(results['U'], results['damping_ratios'][:, 0], 'b-', linewidth=2, label='竖向')
axes1[1].plot(results['U'], results['damping_ratios'][:, 1], 'r-', linewidth=2, label='扭转')
axes1[1].axhline(y=0, color='k', linestyle='--', alpha=0.5)
axes1[1].set_xlabel('风速 (m/s)')
axes1[1].set_ylabel('阻尼比')
axes1[1].set_title('颤振阻尼曲线')
axes1[1].legend()
axes1[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case3_flutter_analysis.png', dpi=150, bbox_inches='tight')
print("\n颤振分析图形已保存")
plt.close()
# 计算临界颤振风速
U_cr = find_critical_flutter_speed(m, I_alpha, xi_h, xi_alpha, omega_h, omega_alpha, B, rho)
print(f"\n临界颤振风速：{U_cr:.2f} m/s")
# 2. 涡激振动分析
print("\n" + "=" * 70)
print("2. 涡激振动分析")
print("=" * 70)
# 计算涡激脱落频率
U_viv_range = np.linspace(1, 30, 100)
f_v_range = vortex_shedding_frequency(U_viv_range, D)
print(f"\n涡激振动参数：")
print(f" 斯特劳哈尔数：0.2")
print(f" 特征尺寸：{D} m")
print(f" 结构固有频率：{f_h} Hz")
# 锁定风速
St = 0.2
U_lock_in = f_h * D / St
print(f" 锁定风速：{U_lock_in:.2f} m/s")
# 绘制涡激频率曲线
fig2, ax2 = plt.subplots(figsize=(10, 6))
ax2.plot(U_viv_range, f_v_range, 'b-', linewidth=2, label='涡激脱落频率')
ax2.axhline(y=f_h, color='r', linestyle='--', label=f'结构频率 {f_h} Hz')
ax2.axvline(x=U_lock_in, color='g', linestyle=':', label=f'锁定风速 {U_lock_in:.1f} m/s')
# 标记锁定区域
ax2.axvspan(U_lock_in*0.9, U_lock_in*1.1, alpha=0.2, color='yellow', label='锁定区域')
ax2.set_xlabel('风速 (m/s)')
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('涡激脱落频率与结构频率关系')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim([0, 30])
ax2.set_ylim([0, 2])
plt.tight_layout()
plt.savefig('case3_vortex_shedding.png', dpi=150, bbox_inches='tight')
print("\n涡激脱落频率图形已保存")
plt.close()
# 3. VIV振幅分析
print("\n" + "=" * 70)
print("3. VIV振幅分析")
print("=" * 70)
# 计算不同风速下的VIV振幅
A_viv = np.zeros(len(U_viv_range))
for i, U in enumerate(U_viv_range):
A_viv[i], _ = viv_response(U, D, m, 2*xi_h*m*omega_h, m*omega_h**2)
# 绘制VIV振幅曲线
fig3, ax3 = plt.subplots(figsize=(10, 6))
ax3.plot(U_viv_range, A_viv * 1000, 'b-', linewidth=2)
ax3.axvline(x=U_lock_in, color='r', linestyle='--', label=f'锁定风速 {U_lock_in:.1f} m/s')
ax3.axvspan(U_lock_in*0.9, U_lock_in*1.1, alpha=0.2, color='yellow', label='锁定区域')
ax3.set_xlabel('风速 (m/s)')
ax3.set_ylabel('振幅 (mm)')
ax3.set_title('涡激振动振幅随风速变化')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim([0, 30])
plt.tight_layout()
plt.savefig('case3_viv_amplitude.png', dpi=150, bbox_inches='tight')
print("VIV振幅图形已保存")
plt.close()
# 最大振幅
A_max = np.max(A_viv)
U_max_viv = U_viv_range[np.argmax(A_viv)]
print(f"\n最大VIV振幅：{A_max*1000:.2f} mm")
print(f"最大振幅对应风速：{U_max_viv:.2f} m/s")
# 4. VIV时程模拟
print("\n" + "=" * 70)
print("4. VIV时程模拟")
print("=" * 70)
# 在锁定风速下模拟
t, u_viv = simulate_viv(U_lock_in, D, m, 2*xi_h*m*omega_h, m*omega_h**2, duration=30)
# 在非锁定风速下模拟
t, u_non_viv = simulate_viv(U_lock_in*0.7, D, m, 2*xi_h*m*omega_h, m*omega_h**2, duration=30)
print(f"\n时程模拟结果：")
print(f" 锁定风速下最大振幅：{np.max(np.abs(u_viv))*1000:.2f} mm")
print(f" 非锁定风速下最大振幅：{np.max(np.abs(u_non_viv))*1000:.2f} mm")
# 绘制时程
fig4, axes4 = plt.subplots(2, 1, figsize=(12, 8))
axes4[0].plot(t, u_viv * 1000, 'b-', linewidth=0.8)
axes4[0].set_ylabel('位移 (mm)')
axes4[0].set_title(f'涡激振动时程（锁定风速 {U_lock_in:.1f} m/s）')
axes4[0].grid(True, alpha=0.3)
axes4[1].plot(t, u_non_viv * 1000, 'g-', linewidth=0.8)
axes4[1].set_xlabel('时间 (s)')
axes4[1].set_ylabel('位移 (mm)')
axes4[1].set_title(f'非锁定风速时程（风速 {U_lock_in*0.7:.1f} m/s）')
axes4[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case3_viv_time_history.png', dpi=150, bbox_inches='tight')
print("VIV时程图形已保存")
plt.close()
# 5. 综合对比
print("\n" + "=" * 70)
print("5. 风致振动综合对比")
print("=" * 70)
fig5, ax5 = plt.subplots(figsize=(12, 6))
# 颤振边界
ax5.axvline(x=U_cr, color='r', linestyle='--', linewidth=2, label=f'颤振临界风速 {U_cr:.1f} m/s')
# VIV区域
ax5.axvspan(U_lock_in*0.9, U_lock_in*1.1, alpha=0.3, color='yellow', label='VIV锁定区域')
# 抖振响应（简化）
U_buffeting = np.linspace(5, 60, 100)
sigma_buffeting = 0.001 * (U_buffeting / 20)**2 # 简化模型
ax5.plot(U_buffeting, sigma_buffeting * 1000, 'b-', linewidth=2, label='抖振响应')
# VIV响应
ax5.plot(U_viv_range, A_viv * 1000, 'g-', linewidth=2, label='VIV响应')
ax5.set_xlabel('风速 (m/s)')
ax5.set_ylabel('振幅 (mm)')
ax5.set_title('斜拉桥风致振动综合对比')
ax5.legend()
ax5.grid(True, alpha=0.3)
ax5.set_xlim([0, 60])
plt.tight_layout()
plt.savefig('case3_wind_induced_vibration_summary.png', dpi=150, bbox_inches='tight')
print("\n综合对比图形已保存")
plt.close()
print("\n" + "=" * 70)
print("案例3运行完成！")
print("=" * 70)
if __name__ == "__main__":
main()
# ============================================
# 生成GIF动画
# ============================================
def create_animation():
"""创建仿真结果动画"""
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10)
ax.set_ylim(-2, 2)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Response', fontsize=12)
ax.set_title('Dynamic Response Animation', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
line, = ax.plot([], [], 'b-', linewidth=2)
def init():
line.set_data([], [])
return line,
def update(frame):
x = np.linspace(0, 10, 100)
y = np.sin(x - frame * 0.2) * np.exp(-frame * 0.01)
line.set_data(x, y)
return line,
anim = FuncAnimation(fig, update, init_func=init, frames=50, interval=100, blit=True)
output_dir = os.path.dirname(os.path.abspath(__file__))
anim.save(f'{output_dir}/simulation_animation.gif', writer='pillow', fps=10)
print(f"动画已保存到: {output_dir}/simulation_animation.gif")
plt.close()
if __name__ == '__main__':
create_animation()
import matplotlib
matplotlib.use('Agg')
"""
案例4：斜拉桥地震响应分析
包含：
1. 地震动生成（人工地震波）
2. 时程分析
3. 响应谱分析
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.size'] = 10
plt.rcParams['axes.unicode_minus'] = False
def kanai_tajimi_spectrum(omega, omega_g, xi_g, S0):
"""
Kanai-Tajimi功率谱
参数:
omega: 圆频率 (rad/s)
omega_g: 场地特征频率 (rad/s)
xi_g: 场地阻尼比
S0: 谱强度
返回:
S: 功率谱密度
"""
S = S0 * (1 + (2 * xi_g * omega / omega_g)**2) / \
((1 - (omega / omega_g)**2)**2 + (2 * xi_g * omega / omega_g)**2)
return S
def generate_artificial_earthquake(duration, dt, omega_g, xi_g, S0, seed=None):
"""
生成人工地震波（谐波叠加法）
参数:
duration: 持续时间 (s)
dt: 时间步长 (s)
omega_g: 场地特征频率 (rad/s)
xi_g: 场地阻尼比
S0: 谱强度
seed: 随机种子
返回:
t: 时间数组
acc: 加速度时程 (m/s²)
"""
if seed is not None:
np.random.seed(seed)
t = np.arange(0, duration, dt)
n_points = len(t)
# 频率范围
omega_max = 50 # 最大频率 (rad/s)
omega_min = 0.1 # 最小频率 (rad/s)
n_freq = 500
omega = np.linspace(omega_min, omega_max, n_freq)
d_omega = omega[1] - omega[0]
# 计算功率谱
S = kanai_tajimi_spectrum(omega, omega_g, xi_g, S0)
# 包络函数（模拟地震动强度变化）
t_rise = 2 # 上升段时间
t_flat = 8 # 平稳段时间
t_decay = 5 # 衰减段时间
envelope = np.ones_like(t)
for i, ti in enumerate(t):
if ti < t_rise:
envelope[i] = (ti / t_rise)**2
elif ti < t_rise + t_flat:
envelope[i] = 1.0
elif ti < t_rise + t_flat + t_decay:
envelope[i] = np.exp(-0.5 * (ti - t_rise - t_flat))
else:
envelope[i] = 0.0
# 谐波叠加
acc = np.zeros(n_points)
for i in range(n_freq):
phi = np.random.uniform(0, 2 * np.pi)
amplitude = np.sqrt(2 * S[i] * d_omega)
acc += amplitude * np.cos(omega[i] * t + phi)
# 应用包络函数
acc *= envelope
return t, acc
def generate_el_centro_record(duration=30, dt=0.02):
"""
生成简化的El Centro地震波
参数:
duration: 持续时间
dt: 时间步长
返回:
t: 时间数组
acc: 加速度时程 (m/s²)
"""
t = np.arange(0, duration, dt)
n_points = len(t)
# 简化的El Centro特征
acc = np.zeros(n_points)
# 主要频率成分
freqs = [1.5, 3.0, 5.0, 8.0]
amps = [3.0, 2.0, 1.0, 0.5]
phases = [0, np.pi/4, np.pi/2, np.pi]
for f, a, p in zip(freqs, amps, phases):
acc += a * np.sin(2 * np.pi * f * t + p)
# 添加随机成分
acc += 0.5 * np.random.randn(n_points)
# 包络函数
envelope = np.exp(-t / 5) * (1 - np.exp(-t / 1))
acc *= envelope
# 归一化到最大加速度约3.4 m/s² (0.35g)
acc = acc / np.max(np.abs(acc)) * 3.4
return t, acc
def newmark_beta(M, C, K, F, dt, beta=0.25, gamma=0.5):
"""
Newmark-beta法求解动力方程
参数:
M: 质量矩阵
C: 阻尼矩阵
K: 刚度矩阵
F: 外力时程 (n_dof x n_steps)
dt: 时间步长
beta, gamma: Newmark参数
返回:
u: 位移响应
v: 速度响应
a: 加速度响应
"""
n_dof = len(M)
n_steps = F.shape[1]
u = np.zeros((n_dof, n_steps))
v = np.zeros((n_dof, n_steps))
a = np.zeros((n_dof, n_steps))
# 等效刚度
K_eff = K + gamma / (beta * dt) * C + 1 / (beta * dt**2) * M
for i in range(1, n_steps):
# 等效载荷
F_eff = F[:, i] + \
(1 / (beta * dt**2) * M + gamma / (beta * dt) * C) @ u[:, i-1] + \
(1 / (beta * dt) * M + (gamma/beta - 1) * C) @ v[:, i-1] + \
(0.5/beta - 1) * M @ a[:, i-1]
# 求解
u[:, i] = np.linalg.solve(K_eff, F_eff)
# 更新加速度和速度
a[:, i] = (u[:, i] - u[:, i-1]) / (beta * dt**2) - \
v[:, i-1] / (beta * dt) - (0.5/beta - 1) * a[:, i-1]
v[:, i] = v[:, i-1] + dt * ((1 - gamma) * a[:, i-1] + gamma * a[:, i])
return u, v, a
def create_simplified_bridge_model():
"""
创建简化桥梁模型（单塔或主跨简化）
返回:
M, C, K: 质量、阻尼、刚度矩阵
phi: 振型矩阵
omega_n: 固有频率
"""
# 简化模型：3自由度系统
# 模拟主跨的竖向、横向和扭转振动
m = 20000 # 等效质量 (kg)
I = 500000 # 等效转动惯量 (kg·m²)
M = np.diag([m, m, I])
# 刚度
k_v = 1.0e7 # 竖向刚度 (N/m)
k_h = 8.0e6 # 横向刚度 (N/m)
k_t = 5.0e8 # 扭转刚度 (N·m/rad)
K = np.diag([k_v, k_h, k_t])
# 阻尼（瑞利阻尼）
omega_1 = np.sqrt(k_v / m)
omega_2 = np.sqrt(k_h / m)
xi = 0.02 # 阻尼比
a0 = xi * 2 * omega_1 * omega_2 / (omega_1 + omega_2)
a1 = xi * 2 / (omega_1 + omega_2)
C = a0 * M + a1 * K
# 计算模态
eigenvalues, eigenvectors = np.linalg.eig(np.linalg.solve(M, K))
idx = np.argsort(eigenvalues)
omega_n = np.sqrt(eigenvalues[idx])
phi = eigenvectors[:, idx]
return M, C, K, phi, omega_n
def response_spectrum(omega, xi, acc, dt):
"""
计算反应谱
参数:
omega: 固有频率数组 (rad/s)
xi: 阻尼比
acc: 加速度时程
dt: 时间步长
返回:
Sd: 位移谱
Sv: 速度谱
Sa: 加速度谱
"""
n_omega = len(omega)
Sd = np.zeros(n_omega)
Sv = np.zeros(n_omega)
Sa = np.zeros(n_omega)
for i, w in enumerate(omega):
m = 1.0
c = 2 * xi * w * m
k = w**2 * m
# Newmark-beta求解
F = -m * acc
u, v, a = newmark_beta(np.array([[m]]), np.array([[c]]), np.array([[k]]),
F.reshape(1, -1), dt)
Sd[i] = np.max(np.abs(u))
Sv[i] = np.max(np.abs(v))
Sa[i] = np.max(np.abs(a + acc))
return Sd, Sv, Sa
def modal_analysis_seismic(M, C, K, phi, omega_n, acc, dt):
"""
模态叠加法地震响应分析
参数:
M, C, K: 质量、阻尼、刚度矩阵
phi: 振型矩阵
omega_n: 固有频率
acc: 地震加速度时程
dt: 时间步长
返回:
u: 位移响应
"""
n_dof = len(M)
n_modes = len(omega_n)
n_steps = len(acc)
# 模态坐标响应
q = np.zeros((n_modes, n_steps))
for i in range(n_modes):
# 模态质量
M_i = phi[:, i] @ M @ phi[:, i]
# 参与系数
gamma_i = -phi[:, i] @ M @ np.ones(n_dof) / M_i
# 模态阻尼
xi_i = (phi[:, i] @ C @ phi[:, i]) / (2 * omega_n[i] * M_i)
# 单自由度响应
m = 1.0
c = 2 * xi_i * omega_n[i] * m
k = omega_n[i]**2 * m
F = gamma_i * acc
u_i, _, _ = newmark_beta(np.array([[m]]), np.array([[c]]), np.array([[k]]),
F.reshape(1, -1), dt)
q[i, :] = u_i[0, :]
# 叠加模态响应
u = phi @ q
return u
def main():
"""
主函数
"""
print("=" * 70)
print("案例4：斜拉桥地震响应分析")
print("=" * 70)
# 创建桥梁模型
print("\n创建桥梁模型...")
M, C, K, phi, omega_n = create_simplified_bridge_model()
f_n = omega_n / (2 * np.pi)
print(f"\n桥梁固有频率：")
for i in range(len(f_n)):
print(f" 第{i+1}阶：f = {f_n[i]:.3f} Hz，T = {1/f_n[i]:.3f} s")
# 1. 生成地震动
print("\n" + "=" * 70)
print("1. 地震动生成")
print("=" * 70)
duration = 20 # 持续时间
dt = 0.01 # 时间步长
# 生成人工地震波
omega_g = 15 # 硬土场地
xi_g = 0.6
S0 = 0.01
t, acc_artificial = generate_artificial_earthquake(duration, dt, omega_g, xi_g, S0, seed=42)
# 生成简化El Centro波
t_el, acc_el = generate_el_centro_record(duration, dt)
# PGA
pga_artificial = np.max(np.abs(acc_artificial))
pga_el = np.max(np.abs(acc_el))
print(f"\n人工地震波：")
print(f" PGA = {pga_artificial:.3f} m/s² ({pga_artificial/9.81*100:.1f}%g)")
print(f"\n简化El Centro波：")
print(f" PGA = {pga_el:.3f} m/s² ({pga_el/9.81*100:.1f}%g)")
# 绘制地震波
fig1, axes1 = plt.subplots(2, 1, figsize=(12, 8))
axes1[0].plot(t, acc_artificial, 'b-', linewidth=0.8)
axes1[0].set_ylabel('加速度 (m/s²)')
axes1[0].set_title(f'人工地震波 (PGA = {pga_artificial:.2f} m/s²)')
axes1[0].grid(True, alpha=0.3)
axes1[1].plot(t_el, acc_el, 'r-', linewidth=0.8)
axes1[1].set_xlabel('时间 (s)')
axes1[1].set_ylabel('加速度 (m/s²)')
axes1[1].set_title(f'简化El Centro波 (PGA = {pga_el:.2f} m/s²)')
axes1[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case4_earthquake_records.png', dpi=150, bbox_inches='tight')
print("\n地震波图形已保存")
plt.close()
# 2. 反应谱分析
print("\n" + "=" * 70)
print("2. 反应谱分析")
print("=" * 70)
# 周期范围
T_range = np.linspace(0.1, 5, 100)
omega_range = 2 * np.pi / T_range
xi_values = [0.02, 0.05, 0.10]
fig2, axes2 = plt.subplots(1, 2, figsize=(14, 5))
for xi in xi_values:
Sd, Sv, Sa = response_spectrum(omega_range, xi, acc_el, dt)
axes2[0].plot(T_range, Sa / 9.81, label=f'ξ = {xi*100:.0f}%')
axes2[1].plot(T_range, Sd * 1000, label=f'ξ = {xi*100:.0f}%')
# 标记桥梁周期
for i, f in enumerate(f_n):
T = 1 / f
axes2[0].axvline(x=T, color='k', linestyle='--', alpha=0.3)
axes2[1].axvline(x=T, color='k', linestyle='--', alpha=0.3)
axes2[0].set_xlabel('周期 (s)')
axes2[0].set_ylabel('加速度谱 (g)')
axes2[0].set_title('加速度反应谱')
axes2[0].legend()
axes2[0].grid(True, alpha=0.3)
axes2[1].set_xlabel('周期 (s)')
axes2[1].set_ylabel('位移谱 (mm)')
axes2[1].set_title('位移反应谱')
axes2[1].legend()
axes2[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case4_response_spectrum.png', dpi=150, bbox_inches='tight')
print("反应谱图形已保存")
plt.close()
# 3. 时程分析
print("\n" + "=" * 70)
print("3. 地震时程分析")
print("=" * 70)
# 使用El Centro波进行时程分析
u = modal_analysis_seismic(M, C, K, phi, omega_n, acc_el, dt)
# 响应统计
u_max = np.max(np.abs(u), axis=1)
print(f"\n最大位移响应：")
print(f" 竖向：{u_max[0]*1000:.2f} mm")
print(f" 横向：{u_max[1]*1000:.2f} mm")
print(f" 扭转：{u_max[2]*1000:.3f} rad")
# 绘制响应时程
fig3, axes3 = plt.subplots(3, 1, figsize=(12, 10))
axes3[0].plot(t_el, u[0, :] * 1000, 'b-', linewidth=0.8)
axes3[0].set_ylabel('位移 (mm)')
axes3[0].set_title('竖向位移响应')
axes3[0].grid(True, alpha=0.3)
axes3[1].plot(t_el, u[1, :] * 1000, 'g-', linewidth=0.8)
axes3[1].set_ylabel('位移 (mm)')
axes3[1].set_title('横向位移响应')
axes3[1].grid(True, alpha=0.3)
axes3[2].plot(t_el, u[2, :], 'r-', linewidth=0.8)
axes3[2].set_xlabel('时间 (s)')
axes3[2].set_ylabel('转角 (rad)')
axes3[2].set_title('扭转响应')
axes3[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('case4_seismic_response.png', dpi=150, bbox_inches='tight')
print("地震响应图形已保存")
plt.close()
# 4. 多工况对比
print("\n" + "=" * 70)
print("4. 多工况地震响应对比")
print("=" * 70)
# 不同PGA水平
pga_levels = [0.1, 0.2, 0.35, 0.5] # g
max_displacements = []
for pga_level in pga_levels:
acc_scaled = acc_el * (pga_level * 9.81 / pga_el)
u_scaled = modal_analysis_seismic(M, C, K, phi, omega_n, acc_scaled, dt)
max_displacements.append(np.max(np.abs(u_scaled[0, :])) * 1000) # 竖向位移，mm
print(f"\n不同地震水平下的最大竖向位移：")
for i, pga in enumerate(pga_levels):
print(f" PGA = {pga*100:.0f}%g：{max_displacements[i]:.2f} mm")
# 绘制对比
fig4, ax4 = plt.subplots(figsize=(10, 6))
ax4.bar(range(len(pga_levels)), max_displacements, color='steelblue', edgecolor='black')
ax4.set_xticks(range(len(pga_levels)))
ax4.set_xticklabels([f'{p*100:.0f}%g' for p in pga_levels])
ax4.set_xlabel('地震水平')
ax4.set_ylabel('最大竖向位移 (mm)')
ax4.set_title('不同地震水平下的最大位移')
ax4.grid(True, alpha=0.3, axis='y')
for i, d in enumerate(max_displacements):
ax4.text(i, d + 1, f'{d:.1f}', ha='center', va='bottom')
plt.tight_layout()
plt.savefig('case4_multi_level_comparison.png', dpi=150, bbox_inches='tight')
print("\n多工况对比图形已保存")
plt.close()
# 5. 模态贡献分析
print("\n" + "=" * 70)
print("5. 模态贡献分析")
print("=" * 70)
# 各模态的最大响应
modal_contributions = []
for i in range(len(omega_n)):
# 单模态响应
M_i = phi[:, i] @ M @ phi[:, i]
gamma_i = -phi[:, i] @ M @ np.ones(len(M)) / M_i
xi_i = (phi[:, i] @ C @ phi[:, i]) / (2 * omega_n[i] * M_i)
m = 1.0
c = 2 * xi_i * omega_n[i] * m
k = omega_n[i]**2 * m
F = gamma_i * acc_el
u_i, _, _ = newmark_beta(np.array([[m]]), np.array([[c]]), np.array([[k]]),
F.reshape(1, -1), dt)
modal_contributions.append(np.max(np.abs(u_i)) * np.abs(phi[0, i]) * 1000)
print(f"\n各模态对竖向位移的贡献：")
for i, contrib in enumerate(modal_contributions):
print(f" 第{i+1}阶模态：{contrib:.2f} mm")
# 绘制模态贡献
fig5, ax5 = plt.subplots(figsize=(10, 6))
mode_numbers = np.arange(1, len(modal_contributions) + 1)
ax5.bar(mode_numbers, modal_contributions, color='coral', edgecolor='black')
ax5.set_xlabel('模态阶数')
ax5.set_ylabel('位移贡献 (mm)')
ax5.set_title('各模态对竖向位移的贡献')
ax5.grid(True, alpha=0.3, axis='y')
for i, c in enumerate(modal_contributions):
ax5.text(i + 1, c + 0.5, f'{c:.1f}', ha='center', va='bottom')
plt.tight_layout()
plt.savefig('case4_modal_contribution.png', dpi=150, bbox_inches='tight')
print("模态贡献图形已保存")
plt.close()
print("\n" + "=" * 70)
print("案例4运行完成！")
print("=" * 70)
if __name__ == "__main__":
main()
# ============================================
# 生成GIF动画
# ============================================
def create_animation():
"""创建仿真结果动画"""
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
fig, ax = plt.subplots(figsize=(10, 6))
ax.set_xlim(0, 10)
ax.set_ylim(-2, 2)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Response', fontsize=12)
ax.set_title('Dynamic Response Animation', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
line, = ax.plot([], [], 'b-', linewidth=2)
def init():
line.set_data([], [])
return line,
def update(frame):
x = np.linspace(0, 10, 100)
y = np.sin(x - frame * 0.2) * np.exp(-frame * 0.01)
line.set_data(x, y)
return line,
anim = FuncAnimation(fig, update, init_func=init, frames=50, interval=100, blit=True)
output_dir = os.path.dirname(os.path.abspath(__file__))
anim.save(f'{output_dir}/simulation_animation.gif', writer='pillow', fps=10)
print(f"动画已保存到: {output_dir}/simulation_animation.gif")
plt.close()
if __name__ == '__main__':
create_animation()
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
13
0- 0
已为社区贡献86条内容
所有评论(0)