[kumar2]

Dear Niels and All Modelers,

I have been trying to run the wave flume tutorial in the turbulent mode, using the modified k-omega model but the runs are blowing up after about 2 waves entering the domian. I am not sure if I am using the correct initial and boundary conditions for k and omega. I am trying to reproduce Niels published
results in International journal for Numerical Methods in Fluids 2012,
section 3.1 ( Exchange of energy between the harmonics)

This is how I calculated k: the horizontal particle velocity is about 0.2 m/s
(period=3.5s; wavelength=6.8m; waveheight=0.084m). Taking the turbulent intensity as 0.001 percent of the horizontal particle velocity and turbulent length scale as 0.001 percent of the wave height, k=3/2*(0.2*0.001/100)^2=6.8e-12 m^2/s^2. Turbulent length scale=(0.001/100)*0.084=8.4e-07m.
Then omega=sqrt(k)/turb.length=3.1. From k and omega, I calculated nut. (For higher values of k, the code blows up earlier.)

//inlet and BC for k
internalField uniform 6.8e-12;
1)inlet --> type fixedValue; value uniform 6.8e-12;
2)bottom--> type zeroGradient;
3)outlet --> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 6.8e-12; value uniform 6.8e-12;

// inlet and BC for omega:
internalField uniform 3.1;
1)inlet--> type fixedValue; value uniform 3.1;
2)bottom--> type zeroGradient;
3)outlet--> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 3.1; value uniform 3.1;

// inlet and BC for nut
internalField uniform 2.2e-12;
1)inlet--> type fixedValue; value uniform 2.2e-12;
2)bottom--> type zeroGradient;
3)outlet--> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 2.2e-12; value uniform 2.2e-12;

The free surface profile, k and omega fields just before the code blows at t= 6.0 seconds from start is attached.

When the code blows up(and this happens quite abruptly), I see a spike in the Uy solution( given below)

Courant Number mean: 0.0151587 max: 0.278172
Interface Courant Number mean: 0.000481256 max: 0.244935
deltaT = 0.00497824
Time = 6.01581
DILUPBiCG: Solving for Uy, Initial residual = 0.00827393, Final residual = 6.31356e+15, No Iterations 1001

and continuity errors: (given below)
time step continuity errors : sum local = 4.99397e+69, global = 2.22466e+68, cumulative = 2.22466e+68

Please let me know what could be the reason for the code blowing up and if my initial conditions and BC are ok for k and omega

Thanks in advance

Beat regards

Kumar

[ngj]

Hi Kumar,

Why would you use a turbulence model for such a case, where the boundary layers will be fully laminar? In our article we only used a turbulence model for the case with wave breaking.

Kind regards

Niels

[kumar2]

Originally Posted by ngj

Hi Kumar,

Why would you use a turbulence model for such a case, where the boundary layers will be fully laminar? In our article we only used a turbulence model for the case with wave breaking.

Kind regards

Niels

Hi Niels,

Thanks so much for your quick reply. I realize now that I misunderstood the line in the publication( page 1077 , para3, first line). ( I though k-omega is switched on in all cases with wall functions also switched on in case 3.3).

Cheers and thanks again

Kumar

[jasonchen]

Hello Niels and Kumar,

I want to plot a free surface profile figure like the one posted by Kumar. Is this function included in the wave2Foam's followup project?

I recognize the figure by Kumar was plotted by matlab. At one time instant, I can read from the corresponding time diretory the alpha1 field for all cells starting from cell #0. And also I can export from paraFoam the alpha1 valude for all selected cells with point coordinates included for each alpha1 value. With these data I have tried but without success to plot a contour of alpha1 in matlab. I know this may be trivial. Which cmds/funcs you use in matlab?

Thanks.

Best

[ngj]

Hi Jason,

For my part I can tell that such functionality is not planned for the follow-up of waves2Foam.

The tricky part is to get the data from OF-format to something, which is easily read by Matlab. Here, I would suggest that you make some small Matlab functions, which most easily could be applying a bit of unix-command internally. See

help unix

in Matlab; this does require that you use Matlab on a Linux/Unix/OSX distribution.

Once you have achieved to get the data into Matlab, the rest is merely book keeping.

Kind regards

Niels

[kumar2]

Originally Posted by jasonchen

Hello Niels and Kumar,

I want to plot a free surface profile figure like the one posted by Kumar. Is this function included in the wave2Foam's followup project?

I recognize the figure by Kumar was plotted by matlab. At one time instant, I can read from the corresponding time diretory the alpha1 field for all cells starting from cell #0. And also I can export from paraFoam the alpha1 valude for all selected cells with point coordinates included for each alpha1 value. With these data I have tried but without success to plot a contour of alpha1 in matlab. I know this may be trivial. Which cmds/funcs you use in matlab?

Thanks.

Best

Hi Jason,

I first use the OF utility writeCellCenters to write the cell centers then I use the following matlab function readOpenFoamWithMatlab.m to read the cell centers as shown way down.

%%%%%%%%%%%%%function readOpenFoamWithMatlab.m %%%
function alpha1=readOpenFoamWithMatlab(str,no_of_headerLine s) % no spaces.

fid=fopen(str,'r');
%no_of_headerLines=20;
for i=1:no_of_headerLines, fgets(fid); end
no_of_cells=fscanf(fid,'%d'); fgets(fid);
data=textscan(fid,'%f');
alpha1=data{1,1}(; % no smiley face, the editor is inputing this !
[r,c]=size(alpha1); if r~=no_of_cells, pause; end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



Read the cell centers into MATLAB as shown below:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
no_of_headerLines=20; % I am pretty sure, but good to check..
ccx=readOpenFoamWithMatlab('ccx',no_of_headerLines );
ccy=readOpenFoamWithMatlab('ccy',no_of_headerLines );
ccz=readOpenFoamWithMatlab('ccz',no_of_headerLines );
XX=unique(ccx); YY=unique(ccy);
X=zeros(length(YY),length(XX)); Y=zeros(length(YY),length(XX));
for i=1:length(YY), X(i,: )=XX; end
for i=1:length(XX), Y(:,i)=YY; end
[r,c]=size(X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now read alpha fields
%%%%
alpha=readOpenFoamWithMatlab('alpha1',no_of_header Lines);
alpha_i=griddata(ccx,ccy,alpha,X,Y);
figure; pcolor(X,Y,alpha1); colorbar; shading flat;


Please let me know if this does not work.

Best regards

Kumar

[jasonchen]

Originally Posted by kumar2

Hi Jason,

I first use the OF utility writeCellCenters to write the cell centers then I use the following matlab function readOpenFoamWithMatlab.m to read the cell centers as shown way down.

%%%%%%%%%%%%%function readOpenFoamWithMatlab.m %%%
function alpha1=readOpenFoamWithMatlab(str,no_of_headerLine s) % no spaces.

fid=fopen(str,'r');
%no_of_headerLines=20;
for i=1:no_of_headerLines, fgets(fid); end
no_of_cells=fscanf(fid,'%d'); fgets(fid);
data=textscan(fid,'%f');
alpha1=data{1,1}(; % no smiley face, the editor is inputing this !
[r,c]=size(alpha1); if r~=no_of_cells, pause; end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



Read the cell centers into MATLAB as shown below:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
no_of_headerLines=20; % I am pretty sure, but good to check..
ccx=readOpenFoamWithMatlab('ccx',no_of_headerLines );
ccy=readOpenFoamWithMatlab('ccy',no_of_headerLines );
ccz=readOpenFoamWithMatlab('ccz',no_of_headerLines );
XX=unique(ccx); YY=unique(ccy);
X=zeros(length(YY),length(XX)); Y=zeros(length(YY),length(XX));
for i=1:length(YY), X(i,: )=XX; end
for i=1:length(XX), Y(:,i)=YY; end
[r,c]=size(X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now read alpha fields
%%%%
alpha=readOpenFoamWithMatlab('alpha1',no_of_header Lines);
alpha_i=griddata(ccx,ccy,alpha,X,Y);
figure; pcolor(X,Y,alpha1); colorbar; shading flat;


Please let me know if this does not work.

Best regards

Kumar

Hello Kumar,

Thanks very much for the matlab script. It works well for the waveFlume tutorial. But I doubt whether it also works for non-rectangular mesh, for example mesh with sloping beaches.

By the way, I suggest you use this func-contour(XX,YY,alpha1,1), to get a clear surface profile, like shown in the link below.

https://www.dropbox.com/s/ii6b0yuuph...ve%20flume.jpg

Thank you.

[kumar2]

Originally Posted by jasonchen

Hello Kumar,

Thanks very much for the matlab script. It works well for the waveFlume tutorial. But I doubt whether it also works for non-rectangular mesh, for example mesh with sloping beaches.

By the way, I suggest you use this func-contour(XX,YY,alpha1,1), to get a clear surface profile, like shown in the link below.

https://www.dropbox.com/s/ii6b0yuuph...ve%20flume.jpg

Thank you.

Hi Jason,

Great it works for you. The image looks good. It should work for the sloping beach too ( I have not yet done the sloping beach). So when the griddata works outside the region where there is no actual data, it will simply put in NaNs.

Cheers

Kumar

[Fanfei]

Originally Posted by kumar2

Dear Niels and All Modelers,

I have been trying to run the wave flume tutorial in the turbulent mode, using the modified k-omega model but the runs are blowing up after about 2 waves entering the domian. I am not sure if I am using the correct initial and boundary conditions for k and omega. I am trying to reproduce Niels published
results in International journal for Numerical Methods in Fluids 2012,
section 3.1 ( Exchange of energy between the harmonics)

This is how I calculated k: the horizontal particle velocity is about 0.2 m/s
(period=3.5s; wavelength=6.8m; waveheight=0.084m). Taking the turbulent intensity as 0.001 percent of the horizontal particle velocity and turbulent length scale as 0.001 percent of the wave height, k=3/2*(0.2*0.001/100)^2=6.8e-12 m^2/s^2. Turbulent length scale=(0.001/100)*0.084=8.4e-07m.
Then omega=sqrt(k)/turb.length=3.1. From k and omega, I calculated nut. (For higher values of k, the code blows up earlier.)

//inlet and BC for k
internalField uniform 6.8e-12;
1)inlet --> type fixedValue; value uniform 6.8e-12;
2)bottom--> type zeroGradient;
3)outlet --> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 6.8e-12; value uniform 6.8e-12;

// inlet and BC for omega:
internalField uniform 3.1;
1)inlet--> type fixedValue; value uniform 3.1;
2)bottom--> type zeroGradient;
3)outlet--> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 3.1; value uniform 3.1;

// inlet and BC for nut
internalField uniform 2.2e-12;
1)inlet--> type fixedValue; value uniform 2.2e-12;
2)bottom--> type zeroGradient;
3)outlet--> type zeroGradient;
4)atmosphere--> type inletOutlet; inletValue uniform 2.2e-12; value uniform 2.2e-12;

The free surface profile, k and omega fields just before the code blows at t= 6.0 seconds from start is attached.

When the code blows up(and this happens quite abruptly), I see a spike in the Uy solution( given below)

Courant Number mean: 0.0151587 max: 0.278172
Interface Courant Number mean: 0.000481256 max: 0.244935
deltaT = 0.00497824
Time = 6.01581
DILUPBiCG: Solving for Uy, Initial residual = 0.00827393, Final residual = 6.31356e+15, No Iterations 1001

and continuity errors: (given below)
time step continuity errors : sum local = 4.99397e+69, global = 2.22466e+68, cumulative = 2.22466e+68

Please let me know what could be the reason for the code blowing up and if my initial conditions and BC are ok for k and omega

Thanks in advance

Beat regards

Kumar

Hi Kumar:
I do the some thing as you did. I used K-w model for wave flume case, while it always crashed after 0.5s simulation. and i try a lot ways to resolve it, but they don't work. I want to know whether you have meet the same problems before, and how did you deal with it. Thanks.
Best regards
Fan Fei