SJ-G12-Sundry-聚焦波-波浪参数-ITTC波谱-代码
FLUENT UDF 代码
freakBAN75.C
// freakBAN75.c
/*
inlet moving mesh : paddles
outlet damping
fluent udf
*/
#include<stdio.h>
#include"udf.h"
#define N 75
#define Xp 10.0
#define Tp 15.0
#define pi M_PI
#define g 9.81
#define h 3
#define L0 20
#define L1 30
#define dw (4.83459-1.17196)/N
#define LW 5
DEFINE_CG_MOTION(ban_moving, dt, cg_vel, cg_omega, time, dtime)
{
real Sn[75] = { 3.98E-05,7.69E-05,0.000128005,0.000189661,0.000256353,0.000322136,0.0003819,0.000432067,0.000470719,0.000497378,0.000512633,0.000517739,0.000514288,0.000503961,0.000488366,0.000468946,0.000446936,0.000423355,0.000399016,0.000374549,0.000350429,0.000327,0.000304501,0.000283088,0.000262853,0.000243841,0.000226057,0.000209483,0.00019408,0.0001798,0.000166583,0.000154367,0.000143089,0.000132685,0.000123091,0.000114247,0.000106097,9.86E-05,9.17E-05,8.53E-05,7.94E-05,7.40E-05,6.90E-05,6.43E-05,6.00E-05,5.61E-05,5.24E-05,4.91E-05,4.59E-05,4.30E-05,4.03E-05,3.78E-05,3.55E-05,3.34E-05,3.14E-05,2.95E-05,2.78E-05,2.62E-05,2.46E-05,2.32E-05,2.19E-05,2.07E-05,1.96E-05,1.85E-05,1.75E-05,1.66E-05,1.57E-05,1.48E-05,1.41E-05,1.33E-05,1.27E-05,1.20E-05,1.14E-05,1.09E-05,1.03E-05 };
real k[75] = { 0.232371455840487,0.249786814149107,0.254149656827340,0.274659964303511,0.285474038790606,0.297912164299048,0.313923158777687,0.321324277725442,0.332577950485603,0.345229949359504,0.360209632718581,0.373869700821059,0.386704443897790,0.400905988782487,0.427237180995047,0.431434306971425,0.450262412027362,0.476991100267958,0.482080568587278,0.509089000706641,0.523833411249025,0.542565682494764,0.557217360473855,0.575792813470086,0.608084429196712,0.617619269403596,0.656201811860559,0.679729793272655,0.686334007551193,0.718583444405394,0.748198795453447,0.768706480317268,0.785281925083144,0.830053292544376,0.837797839083065,0.872058923626267,0.908616475332218,0.934559436401070,0.978228744079897,0.998367031019269,1.03032581517445,1.05822982277296,1.08519853842522,1.12335645542617,1.16763283687379,1.17385947072239,1.21893049141042,1.24385049120768,1.29256614364529,1.32825277957608,1.35830589618821,1.39094138295851,1.44126503331233,1.48426098449304,1.52488396596574,1.57111612906692,1.60899440375169,1.63502351616304,1.67799628996486,1.70742018650009,1.74643523019032,1.79116943440899,1.83888251717798,1.89037051470502,1.91890224124605,1.96465747202773,2.01779476405335,2.06291255223232,2.11822833687560,2.13575483404363,2.19861129869732,2.26221147909793,2.29036169441906,2.35191263607766,2.38270458677882 };
real w, Tr, A, veltemp, T1;
real f = 0;
int i;
for (i = 0; i < N; i++)
{
w = pow((g * k[i] * tanh(k[i] * h)), 0.5);
A = pow((2 * Sn[i] * dw), 0.5);
Tr = 4 * sinh(k[i] * h) * sinh(k[i] * h) / (sinh(2 * k[i] * h) + (2 * k[i] * h));
veltemp = w * A / Tr * cos(k[i] * (-Xp) - w * (time - Tp));
f = f + veltemp;
}
w = pow((g * k[0] * tanh(k[0] * h)), 0.5);
T1 = 2 * pi / w;
if (time <= (2 * T1))
{
cg_vel[0] = f * time / (2 * T1);
}
else
{
cg_vel[0] = f;
}
}
/**********************************/
DEFINE_SOURCE(xom_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;
real mu;
C_CENTROID(x,c,t);
if(C_U(c,t)<=0)
{
mu=2*pow(LW,0.5)*pow(((x[0]-L0)/(L1-L0)),2)*C_R(c,t);
source=-mu*C_U(c,t);
dS[eqn]=-mu;
}
else
{
mu=0;
source=0;
dS[eqn]=0;
}
return source;
}
/**********************************/
DEFINE_SOURCE(zom_source,c,t,dS,eqn)
{
real x[ND_ND];
real source;
real mu;
C_CENTROID(x,c,t);
mu=10*pow(LW,0.5)*pow(((x[0]-L0)/(L1-L0)),2)*C_R(c,t);
source=-mu*C_W(c,t);
dS[eqn]=-mu;
return source;
}
/**********************************/
利用ITTC波谱，生成随机波浪参数
getWSn.C
#include <iostream>
#include<math.h>
#include<fstream>
#include <iomanip>
#include <random> //产生随机数功能
#include <stdlib.h> //使用暂停功能
#include <ctime>//设置种子
#define PI 3.1415926
#define M 150 //叠加波浪的个数
using namespace std;
//生成ITTC波浪谱
void ITTC(double H, double Tz, double* w, double* s_w)
{ // 通过H和T, 得到w(波浪圆频率)和sw(对应谱密度)
double As, Bs;//中间值
ofstream ofs1;
ofstream ofs2;
ofs1.open("ITTC_x轴.txt", ios::out);
ofs2.open("ITTC_y轴.txt", ios::out);
As = 173 * pow(H, 2) / pow(Tz, 4);
Bs = 691 / pow(Tz,4);
for (int i = 0; i < 40000; i++)
{
*w = 0.0005 *(i+1);
*s_w = As / pow(*w, 5)*exp(-(Bs / pow(*w, 4)));
ofs1 << fixed << setprecision(4) << *w << endl;
ofs2 << fixed << setprecision(10) << *s_w << endl;
w++;
s_w++;
}
ofs1.close();
ofs2.close();
}
//确定波浪谱上下限的圆频率
void frequency(double low, double high, double* w_low, double* w_high, double* w, double* s_w)
{ // low 最小能量比例; high 最大能量比例; 为了找到频率的上下限
double temp_area[40000];
double T_area=0; // 总面积 (能量)
double temp;
for (int i = 0; i < 40000; i++)
{
temp = 0.0005*(*s_w);
T_area = T_area + temp;
temp_area[i] = T_area;
w++;
s_w++;
}
cout << "总面积为："<< T_area << endl;
for (int i = 0; i < 40000; i++)
{
if (temp_area[i] >= T_area*low)
{ // 找到小于最小能量的
cout << "当前面积为"<< temp_area [i]<< endl;
*w_low = 0.0005 * (i + 1);
break;
}
}
for (int i = 0; i < 40000; i++)
{
if (temp_area[i] >= T_area*high)
{ // 找到大于最大能量的
cout << "当前面积为" << temp_area[i] << endl;
*w_high = 0.0005 *(i + 1);
break;
}
}
}
//生成随机组成波浪的代表频率
void random_frequency(double* random_w, double w_low, double w_high, double* w1)
{ //得到新的圆频率序列w1, 其中w1是间隔均匀的N个圆频率, random_w是频谱上随机的N个频率
//ofstream ofs1;
//ofs1.open("组成波浪的代表频率.txt", ios::out);
double temp[500] = {0};
double random[500] = {0};
double space_w = (w_high - w_low) / M;
temp[0] = w_low;
for (int i = 0; i < M+1; i++)
{
temp[i+1] = temp[i] + space_w;
}
for (int i = 0; i <M+1; i++)
{
*w1 = temp[i];
w1++;
}
srand((int)time(0));
for (int i = 0; i < M+1; i++)
{
random[i] = temp[i] + rand() / double(RAND_MAX)*space_w;
*random_w = random[i];
//ofs1 << fixed << setprecision(4) << random[i] << endl;
random_w++;
}
//ofs1.close();
}
//生成随机相位
void random_phase(double* random_p)
{
double temp[500] = { 0 };
srand((int)time(0));
for (int i = 0; i < M+1; i++)
{
temp[i] = rand() / double(RAND_MAX) *2*PI;
//cout << temp[i]<< endl;
*random_p = temp[i];
random_p++;
}
}
//组成波浪参数
void wave_h(double H, double Tz, double * A, double * w1, double w_low, double w_high, double* s_w1, double* w2)
{
double As, Bs;//中间值
double temp1[500];//接收w1的数组
double temp[500] = {0};//接收(w(i-1)+w(i))/2的数值
for (int i = 0; i < M+1; i++)
{
temp1[i] = *w1;
w1++;
}
for (int i = 0; i < M+1; i++)
{
temp[i] = (temp1[i] + temp1[i + 1])/2;
*w2 = temp[i];
w2++;
}
As = 173 * pow(H, 2) / pow(Tz, 4);
Bs = 691 / pow(Tz, 4);
for (int i = 0; i < M+1; i++)
{
*s_w1= As / pow(temp[i], 5)*exp(-(Bs / pow(temp[i], 4)));
*A = sqrt(2 * As / pow(temp[i], 5)*exp(-(Bs / pow(temp[i], 4)))*((w_high - w_low) / M));
// 此处修改了M
A++;
s_w1++;
}
}
int main()
{
double H=0;//有义波高 0.1
double Tz=0;//平均周期 2.8
double low = 0.002;//略去能量的下限
double w_low = 0.0;//频率下限
double high = 0.998;//略去能量的上限
double w_high = 0.0;//频率上限
double w[40000] = {0.0};//x轴频率坐标
double s_w[40000] = {0.0};//y轴谱密度坐标
double random_w[500] = { 0 };//组成波浪代表频率
double random_p[500] = { 0 };//组成波浪随机相位
double A[500] = { 0 };//组成波浪幅值
double w1[500] = { 0 };//w(i+1)与w(i)的数值
double s_w1[500] = { 0 };//保存S(w^)的数值
double out_f[500][5] = { 0 };
double w2[500];接收(w(i-1)+w(i))/2的数值
cout<<"请输入有义波高："<< endl;
cin >> H;
cout << "请输入平均周期：" << endl;
cin >> Tz;
//ITTC波浪谱计算
ITTC(H, Tz, w, s_w);
//确定波浪谱上下限的圆频率
frequency(low, high, &w_low, &w_high, w, s_w);
cout <<"频率下限为："<< w_low << endl;
cout << "频率上限为：" << w_high << endl;
//生成随机频率
random_frequency(random_w, w_low, w_high, w1);
//生成随机相位
random_phase(random_p);
//生成波浪幅值
wave_h(H, Tz, A, w1, w_low, w_high, s_w1, w2);
//输出到一个文件中
ofstream ofs1;
ofs1.open("组成波浪参数表.txt", ios::out);
for (int i = 0; i < M; i++)
{
out_f[i][0] = random_w[i];
out_f[i][1] = s_w1[i];
ofs1 << out_f[i][0] << " " << out_f[i][1] << endl;
}
system("pause");
}
利用给定的参数生成聚焦波参数
focus_wave.py
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation
class mesh:
__x_min__ = None
__x_max__ = None
__y_min__ = None
__y_max__ = None
__x_num__ = None
def __init__(self, dictionary_arg):
self.__x_min__ = dictionary_arg['x_min']
self.__x_max__ = dictionary_arg['x_max']
self.__y_min__ = dictionary_arg['y_min']
self.__y_max__ = dictionary_arg['y_max']
self.__x_num__ = dictionary_arg['x_num']
def __get_xp__(self):
num = self.__x_num__ + 1
xp = np.linspace(self.__x_min__, self.__x_max__, num)
return xp
def get_xx(self):
xp = self.__get_xp__()
xc = []
for ii_in in range(len(xp) - 1):
xm = (xp[ii_in] + xp[ii_in + 1]) / 2.0
xc.append(xm)
return xc
class wave_para_converter:
__water_depth__ = None
__wave_dw__ = None
__x_focus__ = None
__t_focus__ = None
__wave_sn__ = None
__wave_number__ = None
def __init__(self, dictionary_arg):
self.__water_depth__ = dictionary_arg['water_depth']
self.__wave_dw__ = dictionary_arg['wave_dw']
self.__x_focus__ = dictionary_arg['x_focus']
self.__t_focus__ = dictionary_arg['t_focus']
self.__wave_sn__ = dictionary_arg['wave_sn']
self.__wave_number__ = dictionary_arg['wave_number']
def wave_height(self):
sn = self.__wave_sn__
dw = self.__wave_dw__
wave_height = 2.0 * pow((2.0 * sn * dw), 0.5)
return wave_height
def wave_period(self):
g = 9.81
k = self.__wave_number__
h = self.__water_depth__
wave_omega = pow((g * k * math.tanh(k * h)), 0.5)
wave_period = (2.0 * math.pi) / wave_omega
return wave_period
def wave_phase(self):
k = self.__wave_number__
xp = self.__x_focus__
tp = self.__t_focus__
omega = (2.0 * math.pi) / (self.wave_period())
phase = - k * xp + omega * tp
return phase
class regular_wave:
__water_depth__ = None
__wave_height__ = None
__wave_phase__ = None
__wave_period__ = None
__wave_length__ = None
def __init__(self, dictionary_arg):
self.__water_depth__ = dictionary_arg['water_depth']
self.__wave_height__ = dictionary_arg['wave_height']
self.__wave_phase__ = dictionary_arg['wave_phase']
self.__wave_period__ = dictionary_arg['wave_period']
self.__wave_length__ = self.wave_length()
def wave_speed(self):
g_mag = 9.81
wave_number = self.wave_number()
water_depth = self.__water_depth__
wave_speed = math.sqrt((g_mag / wave_number) * math.tanh(wave_number * water_depth))
return wave_speed
def wave_length(self):
water_depth = self.__water_depth__
wave_period = self.__wave_period__
g = 9.81
wave_length0 = abs(g) * wave_period * wave_period / (2.0 * math.pi)
wave_length = wave_length0
# stokes I
for ii_in in range(1, 11): # optional: 11 or 101
wave_length = wave_length0 * math.tanh(2.0 * math.pi * water_depth / wave_length)
return wave_length
def wave_number(self):
return (2.0 * math.pi) / (self.wave_length())
def wave_omega(self):
return (2.0 * math.pi) / self.__wave_period__
def wave_surface(self, xx_arg, tt_arg):
wave_number = self.wave_number()
wave_omega = self.wave_omega()
wave_height = self.__wave_height__
wave_phase = self.__wave_phase__
theta00 = wave_number * xx_arg - wave_omega * tt_arg + wave_phase
wave_surface = (wave_height / 2.0) * math.cos(theta00)
wave_speed = self.wave_speed()
distance = wave_speed * tt_arg
if xx_arg > distance:
wave_surface = 0.0
return wave_surface
class focus_wave:
__wave_sn__ = None
__wave_number__ = None
__water_depth__ = None
__wave_dw__ = None
__x_focus__ = None
__t_focus__ = None
__n_focus__ = None
def __init__(self, dictionary_arg):
self.__wave_sn__ = dictionary_arg['wave_sn']
self.__wave_number__ = dictionary_arg['wave_number']
self.__water_depth__ = dictionary_arg['water_depth']
self.__wave_dw__ = dictionary_arg['wave_dw']
self.__x_focus__ = dictionary_arg['x_focus']
self.__t_focus__ = dictionary_arg['t_focus']
self.__n_focus__ = dictionary_arg['n_focus']
def wave_surface(self, xx_arg0, tt_arg0):
wave_surface = 0.0
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'wave_dw': self.__wave_dw__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_sn': self.__wave_sn__[ii_in],
'wave_number': self.__wave_number__[ii_in]
}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
regular_wave_dictionary = {
'water_depth': 3.0,
'wave_height': my_wpc.wave_height(),
'wave_period': my_wpc.wave_period(),
'wave_phase': my_wpc.wave_phase()
}
my_regular_wave = regular_wave(dictionary_arg=regular_wave_dictionary)
wave_surface_tmp = my_regular_wave.wave_surface(xx_arg=xx_arg0, tt_arg=tt_arg0)
wave_surface = wave_surface_tmp + wave_surface
return wave_surface
def write_paras(self):
my_out_file_name = './waveDict.txt'
my_out_file = open(my_out_file_name, 'w')
my_out_file.write('wavePeriods\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'wave_dw': self.__wave_dw__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_sn': self.__wave_sn__[ii_in],
'wave_number': self.__wave_number__[ii_in]
}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
my_out_file.write(str(my_wpc.wave_period()) + '\n')
my_out_file.write(');\n')
my_out_file.write('waveHeights\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'wave_dw': self.__wave_dw__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_sn': self.__wave_sn__[ii_in],
'wave_number': self.__wave_number__[ii_in]
}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
my_out_file.write(str(my_wpc.wave_height()) + '\n')
my_out_file.write(');\n')
my_out_file.write('wavePhases\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'wave_dw': self.__wave_dw__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_sn': self.__wave_sn__[ii_in],
'wave_number': self.__wave_number__[ii_in]
}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
wave_phase = my_wpc.wave_phase() # / 180.0 * math.pi
my_out_file.write(str(wave_phase) + '\n')
my_out_file.write(');\n')
my_out_file.write('waveDirs\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
my_out_file.write('0.0\n')
my_out_file.write(');\n')
my_out_file.flush()
if __name__ == '__main__':
mesh_dictionary = {
'x_min': 0,
'x_max': 30,
'y_min': -0.5,
'y_max': 0.5,
'x_num': 100 # optional: 100 or 3000
}
my_mesh = mesh(dictionary_arg=mesh_dictionary)
xx = my_mesh.get_xx()
focus_dictionary = {'water_depth': 3.0, 'wave_dw': (4.83459 - 1.17196) / 75.0, 'x_focus': 10.0, 't_focus': 15.0,
'n_focus': 75,
'wave_sn': [3.98E-05, 7.69E-05, 0.000128005, 0.000189661, 0.000256353, 0.000322136, 0.0003819,
0.000432067, 0.000470719, 0.000497378, 0.000512633, 0.000517739, 0.000514288,
0.000503961, 0.000488366, 0.000468946, 0.000446936, 0.000423355, 0.000399016,
0.000374549, 0.000350429, 0.000327, 0.000304501, 0.000283088, 0.000262853,
0.000243841, 0.000226057, 0.000209483, 0.00019408, 0.0001798, 0.000166583,
0.000154367, 0.000143089, 0.000132685, 0.000123091, 0.000114247, 0.000106097,
9.86E-05, 9.17E-05, 8.53E-05, 7.94E-05, 7.40E-05, 6.90E-05, 6.43E-05, 6.00E-05,
5.61E-05, 5.24E-05, 4.91E-05, 4.59E-05, 4.30E-05, 4.03E-05, 3.78E-05, 3.55E-05,
3.34E-05, 3.14E-05, 2.95E-05, 2.78E-05, 2.62E-05, 2.46E-05, 2.32E-05, 2.19E-05,
2.07E-05, 1.96E-05, 1.85E-05, 1.75E-05, 1.66E-05, 1.57E-05, 1.48E-05, 1.41E-05,
1.33E-05, 1.27E-05, 1.20E-05, 1.14E-05, 1.09E-05, 1.03E-05],
'wave_number': [0.232371455840487, 0.249786814149107, 0.254149656827340, 0.274659964303511,
0.285474038790606, 0.297912164299048, 0.313923158777687, 0.321324277725442,
0.332577950485603, 0.345229949359504, 0.360209632718581, 0.373869700821059,
0.386704443897790, 0.400905988782487, 0.427237180995047, 0.431434306971425,
0.450262412027362, 0.476991100267958, 0.482080568587278, 0.509089000706641,
0.523833411249025, 0.542565682494764, 0.557217360473855, 0.575792813470086,
0.608084429196712, 0.617619269403596, 0.656201811860559, 0.679729793272655,
0.686334007551193, 0.718583444405394, 0.748198795453447, 0.768706480317268,
0.785281925083144, 0.830053292544376, 0.837797839083065, 0.872058923626267,
0.908616475332218, 0.934559436401070, 0.978228744079897, 0.998367031019269,
1.03032581517445, 1.05822982277296, 1.08519853842522, 1.12335645542617,
1.16763283687379, 1.17385947072239, 1.21893049141042, 1.24385049120768,
1.29256614364529, 1.32825277957608, 1.35830589618821, 1.39094138295851,
1.44126503331233, 1.48426098449304, 1.52488396596574, 1.57111612906692,
1.60899440375169, 1.63502351616304, 1.67799628996486, 1.70742018650009,
1.74643523019032, 1.79116943440899, 1.83888251717798, 1.89037051470502,
1.91890224124605, 1.96465747202773, 2.01779476405335, 2.06291255223232,
2.11822833687560, 2.13575483404363, 2.19861129869732, 2.26221147909793,
2.29036169441906, 2.35191263607766, 2.38270458677882]}
my_focus = focus_wave(dictionary_arg=focus_dictionary)
my_focus.write_paras()
fig, ax = plt.subplots()
# tt = 15.0
# phi = []
# for ii in range(len(xx)):
# phi_tmp = my_focus.wave_surface(xx_arg0=xx[ii], tt_arg0=tt)
# phi.append(phi_tmp)
#
# ax.set_xlim([0.0, 30.0])
# ax.set_ylim([-0.3, 0.3])
# line = ax.plot(xx, phi, '-k')
# plt.show()
tt_max = 30.0 # optional: 30.0 or 50.0
tt_step = 0.3 # optional: 0.3 or 0.1
tt_min = 0.0
# for tt in np.arange(tt_min, tt_max, tt_step):
# plt.cla()
# phi = []
# for ii in range(len(xx)):
# phi_tmp = my_focus.wave_surface(xx_arg0=xx[ii], tt_arg0=tt)
# phi.append(phi_tmp)
#
# ax.set_xlim([0.0, 30.0])
# ax.set_ylim([-0.3, 0.3])
#
# line = ax.plot(xx, phi, '-k')
# plt.pause(0.001)
line, = ax.plot([], [])
ax.set_xlim([0.0, 30.0])
ax.set_ylim([-0.3, 0.3])
ims = []
for tt in np.arange(tt_min, tt_max, tt_step):
phi = []
for ii in range(len(xx)):
phi_tmp = my_focus.wave_surface(xx_arg0=xx[ii], tt_arg0=tt)
phi.append(phi_tmp)
line, = ax.plot(xx, phi, '-k')
ims.append([line])
print(str(tt / tt_max * 100.0) + '%\n', flush=True)
ani = animation.ArtistAnimation(fig, ims, interval=10, repeat=False)
writer = animation.PillowWriter(fps=30)
ani.save('wave.gif', writer=writer)
利用有义波高和平均周期生成聚焦波参数
getFocusWave.py
import matplotlib.pyplot as plt
import random
import time
import numpy as np
import math
import matplotlib.animation as animation
class ITTC_wave_spectral:
__min_wave_omega__ = 0.0
__max_wave_omega__ = 20.0
__num_wave_omega__ = 40000
__wave_energy_low__ = 0.002
__wave_energy_high__ = 0.998
__significant_wave_height__ = None
__mean_wave_period__ = None
__num_wave_out__ = None
def __init__(self, dictionary_arg):
self.__significant_wave_height__ = dictionary_arg['significant_wave_height']
self.__mean_wave_period__ = dictionary_arg['mean_wave_period']
self.__num_wave_out__ = dictionary_arg['num_wave_out']
def __space_wave_omega__(self):
max_wave_omega = self.__max_wave_omega__
min_wave_omega = self.__min_wave_omega__
num_wave_omega = self.__num_wave_omega__
space_wave_omega = (max_wave_omega - min_wave_omega) / num_wave_omega
return space_wave_omega
def __spectral_density__(self):
significant_wave_height = self.__significant_wave_height__
mean_wave_period = self.__mean_wave_period__
num_wave_omega = self.__num_wave_omega__
space_wave_omega = self.__space_wave_omega__()
a_as = 173.0 * pow(significant_wave_height, 2.0) / pow(mean_wave_period, 4.0)
b_bs = 691.0 / pow(mean_wave_period, 4.0)
spectral_density = []
wave_omega = []
for ii in range(num_wave_omega):
wave_omega_tmp = space_wave_omega * (ii + 1.0)
wave_omega.append(wave_omega_tmp)
spectral_density_tmp = a_as / pow(wave_omega_tmp, 5.0) * math.exp(-(b_bs / pow(wave_omega_tmp, 4.0)))
spectral_density.append(spectral_density_tmp)
return spectral_density
def __wave_omega__(self):
wave_omega = []
num_wave_omega = self.__num_wave_omega__
space_wave_omega = self.__space_wave_omega__()
for ii in range(num_wave_omega):
wave_omega_tmp = space_wave_omega * (ii + 1.0)
wave_omega.append(wave_omega_tmp)
return wave_omega
def plot_wave_spectral_full(self):
wave_omega = self.__wave_omega__()
spectral_density = self.__spectral_density__()
plt.plot(wave_omega, spectral_density)
plt.show()
def test_func(self):
sth = self.__wave_omega_random__()
print(sth)
def __total_energy__(self):
total_energy = 0.0
spectral_density = self.__spectral_density__()
space_wave_omega = self.__space_wave_omega__()
num_wave_omega = self.__num_wave_omega__
for ii in range(num_wave_omega):
total_energy_tmp = space_wave_omega * (spectral_density[ii])
total_energy += total_energy_tmp
return total_energy
def __wave_omega_low__(self):
wave_omega_low = 0.0
total_energy = self.__total_energy__()
spectral_density = self.__spectral_density__()
space_wave_omega = self.__space_wave_omega__()
num_wave_omega = self.__num_wave_omega__
wave_energy_low = self.__wave_energy_low__
energy_now = 0.0
for ii in range(num_wave_omega):
energy_tmp = space_wave_omega * (spectral_density[ii])
energy_now += energy_tmp
if energy_now >= total_energy * wave_energy_low:
wave_omega_low = space_wave_omega * (ii + 1)
break
return wave_omega_low
def __wave_omega_high__(self):
wave_omega_high = 0.0
total_energy = self.__total_energy__()
wave_energy_high = self.__wave_energy_high__
spectral_density = self.__spectral_density__()
space_wave_omega = self.__space_wave_omega__()
num_wave_omega = self.__num_wave_omega__
energy_now = 0.0
for ii in range(num_wave_omega):
energy_tmp = space_wave_omega * (spectral_density[ii])
energy_now += energy_tmp
if energy_now >= total_energy * wave_energy_high:
wave_omega_high = space_wave_omega * (ii + 1)
break
return wave_omega_high
def __space_wave_omega_random__(self):
wave_omega_low = self.__wave_omega_low__()
wave_omega_high = self.__wave_omega_high__()
num_wave_out = self.__num_wave_out__
space_wave_omega_random = (wave_omega_high - wave_omega_low) / num_wave_out
return space_wave_omega_random
def __wave_omega_homogeneous__(self):
wave_omega_low = self.__wave_omega_low__()
num_wave_out = self.__num_wave_out__
wave_omega_homogeneous = [wave_omega_low]
space_wave_omega_random = self.__space_wave_omega_random__()
for ii in range(num_wave_out + 1):
wave_omega_homogeneous_tmp = wave_omega_homogeneous[ii] + space_wave_omega_random
wave_omega_homogeneous.append(wave_omega_homogeneous_tmp)
return wave_omega_homogeneous
def __wave_omega_random__(self):
wave_omega_homogeneous = self.__wave_omega_homogeneous__()
space_wave_omega_random = self.__space_wave_omega_random__()
num_wave_out = self.__num_wave_out__
wave_omega_random = []
random.seed(time.time())
for ii in range(num_wave_out):
wave_omega_random_tmp = wave_omega_homogeneous[ii] + random.random() * space_wave_omega_random
wave_omega_random.append(wave_omega_random_tmp)
return wave_omega_random
def wave_omega(self):
wave_omega = self.__wave_omega_random__()
return wave_omega
def wave_height(self):
wave_omega_random = self.__wave_omega_random__()
wave_omega_low = self.__wave_omega_low__()
wave_omega_high = self.__wave_omega_high__()
significant_wave_height = self.__significant_wave_height__
mean_wave_period = self.__mean_wave_period__
a_as = 173.0 * pow(significant_wave_height, 2.0) / pow(mean_wave_period, 4.0)
b_bs = 691.0 / pow(mean_wave_period, 4.0)
num_wave_out = self.__num_wave_out__
wave_height = []
for ii in range(num_wave_out):
wave_height_tmp = math.sqrt(
2.0 * a_as / pow(wave_omega_random[ii], 5.0) * math.exp(-(b_bs / pow(wave_omega_random[ii], 4.0))) * (
(wave_omega_high - wave_omega_low) / num_wave_out))
wave_height.append(2.0 * wave_height_tmp)
return wave_height
class focus_wave:
__water_depth__ = None
__x_focus__ = None
__t_focus__ = None
__n_focus__ = None
__significant_wave_height__ = None
__mean_wave_period__ = None
def __init__(self, dictionary_arg):
self.__water_depth__ = dictionary_arg['water_depth']
self.__x_focus__ = dictionary_arg['x_focus']
self.__t_focus__ = dictionary_arg['t_focus']
self.__n_focus__ = dictionary_arg['n_focus']
self.__significant_wave_height__ = dictionary_arg['significant_wave_height']
self.__mean_wave_period__ = dictionary_arg['mean_wave_period']
def ittc_wave_spectral(self):
ittc_dictionary = {
'significant_wave_height': self.__significant_wave_height__,
'mean_wave_period': self.__mean_wave_period__,
'num_wave_out': self.__n_focus__}
my_ittc_in = ITTC_wave_spectral(dictionary_arg=ittc_dictionary)
return my_ittc_in
def wave_para_convertor(self, wave_height_arg, wave_omega_arg):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_height': wave_height_arg,
'wave_omega': wave_omega_arg}
my_wpc_in = wave_para_converter(dictionary_arg=wpc_dictionary)
return my_wpc_in
def write_paras(self):
my_out_file_name = './waveDict.txt'
my_out_file = open(my_out_file_name, 'w')
my_out_file.write('wavePeriods\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
ittc_dictionary = {
'significant_wave_height': self.__significant_wave_height__,
'mean_wave_period': self.__mean_wave_period__,
'num_wave_out': self.__n_focus__}
my_ittc = ITTC_wave_spectral(dictionary_arg=ittc_dictionary)
wave_height = my_ittc.wave_height()
wave_omega = my_ittc.wave_omega()
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_height': wave_height[ii_in],
'wave_omega': wave_omega[ii_in]}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
my_out_file.write(str(my_wpc.wave_period()) + '\n')
my_out_file.write(');\n')
my_out_file.write('waveHeights\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
my_out_file.write(str(wave_height[ii_in]) + '\n')
my_out_file.write(');\n')
my_out_file.write('wavePhases\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
wpc_dictionary = {
'water_depth': self.__water_depth__,
'x_focus': self.__x_focus__,
't_focus': self.__t_focus__,
'wave_height': wave_height[ii_in],
'wave_omega': wave_omega[ii_in]}
my_wpc = wave_para_converter(dictionary_arg=wpc_dictionary)
wave_phase = my_wpc.wave_phase() # / 180.0 * math.pi
my_out_file.write(str(wave_phase) + '\n')
my_out_file.write(');\n')
my_out_file.write('waveDirs\n')
my_out_file.write(str(self.__n_focus__) + '\n')
my_out_file.write('(\n')
for ii_in in range(self.__n_focus__):
my_out_file.write('0.0\n')
my_out_file.write(');\n')
my_out_file.flush()
print('write paras finished \n')
class wave_para_converter:
__water_depth__ = None
__x_focus__ = None
__t_focus__ = None
__wave_height__ = None
__wave_omega__ = None
def __init__(self, dictionary_arg):
self.__water_depth__ = dictionary_arg['water_depth']
self.__x_focus__ = dictionary_arg['x_focus']
self.__t_focus__ = dictionary_arg['t_focus']
self.__wave_height__ = dictionary_arg['wave_height']
self.__wave_omega__ = dictionary_arg['wave_omega']
def wave_period(self):
wave_omega = self.__wave_omega__
wave_period = (2.0 * math.pi) / wave_omega
return wave_period
def __wave_number__(self):
wave_length = self.__wave_length__()
wave_number = (2.0 * math.pi) / wave_length
return wave_number
def __wave_length__(self):
water_depth = self.__water_depth__
wave_period = self.wave_period()
g = 9.81
wave_length0 = abs(g) * wave_period * wave_period / (2.0 * math.pi)
wave_length = wave_length0
# stokes I
for ii_in in range(1, 11): # optional: 11 or 101
wave_length = wave_length0 * math.tanh(2.0 * math.pi * water_depth / wave_length)
return wave_length
def wave_phase(self):
wave_number = self.__wave_number__()
x_focus = self.__x_focus__
t_focus = self.__t_focus__
wave_omega = self.__wave_omega__
wave_phase = - wave_number * x_focus + wave_omega * t_focus
return wave_phase
def wave_surface(self, xx_arg, tt_arg):
wave_number = self.__wave_number__()
wave_omega = self.__wave_omega__
wave_height = self.__wave_height__
wave_phase = self.wave_phase()
theta00 = wave_number * xx_arg - wave_omega * tt_arg + wave_phase
wave_surface = (wave_height / 2.0) * math.cos(theta00)
wave_speed = self.__wave_speed__()
distance = wave_speed * tt_arg
if xx_arg > distance:
wave_surface = 0.0
return wave_surface
def __wave_speed__(self):
g_mag = 9.81
wave_number = self.__wave_number__()
water_depth = self.__water_depth__
wave_speed = math.sqrt((g_mag / wave_number) * math.tanh(wave_number * water_depth))
return wave_speed
class mesh:
__x_min__ = None
__x_max__ = None
__y_min__ = None
__y_max__ = None
__x_num__ = None
def __init__(self, dictionary_arg):
self.__x_min__ = dictionary_arg['x_min']
self.__x_max__ = dictionary_arg['x_max']
self.__y_min__ = dictionary_arg['y_min']
self.__y_max__ = dictionary_arg['y_max']
self.__x_num__ = dictionary_arg['x_num']
def __get_xp__(self):
num = self.__x_num__ + 1
xp = np.linspace(self.__x_min__, self.__x_max__, num)
return xp
def get_xx(self):
xp = self.__get_xp__()
xc = []
for ii_in in range(len(xp) - 1):
xm = (xp[ii_in] + xp[ii_in + 1]) / 2.0
xc.append(xm)
return xc
if __name__ == '__main__':
mesh_dictionary = {
'x_min': 0,
'x_max': 30,
'y_min': -0.5,
'y_max': 0.5,
'x_num': 100 # optional: 100 or 3000
}
my_mesh = mesh(dictionary_arg=mesh_dictionary)
xx = my_mesh.get_xx()
t_focus = 15.0
x_focus = 13.0
focus_wave_dictionary = {'water_depth': 3.0, 'x_focus': x_focus, 't_focus': t_focus,
'n_focus': 75, 'significant_wave_height': 0.1, 'mean_wave_period': 2.8}
my_focus = focus_wave(dictionary_arg=focus_wave_dictionary)
my_focus.write_paras()
my_ittc = my_focus.ittc_wave_spectral()
wave_height_out = my_ittc.wave_height()
wave_omega_out = my_ittc.wave_omega()
fig, ax = plt.subplots()
tt_max = 30.0 # optional: 30.0 or 50.0
tt_step = 0.3 # optional: 0.3 or 0.1
tt_min = 0.0
line, = ax.plot([], [])
ax.set_xlim([0.0, 30.0])
ax.set_ylim([-0.3, 0.3])
ims = []
for tt in np.arange(tt_min, tt_max, tt_step):
phi = []
for ii_xx in range(len(xx)):
phi_tmp = 0.0
for ii_inner in range(len(wave_height_out)):
my_wpc = my_focus.wave_para_convertor(wave_height_arg=wave_height_out[ii_inner],
wave_omega_arg=wave_omega_out[ii_inner])
wave_surface_tmp = my_wpc.wave_surface(xx_arg=xx[ii_xx], tt_arg=tt)
phi_tmp += wave_surface_tmp
phi.append(phi_tmp)
line, = ax.plot(xx, phi, '-k')
ims.append([line])
print(str(tt / tt_max * 100.0) + '%\n', flush=True)
ani = animation.ArtistAnimation(fig, ims, interval=10, repeat=False)
writer = animation.PillowWriter(fps=30)
ani.save('wave.gif', writer=writer)
plt.cla()
ax.set_xlim([0.0, 30.0])
ax.set_ylim([-0.3, 0.3])
phi = []
for ii_xx in range(len(xx)):
phi_tmp = 0.0
for ii_inner in range(len(wave_height_out)):
my_wpc = my_focus.wave_para_convertor(wave_height_arg=wave_height_out[ii_inner],
wave_omega_arg=wave_omega_out[ii_inner])
wave_surface_tmp = my_wpc.wave_surface(xx_arg=xx[ii_xx], tt_arg=t_focus)
phi_tmp += wave_surface_tmp
phi.append(phi_tmp)
plt.plot(xx, phi, '-k')
plt.savefig('wave.png')
参考文献
畸形波的实验研究和数值模拟_赵西增.pdf
400
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
