波浪谱荷载的Python实现

最新推荐文章于 2023-07-14 11:36:15 发布

海苏-

最新推荐文章于 2023-07-14 11:36:15 发布

阅读量999

点赞数

文章标签： python 开发语言

本文链接：https://blog.csdn.net/weixin_46258164/article/details/128009623

版权

众所周知，在海工结构分析中，波浪荷载是的表达和施加是研究的热点。根据平稳随机过程的谐波合成法，波浪荷载F(t)可按如下方式生成：

$F(t)= \sum_{n=1}^{m}\sqrt{2S_{ff}(w_{n})\Delta w}cos(\tilde{w}t+\varepsilon _{n})$

其中 $S_{ff}$ 被称为波浪荷载谱， $S_{ff}$ 由总速度力谱 $S_{FD}$ 和总惯性力谱 $S_{FM}$ 相加得

$S_{FF}(w)=S_{FD}(w)+S_{FM}(w)$

$S_{FD}$ 的表达式如下：

$S_{FD}(w)=\frac{1}{8\pi }(\frac{1}{2}\rho C_{D}DgkH_{s})^{2}[\frac{sinh2kd+2kd}{sin2kd}]^{2}S_{\eta \eta }(w)$

$S_{FM}$ 的表达式如下：

$S_{FM}(w)=[\frac{\rho g }{4}C_{M}\pi D^{2}tanhkd]^{2}S_{\eta \eta }(w)$

$S_{\eta \eta }$ 为谱密度函数，w为原频率， $\varepsilon _{n}$ 为随机初相位，d为水深，k为波数

w与k满足色散关系式如下：

$w_{n}=gk_{n}tanhk_{n}d$

实际海洋环境中波浪不是固定值，无数个圆频率叠加。此次分析采用等分频率法模拟海浪。确定 $w_{1}$ 和 $w_{m+1}$ ，即频率最小值和最大值。再确定m值，一般取50—100。

$w_{1}$ 和 $w_{m+1}$ 的表达式如下：

$w_{1}=(-\frac{3.11}{H_{s}^{2}ln(\mu )})^{1/4}$

$w_{m+1}=(-\frac{3.11}{H_{s}^{2}ln(1-\mu )})^{1/4}$

接下来利用python进行实现，因为需要用到函数，需要导入math模块

import math
import random
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

然后定义函数，由公式可知，函数需要定义的参数较多。

首先将w0和wmplus1用公式直接定谱密度函数、总速度力谱、总惯性力谱等与圆频率值有关，先定义零值数组。引入for循环。

def hanshu(z1,z2,d,rho,CM,CD,D,Hs,m):

    w0 = (-3.11 / (Hs ** 2 * math.log(0.002))) ** (1 / 4)
    wmplus1 = (-3.11 / (Hs ** 2 * math.log(1 - 0.002))) ** (1 / 4)
    w = np.zeros(m)
    Seta = np.zeros(m)
    Svv = np.zeros(m)
    Saa = np.zeros(m)
    F = np.zeros(m)
    eps = np.zeros(m)

    w[0] = w0

    dw = (wmplus1 - w0) / m
    for i in range(m):
       w[i] = w0+(i-1)*dw
       Seta[i] = 0.78 * w[i] ** 5 * math.exp(-3.11 / (w[i] ** 4 * Hs ** 2))
       Svv[i] =1 / (8 * math.pi) * ((0.5 * CD * rho * D) ** 2) * (9.8 ** 2) * w[i] ** 2 / (9.8 * d) * Hs ** 2 * ((2 * math.sqrt(
           w[i] ** 2 / (9.8 * d)) * z2 - 2 * math.sqrt(w[i] ** 2 / (9.8 * d)) * z1 + math.sinh(2 * math.sqrt(w[i] ** 2 / (9.8 * d)) * z2) - math.sinh(
           2 * math.sqrt(w[i] ** 2 / (9.8 * d)) * z1)) / (math.sinh(math.sqrt(w[i] ** 2 / (9.8 * d)) * d)) ** 2) ** 2 * Seta[i]
       Saa[i] =(rho / 4 * 9.8 * CM * math.pi * (D ** 2) * (math.sinh(math.sqrt(w[i] ** 2 / (9.8 * d)) * z2) - math.sinh(math.sqrt(w[i] ** 2 / (9.8 * d)) * z1)) / math.cosh(
           math.sqrt(w[i] ** 2 / (9.8 * d)) * d)) ** 2 * Seta[i]
       F[i] = Svv[i] + Saa[i]

  

    return w,Seta,Svv,Saa,F

定义完函数后，可以将初值代入，求得不同圆频率下的谱密度。根据公式知，波浪荷载是关于时间的函数，是所有原频率下时间函数的和。因为时间t是等差数列，所以可以不用for循环嵌套，而是直接采用np.linspace等差函数实现。

步骤分两步：1.先求出所有圆频率下t时刻的函数值 2.将对应的时间下的荷载值求和

第一步的具体操作是：先定义一个空数组y，采用y.append将圆频率数组扩充。

第二步直接采用sum函数求和。

代码中eps数组调用了随机函数，是为了模拟随机初相位 $\varepsilon _{n}$

if __name__ == '__main__':
    w0 = (-3.11 / (11.6 ** 2 * math.log(0.002))) ** (1 / 4)
    wmplus1 = (-3.11 / (11.6 ** 2 * math.log(1 - 0.002))) ** (1 / 4)
    m = 50
    dw = (wmplus1 - w0) / m
    Hs = 11.6
    z1 = 0
    z2 = 61
    d = 61
    rho = 1031
    CM = 2
    CD = 1.2
    D = 5.5
    w, Seta, Svv, Saa, F  = hanshu(z1, z2, d, rho, CM, CD, D, Hs, m)
    y = []
    Ft = []
    t = np.linspace(0.1, 20, 200)
    Farray = np.array(F)
    eps = np.zeros(m)
    for j in range(m):
        eps[j] = 2 * math.pi * random.uniform(1, m)
        y.append(math.sqrt(2*dw*F[j]/2)*np.cos(w[j]*t+eps[j]*t))
        Ft = sum(y)
print(Ft)
plt.plot(t,Ft)
plt.show()

随机波浪谱图形如下图所示：