Here are some equations that prove useful in modeling the VBP's behavior.

Approximating Atmospheric Conditions

This approximation describes atmospheric pressure:

P(h) = Po e^{-r h}

where:

P(h) = Mean Atmospheric Pressure at Altitude h

h = Altitude Above Sea Level

Po = Mean Atmospheric Pressure at Sea Level = 1013 mb

r = Atmospheric Density Decay Rate = 0.00014 ¹/_m

Note that atmospheric pressure decays exponentially with altitude. The actual empirical atmospheric decay rate is given by an exponential of the form P(h) = Po e^{a h ^2 + b h + c} , but this approximation holds to within +/-5mb up to 25000km.

Air Density behaves similarly, following this approximation:

r(h) = ro e^{-r h}

where:

r(h) = Atmospheric Density at Altitude h

ro = Mean Atmospheric Density at Sea Level = 1.206 ^kg/_m3

Air temperature is more variable with location and time than either air density or pressure. However, on the mean it follows this approximation:

T(h) = To + mtroposphere h, h < h thermocline

= Tthermocline + mstratosphere h, h > hthermocline

where:

T(h) = Mean Atmospheric Temperature at Altitude h

To = Mean Atmospheric Temperature at Sea Level

Tthermocline = Mean Atmospheric Temperature at Thermocline

mtroposphere = Rate of Temperature Decrease in Troposphere

mstratosphere = Rate of Temperature Decrease in Stratosphere

hthermocline = Mean Altitude of Thermocline

All of these variables change with location and time of year, but some useful averages are To = 288 K, hthermocline = 16500 m, Tthermocline = 183 K, mtroposphere = -0.00636 ^K/_m, mstratopshere = +0.002 ^K/_m. Note that tropospheric temperatures generally decrease with altitude and stratospheric temperatures generally increase with altitude, making the thermocline -- the inflection point in this equation -- the coldest point in the lower atmospheric regions.

Wind speed is nearly unpredictable, and must be determined empirically using atmospheric soundings for each available data point. However, maximum wind speed expected at altitude can be approximated as:

vmax(h) = a vmax,o e ^{+0.5 r h}

where:

vmax(h) = Maximum Expected Windspeed at Altitude h

vmax,o = Maximum Expected Windspeed at Sea Level

a = Scaling Constant (varies with location)

Note that this last equation is nothing more than a guess regarding the worst case scenario, and is no reflection of what upper level winds will actually be doing at any given time. In fact, this smooth distribution is NEVER found in nature. It is mainly useful for computing maximum loads.

Approximating Virtual Beanstalk Behavior -- The Platform

The Universal Gas Equation describes the lift gas in the lift cells under static conditions:

Pgas Vgas = ngas R Tgas

where:

Pgas = Gas Pressure Inside Balloon

Vgas = Volume of Gas Inside Balloon

ngas = Number of Moles of Lift Gas

R = Universal Gas Constant = 8.314 ^{N m}/_{K mol}

Tgas = Temperature of Lift Gas

The mass of the gas inside the balloon is:

mgas = ngas Wgas

where

mgas = Mass of Lift Gas Inside Balloon

Wgas = Molecular Weight of Lift Gas Inside Balloon = 2 ^g/_mol for Hydrogen

Since the air displaced by the balloon also obeys the Universal Gas Equation, we can derive a formula for the difference between the mass of the lift gas and the mass of the air it displaces. This is the equation for balloon lift:

Flift = g Vgas / R ( Pair Wair / Tair - Pgas Wgas / Tgas ) j

where:

Flift = Vector Lift Force Exerted by Platform Lift Cells

g = Gravitational Acceleration Constant = 9.81 m/s2

Pair = Ambient Air Pressure

Wair = Mean Molecular Weight of Air = 28.6 ^g/_mol

Tair = T(h)

j = y-axis Unit Vector

The derivation of this equation assumes Vgas = Vair = Vballoon. Also, in practice, Pair = Pgas for a zero-pressure balloon. However, Tair <> Tgas. Whenever the balloon changes altitude, it must cool or warm up to the temperature at its new altitude. The balloon changes temperature by two mechanisms: Newtonian Cooling and Adiabatic Cooling.

Newtonian Cooling is caused by the exchange of heat with the balloon's surroundings, and follows the equation:

DTNewtonian = ( Ti - Tair ) e^{-q t}

where:

DTNewtonian = Change in Balloon Temperature due to Newtonian Cooling

Ti = Balloon Initial Temperature

q = Newtonian Cooling Constant of Balloon

t = Elapsed Time

The conditions for strict Newtonian Cooling are that the temperature of the air be constant and that there be no adiabatic cooling.

Adiabatic cooling occurs when the volume of the lift gas changes due to a change in pressure, and follows the equation:

Ti ( Vgas,i ) ^{ggas - 1} = Tf ( Vgas,f ) ^{ggas - 1}

where:

Tf = Balloon Final Temperature

Vgas,i = Balloon Initial Volume

Vgas,f = Balloon Final Volume

ggas = Ratio of Specific Heats of Lift Gas = 1.45 for Hydrogen

The condition for strict adiabatic cooling are that the gas does not exchange any heat with its surroundings, i.e., no Newtonian Cooling.

Mathematically, the conditions for strict Newtonian and strict adiabatic cooling are mutually exclusive. However, both effects do occur during operation of the balloon, typically at the same time. The effects can be averaged in simulation.

Wind drag force on the platform can be approximated as:

F_D = ¹/₂ r(h) s_platform c_D v_rel v_rel

where:

F_D = Vector Force on Platform Due to Drag

s_platform = Cross Section Area of Platform

v_rel = Speed of Wind relative to Platform (Magnitude of v_rel)

v_rel = Velocity Vector of Wind Relative to Platform

Note that:

v_rel = v_wind - v_platform

where:

v_wind = Wind Velocity Vector

v_platform = Platform Velocity Vector

If the platform can be approximated as a cylinder, it's cross section can be estimated:

s_platform = p r_cyl² ( z_cyl ^. v_rel ) / ( z_cyl v_rel ) + r_cyl ( z_cyl ^x v_rel) / v_rel

where:

r_cyl = Radius of Platform in Cylindrical Approximation

z_cyl = Length of Platform (Magnitude of zcyl)

z_cyl = Platform Axis Vector

(Note that z_cyl ^. v_rel and z_cyl ^x v_rel are the vector scalar and cross products, respectively.) This approximation assumes that the lift cells are the primary source of drag for the platform, and thus that only lift cell orientation is relavent to this computation.

The gondola restoring torque on the platform is:

t_restore = ( F_tension + W_gondola ) ^x r_cyl

where:

t_restore = Restoring Torque Vector on Platform

F_tension = Vector Force on Platform due to Tension

W_gondola = Vector Weight of Platform Gondola

r_cyl = Radius Vector Relative to Platform Center of Pressure

Note that this approximation assumes that the tether mount and gondola are relatively close together, most of the platform mass resides in the gondola, and that the platform will tend to revolve about its center of pressure rather than its center of mass. This means that r_cyl is effectively the magnitude of r_cyl, even though r_cyl is an integral rather than a geometric radius.

Approximating Virtual Beanstalk Behavior -- The Tether

The shape of the tether can be approximated by the catenary equation:

y = a cosh( ^x / _a ) - a

where:

y = Vertical Displacement of a Point on the Tether Relative to the Tether's Mathematical Origin

x = Horizontal Displacement of a Point on the Tether Relative to the Tether's Mathematical Origin

a = Ratio of Tether Horizontal to Vertical Tension Components at the Tether's Mathematical Origin

Note that cosh( x ) is the hyperbolic cosine of x. Likewise, sinh( x ) is the hyperbolic sine of x.

Note that a = F_tension,x / F_tension,y at the the tether's mathematical origin. However, the virtual beanstalk tether will not extend all the way to its mathematical origin. Rather, it will behave as though it is "clamped" at the ground station, following this approximation only until that point. (The mathematical origin is where the tether weight and platform bouyancy exactly cancel out. This is not the case at the ground station, but the tether doesn't extend past it into the ground.) Note also that this approximation assumes strictly static behavior -- no motion -- and that the only forces on the tether other than tension are exerted at the platform and ground station.

For a strictly static tether extending from its origin, F_tension,y, which is the vertical component of the tether tension at a vertical position y, is due entirely to the weight of the length of tether theoretically hanging at or below y (reaching all the way down to the origin). Thus, F_tension,o (vertical tension at the origin) is roughly equal to the linear density of the tether because there is nothing below the origin. By contrast, the only horizontal tension source of consequence is that exerted on the tether mount and portions extending above y. So only the wind drag on the platform and tether portions overhead are of importance. Thus:

a = ( F_D + SF_D,tether ) / l_tether g

where:

SF_D,tether = Magnitude of Vector Sum of Wind Drag on Tether

l_tether = Linear Density of Tether

The length of the tether curve from its mathematical origin is given by the integral of the derivative of the catenary equation:

l_curve = a sinh ( ^x / _a )

where:

l_curve = Tether Curve Length from its Mathematical Origin.

Thus:

l_tether = a ( sinh( ^yg / _a ) - sinh ( ^{( yg + h )} / _a ) )

where

l_tether = Length of Tether Between Ground Station and Platform

y_g = y Position of Ground Station, relative to origin

The mass of the tether is:

m_tether = l_tether l_tether

The length of rope between the tether and the ground station is what we will actually deploy. l_curve is just a mathematical construct to provide a description consistent with the catenary equation, and would only be the length of the tether if the platform lift force fell so low that the tether began to drag the ground.

This is only an approximation because there is an additional force on the tether besides tension and the force applied at the end points: wind drag. Though small, the tether does experience horizontal wind drag. The wind changes force and direction over the length of the tether, which averages its force out to an even smaller amount over the length of the entire rope, allowing us to use the catenary approximation. However, over short distances wind can deflect the tether by an amount roughly equal to Dx = c_a Dy2 + c_b Dy + c_c, where Dx is the parabolic wind deflection and the constants c_i will vary from wind layer to wind layer. Because the wind forces on the tether are so much smaller than those on the platform balloons, the wind deflection at any given altitude will be quite small compared to the overall behavior of the tether.

For very large x, the approximation cosh( ^x / _a) = l_curve + a holds true, providing a rough guess for the tether length required for a known altitude and origin.

If you break the tether down into a bunch of segments of length Dl, then the wind drag on each is:

DF_D,tether = ¹/₂ r(h) c_D,tether i v_rel,tether v_rel,tether Dl sin(q)

where:

DF_D,tether = Wind Drag Vector Force on Tether Segment

c_D,tether = Coefficient of Drag of Tether Segment

v_rel,tether = Relative Velocity of Wind to Tether

q = Tether Segment Angle from Horizontal

i = Tether Diameter

Since Dh = Dl sin(q), if you know the height of a wind layer (of constant wind velocity) that a segment passes through, you can arbitrarily select Dl so that DF_D,tether is the force exerted on the tether by that wind layer. You can thus add together the forces exerted by all wind layers through which the tether passes, arriving at a value for SF_D,tether.

When in motion, each of these arbitrary tether segments also has angular momentum, and transmits angular momentum to neighboring segments.

All of this takes a bit of simulation.

Elementary Numerical Integration

If you have some formula for a particle's acceleration at any given time, a(t), and you want to simulate its behavior over time, you may use either continuous integration or numerical integration.

Continuous integration is known to every student of calculus.

a(t) = v'(t) = s''(t)

The acceleration is the time derivative of the velocity, which is in turn the time derivative of the position. This is true so long as the acceleration is continuous. If the acceleration is constant, the following integration steps hold true.

a(t) = a_o

v(t) = a_o t + v_o

s(t) = ¹/₂ a_o t² + v_o t + s_o

However, if the acceleration is not continuous, has a random component, or is simply to difficult to integrate, this method becomes impractical. It can be approximated, though, by considering only very small timesteps and treating acceleration as a finite step in velocity occuring during that time. Acceleration is then the discrete difference in velocity during a discrete time, rather than the derivative of a continuous velocity over all time. Velocity can be treated the exact same way.

Dv_n = a_n

Ds_n = v_n = v_n-1 + a_n Dt

s_n = s_n-1 + Ds_n Dt = s_n-1 + v_n-1 Dt + a_n Dt²

This is Euler's method of integration.

At each step in Euler's Method, the difference is taken at face value, without further attempt at continuous integration. The rates of change are assumed constant over the entier timestep. Note that this method of integration -- using differentials rather than derivatives -- grows more inexact as the number of steps increases. Assuming that the acceleration was constant, the numerical integration step to determine velocity from acceleration would still be quite accurate, but the value arrived at for displacement would be off by ¹/₂ a_n Dt². This can be a significant error. To avoid it, one can apply the integration twice and average the displacement differential over two steps.

Dv_n = a_n

v_n = v_n-1 + Dv_n Dt

s_n = s_n-1 + v_n-1 Dt + ¹/₂ Dv_n Dt²

This is the second order Runge-Kutte Method of Numerical integration. The first order Runge-Kutte Method is Euler's Method. The order of the integration step is determined by the number of steps over which the integration is averaged.

These techniques are applicable not only to linear motion of particles but to cooling processes and wind changes as well. They are very useful in simulating the behavior of the Virtual Beanstalk.