1 package com.billkuker.rocketry.motorsim;
\r
5 import java.util.SortedMap;
\r
6 import java.util.TreeMap;
\r
8 import javax.measure.quantity.Area;
\r
9 import javax.measure.quantity.Dimensionless;
\r
10 import javax.measure.quantity.Duration;
\r
11 import javax.measure.quantity.Force;
\r
12 import javax.measure.quantity.Length;
\r
13 import javax.measure.quantity.Mass;
\r
14 import javax.measure.quantity.MassFlowRate;
\r
15 import javax.measure.quantity.Pressure;
\r
16 import javax.measure.quantity.Temperature;
\r
17 import javax.measure.quantity.Velocity;
\r
18 import javax.measure.quantity.Volume;
\r
19 import javax.measure.quantity.VolumetricDensity;
\r
20 import javax.measure.unit.SI;
\r
22 import org.apache.log4j.Logger;
\r
23 import org.jscience.physics.amount.Amount;
\r
24 import org.jscience.physics.amount.Constants;
\r
26 import com.billkuker.rocketry.motorsim.fuel.KNSB;
\r
27 import com.billkuker.rocketry.motorsim.fuel.KNSU;
\r
28 import com.billkuker.rocketry.motorsim.grain.CompoundGrain;
\r
29 import com.billkuker.rocketry.motorsim.grain.CoredCylindricalGrain;
\r
30 import com.billkuker.rocketry.motorsim.grain.MultiGrain;
\r
31 import com.billkuker.rocketry.motorsim.visual.BurnPanel;
\r
35 private static Logger log = Logger.getLogger(Burn.class);
\r
36 protected final Motor motor;
\r
38 private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);
\r
41 private static double combustionEfficency = 0.97;
\r
43 private static double densityRatio = 0.96;
\r
45 public class Interval{
\r
46 Amount<Duration> time;
\r
47 public Amount<Length> regression;
\r
48 public Amount<Pressure> chamberPressure;
\r
49 Amount<Mass> chamberProduct;
\r
50 Amount<Force> thrust;
\r
52 public String toString(){
\r
53 return time + " " + regression + " " + chamberPressure + " " + chamberProduct;
\r
57 protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();
\r
59 public SortedMap<Amount<Duration>,Interval> getData(){
\r
63 public Motor getMotor(){
\r
67 public Amount<Duration> burnTime(){
\r
68 return data.lastKey();
\r
71 public Burn(Motor m){
\r
75 private void burn(){
\r
76 log.info("Starting burn...");
\r
78 Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);
\r
80 //if ( motor.getGrain() instanceof Grain.DiscreteRegression )
\r
81 //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();
\r
83 Interval initial = new Interval();
\r
84 initial.time = Amount.valueOf(0, SI.SECOND);
\r
85 initial.regression = Amount.valueOf(0, SI.MILLIMETER);
\r
86 initial.chamberPressure = atmosphereicPressure;
\r
87 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);
\r
88 initial.thrust = Amount.valueOf(0, SI.NEWTON);
\r
90 data.put(Amount.valueOf(0, SI.SECOND), initial);
\r
92 for ( int i = 0; i < 5000; i++ ){
\r
94 Interval prev = data.get(data.lastKey());
\r
96 log.debug("Step " + i + " ==============================");
\r
97 Interval next = new Interval();
\r
99 Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);
\r
101 log.debug("Burn Rate: " + burnRate);
\r
103 Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);
\r
105 data.put(data.lastKey().plus(dt), next);
\r
107 log.debug("Dt: " + dt);
\r
109 next.regression = prev.regression.plus(regStep);
\r
111 log.info("Regression: " + next.regression);
\r
113 next.time = prev.time.plus(dt);
\r
115 log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));
\r
117 log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));
\r
119 //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));
\r
120 Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);
\r
122 log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));
\r
124 Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().idealDensity().times(densityRatio)).divide(dt).to(MassFlowRate.UNIT);
\r
126 log.debug("Mass Gen Rate: " + mGenRate);
\r
128 Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().effectiveMolarWeight());
\r
129 Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().idealCombustionTemperature().times(combustionEfficency);
\r
131 Amount<MassFlowRate> mNozzle;
\r
133 Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);
\r
135 //pDiff = Amount.valueOf(.7342, MPA).minus(atmosphereicPressure);
\r
137 log.debug("Pdiff: " + pDiff);
\r
139 Amount<Area> aStar = motor.getNozzle().throatArea();
\r
141 double k = motor.getFuel().getCombustionProduct().ratioOfSpecificHeats();
\r
143 log.debug("K: " + k);
\r
145 double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1))); //Math.pow(2/k+1, (k+1)/(2*(k-1)));
\r
147 log.debug("K-Part: (good)" + kSide);
\r
151 //This unit conversion helps JScience to convert nozzle flow rate to
\r
152 //kg/s a little later on I verified the conversion by hand and
\r
153 //JScience checks it too.
\r
154 specificGasConstant = specificGasConstant.to(
\r
155 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));
\r
157 log.debug("Specific Gas Constant: (good)" + specificGasConstant);
\r
159 Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt();
\r
161 //Unit x = SI.JOULE.divide(SI.KILOGRAM).root(2);
\r
163 //sqrtPart = sqrtPart.times(Amount.valueOf(1, x));
\r
165 log.debug("Square Root Part: " + sqrtPart);
\r
167 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);
\r
169 log.debug("Nozzle Flow: " + mNozzle);
\r
171 log.debug("Nozzle Flow: " + mNozzle.to(MassFlowRate.UNIT));
\r
178 Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);
\r
180 log.debug("Chamber Product rate: " + massStorageRate);
\r
182 next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));
\r
184 log.debug("Chamber Product: " + next.chamberProduct);
\r
186 Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);
\r
188 log.debug("Product Density: " + combustionProductDensity);
\r
190 next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);
\r
192 next.chamberPressure = Amount.valueOf(
\r
193 next.chamberPressure.doubleValue(SI.PASCAL),
\r
196 next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().ratioOfSpecificHeats2Phase());
\r
198 if ( next.chamberPressure.approximates(atmosphereicPressure)){
\r
199 log.info("Pressure at Patm on step " + i);
\r
206 public Amount<Pressure> pressure(Amount<Duration> time){
\r
207 return data.get(time).chamberPressure;
\r
210 public Amount<Force> thrust(Amount<Duration> time){
\r
211 return data.get(time).thrust;
\r
214 public Amount<Dimensionless> kn(Amount<Length> regression){
\r
215 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);
\r
218 public static void main( String args[]) throws Exception{
\r
219 Motor m = new Motor();
\r
220 m.setFuel(new KNSU());
\r
222 CylindricalChamber c = new CylindricalChamber();
\r
223 c.setLength(Amount.valueOf(400, SI.MILLIMETER));
\r
224 c.setID(Amount.valueOf(30, SI.MILLIMETER));
\r
227 CoredCylindricalGrain g = new CoredCylindricalGrain();
\r
228 g.setLength(Amount.valueOf(70, SI.MILLIMETER));
\r
229 g.setOD(Amount.valueOf(30, SI.MILLIMETER));
\r
230 g.setID(Amount.valueOf(10, SI.MILLIMETER));
\r
233 CoredCylindricalGrain g1 = new CoredCylindricalGrain();
\r
234 g1.setLength(Amount.valueOf(70, SI.MILLIMETER));
\r
235 g1.setOD(Amount.valueOf(30, SI.MILLIMETER));
\r
236 g1.setID(Amount.valueOf(18, SI.MILLIMETER));
\r
237 g1.inhibit(false, true, true);
\r
239 CoredCylindricalGrain g2 = new CoredCylindricalGrain();
\r
240 g2.setLength(Amount.valueOf(70, SI.MILLIMETER));
\r
241 g2.setOD(Amount.valueOf(12, SI.MILLIMETER));
\r
242 g2.setID(Amount.valueOf(0, SI.MILLIMETER));
\r
243 g2.inhibit(true, false, true);
\r
245 CompoundGrain cg = new CompoundGrain();
\r
249 //m.setGrain( new MultiGrain(cg, 2) );
\r
251 m.setGrain(new MultiGrain(g,4));
\r
253 //m.setGrain(new ExtrudedGrain());
\r
255 ConvergentDivergentNozzle n = new ConvergentDivergentNozzle();
\r
256 n.setThroatDiameter(Amount.valueOf(5.500, SI.MILLIMETER));
\r
257 n.setExitDiameter(Amount.valueOf(20.87, SI.MILLIMETER));
\r
258 n.setEfficiency(.87);
\r
261 Burn b = new Burn(m);
\r
265 new BurnPanel(b).show();
\r
267 Chart<Duration, Pressure> r = new Chart<Duration, Pressure>(
\r
269 SI.MEGA(SI.PASCAL),
\r
272 r.setDomain(b.data.keySet());
\r
275 Chart<Duration, Force> t = new Chart<Duration, Force>(
\r
280 t.setDomain(b.data.keySet());
\r
283 new GrainPanel( m.getGrain() ).show();*/
\r