From 21b4d80ab190a5d42e9a041bd0fb2480e8645aae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?China=20Hagstr=C3=B6m?= Date: Thu, 13 Oct 2022 16:17:27 -0400 Subject: [PATCH] decreased step sizes, testing to find optimal domain size to avoid BC interactions --- SSME.py | 2 +- plumeSSMEJP.ipynb | 172 ++++++++++++++++++++++++++-------------------- 2 files changed, 98 insertions(+), 76 deletions(-) diff --git a/SSME.py b/SSME.py index 06d678b..c7c528f 100644 --- a/SSME.py +++ b/SSME.py @@ -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 diff --git a/plumeSSMEJP.ipynb b/plumeSSMEJP.ipynb index d86ba4b..07f747d 100644 --- a/plumeSSMEJP.ipynb +++ b/plumeSSMEJP.ipynb @@ -1,5 +1,26 @@ { "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, @@ -7,19 +28,6 @@ "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", @@ -490,6 +498,7 @@ "metadata": {}, "outputs": [], "source": [ + "println(\"started\")\n", "n_species = 53\n", "upper = 20000 #[m]\n", "lower = 16000 #[m]\n", @@ -497,7 +506,7 @@ "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", @@ -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", @@ -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", @@ -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", @@ -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", @@ -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", @@ -655,7 +666,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "d011e69b", "metadata": {}, "outputs": [], @@ -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", @@ -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", @@ -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", @@ -725,12 +736,12 @@ " 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", @@ -738,7 +749,7 @@ " ϕ = 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", @@ -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", @@ -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", @@ -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", @@ -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" ] }, { @@ -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": {