这是我们课程的大作业，题目如下：

设某海域水深30米，其波面适用Stokes五阶波理论。

1）试用Morison方程计算1米直径圆立柱的波浪荷载，包括水平波浪力和波浪力矩；

2）讨论波数、波高等波浪参数对波浪荷载峰值的影响。

这里需要用到Stokes五阶波理论与Morison方程，具体理论这里不做过多介绍，可参考书籍与文献。

我采用了python编程来解决这个问题，并且给出我的做题过程，其中存在的错误，欢迎与讨论交流：

1）参数基本假设：海水密度：ρs=1.05g/cm3,惯性系数：Cm=2.0,拖曳系数Cd=1.0,立柱直径：D=1m,水深:h=30m,假设波高H=5m,周期T=10s。

Stokes5阶波:

                      （1）

                      （2）

(1)、(2)超越方程求解困难，这里使用了经验公式（3）[该经验公式参考：李炎保,杨鹏.五阶Stokes波计算的几个注意点[J].海洋科学,1997,(02):65-67.]：

   （3）

求得L=140.93，波数k=2π/L=0.0445,接着对式（1）进行变换,如式（4）所示，进行迭代求解，求出λ=0.11

                      （4）

根据式（5）计算波面抬高，如图（1）所示：

                     （5）

根据式（6）计算瞬时波浪速度，如图2所示：

(6)

根据式（7）计算瞬时波浪加速度，如图3所示：

(7)

根据Morison方程，水平波浪力如式（8），波浪力矩如式（9）：

                       (8)

                  (9)

直接求解十分困难，这里采用了数值求解方法，详见代码的calculate_force_and_moment函数。根据计算，最大总波浪力:25.62 kN，最大拖曳力:10.21 kN(39.84%)，最大惯性力:24.40 kN(60.16%)，水平波浪力如图4所示：

根据计算，最大总波浪力矩:390.28 kN·m，最大拖曳力矩:161.47 kN·m (41.37%)，最大惯性力矩:375.92 kN·m (58.63%)，水平波浪力矩如图5所示：

2）这里保持周期为T=10对比了不同波高对峰值力与峰值力矩的影响，如图6所示。结果表明，波高越高峰值载荷越大

这里保持波高H=5对比了不同周期对峰值力与峰值力矩的影响，如图7所示。结果表明，周期对高峰值载荷的影响并不是简单的线性关系。

python代码如下：

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

class WaveForceCalculator:
    def __init__(self, depth=30.0, diameter=1.0, rho=1025.0, g=9.81, cm=2.0, cd=1.0):
        """初始化波浪荷载计算器"""
        self.depth = depth  # 水深(m)
        self.diameter = diameter  # 立柱直径(m)
        self.radius = diameter/2  # 圆柱半径(m)
        self.rho = rho  # 海水密度(kg/m³)
        self.g = g  # 重力加速度(m/s²)
        self.cm = cm  # 惯性力系数
        self.cd = cd  # 拖曳力系数
        
        # 横截面面积和周长
        self.area = np.pi * self.radius**2
        self.perimeter = np.pi * diameter

    def calculate_wave_parameters(self, T, d, H):
        """计算波浪参数"""
        # 计算波长L（经验公式）
        L0 = self.g * T**2 / (2 * np.pi)
        L = d / (0.0432 + (0.9599 + 0.1284*H/d - 0.3448*(H/d)**2) * d/L0 + 
                (-0.2918 - 1.1732*H/d + 0.9292*(H/d)**2) * (d/L0)**2)
        
        # 计算波数k
        k = 2 * np.pi / L
        
        # 计算双曲余弦和正弦
        c = np.cosh(k * d)
        s = np.sinh(k * d)
        
        return L, k, c, s

    def calculate_coefficients(self, c, s):
        """计算各种系数"""
        # A系数
        A11 = 1 / s
        A13 = -c**2 * (5 * c**2 + 1) / (8 * s**5)
        A15 = -(1184 * c**10 - 1440 * c**8 - 1992 * c**6 + 2641 * c**4 - 249 * c**2 + 18) / (1536 * s**11)
        A22 = 3 / (8 * s**4)
        A24 = (192 * c**8 - 424 * c**6 - 312 * c**4 + 480 * c**2 - 17) / (768 * s**10)
        A33 = (13 - 4 * c**2) / (64 * s**7)
        A35 = (512 * c**12 + 4224 * c**10 - 6800 * c**8 - 12808 * c**6 + 16704 * c**4 - 3154 * c**2 + 107) / (4096 * s**13) / (6 * c**2 - 1)
        A44 = (80 * c**6 - 816 * c**4 + 1338 * c**2 - 197) / (1536 * s**10) / (6 * c**2 - 1)
        A55 = -(2880 * c**10 - 72480 * c**8 + 324000 * c**6 - 432000 * c**4 + 163470 * c**2 - 16245) / (61440 * s**11) / (6 * c**2 - 1) / (8 * c**4 - 11 * c**2 + 3)

        # B系数
        B22 = (2 * c**2 + 1) * c / (4 * s**3)
        B24 = c * (272 * c**8 - 504 * c**6 - 192 * c**4 + 322 * c**2 + 21) / (384 * s**9)
        B44 = c * (768 * c**10 - 488 * c**8 - 48 * c**6 + 48 * c**4 + 106 * c**2 - 21) / (384 * s**9) / (6 * c**2 - 1)
        B33 = 3 * (8 * c**6 + 1) / (64 * s**6)
        B35 = (88128 * c**14 - 208224 * c**12 + 70848 * c**10 + 54000 * c**8 - 21816 * c**6 + 6264 * c**4 - 54 * c**2 - 81) / (12288 * s**12) / (6 * c**2 - 1)
        B55 = (192000 * c**16 - 262720 * c**14 + 83680 * c**12 + 20160 * c**10 - 7280 * c**8 + 7160 * c**6 - 1800 * c**4 - 1050 * c**2 + 225) / (12288 * s**10) / (6 * c**2 - 1) / (8 * c**4 - 11 * c**2 + 3)

        return {
            'A': [A11, A13, A15, A22, A24, A33, A35, A44, A55],
            'B': [B22, B24, B33, B35, B44, B55]
        }

    def get_lamda(self, k, d, H, L, B33, B35, B55):
        """计算lamda值"""
        lamda = np.pi * H / (L * d)  # 初始估计
        for _ in range(10):  # 迭代10次
            lamda_new = np.pi * H / (L * (1 + lamda**2 * B33 + lamda**4 * (B35 + B55)))
            if abs(lamda_new - lamda) < 1e-6:
                break
            lamda = lamda_new
        return lamda

    def calculate_wave_kinematics(self, t, z, T, H, d):
        """计算波浪运动学参数"""
        # 计算基本参数
        L, k, c, s = self.calculate_wave_parameters(T, d, H)
        coefs = self.calculate_coefficients(c, s)
        A11, A13, A15, A22, A24, A33, A35, A44, A55 = coefs['A']
        B22, B24, B33, B35, B44, B55 = coefs['B']
        
        # 计算lamda
        lamda = self.get_lamda(k, d, H, L, B33, B35, B55)
        
        # 计算波速
        C0s = L/T
        omega = 2 * np.pi / T
        theta = (k * z - omega * t)

        # 计算水质点水平速度
        u = C0s * ((lamda*A11+lamda**3*A13+lamda**5*A15)*np.cosh(s*k)*np.cos(theta)
                   +2*(lamda**2*A22+lamda**4*A24)*np.cosh(2*s*k)*np.cos(2*theta)
                   +3*(lamda**3*A33+lamda**5*A35)*np.cosh(3*s*k)*np.cos(3*theta)
                   +4*(lamda**4*A44)*np.cosh(4*s*k)*np.cos(4*theta)
                   +5*(lamda**5*A55)*np.cosh(5*s*k)*np.cos(5*theta))

        # 计算水质点水平加速度
        a = C0s*(omega*(lamda*A11+lamda**3*A13+lamda**5*A15)*np.cosh(s*k)*np.sin(theta)
                 +4*omega*(lamda**2*A22+lamda**4*A24)*np.cosh(2*s*k)*np.sin(2*theta)
                 +9*omega*(lamda**3*A33+lamda**5*A35)*np.cosh(3*s*k)*np.sin(3*theta)
                 +16*omega*(lamda**4*A44)*np.cosh(4*s*k)*np.sin(4*theta)
                 +25*omega*(lamda**5*A55)*np.cosh(5*s*k)*np.sin(5*theta))

        return u, a, L, k

    def drag_force_integrand(self, z, t, wave_period, wave_height, d=None):
        """拖曳力积分项 |u|u 或 (d+z)|u|u"""
        u, _, _, _ = self.calculate_wave_kinematics(t, z, wave_period, wave_height, self.depth)
        
        if d is None:
            return abs(u) * u
        else:
            return (d + z) * abs(u) * u
    
    def inertia_force_integrand(self, z, t, wave_period, wave_height, d=None):
        """惯性力积分项 a 或 (d+z)a"""
        _, a, _, _ = self.calculate_wave_kinematics(t, z, wave_period, wave_height, self.depth)
        
        if d is None:
            return a
        else:
            return (d + z) * a

    def calculate_force_and_moment(self, t, wave_period, wave_height, d=None):
        """计算总波浪力和波浪力矩
        F = Cd·ρ·r·∫|u|u dz + Cm·ρ·π·r²·∫a dz
        M = Cd·ρ·r·∫(d+z)|u|u dz + Cm·ρ·π·r²·∫(d+z)a dz
        """
        # 计算拖曳力积分
        drag_integral, _ = quad(
            self.drag_force_integrand, 
            -self.depth, 0, 
            args=(t, wave_period, wave_height, d)
        )
        
        # 计算惯性力积分
        inertia_integral, _ = quad(
            self.inertia_force_integrand, 
            -self.depth, 0, 
            args=(t, wave_period, wave_height, d)
        )
        
        # 计算总波浪力和力矩
        drag_component = self.cd * self.rho * self.radius * drag_integral
        inertia_component = self.cm * self.rho * self.area * inertia_integral
        
        return drag_component, inertia_component

    def calculate_max_values(self, wave_period, wave_height, time_points=100, d=None):
        """计算波浪力和力矩的最大值"""
        t_values = np.linspace(0, wave_period, time_points)
        
        forces = np.zeros_like(t_values)
        moments = np.zeros_like(t_values)
        drag_forces = np.zeros_like(t_values)
        inertia_forces = np.zeros_like(t_values)
        drag_moments = np.zeros_like(t_values)
        inertia_moments = np.zeros_like(t_values)
        
        for i, t in enumerate(t_values):
            # 计算波浪力
            drag_force, inertia_force = self.calculate_force_and_moment(t, wave_period, wave_height)
            total_force = drag_force + inertia_force
            
            # 计算波浪力矩
            drag_moment, inertia_moment = self.calculate_force_and_moment(t, wave_period, wave_height, self.depth)
            total_moment = drag_moment + inertia_moment
            
            # 存储结果
            forces[i] = total_force
            moments[i] = total_moment
            drag_forces[i] = drag_force
            inertia_forces[i] = inertia_force
            drag_moments[i] = drag_moment
            inertia_moments[i] = inertia_moment
        
        # 计算最大值
        max_force = max(abs(forces))
        max_moment = max(abs(moments))
        max_drag_force = max(abs(drag_forces))
        max_inertia_force = max(abs(inertia_forces))
        max_drag_moment = max(abs(drag_moments))
        max_inertia_moment = max(abs(inertia_moments))
        
        return {
            'max_force': max_force,
            'max_moment': max_moment,
            'max_drag_force': max_drag_force,
            'max_inertia_force': max_inertia_force,
            'max_drag_moment': max_drag_moment,
            'max_inertia_moment': max_inertia_moment,
            'force_ratio': max_drag_force / max_force * 100,
            'moment_ratio': max_drag_moment / max_moment * 100,
            't_values': t_values,
            'forces': forces,
            'moments': moments,
            'drag_forces': drag_forces,
            'inertia_forces': inertia_forces,
            'drag_moments': drag_moments,
            'inertia_moments': inertia_moments
        }

    def plot_results(self, wave_period, wave_height):
        """绘制波浪力和力矩随时间的变化"""
        results = self.calculate_max_values(wave_period, wave_height, d=self.depth)
        
        t_values = results['t_values']
        forces = results['forces']
        moments = results['moments']
        drag_forces = results['drag_forces']
        inertia_forces = results['inertia_forces']
        drag_moments = results['drag_moments']
        inertia_moments = results['inertia_moments']
        
        # 转换为kN和kN·m
        forces_kN = forces / 1000
        moments_kNm = moments / 1000
        drag_forces_kN = drag_forces / 1000
        inertia_forces_kN = inertia_forces / 1000
        drag_moments_kNm = drag_moments / 1000
        inertia_moments_kNm = inertia_moments / 1000
        
        # 绘制波浪力随时间的变化
        plt.figure(figsize=(15, 10))
        
        plt.subplot(2, 1, 1)
        plt.plot(t_values, forces_kN, 'b-', label='总波浪力')
        plt.plot(t_values, drag_forces_kN, 'r--', label='拖曳力')
        plt.plot(t_values, inertia_forces_kN, 'g--', label='惯性力')
        plt.xlabel('时间 (s)')
        plt.ylabel('波浪力 (kN)')
        plt.title(f'波浪力随时间的变化 (波高={wave_height}m, 周期={wave_period}s)')
        plt.grid(True)
        plt.legend()
        
        # 绘制波浪力矩随时间的变化
        plt.subplot(2, 1, 2)
        plt.plot(t_values, moments_kNm, 'b-', label='总波浪力矩')
        plt.plot(t_values, drag_moments_kNm, 'r--', label='拖曳力矩')
        plt.plot(t_values, inertia_moments_kNm, 'g--', label='惯性力矩')
        plt.xlabel('时间 (s)')
        plt.ylabel('波浪力矩 (kN·m)')
        plt.title(f'波浪力矩随时间的变化 (波高={wave_height}m, 周期={wave_period}s)')
        plt.grid(True)
        plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        # 输出最大值
        print(f"最大总波浪力: {results['max_force']/1000:.2f} kN")
        print(f"最大拖曳力: {results['max_drag_force']/1000:.2f} kN ({results['force_ratio']:.2f}%)")
        print(f"最大惯性力: {results['max_inertia_force']/1000:.2f} kN ({100-results['force_ratio']:.2f}%)")
        print(f"最大总波浪力矩: {results['max_moment']/1000:.2f} kN·m")
        print(f"最大拖曳力矩: {results['max_drag_moment']/1000:.2f} kN·m ({results['moment_ratio']:.2f}%)")
        print(f"最大惯性力矩: {results['max_inertia_moment']/1000:.2f} kN·m ({100-results['moment_ratio']:.2f}%)")

    def analyze_parameter_effects(self, base_height=5.0, base_period=10.0):
        """分析波高和周期对波浪荷载的影响"""
        # 波高变化范围
        heights = np.linspace(1, 20, 20)  # 2m到8m，7个点
        # 周期变化范围
        periods = np.linspace(1, 20, 20)  # 6s到14s，9个点
        
        # 存储结果
        height_effects = {
            'forces': [],
            'moments': [],
            'drag_forces': [],
            'inertia_forces': [],
            'drag_moments': [],
            'inertia_moments': []
        }
        
        period_effects = {
            'forces': [],
            'moments': [],
            'drag_forces': [],
            'inertia_forces': [],
            'drag_moments': [],
            'inertia_moments': []
        }
        
        # 分析波高的影响
        for H in heights:
            results = self.calculate_max_values(base_period, H, d=self.depth)
            height_effects['forces'].append(results['max_force']/1000)  # 转换为kN
            height_effects['moments'].append(results['max_moment']/1000)  # 转换为kN·m
            height_effects['drag_forces'].append(results['max_drag_force']/1000)
            height_effects['inertia_forces'].append(results['max_inertia_force']/1000)
            height_effects['drag_moments'].append(results['max_drag_moment']/1000)
            height_effects['inertia_moments'].append(results['max_inertia_moment']/1000)
        
        # 分析周期的影响
        for T in periods:
            results = self.calculate_max_values(T, base_height, d=self.depth)
            period_effects['forces'].append(results['max_force']/1000)
            period_effects['moments'].append(results['max_moment']/1000)
            period_effects['drag_forces'].append(results['max_drag_force']/1000)
            period_effects['inertia_forces'].append(results['max_inertia_force']/1000)
            period_effects['drag_moments'].append(results['max_drag_moment']/1000)
            period_effects['inertia_moments'].append(results['max_inertia_moment']/1000)
        
        # 绘制波高影响图
        plt.figure(figsize=(15, 10))
        plt.subplot(2, 1, 1)
        plt.plot(heights, height_effects['forces'], 'b-', label='总波浪力')
        plt.plot(heights, height_effects['drag_forces'], 'r--', label='拖曳力')
        plt.plot(heights, height_effects['inertia_forces'], 'g--', label='惯性力')
        plt.xlabel('波高 (m)')
        plt.ylabel('波浪力峰值 (kN)')
        plt.title(f'波高对波浪力峰值的影响 (周期={base_period}s)')
        plt.grid(True)
        plt.legend()
        
        plt.subplot(2, 1, 2)
        plt.plot(heights, height_effects['moments'], 'b-', label='总波浪力矩')
        plt.plot(heights, height_effects['drag_moments'], 'r--', label='拖曳力矩')
        plt.plot(heights, height_effects['inertia_moments'], 'g--', label='惯性力矩')
        plt.xlabel('波高 (m)')
        plt.ylabel('波浪力矩峰值 (kN·m)')
        plt.title(f'波高对波浪力矩峰值的影响 (周期={base_period}s)')
        plt.grid(True)
        plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        # 绘制周期影响图
        plt.figure(figsize=(15, 10))
        plt.subplot(2, 1, 1)
        plt.plot(periods, period_effects['forces'], 'b-', label='总波浪力')
        plt.plot(periods, period_effects['drag_forces'], 'r--', label='拖曳力')
        plt.plot(periods, period_effects['inertia_forces'], 'g--', label='惯性力')
        plt.xlabel('周期 (s)')
        plt.ylabel('波浪力峰值 (kN)')
        plt.title(f'周期对波浪力峰值的影响 (波高={base_height}m)')
        plt.grid(True)
        plt.legend()
        
        plt.subplot(2, 1, 2)
        plt.plot(periods, period_effects['moments'], 'b-', label='总波浪力矩')
        plt.plot(periods, period_effects['drag_moments'], 'r--', label='拖曳力矩')
        plt.plot(periods, period_effects['inertia_moments'], 'g--', label='惯性力矩')
        plt.xlabel('周期 (s)')
        plt.ylabel('波浪力矩峰值 (kN·m)')
        plt.title(f'周期对波浪力矩峰值的影响 (波高={base_height}m)')
        plt.grid(True)
        plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        return height_effects, period_effects

# 使用示例
if __name__ == "__main__":
    # 创建计算器实例
    calculator = WaveForceCalculator(
        depth=30.0,      # 水深(m)
        diameter=1.0,    # 直径(m)
        rho=1025.0,      # 海水密度(kg/m³)
        cd=1.0,          # 拖曳力系数
        cm=2.0           # 惯性力系数
    )
    
    # 分析波高和周期的影响
    calculator.analyze_parameter_effects(base_height=5.0, base_period=10.0) 

代码中参考了yl_YANG Liu所分享的五阶stokes波MATLAB代码，这里附上原文链接
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。原文链接：https://blog.csdn.net/yl_shqy/article/details/143161268