主题049: 海洋工程流动

1. 理论基础

1.1 波浪理论

线性波浪理论(Airy波):
$$\eta =a\cos (kx-\omega t)$$

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

深水近似 ($kh\gg 1$):
$${\omega }^2=gk,c=\sqrt{\frac{g}{k}}=\frac{gT}{2\pi }$$

浅水近似 ($kh\ll 1$):
$${\omega }^2=gh{k}^2,c=\sqrt{gh}$$

1.2 莫里森方程

作用在柱体上的波浪力：
$$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}$$

其中：

• $C_D$: 阻力系数
• $C_M$: 惯性力系数 ($C_M=C_a+1$)
• $C_a$: 附加质量系数

1.3 弗劳德数相似

$$Fr=\frac{U}{\sqrt{gL}}$$

模型试验需保持弗劳德数相等：
$${\lambda }_U=\sqrt{{\lambda }_L},{\lambda }_T=\sqrt{{\lambda }_L}$$

2. 数值方法

2.1 波浪模拟

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

class OceanEngineeringFlow:
    """
    海洋工程流动模拟
    """
    
    def __init__(self, g=9.81, rho=1025):
        """
        参数:
            g: 重力加速度
            rho: 海水密度
        """
        self.g = g
        self.rho = rho
    
    def airy_wave(self, x, z, t, a, T, h):
        """
        线性Airy波
        
        参数:
            x: 水平位置
            z: 垂直位置 (向上为正，自由面为0)
            t: 时间
            a: 波幅
            T: 周期
            h: 水深
        """
        omega = 2 * np.pi / T
        
        # 求解波数 (色散关系)
        k = self.wavenumber(omega, h)
        
        # 自由面升高
        eta = a * np.cos(k * x - omega * t)
        
        # 速度势
        phi = self.g * a / omega * np.cosh(k * (z + h)) / np.cosh(k * h) * np.sin(k * x - omega * t)
        
        # 速度
        u = self.g * a * k / omega * np.cosh(k * (z + h)) / np.cosh(k * h) * np.cos(k * x - omega * t)
        w = self.g * a * k / omega * np.sinh(k * (z + h)) / np.cosh(k * h) * np.sin(k * x - omega * t)
        
        # 压力 (线性化)
        p = -self.rho * self.g * z + self.rho * self.g * a * np.cosh(k * (z + h)) / np.cosh(k * h) * np.cos(k * x - omega * t)
        
        return {
            'eta': eta, 'phi': phi, 'u': u, 'w': w, 'p': p,
            'k': k, 'omega': omega, 'L': 2 * np.pi / k
        }
    
    def wavenumber(self, omega, h):
        """
        求解波数 (迭代法)
        """
        # 初始猜测 (深水)
        k = omega**2 / self.g
        
        # 迭代求解
        for _ in range(100):
            k_new = omega**2 / (self.g * np.tanh(k * h))
            if abs(k_new - k) < 1e-10:
                break
            k = k_new
        
        return k
    
    def morison_force(self, u, du_dt, D, C_D=1.0, C_M=2.0):
        """
        莫里森方程计算波浪力
        
        参数:
            u: 水平速度
            du_dt: 水平加速度
            D: 柱体直径
            C_D: 阻力系数
            C_M: 惯性力系数
        """
        # 阻力
        F_D = 0.5 * self.rho * C_D * D * u * np.abs(u)
        
        # 惯性力
        F_I = self.rho * C_M * np.pi * D**2 / 4 * du_dt
        
        # 总力
        F = F_D + F_I
        
        return F, F_D, F_I
    
    def wave_spectrum_jonswap(self, omega, Hs, Tp, gamma=3.3):
        """
        JONSWAP波谱
        
        参数:
            omega: 角频率数组
            Hs: 有效波高
            Tp: 峰值周期
            gamma: 峰升因子
        """
        omega_p = 2 * np.pi / Tp  # 峰值频率
        
        # 阿尔法参数
        alpha = 0.076 * (Hs**2 * omega_p**4 / self.g**2)**0.22
        
        # 西格玛参数
        sigma = np.where(omega <= omega_p, 0.07, 0.09)
        
        # 峰升因子
        A = np.exp(-((omega - omega_p)**2) / (2 * sigma**2 * omega_p**2))
        
        # PM谱
        S_pm = alpha * self.g**2 / omega**5 * np.exp(-1.25 * (omega_p / omega)**4)
        
        # JONSWAP谱
        S = S_pm * gamma**A
        
        return S
    
    def irregular_wave(self, x, z, t, Hs, Tp, n_components=50):
        """
        不规则波模拟 (线性叠加)
        """
        # 频率范围
        omega_min = 0.5 * 2 * np.pi / Tp
        omega_max = 3.0 * 2 * np.pi / Tp
        omega = np.linspace(omega_min, omega_max, n_components)
        domega = omega[1] - omega[0]
        
        # 波谱
        S = self.wave_spectrum_jonswap(omega, Hs, Tp)
        
        # 随机相位
        np.random.seed(42)
        phi = np.random.uniform(0, 2*np.pi, n_components)
        
        # 振幅
        a = np.sqrt(2 * S * domega)
        
        # 叠加
        eta = 0
        u = 0
        w = 0
        
        for i in range(n_components):
            k = self.wavenumber(omega[i], 50)  # 假设水深50m
            
            eta += a[i] * np.cos(k * x - omega[i] * t + phi[i])
            u += a[i] * omega[i] * np.cosh(k * (z + 50)) / np.sinh(k * 50) * np.cos(k * x - omega[i] * t + phi[i])
            w += a[i] * omega[i] * np.sinh(k * (z + 50)) / np.sinh(k * 50) * np.sin(k * x - omega[i] * t + phi[i])
        
        return eta, u, w
    
    def spar_platform_dynamics(self, y, t, M, C, K, F_wave):
        """
        Spar平台运动方程
        """
        # y = [位移, 速度]
        # dy/dt = [速度, 加速度]
        
        displacement = y[0]
        velocity = y[1]
        
        # 波浪力 (简化)
        F = F_wave(t)
        
        # 恢复力
        F_restoring = -K * displacement
        
        # 阻尼力
        F_damping = -C * velocity
        
        # 加速度
        acceleration = (F + F_restoring + F_damping) / M
        
        return [velocity, acceleration]


def simulate_regular_wave():
    """
    规则波模拟
    """
    ocean = OceanEngineeringFlow()
    
    # 波浪参数
    H = 4.0  # 波高
    a = H / 2  # 波幅
    T = 10.0  # 周期
    h = 50.0  # 水深
    
    # 计算波长
    omega = 2 * np.pi / T
    k = ocean.wavenumber(omega, h)
    L = 2 * np.pi / k
    
    print("=" * 60)
    print("规则波参数")
    print("=" * 60)
    print(f"波高 H = {H} m")
    print(f"周期 T = {T} s")
    print(f"水深 h = {h} m")
    print(f"波长 L = {L:.2f} m")
    print(f"波数 k = {k:.4f} rad/m")
    print(f"波速 c = {L/T:.2f} m/s")
    
    # 可视化
    x = np.linspace(0, 2*L, 200)
    z = np.linspace(-h, 0, 50)
    X, Z = np.meshgrid(x, z)
    
    t = 0
    wave = ocean.airy_wave(X, Z, t, a, T, h)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 自由面形状
    ax1 = axes[0, 0]
    eta = a * np.cos(k * x - omega * t)
    ax1.plot(x, eta, 'b-', linewidth=2)
    ax1.fill_between(x, eta, -h, alpha=0.3, color='blue')
    ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
    ax1.set_xlabel('x (m)', fontsize=12)
    ax1.set_ylabel('η (m)', fontsize=12)
    ax1.set_title('Free Surface Elevation', fontsize=12)
    ax1.set_ylim([-h, H])
    ax1.grid(True, alpha=0.3)
    
    # 2. 速度场
    ax2 = axes[0, 1]
    speed = np.sqrt(wave['u']**2 + wave['w']**2)
    contour = ax2.contourf(X, Z, speed, levels=20, cmap='viridis')
    plt.colorbar(contour, ax=ax2, label='Velocity magnitude (m/s)')
    ax2.quiver(X[::5, ::10], Z[::5, ::10], 
               wave['u'][::5, ::10], wave['w'][::5, ::10],
               scale=5, alpha=0.7)
    ax2.set_xlabel('x (m)', fontsize=12)
    ax2.set_ylabel('z (m)', fontsize=12)
    ax2.set_title('Velocity Field', fontsize=12)
    
    # 3. 压力分布
    ax3 = axes[1, 0]
    p = wave['p'] / 1000  # kPa
    contour = ax3.contourf(X, Z, p, levels=20, cmap='RdBu_r')
    plt.colorbar(contour, ax=ax3, label='Pressure (kPa)')
    ax3.set_xlabel('x (m)', fontsize=12)
    ax3.set_ylabel('z (m)', fontsize=12)
    ax3.set_title('Pressure Distribution', fontsize=12)
    
    # 4. 水质点轨迹
    ax4 = axes[1, 1]
    x0 = L / 4
    z0 = -h / 4
    
    T_period = np.linspace(0, T, 100)
    trajectory_x = []
    trajectory_z = []
    
    for t in T_period:
        wave_t = ocean.airy_wave(x0, z0, t, a, T, h)
        # 近似轨迹
        dx = a * np.cosh(k * (z0 + h)) / np.sinh(k * h) * np.sin(k * x0 - omega * t)
        dz = -a * np.sinh(k * (z0 + h)) / np.sinh(k * h) * np.cos(k * x0 - omega * t)
        trajectory_x.append(x0 + dx)
        trajectory_z.append(z0 + dz)
    
    ax4.plot(trajectory_x, trajectory_z, 'b-', linewidth=2)
    ax4.plot(x0, z0, 'ro', markersize=10)
    ax4.set_xlabel('x (m)', fontsize=12)
    ax4.set_ylabel('z (m)', fontsize=12)
    ax4.set_title('Water Particle Trajectory', fontsize=12)
    ax4.set_aspect('equal')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('regular_wave.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n图形已保存为 'regular_wave.png'")


def simulate_morison_force():
    """
    莫里森力计算
    """
    ocean = OceanEngineeringFlow()
    
    # 波浪参数
    H = 4.0
    a = H / 2
    T = 10.0
    h = 50.0
    
    # 柱体参数
    D = 2.0  # 直径
    z = -10.0  # 计算位置 (水面下10m)
    
    # 时间序列
    t = np.linspace(0, 3*T, 300)
    
    # 波浪参数
    omega = 2 * np.pi / T
    k = ocean.wavenumber(omega, h)
    
    # 计算波浪力和速度
    forces = []
    forces_D = []
    forces_I = []
    velocities = []
    accelerations = []
    
    for ti in t:
        wave = ocean.airy_wave(0, z, ti, a, T, h)
        u = wave['u']
        du_dt = a * omega**2 * np.cosh(k * (z + h)) / np.sinh(k * h) * np.sin(-omega * ti)
        
        F, F_D, F_I = ocean.morison_force(u, du_dt, D, C_D=1.2, C_M=2.0)
        
        forces.append(F)
        forces_D.append(F_D)
        forces_I.append(F_I)
        velocities.append(u)
        accelerations.append(du_dt)
    
    forces = np.array(forces)
    forces_D = np.array(forces_D)
    forces_I = np.array(forces_I)
    velocities = np.array(velocities)
    accelerations = np.array(accelerations)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 波浪力时程
    ax1 = axes[0, 0]
    ax1.plot(t, forces/1000, 'b-', linewidth=2, label='Total')
    ax1.plot(t, forces_D/1000, 'r--', linewidth=1.5, label='Drag')
    ax1.plot(t, forces_I/1000, 'g:', linewidth=1.5, label='Inertia')
    ax1.set_xlabel('t (s)', fontsize=12)
    ax1.set_ylabel('F (kN/m)', fontsize=12)
    ax1.set_title('Wave Force Time Series', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 速度时程
    ax2 = axes[0, 1]
    ax2.plot(t, velocities, 'b-', linewidth=2)
    ax2.set_xlabel('t (s)', fontsize=12)
    ax2.set_ylabel('u (m/s)', fontsize=12)
    ax2.set_title('Wave Velocity', fontsize=12)
    ax2.grid(True, alpha=0.3)
    
    # 3. 力-速度关系
    ax3 = axes[1, 0]
    ax3.plot(velocities, forces/1000, 'b.', markersize=3)
    ax3.set_xlabel('u (m/s)', fontsize=12)
    ax3.set_ylabel('F (kN/m)', fontsize=12)
    ax3.set_title('Force-Velocity Relationship', fontsize=12)
    ax3.grid(True, alpha=0.3)
    
    # 4. Keulegan-Carpenter数
    ax4 = axes[1, 1]
    KC = []
    F_max_drag = []
    F_max_inertia = []
    
    D_range = np.linspace(0.5, 5, 20)
    for Di in D_range:
        u_max = a * omega * np.cosh(k * (z + h)) / np.sinh(k * h)
        KC_i = u_max * T / Di
        
        # 最大力
        F_D_max = 0.5 * ocean.rho * 1.2 * Di * u_max**2
        F_I_max = ocean.rho * 2.0 * np.pi * Di**2 / 4 * a * omega**2 * np.cosh(k * (z + h)) / np.sinh(k * h)
        
        KC.append(KC_i)
        F_max_drag.append(F_D_max / 1000)
        F_max_inertia.append(F_I_max / 1000)
    
    ax4.semilogy(KC, F_max_drag, 'r-', linewidth=2, label='Drag')
    ax4.semilogy(KC, F_max_inertia, 'b-', linewidth=2, label='Inertia')
    ax4.set_xlabel('KC number', fontsize=12)
    ax4.set_ylabel('Max Force (kN/m)', fontsize=12)
    ax4.set_title('Drag vs Inertia Dominance', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('morison_force.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n图形已保存为 'morison_force.png'")


def simulate_irregular_wave():
    """
    不规则波模拟
    """
    ocean = OceanEngineeringFlow()
    
    # 海况参数
    Hs = 5.0  # 有效波高
    Tp = 12.0  # 峰值周期
    
    print("=" * 60)
    print("不规则波参数")
    print("=" * 60)
    print(f"有效波高 Hs = {Hs} m")
    print(f"峰值周期 Tp = {Tp} s")
    
    # 频率范围
    omega = np.linspace(0.3, 2.0, 200)
    
    # JONSWAP谱
    S = ocean.wave_spectrum_jonswap(omega, Hs, Tp)
    
    # 时间序列
    t = np.linspace(0, 300, 3000)
    x = 0
    z = -5
    
    eta, u, w = ocean.irregular_wave(x, z, t, Hs, Tp, n_components=100)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 波谱
    ax1 = axes[0, 0]
    ax1.plot(omega, S, 'b-', linewidth=2)
    ax1.fill_between(omega, S, alpha=0.3)
    ax1.set_xlabel('ω (rad/s)', fontsize=12)
    ax1.set_ylabel('S(ω) (m²·s)', fontsize=12)
    ax1.set_title('JONSWAP Wave Spectrum', fontsize=12)
    ax1.grid(True, alpha=0.3)
    
    # 2. 自由面时程
    ax2 = axes[0, 1]
    ax2.plot(t, eta, 'b-', linewidth=1)
    ax2.axhline(y=Hs/2, color='r', linestyle='--', label='H/2')
    ax2.axhline(y=-Hs/2, color='r', linestyle='--')
    ax2.set_xlabel('t (s)', fontsize=12)
    ax2.set_ylabel('η (m)', fontsize=12)
    ax2.set_title('Irregular Wave Elevation', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim([0, 300])
    
    # 3. 速度时程
    ax3 = axes[1, 0]
    ax3.plot(t, u, 'b-', linewidth=1, label='u')
    ax3.plot(t, w, 'r-', linewidth=1, label='w')
    ax3.set_xlabel('t (s)', fontsize=12)
    ax3.set_ylabel('Velocity (m/s)', fontsize=12)
    ax3.set_title('Wave Velocity Time Series', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim([0, 300])
    
    # 4. 统计分布
    ax4 = axes[1, 1]
    ax4.hist(eta, bins=50, density=True, alpha=0.7, color='blue', label='Data')
    
    # 理论Rayleigh分布
    from scipy.stats import rayleigh
    x_rayleigh = np.linspace(0, max(eta), 100)
    sigma = Hs / (2 * np.sqrt(2))
    ax4.plot(x_rayleigh, rayleigh.pdf(x_rayleigh, scale=sigma), 'r-', linewidth=2, label='Rayleigh')
    
    ax4.set_xlabel('η (m)', fontsize=12)
    ax4.set_ylabel('PDF', fontsize=12)
    ax4.set_title('Wave Height Distribution', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('irregular_wave.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n图形已保存为 'irregular_wave.png'")


if __name__ == "__main__":
    simulate_regular_wave()
    simulate_morison_force()
    simulate_irregular_wave()