clear all
close all
clc

a=[0.005 0.005 0.005 0.005 0.005 0.000];
b=[0.005 0.005 0.005 0.005 0.005 0.000];

foil=[1	0.00126
0.987365	0.003184736
0.967647	0.007767849
0.947188	0.014411609
0.926179	0.021342294
0.904775	0.027050448
0.883095	0.031305955
0.861225	0.034685626
0.839225	0.037839926
0.817136	0.041131026
0.794988	0.044631904
0.7728	0.048250888
0.750587	0.051848048
0.728359	0.05530328
0.706124	0.058536669
0.683889	0.061510081
0.661658	0.064215981
0.639437	0.066660754
0.61723	0.068858852
0.595041	0.070824229
0.572875	0.072567976
0.550736	0.074097221
0.52863	0.075413448
0.506562	0.076515475
0.484538	0.077400646
0.462567	0.078064974
0.440656	0.078503844
0.418816	0.0787127
0.397062	0.078689393
0.375412	0.078431267
0.353887	0.077936145
0.332519	0.077203671
0.311348	0.076237145
0.290429	0.075037849
0.269832	0.073614001
0.249649	0.071976052
0.229991	0.070140256
0.210992	0.068131451
0.192796	0.065979514
0.175549	0.06372291
0.159376	0.061403269
0.144366	0.059059732
0.13056	0.056730022
0.11795	0.054441758
0.106485	0.052215645
0.096086	0.050060277
0.086659	0.047980715
0.078106	0.045974699
0.070329	0.044037339
0.06324	0.042160702
0.05676	0.040335856
0.050819	0.038553868
0.045356	0.036803727
0.040318	0.035077221
0.035663	0.033361795
0.031351	0.031649373
0.027354	0.029927231
0.023644	0.028185844
0.020203	0.026411723
0.017015	0.024594265
0.014071	0.022716563
0.011365	0.020766167
0.008901	0.018722882
0.006684	0.01657023
0.004734	0.014289333
0.003076	0.011862976
0.001742	0.009275052
0.000771	0.006528425
0.000194	0.003624466
5.20E-05	0.002083263
0	0
0.000194	-0.00127937
0.000749	-0.003196403
0.001693	-0.005246145
0.00299	-0.007315438
0.004605	-0.009354601
0.006505	-0.011350674
0.008666	-0.013298179
0.01107	-0.015195015
0.013708	-0.017044555
0.016578	-0.018850289
0.019683	-0.020617483
0.023033	-0.022349875
0.02664	-0.024048493
0.030524	-0.025718416
0.034706	-0.02735867
0.039217	-0.028971991
0.044089	-0.03055579
0.049364	-0.032109898
0.055089	-0.033632521
0.061321	-0.035116023
0.068124	-0.036556922
0.07557	-0.037944874
0.083744	-0.039267505
0.092734	-0.040510041
0.102635	-0.041654522
0.113541	-0.042678799
0.125533	-0.043557106
0.138673	-0.044262449
0.152984	-0.044770915
0.168448	-0.045060473
0.184998	-0.045117564
0.202526	-0.044940786
0.220902	-0.044534182
0.239983	-0.043915053
0.259636	-0.043101993
0.279742	-0.042118743
0.300204	-0.040991275
0.320945	-0.039741932
0.341905	-0.038393814
0.363037	-0.036964275
0.384308	-0.035473175
0.405692	-0.033934468
0.427168	-0.03236144
0.448722	-0.030764816
0.470341	-0.029153139
0.492016	-0.027533965
0.513738	-0.025914385
0.535502	-0.024299634
0.557301	-0.022695052
0.57913	-0.021110126
0.600985	-0.019555617
0.62286	-0.018044322
0.644752	-0.016594129
0.666657	-0.015226473
0.68857	-0.013958502
0.710487	-0.012799324
0.732402	-0.011740917
0.754308	-0.010744301
0.776199	-0.009731089
0.798064	-0.008582607
0.819888	-0.007173034
0.841651	-0.005431711
0.863324	-0.003462501
0.884867	-0.001656778
0.906218	-0.000666078
0.927295	-0.001039419
0.947982	-0.002521596
0.968128	-0.003662041
0.987549	-0.002834997
1	-0.00126];

x=(0:0.01:1);
x1=foil(1:71,1);y1=foil(1:71,2);
x2=foil(71:141,1);y2=foil(71:141,2);

% figure
% plot(x1,y1,x2,y2)

z=[0.15 0.3 0.45 0.6 0.75 0.9];
pi=3.1415926;

f1=x1.^0.25.*(1-x1).*exp(-20*x1);
f2=(sin(pi*x1.^(log(0.5)/log(z(2))))).^3;
f3=(sin(pi*x1.^(log(0.5)/log(z(3))))).^3;
f4=(sin(pi*x1.^(log(0.5)/log(z(4))))).^3;
f5=(sin(pi*x1.^(log(0.5)/log(z(5))))).^3;
f6=(sin(pi*x1.^(log(0.5)/log(z(6))))).^3;
yup=y1+a(1)*f1+a(2)*f2+a(3)*f3+a(4)*f4+a(5)*f5+a(6)*f6;

f1=x2.^0.25.*(1-x2).*exp(-20*x2);
f2=(sin(pi*x2.^(log(0.5)/log(z(2))))).^3;
f3=(sin(pi*x2.^(log(0.5)/log(z(3))))).^3;
f4=(sin(pi*x2.^(log(0.5)/log(z(4))))).^3;
f5=(sin(pi*x2.^(log(0.5)/log(z(5))))).^3;
f6=(sin(pi*x2.^(log(0.5)/log(z(6))))).^3;
ydown=y2+b(1)*f1+b(2)*f2+b(3)*f3+b(4)*f4+b(5)*f5+b(6)*f6;

%  figure
%  plot(x1,y1,x2,y2,x1,yup,x2,ydown)

FID=fopen('a.txt','w');
for i=1:71;
fprintf(FID, 'vertex create coordinates %d %d 0\n',x1(i),yup(i));
end 
for i=2:71;
fprintf(FID, 'vertex create coordinates %d %d 0\n',x2(i),ydown(i));
end

fprintf(FID,'edge create nurbs ');
for i=1:61;
fprintf(FID, '"vertex.%d"\\ \n',i);
end 
fprintf(FID,'interpolate\n');

fprintf(FID,'edge create nurbs ');
for i=61:71;
fprintf(FID, '"vertex.%d"\\ \n',i);
end 
fprintf(FID,'interpolate\n');

fprintf(FID,'edge create nurbs ');
for i=71:81;
fprintf(FID, '"vertex.%d"\\ \n',i);
end 
fprintf(FID,'interpolate\n');

fprintf(FID,'edge create nurbs ');
for i=81:141;
fprintf(FID, '"vertex.%d"\\ \n',i);
end 
fprintf(FID,'interpolate\n');

fprintf(FID,'edge create radius %d startangle 90 endangle 270 center "vertex.71" xyplane arc \n',5);
fprintf(FID,'edge create straight "vertex.1" "vertex.141" \n');

fprintf(FID,'edge create translate "vertex.142" "vertex.143" vector %d 0 0 \n',1);
fprintf(FID,'edge create translate "vertex.1" "vertex.141" vector %d 0 0 \n',5);
fprintf(FID,'edge create translate "vertex.144" "vertex.145" vector %d 0 0 \n',5);

fprintf(FID,'edge create straight "vertex.144" "vertex.1" \n');
fprintf(FID,'edge create straight "vertex.141" "vertex.145" \n');
fprintf(FID,'edge create straight "vertex.146" "vertex.147" \n');
fprintf(FID,'edge create straight "vertex.142" "vertex.61" \n');
fprintf(FID,'edge create straight "vertex.81" "vertex.143" \n');
fprintf(FID,'edge create straight "vertex.146" "vertex.148"\n');
fprintf(FID,'edge create straight "vertex.147" "vertex.149"\n');

fprintf(FID,'face create wireframe "edge.16" "edge.5" "edge.17" "edge.2" "edge.3" real \n');
fprintf(FID,'face create wireframe "edge.17" "edge.4" "edge.14" "edge.8" real \n');
fprintf(FID,'face create wireframe "edge.14" "edge.10" "edge.19" "edge.12" real\n');
fprintf(FID,'face create wireframe "edge.9" "edge.10" "edge.6" "edge.15" real\n');
fprintf(FID,'face create wireframe "edge.16" "edge.7" "edge.13" "edge.1" real\n');
fprintf(FID,'face create wireframe "edge.11" "edge.13" "edge.18" "edge.9" real \n');

fprintf(FID,'edge mesh "edge.16" "edge.13" successive ratio1 0.98039216 intervals 150 \n');
fprintf(FID,'edge mesh "edge.17" "edge.14" "edge.19" "edge.18" successive ratio1 1.02 intervals 150 \n');
fprintf(FID,'edge mesh "edge.9" "edge.10" "edge.11" "edge.12" successive ratio1 1.04 intervals 100 \n');

fprintf(FID,'face mesh "face.5" map size %d\n', 0.002);
fprintf(FID,'face mesh "face.2" map size %d\n', 0.002);
fprintf(FID,'face mesh "face.3" map size %d\n', 0.01);
fprintf(FID,'face mesh "face.6" map size %d\n', 0.01);
fprintf(FID,'face mesh "face.4" map size %d\n', 0.01);
fprintf(FID,'face mesh "face.1" map size %d\n', 0.25);

fprintf(FID,'physics create "outlet" btype "PRESSURE_OUTLET" edge "edge.18" "edge.19" "edge.15"\n');
fprintf(FID,'physics create "far" btype "PRESSURE_FAR_FIELD" edge "edge.11" "edge.7" "edge.8" "edge.12" "edge.5"\n');
fprintf(FID,'physics create "wall1" btype "WALL" edge "edge.1" "edge.2" \n');
fprintf(FID,'physics create "wall2" btype "WALL" edge "edge.4" "edge.3"\n');
fprintf(FID,'physics create "wall3" btype "WALL" edge "edge.6" \n');
fprintf(FID,'physics create "flow" ctype "FLUID" face "face.1" "face.2" "face.3" "face.4" "face.5" "face.6"\n');
fprintf(FID,'export fluent5 "blade2.msh" nozval \n');

fclose(FID);

%fluent journel b.txt
FID=fopen('b.txt','w');
%-------------------------read model
fprintf(FID,'file/read-case blade2.msh \n');
fprintf(FID,'grid/reorder/reorder-domain \n');

%-----------materials
fprintf(FID,'define/materials/change-create \n');
fprintf(FID,'air\n');
% material-name> air
fprintf(FID,'\n');
% material name [air] 
fprintf(FID,'y\n');
% change Density? [no] y
fprintf(FID,'ideal-gas\n');
% new method [constant] ideal-gas
fprintf(FID,'n \n');
% change Cp (Specific Heat)? [no] 
fprintf(FID,'n\n');
% change Thermal Conductivity? [no] 
fprintf(FID,'n\n');
% change Viscosity? [no] 
fprintf(FID,'n\n');
% change Molecular Weight? [no] 
fprintf(FID,'n\n');
% change L-J Characteristic Length? [no] 
fprintf(FID,'n\n');
% change L-J Energy Parameter? [no] 
fprintf(FID,'n\n');
% change Thermal Expansion Coefficient? [no] 
fprintf(FID,'n\n');
% change Degrees of Freedom? [no] 
fprintf(FID,'n\n');
% change Speed of Sound? [no] 
fprintf(FID,'define/models/viscous \n');
fprintf(FID,'spalart-allmaras \n');
fprintf(FID,'y \n');
% Enable Spalart-Allmaras turbulence model? [no] y
fprintf(FID,'q \n');

%------------------inlet condition
fprintf(FID,'define/boundary-conditions \n');
fprintf(FID,' pressure-far-field \n');
fprintf(FID,'far \n');
% zone id/name [far] 
fprintf(FID,'n \n');
% Use Profile for Gauge Pressure? [no] 
fprintf(FID,'101325 \n');
% Gauge Pressure (pascal) [0] 101325
fprintf(FID,'n \n');
% Use Profile for Mach Number? [no] 
fprintf(FID,'0.63 \n');
% Mach Number [0.6] 0.63
fprintf(FID,'n \n');
% Use Profile for Temperature? [no] 
fprintf(FID,'300 \n');
% Temperature (k) [300] 
fprintf(FID,'n \n');
% Use Profile for X-Component of Flow Direction? [no] 
fprintf(FID,'1 \n');
% X-Component of Flow Direction [1] 
fprintf(FID,'n \n');
% Use Profile for Y-Component of Flow Direction? [no] 
fprintf(FID,'0 \n');
% Y-Component of Flow Direction [0] 
fprintf(FID,'y \n');
% Turbulent Specification Method: Modified Turbulent Viscosity [yes]
fprintf(FID,'n \n');
% Use Profile for Modified Turbulent Viscosity? [no]
fprintf(FID,'0.001 \n');
% Modified Turbulent Viscosity (m2/s) [0.001]
fprintf(FID,'q \n');

%---------------outlet condition
fprintf(FID,'define/boundary-conditions \n');
fprintf(FID,'pressure-outlet \n');
fprintf(FID,'outlet \n');
% zone id/name [outlet] outlet
fprintf(FID,'n \n');
% Use Profile for Gauge Pressure? [no] 
fprintf(FID,'101325 \n');
% Gauge Pressure (pascal) [0] 101325
fprintf(FID,'n \n');
% Use Profile for Backflow Total Temperature? [no] 
fprintf(FID,'300 \n');
% Backflow Total Temperature (k) [300] 
fprintf(FID,'n \n');
% Backflow Direction Specification Method: Direction Vector [no] 
fprintf(FID,'y \n');
% Backflow Direction Specification Method: Normal to Boundary [yes] 
fprintf(FID,'y \n');
% Turbulent Specification Method: Modified Turbulent Viscosity [yes] 
fprintf(FID,'n \n');
% Use Profile for Backflow Modified Turbulent Viscosity? [no] 
fprintf(FID,'0.001 \n');
% Backflow Modified Turbulent Viscosity (m2/s) [0.001] 
fprintf(FID,'n \n');
% Specify targeted mass flow rate [no] 
fprintf(FID,'q \n');

%%%-----------------operating-condition
fprintf(FID,'define/operating-conditions \n');
fprintf(FID,'operating-pressure \n');
fprintf(FID,'0  \n');
%operating pressure (pascal) [101325] 0
fprintf(FID,'q \n');

fprintf(FID,'solve/initialize/set-defaults/x-velocity \n');
fprintf(FID,'150 \n');
fprintf(FID,'q \n');
fprintf(FID,'solve/initialize/set-defaults/pressure \n');
fprintf(FID,'101325 \n');
fprintf(FID,'q \n');
fprintf(FID,'q \n');
fprintf(FID,'q \n');
fprintf(FID,'q \n');

fprintf(FID,'solve/initialize \n');
fprintf(FID,'initialize-flow \n');
fprintf(FID,'q \n');

fprintf(FID,'solve/monitors/residual/convergence-criteria \n');
fprintf(FID,'0.00001 \n');
%continuity residual convergence criterion [0.001] 0.00001
fprintf(FID,'0.001 \n');
%x-velocity residual convergence criterion [0.001] 0.001
fprintf(FID,'0.001 \n');
%y-velocity residual convergence criterion [0.001] 0.001
fprintf(FID,'1e-6 \n');
%energy residual convergence criterion [1e-06] 1e-6
fprintf(FID,'1e-3 \n');
%nut residual convergence criterion [0.001] 0.001
fprintf(FID,'q \n');

fprintf(FID,'report/reference-values \n'); 
fprintf(FID,'pressure \n'); 
fprintf(FID,'101325 \n');
fprintf(FID,'velocity \n'); 
fprintf(FID,'218 \n'); 
fprintf(FID,'q \n');

fprintf(FID,'solve \n');
fprintf(FID,'iterate %i \n', 10000);
fprintf(FID,'q \n');

fclose(FID);

status = system('del blade*.*');
status = system('C:\Fluent.Inc\ntbin\ntx86\gambit.exe -r2.4.6 -id blade2.dbs -inp a.txt');
status = system('C:\Fluent.Inc\ntbin\win64\fluent.exe 2ddp -t4 -i b.txt');
