1.Fluent中的UDF造波代码
// 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;
}
/**********************************/
2.利用Python实现波浪参数转化,并制作GIF动图
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, 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
}
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
tt_step = 0.3
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])
ani = animation.ArtistAnimation(fig, ims, interval=10, repeat=False)
writer = animation.PillowWriter(fps=30)
ani.save('wave.gif', writer=writer)
参考文献