function foilgmsh(archi,alfa,yplus,eter,Re,M,T0,N,bump) %foilgmsh(archi,alfa,yplus,eter,Re,M,T0,N,bump) % %foilgmsh will create a GEO file with all necessary instructions to mesh an airfoil geometry only with hexahedra in Gmsh, departuring from a set of points defining the airfoil section contour. So far this code can only deal with monoelement airfoils and 2D simulations for finite volume CFD codes. %The on-screen output displays a summary of the input, some statistics of the mesh to be created and useful information for definition of the case in CFD softwares, specially Code_Saturne. % %- archi: string which defines the complete name of the airfoil coordinates file. The coordinates must be given in XFoil format starting from the trailing edge and circling counterclockwise. The coordinates don't need to be adimensionalized, but the program won't do it either, as the airfoil chord will be estimated automatically from it. The author consider's this is useful to compare further results with experimental data. The coordinates can be off-centered by small amounts. %The file must be located in the active directory. The output file will have the same name with the string '.geo' appended. If you wish to include text in your coordinates file, please add a percentage sign on the first column, otherwise an error will occur. % %- alfa: desired angle of attack in degrees. The mesh generated will be so in wind axis system, where wind velocity coincides in sense and direction with the positive X axis, so when defining inlet speed, Y and Z components should be zero. Also, the section plane coincides with the mesh XY plane. % %- yplus: desired y+ value around the airfoil. Spacing between susecquent cells normal to the airfoil surface is done with an authomatically defined geometric progression. % %- eter: string defining material medium where airfoil is submerged. Only three posibilities are allowed, 'a' for normal air, 'h' for distilled water and 'n' for nitrogen. % %- Re: Reynolds number based on airfoil chord and freestream speed. The reference chord and that of the airfoil in the coordinates file must coincide. % %- M: if the medium is air or nitrogen, it is the Mach number. If the medium is water, it is the freestream speed magnitude in m/s. Although the applied formulas hold for transonic and supersonic flows, the resulting mesh might not be suitable for those cases. Compressible fully subsonic flows should not be a problem. % %- T0: if the medium is air or nitrogen, it is the stagnation temperature in Kelvin degrees. If the medium is water, it is the non-stagnated flow temperature in Celsius degrees. % %- N: four integers vector whre N(1) is the number of hexahedra along the airfoil chord, therefore the whole contour of the airfoil will be discretized in 2*N(1) elements. N(2) is the number of elements normal to the chord and inside a semiellipse closely enclosing the airfoil (a value between 25 and 100 should suffice for most applications). N(3) is the number of elements into which the horizontal leading edge-inlet and trailing edge-outlet gaps will be discretized. N(4) is similar to N(3) but applied to the vertical gaps between the airfoil and the top/bottom walls defining the box. % %- bump: a value which defines the mesh concentration degree around the leading and trailing edge. If bump=1 the cell lenghts will be uniform along the airfoil contour, if bump>1 the elements will be concentrated around the midchord. If 0 0.001 Prog_t = Prog + ((Prog^4-2*Prog^3+Prog^2)*L+(Prog^3-Prog^2+Prog^N*(Prog-Prog^2))*Ymin)/((Prog-1)*Prog^N*Ymin*N+((1-2*Prog)*Prog^N+Prog^2)*Ymin); e = abs((Prog_t-Prog)/Prog); Prog = Prog_t; end if Prog < 1 Prog = 1; end end %%%%%%%%%%----%%%%%%%%%%----%%%%%%%%%%----%%%%%%%%%%----%%%%%%%%%%----%%%%%%%%%% function Ymin = ypar (yplus,cuerda,Re,M,T0,eter) %This function computes the minimum cell distance from a wall, according to equations for turbulent flow over a flat plate at zero incidence. It also computes some bonus data to define simulation parameters. switch eter case 'a' T = T0/(1+0.2*M^2); V = M*sqrt(1.4*287.074*T); muT = 1.458e-6*T^1.5/(110.4+T); rho = Re*muT/(V*cuerda); P = rho*287.074*T; P0 = P*(T0/T)^(1.4/0.4); printf ('\n-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_\n\n') printf ('Medium: air\nRe = %g\nChord = %f [m]\nDensity = %g [Kg/m^3]\nDynamic viscosity = %g [Kg/(ms)]\nFreestream speed = %f [m/s]\nFreestream Mach = %f\nStatic pressure (absolute) = %f [Pa]\nStagnation pressure = %f [Pa]\nTemperature = %f [K]\nStagnation temeprature = %f[K]\nY+ = %f\n\n', Re,cuerda,rho,muT,V,M,P,P0,T,T0,yplus) case 'n' T = T0/(1+0.2*M^2); V = M*sqrt(1.4*297*T); muT = 1.781e-5*(111+300.55)/(111+T)*(T/300.55)^1.5; rho = Re*muT/(V*cuerda); P = rho*297*T; P0 = P*(T0/T)^(1.4/0.4); printf ('\n-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_\n\n') printf ('Medium: nitrogen\nRe = %g\nChord = %f [m]\nDensity = %g [Kg/m^3]\nDynamic viscosity = %g [Kg/(ms)]\nFreestream speed = %f [m/s]\nFreestream Mach = %f\nStatic pressure (absolute) = %f [Pa]\nStagnation pressure = %f [Pa]\nTemperature = %f [K]\nStagnation temperature = %f[K]\nY+ = %f\n\n', Re,cuerda,rho,muT,V,M,P,P0,T,T0,yplus) case 'h' V = M; rho = 1000*(1-(T0+288.9414)/(508929.2*(T0+68.12963))*(T0-3.9863)^2); muT = rho*V*cuerda/Re; printf ('\n-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_\n\n') printf ('Medium: water\nRe = %g\nChord = %f [m]\nDensity = %g [Kg/m^3]\nDynamic viscosity = %g [Kg/(ms)]\nFreestream speed = %f [m/s]\nTemperature = %f [K]\nY+ = %f\n\n', Re,cuerda,rho,muT,V,yplus) end Cf = 0.02; while i<10 funcion = 4.15*sqrt(Cf)*log10(Re*Cf)+1.7*sqrt(Cf)-1.0; derfunc = (4.15*log10(exp(1.0))+0.5*4.15*log10(Re*Cf)+1.7/2.0)/sqrt(Cf); fsd = funcion/derfunc; if abs(fsd/Cf) <= exp(-5.0) break end Cfo = Cf-fsd; if Cfo <= 0.0 Cf = 0.5*Cf; else Cf = Cfo; end i=i+1; end %Cfo = (1./(4.15*log10(Re*Cf)+1.7))^2; %tau = 0.5*rho*V*V*Cf %aus = sqrt(tau/rho) Ymin = yplus*muT/(V*sqrt(Cf/2)); printf ('To obtain CFL(max) <= 20 across the whole flowfield, a timestep dt <= %.10gs is recomended for transient simulations.\n\n', 20*Ymin/V) printf ('The following values are recomended to initialize external flowfield variables:\nK = %g [m^2/s^2]\nEpsilon = %g [m^2/s^3]\nOmega = %g [1/s]\nTurbulent intensity = 6.6667e-7 [%%]\nHydraulic diameter = %f [m]\n\n', 1e-6*V^2,4.5e-7*V^3/cuerda,0.45*V/cuerda,0.0052164*cuerda) end