N=500;T31=10;H31=6;gamma=3.3;
Tp=T31/(1-0.132*(gamma+0.2)**(-0.559))
fp=1/Tp
fL=0.01;fH=5*fp;d=25
delatf=(fH-fL)/N
f=np.arange(fL,fH,delatf)
def sigma(f,fp):
sigma=np.piecewise(f,[f<=fp,f>fp],[0.07,0.09])
return sigma
def k_solve(f,d):
T=1/f
mu0=(2*np.pi)**2*d/(9.8*T**2)
mu=mu0*(1+mu0*np.exp(-(1.835+1.225*mu0**1.35)))/np.sqrt(np.tanh(mu0))
k=(2*np.pi)/(2*np.pi*d/mu)
return k
betaJ=0.0624*(1.094-0.01915*np.log(gamma))/(0.23+0.0336*gamma-0.185*(1.9+gamma)**(-1))
sigma1=sigma(f,fp)
sgm=np.exp(-(Tp*f-1)**2/2*sigma1**2)
Sf=betaJ*H31**2*Tp**(-4)*(f**(-5))*np.exp(-1.25*(Tp*f)**(-4))*(gamma**sgm)
h=np.sqrt(2*Sf*delatf)
k=k_solve(f,d)
xx=np.arange(0,600,0.5)
tt=np.arange(800,900,0.5)
x,t=np.meshgrid(xx,tt)
theta=np.random.uniform(0,2*np.pi,N)
C=np.cos(np.multiply.outer(x,k)-np.multiply.outer(t,2*np.pi*f)+theta)/np.cosh(k*d)
print(C.shape)
D=C*h
print(D.shape)
result=np.sum(D,axis=2)*1e3*10/1000
print(result.shape)
import matplotlib.pyplot as plt
######################设置全局字体######################
plt.rcParams['font.family']='Times new roman'
#字体大小
plt.rcParams['font.size']=50
#stix
plt.rcParams['mathtext.fontset'] = 'stix'
%matplotlib inline
fig, ax = plt.subplots(subplot_kw={"projection": "3d"},figsize=(6,6))
surf = ax.plot_surface(x, t, result, cmap='hsv')
# ax.ticklabel_format(style='sci', scilimits=(-1,2), axis='y')
ax.view_init(25,20)
#添加坐标轴标题
ax.set_xlabel(r'$x\ (m)$',labelpad=5,fontsize=15)
ax.set_ylabel(r'$t\ (s)$',labelpad=5,fontsize=15)
ax.set_zlabel(r'$Wave\ pressure\ (kPa)$',labelpad=5,fontsize=15)
ax.tick_params(axis='x',pad=1,labelsize=15)
ax.tick_params(axis='y',pad=1,labelsize=15)
ax.tick_params(axis='z',pad=1,labelsize=15)
from matplotlib.pyplot import MultipleLocator
x_major_locator=MultipleLocator(200)
ax.xaxis.set_major_locator(x_major_locator)
# ax.zlabel.set_rotation(10)
cbar=fig.colorbar(surf, shrink=0.4, aspect=10,orientation='horizontal',pad=0)
cbar.ax.tick_params(labelsize=14)
plt.tight_layout()
# ax.set_box_aspect((3,5,2))
ax.tick_params(which='major',direction='in',axis='x',width=2,length=4)
# fig.savefig(r'三维波压力1.png',bbox_inches='tight')
plt.show(block=True)
