Showing posts with label octave. Show all posts
Showing posts with label octave. Show all posts

Tuesday, February 24, 2015

Electron Cont'd (9)

We're now ready to discuss some results from our JSBSim flight dynamics simulation. In this post, we'll talk mainly about how to survive maximum aerodynamic pressure, or max q.
Let's begin by defining the coordinate systems used for our rocket, as well as major forces and velocities involved.


The engines produce a forward thrust that pushes the rocket along the direction of velocity V. Most liquid rocket engines require a little bit of time to develop enough thrust to be able to lift the combined mass of the rocket plus fuel and oxidizer off the ground. During this short period, a hold down mechanism is typically used to keep the rocket stable on the pad. While the rocket is flying within the atmosphere, a drag force due to air resistance (wind) directly opposes the rocket's thrust.

If the rocket's velocity vector is not aligned with the actual pitch angle with the horizon, the resulting angular difference is called the angle of attack alpha. A large alpha may result in a strong perpendicular lifting force due to aerodynamic lift. This is usually bad, because for a long object such as a rocket, excessive perpendicular forces could cause the rocket to bend and result in structural failure. Thus many rockets aim to fly a zero angle of attack and throttle back during flight through max q. This reduces bending forces to a manageable extent.

The following Octave (Matlab) plots of a JSBSim run demonstrate the gravity turn maneuver in terms of pitch angle variation with time (left) and angle of attack (right). A 90-degree pitch angle means "vertical" and 0-degrees is "horizontal". Dynamic pressure appears as peaks in both plots.

Notice how the gravity turn drops the angle of attack to zero just before peak dynamic pressure. We're also throttling back all main engines to 75% at MET+79 and throttling up again at MET+129. To show how bad things would be like without these countermeasures, here is another JSBSim run, but with no engine throttling and the gravity turn approximated by a straight line. Max q is now about 66% greater (500 psi), and the angle of attack nearly hits 7 degrees (as a rule of thumb, most launch vehicles try not to go beyond 4 degrees while passing through max q):

At this point, it is important to note that zero alpha is not a panacea. In fact, zero alpha must not (and does not have to) be maintained throughout the whole flight. This is because sticking to such a flight profile can result in excessive pitchover, and penalties in altitude gain later in flight (see for example: Turner, 2006. Rocket and Spacecraft Propulsion: Principles, Practice and New Developments). So we switch to a lower rate of pitchover as soon as we can (once the rocket passes max q). This allows the rocket to gain altitude more quickly. The following plot illustrates this nicely. Apart from revealing that our rocket has come a long way since our initial experiments and is now able to reach altitudes >1,500,000 ft (>450 km!), notice how the altitude line becomes steeper after passing through the gravity turn phase, indicating that the rocket is able to gain altitude quicker:


Sideways winds gives rise to side forces that causes the rocket to yaw (swivel and spin sideways). Left unchecked, this can result in oscillations that lead to destabilization and tumbling (as we saw in our earlier experiments with OpenRocket). Thus we have also added controllers for holding the yaw rate steady by actuating the first and second stage engines.

In verifying the performance of our yaw and pitch control, it is useful to look at the 3 different components u, v, w of the total velocity V:
  • u: Velocity directly along the axis of the rocket ("forward")
  • v: Sideways velocity
  • w: Vertical velocity
Our aim is, by the end of our launch trajectory, eventually maximize u, and make v (yaw) and w as small as possible—in a circular orbit, u would be a constant orbital velocity and v, w should be zero.

To be continued...

Electron Cont'd (2)

To simulate the atmospheric phase of an Electron launch, I tried using OpenRocket, a Java GUI originally designed for simulating model rocket launches.
However the limitations of using a model rocket simulator to simulate a 10.5-ton orbital launch vehicle quickly became apparent:


In the above screen, I roughly sketched out the two stages of the Electron vehicle using the measurements in my Blender model.
As Electron will use composite materials for its hull, I specified carbon fibre for my OpenRocket model as well.
There is no publicly available data on the thrust curve of the Rutherford engines powering the Electron, so I took the highest thrust motor O8000 available on ThrustCurve.org and scaled up the maximum thrust to match the peak thrust figure quoted for Rutherford on the Rocket Lab web site (146 kN for the 9-engine 1st stage cluster, or 16.2 kN/motor, and 18 kN for the 2nd stage vacuum motor). I also set the burn times for the engines to 114 seconds, which is about how long ESA's Vega burns for its 1st stage, and increased the propellant masses until the wet mass of the rocket was in the ballpark of the quoted figure of 10.5 metric tons.

This gets Electron up to about 1.6 km altitude and nearly Mach 0.4, but without active GNC (guidance, navigation and control), the LV starts to "dance" around and tumble out of control and hits the ground after only 70 seconds in the air. Here is a GNU Octave (Matlab) plot of the simulated flight profile from OpenRocket:


The problem is that all model rocket simulator programs like OpenRocket assume that you are using appendages such as fins to passively stabilize the rocket and cannot simulate active stabilization using a GNC system. Still, it's quite interesting to see how long such a large rocket can stay aloft even without stabilization.

To be able to simulate active control, a different program is required. Apparently JSBSim, which is also used by the open source flight sim FlightGear, allows you to model a rocket GNC system. This is not surprising, as the designer of JSBSim is Jon S. Berndt, who did rocket dynamics engineering work previously for NASA. My plan is to try this next.

Continue reading...

Friday, February 13, 2015

KARI Lunar Orbiter Cont'd (3)

I think that the smoothed velocities are correct (i.e., the unsmoothed velocities with gaps are wrong).

Here is a plot comparing the smoothed and unsmoothed velocities of the spacecraft at launch time:

For both, the velocity profiles can be broken down into 3 main parts:

  1. Launch: Steep increase in velocity
  2. LEO orbit: Constant velocity (circular orbit)
  3. Trans-lunar Injection (TLI): Delta-v increase to push spacecraft towards the Moon

For the unsmoothed profile, the velocity evens out at only ~5.5 km/s, and the delta-v of TLI is only about +2 km/s.
These numbers are too small, especially considering that a typical orbital velocity of a circular LEO should be ~8 km/s.

By comparison, the smoothed profile flattens out at ~7.8 km/s, which is a reasonable velocity for LEO. Also, TLI delta-v is roughly +3 km/s, which is also in the ballpark for typical TLI delta-v's. :D

Continue reading...

KARI Lunar Orbiter Cont'd (2)

I noticed that with the provided trajectory .xyz file, the spacecraft would periodically speed up, slow down, then speed up again, etc.
To illustrate, here is a Octave (Matlab) plot of the spacecraft velocity:


To see why, let's take a look at the .xyz file:
2458330.207  -2956.863152  4765.437104  3027.647441
2458330.207  -2957.134383  4765.280585  3027.653452
2458330.207  -3042.98079   4756.968872  3041.183781
2458330.208  -3244.743734  4682.744325  3039.987276
2458330.209  -3525.597221  4541.81125   3018.195342
2458330.209  -3856.933379  4338.323259  2973.66722
2458330.21   -4210.611914  4078.911213  2905.378869
It seems that the time (first column) was not written with sufficient number of decimal places, causing several rows to have the same time.
So I tried to simply delete any rows with duplicated times. But this didn't work, because it caused the [b]distance[/b] between the same timestamps to vary.
In the above example, deleting duplicate .207 and .209 entries would cause the distance to jump suddenly between .207 and .208, and .209 to .21.
Since velocity=distance/time, this causes the velocity to jump suddenly also.

The solution was to recalculate all the timestamps. I simply took the total time interval between the first row of the .xyz file and the last row, and divided by the total number of .xyz entries minus 1. Then I made all the timestamps increment by this constant time interval.

This is the resulting plot:


As you can see,  there are no more gaps and moreover the velocity curve is smooth. The motion of the spacecraft in celestia.Sci is also much smoother. However, the velocities are also higher than the original; I'm still investigating why.

Continue reading...

Microlensing and exoplanets

Finally, let’s discuss how to represent amplification due to lensing. Previously we’ve seen that amplification due to microlensing is an important technique used to detect exoplanets (other techniques include radial velocity, transit, etc). But how can this be possible, since planets are so much smaller in mass than stars and galaxy clusters? Actually, even stars by themselves amplify light by focusing it, but the effect is very small. This GNU Octave plot illustrates the amplification factor around our Sun:


Things get interesting however when a planet orbits close to the star. While the planet itself exerts an even smaller gravitational influence than its star, when the two influences combine they can produce extreme, discontinuous jumps in brightness called caustics. Caustics can be seen as bright wavy lines in a pool of water, as illustrated in the excellent Optics Picture of the Day website: http://www.atoptics.co.uk/fz535.htm.

The following Octave light curve plot shows the amplification that is observed when a single exoplanet passes near the Einstein radius of its star. The Einstein radius is defined as \(2R_s\frac{D_{ds}}{D_s}\) where \(R_s\) is the Schwarzschild radius of the lens, a measure of how much the lens curves spacetime (and also the radius of the event horizon in a black hole). The Einstein radius can be thought of as the characteristic size of the lens, and is also the size by which an Einstein ring would appear around the lens.


Here is an actual light curve observed for an extrasolar system 4.1 kpc from Earth called OGLE-2012-BLG-0026.

Han, C. et al., 2013. The Second Multiple-planet System Discovered by Microlensing: OGLE-2012-BLG-0026Lb, c - A Pair of Jovian Planets beyond the Snow Line. The Astrophysical Journal, 762(2), p.L28.
I've implemented a real-time plotting feature in celestia.Sci using the Qwt framework. Here it is in action:


First to briefly describe how this was done, I had to create an add-on for the extrasolar system as stock celestia.Sci doesn't have it. Then I chose View > Plot, selected the star OGLE-2012-BLG-0026L, lined it up with a background star which I arbitrarily named OGLE-2012-BLG-0026L-SRC, and hit Refresh in the plot panel. The microlensing code then starts varying the impact parameter across the lens to produce the time-varying light curve.

The light curve of OGLE-2012-BLG-0026 plotted by celestia.Sci shows several differences such as higher peaks. This however could be explained by the fact that we are approximating lenses and source objects as points of zero area;  this results in singularities in the calculated amplification factors. Real microlenses have non-zero area, and this smooths out the singularities. However, qualitatively we can identify similarities between the simulated and real light curves, such as the presence of multiple peaks. Future improvements could try to increase the realism by taking into account the disc sizes of stars and planets.