|
|
'Free Balloon Simulator
'**************************************
'user-defined variable types
TYPE cartesian 'cartesian (x,y,z) vector
x AS SINGLE
y AS SINGLE
z AS SINGLE
mag AS SINGLE
END TYPE
TYPE spherical 'spherical (radius, altitude, azimuth) vector
R AS SINGLE
theta AS SINGLE
rho AS SINGLE
END TYPE
TYPE platform 'platform parameters
position AS cartesian
velocity AS cartesian
orientation AS spherical
force AS cartesian
radius AS SINGLE
mass AS SINGLE
cd AS SINGLE
volume AS SINGLE
pressure AS SINGLE
temperature AS SINGLE
moles AS SINGLE
max AS SINGLE
cooling AS SINGLE
END TYPE
TYPE segment 'segment parameters
position AS cartesian
velocity AS cartesian
orientation AS spherical
radius AS SINGLE
cd AS SINGLE
END TYPE
'**************************************
'Subroutines & functions
'AIRPRESSURE estimates ambient air pressure at altitude
DECLARE FUNCTION AirPressure (h)
'AIRTEMPERATURE estimates ambient air temperature at altitude
DECLARE FUNCTION AirTemperature (h)
'GASDENSITY estimates the density of a gas, using universal gas equation
DECLARE FUNCTION GasDensity (pressure, temperature, mole.weight)
'WINDDISTRIBUTION creates an upper level wind distribution
DECLARE SUB WindDistribution (wind() AS spherical)
'MAGNITUDE computes the magnitude of a cartesian vector
DECLARE FUNCTION Magnitude (v AS cartesian)
'SPHERETOCART converts a spherical vector to a cartesian vector
DECLARE SUB SphereToCart (v AS spherical, u AS cartesian)
'DOTPRODUCT computes scalar product of two cartesian vectors
DECLARE FUNCTION DotProduct (v1 AS cartesian, v2 AS cartesian)
'CROSSPRODUCT computes the cross product of two cartesian vectors
DECLARE SUB CrossProduct (v1 AS cartesian, v2 AS cartesian, u AS cartesian)
'DRAG computes drag force on a tether segment or the platform
DECLARE SUB drag (b AS segment, wind() AS spherical, D AS cartesian)
'LIFT computes lift force of platform
DECLARE FUNCTION Lift (b AS platform)
'LIFTCELL computes lift cell parameters
DECLARE SUB LiftCell (b AS platform)
'PLATFORMTOSEGMENT converts from platform to segment variable type
DECLARE SUB PlatformToSegment (v1 AS platform, v2 AS segment)
'DISPLAYPARAM provides a readout of platform parameters
DECLARE SUB DisplayParam (balloon AS platform)
'**************************************
'Global and Shared Variables
'**************************************
'Constants
CONST P.0 = 99375 'Mean Sea level air pressure (Pa)
CONST r.h = .00014 'Mean Atmospheric Density Decay Rate (1/m)
CONST T.0 = 288 'Mean Atmospheric Temperature at Surface(K)
CONST h.tropopause = 16500 'Mean Altitude of Tropopause(m)
CONST T.tropopause = 183 'Mean Temperature of Tropopause(K)
CONST day = 86400 'Mean length of day (s)
CONST R.gas = 8.314 'Universal Gas Constant (Nm/K-mole)
CONST m.stratosphere = .005164 'mean rate of temperature change in
'stratosphere (K / m)
CONST pi = 3.1415926# 'Geometric ratio of Circumference to Diameter
CONST W.air = 28.4 'Molecular Weight of dry air (g/mole)
CONST W.lift = 2 'Molecular weight of hydrogen (g/mole)
CONST gamma.lift = 1.45 'ratio of specific heats of hydrogen
CONST g = 9.81 'Acceleration due to gravity (m/s^2)
'**************************************
'shared variables
DIM SHARED Time AS SINGLE 'Simulation time
DIM SHARED timestep AS SINGLE 'calculation timestep
'**************************************
DIM wind(250) AS spherical 'Wind distribution array
DIM balloon AS platform 'platform parameters
DIM ptemp AS segment 'segment holding variable
DIM D AS cartesian 'Drag Force
DIM accel AS cartesian 'linear acceleration
GOSUB parameters 'initialize variables
DO WHILE Time < day
Time = Time + timestep
CALL LiftCell(balloon)
CALL DisplayParam(balloon)
CALL PlatformToSegment(balloon, ptemp)
CALL drag(ptemp, wind(), D)
balloon.force.x = D.x
balloon.force.y = Lift(balloon) + D.y - balloon.mass * g
balloon.force.z = D.z
accel.x = balloon.force.x / balloon.mass
accel.y = balloon.force.y / balloon.mass
accel.z = balloon.force.z / balloon.mass
balloon.position.x = balloon.position.x + balloon.velocity.x * timestep + .5 * accel.x * timestep * timestep
balloon.position.y = balloon.position.y + balloon.velocity.y * timestep + .5 * accel.y * timestep * timestep
balloon.position.z = balloon.position.z + balloon.velocity.z * timestep + .5 * accel.z * timestep * timestep
balloon.velocity.x = balloon.velocity.x + accel.x * timestep
balloon.velocity.y = balloon.velocity.y + accel.y * timestep
balloon.velocity.z = balloon.velocity.z + accel.z * timestep
bx = INT(balloon.position.x / 100) + 320
by = 400 - INT(balloon.position.y / 100)
IF bx > 639 THEN bx = bx - INT(bx / 640) * 640
IF bx < 1 THEN bx = bx - INT(bx / 640) * 640 + 640
PSET (bx, by)
LOOP
END
parameters:
SCREEN 12: CLS
'**************************************
'Initialize variables
'**************************************
Time = 1
timestep = .3
'platform parameters
balloon.orientation.R = 100 'platform length (m)
balloon.orientation.rho = 0 'platform azimuth
balloon.orientation.theta = 0 'platform altitude
balloon.mass = 10000 'platform mass (kg)
balloon.cd = .8 'platform coefficient of drag
balloon.radius = 12 'platform radius (m)
balloon.position.x = 0 'platform position
balloon.position.y = 0
balloon.position.z = 0
balloon.velocity.x = 0 'platform velocity
balloon.velocity.y = 0
balloon.velocity.z = 0
'#moles lift gas in balloon
balloon.moles = (balloon.mass + 1000) / (W.air - W.lift) * 1000
h.target = 25200 'Target Altitude
'max balloon volume
balloon.max = R.gas * balloon.moles * AirTemperature(h.target) / AirPressure(h.target)
balloon.cooling = .000625 'Newtonian Cooling Constant of Balloon
balloon.temperature = AirTemperature(balloon.position.y)
CALL LiftCell(balloon)
'**************************************
'display starting parameters
start1:
LOCATE 1, 1: PRINT "Balloon Mass:"; balloon.mass; "kg"
LOCATE 1, 40: PRINT "Balloon Length:"; balloon.orientation.R; "m"
LOCATE 2, 1: PRINT "Simulation Timestep:"; timestep; "s"
LOCATE 2, 40: PRINT "Newtonian Cooling Constant:"; balloon.cooling
LOCATE 3, 1: PRINT "Balloon Coefficient of Drag"; balloon.cd
LOCATE 3, 40: PRINT "Balloon Cylinder Radius"; balloon.radius
LOCATE 4, 1: PRINT "Balloon Starting Position:";
PRINT " x = "; balloon.position.x; ", y = "; balloon.position.y;
PRINT ", z = "; balloon.position.z
LOCATE 5, 1: PRINT "Target Altitude:"; h.target
LOCATE 7, 1: INPUT "Accept Defaults (y/n)"; a$
'**************************************
IF LEFT$(a$, 1) = "n" OR LEFT$(a$, 1) = "N" THEN
'manual entry of parameters
CLS
PRINT "Balloon Mass = "; balloon.mass
INPUT "Balloon Mass"; v$
IF VAL(v$) > 0 THEN balloon.mass = VAL(v$)
PRINT "Balloon Length = "; balloon.orientation.R
INPUT "Balloon Length"; v$
IF VAL(v$) > 0 THEN balloon.orientation.R = VAL(v$)
PRINT "Simulation Timestep = "; timestep
INPUT "Simulation Timestep"; v$
IF VAL(v$) > 0 THEN timestep = VAL(v$)
PRINT "Newtonian Cooling Constant = "; balloon.cooling
INPUT "Newtonian Cooling Constant"; v$
IF VAL(v$) > 0 THEN balloon.cooling = VAL(v$)
PRINT "Coefficient of Drag = "; balloon.cd
INPUT "Coefficient of Drag"; v$
IF VAL(v$) > 0 THEN balloon.cd = VAL(v$)
PRINT "Cylindrical Radius = "; balloon.radius
INPUT "Cylindrical Radius"; v$
IF VAL(v$) > 0 THEN balloon.radius = VAL(v$)
PRINT "Starting Position = "
PRINT " x = "; balloon.position.x; ", y = "; balloon.position.y;
PRINT ", z = "; balloon.position.z
INPUT " x "; x$
INPUT " y "; y$
INPUT " z "; z$
balloon.position.x = VAL(x$)
balloon.position.y = VAL(y$)
balloon.position.z = VAL(z$)
PRINT "Target Altitude = "; h.target
INPUT "Target Altitude"; v$
IF VAL(v$) > 0 THEN h.target = VAL(v$)
CLS
GOTO start1
END IF
CLS
'**************************************
'generate wind distribution
CALL WindDistribution(wind())
RETURN
FUNCTION AirPressure (h)
'Returns estimate for air pressure at altitude h
'Local variables:
'AirPressure=Ambient Air Pressure
'Global and shared variables:
'h=Altitude
'P.0=Air Pressure at sea level
'r.h=Atmospheric density decay constant
AirPressure = P.0 * EXP(-r.h * h)
END FUNCTION
FUNCTION AirTemperature (h)
'Returns estimate for atmospheric temperature
'Local Variables:
'AirTemperature=Ambient Air Temperature
'T.intersect=y-intersect (starting temperature) of linear temperature
' function
'm=rate of temperature change (K/m)
'height=altitude above boundary layer
'Global and Shared Variables:
'h=Altitude
'time=Time of day
'T.0=Mean Surface Temperature
'T.tropopause=Mean tropopause temperature
'm.stratosphere=mean rate of temperature change in stratosphere (K/m)
IF h < h.tropopause THEN
T.intersect = T.0
m = (T.tropopause - T.0) / h.tropopause
height = h
ELSE
T.intersect = T.tropopause
m = m.stratosphere
height = h - h.tropopause
END IF
variation = 5 * SIN(Time / day * 2 * pi)
AirTemperature = T.intersect + m * height + variation
END FUNCTION
SUB CrossProduct (v1 AS cartesian, v2 AS cartesian, u AS cartesian)
'Computes the cross product of two cartesian vectors
'Global and Shared Variables:
'u = Vector Cross-Product of two cartesian vectors
'v1 & v2 = vector factors
u.x = v1.y * v2.z - v1.z * v2.y
u.y = v1.z * v2.x - v1.x * v2.z
u.z = v1.x * v2.y - v1.y * v2.x
END SUB
SUB DisplayParam (balloon AS platform)
'Displays a readout of the platform parameters
LOCATE 1, 10: PRINT "Simulation Time = "; INT(Time); "s"; SPC(3);
LOCATE 1, 40: PRINT "Real Time = "; TIMER - start.time
LOCATE 2, 10: PRINT "Altitude = "; INT(balloon.position.y); "m"; SPC(5);
Horizontal = SQR(balloon.position.x * balloon.position.x + balloon.position.z * balloon.position.z)
LOCATE 2, 40: PRINT "Drift = "; INT(Horizontal); "m"; SPC(5);
LOCATE 3, 10: PRINT "Vert Velocity = "; INT(balloon.velocity.y); "m/s"; SPC(5);
Horizontal = SQR(balloon.velocity.x * balloon.velocity.x + balloon.velocity.z * balloon.velocity.z)
LOCATE 3, 40: PRINT "Horiz Velocity = "; INT(Horizontal); "m/s"; SPC(5);
END SUB
FUNCTION DotProduct (v1 AS cartesian, v2 AS cartesian)
'Returns the scalar product of two cartesian vectors
'Local Variables:
'DotProduct = Scalar product of two cartesian vectors
'Global and Shared variables
'v1 & v2 = cartesian vectors to multiply
DotProduct = v1.x * v2.x + v1.y * v2.y + v1.z * v2.z
END FUNCTION
SUB drag (b AS segment, wind() AS spherical, D AS cartesian)
'Computes drag force on a tether segment or on the platform. The segment
'is assumed cylindrical.
'Local Variables:
'vrel=relative velocity vector
'sigma=cross-section area
'u=spherical vector holding variable
'wv=cartesion vector holding variable
'dp=scalar holding variable
'hy=wind layer
'orient=cartesian orientation vector
'Density.air=ambient air density
'T.air=Ambient Air temperature
'P.air=ambient air pressure
'Global and Shared Variables:
'b = segment parameters
'wind() = wind distribution
'D = Drag force
DIM vrel AS cartesian, wv AS cartesian, u AS spherical, orient AS cartesian
hy = b.position.y / 100 'current wind layer
IF hy < 0 THEN hy = 0
u.R = wind(hy).R
u.rho = wind(hy).rho
u.theta = wind(hy).theta
CALL SphereToCart(u, wv)
CALL SphereToCart(b.orientation, orient)
vrel.x = wv.x - b.velocity.x
vrel.y = wv.y - b.velocity.y
vrel.z = wv.z - b.velocity.z
vrel.mag = Magnitude(vrel)
IF vrel.mag > 0 THEN
dp = DotProduct(orient, vrel)
CALL CrossProduct(orient, vrel, wv)
wv.mag = Magnitude(wv)
sigma = pi * b.radius * b.radius * dp / b.orientation.R / vrel.mag
sigma = sigma + b.radius / vrel.mag * wv.mag
P.air = AirPressure(b.position.y)
T.air = AirTemperature(b.position.y)
density.air = GasDensity(P.air, T.air, W.air)
D.mag = .5 * density.air * sigma * b.cd * vrel.mag
D.x = D.mag * vrel.x
D.y = D.mag * vrel.y
D.z = D.mag * vrel.z
D.mag = D.mag * vrel.mag
ELSE 'no relative velocity
D.mag = 0
D.x = 0
D.y = 0
D.z = 0
END IF
END SUB
FUNCTION GasDensity (pressure, temperature, mole.weight)
'Returns estimate for gas density
'Local Variables:
'Gas Density = Density of Gas (kg/m^3)
'Global and Shared Variables
'pressure = gas pressure (Pa)
'temperature = gas temperature (K)
'R.gas=Universal Gas Constant (N m/K-Mole)
'mole.weight=Molar weight of gas (g/mole)
GasDensity = pressure / temperature / R.gas * mole.weight / 1000
END FUNCTION
FUNCTION Lift (b AS platform)
'Returns lift force of platform
'Local variables:
'P.air=Ambient air pressure
'T.air=Ambient Air Temperature
'l=holding variable
'Global and Shared Variables:
'b = platform parameters
P.air = AirPressure(b.position.y)
T.air = AirTemperature(b.position.y)
l = (P.air * W.air / T.air - b.pressure * W.lift / b.temperature)
Lift = l * g * b.volume / R.gas / 1000
END FUNCTION
SUB LiftCell (b AS platform)
'Computes the current properties of the lift cells.
'Cell temperature is computed using Euler's method.
'Local Variables:
'delta.T=Temperature difference between balloon and outside air
'dT.Newtonian = Newtonian cooling rate
'Previous.volume = previously recorded cell volume
'Global and shared variables:
'b = platform parameters
'timestep= calculation timestep
'gamma.lift = ratio of specific heats of lift gas
STATIC previous.volume
y = b.position.y
b.pressure = AirPressure(y) 'lift cell pressure
'lift cell volume
b.volume = b.temperature * R.gas * b.moles / b.pressure
IF b.volume > b.max THEN 'lift cell venting
'Note: Venting is assumed instantaneous
b.moles = b.max * b.pressure / b.temperature / R.gas
b.volume = b.max
END IF
'Newtonian cooling of lift cells
T.air = AirTemperature(y)
delta.t = b.temperature - T.air
dT.Newtonian = -delta.t * b.cooling
'Adiabatic Cooling of Lift Cells
IF previous.volume = 0 THEN
'first time
previous.volume = b.volume
ELSE
dT.Adiabatic = b.temperature * (previous.volume / b.volume) ^ (gamma.lift - 1)
dT.Adiabatic = dT.Adiabatic - b.temperature
END IF
previous.volume = b.volume
dT = dT.Newtonian + dT.Adiabatic
b.temperature = b.temperature + dT * timestep
END SUB
FUNCTION Magnitude (v AS cartesian)
'Computes magnitude of cartesian vector and returns as v.mag
Magnitude = SQR(v.x * v.x + v.y * v.y + v.z * v.z)
END FUNCTION
SUB PlatformToSegment (v1 AS platform, v2 AS segment)
'Converts a platform variable to a segment variable by copying the
'common information
v2.position.x = v1.position.x
v2.position.y = v1.position.y
v2.position.z = v1.position.z
v2.velocity.x = v1.velocity.x
v2.velocity.y = v1.velocity.y
v2.velocity.z = v1.velocity.z
v2.orientation.R = v1.orientation.R
v2.orientation.rho = v1.orientation.rho
v2.orientation.theta = v1.orientation.theta
v2.radius = v1.radius
v2.cd = v1.cd
END SUB
SUB SphereToCart (v AS spherical, u AS cartesian)
'Converts a spherical vector to a cartesian vector
'Global and shared variables:
'v=spherical (radius, altitude and azimuth) vector
'u=cartesian (x,y,z) vector
u.mag = v.R
u.x = v.R * COS(v.theta) * COS(v.rho)
u.y = v.R * SIN(v.theta)
u.z = v.R * COS(v.theta) * SIN(v.rho)
END SUB
SUB WindDistribution (wind() AS spherical)
'Randomly generates a wind distribution to 25000m, assuming only horizontal
'winds. Wind array is arranged as a set of 250 sections, each 100m thick.
'Local Variables:
'h=altitude
'P.h=air pressure at altitude
'hy=position in array
'prob=probability of wind change
'wv=cartesian holding vector
'Global and Shared Variables:
'wind() = Array of wind magnitudes and directions
'P.0 = Sea level air pressure
DIM wv AS cartesian
RANDOMIZE TIMER 'initialize random number generator
wind(0).R = INT(RND(1) * 10) 'ground level winds
'circle grid for x-z wind speed chart
CIRCLE (40, 80), 10 'each circle measures 10m/s
CIRCLE (40, 80), 20
CIRCLE (40, 80), 30
'ground line
LINE (0, 400)-(639, 400), 2
FOR hy = 1 TO 250 'reset wind orientation
wind(hy).rho = 0
NEXT hy
'establish wind distribution
FOR h = 1 TO 25000 STEP 100
hy = INT(h / 100) + 1
P.h = AirPressure(h)
wind(hy).R = wind(0).R * SQR(P.0 / P.h) 'magnitude of wind in layer
IF h > h.tropopause AND h < 23000 THEN
wind(hy).R = wind(hy).R * .25 'less intense wind
wind(hy).rho = wind(hy).rho * .5 'less angular variation
END IF
prob = INT(RND(1) * 100)
IF prob > 94 THEN 'wind change
wind(0).R = INT(RND(1) * 10) 'change in magnitude
'change in direction
IF prob = 95 OR prob = 96 THEN
wind(hy).rho = wind(hy - 1).rho + pi / (INT(RND(1) * 10) + 1) + pi / 10
ELSEIF prob = 97 OR prob = 98 THEN
wind(hy).rho = wind(hy - 1).rho - pi / (INT(RND(1) * 10) + 1) - pi / 10
END IF
END IF
'display wind distribution
CALL SphereToCart(wind(hy), wv)
col% = INT(hy / 16) + 1
LINE (40, 400 - hy)-(40 + wv.x, 400 - hy), col% 'x projection of wind
LINE (40, 80)-(40 + wv.x, 80 + wv.z), col% 'x-z projection
NEXT h
END SUB