一、仿真结果

1. 二维：波浪谱

2.二维：海浪波高随着时间的变化

3.三维：波浪谱

4. 三维：t=0时海浪波高

5.三维：海浪波高随时间的变化

二、代码

#%%

"""
%三维不规则短峰波的仿真
%参考论文：https://www.docin.com/p-1756913729.html
%仿真过程：1.根据风速确定P-M谱和仿真频率范围
%          2.求各个谐波的幅值
%          3.求0-2π之间均匀分布的随机数作为各个谐波的初相位
%          4.将各个谐波叠加
%参数：1.风速：10m/s;2.仿真频段：0.35-2.75；3.频率增量：0.08

"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
# ###############  1.设置参数  ################


# ######  v,g  ######

v=10
g=9.8

# ######  w  ######

dw=0.01
w_min = 0.35
w_max = 2.75
dw = 0.1
w = [w_min + dw*i for i in range(int((w_max - w_min) / dw) + 1)]
w = np.array([w])
print('w=')
print(w)

# ######  t  #######

t_min = 0
dt = 0.01
t_max = 50
T = [t_min + dt*i for i in range(int((t_max - t_min) / dt) + 1)]
T = np.array([T])
print('T=')
print(T)

# ######  X  #######

x_min = 0
dx = 1
x_max = 70
X = [t_min + dx*i for i in range(int((x_max - x_min) / dx) + 1)]
X = np.array([X])
print('X=')
print(X)

# ######  Y  #######

y_min = 0
dy = 1
y_max = 70
Y= [t_min + dy*i for i in range(int((y_max - y_min) / dy) + 1)]
Y = np.array([Y])
print('Y=')
print(Y)

# ######  num i,j  ######

num_i=np.size(w)
num_j=20
num_t=np.size(T)
num_x=np.size(X)
num_y=np.size(Y)
print('num_i=', num_i)
print('num_j=', num_j)
print('num_t=', num_t)
print('num_x=', num_x)
print('num_y=', num_y)

# ######  u  #######

u = np.random.rand(1,num_j) * np.pi - np.pi / 2 
print('u=')
print(u)
print(np.shape(u))

def du(uu):
    du = np.zeros((1,num_j))
    for j in range(num_j-1):
        du[0,j+1] = abs(uu[0,j+1] - uu[0,j])
    return du
du_ = du(u)
print('du=')
print(du_)
print(np.shape(du_))

#%%

# ###########  2.定义子函数  #################


# ######  扩散函数  #######

def phi(u, num_j):
    phi = np.zeros((1, num_j))
    for j in range(num_j - 1):
        phi[0, j] = 2 * ((np.cos(u[0, j])) ** 2) /np.pi
    return phi
phi_ = phi(u, num_j)
print('phi=')
print(phi_)

#%%

# ######  波数 k_i  #######


def k_i(num_i, w, g):
    k_i = np.zeros((1, num_i))
    for i in range(num_i - 1):
        k_i[0, i] = (w[0, i] ** 2) / g
    return k_i
k_i_ = k_i(num_i, w, g)
print('k_i=')
print(k_i_)

#%%

# ######  初相位  ######

def epsillon(num_i, num_j):
    epsillon = np.random.rand(num_i, num_j) * 2 * np.pi
    return epsillon
epsillon_ = epsillon(num_i, num_j)
print('epsillon=')
print(epsillon_)

#%%

# ######  P-M谱  #######

# ###  长峰波  ###
def S_zeta_l(v, w, g):
    S_zeta_l = np.zeros((1, np.size(w)))
    for i in range(np.size(w) - 1):
        S_zeta_l[0, i] = 8.1 * (10 ** (-3)) * g * g * np.exp(-0.74 * ((g /(v * w[0, i])) ** 4)) / (w[0, i] ** 5)
    return S_zeta_l
S_zeta_long = S_zeta_l(v, w, g)
print('长峰波波谱，S_zeta_long=')
print(S_zeta_long)
print(np.shape(S_zeta_long))

#%%

# ###  短峰波  ###

def S_zeta_s(S_zeta_long, phi_, num_i, num_j):
    S_zeta_s = np.zeros((num_i, num_j))
    for i in range(num_i - 1):
        for j in range(num_j - 1):
            S_zeta_s[i, j] = S_zeta_long[0, i] * phi_[0, j]
    return S_zeta_s
S_zeta_short = S_zeta_s(S_zeta_long, phi_, num_i, num_j)
print('短峰波，S_zeta_short=')
print(S_zeta_short)
print(np.shape(S_zeta_short))

#%%

# ######  时域波高  ######

# ### 长峰波 ###
def zeta_ai(S_zeta_long, dw, num_i):
    zeta_ai = np.zeros((1, num_i))
    for i in range(num_i - 1):
        zeta_ai[0, i] = np.sqrt(2 * S_zeta_long[0, i] * dw)
    return zeta_ai
zeta_ai_ = zeta_ai(S_zeta_long, dw, num_i)
print('长峰波，波高，zeta_ai=')
print(zeta_ai_)
print(np.shape(zeta_ai_))

#%%

# ### 短峰波 ###

def zeta_aij(S_zeta_short, dw, du_, num_i, num_j):
    zeta_aij = np.zeros((num_i, num_j))
    for i in range(num_i - 1):
        for j in range(num_j - 1):
            zeta_aij[i, j] = np.sqrt(2 * S_zeta_short[i, j] * dw * du_[0, j])
    return zeta_aij
zeta_aij_ = zeta_aij(S_zeta_short, dw, du_, num_i, num_j)
print('短峰波，波高，zeta_aij=')
print(zeta_aij_)
print(np.shape(zeta_aij_))

#%%

# ######  长峰波波高与时间的关系  ######

def zeta_t(zeta_ai_, epsillon_, T, num_i, num_t):
    zeta_t = np.zeros((1,num_t))
    for t in range(num_t - 1):
        for i in range(num_i -1):
            zeta_t[0, t] = zeta_t[0, t] + zeta_ai_[0, i] * np.cos(w[0, i] * T[0, t] + epsillon_[i, 0])
    return zeta_t
zeta_t_ = zeta_t(zeta_ai_, epsillon_, T, num_i, num_t)
print('zeta_t=')
print(zeta_t_)
print(np.shape(zeta_t_))

#%%

# ######  短峰波波高在不同位置（x, y）与时间的关系  ######

def zeta_xyt(X, Y, t, zeta_aij_, k_i_, w, epsillon_, num_x, num_y, num_t, num_i, num_j):
    zeta_xyt = np.zeros((num_x, num_y))
    
    for x in range(num_x -1):
        for y in range(num_y -1):
            for i in range(num_i -1):
                for j in range(num_j -1):
                    zeta_xyt[x, y] = zeta_xyt[x, y] + zeta_aij_[i, j] * np.cos(k_i_[0,i] * X[0,x] * np.cos(u[0, j]) + k_i_[0, i] * Y[0, y] * np.sin(u[0, j]) - w[0, i] * t + epsillon_[i, j])
    return zeta_xyt
zeta_xyt_ = zeta_xyt(X, Y, 5.001, zeta_aij_, k_i_, w, epsillon_, num_x, num_y, num_t, num_i, num_j)
print('zeta_xyt=')
print(zeta_xyt_)
print(np.shape(zeta_xyt_))

#%%

# ######  三、画图  ######


# ### 1. 2d ###
print(np.shape(w))
print(np.shape(S_zeta_long))
ww = np.linspace(w_min, w_max, num_i)
print(np.shape(ww))
f = interp1d(ww, S_zeta_long, kind='cubic')
S_zeta_long_new = f(ww)
plt.plot(ww, S_zeta_long_new.T)
plt.show()

#%%

print(np.shape(T))
print(np.shape(zeta_t_))
tt = np.linspace(t_min, t_max, num_t)
print(np.shape(tt))
f = interp1d(tt, zeta_t_, kind='cubic')
zeta_t_new = f(tt)
plt.plot(tt, zeta_t_new.T)
plt.show()

#%%

# ### 2. 3d ###

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

W, U = np.meshgrid(w, u.T)
Z = S_zeta_short

ax.plot_surface(W, U, Z.T)

plt.show()

#%%

# fig = plt.figure(figsize=(200, 200))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

XX, YY = np.meshgrid(X, Y.T)
Z = zeta_xyt_

print(np.shape(Z))
ax.plot_surface(XX, YY, Z.T)
plt.savefig("filename2.png") # 保存图片
plt.show()


#%%

# ######  短峰波波高在不同位置（x, y）与时间的关系  ######

def zeta_xyt(X, Y, T, zeta_aij_, k_i_, w, epsillon_, num_x, num_y, num_t, num_i, num_j):
    zeta_xyt = np.zeros((num_t, num_x, num_y))
    for t in range(num_t -1):
        for x in range(num_x -1):
            for y in range(num_y -1):
                for i in range(num_i -1):
                    for j in range(num_j -1):
                        zeta_xyt[t, x, y] = zeta_xyt[t, x, y] + zeta_aij_[i, j] * np.cos(k_i_[0,i] * X[0,x] * np.cos(u[0, j]) + k_i_[0, i] * Y[0, y] * np.sin(u[0, j]) - w[0, i] * T[0, t] + epsillon_[i, j])
    return zeta_xyt
zeta_xyt_ = zeta_xyt(X, Y, T, zeta_aij_, k_i_, w, epsillon_, num_x, num_y, num_t, num_i, num_j)
print('zeta_xyt=')
print(zeta_xyt_)
print(np.shape(zeta_xyt_))

#%%

np.save("zeta_xyt.npy",zeta_xyt_)
zeta_xyt_load = np.load("zeta_xyt.npy")
print(np.shape(zeta_xyt_load))

#%%

from PIL import Image
import matplotlib.pyplot as plt
from PIL import ImageFile
ImageFile.LOAD_TRUNCATED_IMAGES = True
Image.MAX_IMAGE_PIXELS = None



XX, YY = np.meshgrid(X, Y.T)
num_test = 1000

# 画图并保存
for t in range(num_test - 1):
    fig = plt.figure(figsize=(50, 50))
    ax = fig.add_subplot(111, projection='3d')
    Z_ = zeta_xyt_[t,:,:]
    print(Z_)
    surf = ax.plot_surface(XX, YY, Z_.T)
    filename = "images/t={}.png".format(T[0, t])
    
    plt.savefig(filename) # 保存图片
    plt.show()
    



#%%
import numpy as np
import cv2

# 读取一张图片
size = (5000, 5000)
num_test = 300 #这里因为内存问题没有办法将所有的图片都做成视频
fourcc = cv2.VideoWriter_fourcc('M','J','P','G')
print(size)
# 完成写入对象的创建，第一个参数是合成之后的视频的名称，第二个参数是可以使用的编码器，第三个参数是帧率即每秒钟展示多少张图片，第四个参数
videowrite = cv2.VideoWriter("movie_10s.avi", fourcc, 20, size)  # 20是帧数，size是图片尺寸
img_array = []
for filename in ["images/t={}.png".format(T[0, t]) for t in range(num_test - 1)]:
    img = cv2.imread(filename)
    img_array.append(img)

for t in range(num_test - 1):
    videowrite.write(img_array[t])
print('end!')