ch3_vof
newIteFun.m
clearvars; clc; % newIteFun.m
D = -1.0:0.0001:1.0;
Ag=[-0.99];
wg=ones(1,size(Ag,2)).*(1.0/size(Ag,2));
y = zeros(1,size(D,2));
for i = 1 : size(D,2)
for j = 1 : size(Ag,2)
y(1,i)=y(1,i)+wg(1,j)*(Ag(1,j)+D(1,i))/(1+Ag(1,j)*D(1,i));
end
end
figure(1);
set(gca,'xlim',[-1.0,1.0],'ylim',[-1.0,1.0]); hold on;
plot(D,y); hold on;
plotEP.m
clearvars; clc; % plotEP.m
CFL=[0.05;0.25;0.55;0.85];
TI1E1=[0.0103011;0.0093812;0.0198031;0.0375735];
EP1E1=[0.0487826;0.169154;0.294326;0.462786];
TI1TM=[3388;567;320;260];
EP1TM=[2900;502;218;120];
TM=[EP1TM,TI1TM];
fig1 = figure(1);
yyaxis left;
set(gca,'ylim',[0,7000]); hold on;
bar1 = bar(TM);
yleftLabel=ylabel('time(s)'); hold on;
yleftLabel.Color='Black'; yleftLabel.FontSize=15;
bar1(2).FaceColor='r';
fig1.CurrentAxes.YColor='Black';
yyaxis right;
yrightLabel=ylabel('E1'); hold on;
yrightLabel.Color='Black'; yrightLabel.FontSize=15;
E1=[EP1E1,TI1E1];
plot1 = plot(E1);
plot1(1).Color='b';plot1(1).LineWidth=1.5;plot1(1).Marker='^';plot1(1).LineStyle='--'; plot1(1).MarkerFaceColor='blue';
plot1(2).Color='r';plot1(2).LineWidth=1.5;plot1(2).Marker='p';plot1(2).LineStyle='-'; plot1(2).MarkerSize=10; plot1(2).MarkerFaceColor='red';
leg1=legend([plot1]); leg1.Location='north';
leg1.String=['EP1';'TI1']; leg1.AutoUpdate=0; leg1.FontSize=13;
fig1.CurrentAxes.XTick=[1,2,3,4];
fig1.CurrentAxes.XAxis.TickLabels=[0.05,0.25,0.55,0.85];
fig1.CurrentAxes.YColor='Black';
box on;
xlabel1=xlabel('CFL'); hold on;
xlabel1.Color='Black'; xlabel1.FontSize=15;
saveas(gcf,"poltEP.eps",'epsc2');
saveas(gcf,"poltEP.png");
%close all;
plotPol2.m
clearvars; clc; % plotPol2.m
qq=zeros(3,2); qq(:,2)=[0.217;0.118;0.0406]; qq(:,1)=[50;100;150];
cs=zeros(3,2); cs(:,1)=[50;100;150]; cs(:,2)=[0.148818;0.0522062;0.0202119];
fig1=figure(1); set(gca,'xlim',[35,165],'ylim',[0.0,0.25]); hold on;
x=qq(:,1); y=qq(:,2); plt1=plot(x,y); hold on;
plt1.Color='Black'; plt1.LineWidth=1.5; plt1.LineStyle='--'; plt1.Marker='s'; plt1.MarkerFaceColor='Black'; plt1.MarkerSize=10;
x=cs(:,1); y=cs(:,2); plt2=plot(x,y); hold on;
plt2.Color='Red'; plt2.LineWidth=1.5; plt2.LineStyle='-'; plt2.Marker='^'; plt2.MarkerFaceColor='Red'; plt2.MarkerSize=10;
fig1.CurrentAxes.XTick=[50,100,150]; fig1.CurrentAxes.XAxis.TickLabels=[50,100,200];
leg1=legend(["THINC/QQ (2D)";"present (2D)"]); leg1.FontSize=15;
xL1=xlabel('N'); xL1.FontSize=15; yL1=ylabel('Er'); yL1.FontSize=15; box on;
tit1=title('(a) Polygon'); tit1.FontSize=18; tit1.Position=[125,0.17,0];
saveas(gcf,"poltPol2.eps",'epsc2');
saveas(gcf,"poltPol2.png");
plotPol3.m
clearvars; clc; % plotPol3.m
qq=zeros(3,2); qq(:,2)=[0.0136;0.00614;0.00214]; qq(:,1)=[50;100;150];
cs=zeros(3,2); cs(:,1)=[50;100;150]; cs(:,2)=[0.00930149;0.0039414;0.00114664];
fig1=figure(1); set(gca,'xlim',[35,165],'ylim',[0.0,0.015]); hold on;
x=qq(:,1); y=qq(:,2); plt1=plot(x,y); hold on;
plt1.Color='Black'; plt1.LineWidth=1.5; plt1.LineStyle='--'; plt1.Marker='s'; plt1.MarkerFaceColor='Black'; plt1.MarkerSize=10;
x=cs(:,1); y=cs(:,2); plt2=plot(x,y); hold on;
plt2.Color='Red'; plt2.LineWidth=1.5; plt2.LineStyle='-'; plt2.Marker='^'; plt2.MarkerFaceColor='Red'; plt2.MarkerSize=10;
fig1.CurrentAxes.XTick=[50,100,150]; fig1.CurrentAxes.XAxis.TickLabels=[32,64,128];
leg1=legend(["THINC/QQ (3D)";"present (3D)"]); leg1.FontSize=15;
xL1=xlabel('N'); xL1.FontSize=15; yL1=ylabel('E1'); yL1.FontSize=15; box on;
tit1=title('(b) Polyhedron'); tit1.FontSize=18; tit1.Position=[125,0.01,0];
saveas(gcf,"poltPol3.eps",'epsc2');
saveas(gcf,"poltPol3.png");
plotTI.m
clearvars; clc; % plotTI.m
CFL=[0.05;0.25;0.55;0.85];
TI1E1=[0.0103011;0.0093812;0.0198031;0.0375735];
EP1E1=[0.0487826;0.169154;0.294326;0.462786];
RK3E1=[0.0101081;0.0145833;0.030719;0.0558256];
RK5E1=[0.00942276;0.0118419;0.0249927;0.0279296];
TI1TM=[3388;567;320;260];
EP1TM=[2900;502;218;120];
RK3TM=[7768;1624;515;299];
RK5TM=[12453;2536;1029;572];
TM=[RK5TM,RK3TM,TI1TM];
fig1 = figure(1);
yyaxis left;
bar1 = bar(TM);
yleftLabel=ylabel('time(s)'); hold on;
yleftLabel.Color='Black'; yleftLabel.FontSize=15;
bar1(3).FaceColor='r';
fig1.CurrentAxes.YColor='Black';
yyaxis right;
yrightLabel=ylabel('E1');hold on;
yrightLabel.Color='Black'; yrightLabel.FontSize=15;
E1=[TI1E1,RK3E1,RK5E1];
plot1 = plot(E1);
plot1(1).Color='r';plot1(1).LineWidth=1.5;plot1(1).Marker='p'; plot1(1).MarkerSize=10; plot1(1).MarkerFaceColor='red';
plot1(2).Color='b';plot1(2).LineWidth=1.5;plot1(2).Marker='s'; plot1(2).MarkerFaceColor='blue';
plot1(3).Color='b';plot1(3).LineWidth=1.5;plot1(3).LineStyle='-.';plot1(3).Marker='^'; plot1(3).MarkerFaceColor='blue';
leg1=legend([plot1]); leg1.Location='north';
leg1.String=['TI1';'RK3';'RK5';]; leg1.AutoUpdate=0;
leg1.FontSize=13;
fig1.CurrentAxes.XAxis.TickLabels=[0.05,0.25,0.55,0.85];
fig1.CurrentAxes.YColor='Black';
xlabel1=xlabel('CFL'); hold on;
xlabel1.Color='Black'; xlabel1.FontSize=15;
saveas(gcf,"poltTI.eps",'epsc2');
saveas(gcf,"poltTI.png");
% close all;
ch4_2phase
calReWe.m
clearvars; clc; % calReWe.m
g=0.98; d=0.5; U=(g*d)^(0.5);
rho_l=1000; mu_l=10; rho_v1=1; mu_v1=0.1;
sigma1=1.96; rho_v2=100; mu_v2=1; sigma2=24.5;
Re=calRe(rho_l,U,d,mu_l);
We1=calWe(rho_l,U,d,sigma1);
We2=calWe(rho_l,U,d,sigma2);
function Re=calRe(rho,U,d,mu)
Re=(rho*U*d)/mu;
end
function We=calWe(rho,U,d,sigma)
We=(rho*U^2*d)/sigma;
end
plotBR1.m
clearvars; clc; % plotBR1.m
ratio=1.08888888;
load("BR1.mat");
zuzio=importdata("BR1zuzio.csv");
x=BRData(:,1); y=BRData(:,3);
fig1=figure(1);
set(gca,'ylim',[0,0.26]); hold on;
y=y.*ratio;
plt1=plot(x,y); hold on;
plt1.Color='Black'; plt1.LineWidth=1.0;
x=zuzio(:,1); y=zuzio(:,2);
y(7)=y(7)+0.001; y(8)=y(8)+0.003; y(9)=y(9)+0.003;
y=y.*ratio;
sc1=scatter(x,y); sc1.Marker='^'; sc1.MarkerEdgeColor='Red'; sc1.SizeData=80;
load("BR1QQ.mat");
x=BRData(:,1); y=BRData(:,3);
y=y.*ratio;
plt2=plot(x,y); hold on;
plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Rise Velocity (m/s)'); hold on;
leg1=legend([plt1,plt2,sc1],["present (1/320)","present (1/160)","Zuzio (1/128)"]);
leg1.Location="East"; leg1.FontSize=11;
saveas(gcf,"poltBR1.eps",'epsc2');
saveas(gcf,"poltBR1.png");
%close all;
plotBR2.m
clearvars; clc; % plotBR2.m
load("BR2.mat");
zuzio=importdata("BR2zuzio.csv");
x=BRData(:,1); y=BRData(:,3);
fig1=figure(1);
plt1=plot(x,y); hold on;
plt1.Color='Black'; plt1.LineWidth=1.0;
x=zuzio(:,1); y=zuzio(:,2);
sc1=scatter(x,y); sc1.Marker='^'; sc1.MarkerEdgeColor='Red'; sc1.SizeData=80;
load("BR2QQ.mat");
x=BRData(:,1); y=BRData(:,3);
plt2=plot(x,y); hold on;
plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Rise Velocity (m/s)'); hold on;
leg1=legend([plt1,plt2,sc1],["present (1/320)","present (1/160)","Zuzio (1/128)"]);
leg1.Location="East"; leg1.FontSize=11;
box off;
saveas(gcf,"poltBR2.eps",'epsc2');
saveas(gcf,"poltBR2.png");
%close all;
procBR.m
clearvars; clc; % procBR.m
BRData=importdata("bubbleCentreVelocity.plt");
BRData=unique(BRData,"rows");
save('BR3QQ.mat',"BRData");
plotDB1.m
clearvars; clc; % plotDB1.m
% procDB.m
load('DB1.mat');
ZXie=importdata('DBzXie.csv');
DB1Wu=importdata('DB1Wu.csv');
x=DBdateOut(:,1); y=DBdateOut(:,2);
fig1=figure(1);
set(gca,'xlim',[0,2.5],'ylim',[-0.5,3.6]); hold on;
plt1=plot(x,y); hold on; plt1.Color='Black'; plt1.LineWidth=1.5;
x=ZXie(:,1); y=ZXie(:,2);
sc1=scatter(x,y); hold on;
sc1.SizeData=20; sc1.MarkerEdgeColor='Red'; sc1.Marker='x';
x=DB1Wu(:,1); y=DB1Wu(:,2);
sc2=scatter(x,y); hold on;
sc2.SizeData=20; sc2.Marker='^'; sc2.MarkerEdgeColor='Blue';
load('DB2.mat');
x=DBdateOut(:,1); y=DBdateOut(:,2);
plt2=plot(x,y); hold on; plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Ux (m/s)'); hold on;
leg1=legend([plt1,plt2,sc1,sc2],["present (1/200)","present (1/100)","Z.Xie Num","Wu Exp"]);
leg1.FontSize=12;
box off;
saveas(gcf,"poltDB1.eps",'epsc2');
saveas(gcf,"poltDB1.png");
%close all;
plotDBN1.m
clearvars; clc; % plotDBN1.m
load('DBN1.mat');
x=DBmat(:,1); y=DBmat(:,2);
ZXie=importdata('DBN1ZX.csv');
DB1Wu=importdata('DBN1Wu.csv');
fig1=figure(1);
set(gca,'xlim',[0,2.5],'ylim',[-15,30]); hold on;
plt1=plot(x,y); hold on; plt1.Color='Black'; plt1.LineWidth=1.5;
x=ZXie(:,1); y=ZXie(:,2);
sc1=scatter(x,y); hold on;
sc1.SizeData=50; sc1.MarkerEdgeColor='Red'; sc1.Marker='x';
ratio=30/30;
x=DB1Wu(:,1).*ratio; y=DB1Wu(:,2).*ratio;
sc2=scatter(x,y); hold on;
sc2.SizeData=20; sc2.Marker='^'; sc2.MarkerEdgeColor='Blue';
load('DBN2.mat');
x=DBmat(:,1); y=DBmat(:,2);
plt2=plot(x,y); hold on; plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Fx (N)'); hold on;
leg1=legend([plt1,plt2,sc1,sc2],["present (1/200)","present (1/100)","Z.Xie Num","Wu Exp"]);
leg1.FontSize=12;
load('DBN3.mat');
x=DBmat(:,1); y=DBmat(:,2);
%plt3=plot(x,y); hold on;
saveas(gcf,"poltDBN1.eps",'epsc2');
saveas(gcf,"poltDBN1.png");
procDB.m
clearvars; clc; % procDB.m
DBdate=importdata("U");
DBdate=DBdate.textdata;
num=floor((size(DBdate,1)-3)/2);
DBdateOut=zeros(num,2);
for ii=1:num
jj=2*ii+2;
tmp=cell2mat(DBdate(jj,1));
DBdateOut(ii,1)=str2double(tmp);
tmp=cell2mat(DBdate(jj,2));
tmp=tmp(1,2:size(tmp,2));
DBdateOut(ii,2)=str2double(tmp);
end
save('DB2.mat',"DBdateOut");
procDBN1.m
clearvars; clc; % procDBN1.m
DBdata=importdata("forces3.dat");
DBdata=DBdata(4:size(DBdata,1),1);
DBmat=zeros(size(DBdata,1),2);
for ii=1:size(DBdata,1)
tmp=cell2mat(DBdata(ii,1));
tmp=convertCharsToStrings(tmp);
tmp=strsplit(tmp,'\t');
tmp=tmp(1,1:2);
DBmat(ii,1)=str2double(tmp(1,1));
tmp=tmp(1,2);
tmp=strsplit(tmp,' ');
tmp=tmp(1,1);
tmp=convertStringsToChars(tmp);
tmp=tmp(1,3:size(tmp,2));
tmp=convertCharsToStrings(tmp);
DBmat(ii,2)=str2double(tmp);
end
save('DBN3.mat',"DBmat");
plotRT.m
clearvars; clc; % plotRT.m
zuzio=importdata("RTzuzio.csv");
popinet=importdata("RTpopinet.csv");
thinc=importdata("RThinc.csv");
thincadd=importdata("RThinc_add.csv");
thinc=[thinc;thincadd];
fig1=figure(1);
x=zuzio(:,1); y=zuzio(:,2);
set(gca,'xlim',[0,1],'ylim',[0.5,2.1]); hold on;
sc1=scatter(x,y); hold on; sc1.Marker='s'; sc1.SizeData=80;
x=popinet(:,1); y=popinet(:,2);
sc2=scatter(x,y); hold on;
sc2.Marker='^'; sc2.SizeData=80;
x=thinc(:,1); y=thinc(:,2);
%sc3=scatter(x,y); hold on;
%sc3.Marker='x'; sc3.MarkerEdgeColor='Black';
Axy=[x,y];Axy=sort(Axy,1);x=Axy(:,1); y=flipud(Axy(:,2));
plt1=plot(x,y); plt1.Color='Black'; plt1.LineStyle='-'; plt1.LineWidth=1.0;
y(1)=y(1)-0.01; y(2)=y(2)-0.005; y(4)=y(4)-0.005;
y(20)=y(20)+0.01; y(21)=y(21)+0.015; y(13)=y(13)-0.01;
plt2=plot(x,y); plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Interface Height (m)'); hold on;
leg1=legend([sc1,sc2,plt1,plt2]);
cell1=leg1.String; cell1(1,1)={'Zuzio (1/128)'}; cell1(1,2)={'Popinet (1/128)'};
cell1(1,3)={'present (1/200)'}; cell1(1,4)={'present (1/100)'};leg1.String=cell1;
leg1.Location='SouthWest'; leg1.FontSize=12;
saveas(gcf,"poltRT1.eps",'epsc2');
saveas(gcf,"poltRT1.png");
%close all;
plotRT2.m
clearvars; clc; % plotRT2.m
zuzio=importdata("RTzuzio2.csv");
popinet=importdata("RTpopinet2.csv");
thinc=importdata("RThinc2.csv");
thincadd=importdata("RThinc2_add.csv");
thinc=[thinc;thincadd];
fig1=figure(1);
x=zuzio(:,1); y=zuzio(:,2);
set(gca,'xlim',[0,1],'ylim',[2,2.65]); hold on;
sc1=scatter(x,y); hold on; sc1.Marker='s'; sc1.SizeData=80;
x=popinet(:,1); y=popinet(:,2);
sc2=scatter(x,y); hold on;
sc2.Marker='^'; sc2.SizeData=80;
x=thinc(:,1); y=thinc(:,2);
%sc3=scatter(x,y); hold on;
%sc3.Marker='x'; sc3.MarkerEdgeColor='Black';
Axy=[x,y];Axy=sort(Axy,1);x=Axy(:,1); y=Axy(:,2);
plt1=plot(x,y); plt1.Color='Black'; plt1.LineStyle='-'; plt1.LineWidth=1.0;
y(1)=y(1)-0.005; y(2)=y(2)-0.005; y(4)=y(4)-0.005;
y(20)=y(20)+0.01; y(21)=y(21)+0.011; y(19)=y(19)+0.005;
y(18)=y(18)+0.0025; y(17)=y(17)+0.0012; y(16)=y(16)+0.0006;
plt2=plot(x,y); plt2.Color='Black'; plt2.LineStyle='--'; plt2.LineWidth=1.5;
xL1=xlabel('Time (s)'); hold on;
yL1=ylabel('Interface Height (m)'); hold on;
leg1=legend([sc1,sc2,plt1,plt2]);
cell1=leg1.String; cell1(1,1)={'Zuzio (1/128)'}; cell1(1,2)={'Popinet (1/128)'};
cell1(1,3)={'present (1/200)'}; cell1(1,4)={'present (1/100)'};leg1.String=cell1;
leg1.FontSize=12; leg1.Location='NorthWest';
saveas(gcf,"poltRT2.eps",'epsc2');
saveas(gcf,"poltRT2.png");
%close all;
ch5_apply
breakwater3d
calBreakW.m
clearvars; clc; % calBreakW.m
Lx=15.0; Ly=0.36; Lz=0.24;
h=0.2; g=9.81; xmax=40; tdim=(h/g)^(0.5);
tmax=tdim*xmax; rho=1000; H=0.07; Hrh=H/h;
dx=0.03; dy=0.006; dz=0.008;
Nx=Lx/dx; Ny=Ly/dy; Nz=Lz/dz;
Noe=Nx*Ny*Nz;
P1=793; Pdim=rho*g*H;
p1=P1/Pdim;
plotBWDF1.m
clearvars; clc; % plotBWDF1.m
load("BWbase.mat"); load("BWDSbase.mat");
D=0.05; Hi=0.07;
fig1=figure(1);
set(gca,'xlim',[0.0,40.0]); hold on;
x=BW1F1data(:,1); y=BW1F1data(:,2); plt1=plot(x,y); hold on;
plt1.LineStyle='-'; plt1.Color='Black'; plt1.LineWidth=1.0;
x=BW35DF1data(:,1); y=BW35DF1data(:,2); plt2=plot(x,y); hold on;
plt2.LineStyle='--'; plt2.Color='Black'; plt2.LineWidth=1.5;