Skip to content

Commit

Permalink
decreased step sizes, testing to find optimal domain size to avoid BC…
Browse files Browse the repository at this point in the history 
… interactions
  • Loading branch information
@chinahg
chinahg committed Oct 20, 2022
1 parent 7c6433e commit 21b4d80
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 76 deletions.
2 changes: 1 addition & 1 deletion SSME.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@
u = results_u[1] #velocity at exit of nozzle (m/s)
gamma = 1.1
print(u)
P1 = P_e #Noz_states.P[n] PLACEHOLDER, need design exit pressure or design alt
P1 = P_e #TODO Noz_states.P[n] PLACEHOLDER, need design exit pressure or design alt
T1 = Noz_states.T[n]
M1 = u/math.sqrt(gamma*P1/Noz_states.density[n])
P2 = P_atm #Pa
Expand Down
172 changes: 97 additions & 75 deletions plumeSSMEJP.ipynb
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,33 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "0f5cbf93",
"metadata": {},
"outputs": [],
"source": [
"#import Pkg\n",
"#\n",
"#Pkg.add(\"Interpolations\")\n",
"#Pkg.add(\"PyPlot\")\n",
"#Pkg.add(\"PyCall\")\n",
"#Pkg.add(\"OrdinaryDiffEq\")\n",
"#Pkg.add(\"YAML\")\n",
"#Pkg.add(\"DelimitedFiles\")\n",
"#Pkg.add(\"CSV\")\n",
"#Pkg.add(\"HDF5\")\n",
"#Pkg.add(\"StructArrays\")\n",
"#Pkg.add(\"Random\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f73d713f",
"metadata": {},
"outputs": [],
"source": [
"import Pkg\n",
"\n",
"Pkg.add(\"Interpolations\")\n",
"Pkg.add(\"PyPlot\")\n",
"Pkg.add(\"PyCall\")\n",
"Pkg.add(\"OrdinaryDiffEq\")\n",
"Pkg.add(\"YAML\")\n",
"Pkg.add(\"DelimitedFiles\")\n",
"Pkg.add(\"CSV\")\n",
"Pkg.add(\"HDF5\")\n",
"Pkg.add(\"StructArrays\")\n",
"Pkg.add(\"Random\")\n",
"\n",
"using Interpolations, PyCall, OrdinaryDiffEq,\n",
"YAML, DelimitedFiles, CSV, HDF5, StructArrays, Random, \n",
"NBInclude, PyPlot\n",
Expand Down Expand Up @@ -490,14 +498,15 @@
"metadata": {},
"outputs": [],
"source": [
"println(\"started\")\n",
"n_species = 53\n",
"upper = 20000 #[m]\n",
"lower = 16000 #[m]\n",
"space = convert(Int, (upper - lower) / 2000 + 1) #250\n",
"h = [16000] #Int.(LinRange(16000, 20000, space))\n",
"g = 1\n",
"\n",
"n = 300 #90 #n steps in x dir\n",
"n = 200 #90 #n steps in x dir\n",
"s = 900 #80 #s steps in y direction\n",
"T_save = zeros(length(h), s, n)\n",
"u_save = zeros(length(h), s, n)\n",
Expand All @@ -507,6 +516,7 @@
"gas_g = StructArray{gas_type}(undef,s,n,length(h))\n",
"gas_g .= [gas_type(0)]\n",
"\n",
"println(\"started altitude for loop\")\n",
"for m = 1:length(h)\n",
" h_string = string(h[m])\n",
"\n",
Expand All @@ -526,18 +536,19 @@
" #curve fit #a = ambient vel [m/s] (speed of rocket) \n",
" T_a = HDF5.h5read(\"/home/chinahg/GCresearch/rocketemissions/plot_data.h5\", h_string * \"m/T_a\")\n",
" #T0 = T_a #for testing see vincent messages\n",
" u0 = u_a #for testing see vincent messages\n",
" #u0 = u_a #for testing see vincent messages\n",
" \n",
" ### CALCULATE VELOCITY AND TEMPERATURE FIELDS (NO CHEMISTRY) ###\n",
"\n",
" Δϕ = 0.1 * ones(n) #step size in phi\n",
" println(\"started populating phi and psi\")\n",
" Δϕ = 0.01 * ones(n) #step size in phi\n",
" for i = 2:n\n",
" Δϕ[i] = 1.05 * Δϕ[i-1] #enlarge with each step by 1.1\n",
" end\n",
"\n",
" Δψ = 0.1 * ones(s) #s vertical grid points\n",
" Δψ = 0.001 * ones(s) #s vertical grid points\n",
" for i = 2:s\n",
" Δψ[i] = 1.05 * Δψ[i-1] #enlarge with each step by 1.1\n",
" Δψ[i] = 1.01 * Δψ[i-1] #enlarge with each step by 1.1\n",
" end\n",
"\n",
" Pr = 1.0 #prandtl number\n",
Expand All @@ -549,7 +560,7 @@
" y_init = compute_y(u_init, T_init, Δψ, R, p)\n",
"\n",
" #Geometry of plume\n",
" radius = 0.2\n",
" radius = 1.147\n",
" u_init[y_init.>radius] .= convert(AbstractFloat, u_a)\n",
" T_init[y_init.>radius] .= convert(AbstractFloat, T_a)\n",
" \n",
Expand Down Expand Up @@ -594,7 +605,7 @@
" χ_init = zeros(s, n_species)\n",
"\n",
" ### COMPUTE MOLE FRACTIONS IN PLUME ###\n",
"\n",
" println(\"at radius allocation\")\n",
" i = y_init .>= radius\n",
" for j = 1:length(i) \n",
" if i[j] == 0 #if y position is less than the initial plume radius, assign plume conditions\n",
Expand Down Expand Up @@ -622,11 +633,11 @@
" χ_1 = zeros(size(χ_init))\n",
" \n",
" # Create gas object to store Reactor output gas object state\n",
" gas = ct.Solution(\"gri30.yaml\")\n",
" #gas = ct.Solution(\"gri30.yaml\")\n",
"\n",
" # Create a dummy reactor to ????\n",
" dummy_reactor = ct.IdealGasReactor(gas)\n",
"\n",
" #dummy_reactor = ct.IdealGasReactor(gas)\n",
" println(\"starting splitting\")\n",
" for i = 1:n-1 #x\n",
" \n",
" for j = 1:n_species #species\n",
Expand Down Expand Up @@ -655,7 +666,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": null,
"id": "d011e69b",
"metadata": {},
"outputs": [],
Expand All @@ -678,7 +689,7 @@
" R = 287 #[J/kgK] placeholder\n",
"\n",
" #REGRID SOLUTION\n",
" #xx, yy, u_g, T_g, χ_gO2 = regrid_solution(x, y, u, T, χ[:,:,4], 0.01)\n",
" xx, yy, u_g, T_g, χ_gO2 = regrid_solution(x, y, u, T, χ[:,:,4], 0.01)\n",
" #xx, yy, u_g, T_g, χ_gN2 = regrid_solution(x, y, u, T, χ[:,:,48], 0.01)\n",
" #xx, yy, u_g, T_g, χ_gNO2 = regrid_solution(x, y, u, T, χ[:, :, 37], 0.01)\n",
" #xx, yy, u_g, T_g, χ_gN2O = regrid_solution(x, y, u, T, χ[:, :, 38], 0.01)\n",
Expand All @@ -690,7 +701,7 @@
"\n",
" #CALCULATE NO EI FOR ALTITUDE\n",
" # Sum of all NOx species as a function of x\n",
" n = 300\n",
" n = 200\n",
" s = 900\n",
" sumNO = zeros(n)\n",
" sumNO2 = zeros(n)\n",
Expand All @@ -713,9 +724,9 @@
" mdot_Ar_sum = zeros(n)\n",
" mdot_N2_sum = zeros(n)\n",
" mdot_O2_sum = zeros(n)\n",
" mdot_Ar_sumTrunc = zeros(n)\n",
" mdot_N2_sumTrunc = zeros(n)\n",
" mdot_O2_sumTrunc = zeros(n)\n",
" #mdot_Ar_sumTrunc = zeros(n)\n",
" #mdot_N2_sumTrunc = zeros(n)\n",
" #mdot_O2_sumTrunc = zeros(n)\n",
" mdot_Ar = zeros(s,n)\n",
" mdot_N2 = zeros(s,n)\n",
" mdot_O2 = zeros(s,n)\n",
Expand All @@ -725,20 +736,20 @@
" mdot_tot = zeros(n)\n",
" ppm_tot = zeros(n)\n",
"\n",
" Δψ = 0.1 * ones(s) #s vertical grid points\n",
" Δψ = 0.001 * ones(s) #s vertical grid points\n",
" for i = 2:s\n",
" Δψ[i] = 1.05 * Δψ[i-1] #enlarge with each step by 1.1\n",
" Δψ[i] = 1.01 * Δψ[i-1] #enlarge with each step by 1.1\n",
" end\n",
"\n",
" Δϕ = 0.1 * ones(n) #step size in phi\n",
" Δϕ = 0.01 * ones(n) #step size in phi\n",
" for i = 2:n\n",
" Δϕ[i] = 1.05 * Δϕ[i-1] #enlarge with each step by 1.1\n",
" end\n",
"\n",
" ϕ = cumsum(Δϕ)\n",
" ψ = cumsum(Δψ)\n",
"\n",
" for i = 1:n, j = 1:s-300\n",
" for i = 1:n, j = 1:s-1\n",
" rho_tot[j,i] = P_atm/(R*T[j,i]) #[kg/m^3] #need to convert to xx yy coordinates\n",
" MFar[j,i] = (χ[j,i,49]*(1/28.97)*(39.9))/1000000 #[kg/kg tot]\n",
" MFn2[j,i] = (χ[j,i,48]*(1/28.97)*(28))/1000000 #[kg/kg tot]\n",
Expand All @@ -753,15 +764,16 @@
" ppm_n2[i] += χ[j,i,48]\n",
" ppm_o2[i] += χ[j,i,4]\n",
"\n",
" #TODO: populate the last entry of sum arrays (stops at s-1)\n",
" mdot_O2_sum[i] += mdot_O2[j,i]\n",
" mdot_N2_sum[i] += mdot_N2[j,i]\n",
" mdot_Ar_sum[i] += mdot_Ar[j,i]\n",
" \n",
" if j < 250\n",
" mdot_O2_sumTrunc[i] += mdot_O2[i,j]\n",
" mdot_N2_sumTrunc[i] += mdot_N2[i,j]\n",
" mdot_Ar_sumTrunc[i] += mdot_Ar[i,j]\n",
" end\n",
" #if j < 150\n",
" # mdot_O2_sumTrunc[i] += mdot_O2[i,j]\n",
" # mdot_N2_sumTrunc[i] += mdot_N2[i,j]\n",
" # mdot_Ar_sumTrunc[i] += mdot_Ar[i,j]\n",
" #end\n",
" mdot_tot[i] += mdot_Ar[j,i] + mdot_N2[j,i] + mdot_O2[j,i]\n",
" end\n",
"\n",
Expand All @@ -779,7 +791,7 @@
" axc.plot(ϕ, u[s-200,:], label=\"100/300\")\n",
" axc.plot(ϕ, u[s-250,:], label=\"50/300\")\n",
" axc.plot(ϕ, u[s-299,:], label=\"1/300\")\n",
"\n",
"#\n",
" axc.set_xlabel(\"ϕ\")\n",
" axc.set_ylabel(\"m/s at ψ = 0.1\")\n",
" axc.set_xscale(\"log\")\n",
Expand Down Expand Up @@ -823,53 +835,53 @@
" savefig(\"/home/chinahg/GCresearch/rocketemissions/rockettests/\" * h_string * \"m/u.png\")\n",
"\n",
"\n",
" fig,axs = plt.subplots(2,2)\n",
" axs[1,1].plot(ϕ, mdot_tot)\n",
" axs[1,1].set_xlabel(\"ϕ\")\n",
" axs[1,1].set_ylabel(\"Total mass flow rate [kg/s]\")\n",
" axs[1,1].set_title(\"Total mdot\")\n",
" axs[1,1].set_xscale(\"log\")\n",
"\n",
" axs[1,2].plot(ϕ, mdot_N2_sumTrunc)\n",
" axs[1,2].set_xlabel(\"ϕ\")\n",
" axs[1,2].set_ylabel(\"N2 mass flow rate [kg/s]\")\n",
" axs[1,2].set_title(\"N2\")\n",
" axs[1,2].set_xscale(\"log\")\n",
"\n",
" axs[2,1].plot(ϕ, mdot_O2_sumTrunc)\n",
" axs[2,1].set_xlabel(\"ϕ\")\n",
" axs[2,1].set_ylabel(\"O2 mass flow rate [kg/s]\")\n",
" axs[2,1].set_xscale(\"log\")\n",
"\n",
" axs[2,2].plot(ϕ, mdot_Ar_sumTrunc)\n",
" axs[2,2].set_xlabel(\"ϕ\")\n",
" axs[2,2].set_ylabel(\"Ar mass flow rate [kg/s]\")\n",
" axs[2,2].set_title(\"Ar\")\n",
" axs[2,2].set_xscale(\"log\")\n",
" \n",
" fig.suptitle(\"Truncated Mass Flow Rates, s = 250/300\")\n",
" fig.tight_layout()\n",
" savefig(\"/home/chinahg/GCresearch/rocketemissions/rockettests/\" * h_string * \"m/mdot_trunc.png\")\n",
" #fig,axs = plt.subplots(2,2)\n",
" #axs[1,1].plot(ϕ, mdot_tot)\n",
" #axs[1,1].set_xlabel(\"ϕ\")\n",
" #axs[1,1].set_ylabel(\"Total mass flow rate [kg/s]\")\n",
" #axs[1,1].set_title(\"Total mdot\")\n",
" #axs[1,1].set_xscale(\"log\")\n",
"#\n",
" #axs[1,2].plot(ϕ, mdot_N2_sumTrunc)\n",
" #axs[1,2].set_xlabel(\"ϕ\")\n",
" #axs[1,2].set_ylabel(\"N2 mass flow rate [kg/s]\")\n",
" #axs[1,2].set_title(\"N2\")\n",
" #axs[1,2].set_xscale(\"log\")\n",
"#\n",
" #axs[2,1].plot(ϕ, mdot_O2_sumTrunc)\n",
" #axs[2,1].set_xlabel(\"ϕ\")\n",
" #axs[2,1].set_ylabel(\"O2 mass flow rate [kg/s]\")\n",
" #axs[2,1].set_xscale(\"log\")\n",
"#\n",
" #axs[2,2].plot(ϕ, mdot_Ar_sumTrunc)\n",
" #axs[2,2].set_xlabel(\"ϕ\")\n",
" #axs[2,2].set_ylabel(\"Ar mass flow rate [kg/s]\")\n",
" #axs[2,2].set_title(\"Ar\")\n",
" #axs[2,2].set_xscale(\"log\")\n",
" #\n",
" #fig.suptitle(\"Truncated Mass Flow Rates, s = 250/300\")\n",
" #fig.tight_layout()\n",
" #savefig(\"/home/chinahg/GCresearch/rocketemissions/rockettests/\" * h_string * \"m/mdot_trunc.png\")\n",
"\n",
" fig,axX = plt.subplots(2,2)\n",
" axX[1,1].plot(ϕ[1:10], mdot_tot[1:10])\n",
" axX[1,1].plot(ϕ, mdot_tot)\n",
" axX[1,1].set_xlabel(\"ϕ\")\n",
" axX[1,1].set_ylabel(\"Total mass flow rate [kg/s]\")\n",
" axX[1,1].set_title(\"Total mdot\")\n",
" axX[1,1].set_xscale(\"log\")\n",
"\n",
" axX[1,2].plot(ϕ[1:10], mdot_N2_sum[1:10])\n",
" axX[1,2].plot(ϕ, mdot_N2_sum)\n",
" axX[1,2].set_xlabel(\"ϕ\")\n",
" axX[1,2].set_ylabel(\"N2 mass flow rate [kg/s]\")\n",
" axX[1,2].set_title(\"N2\")\n",
" axX[1,2].set_xscale(\"log\")\n",
"\n",
" axX[2,1].plot(ϕ[1:10], mdot_O2_sum[1:10])\n",
" axX[2,1].plot(ϕ, mdot_O2_sum)\n",
" axX[2,1].set_xlabel(\"ϕ\")\n",
" axX[2,1].set_ylabel(\"O2 mass flow rate [kg/s]\")\n",
" axX[2,1].set_xscale(\"log\")\n",
"\n",
" axX[2,2].plot(ϕ[1:10], mdot_Ar_sum[1:10])\n",
" axX[2,2].plot(ϕ, mdot_Ar_sum)\n",
" axX[2,2].set_xlabel(\"ϕ\")\n",
" axX[2,2].set_ylabel(\"Ar mass flow rate [kg/s]\")\n",
" axX[2,2].set_title(\"Ar\")\n",
Expand Down Expand Up @@ -907,15 +919,17 @@
" fig.tight_layout()\n",
" savefig(\"/home/chinahg/GCresearch/rocketemissions/rockettests/\" * h_string * \"m/ppm_tot.png\")\n",
"\n",
" println(\"Δ mdot Ar total is \", mdot_Ar_sum[n]-mdot_Ar_sum[1])\n",
" println(\"Δ mdot Ar/Δ mdot total is \", (mdot_Ar_sum[n]-mdot_Ar_sum[1])/(mdot_tot[n]-mdot_tot[1]))\n",
" println(mdot_Ar_sum[300], mdot_tot[300])\n",
" println(mdot_Ar_sum[n], mdot_tot[n])\n",
"\n",
" println(\"plotted \", h_string[m], \" altitude!\\n\")\n",
" println(mdot_Ar_sum[1])\n",
" println(mdot_Ar[:,1])\n",
" #println(mdot_Ar_sum[1])\n",
" #println(mdot_Ar[:,1])\n",
" println(yy[473])\n",
"end\n",
"\n",
"println(\"\\n done plotting!\")\n"
"println(\"done plotting!\")\n"
]
},
{
Expand All @@ -926,8 +940,16 @@
"outputs": [],
"source": [
"#could be due to error from interpolating between points \n",
"#(areas increase exponentially so error increases as well)"
"#(areas increase exponentially so error increases as well)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3e6aeb8e",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down

0 comments on commit 21b4d80

Please sign in to comment.