X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=src%2Fcom%2Fbillkuker%2Frocketry%2Fmotorsim%2FBurn.java;h=0416d1ec0892a087710627ceb61555854534f99c;hb=6a3d94f706ba306d48d96cbad16c3bd3004e2c6e;hp=39933d2d91e05d0cad36663569449d1bc834e96c;hpb=98b8e334f9300f6696abf74a0c03be748bbcc37f;p=sw%2Fmotorsim diff --git a/src/com/billkuker/rocketry/motorsim/Burn.java b/src/com/billkuker/rocketry/motorsim/Burn.java index 39933d2..0416d1e 100644 --- a/src/com/billkuker/rocketry/motorsim/Burn.java +++ b/src/com/billkuker/rocketry/motorsim/Burn.java @@ -2,6 +2,9 @@ package com.billkuker.rocketry.motorsim; +import java.util.Date; +import java.util.HashSet; +import java.util.Set; import java.util.SortedMap; import java.util.TreeMap; @@ -13,6 +16,7 @@ import javax.measure.quantity.Length; import javax.measure.quantity.Mass; import javax.measure.quantity.MassFlowRate; import javax.measure.quantity.Pressure; +import javax.measure.quantity.Quantity; import javax.measure.quantity.Temperature; import javax.measure.quantity.Velocity; import javax.measure.quantity.Volume; @@ -23,40 +27,103 @@ import org.apache.log4j.Logger; import org.jscience.physics.amount.Amount; import org.jscience.physics.amount.Constants; -import com.billkuker.rocketry.motorsim.fuel.KNSB; -import com.billkuker.rocketry.motorsim.fuel.KNSU; -import com.billkuker.rocketry.motorsim.grain.CompoundGrain; -import com.billkuker.rocketry.motorsim.grain.CoredCylindricalGrain; -import com.billkuker.rocketry.motorsim.grain.MultiGrain; -import com.billkuker.rocketry.motorsim.visual.BurnPanel; +import com.billkuker.rocketry.motorsim.Validating.ValidationException; public class Burn { - private static Logger log = Logger.getLogger(Burn.class); - protected final Motor motor; - private static final Amount atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL); + private static final BurnSettings settings = new BurnSettings(); + public static final BurnSettings getBurnSettings(){ + return settings; + } + /** + * A class representing all the settigns one can change on a burn + * @author bkuker + */ + public static class BurnSettings { + private BurnSettings(){}; + + public enum BurnVolumeMethod { + DeltaVolume, + SurfaceTimesRegression; + } + + private BurnVolumeMethod volumeMethod = BurnVolumeMethod.SurfaceTimesRegression; + private double regStepIncreaseFactor = 1.01; + private double regStepDecreaseFactor = .5; + private Amount chamberPressureMaxDelta = Amount.valueOf(.5, SI.MEGA(SI.PASCAL)); + private Amount endPressure = Amount.valueOf(.1, RocketScience.PSI); + private Amount atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL); + + public void setVolumeMethod(BurnVolumeMethod volumeMethod) { + this.volumeMethod = volumeMethod; + } + public BurnVolumeMethod getVolumeMethod() { + return volumeMethod; + } + public double getRegStepIncreaseFactor() { + return regStepIncreaseFactor; + } + public void setRegStepIncreaseFactor(double regStepIncreaseFactor) { + this.regStepIncreaseFactor = regStepIncreaseFactor; + } + public double getRegStepDecreaseFactor() { + return regStepDecreaseFactor; + } + public void setRegStepDecreaseFactor(double regStepDecreaseFactor) { + this.regStepDecreaseFactor = regStepDecreaseFactor; + } + public Amount getChamberPressureMaxDelta() { + return chamberPressureMaxDelta; + } + public void setChamberPressureMaxDelta(Amount chamberPressureMaxDelta) { + this.chamberPressureMaxDelta = chamberPressureMaxDelta; + } + public Amount getEndPressure() { + return endPressure; + } + public void setEndPressure(Amount endPressure) { + this.endPressure = endPressure; + } + public Amount getAtmosphereicPressure() { + return atmosphereicPressure; + } + public void setAtmosphereicPressure(Amount atmosphereicPressure) { + this.atmosphereicPressure = atmosphereicPressure; + } + } - private static double combustionEfficency = 0.97; + protected final Motor motor; - private static double densityRatio = 0.96; + private boolean burning = false; + private boolean done = false; + public interface BurnProgressListener{ + public void setProgress(float p); + public void burnComplete(); + } + + private Set bpls = new HashSet(); + public class Interval{ - Amount time; + public Amount time; + public Amount dt; public Amount regression; public Amount chamberPressure; Amount chamberProduct; - Amount thrust; + public Amount thrust; public String toString(){ - return time + " " + regression + " " + chamberPressure + " " + chamberProduct; + return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct; } } protected SortedMap,Interval> data = new TreeMap, Interval>(); public SortedMap,Interval> getData(){ + if ( !done ) + throw new IllegalStateException("Burn not complete!"); return data; } @@ -65,225 +132,211 @@ public class Burn { } public Amount burnTime(){ - return data.lastKey(); + return getData().lastKey(); } public Burn(Motor m){ + try { + m.validate(); + } catch (ValidationException e) { + throw new IllegalArgumentException("Invalid Motor: " + e.getMessage()); + } motor = m; - burn(); } - private void burn(){ + public void addBurnProgressListener( BurnProgressListener bpl ){ + bpls.add(bpl); + if ( done ) + bpl.burnComplete(); + } + + public void burn(){ + synchronized(this){ + if ( burning ) + throw new IllegalStateException("Already burning!"); + burning = true; + } log.info("Starting burn..."); + for (BurnProgressListener bpl : bpls ){ + bpl.setProgress(0); + } - Amount regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER); - - //if ( motor.getGrain() instanceof Grain.DiscreteRegression ) - //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep(); + int endPressureSteps = 0; + long start = new Date().getTime(); + Amount regStep = Amount.valueOf(0.01, SI.MILLIMETER); + Interval initial = new Interval(); initial.time = Amount.valueOf(0, SI.SECOND); + initial.dt = Amount.valueOf(0, SI.SECOND); initial.regression = Amount.valueOf(0, SI.MILLIMETER); - initial.chamberPressure = atmosphereicPressure; + initial.chamberPressure = settings.getAtmosphereicPressure(); initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM); initial.thrust = Amount.valueOf(0, SI.NEWTON); data.put(Amount.valueOf(0, SI.SECOND), initial); - for ( int i = 0; i < 5000; i++ ){ - + step: + for ( int i = 0; i < 5000; i++ ) { + assert(positive(regStep)); + regStep = regStep.times(settings.getRegStepIncreaseFactor()); + Interval prev = data.get(data.lastKey()); log.debug(prev); log.debug("Step " + i + " =============================="); Interval next = new Interval(); Amount burnRate = motor.getFuel().burnRate(prev.chamberPressure); + assert(positive(burnRate)); log.debug("Burn Rate: " + burnRate); Amount dt = regStep.divide(burnRate).to(Duration.UNIT); + assert(positive(dt)); + next.dt = dt; - data.put(data.lastKey().plus(dt), next); + log.debug("Dt: " + dt); next.regression = prev.regression.plus(regStep); + assert(positive(next.regression)); - log.info("Regression: " + next.regression); + log.debug("Regression: " + next.regression); + + //Update BurnProgressListeners + Amount a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT); + for (BurnProgressListener bpl : bpls ){ + bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT)); + } + next.time = prev.time.plus(dt); - log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3))); + //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3))); - log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3))); + //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3))); - //TODO Amount volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression)); - Amount volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT); + Amount volumeBurnt; + if ( settings.getVolumeMethod() == BurnSettings.BurnVolumeMethod.DeltaVolume ){ + volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression)); + } else { + volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT); + } + assert(positive(volumeBurnt)); + //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3))); - log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3))); + Amount mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT); + assert(positive(mGenRate)); - Amount mGenRate = volumeBurnt.times(motor.getFuel().idealDensity().times(densityRatio)).divide(dt).to(MassFlowRate.UNIT); + //log.debug("Mass Gen Rate: " + mGenRate); - log.debug("Mass Gen Rate: " + mGenRate); + //Calculate specific gas constant + Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight()); + //This unit conversion helps JScience to convert nozzle flow rate to + //kg/s a little later on I verified the conversion by hand and + //JScience checks it too. + specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant); - Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().effectiveMolarWeight()); - Amount chamberTemp = motor.getFuel().getCombustionProduct().idealCombustionTemperature().times(combustionEfficency); + //Calculate chamber temperature + Amount chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency()); Amount mNozzle; { - Amount pDiff = prev.chamberPressure.minus(atmosphereicPressure); - - //pDiff = Amount.valueOf(.7342, MPA).minus(atmosphereicPressure); - - log.debug("Pdiff: " + pDiff); - + Amount pDiff = prev.chamberPressure.minus(settings.getAtmosphereicPressure()); + //log.debug("Pdiff: " + pDiff); Amount aStar = motor.getNozzle().throatArea(); - - double k = motor.getFuel().getCombustionProduct().ratioOfSpecificHeats(); - - log.debug("K: " + k); - - 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))); - - log.debug("K-Part: (good)" + kSide); - - - - //This unit conversion helps JScience to convert nozzle flow rate to - //kg/s a little later on I verified the conversion by hand and - //JScience checks it too. - specificGasConstant = specificGasConstant.to( - SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN))); - - log.debug("Specific Gas Constant: (good)" + specificGasConstant); - - Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt(); - - //Unit x = SI.JOULE.divide(SI.KILOGRAM).root(2); - - //sqrtPart = sqrtPart.times(Amount.valueOf(1, x)); - - log.debug("Square Root Part: " + sqrtPart); - + double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats(); + double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1))); + Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt(); mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT); - - log.debug("Nozzle Flow: " + mNozzle); - - log.debug("Nozzle Flow: " + mNozzle.to(MassFlowRate.UNIT)); - - - - + //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT)); } + assert(positive(mNozzle)); Amount massStorageRate = mGenRate.minus(mNozzle); - log.debug("Chamber Product rate: " + massStorageRate); - + //log.debug("Mass Storage Rate: " + massStorageRate); + next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt)); - log.debug("Chamber Product: " + next.chamberProduct); + //Product can not go negative! + if ( !positive(next.chamberProduct) ){ + log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!"); + regStep = regStep.times(settings.getRegStepDecreaseFactor()); + continue step; + } + assert(positive(next.chamberProduct)); + if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) ) + next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM); + + //log.debug("Chamber Product: " + next.chamberProduct); Amount combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT); log.debug("Product Density: " + combustionProductDensity); - next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT); + next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(settings.getAtmosphereicPressure()).to(Pressure.UNIT); + assert(positive(next.chamberPressure)); next.chamberPressure = Amount.valueOf( next.chamberPressure.doubleValue(SI.PASCAL), SI.PASCAL); - next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().ratioOfSpecificHeats2Phase()); + Amount dp = next.chamberPressure.minus(prev.chamberPressure); + if ( dp.abs().isGreaterThan(settings.getChamberPressureMaxDelta())){ + log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!"); + regStep = regStep.times(settings.getRegStepDecreaseFactor()); + continue step; + } + + next.thrust = motor.getNozzle().thrust(next.chamberPressure, settings.getAtmosphereicPressure(), settings.getAtmosphereicPressure(), motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase()); + assert(positive(next.thrust)); - if ( next.chamberPressure.approximates(atmosphereicPressure)){ - log.info("Pressure at Patm on step " + i); - break; + if ( i > 100 && next.chamberPressure.minus(settings.getAtmosphereicPressure()).abs().isLessThan(settings.getEndPressure())){ + log.info("Pressure at ~Patm on step " + i); + endPressureSteps++; + if ( endPressureSteps > 5 ) + break; } + + data.put(data.lastKey().plus(dt), next); } - + + long time = new Date().getTime() - start; + log.info("Burn took " + time + " millis."); + done = true; + for (BurnProgressListener bpl : bpls ){ + bpl.burnComplete(); + } + } + + @SuppressWarnings("unchecked") + /* + * This converts the units of this constant to something JScience is able + * to work from. This conversion is unchecked at compile time, but + * JScience keeps me honest at runtime. + */ + private Amount convertSpecificGasConstantUnits(Amount a){ + return a.to( + SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN))); } public Amount pressure(Amount time){ - return data.get(time).chamberPressure; + return getData().get(time).chamberPressure; } public Amount thrust(Amount time){ - return data.get(time).thrust; + return getData().get(time).thrust; } public Amount kn(Amount regression){ return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT); } - public static void main( String args[]) throws Exception{ - Motor m = new Motor(); - m.setFuel(new KNSB()); - - CylindricalChamber c = new CylindricalChamber(); - c.setLength(Amount.valueOf(300, SI.MILLIMETER)); - c.setID(Amount.valueOf(30, SI.MILLIMETER)); - m.setChamber(c); - - CoredCylindricalGrain g = new CoredCylindricalGrain(); - g.setLength(Amount.valueOf(35, SI.MILLIMETER)); - g.setOD(Amount.valueOf(30, SI.MILLIMETER)); - g.setID(Amount.valueOf(10, SI.MILLIMETER)); - m.setGrain(g); - - CoredCylindricalGrain g1 = new CoredCylindricalGrain(); - g1.setLength(Amount.valueOf(70, SI.MILLIMETER)); - g1.setOD(Amount.valueOf(30, SI.MILLIMETER)); - g1.setID(Amount.valueOf(18, SI.MILLIMETER)); - g1.inhibit(false, true, true); - - CoredCylindricalGrain g2 = new CoredCylindricalGrain(); - g2.setLength(Amount.valueOf(70, SI.MILLIMETER)); - g2.setOD(Amount.valueOf(12, SI.MILLIMETER)); - g2.setID(Amount.valueOf(0, SI.MILLIMETER)); - g2.inhibit(true, false, true); - - CompoundGrain cg = new CompoundGrain(); - cg.add(g1); - cg.add(g2); - - //m.setGrain( new MultiGrain(cg, 2) ); - - //g.setAftEndInhibited(true); - //g.setForeEndInhibited(true); - m.setGrain(new MultiGrain(g,3)); - - //m.setGrain(new ExtrudedGrain()); - - ConvergentDivergentNozzle n = new ConvergentDivergentNozzle(); - n.setThroatDiameter(Amount.valueOf(5.500, SI.MILLIMETER)); - n.setExitDiameter(Amount.valueOf(20.87, SI.MILLIMETER)); - n.setEfficiency(.87); - m.setNozzle(n); - - Burn b = new Burn(m); - - b.burn(); - - new BurnPanel(b).show(); - /* - Chart r = new Chart( - SI.SECOND, - SI.MEGA(SI.PASCAL), - b, - "pressure"); - r.setDomain(b.data.keySet()); - r.show(); - - Chart t = new Chart( - SI.SECOND, - SI.NEWTON, - b, - "thrust"); - t.setDomain(b.data.keySet()); - t.show(); - - new GrainPanel( m.getGrain() ).show();*/ - + + private boolean positive(Amount a){ + return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a))); } + }