Merge pull request #2570 from SiboVG/improve-beta

Clamp beta factor (compressibility correction) in transonic region + apply it in supersonic region

core/src/main/java/info/openrocket/core/aerodynamics/FlightConditions.java

@@ -22,6 +22,7 @@ import info.openrocket.core.util.ModID;
  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  */
 public class FlightConditions implements Cloneable, ChangeSource, Monitorable {
+	private static final double MIN_BETA = 0.25;
 
 	private List<EventListener> listenerList = new ArrayList<>();
 	private EventObject event = new EventObject(this);
@@ -55,7 +56,7 @@ public class FlightConditions implements Cloneable, ChangeSource, Monitorable {
 	 * Sqrt(1 - M^2) for M<1
 	 * Sqrt(M^2 - 1) for M>1
 	 */
-	private double beta = MathUtil.safeSqrt(1 - mach * mach);
+	private double beta = calculateBeta(mach);
 
 	/** Current roll rate. */
 	private double rollRate = 0;
@@ -243,10 +244,7 @@ public class FlightConditions implements Cloneable, ChangeSource, Monitorable {
 			return;
 
 		this.mach = mach;
-		if (mach < 1)
+		this.beta = calculateBeta(mach);
-			this.beta = MathUtil.safeSqrt(1 - mach * mach);
-		else
-			this.beta = MathUtil.safeSqrt(mach * mach - 1);
 
 		fireChangeEvent();
 	}
@@ -288,6 +286,19 @@ public class FlightConditions implements Cloneable, ChangeSource, Monitorable {
 		return beta;
 	}
 
+	/**
+	 * Calculate the beta value (compressibility factor/Prandtl-Glauert correction factor) for the given Mach number.
+	 * @param mach the Mach number.
+	 * @return the beta value.
+	 */
+	private static double calculateBeta(double mach) {
+		if (mach < 1) {
+			return MathUtil.max(MIN_BETA, MathUtil.safeSqrt(1 - mach * mach));
+		} else {
+			return MathUtil.max(MIN_BETA, MathUtil.safeSqrt(mach * mach - 1));
+		}
+	}
+
 	/**
 	 * @return the current roll rate.
 	 */

core/src/main/java/info/openrocket/core/aerodynamics/barrowman/FinSetCalc.java

@@ -162,9 +162,6 @@ public class FinSetCalc extends RocketComponentCalc {
 		// TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25
 		// Calculate CP position
 		double x = macLead + calculateCPPos(conditions) * macLength;
-		//		logger.debug("Component macLead = {}", macLead);
-		//		logger.debug("Component macLength = {}", macLength);
-		//		logger.debug("Component x = {}", x);
 		
 		
 		// Calculate roll forces, reduce forcing above stall angle
@@ -176,17 +173,12 @@ public class FinSetCalc extends RocketComponentCalc {
 		forces.setCrollForce((macSpan + r) * cna1 * (1 + tau) * cantAngle / conditions.getRefLength());
 		
 		if (conditions.getAOA() > STALL_ANGLE) {
-			//			System.out.println("Fin stalling in roll");
 			forces.setCrollForce(forces.getCrollForce() * MathUtil.clamp(
 					1 - (conditions.getAOA() - STALL_ANGLE) / (STALL_ANGLE / 2), 0, 1));
 		}
 		forces.setCrollDamp(calculateDampingMoment(conditions));
 		forces.setCroll(forces.getCrollForce() - forces.getCrollDamp());
 		
-		//		System.out.printf(component.getName() + ":  roll rate:%.3f  force:%.3f  damp:%.3f  " +
-		//				"total:%.3f\n",
-		//				conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll);
-		
 		forces.setCNa(cna);
 		forces.setCN(cna * MathUtil.min(conditions.getAOA(), STALL_ANGLE));
 		forces.setCP(new Coordinate(x, 0, 0, cna));
@@ -349,7 +341,6 @@ public class FinSetCalc extends RocketComponentCalc {
 			double y = i * dy;
 			
 			macLength += length * length;
-			//logger.debug("macLength = {}, length = {}, i = {}", macLength, length, i);
 			macSpan += y * length;
 			macLead += chordLead[i] * length;
 			area += length;
@@ -408,10 +399,6 @@ public class FinSetCalc extends RocketComponentCalc {
 		K1 = new LinearInterpolator(x, k1);
 		K2 = new LinearInterpolator(x, k2);
 		K3 = new LinearInterpolator(x, k3);
-		
-		//		System.out.println("K1[m="+CNA_SUPERSONIC+"] = "+k1[0]);
-		//		System.out.println("K2[m="+CNA_SUPERSONIC+"] = "+k2[0]);
-		//		System.out.println("K3[m="+CNA_SUPERSONIC+"] = "+k3[0]);
 	}
 	
 	protected double calculateFinCNa1(FlightConditions conditions) {
@@ -445,8 +432,6 @@ public class FinSetCalc extends RocketComponentCalc {
 				K3.getValue(CNA_SUPERSONIC) * pow2(alpha)) / ref;
 		superD = -finArea / ref * 2 * CNA_SUPERSONIC / CNA_SUPERSONIC_B;
 		
-		//		System.out.println("subV="+subV+" superV="+superV+" subD="+subD+" superD="+superD);
-		
 		return cnaInterpolator.interpolate(mach, subV, superV, subD, superD, 0);
 	}
 	
@@ -473,24 +458,16 @@ public class FinSetCalc extends RocketComponentCalc {
 			}
 			sum = sum * (span / DIVISIONS) * 2 * Math.PI / conditions.getBeta() /
 					(conditions.getRefArea() * conditions.getRefLength());
-			
+
-			//			System.out.println("SPECIAL: " + 
-			//					(MathUtil.sign(rollRate) * sum));
 			return MathUtil.sign(rollRate) * sum;
 		}
 		
 		if (mach <= CNA_SUBSONIC) {
-			//			System.out.println("BASIC:   "+
-			//					(2*Math.PI * rollRate * rollSum / 
-			//			(conditions.getRefArea() * conditions.getRefLength() * 
-			//					conditions.getVelocity() * conditions.getBeta())));
-			
 			return 2 * Math.PI * rollRate * rollSum /
 					(conditions.getRefArea() * conditions.getRefLength() *
 							conditions.getVelocity() * conditions.getBeta());
 		}
 		if (mach >= CNA_SUPERSONIC) {
-			
 			double vel = conditions.getVelocity();
 			double k1 = K1.getValue(mach);
 			double k2 = K2.getValue(mach);
@@ -511,7 +488,6 @@ public class FinSetCalc extends RocketComponentCalc {
 		}
 		
 		// Transonic, do linear interpolation
-		
 		FlightConditions cond = conditions.clone();
 		cond.setMach(CNA_SUBSONIC - 0.01);
 		double subsonic = calculateDampingMoment(cond);
@@ -532,7 +508,7 @@ public class FinSetCalc extends RocketComponentCalc {
 	 */
 	private double calculateCPPos(FlightConditions cond) {
 		double m = cond.getMach();
-		//		logger.debug("m = {} ", m);
+
 		if (m <= 0.5) {
 			// At subsonic speeds CP at quarter chord
 			return 0.25;
@@ -551,7 +527,7 @@ public class FinSetCalc extends RocketComponentCalc {
 			val += v * x;
 			x *= m;
 		}
-		//		logger.debug("val = {}", val);
+
 		return val;
 	}
 	

core/src/test/java/info/openrocket/core/aerodynamics/FlightConditionsTest.java

@@ -0,0 +1,213 @@
+package info.openrocket.core.aerodynamics;
+
+import static org.junit.jupiter.api.Assertions.*;
+
+import com.google.inject.Guice;
+import com.google.inject.Injector;
+import com.google.inject.Module;
+import info.openrocket.core.ServicesForTesting;
+import info.openrocket.core.document.OpenRocketDocument;
+import info.openrocket.core.document.OpenRocketDocumentFactory;
+import info.openrocket.core.plugin.PluginModule;
+import info.openrocket.core.rocketcomponent.BodyTube;
+import info.openrocket.core.rocketcomponent.Rocket;
+import info.openrocket.core.startup.Application;
+import org.junit.jupiter.api.BeforeEach;
+import org.junit.jupiter.api.Test;
+
+import info.openrocket.core.models.atmosphere.AtmosphericConditions;
+import info.openrocket.core.rocketcomponent.FlightConfiguration;
+import info.openrocket.core.util.Coordinate;
+
+class FlightConditionsTest {
+	private FlightConditions conditions;
+	private static final double EPSILON = 1e-6;
+
+	private Rocket rocket;
+
+	@BeforeEach
+	void setUp() {
+		com.google.inject.Module applicationModule = new ServicesForTesting();
+		Module pluginModule = new PluginModule();
+
+		Injector injector = Guice.createInjector(applicationModule, pluginModule);
+		Application.setInjector(injector);
+
+		OpenRocketDocument document = OpenRocketDocumentFactory.createNewRocket();
+		this.rocket = document.getRocket();
+		BodyTube bodyTube = new BodyTube();
+		bodyTube.setLength(1.0);
+		bodyTube.setOuterRadius(0.05);
+		this.rocket.getChild(0).addChild(bodyTube);
+
+		FlightConfiguration mockConfig = new FlightConfiguration(rocket);
+		this.conditions = new FlightConditions(mockConfig);
+	}
+
+	@Test
+	void testSetAndGetRefLength() {
+		double expectedLength = ((BodyTube) rocket.getChild(0).getChild(0)).getOuterRadius() * 2;
+		assertEquals(expectedLength, conditions.getRefLength(), EPSILON);
+
+		conditions.setRefLength(2.0);
+		assertEquals(2.0, conditions.getRefLength(), EPSILON);
+		assertEquals(Math.PI, conditions.getRefArea(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetRefArea() {
+		// Get the actual reference area from FlightConditions
+		double actualRefArea = conditions.getRefArea();
+
+		// Test that setting this area results in the same value
+		conditions.setRefArea(actualRefArea);
+		assertEquals(actualRefArea, conditions.getRefArea(), EPSILON);
+
+		// Test that setting a new area works correctly
+		double newArea = 4.0;
+		conditions.setRefArea(newArea);
+		assertEquals(newArea, conditions.getRefArea(), EPSILON);
+		assertEquals(Math.sqrt(newArea / Math.PI) * 2, conditions.getRefLength(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetAOA() {
+		conditions.setAOA(Math.PI / 4);
+		assertEquals(Math.PI / 4, conditions.getAOA(), EPSILON);
+		assertEquals(Math.sin(Math.PI / 4), conditions.getSinAOA(), EPSILON);
+		assertEquals(Math.sin(Math.PI / 4) / (Math.PI / 4), conditions.getSincAOA(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetTheta() {
+		conditions.setTheta(Math.PI / 3);
+		assertEquals(Math.PI / 3, conditions.getTheta(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetMach() {
+		conditions.setMach(0.2);
+		assertEquals(0.2, conditions.getMach(), EPSILON);
+		assertEquals(0.9797958971, conditions.getBeta(), EPSILON);
+
+		conditions.setMach(0.8);
+		assertEquals(0.8, conditions.getMach(), EPSILON);
+		assertEquals(0.6, conditions.getBeta(), EPSILON);
+
+		conditions.setMach(0.9999999999);
+		assertEquals(0.9999999999, conditions.getMach(), EPSILON);
+		assertEquals(0.25, conditions.getBeta(), EPSILON);
+
+		conditions.setMach(1.00000000001);
+		assertEquals(1.00000000001, conditions.getMach(), EPSILON);
+		assertEquals(0.25, conditions.getBeta(), EPSILON);
+
+		conditions.setMach(1.3);
+		assertEquals(1.3, conditions.getMach(), EPSILON);
+		assertEquals(0.8306623863, conditions.getBeta(), EPSILON);
+
+		conditions.setMach(3);
+		assertEquals(3, conditions.getMach(), EPSILON);
+		assertEquals(2.8284271247, conditions.getBeta(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetVelocity() {
+		AtmosphericConditions atm = new AtmosphericConditions();
+		conditions.setAtmosphericConditions(atm);
+
+		double expectedMachSpeed = atm.getMachSpeed();
+		conditions.setVelocity(expectedMachSpeed / 2);
+
+		assertEquals(0.5, conditions.getMach(), EPSILON);
+		assertEquals(expectedMachSpeed / 2, conditions.getVelocity(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetRollRate() {
+		conditions.setRollRate(5.0);
+		assertEquals(5.0, conditions.getRollRate(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetPitchRate() {
+		conditions.setPitchRate(2.5);
+		assertEquals(2.5, conditions.getPitchRate(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetYawRate() {
+		conditions.setYawRate(1.5);
+		assertEquals(1.5, conditions.getYawRate(), EPSILON);
+	}
+
+	@Test
+	void testSetAndGetPitchCenter() {
+		Coordinate center = new Coordinate(1.0, 2.0, 3.0);
+		conditions.setPitchCenter(center);
+		assertEquals(center, conditions.getPitchCenter());
+	}
+
+	@Test
+	void testClone() {
+		conditions.setAOA(Math.PI / 6);
+		conditions.setMach(0.7);
+		conditions.setRollRate(3.0);
+		AtmosphericConditions atm = new AtmosphericConditions(280, 90000);
+		conditions.setAtmosphericConditions(atm);
+
+		FlightConditions cloned = conditions.clone();
+
+		assertNotSame(conditions, cloned);
+		assertEquals(conditions.getAOA(), cloned.getAOA(), EPSILON);
+		assertEquals(conditions.getMach(), cloned.getMach(), EPSILON);
+		assertEquals(conditions.getRollRate(), cloned.getRollRate(), EPSILON);
+		assertEquals(conditions.getAtmosphericConditions().getTemperature(),
+				cloned.getAtmosphericConditions().getTemperature(), EPSILON);
+		assertEquals(conditions.getAtmosphericConditions().getPressure(),
+				cloned.getAtmosphericConditions().getPressure(), EPSILON);
+	}
+
+	@Test
+	void testEquals() {
+		FlightConditions conditions1 = new FlightConditions(null);
+		FlightConditions conditions2 = new FlightConditions(null);
+
+		conditions1.setAOA(Math.PI / 6);
+		conditions1.setMach(0.7);
+		conditions2.setAOA(Math.PI / 6);
+		conditions2.setMach(0.7);
+
+		assertTrue(conditions1.equals(conditions2));
+
+		conditions2.setMach(0.8);
+		assertFalse(conditions1.equals(conditions2));
+	}
+
+	@Test
+	void testSetAndGetAtmosphericConditions() {
+		AtmosphericConditions atm = new AtmosphericConditions(280, 90000);
+		conditions.setAtmosphericConditions(atm);
+
+		assertEquals(280, conditions.getAtmosphericConditions().getTemperature(), EPSILON);
+		assertEquals(90000, conditions.getAtmosphericConditions().getPressure(), EPSILON);
+	}
+
+	@Test
+	void testGetVelocityWithChangedAtmosphere() {
+		AtmosphericConditions atm = new AtmosphericConditions(280, 90000);
+		conditions.setAtmosphericConditions(atm);
+		conditions.setMach(0.5);
+
+		double expectedVelocity = 0.5 * atm.getMachSpeed();
+		assertEquals(expectedVelocity, conditions.getVelocity(), EPSILON);
+
+		// Change atmospheric conditions
+		atm.setTemperature(300);
+		conditions.setAtmosphericConditions(atm);
+
+		// Velocity should change with new atmospheric conditions
+		expectedVelocity = 0.5 * atm.getMachSpeed();
+		assertEquals(expectedVelocity, conditions.getVelocity(), EPSILON);
+	}
+}