A list,

Matlab based 3d wave model simulation: floating structure design in the field of ocean engineering needs to accurately calculate wave load to ensure stability and safety. In this paper, longuet-Higgins long peak wave model and 3d random wave model of irregular short peak wave are simulated by MATLAB. The results show that using wave spectrum to simulate 3d random wave can get more accurate wave surface and wave height, and the wave load on the structure can be analyzed and calculated according to the potential flow theory of fluid, which provides a necessary basis for verifying the structural strength of the structure.

Ii. Source code

clear all;
close all;

nhFig = 0; % Figure Number;
SeaRegLx = 40e+3; % Sea region length, unit: m
SeaRegLy = 40e+3; % Sea region length, unit: M % According to the sampling theorem in the frequency domain, after determining the sampling period, that is, the footprint width of the beam, % can realize sampling in the frequency domain by using the number of sampling points in the time domain sampling (the sampling points can be calculated according to the width of the spatial spectrum). Therefore, its spectral width generally refers to the incomplete spectral width % that includes most of the energy (or the main energy part). Therefore, the spectral width corresponding to the main energy part is calculated according to the spectral function in advance. XSampleStep =50; % unit: m
ySampleStep = 50;
xNs = round( SeaRegLx / xSampleStep );
yNs = round( SeaRegLy / ySampleStep );

x = linspace( -SeaRegLx / 2, SeaRegLx / 2, xNs );
y = linspace( -SeaRegLy/ 2, SeaRegLy / 2, yNs ); % wave spectrum simulation g =9.8; % gravity acceleration
% swell wave spectrum parameter
SwellWaveLength = 1000; % wave length KswellWavePeak =2 * pi / SwellWaveLength;
SwellAngled=0;
SwellAngle = SwellAngled/ 180* pi; % Angle between surge and observed direction KxswellWavePeak = KswellWavePeak *cos( SwellAngle );
KyswellWavePeak = KswellWavePeak * sin( SwellAngle );
SigmaKx = 2.5 e-3; % wave width SigmaKy =2.5 e-3;
SigmaHSwell = 2; % Wind wave spectrum parameter WindAngle =45 / 180* pi; % Angle between wind direction and observation direction U10 =12; % 10Sea surface wind speed at meters5m/s, 10m/s, 15The width of m/s wind wave spectrum is respectively1.5.0.4.0.15The corresponding maximum spatial sampling interval is2m, 10m, 20m
KwindPeak = g / ( 1.2 * U10 )^2; % sea wave spectrum parameter NxSeaWave = xNs; % frequency domain sampling points, convenient Fourier transform NySeaWave = yNs; KxSeaWave =2 * pi / SeaRegLx * ( -xNs / 2 : 1 : xNs / 2 - 1 );
KySeaWave = 2 * pi / SeaRegLy * ( -yNs / 2 : 1 : yNs / 2 - 1 );
KxSeaWaveTicks = KxSeaWave;
KySeaWaveTicks = KySeaWave;
KxSeaWave = ( KxSeaWave == 0 ) * ( max( KxSeaWave ) * 1e-16 ) + KxSeaWave; % avoid divided by 0
KySeaWave = ( KySeaWave == 0 ) * ( max( KySeaWave ) * 1e-16 ) + KySeaWave; % avoid divided by 0

% swell wave spectrum
SpectrumSwell = zeros( NySeaWave, NxSeaWave );
Temp1 = ones( NySeaWave, 1 ) * ( KxSeaWave - KxswellWavePeak ) / SigmaKx;
Temp2 = ( KySeaWave - KyswellWavePeak )' / SigmaKy * ones( 1, NxSeaWave ); SpectrumSwell = SigmaHSwell^ 2/2 / PI/SigmaKx/SigmaKy * exp(-0.5 * (Temp1.^2 + Temp2.^2)); clear Temp1; clear Temp2; figure; colormap(gray(256)); image( KxSeaWave, KySeaWave, 256 - 255 / ( max( max( abs( SpectrumSwell ) ) ) - min( min ( abs( SpectrumSwell ) ) ) ) * ( abs( SpectrumSwell ) - min(  min ( abs( SpectrumSwell ) ) ) ) ); axis('xy');
xlabel( 'Kx: wave number in the X direction');
ylabel( 'KY: wave number in Y direction');
title( 'Surge spectrum');

KxSeaWaveMatrix = ones( NySeaWave, 1 ) * KxSeaWave;
KySeaWaveMatrix = KySeaWave' * ones( 1, NxSeaWave );
KSeaWaveTemp1 = ( sqrt( KxSeaWaveMatrix.^2 + KySeaWaveMatrix.^2)); Fwk1 =exp(-1.22 * ( sqrt( KSeaWaveTemp1 ./ KwindPeak ) - 1). ^2  );
HKKpx1 = 1.24 * ( ( KSeaWaveTemp1 / KwindPeak < 0.31 ) & ( KSeaWaveTemp1 / KwindPeak >= 0)); HKKpx2 =2.61 * ( ( KSeaWaveTemp1 / KwindPeak ).^0.65 ) .* ( ( KSeaWaveTemp1 / KwindPeak < 0.9 ) & ( KSeaWaveTemp1 / KwindPeak >= 0.31)); HKKpx3 =2.28 * ( ( KSeaWaveTemp1 / KwindPeak ).^( 0.65 ) ) .* ( KSeaWaveTemp1 / KwindPeak >= 0.9 );
HKKp1 = HKKpx1 + HKKpx2 + HKKpx3;
Temp1x = 1.62 * 1e-3 * U10 / ( g^0.5 ) ./ ( KSeaWaveTemp1 ).^3.5;
Temp2x = exp( -( KwindPeak ./ KSeaWaveTemp1 ).^2* ().1.7 .^ Fwk1 );
Temp3x = ( HKKp1 .* ( sech( ( HKKp1 .* ( atan( ( KySeaWaveMatrix ./ KxSeaWaveMatrix ) ) - WindAngle ) ) ) ).^2 );
% SpectrumWind1 = ( Temp1x .* Temp2x .* Temp3x )';
fid=fopen('sea_top_wl1000_wh50m.dat'.'w')
fprintf(fid,'% 10.5 f \ n',real(wh));
status=fclose(fid);

fid=fopen('sea_top_wl1000_wh50m.dat'.'r');
sea_dem=fscanf(fid,'%g'[800.800]);
fclose(fid);
%image( x, y, 256 - 255/ ( max( max( real( sea_dem ) ) ) - min( min ( real( sea_dem ) ) ) ) * ( real( sea_dem ) - min( min ( real( sea_dem ) ) ))); % mesh(sea_dem) colorbar % view(0.90)
%colormap(gray)
% This function computes the non-coherent integration improvment
% factor using the empirical formula defind in Eq. (2.48)
fact1 = 1.0 + log10( 1.0 / pfa) / 46.6;
fact2 = 6.79 * (1.0 + 0.235 * pd);
fact3 = 1.0 - 0.14 * log10(np) + 0.0183 * (log10(np))^2;
impr_of_np = fact1 * fact2 * fact3 * log10(np);
return

% SpectrumWind = SpectrumWind1';
SpectrumWind = Temp1x .* Temp2x .* Temp3x;
clear KSeaWaveTemp1 Fwk1 HKKpx1 HKKpx2 HKKpx3 HKKp1 Temp1x Temp2x Temp3x;

Copy the code

3. Operation results

Fourth, note

Version: 2014 a