Flow Through a de Laval Nozzle Appendix

7. Appendix Included in this section are the parameter files that were used to run the supersonic de Laval nozzle in addition to the custom python script used to generate the structured mesh for the grid. The actual grid file is not included as it is remarkably large. The files are listed here in the order in which they are presented above. import numpy as np import matplotlib.pyplot as plt import pdb ''' makewedge.py PURPOSE: Script that writes instructions to create a wedge shape for the CFD openfoam project. For use with Gmsh AUTHOR: Connor Robinson, Dec 6th, 2016 ''' #Set up paths for saving path = '/Users/Connor/Desktop/Coursework/CFD/project/' outfile = 'mesh_shape.txt' #Define the shape of the nozzle in 1D #Set the number of cells along the cylinder length = 400 #Set up cells between the center of the cylinder and the axis #Specify number of cells hnum = 200 #Specify the transfinite progressive scaling for those cells hprog = .97 #Define r as a function of Z zbase = np.arange(length) r = (45-zbase**2/(.08*length)**2 + zbase**3/(.186*length)**3) r = (r/np.max(r))[::-1] #Reverse r (Better nozzle shape, not necessary) #Scale zbase back to something reasonable. zbase = zbase/(length/2) #Define the fraction of a circle for the wedge to include (1/144 = 2.5 degrees) anglefrac = 1/144 #Actually go about generating the mesh, unlikely necessary to change anything below here. anglepoints = 3 phi = [-np.pi*anglefrac,np.pi*anglefrac] x = np.array([r * np.cos(p) for p in phi]).flatten() y = np.array([r * np.sin(p) for p in phi]).flatten() z = np.hstack([zbase, zbase, zbase]) x = np.hstack([x,np.zeros(length)]) y = np.hstack([y,np.zeros(length)]) vlines = [] #Give a number to each point points = np.arange(length*anglepoints)+1 counter = 1 #Connect the lines along the flow axis, start at bottom and work around with increasing phi for angle in np.arange(anglepoints): for xval in np.arange(length-1): vlines.append([counter, xval+ angle*length+1, xval+angle*length+2]) counter = counter +1 hlines = [] nvert = counter -1 #Connect the lines around the tube, start at bottom and work way up. for zval in np.arange(length)+1: for angle in np.arange(anglepoints): if (zval == length) and (angle == anglepoints -2): #hlines.append([counter, zval+angle*length, length*anglepoints]) hlines.append([counter, length*anglepoints, zval+angle*length]) else: if angle == 2 or angle == 0: hlines.append([counter, zval+angle*length, np.mod(zval+(angle+1)*length, length*anglepoints)]) if angle == 1: hlines.append([counter, np.mod(zval+(angle+1)*length, length*anglepoints), zval+angle*length]) counter = counter+1 nhorz = counter - 1 -nvert #Now form the line loops, start at bottom and work around with increasing phi loop = [] loopcounter = 1 for zval in np.arange(length-1): for angle in np.arange(anglepoints): loop.append([loopcounter,\ vlines[zval+(angle*(length-1))][0], \ hlines[(zval+1)*anglepoints + angle][0], \ vlines[np.mod(zval + (angle+1)*(length-1), (anglepoints)*(length-1))][0], \ hlines[zval*anglepoints + angle][0]]) loopcounter = loopcounter +1 #Write out the lines and line loops to a text file. file = open(path+outfile, 'w') #Write each set of points file.write('// Points \n') for p in np.arange(length*anglepoints): file.write('Point ('+str(np.int(p+1))+') = {'+str(x[p])+','+str(y[p])+','+str(z[p])+',1};\n') #Write the vertical lines file.write('// Verticle lines \n') for v in vlines: file.write('Line('+str(np.int(v[0]))+') = {'+str(np.int(v[1]))+','+str(np.int(v[2]))+'};\n') #Write the horizontal lines file.write('// Horizontal lines \n') for h in hlines: file.write('Line('+str(np.int(h[0]))+') = {'+str(np.int(h[1]))+','+str(np.int(h[2]))+'};\n') #Write the line loops surfacecounter = 1 file.write('// Line loops \n') for i, L in enumerate(loop): if i != len(loop)-2: if np.mod(i+1, 3) == 0: file.write('Line Loop('+str(np.int(L[0]))+') = {'+str(np.int(L[1]))+','+str(np.int(L[2]))+','+str(np.int(-L[3]))+','+str(np.int(-L[4] ))+'};\n') if np.mod(i+2, 3) == 0: file.write('Line Loop('+str(np.int(L[0]))+') = {'+str(np.int(L[1]))+','+str(-np.int(L[2]))+','+str(np.int(-L[3]))+','+str(np.int(L[4] ))+'};\n') if np.mod(i+3, 3) == 0: file.write('Line Loop('+str(np.int(L[0]))+') = {'+str(np.int(L[1]))+','+str(np.int(L[2]))+','+str(np.int(-L[3]))+','+str(np.int(-L[4] ))+'};\n') #Handle the final loop... else: file.write('Line Loop('+str(np.int(L[0]))+') = {'+str(np.int(L[1]))+','+str(-np.int(L[2]))+','+str(np.int(-L[3]))+','+str(np.int(L[4] ))+'};\n') file.write('Ruled Surface('+str(surfacecounter)+') = {'+str(surfacecounter)+'};\n') surfacecounter = surfacecounter +1 #Set up side1, side2, and the outer wall #Side1 string ='' file.write('Physical Surface("Left") = {') for val in np.arange(loop[-1][0])+1: if np.mod(val, 3) == 0: string = string+(np.str(np.int(val))+',') file.write(string[:-1]+'};\n') #Side2 string ='' file.write('Physical Surface("Right") = {') for val in np.arange(loop[-1][0])+1: if np.mod(val+1, 3) == 0: string = string+(np.str(np.int(val))+',') file.write(string[:-1]+'};\n') #outerwall string ='' file.write('Physical Surface("Outerwall") = {') for val in np.arange(loop[-1][0])+1: if np.mod(val+2, 3) == 0: string = string+(np.str(np.int(val))+',') file.write(string[:-1]+'};\n') #Set up line loop and empty surface for axis: string = '' file.write('Physical Line("Axis") = {') for val in np.arange(length-1)+(length)*2-1: string = string+(np.str(np.int(val))+',') file.write(string[:-1]+'};\n') #Form a volume for each cell in the center #Write the internal line interiorloopstart = surfacecounter for val in np.arange(length): string = 'Line loop('+str(surfacecounter)+') = {' baseval = (val+length-1)*3 string = string + str(-int(baseval+3))+','+str(int(baseval+2))+','+str(-int(baseval+1)) file.write(string+'};\n') file.write('Ruled Surface('+str(surfacecounter)+') = {'+str(surfacecounter)+'};\n') surfacecounter = surfacecounter+1 file.write('Physical Surface("Input") = {'+str(interiorloopstart)+'};\n') file.write('Physical Surface("Output") = {'+str(surfacecounter-1)+'};\n') volcounter = 1 for val in (np.arange(length-1)): string string string string string string = = = = = = 'Surface loop('+str(volcounter)+') = {' string+(np.str(np.int(val*3+1))+',') string+(np.str(np.int(val*3+2))+',') string+(np.str(np.int(val*3+3))+',') string+(np.str(np.int(interiorloopstart+val))+',') string+(np.str(np.int(interiorloopstart+val+1))) file.write(string+'};\n') file.write('Volume('+str(volcounter)+') = {'+str(volcounter)+'};\n') volcounter = volcounter+1 #Set up the mesh structure for the horizontal lines string = '' file.write('Transfinite Line {') for val in np.arange(nhorz): #if val != nhorz-1: if np.mod(val, 3) != 0: string = string+str(int(val+nvert+1))+',' file.write(string[:-1]+'} = '+str(hnum)+' Using Progression '+str(hprog)+';\n') #Set up the mesh for the outer wall string = '' file.write('Transfinite Line {') for val in np.arange(nhorz): #if val != nhorz-1: if np.mod(val, 3) == 0: string = string+str(int(val+nvert+1))+',' file.write(string[:-1]) file.write('} = 1;\n') #Set up the mesh structure for the vertical lines file.write('Transfinite Line {') for val in np.arange(nvert): if val != nvert-1: #if np.mod(val, 3) != 0 and val != nhorz-1: file.write(str(int(val+1))+',') if val == nvert-1: #if (np.mod(val, 3) != 0) and val == nhorz-1: file.write(str(int(val+1))) file.write('} = 1;\n') file.write('Transfinite Surface "*";\n') file.write('Recombine Surface "*";\n') file.write('Transfinite Volume "*";\n') string = 'Physical Volume("internal") = {' for val in np.arange(volcounter-1)+1: string = string+str(int(val))+',' file.write(string[:-1]+'};\n') file.close() /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // fluxScheme ddtSchemes { default } gradSchemes { default } divSchemes { default div(tauMC) } laplacianSchemes { default } Kurganov; Euler; Gauss linear; none; Gauss linear; Gauss linear corrected; interpolationSchemes { default linear; reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; } snGradSchemes { default corrected; } // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { "(rho|rhoU|rhoE)" { solver } diagonal; U { solver smoother nSweeps tolerance relTol smoothSolver; GaussSeidel; 2; 1e-10; 0; $U; tolerance relTol 1e-10; 0; } e { } } // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object thermophysicalProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // thermoType { type mixture transport thermo equationOfState specie energy } mixture { specie { nMoles molWeight } thermodynamics { Cp Hf } transport { As Ts Pr } } hePsiThermo; pureMixture; sutherland; hConst; perfectGas; specie; sensibleInternalEnergy; 1; 28.96; 1004.5; 0; 1.458e-06; 110.4; 1; // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType laminar; // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class polyBoundaryMesh; location "constant/polyMesh"; object boundary; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 5 ( Input { type physicalType nFaces startFace } Output { type physicalType nFaces startFace } Outerwall { type physicalType nFaces startFace } Right { type physicalType nFaces startFace } Left { type patch; patch; 199; 158204; patch; patch; 199; 158403; patch; patch; 399; 158602; wedge; wedge; 79401; 159001; wedge; physicalType nFaces startFace wedge; 79401; 238402; } ) // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [1 -1 -2 0 0 0 0]; internalField uniform 15000; boundaryField { Output { type waveTransmissive; field p; phi phi; rho rho; psi thermo:psi; gamma 1.4; fieldInf 2533.33; lInf 1; value uniform 2533.33; } Outerwall { type } slip; Input { type U phi rho totalPressure; U; phi; none; psi gamma p0 value none; 1.4; uniform 101325; uniform 101325; } Left {type wedge;} Right {type wedge;} Axis {type empty;} } // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volVectorField; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 -1 0 0 0 0]; internalField uniform (5 0 0); boundaryField { Input { type } Output { type } Outerwall { type } zeroGradient; zeroGradient; slip; Left {type wedge;} Right {type wedge;} Axis {type empty;} } // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object T; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 0 1 0 0 0]; internalField uniform 298.0; boundaryField { Input { type value } Output { type } Outerwall { type } fixedValue; uniform 298.0; zeroGradient; slip; Left {type wedge;} Right {type wedge;} Axis {type empty;} } // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application rhoCentralFoam; startFrom latestTime; startTime 0.0161; stopAt endTime; endTime 2e-02; deltaT 1e-10; writeControl adjustableRunTime; writeInterval 1e-04; cycleWrite 0; writeFormat ascii; writePrecision 15; writeCompression off; timeFormat general; timePrecision 6; adjustTimeStep yes; maxCo 0.5; maxDeltaT 1; // ************************************************************************* // /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object sampleDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // interpolationScheme cellPoint; setFormat raw; sets ( horizontalLine { type uniform; axis distance; start ( 0.1 0.0 0.005 ); end ( 0.1 0.0 1.995 ); nPoints 100; } ); fields ( p U T); // ************************************************************************* //