/************************************************
 *
 * Title	: Pintle Injctor Geometry 1.5tonf
 * Date		: 2019.02.01
 * Writer	: Dokeun Hwang(dokeun@kari.re.kr)
 * Version	: 1.0
 * Geometric Unit : mm, sec, kg
 ************************************************/

lc = 1e0; //prescribed mesh element size of the points
m = 1e3; //millimetre to meter
bar2Pascal = 1e5; //bar to pascal

/************************************************/
/* Pintle Injector Design Parameter Calculation */
/************************************************/

/*Physical properties defined*/
rho_annulus	= 846.5;		// density of annulus fluid, kg/m3
rho_pintle	= 1263.4;		// density of pintle fluid, kg/m3
mfr_annulus	= 1.1;			// mass flow rate of annulus gap, kg/s
mfr_pintle	= 3.08;			// mass flow rate of pintle orifice, kg/s
mfrPalloc	= 85;			// mass flow rate allocation on primary slots, %
cd_annulus	= 0.7;			// flow coefficient of annulus gap
cd_pintle	= 0.65;			// flow coefficient of annulus gap
delP_pintle	= 10;			// pressure drop through pintle slot, bar
delP_annulus= 10;			// pressure drop through annulus gap, bar

/*Geometric parameters defined*/
Dc			= 78.67;							// Chamber diameter, mm
DcDp		= 4.4;								// Chamber/Pintle diameter ratio
lsDP		= 0.559284116331;			// Pintle length/diameter ratio
n			= 36;								    // Number of primary, Secondary orifice
AR1			= 8;								  // Aspect ratio of primary orifice
AR2			= 3;							   	// Aspect ratio of secondary orifice
d12			= 1;  								// Distance from Primary to Secondary slot, mm
d23			= 0.2;  								// Distance from Secondary slot to Pintle Cylinder end, mm
Eccent		= 0.866;						// Eccentricity of elliptic Pintle Tip, E=0.866 -> A/B~2

/*Geometric parameters derived*/
Dp			= Dc/DcDp;																				// Pintle Diameter, mm
ls			= Dp*lsDP;																				// Skip distance, mm
ellipseA	= Dp/2.0;							// vertex for ellipse shape at Pintle Tip
ellipseB	= ellipseA*Sqrt((1-Eccent));		// co-vertex for ellipse shape at Pintle Tip
ellipseC	= Sqrt(ellipseA^2 - ellipseB^2);	// focus for ellipse shape at Pintle Tip

AreaPrim	= (mfrPalloc/100)*mfr_pintle/(n*cd_pintle*Sqrt(2.0*rho_pintle*delP_pintle*bar2Pascal))*m^2;		// Primary pintle slot area,mm^2
WPrim		= Sqrt(AreaPrim/AR1);																// Primary pintle slot width, mm
HPrim		= AreaPrim/WPrim;																	// Primary pintle slot height, mm
WPrimg		= Pi*Dp/n-WPrim;																		// Primary pintle slot vent, mm

AreaSec		= AreaPrim*(1.0-mfrPalloc/100);										// Secondary pintle slot area, mm^2
WSec		= Sqrt(AreaSec/AR2);																// Secondary pintle slot width, mm
HSec		= AreaSec/WSec;																	// Secondary pintle slot height, mm
WSecg		= Pi*Dp/n-WSec;																			// Primary pintle slot vent, mm

gap			= (-Dp/m+Sqrt((Dp/m)^2+4*(mfr_annulus/(Pi*cd_annulus*Sqrt(2*rho_annulus*delP_annulus*bar2Pascal)))))/2*m;						// annulus gap derived from Cd & target mfr, mm

BF 			= n*(WPrim+WSec)/(Pi*Dp);																// Blockage Factor

Printf("Dc = %g mm, Dp = %g mm, ls = %g mm",Dc, Dp, ls);
Printf("WPrim = %g mm, WSec = %g mm, gap = %g mm, BF = %g",WPrim,WSec,gap,BF);

/*Calculation of Extrude angle*/
angle =2*Pi/n;																	// Extrude(Sweep) angle decision
angle_PrimHalf = Asin(WPrim/Dp);
angle_SecHalf  = Asin(WSec/Dp);

/*Physical parameters derived*/
mfr_Prim	= mfr_pintle*(mfrPalloc/100)/n;															// mass flow rate per one Primary slot, kg/s
mfr_Sec		= mfr_pintle*(1-mfrPalloc/100)/n;														// mass flow rate per one Secondary slot, kg/s
UPrim		= mfr_Prim/(rho_pintle*AreaPrim);														// velocity of jets from Primary slot, m/s
USec		= mfr_Sec/(rho_pintle*AreaSec);															// velocity of jets from Secondary slot = UPrim, m/s


If (d23 > 0)
	/************************************************/
	/*               Geometry Drawing               */
	/************************************************/

	/*Geometric Points & Line drawing*/
	p1 = newp; Point(p1)  	 = {0,									Dp/2,	0,lc};
	p2 = newp; Point(p2)  	 = {ls,									Dp/2,	0,lc};
	p3 = newp; Point(p3)  	 = {ls+HPrim,							Dp/2,	0,lc};
	p4 = newp; Point(p4)  	 = {ls+HPrim+d12,						Dp/2,	0,lc};
	p5 = newp; Point(p5)   	 = {ls+HPrim+d12+HSec,					Dp/2,	0,lc};
	p6 = newp; Point(p6)  	 = {ls+HPrim+d12+HSec+d23,				Dp/2,	0,lc};
	p7 = newp; Point(p7)   	 = {ls+HPrim+d12+HSec+d23,				0, 		0,lc};
	p8 = newp; Point(p8)  	 = {ls+HPrim+d12+HSec+d23+ellipseB,		0, 		0,lc};
	p9 = newp; Point(p9)  	 = {ls+HPrim+d12+HSec+d23+ellipseB+gap,	0, 		0,lc}; 
	p10 = newp; Point(p10)   = {ls+HPrim+d12+HSec+d23+Dc/2,	0, 0,lc}; 
		
	l1=newl; Line(l1) = {p1, p2}; l2=newl; Line(l2) = {p2, p3};
	l3=newl; Line(l3) = {p3, p4}; l4=newl; Line(l4) = {p4, p5};
	l5=newl; Line(l5) = {p5, p6};

	//Extrude {0, gap, 	0} { Curve{l1, l2, l3, l4, l5}; Layers{5}; Recombine; }
	//Extrude {0, Dc/2-Dp/2-gap, 0} { Curve{l1+5, l2+8, l3+11, l4+14, l5+17}; Layers{5}; Recombine; }
	Extrude {0, gap, 	0} { Curve{l1, l2, l3, l4, l5}; Recombine;}
	Extrude {0, Dc/2-Dp/2-gap, 0} { Curve{l1+5, l2+8, l3+11, l4+14, l5+17}; Recombine;}

	l5 = newl; Ellipse(l5) 	= {p6, p7, p6, p8};
	l6 = newl; Ellipse(l6) 	= {p6+14, p7, p6+14, p9};
	l7 = newl; Line(l7)    	= {p8, p9};
	l8 = newl; Line(l8)    	= {p9, p10};
	l9 = newl; Circle(l9)	= {p10, p7, p6+24};

	Curve Loop(1) = {24, 47, -48, -46}; 	s1 = news; Surface(s1) = {1};
	Curve Loop(2) = {44, -50, -49, -47};	s2 = news; Surface(s2) = {2};
	
	//Extrude {{1, 0, 0}, {0, 0, 0}, angle_PrimHalf} { Surface{:}; Layers{5}; Recombine;}
    //Extrude {{1, 0, 0}, {0, 0, 0}, angle/2-(angle_PrimHalf+angle_SecHalf)} 	{Surface{74, 96, 118,140,162,184,206,228,250,272,289,306}; Layers{5}; Recombine;}
	//Extrude {{1, 0, 0}, {0, 0, 0}, angle_SecHalf} 							{Surface{328,350,372,394,416,438,460,482,504,526,543,560}; Layers{5}; Recombine;}
	Extrude {{1, 0, 0}, {0, 0, 0}, angle_PrimHalf} { Surface{:}; Recombine;}
    Extrude {{1, 0, 0}, {0, 0, 0}, angle/2-(angle_PrimHalf+angle_SecHalf)} 	{Surface{74, 96, 118,140,162,184,206,228,250,272,289,306}; Recombine;}
	Extrude {{1, 0, 0}, {0, 0, 0}, angle_SecHalf} 							{Surface{328,350,372,394,416,438,460,482,504,526,543,560}; Recombine;}
	//angle_PrimHalf
	//angle/2-(angle_PrimHalf+angle_SecHalf)
	//2*angle_SecHalf
	//angle/2-(angle_PrimHalf+angle_SecHalf)
	//angle_PrimHalf
	
	/************************************************/
	/*        Physical Boundary Definition          */
	/************************************************/

	Physical Surface("Wall_Pintle") = {61, 315, 569, 337, 591, 105, 359, 613, 381, 127, 149, 403, 657, 288, 542, 796};
	Physical Surface("Inlet_PrimSlot") = {83};
	Physical Surface("Inlet_SecSlot") = {635};
	Physical Surface("Inlet_Annul") = {73, 327, 581};
	Physical Surface("Wall_Head") = {183, 437, 691};
	Physical Surface("Farfield") = {687, 433, 179, 201, 455, 709, 223, 477, 731, 753, 499, 245, 267, 521, 775, 301, 555, 809};
	Physical Surface("Periodic_L") = {52,  29,   33,   37,  41, 45,  51,  25,  21,  17,  13,   9};
	Physical Surface("Periodic_R") = {814, 692, 714, 736, 758, 780, 797, 670, 648, 626, 604, 582};
	
	v1= newv; Physical Volume("Inner_volume") = Volume{:};
	
	
	/************************************************/
	/*         Grid point calculation with y+       */
	/************************************************/
	// Calulation will be performed later version
	G_WPrim = 20; G_WPS = 20; G_WSecHalf  = 10; 
	G_HPrim = 30; G_HSec  = 10;
	G_Hls   = 30; G_Hd12  = 10; G_Hd23  = 10; 
	G_Rgap = 20; G_RHead = 80; G_RTip  = 160;
	//G_WRSec = G_WPrim + G_WSecHalf - 2 + 5; //other equation needs for automatic calculation
	//G_HP2CE = G_Hd12 + G_HSec - 1;
	//G_HCY = G_Hls + G_HPrim + G_HP2CE - 2;
	//G_Wgap  = 2*G_WRSec + G_WSec - 2; G_Rgap  = 10;
	
	
	//Transfinite Curves
	//Vertical line of ls
	Transfinite Curve {1, 54, 308, 562}		= G_Hls Using Bump 0.01;
	Transfinite Curve {6, 56, 310, 564}		= G_Hls Using Bump 0.01;
	Transfinite Curve {26, 166, 420, 674}	= G_Hls Using Bump 0.01;
	
	//Vertical line of Prim Slot	
	Transfinite Curve {2, 76, 330, 584} 	= G_HPrim Using Bump 1; 
	Transfinite Curve {10, 78, 332, 586} 	= G_HPrim Using Bump 1; 
	Transfinite Curve {30,188, 442, 696}	= G_HPrim Using Bump 1;	
	
	//Vertical line of d12
	Transfinite Curve {3, 98, 352, 606}		= G_Hd12 Using Bump 1;
	Transfinite Curve {14, 100, 354, 608}	 = G_Hd12 Using Bump 1;
	Transfinite Curve {34, 210, 464, 718}	= G_Hd12 Using Bump 1;
	
	//Vertical line of Secondary Slot
	Transfinite Curve {4, 120, 374, 628}	= G_HSec Using Bump 1;
	Transfinite Curve {18,122,376,630}	= G_HSec Using Bump 1;
	Transfinite Curve {38,232,486,740}	= G_HSec Using Bump 1;
	
	//Vertical line of d23
	Transfinite Curve {5, 142, 396, 650} = G_Hd23 Using Bump 1;
	Transfinite Curve {22,144, 398, 652} = G_Hd23 Using Bump 1;
	Transfinite Curve {42, 254, 508, 762} = G_Hd23 Using Bump 1;
		
	//Tangential line of Prim Slot Half 
	Transfinite Curve {59, 60, 82, 104, 126,  148}		= G_WPrim	 Using Bump 1;	
	Transfinite Curve {68, 64, 86, 108, 130, 152 } 		= G_WPrim	 Using Bump 1;
	Transfinite Curve {178, 174, 196, 218, 240, 262}	= G_WPrim	 Using Bump 1;
	
	//Tangential line of Prim-Sec vent
	Transfinite Curve {313, 314, 336, 358, 380, 402}	= G_WPS		 Using Bump 1;
	Transfinite Curve {322, 318, 340, 362, 384, 406}	= G_WPS		 Using Bump 1;
	Transfinite Curve {432, 428, 450, 472, 494, 516}	= G_WPS		 Using Bump 1;
	
	//Tangential line of Sec Slot Half 
	Transfinite Curve {567, 568, 590, 612, 634, 656}	= G_WSecHalf Using Bump 1;
	Transfinite Curve {576, 572, 594, 616, 638, 660}	= G_WSecHalf Using Bump 1;
	Transfinite Curve {686, 682, 704, 726, 748, 770}	= G_WSecHalf Using Bump 1;
	
	//Radial line of annulus gap
	Transfinite Curve {7,57,311,565}	= G_Rgap Using Bump 0.01;
	Transfinite Curve {8,55,309,563}	= G_Rgap Using Bump 0.01; 
	Transfinite Curve {12,77,311,585}	= G_Rgap Using Bump 0.01;
	Transfinite Curve {16,99,353,607}	= G_Rgap Using Bump 0.01;
	Transfinite Curve {20,121,375,629}	= G_Rgap Using Bump 0.01;
	Transfinite Curve {24,143,397,651}	= G_Rgap Using Bump 0.01;
	Transfinite Curve {48}				= G_Rgap Using Bump 0.01;
	
	//Radial line of Head face plate
	Transfinite Curve {27, 167, 421, 675} = G_RHead Using Bump 0.1;
	Transfinite Curve {28, 165, 419, 673} = G_RHead Using Bump 0.1;
	Transfinite Curve {32, 187, 441, 695} = G_RHead Using Bump 0.1;
	Transfinite Curve {36, 209, 463, 717} = G_RHead Using Bump 0.1;
	Transfinite Curve {40, 231, 485, 739} = G_RHead Using Bump 0.1;
	Transfinite Curve {44, 253, 507, 761} = G_RHead Using Bump 0.1;
	Transfinite Curve {49} 				  = G_RHead Using Bump 0.01;
	
	//Pintle Tip and Farfield 
	Transfinite Curve {46, 277, 531, 785} = G_RTip Using Bump 0.01;
	Transfinite Curve {47, 275, 529, 783} = G_RTip Using Bump 0.01;
	Transfinite Curve {50, 292, 546, 800} = G_RTip Using Bump 0.01;
	
	
	//Problems on surface 341 350 372 595
	//Surface cannot be meshed using the transfinite algo
	//Transfinite Surfaces	
	//annulus gap
	Transfinite Surface {73,327,581};	Transfinite Surface {65,319,573};	
	Transfinite Surface {87,341,595};	Transfinite Surface {109,363,617};	
	Transfinite Surface {131,385,639};	Transfinite Surface {153,407,661};	
	
	// Head face
	Transfinite Surface {183,437,691};	Transfinite Surface {175,429,683};
	Transfinite Surface {197,451,705};	Transfinite Surface {219,473,727};
	Transfinite Surface {749,495,241};	Transfinite Surface {517,526,771};
	
	//ls
	Transfinite Surface {61,315,569};	Transfinite Surface {69,323,577};
	Transfinite Surface {179,433,687};	
	
	// PrimSlot	
	Transfinite Surface {591,337,83}; 	Transfinite Surface {599,345,91};
	Transfinite Surface {201,455,709}; 	
	
	// d12	
	Transfinite Surface {613,359,105}; 	Transfinite Surface {621,367,113};
	Transfinite Surface {223,477,731};
	
	// SecSlot
	Transfinite Surface {635,381,127}; 	Transfinite Surface {643,389,135};
	Transfinite Surface {245,499,753}; 
	
	// d23
	Transfinite Surface {657,403,149}; 	Transfinite Surface {665,411,157};
	Transfinite Surface {775,521,267};
	
	//Periodic
	Transfinite Surface {9,74,328,582};  	Transfinite Surface {29,184,438,692};	
	Transfinite Surface {13,96,350,604}; 	Transfinite Surface {33,206,460,714};
	Transfinite Surface {626,372,118,17};	Transfinite Surface {736,482,228,37};
	Transfinite Surface {21,140,394,648}; 	Transfinite Surface {758,504,250,41};
	Transfinite Surface {25,162,416,670}; 	Transfinite Surface {780,526,272,45};
	Transfinite Surface {51,289,543,797};	Transfinite Surface {814,560,306,52};

	// Curve Tip
	Transfinite Surface {288,542,796};		Transfinite Surface {274,538,792};	
	Transfinite Surface {301,555,809};
	
	//Transfinite volume{1,13,25};	Transfinite volume{6,18,30};
	//Transfinite volume{2,14,26};	Transfinite volume{7,19,31};
	//Transfinite volume{3,15,27};	Transfinite volume{8,20,32};
	//Transfinite volume{4,16,28};	Transfinite volume{9,21,33};
	//Transfinite volume{5,17,29};	Transfinite volume{10,22,34};
	//Transfinite volume{};	Transfinite volume{};
	
	//Recombine surface 
	//annulus gap
	Recombine Surface {73,327,581};	Recombine Surface {65,319,573};	
	Recombine Surface {87,341,595};	Recombine Surface {109,363,617};	
	Recombine Surface {131,385,639};	Recombine Surface {153,407,661};	
	
	// Head face
	Recombine Surface {183,437,691};	Recombine Surface {175,429,683};
	Recombine Surface {197,451,705};	Recombine Surface {219,473,727};
	Recombine Surface {749,495,241};	Recombine Surface {517,526,771};
	
	//ls
	Recombine Surface {61,315,569};	Recombine Surface {69,323,577};
	Recombine Surface {179,433,687};	
	
	// PrimSlot	
	Recombine Surface {591,337,83}; 	Recombine Surface {599,345,91};
	Recombine Surface {201,455,709}; 	
	
	// d12	
	Recombine Surface {613,359,105}; 	Recombine Surface {621,367,113};
	Recombine Surface {223,477,731};
	
	// SecSlot
	Recombine Surface {635,381,127}; 	Recombine Surface {643,389,135};
	Recombine Surface {245,499,753}; 
	
	// d23
	Recombine Surface {657,403,149}; 	Recombine Surface {665,411,157};
	Recombine Surface {775,521,267};
	
	//Periodic
	Recombine Surface {9,74,328,582};  	Recombine Surface {29,184,438,692};	
	Recombine Surface {13,96,350,604}; 	Recombine Surface {33,206,460,714};
	Recombine Surface {626,372,118,17};	Recombine Surface {736,482,228,37};
	Recombine Surface {21,140,394,648}; 	Recombine Surface {758,504,250,41};
	Recombine Surface {25,162,416,670}; 	Recombine Surface {780,526,272,45};
	Recombine Surface {51,289,543,797};	Recombine Surface {814,560,306,52};
//ElseIF(d23==0)
//Else()
	//Print error message
EndIf//+
//+
