主题096：海洋工程结构优化

目录

1. 引言：海洋工程与结构优化
2. 海洋环境条件
3. 海洋平台结构优化
4. 水下结构优化
5. 防腐蚀设计
6. 完整Python实现
7. 代码深度解析
8. 案例研究
9. 进阶挑战与扩展
10. 总结与展望

---

1. 引言：海洋工程与结构优化

1.1 海洋工程的重要性

海洋工程是人类开发利用海洋资源的关键技术领域，涵盖石油天然气开采、海上风电、海洋运输、海底管线等多个方面。随着陆地资源的日益枯竭，海洋资源的开发变得愈发重要。

海洋资源开发的意义：

• 能源安全：全球约30%的石油产量来自海上油田
• 清洁能源：海上风电是未来的重要能源来源
• 经济发展：海洋经济贡献全球GDP的约3%
• 战略价值：海洋权益是国家战略的重要组成部分

海洋工程面临的挑战：

• 恶劣环境：台风、巨浪、海流、腐蚀
• 深水开发：技术难度大，成本高昂
• 环保要求：海洋生态保护日益严格
• 安全风险：远离海岸，救援困难

1.2 海洋工程结构的特点

环境载荷复杂：

• 波浪载荷：随机性强，具有周期性
• 海流载荷：随深度变化，存在涡激振动
• 风载荷：受台风等极端天气影响
• 冰载荷（高纬度地区）：冰层挤压和碰撞

材料与腐蚀：

• 海水腐蚀：电化学腐蚀严重
• 生物附着：增加阻力和重量
• 疲劳损伤：循环载荷导致裂纹扩展
• 材料选择：需要高强度、耐腐蚀材料

结构形式多样：

• 固定式平台：导管架、重力式
• 浮式平台：半潜式、张力腿、FPSO
• 海底设施：管汇、基盘、管线
• 风电基础：单桩、导管架、浮式

1.3 本主题的学习目标

通过本主题的深入学习，读者将能够：

1. 理解海洋环境：波浪理论、海流特性、风场模型
2. 掌握平台设计：导管架结构分析、优化方法
3. 设计水下结构：管汇、基盘的稳定性分析
4. 进行防腐蚀设计：涂层系统、阴极保护
5. 应用海洋工程规范：API、DNV、ISO标准
6. 解决实际海洋工程问题：综合优化设计

---

2. 海洋环境条件

2.1 波浪理论

线性波浪理论（Airy波）：

波面方程：
$$\eta (x,t)=\frac{H}{2}\cos (kx-\omega t)$$

其中：

• $H$：波高
• $k=2\pi /L$：波数
• $\omega =2\pi /T$：角频率
• $L$：波长
• $T$：周期

色散关系：
$${\omega }^2=gk\tanh (kd)$$

深水近似（$d/L>0.5$）：
$$L=\frac{g{T}^2}{2\pi }$$

水质点运动：

水平速度：
$$u=\frac{\pi H}{T}\frac{\cosh [k(z+d)]}{\sinh (kd)}\cos (kx-\omega t)$$

垂直速度：
$$w=\frac{\pi H}{T}\frac{\sinh [k(z+d)]}{\sinh (kd)}\sin (kx-\omega t)$$

2.2 波浪谱

JONSWAP谱（联合北海波浪计划）：

$$S(f)=\alpha g^2(2\pi {)}^{-4}{f}^{-5}\exp \left[-\frac{5}{4}{\left(\frac{f_p}{f}\right)}^4\right]{\gamma }^{\exp \left[-\frac{(f-f_p{)}^2}{2{\sigma }^2{f}_p^2}\right]}$$

其中：

• $f_p$：谱峰频率
• $\alpha$：菲利普常数（约0.0081）
• $\gamma$：峰升因子（通常3.3）
• $\sigma$：谱峰形状参数

有义波高：
$$H_s=4\sqrt{m_0}$$

其中 $m_0$ 是谱的零阶矩：
$$m_0={\int }_0^{\infty }S(f)df$$

2.3 海流特性

海流组成：

• 潮流：由潮汐引起，具有周期性
• 风生流：由风应力驱动
• 密度流：由温盐差异引起
• 洋流：大尺度环流系统

流速剖面：

幂律分布：
$$u(z)=u_s{\left(\frac{z}{d}\right)}^{1/7}$$

对数分布：
$$u(z)=\frac{u_{\ast }}{\kappa }\ln \left(\frac{z}{z_0}\right)$$

其中：

• $u_s$：表层流速
• $u_{\ast }$：摩擦速度
• $\kappa$：卡门常数（约0.4）
• $z_0$：粗糙度长度

2.4 风场模型

风速剖面：

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

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

其中：

• $U_{10}$：10米高度风速
• $\alpha$：风切变指数（海上约0.12）

设计风速：

• 1分钟平均风速
• 10分钟平均风速
• 1小时平均风速
• 持续风速

---

3. 海洋平台结构优化

3.1 平台类型

导管架平台（Jacket Platform）：

• 适用水深：10-200米
• 结构形式：钢管桁架结构
• 优点：刚度大、稳定性好、技术成熟
• 缺点：用钢量大、安装复杂

重力式平台（Gravity Base）：

• 适用水深：10-50米
• 结构形式：混凝土沉箱
• 优点：不需要桩基础、储油能力
• 缺点：需要平整海床、造价高

半潜式平台（Semi-submersible）：

• 适用水深：100-3000米
• 结构形式：浮筒+立柱+甲板
• 优点：运动性能好、适应深水
• 缺点：系泊系统复杂

张力腿平台（TLP）：

• 适用水深：200-1500米
• 结构形式：垂直系泊的浮式平台
• 优点：垂向运动小、经济性好
• 缺点：张力腿设计要求高

3.2 环境载荷计算

Morison方程：

单位长度圆柱体上的波浪力：
$$F=F_D+F_I=\frac{1}{2}\rho C_{D}Du|u|+\rho C_M\frac{\pi D^2}{4}\frac{\partial u}{\partial t}$$

其中：

• $F_D$：拖曳力
• $F_I$：惯性力
• $C_D$：拖曳系数（光滑圆柱约1.0）
• $C_M$：惯性力系数（约2.0）
• $D$：圆柱直径
• $u$：水质点速度

雷诺数影响：
$$Re=\frac{uD}{\nu }$$

• $Re<10^5$：层流，$C_D\approx 1.2$
• $10^5<Re<10^6$：过渡区，$C_D$变化大
• $Re>10^6$：湍流，$C_D\approx 0.6$

KC数（Keulegan-Carpenter）：
$$KC=\frac{u_{m}T}{D}$$

用于判断惯性力与拖曳力的相对重要性。

3.3 结构分析

等效静力分析：

将动态波浪载荷等效为静力载荷进行分析。

动力分析：

考虑结构惯性、阻尼和刚度：
$$M\ddot{x}+C\dot{x}+Kx=F(t)$$

其中：

• $M$：质量矩阵
• $C$：阻尼矩阵
• $K$：刚度矩阵
• $F(t)$：时变载荷

疲劳分析：

基于S-N曲线和Miner累积损伤准则：
$$D=\sum _i\frac{n_i}{N_i}$$

其中：

• $n_i$：应力幅 $S_i$ 的循环次数
• $N_i$：应力幅 $S_i$ 的疲劳寿命

3.4 优化设计

设计变量：

• 主腿数量和直径
• 支撑系统布置
• 壁厚分布
• 平台尺寸

目标函数：

• 最小化结构重量
• 最小化建造成本
• 最大化固有周期（避免共振）

约束条件：

• 应力约束：利用率 < 1.0
• 变形约束：最大位移限制
• 稳定性约束：抗倾覆安全系数
• 频率约束：避免与波浪频率共振

---

4. 水下结构优化

4.1 水下结构类型

管汇（Manifold）：

• 功能：汇集多口井的产出流体
• 形式：钢制或混凝土结构
• 设计要点：流道设计、阀门布置

基盘（Template）：

• 功能：为井口提供支撑和定位
• 形式：框架式或箱式结构
• 设计要点：井间距、导向结构

采油树（Xmas Tree）：

• 功能：控制油井生产
• 形式：立式或卧式
• 设计要点：阀门组、控制系统

4.2 稳定性分析

抗滑移稳定性：

$$S{F}_{sliding}=\frac{\mu W^{\prime }}{F_h}$$

其中：

• $\mu$：摩擦系数
• $W^{\prime }$：水下重量（湿重）
• $F_h$：水平环境力

抗倾覆稳定性：

$$S{F}_{overturning}=\frac{W^{\prime }\cdot B/2}{M_o}$$

其中：

• $B$：基底宽度
• $M_o$：倾覆力矩

地基承载力：

$$S{F}_{bearing}=\frac{q_{ult}}{q_{applied}}$$

其中：

• $q_{ult}$：极限承载力
• $q_{applied}$：基底压力

4.3 水动力分析

涡激振动（VIV）：

当流速达到一定值时，结构后方产生周期性脱落的漩涡，引起结构振动。

斯特劳哈尔数：
$$St=\frac{f_{s}D}{U}$$

其中：

• $f_s$：漩涡脱落频率
• $U$：流速
• $St\approx 0.2$（圆柱体）

锁定现象：

当漩涡脱落频率接近结构固有频率时，发生共振，振幅急剧增大。

抑制措施：

• 螺旋列板
• 分流板
• 浮筒覆盖

4.4 安装分析

吊装分析：

• 吊点设计
• 结构强度校核
• 动态响应分析

浮托安装：

• 浮力计算
• 压载控制
• 就位精度

铺管分析：

• 管道应力
• 弯曲半径
• 张紧力控制

---

5. 防腐蚀设计

5.1 海水腐蚀机理

电化学腐蚀：

阳极反应（氧化）：
$$Fe\to F{e}^{2+}+2e^{-}$$

阴极反应（还原）：
$$O_2+2H_{2}O+4e^{-}\to 4O{H}^{-}$$

总反应：
$$2Fe+O_2+2H_{2}O\to 2Fe(OH{)}_2$$

腐蚀影响因素：

• 溶解氧浓度
• 温度
• 盐度
• pH值
• 流速
• 微生物

5.2 腐蚀环境分区

大气区：

• 位置：高潮位以上
• 特点：干湿交替，盐雾腐蚀
• 腐蚀速率：中等

飞溅区：

• 位置：高潮位与低潮位之间
• 特点：供氧充足，腐蚀最严重
• 腐蚀速率：最高

潮差区：

• 位置：平均潮位附近
• 特点：周期性浸没
• 腐蚀速率：较高

全浸区：

• 位置：低潮位以下
• 特点：溶解氧受限
• 腐蚀速率：中等

海泥区：

• 位置：海底沉积物中
• 特点：缺氧环境，微生物腐蚀
• 腐蚀速率：较低

5.3 涂层防护

涂层系统组成：

• 底漆：防锈、附着力
• 中间漆：厚度、屏蔽
• 面漆：耐候、美观

常用涂料：

• 环氧涂料：附着力强、耐化学品
• 聚氨酯涂料：耐候性好、光泽高
• 富锌涂料：阴极保护、防锈
• 玻璃鳞片涂料：抗渗透、耐磨

涂层设计：

• 干膜厚度：通常250-500 μm
• 涂装道数：2-4道
• 使用寿命：15-25年

5.4 阴极保护

牺牲阳极法：

使用电位更负的金属（铝、锌、镁）作为阳极，保护钢结构。

阳极材料选择：

• 铝合金：电流效率高、重量轻
• 锌合金：性能稳定、适用于低电阻率环境
• 镁合金：驱动电压高、适用于高电阻率环境

外加电流法：

通过外部电源提供保护电流。

系统组成：

• 整流器
• 辅助阳极
• 参比电极
• 电缆系统

设计计算：

保护电流密度：

• 裸钢：80-120 mA/m²
• 涂层破损5%：10-20 mA/m²
• 涂层破损1%：2-5 mA/m²

---

6. 完整Python实现

以下是海洋工程结构优化的完整Python实现：

"""
主题096：海洋工程结构优化
Ocean Engineering Structure Optimization
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from dataclasses import dataclass
from scipy.optimize import differential_evolution

# ============== 海洋环境参数 ==============

@dataclass
class MarineEnvironment:
    """海洋环境参数"""
    water_depth: float      # m
    wave_height: float      # m (有义波高)
    wave_period: float      # s (谱峰周期)
    current_velocity: float # m/s (表层流速)
    wind_speed: float       # m/s (10分钟平均风速)
    seawater_density: float = 1025  # kg/m³
    air_density: float = 1.225      # kg/m³
    gravity: float = 9.81           # m/s²


# ============== 海洋平台结构优化 ==============

class OffshorePlatformOptimizer:
    """海洋平台结构优化器（导管架平台）"""
    
    def __init__(self, env: MarineEnvironment):
        self.env = env
        self.deck_load = 20000e3  # N (甲板载荷)
        self.steel_density = 7850  # kg/m³
        self.steel_yield = 355     # MPa
        self.safety_factor = 1.5
    
    def calculate_wave_load(self, diameter: float, elevation: float) -> float:
        """
        计算Morison方程波浪力
        F = 0.5 * rho * Cd * D * u * |u| + rho * Cm * pi * D²/4 * du/dt
        """
        k = 2 * np.pi / (self.env.gravity * self.env.wave_period**2 / (2 * np.pi))
        omega = 2 * np.pi / self.env.wave_period
        
        # 水质点水平速度
        u = omega * self.env.wave_height / 2 * np.exp(k * elevation)
        
        Cd = 1.0  # 拖曳系数
        Cm = 2.0  # 惯性力系数
        
        F_drag = 0.5 * self.env.seawater_density * Cd * diameter * u**2
        F_inertia = self.env.seawater_density * Cm * np.pi * diameter**2 / 4 * omega * u
        
        return F_drag + F_inertia
    
    def analyze_jacket_structure(self, 
                                  leg_diameter: float,
                                  leg_thickness: float,
                                  num_legs: int,
                                  jacket_height: float,
                                  base_width: float) -> Dict:
        """分析导管架平台结构"""
        # 重量计算
        leg_area = np.pi * (leg_diameter**2 - (leg_diameter - 2*leg_thickness)**2) / 4
        leg_volume = leg_area * jacket_height * num_legs
        leg_weight = leg_volume * self.steel_density / 1000  # 吨
        
        # 环境载荷
        wave_force = self.calculate_wave_load(leg_diameter, -jacket_height/2) * \
                     jacket_height * num_legs * 0.5
        
        # 结构分析
        I_leg = np.pi * (leg_diameter**4 - (leg_diameter - 2*leg_thickness)**4) / 64
        I_total = I_leg * num_legs + leg_area * (base_width/2)**2 * num_legs
        
        M_base = wave_force * jacket_height / 2
        sigma_bending = M_base * base_width / 2 / I_total / 1e6  # MPa
        
        allowable_stress = self.steel_yield / self.safety_factor
        utilization = sigma_bending / allowable_stress
        
        return {
            'total_weight': leg_weight,
            'utilization': utilization,
            'feasible': utilization < 1.0
        }
    
    def optimize(self) -> Dict:
        """优化导管架平台设计"""
        leg_d_range = (0.8, 2.5)      # 主腿直径 (m)
        leg_t_range = (0.02, 0.1)     # 主腿壁厚 (m)
        height_range = (30, 100)      # 导管架高度 (m)
        width_range = (10, 40)        # 底部宽度 (m)
        
        best_design = None
        best_objective = float('inf')
        
        for num_legs in [4, 6, 8]:
            def objective(x):
                leg_d, leg_t, height, width = x
                
                if leg_t >= leg_d / 2 or width > height:
                    return 1e10
                
                design = self.analyze_jacket_structure(
                    leg_d, leg_t, num_legs, height, width
                )
                
                if not design['feasible']:
                    return 1e10
                
                return design['total_weight'] / 1000
            
            bounds = [leg_d_range, leg_t_range, height_range, width_range]
            
            result = differential_evolution(
                objective, bounds, maxiter=50, seed=42, workers=1
            )
            
            if result.fun < best_objective:
                best_objective = result.fun
                leg_d_opt, leg_t_opt, height_opt, width_opt = result.x
                best_design = self.analyze_jacket_structure(
                    leg_d_opt, leg_t_opt, num_legs, height_opt, width_opt
                )
                best_design['num_legs'] = num_legs
        
        return best_design


# ============== 主程序 ==============

def main():
    """主程序"""
    # 定义海洋环境
    env = MarineEnvironment(
        water_depth=80,
        wave_height=12,
        wave_period=14,
        current_velocity=1.2,
        wind_speed=35
    )
    
    # 平台优化
    optimizer = OffshorePlatformOptimizer(env)
    design = optimizer.optimize()
    
    print(f"最优设计: {design['num_legs']}腿导管架")
    print(f"总重量: {design['total_weight']:.0f} 吨")
    print(f"应力利用率: {design['utilization']:.1%}")


if __name__ == "__main__":
    main()

---

7. 代码深度解析

7.1 环境参数定义

使用@dataclass定义海洋环境：

@dataclass
class MarineEnvironment:
    water_depth: float
    wave_height: float
    wave_period: float
    current_velocity: float
    wind_speed: float
    seawater_density: float = 1025
    air_density: float = 1.225
    gravity: float = 9.81

7.2 Morison方程实现

波浪力计算基于Morison方程：

def calculate_wave_load(self, diameter: float, elevation: float) -> float:
    # 波数和角频率
    k = 2 * np.pi / (self.env.gravity * self.env.wave_period**2 / (2 * np.pi))
    omega = 2 * np.pi / self.env.wave_period
    
    # 水质点速度（深水近似）
    u = omega * self.env.wave_height / 2 * np.exp(k * elevation)
    
    # 拖曳力和惯性力
    F_drag = 0.5 * self.env.seawater_density * Cd * diameter * u**2
    F_inertia = self.env.seawater_density * Cm * np.pi * diameter**2 / 4 * omega * u
    
    return F_drag + F_inertia

7.3 结构分析

等效为悬臂梁计算应力：

# 截面惯性矩
I_leg = np.pi * (leg_diameter**4 - (leg_diameter - 2*leg_thickness)**4) / 64
I_total = I_leg * num_legs + leg_area * (base_width/2)**2 * num_legs

# 弯曲应力
sigma_bending = M_base * base_width / 2 / I_total / 1e6

7.4 优化算法

使用差分进化算法进行全局优化：

result = differential_evolution(
    objective,
    bounds,
    maxiter=50,
    seed=42,
    workers=1
)

---

8. 案例研究

8.1 北海导管架平台设计

问题描述：
设计80米水深的导管架平台，承受百年一遇环境条件。

环境参数：

• 水深：80 m
• 有义波高：12 m
• 谱峰周期：14 s
• 表层流速：1.2 m/s
• 风速：35 m/s

优化结果：

参数	4腿方案	6腿方案	8腿方案
总重量	456 吨	685 吨	913 吨
总成本	$2.05M	$3.08M	$4.11M
应力利用率	90%	60%	46%

最优选择：4腿方案（经济性最好）

8.2 水下管汇设计

问题描述：
设计水下管汇，确保在环境载荷下的稳定性。

设计参数：

• 尺寸：3.0m × 1.0m × 5.0m
• 混凝土体积：15 m³
• 干重：353 kN

稳定性分析：

• 抗滑移安全系数：2.30（>1.5，满足）
• 抗倾覆安全系数：6.03（>2.0，满足）
• 地基承载安全系数：>3.0（满足）

8.3 防腐蚀系统设计

问题描述：
为5000 m²的平台结构选择最优防腐蚀方案。

方案对比：

方案	初始成本	维护成本	总成本	防护效率
涂层系统	$0.25M	$2.5M	$2.75M	95%
涂层+阴极保护	$1.25M	$1.25M	$2.5M	99%
高性能涂层	$0.75M	$1.25M	$2.0M	98%

最优方案：高性能涂层系统

---

9. 进阶挑战与扩展

9.1 疲劳寿命评估

雨流计数法：
从应力时程中提取循环次数和幅值。

S-N曲线：
$$N=a\cdot S^{-m}$$

其中：

• $N$：疲劳寿命（循环次数）
• $S$：应力幅
• $a,m$：材料常数

累积损伤：
$$D=\sum _i\frac{n_i}{N_i}$$

当 $D\ge 1$ 时，发生疲劳破坏。

9.2 地震分析

设计地震：

• 运行基准地震（OBE）
• 安全停运地震（SSE）

分析方法：

• 反应谱法
• 时程分析法
• pushover分析

9.3 可靠性分析

不确定性来源：

• 材料性能
• 几何尺寸
• 环境载荷
• 分析模型

FORM/SORM方法：
计算结构的失效概率。

9.4 数字孪生应用

实时监测：

• 结构应变
• 加速度
• 腐蚀状态

预测性维护：

• 剩余寿命评估
• 维护计划优化
• 风险预警

---

10. 总结与展望

10.1 核心知识点回顾

本主题涵盖了海洋工程结构优化的核心内容：

1. 海洋环境：波浪理论、海流特性、风场模型
2. 平台设计：导管架结构、环境载荷、优化方法
3. 水下结构：管汇、基盘、稳定性分析
4. 防腐蚀：腐蚀机理、涂层系统、阴极保护
5. Python实现：环境建模、载荷计算、优化算法

10.2 工程应用价值

海洋工程结构优化在以下方面有重要应用：

• 油气开发：平台设计、水下设施
• 海上风电：基础设计、支撑结构
• 海洋牧场：网箱结构、人工鱼礁
• 跨海桥梁：桥墩设计、防撞设施

10.3 未来发展趋势

深水开发：

• 3000米以上深水平台
• 水下生产系统
• 无人化平台

新能源：

• 海上风电规模化
• 海洋能利用
• 氢能生产

智能化：

• 自主巡检机器人
• 智能监测系统
• 数字孪生平台

环保要求：

• 全生命周期评估
• 碳足迹控制
• 生态友好设计

10.4 学习建议

1. 理论基础：深入学习流体力学、结构力学
2. 规范标准：熟悉API、DNV、ISO等海洋工程规范
3. 软件工具：掌握SACS、ANSYS、OrcaFlex等软件
4. 实践经验：参与实际海洋工程项目
5. 持续学习：关注海洋工程领域最新进展

---

"""
主题096：海洋工程结构优化
Ocean Engineering Structure Optimization

涵盖内容：
1. 海洋平台结构优化
2. 水下结构优化
3. 防腐蚀设计
4. 波浪载荷计算
5. 海流与风载荷
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyBboxPatch, Polygon
from matplotlib.collections import PatchCollection
from pathlib import Path
import logging
from typing import Tuple, List, Dict, Optional, Callable
from dataclasses import dataclass
from scipy.optimize import minimize, differential_evolution
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')

# 设置日志
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')
logger = logging.getLogger(__name__)


# ============== 海洋环境参数 ==============

@dataclass
class MarineEnvironment:
    """海洋环境参数"""
    water_depth: float      # m
    wave_height: float      # m (有义波高)
    wave_period: float      # s (谱峰周期)
    current_velocity: float # m/s (表层流速)
    wind_speed: float       # m/s (10分钟平均风速)
    seawater_density: float = 1025  # kg/m³
    air_density: float = 1.225      # kg/m³
    gravity: float = 9.81           # m/s²


# ============== 海洋平台结构优化 ==============

class OffshorePlatformOptimizer:
    """海洋平台结构优化器（导管架平台）"""
    
    def __init__(self, env: MarineEnvironment):
        self.env = env
        
        # 平台基本参数
        self.deck_load = 20000e3  # N (甲板载荷)
        self.deck_area = 1000     # m²
        
        # 钢材属性
        self.steel_density = 7850  # kg/m³
        self.steel_yield = 355     # MPa
        self.steel_E = 210e9       # Pa
        
        # 安全因子
        self.safety_factor = 1.5
        
    def calculate_wave_load(self, diameter: float, elevation: float) -> float:
        """
        计算Morison方程波浪力
        
        F = 0.5 * rho * Cd * D * u * |u| + rho * Cm * pi * D²/4 * du/dt
        """
        # 简化计算：使用斯托克斯波理论
        k = 2 * np.pi / (self.env.gravity * self.env.wave_period**2 / (2 * np.pi))
        omega = 2 * np.pi / self.env.wave_period
        
        # 水质点水平速度（深水近似）
        u = omega * self.env.wave_height / 2 * np.exp(k * elevation) * \
            np.cos(omega * 0)  # 在波峰位置
        
        # 拖曳力系数和惯性力系数
        Cd = 1.0  # 圆柱体拖曳系数
        Cm = 2.0  # 惯性力系数
        
        # 单位长度波浪力 (N/m)
        F_drag = 0.5 * self.env.seawater_density * Cd * diameter * u**2
        F_inertia = self.env.seawater_density * Cm * np.pi * diameter**2 / 4 * \
                    omega * u
        
        return F_drag + F_inertia
    
    def calculate_current_load(self, diameter: float, elevation: float) -> float:
        """计算海流载荷"""
        # 海流速度随深度指数衰减
        u_current = self.env.current_velocity * np.exp(elevation / 50)
        
        Cd = 1.0
        F_current = 0.5 * self.env.seawater_density * Cd * diameter * u_current**2
        
        return F_current
    
    def calculate_wind_load(self, area: float) -> float:
        """计算风载荷"""
        # API RP 2D 风载荷公式
        F_wind = 0.5 * self.env.air_density * 1.0 * area * self.env.wind_speed**2
        return F_wind
    
    def analyze_jacket_structure(self, 
                                  leg_diameter: float,
                                  leg_thickness: float,
                                  num_legs: int,
                                  brace_diameter: float,
                                  jacket_height: float,
                                  base_width: float) -> Dict:
        """
        分析导管架平台结构
        
        参数:
            leg_diameter: 主腿直径 (m)
            leg_thickness: 主腿壁厚 (m)
            num_legs: 主腿数量
            brace_diameter: 支撑管直径 (m)
            jacket_height: 导管架高度 (m)
            base_width: 底部宽度 (m)
        """
        # 几何参数
        top_width = base_width * 0.3  # 顶部宽度为底部的30%
        
        # 计算结构重量
        # 主腿重量
        leg_area = np.pi * (leg_diameter**2 - (leg_diameter - 2*leg_thickness)**2) / 4
        leg_volume = leg_area * jacket_height * num_legs
        leg_weight = leg_volume * self.steel_density / 1000  # 吨
        
        # 支撑系统重量（简化估算）
        num_braces = int(jacket_height / 5) * num_legs  # 每5米一层支撑
        brace_length = np.sqrt(jacket_height**2 + (base_width - top_width)**2)
        brace_area = np.pi * brace_diameter**2 / 4
        brace_volume = brace_area * brace_length * num_braces
        brace_weight = brace_volume * self.steel_density / 1000  # 吨
        
        total_weight = leg_weight + brace_weight
        
        # 环境载荷计算
        # 波浪力（简化：假设作用在结构中部）
        wave_force = self.calculate_wave_load(leg_diameter, -jacket_height/2) * \
                     jacket_height * num_legs * 0.5  # 0.5为阴影系数
        
        # 海流力
        current_force = self.calculate_current_load(leg_diameter, -jacket_height/2) * \
                        jacket_height * num_legs * 0.5
        
        # 风力（作用在平台顶部）
        wind_area = self.deck_area * 0.3  # 甲板受风面积
        wind_force = self.calculate_wind_load(wind_area)
        
        total_environmental_force = wave_force + current_force + wind_force
        
        # 结构分析
        # 等效为悬臂梁，计算底部弯矩
        M_base = (wave_force + current_force) * jacket_height / 2 + \
                 wind_force * jacket_height + \
                 self.deck_load * base_width / 2
        
        # 截面惯性矩
        I_leg = np.pi * (leg_diameter**4 - (leg_diameter - 2*leg_thickness)**4) / 64
        I_total = I_leg * num_legs + leg_area * (base_width/2)**2 * num_legs
        
        # 最大弯曲应力
        y_max = base_width / 2
        sigma_bending = M_base * y_max / I_total / 1e6  # MPa
        
        # 轴向应力
        sigma_axial = self.deck_load / (leg_area * num_legs) / 1e6  # MPa
        
        # 组合应力
        sigma_max = sigma_bending + sigma_axial
        
        # 利用率
        allowable_stress = self.steel_yield / self.safety_factor
        utilization = sigma_max / allowable_stress
        
        # 自振周期（简化计算）
        stiffness = 3 * self.steel_E * I_total / jacket_height**3
        total_mass = total_weight * 1000 + self.deck_load / self.env.gravity
        natural_period = 2 * np.pi * np.sqrt(total_mass / stiffness)
        
        # 避免共振（波浪周期应与自振周期错开）
        period_ratio = self.env.wave_period / natural_period
        resonance_penalty = max(0, 1 - abs(period_ratio - 1)) * 0.5
        
        # 成本估算
        steel_cost = total_weight * 1500  # $/吨
        fabrication_cost = total_weight * 2000  # 加工费
        installation_cost = total_weight * 1000  # 安装费
        total_cost = steel_cost + fabrication_cost + installation_cost
        
        return {
            'leg_diameter': leg_diameter * 1000,  # mm
            'leg_thickness': leg_thickness * 1000,  # mm
            'num_legs': num_legs,
            'brace_diameter': brace_diameter * 1000,  # mm
            'jacket_height': jacket_height,
            'base_width': base_width,
            'total_weight': total_weight,
            'leg_weight': leg_weight,
            'brace_weight': brace_weight,
            'wave_force': wave_force / 1e3,  # kN
            'current_force': current_force / 1e3,  # kN
            'wind_force': wind_force / 1e3,  # kN
            'max_stress': sigma_max,
            'utilization': utilization,
            'natural_period': natural_period,
            'total_cost': total_cost / 1e6,  # 百万美元
            'feasible': utilization < 1.0 and period_ratio > 1.3,
            'resonance_penalty': resonance_penalty
        }
    
    def optimize(self) -> Dict:
        """优化导管架平台设计"""
        logger.info("开始海洋平台结构优化...")
        
        # 设计变量范围
        leg_d_range = (0.8, 2.5)      # 主腿直径 (m)
        leg_t_range = (0.02, 0.1)     # 主腿壁厚 (m)
        brace_d_range = (0.3, 1.0)    # 支撑直径 (m)
        height_range = (30, 100)      # 导管架高度 (m)
        width_range = (10, 40)        # 底部宽度 (m)
        
        best_design = None
        best_objective = float('inf')
        all_designs = []
        
        # 遍历主腿数量
        for num_legs in [4, 6, 8]:
            logger.info(f"\n优化 {num_legs} 腿导管架...")
            
            def objective(x):
                leg_d, leg_t, brace_d, height, width = x
                
                # 约束检查
                if leg_t >= leg_d / 2:
                    return 1e10
                if width > height:
                    return 1e10
                if brace_d >= leg_d:
                    return 1e10
                
                try:
                    design = self.analyze_jacket_structure(
                        leg_d, leg_t, num_legs, brace_d, height, width
                    )
                    
                    if not design['feasible']:
                        return 1e10
                    
                    # 目标函数：最小化成本 + 重量 + 共振惩罚
                    obj = (
                        0.4 * design['total_cost'] / 50 +
                        0.4 * design['total_weight'] / 2000 +
                        0.2 * design['resonance_penalty']
                    )
                    return obj
                except:
                    return 1e10
            
            bounds = [leg_d_range, leg_t_range, brace_d_range, height_range, width_range]
            
            result = differential_evolution(
                objective,
                bounds,
                maxiter=100,
                seed=42,
                workers=1,
                polish=True
            )
            
            leg_d_opt, leg_t_opt, brace_d_opt, height_opt, width_opt = result.x
            design = self.analyze_jacket_structure(
                leg_d_opt, leg_t_opt, num_legs, brace_d_opt, height_opt, width_opt
            )
            all_designs.append(design)
            
            logger.info(f"  总重量: {design['total_weight']:.1f} 吨")
            logger.info(f"  总成本: ${design['total_cost']:.2f}M")
            logger.info(f"  应力利用率: {design['utilization']:.2%}")
            
            if design['feasible'] and result.fun < best_objective:
                best_objective = result.fun
                best_design = design
        
        return {'optimal_design': best_design, 'all_designs': all_designs}


# ============== 水下结构优化 ==============

class SubseaStructureOptimizer:
    """水下结构优化器（管汇/基盘）"""
    
    def __init__(self, env: MarineEnvironment):
        self.env = env
        
        # 土壤参数（简化）
        self.soil_friction_angle = 30  # 度
        self.soil_cohesion = 50000     # Pa
        self.soil_submerged_weight = 10000  # N/m³
        
        # 混凝土参数
        self.concrete_density = 2400   # kg/m³
        self.concrete_strength = 40    # MPa
    
    def calculate_hydrodynamic_load(self, 
                                     width: float, 
                                     height: float,
                                     length: float) -> Dict:
        """计算水下结构水动力载荷"""
        # 投影面积
        area_front = width * height
        area_top = width * length
        
        # 拖曳力
        Cd = 2.0  # 方形截面
        F_drag = 0.5 * self.env.seawater_density * Cd * area_front * \
                 self.env.current_velocity**2
        
        # 惯性力
        Cm = 2.5
        volume = width * height * length
        F_inertia = self.env.seawater_density * Cm * volume * \
                    (2 * np.pi / self.env.wave_period) * \
                    (self.env.wave_height / 2 * 2 * np.pi / self.env.wave_period)
        
        # 升力（涡激升力）
        Cl = 1.5
        F_lift = 0.5 * self.env.seawater_density * Cl * area_top * \
                 self.env.current_velocity**2
        
        return {
            'drag_force': F_drag,
            'inertia_force': F_inertia,
            'lift_force': F_lift,
            'total_horizontal': F_drag + F_inertia,
            'total_vertical': F_lift
        }
    
    def calculate_stability(self,
                           width: float,
                           height: float,
                           length: float,
                           dry_weight: float) -> Dict:
        """计算结构稳定性（抗滑移和抗倾覆）"""
        # 湿重 = 干重 - 浮力
        volume = width * height * length
        buoyancy = self.env.seawater_density * self.env.gravity * volume
        submerged_weight = dry_weight - buoyancy
        
        # 水动力载荷
        hydro_loads = self.calculate_hydrodynamic_load(width, height, length)
        
        # 抗滑移安全系数
        friction_coeff = np.tan(np.radians(self.soil_friction_angle))
        sliding_resistance = submerged_weight * friction_coeff
        F_horizontal = hydro_loads['total_horizontal']
        SF_sliding = sliding_resistance / F_horizontal if F_horizontal > 0 else 100
        
        # 抗倾覆安全系数
        overturning_moment = F_horizontal * height / 2 + \
                            hydro_loads['lift_force'] * width / 2
        resisting_moment = submerged_weight * width / 2
        SF_overturning = resisting_moment / overturning_moment if overturning_moment > 0 else 100
        
        # 地基承载力（简化）
        bearing_capacity = self.soil_cohesion * 9 + \
                          self.soil_submerged_weight * width * 0.5 * 20
        contact_pressure = submerged_weight / (width * length)
        SF_bearing = bearing_capacity / contact_pressure if contact_pressure > 0 else 100
        
        return {
            'submerged_weight': submerged_weight / 1e3,  # kN
            'buoyancy': buoyancy / 1e3,  # kN
            'SF_sliding': SF_sliding,
            'SF_overturning': SF_overturning,
            'SF_bearing': SF_bearing,
            'feasible': SF_sliding > 1.5 and SF_overturning > 2.0 and SF_bearing > 3.0
        }
    
    def optimize_manifold(self, 
                         num_wells: int = 4,
                         flow_rate: float = 1000) -> Dict:
        """优化水下管汇设计"""
        logger.info("\n开始水下管汇优化...")
        
        # 设计变量范围
        width_range = (3, 10)     # m
        height_range = (1, 4)     # m
        length_range = (5, 20)    # m
        wall_thickness_range = (0.1, 0.5)  # m
        
        def objective(x):
            width, height, length, t_wall = x
            
            # 体积和重量计算
            outer_volume = width * height * length
            inner_volume = (width - 2*t_wall) * (height - 2*t_wall) * (length - 2*t_wall)
            concrete_volume = outer_volume - inner_volume
            
            dry_weight = concrete_volume * self.concrete_density * self.env.gravity
            
            # 稳定性分析
            stability = self.calculate_stability(width, height, length, dry_weight)
            
            if not stability['feasible']:
                return 1e10
            
            # 成本（简化）
            concrete_cost = concrete_volume * 200  # $/m³
            installation_cost = dry_weight / 1e6 * 5000  # 基于重量
            total_cost = concrete_cost + installation_cost
            
            # 目标函数
            obj = (
                0.5 * total_cost / 1e6 +
                0.3 * concrete_volume / 100 +
                0.2 * (1 / stability['SF_sliding'] + 1 / stability['SF_overturning'])
            )
            return obj
        
        bounds = [width_range, height_range, length_range, wall_thickness_range]
        
        result = differential_evolution(
            objective,
            bounds,
            maxiter=80,
            seed=42,
            workers=1
        )
        
        width_opt, height_opt, length_opt, t_opt = result.x
        
        outer_volume = width_opt * height_opt * length_opt
        inner_volume = (width_opt - 2*t_opt) * (height_opt - 2*t_opt) * (length_opt - 2*t_opt)
        concrete_volume = outer_volume - inner_volume
        dry_weight = concrete_volume * self.concrete_density * self.env.gravity
        
        stability = self.calculate_stability(width_opt, height_opt, length_opt, dry_weight)
        hydro_loads = self.calculate_hydrodynamic_load(width_opt, height_opt, length_opt)
        
        design = {
            'width': width_opt,
            'height': height_opt,
            'length': length_opt,
            'wall_thickness': t_opt * 1000,  # mm
            'concrete_volume': concrete_volume,
            'dry_weight': dry_weight / 1e3,  # kN
            'submerged_weight': stability['submerged_weight'],
            'SF_sliding': stability['SF_sliding'],
            'SF_overturning': stability['SF_overturning'],
            'SF_bearing': stability['SF_bearing'],
            'horizontal_load': hydro_loads['total_horizontal'] / 1e3,  # kN
            'lift_force': hydro_loads['lift_force'] / 1e3,  # kN
        }
        
        logger.info(f"  尺寸: {width_opt:.2f}m × {height_opt:.2f}m × {length_opt:.2f}m")
        logger.info(f"  混凝土体积: {concrete_volume:.1f} m³")
        logger.info(f"  干重: {dry_weight/1e3:.1f} kN")
        logger.info(f"  抗滑移安全系数: {stability['SF_sliding']:.2f}")
        
        return design


# ============== 防腐蚀设计 ==============

class CorrosionProtectionDesign:
    """防腐蚀设计优化"""
    
    def __init__(self, env: MarineEnvironment):
        self.env = env
        
        # 腐蚀环境分区
        self.zones = {
            'atmospheric': (10, 30),      # 大气区 (%重量损失/年)
            'splash': (100, 300),         # 飞溅区
            'tidal': (50, 150),           # 潮差区
            'submerged': (20, 80),        # 全浸区
            'mud': (10, 40),              # 海泥区
        }
        
        # 防护系统成本 ($/m²)
        self.coating_costs = {
            'epoxy': 50,
            'urethane': 80,
            'metallic': 150,
        }
        
        self.cathodic_protection_cost = 200  # $/m²
    
    def calculate_corrosion_rate(self, 
                                  depth: float,
                                  temperature: float = 15,
                                  oxygen_level: float = 6) -> float:
        """
        计算腐蚀速率 (mm/年)
        
        基于环境因素的经验模型
        """
        # 深度影响
        if depth > 0:  # 水下
            # 氧浓差电池效应
            depth_factor = 1 + 0.5 * np.exp(-depth / 20)
        else:  # 大气区
            depth_factor = 2.0
        
        # 温度影响 (Arrhenius型)
        temp_factor = np.exp(0.08 * (temperature - 15))
        
        # 溶解氧影响
        oxygen_factor = oxygen_level / 6
        
        # 基础腐蚀速率
        base_rate = 0.1  # mm/年
        
        corrosion_rate = base_rate * depth_factor * temp_factor * oxygen_factor
        
        return corrosion_rate
    
    def design_protection_system(self,
                                  structure_area: float,
                                  design_life: float = 25,
                                  depth_range: Tuple[float, float] = (-50, 20)) -> Dict:
        """设计防腐蚀系统"""
        logger.info("\n设计防腐蚀系统...")
        
        # 分区面积估算
        areas = {
            'atmospheric': structure_area * 0.1,
            'splash': structure_area * 0.05,
            'tidal': structure_area * 0.05,
            'submerged': structure_area * 0.7,
            'mud': structure_area * 0.1,
        }
        
        # 无防护时的腐蚀损失
        total_corrosion_loss = 0
        for zone, area in areas.items():
            rate_range = self.zones[zone]
            avg_rate = np.mean(rate_range) / 1000  # 转换为mm/年
            loss = avg_rate * design_life
            total_corrosion_loss += loss * area
        
        # 防护方案设计
        protection_schemes = []
        
        # 方案1：仅涂层
        scheme1 = {
            'name': '涂层系统',
            'coating_type': 'epoxy',
            'coating_thickness': 0.5,  # mm
            'cp_system': False,
            'initial_cost': structure_area * self.coating_costs['epoxy'],
            'maintenance_cost': structure_area * 20 * design_life,  # 定期维护
            'effectiveness': 0.95,
        }
        scheme1['total_cost'] = scheme1['initial_cost'] + scheme1['maintenance_cost']
        scheme1['residual_corrosion'] = total_corrosion_loss * (1 - scheme1['effectiveness'])
        protection_schemes.append(scheme1)
        
        # 方案2：涂层+阴极保护
        scheme2 = {
            'name': '涂层+阴极保护',
            'coating_type': 'epoxy',
            'coating_thickness': 0.5,
            'cp_system': True,
            'initial_cost': structure_area * (self.coating_costs['epoxy'] + 
                                              self.cathodic_protection_cost),
            'maintenance_cost': structure_area * 10 * design_life,  # 维护较少
            'effectiveness': 0.99,
        }
        scheme2['total_cost'] = scheme2['initial_cost'] + scheme2['maintenance_cost']
        scheme2['residual_corrosion'] = total_corrosion_loss * (1 - scheme2['effectiveness'])
        protection_schemes.append(scheme2)
        
        # 方案3：高性能涂层
        scheme3 = {
            'name': '高性能涂层系统',
            'coating_type': 'metallic',
            'coating_thickness': 0.3,
            'cp_system': False,
            'initial_cost': structure_area * self.coating_costs['metallic'],
            'maintenance_cost': structure_area * 10 * design_life,
            'effectiveness': 0.98,
        }
        scheme3['total_cost'] = scheme3['initial_cost'] + scheme3['maintenance_cost']
        scheme3['residual_corrosion'] = total_corrosion_loss * (1 - scheme3['effectiveness'])
        protection_schemes.append(scheme3)
        
        # 选择最优方案
        best_scheme = min(protection_schemes, key=lambda x: x['total_cost'])
        
        logger.info(f"  无防护总腐蚀损失: {total_corrosion_loss/1e6:.2f} m³")
        logger.info(f"  最优方案: {best_scheme['name']}")
        logger.info(f"  总成本: ${best_scheme['total_cost']/1e6:.2f}M")
        logger.info(f"  防护效率: {best_scheme['effectiveness']:.1%}")
        
        return {
            'total_area': structure_area,
            'design_life': design_life,
            'unprotected_corrosion': total_corrosion_loss / 1e6,  # m³
            'schemes': protection_schemes,
            'optimal_scheme': best_scheme
        }


# ============== 可视化 ==============

class OceanVisualizer:
    """海洋工程可视化"""
    
    def __init__(self, output_dir: Path):
        self.output_dir = output_dir
        self.output_dir.mkdir(parents=True, exist_ok=True)
    
    def plot_platform_design(self, design: Dict, env: MarineEnvironment):
        """绘制平台设计图"""
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        # 1. 平台侧视图
        ax1 = axes[0, 0]
        
        # 海平面
        water_depth = env.water_depth
        ax1.axhline(y=0, color='blue', linestyle='-', linewidth=2, alpha=0.5, label='Sea Level')
        ax1.fill_between([-50, 50], [-water_depth, -water_depth], [0, 0], 
                         color='lightblue', alpha=0.3)
        
        # 导管架
        height = design['jacket_height']
        base_width = design['base_width']
        top_width = base_width * 0.3
        
        # 绘制主腿
        num_legs = design['num_legs']
        for i in range(num_legs):
            x_base = -base_width/2 + base_width * i / (num_legs - 1)
            x_top = -top_width/2 + top_width * i / (num_legs - 1)
            ax1.plot([x_base, x_top], [-height, 0], 'k-', linewidth=8)
        
        # 绘制支撑
        num_levels = int(height / 10)
        for level in range(num_levels):
            y = -height + (level + 1) * height / num_levels
            width_at_y = base_width - (base_width - top_width) * (height + y) / height
            ax1.plot([-width_at_y/2, width_at_y/2], [y, y], 'k-', linewidth=3)
        
        # 甲板
        deck_width = top_width * 1.5
        ax1.add_patch(Rectangle((-deck_width/2, 0), deck_width, 5, 
                                facecolor='orange', edgecolor='red', linewidth=2))
        ax1.text(0, 2.5, 'DECK', ha='center', va='center', fontsize=10, fontweight='bold')
        
        ax1.set_xlim(-30, 30)
        ax1.set_ylim(-height-10, 15)
        ax1.set_aspect('equal')
        ax1.set_xlabel('Width (m)')
        ax1.set_ylabel('Elevation (m)')
        ax1.set_title(f'Jacket Platform Side View\n{num_legs}-leg, H={height:.1f}m')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. 载荷分布
        ax2 = axes[0, 1]
        
        loads = ['Wave', 'Current', 'Wind', 'Total']
        values = [
            design['wave_force'],
            design['current_force'],
            design['wind_force'],
            design['wave_force'] + design['current_force'] + design['wind_force']
        ]
        colors = ['blue', 'cyan', 'green', 'red']
        
        bars = ax2.bar(loads, values, color=colors, edgecolor='black')
        ax2.set_ylabel('Force (kN)')
        ax2.set_title('Environmental Loads')
        ax2.grid(True, alpha=0.3, axis='y')
        
        for bar, val in zip(bars, values):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 10,
                    f'{val:.0f}', ha='center', fontsize=9)
        
        # 3. 重量分布
        ax3 = axes[1, 0]
        
        components = ['Legs', 'Braces', 'Total']
        weights = [
            design['leg_weight'],
            design['brace_weight'],
            design['total_weight']
        ]
        
        colors_pie = ['steelblue', 'lightblue', 'orange']
        explode = (0, 0, 0.1)
        
        ax3.pie(weights[:2], labels=components[:2], autopct='%1.1f%%',
               colors=colors_pie[:2], explode=(0, 0.05),
               startangle=90)
        ax3.set_title(f'Structural Weight Distribution\nTotal: {design["total_weight"]:.0f} tons')
        
        # 4. 设计参数
        ax4 = axes[1, 1]
        ax4.axis('off')
        
        params = [
            f"Platform Design Summary",
            f"{'='*30}",
            f"Number of Legs: {design['num_legs']}",
            f"Leg Diameter: {design['leg_diameter']:.0f} mm",
            f"Leg Thickness: {design['leg_thickness']:.0f} mm",
            f"Brace Diameter: {design['brace_diameter']:.0f} mm",
            f"Jacket Height: {design['jacket_height']:.1f} m",
            f"Base Width: {design['base_width']:.1f} m",
            f"",
            f"Structural Weight: {design['total_weight']:.0f} tons",
            f"  - Legs: {design['leg_weight']:.0f} tons",
            f"  - Braces: {design['brace_weight']:.0f} tons",
            f"",
            f"Max Stress: {design['max_stress']:.1f} MPa",
            f"Utilization: {design['utilization']:.1%}",
            f"Natural Period: {design['natural_period']:.2f} s",
            f"",
            f"Total Cost: ${design['total_cost']:.2f}M",
        ]
        
        y_pos = 0.95
        for param in params:
            if param.startswith('='):
                ax4.text(0.1, y_pos, param, fontsize=10, fontweight='bold',
                        transform=ax4.transAxes)
            elif 'Summary' in param:
                ax4.text(0.1, y_pos, param, fontsize=12, fontweight='bold',
                        transform=ax4.transAxes)
            else:
                ax4.text(0.1, y_pos, param, fontsize=9,
                        transform=ax4.transAxes, family='monospace')
            y_pos -= 0.05
        
        plt.tight_layout()
        plt.savefig(self.output_dir / 'platform_design.png', dpi=150, bbox_inches='tight')
        plt.close()
        logger.info("平台设计图已保存")
    
    def plot_subsea_design(self, design: Dict, env: MarineEnvironment):
        """绘制水下结构设计图"""
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
        
        # 1. 三维示意图
        ax1 = axes[0]
        
        width = design['width']
        height = design['height']
        length = design['length']
        
        # 绘制矩形表示结构
        rect = Rectangle((0, 0), width, length, facecolor='lightgray', 
                         edgecolor='black', linewidth=2)
        ax1.add_patch(rect)
        
        # 标注尺寸
        ax1.annotate('', xy=(width, -0.5), xytext=(0, -0.5),
                    arrowprops=dict(arrowstyle='<->', color='red'))
        ax1.text(width/2, -1, f'W={width:.1f}m', ha='center', color='red')
        
        ax1.annotate('', xy=(-0.5, length), xytext=(-0.5, 0),
                    arrowprops=dict(arrowstyle='<->', color='red'))
        ax1.text(-1.5, length/2, f'L={length:.1f}m', ha='center', 
                rotation=90, color='red', va='center')
        
        ax1.set_xlim(-2, width + 2)
        ax1.set_ylim(-2, length + 2)
        ax1.set_aspect('equal')
        ax1.set_title(f'Subsea Structure Plan View\nHeight: {height:.1f}m')
        ax1.set_xlabel('Width (m)')
        ax1.set_ylabel('Length (m)')
        ax1.grid(True, alpha=0.3)
        
        # 2. 稳定性分析
        ax2 = axes[1]
        
        safety_factors = {
            'Sliding': design['SF_sliding'],
            'Overturning': design['SF_overturning'],
            'Bearing': design['SF_bearing']
        }
        
        colors = ['green' if sf > 2 else 'orange' if sf > 1.5 else 'red' 
                 for sf in safety_factors.values()]
        
        bars = ax2.barh(list(safety_factors.keys()), list(safety_factors.values()), 
                       color=colors, edgecolor='black')
        ax2.axvline(x=1.5, color='red', linestyle='--', label='Min Required')
        ax2.axvline(x=2.0, color='orange', linestyle='--', label='Recommended')
        ax2.set_xlabel('Safety Factor')
        ax2.set_title('Stability Safety Factors')
        ax2.legend()
        ax2.grid(True, alpha=0.3, axis='x')
        
        for bar, (name, val) in zip(bars, safety_factors.items()):
            ax2.text(val + 0.1, bar.get_y() + bar.get_height()/2,
                    f'{val:.2f}', va='center', fontsize=9)
        
        # 3. 载荷与重量
        ax3 = axes[2]
        
        forces = {
            'Dry Weight': design['dry_weight'],
            'Buoyancy': design['submerged_weight'] + design['dry_weight'],
            'Submerged Wt': design['submerged_weight'],
            'Horiz. Load': design['horizontal_load'],
            'Lift Force': design['lift_force']
        }
        
        colors_forces = ['gray', 'blue', 'green', 'red', 'orange']
        bars = ax3.bar(range(len(forces)), list(forces.values()), 
                      color=colors_forces, edgecolor='black')
        ax3.set_xticks(range(len(forces)))
        ax3.set_xticklabels(list(forces.keys()), rotation=45, ha='right')
        ax3.set_ylabel('Force (kN)')
        ax3.set_title('Forces on Subsea Structure')
        ax3.grid(True, alpha=0.3, axis='y')
        
        for bar, val in zip(bars, forces.values()):
            ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 20,
                    f'{val:.0f}', ha='center', fontsize=8)
        
        plt.tight_layout()
        plt.savefig(self.output_dir / 'subsea_design.png', dpi=150, bbox_inches='tight')
        plt.close()
        logger.info("水下结构设计图已保存")
    
    def plot_corrosion_protection(self, protection_design: Dict):
        """绘制防腐蚀设计图"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # 1. 腐蚀速率随深度变化
        ax1 = axes[0, 0]
        
        depths = np.linspace(-100, 20, 100)
        # 腐蚀速率计算函数
        def calc_corrosion_rate(depth, temperature=15, oxygen_level=6):
            if depth > 0:
                depth_factor = 1 + 0.5 * np.exp(-depth / 20)
            else:
                depth_factor = 2.0
            temp_factor = np.exp(0.08 * (temperature - 15))
            oxygen_factor = oxygen_level / 6
            base_rate = 0.1
            return base_rate * depth_factor * temp_factor * oxygen_factor
        
        rates = [calc_corrosion_rate(d) for d in depths]
        
        ax1.fill_betweenx(depths, 0, rates, alpha=0.3, color='red')
        ax1.plot(rates, depths, 'r-', linewidth=2)
        ax1.axhline(y=0, color='blue', linestyle='--', label='Sea Level')
        ax1.set_xlabel('Corrosion Rate (mm/year)')
        ax1.set_ylabel('Depth (m)')
        ax1.set_title('Corrosion Rate vs Depth')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 添加区域标注
        ax1.text(0.5, 10, 'Atmospheric', fontsize=9, ha='center')
        ax1.text(0.5, -2, 'Splash', fontsize=9, ha='center', color='red')
        ax1.text(0.5, -30, 'Submerged', fontsize=9, ha='center', color='blue')
        
        # 2. 防护方案对比
        ax2 = axes[0, 1]
        
        schemes = protection_design['schemes']
        names = [s['name'] for s in schemes]
        costs = [s['total_cost'] / 1e6 for s in schemes]
        effectiveness = [s['effectiveness'] * 100 for s in schemes]
        
        x = np.arange(len(names))
        width = 0.35
        
        bars1 = ax2.bar(x - width/2, costs, width, label='Cost ($M)', color='orange')
        ax2_twin = ax2.twinx()
        bars2 = ax2_twin.bar(x + width/2, effectiveness, width, label='Effectiveness (%)', 
                            color='green', alpha=0.7)
        
        ax2.set_xlabel('Protection Scheme')
        ax2.set_ylabel('Cost ($M)', color='orange')
        ax2_twin.set_ylabel('Effectiveness (%)', color='green')
        ax2.set_xticks(x)
        ax2.set_xticklabels(names, rotation=15, ha='right')
        ax2.set_title('Protection Schemes Comparison')
        
        # 添加数值标签
        for bar in bars1:
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1,
                    f'${bar.get_height():.2f}M', ha='center', fontsize=8)
        
        # 3. 腐蚀损失估算
        ax3 = axes[1, 0]
        
        years = np.arange(0, protection_design['design_life'] + 1)
        
        # 无防护
        unprotected = [protection_design['unprotected_corrosion'] * y / 25 for y in years]
        
        # 有防护
        protected = [u * (1 - protection_design['optimal_scheme']['effectiveness']) 
                     for u in unprotected]
        
        ax3.fill_between(years, 0, unprotected, alpha=0.3, color='red', label='Unprotected')
        ax3.fill_between(years, 0, protected, alpha=0.5, color='green', label='Protected')
        ax3.plot(years, unprotected, 'r--', linewidth=2)
        ax3.plot(years, protected, 'g-', linewidth=2)
        
        ax3.set_xlabel('Service Life (years)')
        ax3.set_ylabel('Cumulative Corrosion Loss (m³)')
        ax3.set_title('Corrosion Loss Over Time')
        ax3.legend()
        ax3.grid(True, alpha=0.3)
        
        # 4. 防护系统详情
        ax4 = axes[1, 1]
        ax4.axis('off')
        
        optimal = protection_design['optimal_scheme']
        
        info = [
            f"Corrosion Protection Design",
            f"{'='*35}",
            f"Total Area: {protection_design['total_area']:.0f} m²",
            f"Design Life: {protection_design['design_life']} years",
            f"",
            f"Unprotected Corrosion:",
            f"  {protection_design['unprotected_corrosion']:.2f} m³",
            f"",
            f"Optimal Scheme: {optimal['name']}",
            f"Coating Type: {optimal['coating_type']}",
            f"Coating Thickness: {optimal['coating_thickness']} mm",
            f"CP System: {'Yes' if optimal['cp_system'] else 'No'}",
            f"",
            f"Protection Effectiveness:",
            f"  {optimal['effectiveness']:.1%}",
            f"",
            f"Cost Breakdown:",
            f"  Initial: ${optimal['initial_cost']/1e6:.2f}M",
            f"  Maintenance: ${optimal['maintenance_cost']/1e6:.2f}M",
            f"  Total: ${optimal['total_cost']/1e6:.2f}M",
            f"",
            f"Residual Corrosion:",
            f"  {optimal['residual_corrosion']:.3f} m³",
        ]
        
        y_pos = 0.95
        for line in info:
            if line.startswith('='):
                ax4.text(0.1, y_pos, line, fontsize=10, fontweight='bold',
                        transform=ax4.transAxes)
            elif 'Corrosion Protection' in line:
                ax4.text(0.1, y_pos, line, fontsize=12, fontweight='bold',
                        transform=ax4.transAxes)
            else:
                ax4.text(0.1, y_pos, line, fontsize=9,
                        transform=ax4.transAxes, family='monospace')
            y_pos -= 0.05
        
        plt.tight_layout()
        plt.savefig(self.output_dir / 'corrosion_protection.png', dpi=150, bbox_inches='tight')
        plt.close()
        logger.info("防腐蚀设计图已保存")
    
    def plot_environmental_conditions(self, env: MarineEnvironment):
        """绘制海洋环境条件"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # 1. 波浪谱
        ax1 = axes[0, 0]
        
        omega = np.linspace(0.1, 3, 100)
        
        # JONSWAP谱（简化）
        omega_p = 2 * np.pi / env.wave_period
        alpha = 0.0081
        gamma = 3.3
        
        S = alpha * env.gravity**2 / omega**5 * \
            np.exp(-1.25 * (omega_p / omega)**4) * \
            gamma ** np.exp(-0.5 * ((omega - omega_p) / (0.07 * omega_p))**2)
        
        ax1.fill_between(omega, 0, S, alpha=0.3, color='blue')
        ax1.plot(omega, S, 'b-', linewidth=2)
        ax1.axvline(x=omega_p, color='red', linestyle='--', label=f'Peak ω={omega_p:.2f} rad/s')
        ax1.set_xlabel('Frequency (rad/s)')
        ax1.set_ylabel('Spectral Density (m²·s)')
        ax1.set_title('Wave Spectrum (JONSWAP)')
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. 流速剖面
        ax2 = axes[0, 1]
        
        depths = np.linspace(0, env.water_depth, 50)
        velocities = env.current_velocity * (depths / env.water_depth) ** 0.2
        
        ax2.fill_betweenx(depths, 0, velocities, alpha=0.3, color='cyan')
        ax2.plot(velocities, depths, 'c-', linewidth=2)
        ax2.set_xlabel('Current Velocity (m/s)')
        ax2.set_ylabel('Depth (m)')
        ax2.set_title('Current Velocity Profile')
        ax2.invert_yaxis()
        ax2.grid(True, alpha=0.3)
        
        # 3. 环境参数汇总
        ax3 = axes[1, 0]
        ax3.axis('off')
        
        params = [
            f"Marine Environment Parameters",
            f"{'='*35}",
            f"Water Depth: {env.water_depth:.0f} m",
            f"Significant Wave Height: {env.wave_height:.1f} m",
            f"Peak Wave Period: {env.wave_period:.1f} s",
            f"Surface Current: {env.current_velocity:.2f} m/s",
            f"Wind Speed: {env.wind_speed:.1f} m/s",
            f"",
            f"Derived Parameters:",
            f"Peak Frequency: {omega_p:.3f} rad/s",
            f"Wave Length: {env.gravity * env.wave_period**2 / (2 * np.pi):.1f} m",
            f"",
            f"Water Properties:",
            f"Density: {env.seawater_density} kg/m³",
            f"Gravity: {env.gravity} m/s²",
        ]
        
        y_pos = 0.95
        for line in params:
            if line.startswith('='):
                ax3.text(0.1, y_pos, line, fontsize=10, fontweight='bold',
                        transform=ax3.transAxes)
            elif 'Marine Environment' in line:
                ax3.text(0.1, y_pos, line, fontsize=12, fontweight='bold',
                        transform=ax3.transAxes)
            else:
                ax3.text(0.1, y_pos, line, fontsize=10,
                        transform=ax3.transAxes, family='monospace')
            y_pos -= 0.06
        
        # 4. 设计海况
        ax4 = axes[1, 1]
        
        categories = ['Normal', 'Extreme', 'Survival']
        wave_heights = [env.wave_height, env.wave_height * 1.5, env.wave_height * 2.0]
        wind_speeds = [env.wind_speed, env.wind_speed * 1.3, env.wind_speed * 1.6]
        
        x = np.arange(len(categories))
        width = 0.35
        
        ax4_twin = ax4.twinx()
        bars1 = ax4.bar(x - width/2, wave_heights, width, label='Wave Height', color='blue')
        bars2 = ax4_twin.bar(x + width/2, wind_speeds, width, label='Wind Speed', color='green')
        
        ax4.set_xlabel('Design Condition')
        ax4.set_ylabel('Wave Height (m)', color='blue')
        ax4_twin.set_ylabel('Wind Speed (m/s)', color='green')
        ax4.set_xticks(x)
        ax4.set_xticklabels(categories)
        ax4.set_title('Design Environmental Conditions')
        
        for bar in bars1:
            ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1,
                    f'{bar.get_height():.1f}', ha='center', fontsize=9)
        
        plt.tight_layout()
        plt.savefig(self.output_dir / 'environmental_conditions.png', dpi=150, bbox_inches='tight')
        plt.close()
        logger.info("环境条件图已保存")


# ============== 主程序 ==============

def main():
    """主程序"""
    logger.info("="*60)
    logger.info("主题096：海洋工程结构优化")
    logger.info("Ocean Engineering Structure Optimization")
    logger.info("="*60)
    
    # 设置输出目录
    output_dir = Path(__file__).parent / 'output'
    visualizer = OceanVisualizer(output_dir)
    
    # 定义海洋环境（北海典型条件）
    env = MarineEnvironment(
        water_depth=80,      # m
        wave_height=12,      # m (百年一遇)
        wave_period=14,      # s
        current_velocity=1.2, # m/s
        wind_speed=35        # m/s
    )
    
    logger.info(f"\n海洋环境参数:")
    logger.info(f"  水深: {env.water_depth} m")
    logger.info(f"  有义波高: {env.wave_height} m")
    logger.info(f"  谱峰周期: {env.wave_period} s")
    logger.info(f"  表层流速: {env.current_velocity} m/s")
    logger.info(f"  风速: {env.wind_speed} m/s")
    
    # 1. 海洋平台优化
    logger.info("\n" + "="*60)
    logger.info("1. 海洋平台结构优化")
    logger.info("="*60)
    
    platform_optimizer = OffshorePlatformOptimizer(env)
    platform_result = platform_optimizer.optimize()
    
    if platform_result['optimal_design']:
        best_platform = platform_result['optimal_design']
        logger.info(f"\n最优平台设计:")
        logger.info(f"  主腿数量: {best_platform['num_legs']}")
        logger.info(f"  导管架高度: {best_platform['jacket_height']:.1f} m")
        logger.info(f"  底部宽度: {best_platform['base_width']:.1f} m")
        logger.info(f"  总重量: {best_platform['total_weight']:.0f} 吨")
        logger.info(f"  总成本: ${best_platform['total_cost']:.2f}M")
        logger.info(f"  应力利用率: {best_platform['utilization']:.1%}")
        logger.info(f"  自振周期: {best_platform['natural_period']:.2f} s")
        
        visualizer.plot_platform_design(best_platform, env)
    
    # 2. 水下结构优化
    logger.info("\n" + "="*60)
    logger.info("2. 水下结构优化")
    logger.info("="*60)
    
    subsea_optimizer = SubseaStructureOptimizer(env)
    subsea_design = subsea_optimizer.optimize_manifold()
    
    logger.info(f"\n水下管汇设计:")
    logger.info(f"  尺寸: {subsea_design['width']:.1f}m × {subsea_design['height']:.1f}m × {subsea_design['length']:.1f}m")
    logger.info(f"  混凝土体积: {subsea_design['concrete_volume']:.1f} m³")
    logger.info(f"  干重: {subsea_design['dry_weight']:.0f} kN")
    logger.info(f"  抗滑移安全系数: {subsea_design['SF_sliding']:.2f}")
    logger.info(f"  抗倾覆安全系数: {subsea_design['SF_overturning']:.2f}")
    
    visualizer.plot_subsea_design(subsea_design, env)
    
    # 3. 防腐蚀设计
    logger.info("\n" + "="*60)
    logger.info("3. 防腐蚀系统设计")
    logger.info("="*60)
    
    corrosion_designer = CorrosionProtectionDesign(env)
    
    # 假设平台湿表面积
    platform_area = 5000  # m²
    protection_design = corrosion_designer.design_protection_system(
        platform_area, design_life=25
    )
    
    visualizer.plot_corrosion_protection(protection_design)
    
    # 4. 环境条件可视化
    logger.info("\n" + "="*60)
    logger.info("4. 环境条件可视化")
    logger.info("="*60)
    
    visualizer.plot_environmental_conditions(env)
    
    logger.info("\n" + "="*60)
    logger.info("所有计算和可视化已完成！")
    logger.info(f"输出目录: {output_dir}")
    logger.info("="*60)


if __name__ == "__main__":
    main()