第五十一篇：海洋平台结构动力响应

摘要

海洋平台作为海洋油气开发的核心设施，长期处于波浪、风、流、地震等复杂环境荷载作用下，其结构动力响应分析是确保平台安全运营的关键技术。本主题系统介绍海洋平台结构动力响应的分析理论与仿真方法，涵盖导管架平台、半潜式平台、张力腿平台等典型结构形式。首先阐述海洋环境荷载特性，包括波浪理论、风荷载模型、海流作用及地震荷载；然后建立海洋平台的结构动力学模型，考虑流固耦合、桩土相互作用等非线性因素；接着介绍时域和频域动力响应分析方法，重点讲解Morison方程、绕射理论、辐射理论等水动力计算方法；最后通过Python仿真实现固定式平台和浮式平台的动力响应分析，为海洋平台结构设计与安全评估提供理论基础和工程实践指导。

关键词

海洋平台；波浪荷载；流固耦合；动力响应；Morison方程；桩土相互作用；随机振动；疲劳分析

---

1. 海洋平台结构概述

1.1 海洋平台的分类与特点

海洋平台是用于海洋资源开发的大型工程结构，根据结构形式和工作原理可分为以下几类：

固定式平台主要包括导管架平台（Jacket Platform）和重力式平台（Gravity Platform）。导管架平台由钢管构件组成的空间框架结构，通过桩基础固定于海底，适用于水深不超过300米的海域。其优点是结构刚度大、动力响应小、适应恶劣海况能力强；缺点是随着水深增加，结构用钢量急剧增加，经济性下降。重力式平台依靠自身重量和压载保持稳定，不需要桩基础，适用于软土地基，但造价较高。

浮式平台包括半潜式平台（Semi-submersible Platform）、张力腿平台（TLP, Tension Leg Platform）、 Spar平台和浮式生产储油船（FPSO）。半潜式平台由浮箱、立柱和上层甲板组成，通过锚泊系统定位，具有良好的运动性能和较大的甲板面积。张力腿平台通过张力腱与海底连接，在垂直方向具有刚性约束，水平方向可随波运动，适用于深水开发。Spar平台为单柱式结构，重心低于浮心，具有良好的稳定性。FPSO是改装或新建的船舶式平台，具有生产、储存和卸载功能。

顺应式平台如顺应塔平台（Compliant Tower），通过柔性结构适应海洋环境，允许在一定范围内运动，从而减小环境荷载。

1.2 海洋环境荷载特性

海洋平台面临的环境荷载复杂多变，主要包括：

波浪荷载是海洋平台设计的主要控制荷载。波浪由风作用于海面产生，具有随机性和非线性特征。根据波浪理论，可分为线性波浪理论（Airy波）、斯托克斯波浪理论、椭圆余弦波理论等。实际海况中，波浪可视为由多个不同频率、不同方向的简谐波叠加而成的随机过程，常用JONSWAP谱或PM谱描述。

风荷载包括平均风和脉动风两部分。平均风产生静力作用，脉动风引起结构振动。风荷载与风速、风剖面、结构体型系数等因素相关。对于海洋平台，还需考虑风与波浪的联合作用。

海流荷载由潮汐、风生流和密度流引起。海流对平台产生拖曳力和惯性力，同时影响波浪特性。在深水区域，海流剖面随深度变化，表层流速大，底层流速小。

地震荷载对于位于地震活跃区域的海洋平台尤为重要。海底地震通过地基土体传递到平台基础，引起结构振动。海洋平台的地震响应分析需考虑水-结构-地基的相互作用。

冰荷载和海生物附着也是影响平台安全的重要因素，特别是在极地海域和热带海域。

1.3 海洋平台动力响应的特点

与陆地结构相比，海洋平台的动力响应具有以下特点：

流固耦合效应显著：平台结构与水体的相互作用复杂，包括附加质量、辐射阻尼、绕射效应等。对于小构件，可采用Morison方程计算波浪力；对于大尺度构件，需考虑绕射和辐射效应。

非线性因素多：海洋平台的非线性来源包括几何非线性（大位移、大转动）、材料非线性（屈服、疲劳）、边界非线性（桩土相互作用、锚链非线性刚度）和荷载非线性（波浪力的非线性）。

随机振动特征明显：波浪、风等环境荷载本质上是随机过程，平台的响应也是随机过程，需采用随机振动理论进行分析。

多体耦合运动：浮式平台的运动包括纵荡、横荡、垂荡、横摇、纵摇和艏摇六个自由度，各自由度之间存在耦合效应。

疲劳问题突出：海洋平台在长期循环荷载作用下，焊接节点等应力集中部位容易产生疲劳裂纹，疲劳寿命评估是平台设计的重要内容。

2. 海洋环境荷载建模

2.1 波浪理论

2.1.1 线性波浪理论

线性波浪理论（微幅波理论）假设波高与波长相比很小，波浪运动可视为简谐运动。对于深水情况，波面方程为：

$$\eta (x,t)=a\cos (kx-\omega t)$$

其中，$a$为波幅，$k=2\pi /\lambda$为波数，$\omega =2\pi /T$为圆频率，$\lambda$为波长，$T$为周期。

波数与频率满足色散关系：

$${\omega }^2=gk\tanh (kd)$$

其中，$g$为重力加速度，$d$为水深。对于深水（$d/\lambda >0.5$），$\tanh (kd)\approx 1$，色散关系简化为：

$${\omega }^2=gk$$

水质点速度势为：

$$\varphi (x,z,t)=\frac{ag}{\omega }{e}^{kz}\sin (kx-\omega t)$$

水质点水平速度和垂直速度分别为：

$$u=\frac{\partial \varphi }{\partial x}=a\omega e^{kz}\cos (kx-\omega t)$$

$$w=\frac{\partial \varphi }{\partial z}=a\omega e^{kz}\sin (kx-\omega t)$$

水质点加速度为：

$$\dot{u}=a{\omega }^2{e}^{kz}\sin (kx-\omega t)$$

$$\dot{w}=-a{\omega }^2{e}^{kz}\cos (kx-\omega t)$$

2.1.2 非线性波浪理论

当波陡（波高与波长之比）较大时，线性理论不再适用，需采用非线性波浪理论。斯托克斯波浪理论是常用的非线性理论，将速度势展开为波陡的幂级数：

$$\varphi =\sum _{n=1}^N{\varphi }_n$$

其中，${\varphi }_n$为第$n$阶速度势。斯托克斯五阶波理论可较好地描述有限振幅波的波形、水质点运动和波浪破碎前的极限状态。

椭圆余弦波（Cnoidal波）理论适用于浅水情况，其波形为椭圆余弦函数。孤立波（Solitary波）是Cnoidal波在无限波长时的极限情况，波形为双曲正割函数。

2.1.3 随机波浪模型

实际海况中，波浪是随机过程，可用波浪谱描述其统计特性。常用的波浪谱有：

P-M谱（Pierson-Moskowitz谱）：适用于充分发展的风浪

$$S_{\eta \eta }(\omega )=\frac{5.48}{{\omega }^5}{H}_s^2{T}_z^2\exp \left(-1.25{\left(\frac{T_z\omega }{2\pi }\right)}^{-4}\right)$$

其中，$H_s$为有效波高，$T_z$为平均周期。

JONSWAP谱：适用于有限风区的风浪，在P-M谱基础上增加峰值提升因子

$$S_{\eta \eta }(\omega )=\alpha g^2{\omega }^{-5}\exp \left(-1.25{\left(\frac{{\omega }_p}{\omega }\right)}^4\right){\gamma }^{\exp \left(-\frac{(\omega -{\omega }_p{)}^2}{2{\sigma }^2{\omega }_p^2}\right)}$$

其中，${\omega }_p$为谱峰频率，$\gamma$为峰值提升因子（通常取1.6-3.3），$\sigma$为峰形参数。

方向谱：考虑波浪传播方向的分布

$$S(f,\theta )=S(f)G(f,\theta )$$

其中，$G(f,\theta )$为方向分布函数，常用${\cos }^2$或${\cos }^4$形式。

2.2 波浪荷载计算

2.2.1 Morison方程

对于直径与波长比小于0.2的小尺度构件，可采用Morison方程计算波浪力。Morison方程认为波浪力由惯性力和拖曳力两部分组成：

$$F=F_I+F_D=C_M\rho V\dot{u}+\frac{1}{2}{C}_D\rho A|u|u$$

其中，$C_M$为惯性力系数（通常取2.0），$C_D$为拖曳力系数（通常取0.6-1.2），$\rho$为水密度，$V$为单位长度构件排水体积，$A$为单位长度构件投影面积，$u$和$\dot{u}$分别为水质点水平和加速度。

对于圆柱体，单位长度波浪力为：

$$F=C_M\rho \frac{\pi D^2}{4}\dot{u}+\frac{1}{2}{C}_D\rho D|u|u$$

Morison方程中的惯性力系数和拖曳力系数与雷诺数、KC数（Keulegan-Carpenter数）和表面粗糙度有关，需通过试验确定。

2.2.2 绕射理论

对于大尺度构件（直径与波长比大于0.2），波浪遇到结构物时发生绕射，需采用绕射理论计算波浪力。基于线性势流理论，速度势可分解为入射势、绕射势和辐射势：

$$\varphi ={\varphi }_I+{\varphi }_D+{\varphi }_R$$

入射势为已知，绕射势满足物面边界条件（法向速度为零），辐射势与结构运动相关。

通过求解拉普拉斯方程和边界条件，得到速度势后，压力由线性伯努利方程计算：

$$p=-\rho \frac{\partial \varphi }{\partial t}-\rho gz$$

波浪力通过对物面积分得到：

$$F=-{\int }_{S}p\mathbf{n}dS$$

对于简单几何形状（圆柱、球体），有解析解；对于复杂结构，需采用边界元法或有限元法数值求解。

2.2.3 辐射理论

当结构物运动时，会向外辐射波浪，产生附加质量和辐射阻尼。附加质量表示结构带动周围水体运动所需的等效质量，辐射阻尼表示波浪辐射引起的能量耗散。

对于浮式平台，附加质量和阻尼系数是频率的函数，可通过势流理论计算或模型试验确定。在时域分析中，需采用卷积积分考虑频率相关的附加质量和阻尼：

$$F_{rad}(t)=-A_{\infty }\ddot{x}(t)-{\int }_0^{t}K(t-\tau )\dot{x}(\tau )d\tau$$

其中，$A_{\infty }$为无限频率时的附加质量，$K(t)$为脉冲响应函数（延迟函数）。

2.3 风荷载模型

海洋平台的风荷载包括平均风荷载和脉动风荷载。平均风荷载为：

$$F_w=\frac{1}{2}{\rho }_a{C}_{s}A{U}^2$$

其中，${\rho }_a$为空气密度，$C_s$为体型系数，$A$为受风面积，$U$为设计风速。

风速沿高度的分布可用对数律或幂律描述：

对数律：

$$U(z)=\frac{u_{\ast }}{\kappa }\ln \left(\frac{z}{z_0}\right)$$

其中，$u_{\ast }$为摩擦速度，$\kappa$为卡门常数（约0.4），$z_0$为粗糙长度。

幂律：

$$U(z)=U_{10}{\left(\frac{z}{10}\right)}^{\alpha }$$

其中，$U_{10}$为10米高度处风速，$\alpha$为风剖面指数（海上通常取0.12-0.16）。

脉动风引起结构振动，其功率谱常用Kaimal谱或Davenport谱描述。对于海洋平台，还需考虑风与波浪的联合分布，通常采用环境等值线法或载荷组合系数法。

2.4 海流与地震荷载

海流荷载与波浪荷载类似，可采用Morison方程计算。海流速度剖面随深度变化，常用指数分布或线性分布描述。海流与波浪联合作用时，水质点速度为波浪速度与海流速度的矢量和。

地震荷载通过地震波传播到平台基础。海底地震动特性与陆地有所不同，需考虑海底土层的放大效应和水体的附加质量。海洋平台的地震响应分析可采用反应谱法或时程分析法，需考虑水-结构-地基的相互作用。

3. 海洋平台结构动力学模型

3.1 导管架平台动力学模型

导管架平台可简化为空间杆系结构，采用有限元法建立动力学模型。结构质量矩阵包括结构质量、附加质量和设备质量：

$$\mathbf{M}={\mathbf{M}}_s+{\mathbf{M}}_a+{\mathbf{M}}_e$$

其中，附加质量包括杆件周围水体的附加质量和桩基础的附加质量。

刚度矩阵考虑杆件轴向刚度、弯曲刚度和剪切刚度。对于大位移情况，需考虑几何刚度（初应力刚度）矩阵。

阻尼矩阵通常采用Rayleigh阻尼：

$$\mathbf{C}=\alpha \mathbf{M}+\beta \mathbf{K}$$

其中，$\alpha$和$\beta$由两阶模态阻尼比确定。

3.2 桩土相互作用模型

桩基础是导管架平台的重要组成部分，桩土相互作用对平台动力特性有显著影响。常用的桩土相互作用模型有：

p-y曲线法：描述桩侧土抗力与桩水平位移的关系

$$p=k_{h}y$$

其中，$p$为单位长度土抗力，$y$为桩位移，$k_h$为地基水平反力系数，随深度和土性变化。API规范给出了砂土和粘土中p-y曲线的建议公式。

t-z曲线法和q-z曲线法：分别描述桩侧摩阻力和桩端阻力与桩竖向位移的关系。

在动力分析中，桩土相互作用可用等效弹簧-阻尼器模拟，弹簧刚度由p-y曲线割线模量确定，阻尼包括辐射阻尼和材料阻尼。

3.3 浮式平台动力学模型

浮式平台的运动方程为：

$$(\mathbf{M}+\mathbf{A})\ddot{\mathbf{x}}+\mathbf{B}\dot{\mathbf{x}}+\mathbf{Kx}={\mathbf{F}}_{wave}+{\mathbf{F}}_{wind}+{\mathbf{F}}_{current}+{\mathbf{F}}_{moor}$$

其中，$\mathbf{M}$为结构质量矩阵，$\mathbf{A}$为附加质量矩阵，$\mathbf{B}$为阻尼矩阵（包括辐射阻尼和粘性阻尼），$\mathbf{K}$为恢复力矩阵（包括静水恢复力和锚泊恢复力），$\mathbf{x}$为六自由度位移向量。

附加质量：浮体在流体中运动时带动周围水体运动，产生附加质量。附加质量是频率的函数，对于规则形状有理论解，复杂形状需通过势流计算或试验确定。

辐射阻尼：浮体运动辐射波浪，产生能量耗散。辐射阻尼也是频率的函数，在共振频率附近较大。

恢复力：包括静水恢复力和锚泊恢复力。静水恢复力由浮心和重心的相对位置产生，锚泊恢复力由锚链张力产生。对于TLP，张力腱提供主要的垂向恢复力。

锚泊系统：锚泊线（链或缆）具有非线性刚度特性，其恢复力与位移的关系为：

$$T=T_0+k_{eq}\Delta x$$

其中，$T_0$为预张力，$k_{eq}$为等效刚度，与锚泊线长度、水深、弹性模量等因素有关。

3.4 流固耦合分析

流固耦合是海洋平台动力分析的核心问题。根据结构尺度与波长的比值，可采用不同的分析方法：

小尺度构件（$D/\lambda <0.2$）：采用Morison方程，考虑惯性力和拖曳力，忽略绕射效应。

大尺度构件（$D/\lambda >0.2$）：采用绕射理论，求解速度势，计算波浪力。

浮式平台：采用势流理论计算附加质量、阻尼和波浪激励力，考虑辐射和绕射效应。

在时域分析中，对于频率相关的附加质量和阻尼，可采用Cummins方程：

$$(\mathbf{M}+{\mathbf{A}}_{\infty })\ddot{\mathbf{x}}+{\int }_0^t\mathbf{K}(t-\tau )\dot{\mathbf{x}}(\tau )d\tau +\mathbf{Kx}={\mathbf{F}}_{exc}(t)$$

其中，$\mathbf{K}(t)$为延迟函数，可通过频率域的附加质量和阻尼经傅里叶变换得到。

4. 动力响应分析方法

4.1 频域分析法

频域分析法适用于线性系统和规则波情况，计算效率高，可方便地得到统计特性。

对于线性系统，传递函数为：

$$H(\omega )=\frac{1}{-{\omega }^2(\mathbf{M}+\mathbf{A})+i\omega \mathbf{B}+\mathbf{K}}$$

响应谱为：

$$S_{xx}(\omega )=|H(\omega ){|}^2{S}_{FF}(\omega )$$

其中，$S_{FF}(\omega )$为波浪力谱。对于规则波，波浪力谱为离散谱；对于不规则波，波浪力谱由波浪谱和传递函数确定。

响应的统计特性可由谱矩计算：

$$m_n={\int }_0^{\infty }{\omega }^n{S}_{xx}(\omega )d\omega$$

均方根响应：

$${\sigma }_x=\sqrt{m_0}$$

平均周期：

$$T_z=2\pi \sqrt{\frac{m_0}{m_2}}$$

极值分布可用Rayleigh分布或Weibull分布描述。

4.2 时域分析法

时域分析法可考虑非线性因素，适用于极端海况和疲劳分析。运动方程为：

$$\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{Kx}=\mathbf{F}(t,\mathbf{x},\dot{\mathbf{x}})$$

其中，荷载$\mathbf{F}$可能与位移和速度有关（如拖曳力）。

常用的时间积分方法有Newmark-β法、Wilson-θ法和中心差分法。对于海洋平台，需采用较小的时间步长（通常T/100-T/200，T为波浪周期）以保证精度和稳定性。

时域分析可得到完整的响应时程，便于进行疲劳分析和极端响应统计。

4.3 随机振动分析

海洋环境荷载本质上是随机过程，平台的响应也是随机过程。随机振动分析包括：

线性随机振动：对于线性系统，输入输出的功率谱关系为：

$$S_{yy}(\omega )=|H(\omega ){|}^2{S}_{xx}(\omega )$$

响应的均值、方差、自相关函数等统计特性可由谱分析得到。

非线性随机振动：对于非线性系统，可采用等效线性化法、摄动法或Monte Carlo模拟法。等效线性化法将非线性系统用等效线性系统近似，使误差最小。

首次穿越问题：关心结构响应首次超过安全阈值的概率，与结构的可靠性和安全性评估相关。

4.4 疲劳分析

海洋平台在长期循环荷载作用下，焊接节点等部位容易产生疲劳裂纹。疲劳分析包括：

疲劳荷载：由波浪引起的应力循环，可用应力谱或应力时程描述。

S-N曲线：描述应力幅与疲劳寿命的关系，常用对数坐标下的直线表示：

$$\log N=\log a-m\log \Delta \sigma$$

其中，$N$为疲劳寿命（循环次数），$\Delta \sigma$为应力幅，$a$和$m$为材料常数。

Miner线性累积损伤理论：

$$D=\sum _i\frac{n_i}{N_i}$$

当累积损伤$D\ge 1$时，认为发生疲劳破坏。

疲劳寿命计算：根据应力谱和S-N曲线，计算年疲劳损伤，进而预测平台寿命。

5. Python仿真实现

5.1 案例1：海洋平台环境荷载建模

本案例实现海洋环境荷载的建模与仿真，包括波浪、风、流的模拟。

---

补充内容：深入理解与工程应用

工程实践中的关键考虑因素

在实际工程应用中，理论分析结果需要结合工程经验进行综合考虑。以下是几个关键考虑因素：

1. 材料非线性效应

实际工程材料往往表现出复杂的非线性行为，包括：

• 几何非线性：大变形条件下，结构的刚度矩阵会发生变化，需要考虑几何刚度效应。
• 材料非线性：钢材的屈服、混凝土的开裂都会导致材料本构关系的非线性变化。
• 接触非线性：结构连接部位的接触状态变化会影响整体动力特性。

在动力分析中，这些非线性效应可能导致：

• 固有频率随振幅变化
• 振动响应出现谐波成分
• 能量耗散机制复杂化

2. 环境因素的影响

环境因素对结构动力特性有显著影响：

温度效应：

• 材料弹性模量随温度变化
• 热应力改变结构的初始应力状态
• 温度梯度引起的热变形

湿度与腐蚀：

• 混凝土含水率影响刚度
• 钢材腐蚀降低截面特性
• 长期性能退化

风荷载与气动效应：

• 涡激振动
• 颤振失稳
• 气动阻尼

3. 建模不确定性处理

工程分析中存在多种不确定性来源：

参数不确定性：

• 材料参数的离散性
• 几何尺寸的制造误差
• 边界条件的简化假设

模型不确定性：

• 理论模型的简化
• 数值方法的近似
• 边界条件的理想化

处理方法：

• 灵敏度分析识别关键参数
• 蒙特卡洛模拟量化不确定性
• 模型修正技术提高精度

常见问题与解决方案

问题1：计算结果与试验结果偏差较大

可能原因：

• 边界条件假设不合理
• 材料参数取值不准确
• 模型简化过度
• 数值误差累积

解决方案：

• 进行边界条件敏感性分析
• 通过试验数据修正模型参数
• 逐步细化模型，对比不同复杂度模型的结果
• 检查数值方法的稳定性和收敛性

问题2：数值计算不收敛

可能原因：

• 时间步长过大
• 刚度矩阵病态
• 非线性程度过高
• 初始条件设置不当

解决方案：

• 减小时间步长或采用自适应步长
• 检查并改善模型网格质量
• 采用分步加载或弧长法
• 调整初始条件，避免突变

问题3：高阶模态计算困难

可能原因：

• 网格密度不足
• 数值算法精度限制
• 刚体模态干扰
• 局部振动模态混杂

解决方案：

• 在关键区域加密网格
• 采用高精度特征值算法
• 使用刚体约束消除零频模态
• 应用模态选择技术

进阶案例分析

案例：复杂结构的模态分析策略

考虑一个具有多种构件类型的复杂结构，需要采用分层建模策略：

第一层：整体模型

• 建立简化整体模型
• 确定整体振动特性
• 识别关键振动方向

第二层：子结构细化

• 对关键部位建立详细模型
• 分析局部振动特性
• 评估局部-整体相互作用

第三层：连接部位分析

• 详细建模连接节点
• 分析连接刚度影响
• 验证传力路径

分析流程：

1. 建立整体简化模型，进行初步分析
2. 根据初步结果识别关键区域
3. 对关键区域进行细化建模
4. 通过子结构技术整合分析结果
5. 验证并迭代优化模型

案例：动力试验与仿真对比

试验设计要点：

• 激励方式选择（锤击、激振器、环境激励）
• 测点布置优化
• 数据采集参数设置
• 信号处理与模态参数识别

对比分析方法：

• MAC矩阵评估振型相关性
• 频率误差分析
• 阻尼比对比（考虑阻尼机理差异）
• 响应时程对比

模型修正策略：

• 基于灵敏度分析的参数选择
• 贝叶斯方法量化不确定性
• 迭代修正直至满足精度要求

行业发展趋势

1. 智能化仿真技术

机器学习应用：

• 代理模型加速优化
• 数据驱动的模型修正
• 异常检测与故障诊断
• 响应预测与实时控制

数字孪生技术：

• 实时模型更新
• 虚实交互验证
• 全生命周期管理
• 预测性维护

2. 高性能计算应用

并行计算：

• 分布式内存并行
• GPU加速计算
• 云计算平台应用
• 大规模问题求解

新型算法：

• 无网格方法
• 等几何分析
• 降阶建模技术
• 多尺度分析方法

3. 多物理场耦合分析

流固耦合：

• 风-结构相互作用
• 水-结构相互作用
• 气动弹性问题

热结构耦合：

• 热振动分析
• 热屈曲问题
• 火灾作用下结构响应

机电耦合：

• 智能结构控制
• 压电传感与驱动
• 能量收集系统

学习建议与资源

理论基础深化

推荐教材：

1. 《结构动力学》 - Clough & Penzien
2. 《振动理论及应用》 - Thomson
3. 《有限元方法》 - Bathe
4. 《随机振动》 - 方同

关键数学基础：

• 线性代数与矩阵理论
• 常微分方程与偏微分方程
• 概率论与随机过程
• 数值分析方法

软件技能培养

商业软件：

• ANSYS Mechanical
• ABAQUS
• MSC Nastran
• SAP2000
• ETABS

编程工具：

• Python（科学计算生态）
• MATLAB
• Julia（高性能计算）

工程实践积累

规范标准学习：

• 建筑抗震设计规范
• 桥梁抗震设计规范
• 风荷载规范
• 振动控制相关标准

工程案例研究：

• 典型结构动力特性
• 振动控制工程实例
• 灾害案例分析
• 创新技术应用

---

完整Python代码实现

以下是本主题的完整Python仿真代码：

# -*- coding: utf-8 -*-
"""
案例1：海洋平台环境荷载建模
本案例实现海洋环境荷载的建模与仿真，包括波浪、风、流的模拟
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.integrate import trapezoid
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\结构动力学仿真\主题051'
os.makedirs(output_dir, exist_ok=True)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("案例1：海洋平台环境荷载建模")
print("=" * 60)

# ==================== 参数设置 ====================
g = 9.81  # 重力加速度 (m/s^2)
rho_water = 1025  # 海水密度 (kg/m^3)
rho_air = 1.225  # 空气密度 (kg/m^3)

# 海洋环境参数
Hs = 8.0  # 有效波高 (m)
Tz = 10.0  # 平均周期 (s)
d_water = 100.0  # 水深 (m)
U10 = 25.0  # 10米高度处风速 (m/s)
U_current_surface = 1.5  # 表层海流速度 (m/s)

print("\n【海洋环境参数】")
print(f"有效波高 Hs = {Hs} m")
print(f"平均周期 Tz = {Tz} s")
print(f"水深 d = {d_water} m")
print(f"10米高度风速 U10 = {U10} m/s")
print(f"表层海流速度 = {U_current_surface} m/s")

# ==================== 1. 波浪理论 ====================
print("\n" + "=" * 60)
print("1. 波浪理论仿真")
print("=" * 60)

# 1.1 线性波浪理论
def linear_wave(x, z, t, a, T, d, g=9.81):
    """
    线性波浪理论 (Airy波)
    计算波面、水质点速度和加速度
    """
    omega = 2 * np.pi / T  # 圆频率
    
    # 色散关系求解波数 k
    # omega^2 = g*k*tanh(k*d)
    # 使用迭代法求解
    k = omega**2 / g  # 深水初值
    for _ in range(20):
        k_new = omega**2 / (g * np.tanh(k * d))
        if abs(k_new - k) < 1e-6:
            break
        k = k_new
    
    # 波面
    eta = a * np.cos(k * x - omega * t)
    
    # 水质点速度
    u = a * omega * np.cosh(k * (z + d)) / np.sinh(k * d) * np.cos(k * x - omega * t)
    w = a * omega * np.sinh(k * (z + d)) / np.sinh(k * d) * np.sin(k * x - omega * t)
    
    # 水质点加速度
    u_dot = a * omega**2 * np.cosh(k * (z + d)) / np.sinh(k * d) * np.sin(k * x - omega * t)
    w_dot = -a * omega**2 * np.sinh(k * (z + d)) / np.sinh(k * d) * np.cos(k * x - omega * t)
    
    return eta, u, w, u_dot, w_dot, k, omega

# 计算波浪参数
a_wave = Hs / 2  # 波幅
omega_wave = 2 * np.pi / Tz

# 求解波数
k_wave = omega_wave**2 / g
for _ in range(20):
    k_new = omega_wave**2 / (g * np.tanh(k_wave * d_water))
    if abs(k_new - k_wave) < 1e-6:
        break
    k_wave = k_new

lambda_wave = 2 * np.pi / k_wave  # 波长
print(f"\n波浪参数:")
print(f"  波幅 a = {a_wave:.2f} m")
print(f"  圆频率 omega = {omega_wave:.4f} rad/s")
print(f"  波数 k = {k_wave:.4f} rad/m")
print(f"  波长 lambda = {lambda_wave:.2f} m")
print(f"  水深波长比 d/lambda = {d_water/lambda_wave:.3f}")

# 绘制波浪传播过程
x_wave = np.linspace(0, 2*lambda_wave, 200)
t_wave = np.linspace(0, 2*Tz, 100)
X_wave, T_wave = np.meshgrid(x_wave, t_wave)

Eta_wave = a_wave * np.cos(k_wave * X_wave - omega_wave * T_wave)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：波浪传播时空图
ax1 = axes[0, 0]
im1 = ax1.contourf(X_wave, T_wave, Eta_wave, levels=20, cmap='RdBu_r')
ax1.set_xlabel('Position x (m)', fontsize=11)
ax1.set_ylabel('Time t (s)', fontsize=11)
ax1.set_title('Wave Propagation (Space-Time)', fontsize=12, fontweight='bold')
plt.colorbar(im1, ax=ax1, label='Elevation (m)')

# 图2：不同时刻的波面形状
ax2 = axes[0, 1]
t_snapshots = [0, Tz/4, Tz/2, 3*Tz/4]
colors = ['blue', 'green', 'red', 'purple']
for i, t_snap in enumerate(t_snapshots):
    eta_snap = a_wave * np.cos(k_wave * x_wave - omega_wave * t_snap)
    ax2.plot(x_wave, eta_snap, color=colors[i], linewidth=2, 
             label=f't = {t_snap:.1f}s ({i*90}°)')
ax2.axhline(y=0, color='black', linestyle='--', alpha=0.5)
ax2.set_xlabel('Position x (m)', fontsize=11)
ax2.set_ylabel('Wave Elevation (m)', fontsize=11)
ax2.set_title('Wave Profiles at Different Times', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-a_wave*1.2, a_wave*1.2)

# 图3：水质点运动轨迹
ax3 = axes[1, 0]
x_particle = 0  # 观察点位置
z_particles = np.linspace(-d_water, 0, 10)  # 不同深度
t_particle = np.linspace(0, Tz, 100)

for z_p in z_particles:
    eta_p, u_p, w_p, _, _, _, _ = linear_wave(x_particle, z_p, t_particle, 
                                               a_wave, Tz, d_water)
    # 计算位移 (近似积分)
    dx = integrate.cumulative_trapezoid(u_p, t_particle, initial=0)
    dz = integrate.cumulative_trapezoid(w_p, t_particle, initial=0)
    ax3.plot(x_particle + dx, z_p + dz, alpha=0.7, linewidth=1.5)
    ax3.plot(x_particle, z_p, 'ko', markersize=3)

ax3.set_xlabel('Horizontal Position (m)', fontsize=11)
ax3.set_ylabel('Vertical Position z (m)', fontsize=11)
ax3.set_title('Water Particle Trajectories', fontsize=12, fontweight='bold')
ax3.axhline(y=0, color='blue', linestyle='-', alpha=0.5, label='Still Water Level')
ax3.set_xlim(-a_wave*2, a_wave*2)
ax3.set_ylim(-d_water, a_wave)
ax3.grid(True, alpha=0.3)
ax3.legend(fontsize=9)

# 图4：水质点速度和加速度分布
ax4 = axes[1, 1]
z_profile = np.linspace(-d_water, 0, 50)
t_ref = 0  # 参考时刻

u_profile = []
w_profile = []
u_dot_profile = []

for z_p in z_profile:
    _, u_p, w_p, u_dot_p, w_dot_p, _, _ = linear_wave(0, z_p, t_ref, 
                                                       a_wave, Tz, d_water)
    u_profile.append(u_p)
    w_profile.append(w_p)
    u_dot_profile.append(u_dot_p)

ax4.plot(u_profile, z_profile, 'b-', linewidth=2, label='Horizontal velocity u')
ax4.plot(w_profile, z_profile, 'r-', linewidth=2, label='Vertical velocity w')
ax4.plot(u_dot_profile, z_profile, 'g--', linewidth=2, label='Horizontal accel. du/dt')
ax4.axvline(x=0, color='black', linestyle='-', alpha=0.3)
ax4.axhline(y=0, color='blue', linestyle='-', alpha=0.3)
ax4.set_xlabel('Velocity / Acceleration', fontsize=11)
ax4.set_ylabel('Depth z (m)', fontsize=11)
ax4.set_title('Velocity and Acceleration Profiles', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_wave_theory.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 波浪理论仿真图已保存")

# ==================== 2. 随机波浪谱 ====================
print("\n" + "=" * 60)
print("2. 随机波浪谱仿真")
print("=" * 60)

# 2.1 JONSWAP谱
def jonswap_spectrum(omega, Hs, Tp, gamma=3.3):
    """
    JONSWAP波浪谱
    omega: 圆频率 (rad/s)
    Hs: 有效波高 (m)
    Tp: 谱峰周期 (s)
    gamma: 峰值提升因子
    """
    omega_p = 2 * np.pi / Tp  # 谱峰频率
    
    # 参数
    alpha = 0.0081  # Phillips常数
    sigma = np.where(omega <= omega_p, 0.07, 0.09)
    
    # JONSWAP谱
    S = (alpha * g**2 / omega**5) * np.exp(-1.25 * (omega_p / omega)**4)
    S *= gamma ** np.exp(-0.5 * ((omega - omega_p) / (sigma * omega_p))**2)
    
    # 归一化使 Hs = 4*sqrt(m0)
    m0 = trapezoid(S, omega)
    S = S * (Hs**2 / 16) / m0
    
    return S

# 2.2 Pierson-Moskowitz谱
def pm_spectrum(omega, Hs, Tz):
    """
    Pierson-Moskowitz波浪谱 (充分发展风浪)
    """
    omega_p = 2 * np.pi / Tz
    S = (5.48 / omega**5) * Hs**2 * Tz**2 * np.exp(-1.25 * (omega_p / omega)**4)
    return S

# 频率范围
omega_range = np.linspace(0.1, 3.0, 500)

# 计算不同海况的谱
# 海况1: 轻度海况
Hs1, Tp1 = 3.0, 8.0
S1 = jonswap_spectrum(omega_range, Hs1, Tp1, gamma=2.0)

# 海况2: 中度海况
Hs2, Tp2 = 6.0, 10.0
S2 = jonswap_spectrum(omega_range, Hs2, Tp2, gamma=3.0)

# 海况3: 恶劣海况 (设计海况)
Hs3, Tp3 = 12.0, 14.0
S3 = jonswap_spectrum(omega_range, Hs3, Tp3, gamma=3.3)

# 计算谱矩和特征参数
def spectral_moments(omega, S, n):
    """计算谱矩 m_n = ∫ omega^n * S(omega) domega"""
    return trapezoid(omega**n * S, omega)

print("\n不同海况的波浪谱参数:")
for i, (Hs, Tp, S, name) in enumerate([(Hs1, Tp1, S1, "轻度"), 
                                        (Hs2, Tp2, S2, "中度"), 
                                        (Hs3, Tp3, S3, "恶劣")], 1):
    m0 = spectral_moments(omega_range, S, 0)
    m2 = spectral_moments(omega_range, S, 2)
    m4 = spectral_moments(omega_range, S, 4)
    
    Hs_calc = 4 * np.sqrt(m0)
    Tz_calc = 2 * np.pi * np.sqrt(m0 / m2)
    
    print(f"\n海况{i} ({name}):")
    print(f"  输入有效波高 Hs = {Hs} m, 谱峰周期 Tp = {Tp} s")
    print(f"  计算有效波高 Hs = {Hs_calc:.2f} m")
    print(f"  平均周期 Tz = {Tz_calc:.2f} s")
    print(f"  谱峰频率 omega_p = {2*np.pi/Tp:.4f} rad/s")

# 绘制波浪谱
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：不同海况的JONSWAP谱
ax1 = axes[0, 0]
ax1.plot(omega_range, S1, 'b-', linewidth=2, label=f'Mild: Hs={Hs1}m, Tp={Tp1}s')
ax1.plot(omega_range, S2, 'g-', linewidth=2, label=f'Moderate: Hs={Hs2}m, Tp={Tp2}s')
ax1.plot(omega_range, S3, 'r-', linewidth=2, label=f'Severe: Hs={Hs3}m, Tp={Tp3}s')
ax1.fill_between(omega_range, S1, alpha=0.3, color='blue')
ax1.fill_between(omega_range, S2, alpha=0.3, color='green')
ax1.fill_between(omega_range, S3, alpha=0.3, color='red')
ax1.set_xlabel('Frequency omega (rad/s)', fontsize=11)
ax1.set_ylabel('Spectral Density S(omega) (m^2.s)', fontsize=11)
ax1.set_title('JONSWAP Wave Spectra for Different Sea States', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 3)

# 图2：JONSWAP vs PM谱比较
ax2 = axes[0, 1]
S_pm = pm_spectrum(omega_range, Hs2, Tz)
ax2.plot(omega_range, S2, 'b-', linewidth=2, label='JONSWAP (gamma=3.0)')
ax2.plot(omega_range, S_pm, 'r--', linewidth=2, label='Pierson-Moskowitz')
ax2.set_xlabel('Frequency omega (rad/s)', fontsize=11)
ax2.set_ylabel('Spectral Density S(omega) (m^2.s)', fontsize=11)
ax2.set_title('JONSWAP vs PM Spectrum Comparison', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 图3：基于谱合成随机波浪时程
ax3 = axes[1, 0]

# 使用IFFT方法合成随机波浪
def synthesize_wave_time_series(omega, S, duration, dt, seed=None):
    """
    基于波浪谱合成随机波浪时程
    """
    if seed is not None:
        np.random.seed(seed)
    
    N = int(duration / dt)
    t = np.linspace(0, duration, N)
    
    # 频率分辨率
    df = 1 / duration
    f = np.arange(N) * df
    omega_fft = 2 * np.pi * f
    
    # 插值得到各频率的谱值
    S_interp = np.interp(omega_fft, omega, S, left=0, right=0)
    
    # 振幅谱
    A = np.sqrt(2 * S_interp * 2 * np.pi * df)
    
    # 随机相位
    phi = np.random.uniform(0, 2*np.pi, N)
    
    # 构造频域信号
    X = A * np.exp(1j * phi)
    
    # IFFT得到时域信号
    eta = np.fft.ifft(X).real * N
    
    return t, eta

t_wave_syn, eta_wave_syn = synthesize_wave_time_series(omega_range, S2, 
                                                        duration=300, dt=0.5, seed=42)

ax3.plot(t_wave_syn, eta_wave_syn, 'b-', linewidth=0.8, alpha=0.8)
ax3.axhline(y=Hs2/2, color='r', linestyle='--', label=f'Hs/2 = {Hs2/2}m')
ax3.axhline(y=-Hs2/2, color='r', linestyle='--')
ax3.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax3.set_xlabel('Time (s)', fontsize=11)
ax3.set_ylabel('Wave Elevation (m)', fontsize=11)
ax3.set_title('Synthesized Random Wave Time Series', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 300)

# 图4：波浪统计特性
ax4 = axes[1, 1]

# 计算波高统计
from scipy.signal import find_peaks
peaks, _ = find_peaks(eta_wave_syn, height=0)
wave_heights = []
for peak in peaks:
    # 找到相邻的波谷
    if peak > 10 and peak < len(eta_wave_syn) - 10:
        valley_before = np.min(eta_wave_syn[peak-10:peak])
        valley_after = np.min(eta_wave_syn[peak:peak+10])
        wave_height = eta_wave_syn[peak] - min(valley_before, valley_after)
        wave_heights.append(wave_height)

wave_heights = np.array(wave_heights)
H_max = np.max(wave_heights)
H_avg = np.mean(wave_heights)
H_s_calc = 4 * np.std(eta_wave_syn)

# 绘制波高分布
ax4.hist(wave_heights, bins=20, density=True, alpha=0.7, color='steelblue', edgecolor='black')
ax4.axvline(x=H_s_calc, color='r', linestyle='--', linewidth=2, label=f'Hs calc = {H_s_calc:.2f}m')
ax4.axvline(x=H_avg, color='g', linestyle='--', linewidth=2, label=f'H avg = {H_avg:.2f}m')
ax4.set_xlabel('Wave Height (m)', fontsize=11)
ax4.set_ylabel('Probability Density', fontsize=11)
ax4.set_title('Wave Height Distribution', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_wave_spectrum.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 随机波浪谱仿真图已保存")

# ==================== 3. 风荷载模型 ====================
print("\n" + "=" * 60)
print("3. 风荷载模型仿真")
print("=" * 60)

# 3.1 风速剖面
def wind_profile_log(z, u_star, z0=0.0002, kappa=0.4):
    """
    对数律风速剖面
    z: 高度 (m)
    u_star: 摩擦速度 (m/s)
    z0: 粗糙长度 (m)
    kappa: 卡门常数
    """
    return (u_star / kappa) * np.log(z / z0)

def wind_profile_power(z, U10, alpha=0.14):
    """
    幂律风速剖面
    z: 高度 (m)
    U10: 10米高度处风速 (m/s)
    alpha: 风剖面指数
    """
    return U10 * (z / 10) ** alpha

# 计算摩擦速度 (假设CD = 0.0012)
CD = 0.0012
u_star = np.sqrt(CD) * U10

z_wind = np.linspace(5, 150, 100)  # 高度范围
U_log = wind_profile_log(z_wind, u_star)
U_power = wind_profile_power(z_wind, U10, alpha=0.14)

print(f"\n风荷载参数:")
print(f"  10米高度设计风速 U10 = {U10} m/s")
print(f"  摩擦速度 u* = {u_star:.4f} m/s")
print(f"  50米高度风速 (对数律) = {wind_profile_log(50, u_star):.2f} m/s")
print(f"  50米高度风速 (幂律) = {wind_profile_power(50, U10):.2f} m/s")

# 3.2 Kaimal风谱
def kaimal_spectrum(f, U10, z=10):
    """
    Kaimal纵向风谱
    f: 频率 (Hz)
    U10: 10米高度处风速 (m/s)
    z: 高度 (m)
    """
    u_star = 0.4 * U10  # 简化计算
    f_star = f * z / U10
    
    # Kaimal谱
    S = (105 * u_star**2 * f_star) / (f * (1 + 33 * f_star)**(5/3))
    return S

# 3.3 风压计算
def wind_pressure(U, rho_air=1.225):
    """
    计算风压
    """
    return 0.5 * rho_air * U**2

# 频率范围
f_wind = np.linspace(0.001, 2.0, 500)
S_wind = kaimal_spectrum(f_wind, U10, z=30)

# 绘制风荷载
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：风速剖面
ax1 = axes[0, 0]
ax1.plot(U_log, z_wind, 'b-', linewidth=2, label='Logarithmic Law')
ax1.plot(U_power, z_wind, 'r--', linewidth=2, label='Power Law (alpha=0.14)')
ax1.axvline(x=U10, color='g', linestyle=':', alpha=0.7, label=f'U10 = {U10} m/s')
ax1.set_xlabel('Wind Speed (m/s)', fontsize=11)
ax1.set_ylabel('Height z (m)', fontsize=11)
ax1.set_title('Wind Speed Profile', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 图2：Kaimal风谱
ax2 = axes[0, 1]
ax2.loglog(f_wind, S_wind, 'b-', linewidth=2)
ax2.set_xlabel('Frequency f (Hz)', fontsize=11)
ax2.set_ylabel('Spectral Density S(f) (m^2/s)', fontsize=11)
ax2.set_title('Kaimal Wind Spectrum', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')

# 图3：合成脉动风速时程
def synthesize_wind_time_series(f, S, duration, dt, U_mean, seed=None):
    """
    合成脉动风速时程
    """
    if seed is not None:
        np.random.seed(seed)
    
    N = int(duration / dt)
    
    # 频率分辨率
    df = f[1] - f[0]
    
    # 振幅谱
    A = np.sqrt(2 * S * df)
    
    # 随机相位
    phi = np.random.uniform(0, 2*np.pi, len(f))
    
    # 时间序列
    t = np.linspace(0, duration, N)
    u_fluct = np.zeros(N)
    
    for i in range(len(f)):
        u_fluct += A[i] * np.cos(2 * np.pi * f[i] * t + phi[i])
    
    U_total = U_mean + u_fluct
    
    return t, U_total, u_fluct

t_wind_syn, U_wind_syn, u_fluct = synthesize_wind_time_series(
    f_wind, S_wind, duration=300, dt=0.1, U_mean=U10, seed=123)

ax3 = axes[1, 0]
ax3.plot(t_wind_syn, U_wind_syn, 'b-', linewidth=0.8, alpha=0.8, label='Total Wind Speed')
ax3.axhline(y=U10, color='r', linestyle='--', linewidth=2, label=f'Mean U10 = {U10} m/s')
ax3.fill_between(t_wind_syn, U10, U_wind_syn, alpha=0.3, color='gray', label='Fluctuation')
ax3.set_xlabel('Time (s)', fontsize=11)
ax3.set_ylabel('Wind Speed (m/s)', fontsize=11)
ax3.set_title('Synthesized Wind Speed Time Series', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 300)

# 图4：风压分布
ax4 = axes[1, 1]
U_range = np.linspace(0, 50, 100)
P_wind = wind_pressure(U_range)

ax4.plot(U_range, P_wind, 'b-', linewidth=2)
ax4.fill_between(U_range, P_wind, alpha=0.3)
ax4.axvline(x=U10, color='r', linestyle='--', linewidth=2, label=f'U10 = {U10} m/s')
ax4.axhline(y=wind_pressure(U10), color='r', linestyle='--', alpha=0.5)
ax4.scatter([U10], [wind_pressure(U10)], color='red', s=100, zorder=5)
ax4.set_xlabel('Wind Speed (m/s)', fontsize=11)
ax4.set_ylabel('Wind Pressure (Pa)', fontsize=11)
ax4.set_title('Wind Pressure vs Wind Speed', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)
ax4.text(U10+2, wind_pressure(U10), f'q = {wind_pressure(U10):.1f} Pa', 
         fontsize=10, color='red')

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_wind_load.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 风荷载模型仿真图已保存")

# ==================== 4. 海流荷载 ====================
print("\n" + "=" * 60)
print("4. 海流荷载仿真")
print("=" * 60)

# 4.1 海流速度剖面
def current_profile_exp(z, U0, d, decay_depth=50):
    """
    指数分布海流剖面
    z: 深度 (m, 负值)
    U0: 表层流速 (m/s)
    d: 水深 (m)
    decay_depth: 流速衰减特征深度 (m)
    """
    return U0 * np.exp((z + d) / decay_depth)

def current_profile_linear(z, U0, d):
    """
    线性分布海流剖面
    """
    return U0 * (z + d) / d

z_current = np.linspace(-d_water, 0, 100)
U_current_exp = current_profile_exp(z_current, U_current_surface, d_water)
U_current_linear = current_profile_linear(z_current, U_current_surface, d_water)

print(f"\n海流参数:")
print(f"  表层流速 = {U_current_surface} m/s")
print(f"  50米深度流速 (指数) = {current_profile_exp(-50, U_current_surface, d_water):.3f} m/s")
print(f"  海底流速 (指数) = {current_profile_exp(-d_water, U_current_surface, d_water):.4f} m/s")

# 4.2 Morison方程计算海流力
def morison_force_per_length(D, U, U_dot=0, CM=2.0, CD=1.0, rho=1025):
    """
    Morison方程计算单位长度波浪/海流力
    D: 构件直径 (m)
    U: 水质点速度 (m/s)
    U_dot: 水质点加速度 (m/s^2)
    CM: 惯性力系数
    CD: 拖曳力系数
    """
    # 惯性力
    F_I = CM * rho * (np.pi * D**2 / 4) * U_dot
    # 拖曳力
    F_D = 0.5 * CD * rho * D * np.abs(U) * U
    
    return F_I + F_D

# 示例：计算不同直径构件的海流力
D_pile = 2.0  # 桩直径 (m)
F_current_exp = morison_force_per_length(D_pile, U_current_exp, U_dot=0, CM=2.0, CD=1.2)
F_current_linear = morison_force_per_length(D_pile, U_current_linear, U_dot=0, CM=2.0, CD=1.2)

# 绘制海流荷载
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：海流速度剖面
ax1 = axes[0, 0]
ax1.plot(U_current_exp, z_current, 'b-', linewidth=2, label='Exponential Profile')
ax1.plot(U_current_linear, z_current, 'r--', linewidth=2, label='Linear Profile')
ax1.axvline(x=U_current_surface, color='g', linestyle=':', alpha=0.7, label=f'Surface = {U_current_surface} m/s')
ax1.set_xlabel('Current Velocity (m/s)', fontsize=11)
ax1.set_ylabel('Depth z (m)', fontsize=11)
ax1.set_title('Current Velocity Profile', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 图2：海流力分布
ax2 = axes[0, 1]
ax2.plot(F_current_exp/1000, z_current, 'b-', linewidth=2, label='Exponential Profile')
ax2.plot(F_current_linear/1000, z_current, 'r--', linewidth=2, label='Linear Profile')
ax2.set_xlabel('Current Force per Unit Length (kN/m)', fontsize=11)
ax2.set_ylabel('Depth z (m)', fontsize=11)
ax2.set_title(f'Current Force on D={D_pile}m Pile', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 图3：波浪+海流联合作用
ax3 = axes[1, 0]

# 波浪水质点速度 (水面处)
t_combined = np.linspace(0, 3*Tz, 300)
_, u_wave_surface, _, _, _, _, _ = linear_wave(0, 0, t_combined, a_wave, Tz, d_water)

# 海流速度 (恒定)
U_current_const = U_current_surface

# 联合速度
U_combined = u_wave_surface + U_current_const

# 计算拖曳力 (忽略惯性力)
F_wave_only = 0.5 * 1.2 * rho_water * D_pile * np.abs(u_wave_surface) * u_wave_surface
F_combined = 0.5 * 1.2 * rho_water * D_pile * np.abs(U_combined) * U_combined

ax3.plot(t_combined, F_wave_only/1000, 'b-', linewidth=1.5, alpha=0.7, label='Wave Only')
ax3.plot(t_combined, F_combined/1000, 'r-', linewidth=2, label='Wave + Current')
ax3.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax3.set_xlabel('Time (s)', fontsize=11)
ax3.set_ylabel('Drag Force (kN/m)', fontsize=11)
ax3.set_title('Wave+Current Combined Force', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 图4：不同KC数下的拖曳力系数
ax4 = axes[1, 1]

KC_range = np.linspace(0.1, 50, 100)
# 经验公式：CD随KC数变化 (简化模型)
CD_wave = 1.2 * (1 + 0.5 * np.exp(-KC_range/10))
CM_wave = 2.0 - 0.5 * np.exp(-KC_range/5)

ax4.plot(KC_range, CD_wave, 'b-', linewidth=2, label='CD (Drag Coeff.)')
ax4.plot(KC_range, CM_wave, 'r-', linewidth=2, label='CM (Inertia Coeff.)')
ax4.set_xlabel('KC Number', fontsize=11)
ax4.set_ylabel('Coefficient', fontsize=11)
ax4.set_title('Morison Coefficients vs KC Number', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(0, 50)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_current_load.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 海流荷载仿真图已保存")

# ==================== 5. 环境荷载综合展示 ====================
print("\n" + "=" * 60)
print("5. 环境荷载综合展示")
print("=" * 60)

fig, axes = plt.subplots(3, 1, figsize=(14, 12))

# 时间序列
t_total = np.linspace(0, 300, 600)

# 波浪时程
_, eta_total = synthesize_wave_time_series(omega_range, S2, duration=300, dt=0.5, seed=42)
eta_total = eta_total[:len(t_total)]

# 风速时程
_, U_wind_total, _ = synthesize_wind_time_series(f_wind, S_wind, duration=300, dt=0.5, U_mean=U10, seed=123)
U_wind_total = U_wind_total[:len(t_total)]

# 图1：波浪时程
ax1 = axes[0]
ax1.fill_between(t_total, eta_total, alpha=0.5, color='steelblue')
ax1.plot(t_total, eta_total, 'b-', linewidth=1)
ax1.axhline(y=Hs2/2, color='r', linestyle='--', alpha=0.7, label=f'+Hs/2 = {Hs2/2}m')
ax1.axhline(y=-Hs2/2, color='r', linestyle='--', alpha=0.7, label=f'-Hs/2 = -{Hs2/2}m')
ax1.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax1.set_ylabel('Wave Elevation (m)', fontsize=11)
ax1.set_title('Environmental Loads on Offshore Platform', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 300)
ax1.text(0.02, 0.95, f'Hs = {Hs2}m, Tz = {Tz}s', transform=ax1.transAxes, 
         fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# 图2：风速时程
ax2 = axes[1]
ax2.plot(t_total, U_wind_total, 'g-', linewidth=1, alpha=0.8)
ax2.axhline(y=U10, color='r', linestyle='--', alpha=0.7, label=f'Mean = {U10} m/s')
ax2.fill_between(t_total, U10, U_wind_total, alpha=0.3, color='gray')
ax2.set_ylabel('Wind Speed (m/s)', fontsize=11)
ax2.legend(fontsize=9, loc='upper right')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 300)
ax2.text(0.02, 0.95, f'U10 = {U10} m/s', transform=ax2.transAxes, 
         fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))

# 图3：海流剖面
ax3 = axes[2]
ax3.plot(U_current_exp, z_current, 'purple', linewidth=2.5, label='Current Velocity')
ax3.fill_betweenx(z_current, 0, U_current_exp, alpha=0.3, color='purple')
ax3.set_xlabel('Velocity (m/s)', fontsize=11)
ax3.set_ylabel('Depth (m)', fontsize=11)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.text(0.02, 0.95, f'Surface Current = {U_current_surface} m/s', transform=ax3.transAxes, 
         fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='plum', alpha=0.5))

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_environmental_summary.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 环境荷载综合展示图已保存")

# ==================== 生成GIF动画 ====================
print("\n" + "=" * 60)
print("6. 生成波浪传播GIF动画")
print("=" * 60)

from matplotlib.animation import FuncAnimation

# 创建波浪传播动画
fig, ax = plt.subplots(figsize=(12, 5))

x_anim = np.linspace(0, 3*lambda_wave, 300)
line, = ax.plot([], [], 'b-', linewidth=2)
fill = ax.fill_between([], [], alpha=0.3)
ax.set_xlim(0, 3*lambda_wave)
ax.set_ylim(-a_wave*1.5, a_wave*1.5)
ax.set_xlabel('Position x (m)', fontsize=12)
ax.set_ylabel('Wave Elevation (m)', fontsize=12)
ax.set_title('Wave Propagation Animation', fontsize=14, fontweight='bold')
ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax.grid(True, alpha=0.3)

def init():
    line.set_data([], [])
    return line,

def animate(frame):
    t = frame * Tz / 30  # 30帧一个周期
    eta = a_wave * np.cos(k_wave * x_anim - omega_wave * t)
    line.set_data(x_anim, eta)
    
    # 更新填充区域
    for coll in ax.collections[:]:
        coll.remove()
    ax.fill_between(x_anim, 0, eta, alpha=0.3, color='steelblue')
    
    return line,

anim = FuncAnimation(fig, animate, init_func=init, frames=60, interval=100, blit=False)
anim.save(f'{output_dir}/case1_wave_animation.gif', writer='pillow', fps=15, dpi=100)
plt.close()
print("\n✓ 波浪传播GIF动画已保存")

# ==================== 总结 ====================
print("\n" + "=" * 60)
print("案例1仿真完成总结")
print("=" * 60)
print("\n本案例实现了以下环境荷载模型:")
print("1. ✓ 线性波浪理论 (Airy波)")
print("2. ✓ 随机波浪谱 (JONSWAP谱、PM谱)")
print("3. ✓ 风荷载模型 (对数律、幂律、Kaimal谱)")
print("4. ✓ 海流荷载模型 (指数剖面、线性剖面)")
print("5. ✓ Morison方程计算波浪/海流力")
print("6. ✓ 波浪传播GIF动画")
print("\n生成的文件:")
print(f"  - {output_dir}/case1_wave_theory.png")
print(f"  - {output_dir}/case1_wave_spectrum.png")
print(f"  - {output_dir}/case1_wind_load.png")
print(f"  - {output_dir}/case1_current_load.png")
print(f"  - {output_dir}/case1_environmental_summary.png")
print(f"  - {output_dir}/case1_wave_animation.gif")
print("\n" + "=" * 60)

代码说明

1. 参数设置区：定义系统的质量、刚度、阻尼等基本参数
2. 核心计算模块：实现振动方程的求解算法
3. 可视化模块：生成分析图表和动画
4. 结果输出：保存图片文件供文档使用

运行上述代码将生成本主题所需的所有可视化素材。