Attachment 'plat03.c'

Download

   1 // plat03.c
   2 //  Version 0.3    7-19-2011   KHL  more printout
   3 //  Version 0.2    7-19-2011   KHL  tweaks, magnet radius
   4 //  Version 0.1    7-19-2011   KHL  original version
   5 //
   6 // compile with:
   7 // cc -o plat03 plat03.c -lm
   8 //    using gcc compiler, -lm is the gcc math library
   9 // math:   sin, cos, atan, acos, sqrt
  10 // stdio:  printf, fprintf, fopen, fclose, exit
  11 //
  12 // Inputs via DEFINE statements, too lazy to parse argv
  13 //
  14 // Standard output,  tables suitable for inclusion in moin wiki
  15 //
  16 // plat03.dat output, X/Z coordinate points suitable for gnuplot
  17 //
  18 // Units are km and sec, tonnes and meganewtons (tonne-km/sec2)
  19 //  Z is vertical, X is horizontal
  20 //
  21 // Ignore coriolis effects, just centrifugal effects on apparent gravity
  22 // Ignore earth curvature
  23 //
  24 // No static/repair rotor assumed, all rotors combined in weight
  25 //    With static rotor just floating there in the levitation field,
  26 //    it will be under high tension at the top, break,  slide down
  27 //    to the bottom and pack up, won't it?  Iron is not very strong.
  28 
  29 #include <stdlib.h>
  30 #include <stdio.h>
  31 #include <math.h>
  32 
  33 #define  NAME   "plat03"             // data file name
  34 
  35 #define  PLAT    50.0                // platform bottom altitude, km
  36 #define  M0       3.4796             // rotor tons per kilometer (kg/m)
  37 #define  DEG0    38.0                // Incline starting angle, degrees
  38 #define  V0      3.50                // ground velocity, km/s
  39 #define  TRATIO   2.0                // Ratio of track deadweight to rotor weight
  40 #define  MYURI   (0.5*2.08)          // derated strength of kevlar
  41 
  42 // #define  MYURI   100.0               // make higher to remove tension effect
  43 
  44 // John's superconducting magnets assumed for smaller ground radius.
  45 //   Typical superconducting coils surround superconductor with copper
  46 //   matrix to provide current path for microfluctuations.  However,
  47 //   the copper resists high dB/dT needed for active control.
  48 //   Fast control currents may require much power and dissipate huge
  49 //   currents in low-resistance copper matrix and cause heat and quench.
  50 //   Launch loop uses copper wound magnets with NIB core, inefficient
  51 //   but quick response and no risk of quench.
  52 //   Needs more study.
  53 //
  54 // -----------------------------------------------------------------------
  55 // Magnet diameters from John's paper
  56 //
  57 #define  MAGGND  (0.62/(3.5*3.5))    // ground   magnet radius km per (km/s)^2
  58 #define  MAGTOP  (8.50/(3.5*3.5))    // platform magnet radius km per (km/s)^2
  59 
  60 // arbitrary cost functions per unit length relative to incline. 
  61 // The top magnets will be difficult to power and prevent overheating
  62 // The ground magnets will be supercooled and floating and anchored
  63 
  64 #define  CINCL   1.0                 // incline cost
  65 #define  CTOP    20.0                // top magnet cost
  66 #define  CGND    60.0                // ground magnet cost
  67 
  68 #define  DZ      1.0                 // vertical output step height, km
  69 #define  NSTEP   100                 // compute steps per output step
  70 
  71 #define  MU      398600.4418         // grav constant, km3/s2
  72 #define  SDAY    86164.098903691     // sidereal day, seconds
  73 #define  RE      6378.137            // Earth equatorial radius, km
  74 
  75 int main() {
  76    double pi2   = 8.0*atan(1.0)     ;// 
  77    double d2r   = pi2/360.0         ;// degrees to radians
  78    double om2   = pi2/SDAY          ;// radians/second 
  79           om2  *= om2               ;// (rad/s)^2
  80 
  81    double ztop  = PLAT-0.1*DZ/NSTEP ;// stop the main loop
  82    
  83    double z     = 0.0               ;// vertical height, km
  84    double x     = 0.0               ;// horizontal distance, km
  85    double l     = 0.0               ;// total length, km
  86    double mt                        ;// track/tube mass per length, tonnes/km
  87    double sumt  = 0.0               ;// total track/tube mass, tonnes
  88    double mr    = M0                ;// rotor mass per length, tonnes/km
  89    double sumr  = 0.0               ;// total rotor mass, tonnes
  90 
  91    double time  = 0.0               ;// time, seconds
  92    double dtime                     ;// time increment, seconds
  93 
  94    double T     = 0.0               ;// total track tension MN
  95    double T0                        ;// old total track tension MN
  96    
  97    double dz    = DZ                ;// computing step height km
  98           dz   /= (double) NSTEP    ;// smaller for inner loop
  99 
 100    double rad0  = d2r*DEG0          ;// Incline starting angle radians
 101 
 102    double dxdz  = 1.0/tan(rad0)     ;// starting "inverse slope" 
 103 
 104    double tmp                       ;// working variable computing new dxdz
 105 
 106    double dx                        ;// horizontal distance increment, km
 107    double dl                        ;// length increment, km
 108    double dldz                      ;// differential length per height
 109  
 110    double v     = V0                ;// rotor velocity km/s
 111    double v0                        ;// old velocity km/s
 112 
 113    double vx, vz                    ;// velocity components, km/s printout
 114    double vz0   = V0*sin(rad0)      ;// vertical rotor velocity km/s printout
 115    double vx0   = V0*cos(rad0)      ;// horizontal rotor velocity km/s printout
 116 
 117    double mv    = V0*M0             ;// mass flow rate, constant! tonne-km/sec
 118 
 119    double maggnd = MAGGND*V0*V0     ;// magnet radius near ground
 120    double magtop = MAGTOP*V0*V0     ;// magnet radius near platform
 121 
 122    double r                         ;// radius from center of earth, km
 123    double a                         ;// radial acceleration (down), km/s2
 124    double a0    = MU/(RE*RE)        ;// gravitational accel. at surface km/s2
 125           a0   -= om2*RE            ;// minus rotation centrifugal accel
 126    double eg2                       ;// 2* grav/rot potential diff GJ/tonne
 127    double eg20 =om2*RE*RE-2.0*MU/RE ;// potential at surface GJ/tonne
 128 
 129    int    i                         ;// inner loop index
 130    int    stop = 0                  ;// outer loop stop flag
 131 
 132    FILE   *datfile                  ;// data file
 133    char   name[128]                 ;// working string name
 134 
 135    sprintf( name, "%s.dat", NAME )  ;// data file name
 136    datfile = fopen( name, "wb" )    ;// gnuplot data file, x z pairs
 137 
 138    sprintf( name, "%s.c", NAME )    ;// source file name
 139 
 140    for( ;; ) {                       // break loop with if statements
 141       fprintf( datfile,
 142              "%9.3f%9.3f\n", x, z ) ;// print to gnuplot data file
 143 
 144       for( i=0 ; i<NSTEP ; i++ ) {   // inner loop
 145 
 146          // all this is from the previous iteration:
 147 
 148          v0    = v                  ;// save old velocity
 149          T0    = T                  ;// save old tension
 150 
 151          dldz  = sqrt(1.0+dxdz*dxdz);// l^2 = x^2 + y^2
 152          dx    = dxdz * dz          ;//
 153          dl    = dldz * dz          ;//
 154 	 dtime = dl/v               ;// transit time for dl at velocity v
 155          
 156          mt    = M0*TRATIO+T/MYURI  ;// track/tube weight 
 157          mr    = mv /  v            ;// rotor density increases with height
 158 
 159          x    += dx                 ;// horizontal km
 160          z    += dz                 ;// vertical km
 161          l    += dl                 ;// total length km
 162          time += dtime              ;// total time seconds
 163          sumt += dl * mt            ;// total track mass tonnes
 164          sumr += dl * mr            ;// total rotor mass tonnes
 165 
 166          // start of new iteration, dz higher.  What is new acceleration
 167          // and velocity and slope at the top of the new segment?
 168 
 169          r     = RE + z             ;// earth radius from altitude
 170          a     = MU/(r*r)           ;// gravitational acceleration at altit. z
 171          a    -= om2*r              ;// minus centrifugal acceleration
 172          eg2   = om2*r*r -2.0*MU/r  ;// twice the potential difference
 173          eg2  -= eg20               ;// minus the ground potential difference
 174          v     = sqrt(V0*V0-eg2)    ;// velocity, reduced by potential difference
 175 
 176          // we know from Archimedes' "chain hanging down a slope"
 177          // thought experiment that the track tension increases by the
 178          // extra vertical weight of track, not the horizontal weight
 179 
 180          T    += dz*a*mt            ;//
 181  
 182          // The total force pushing this segment horizontally (in the
 183          // X direction) is zero.
 184          //
 185          // A force of mv*v0-T0 is pushing on the segment from below.
 186          // The horizontal portion of that, pushing in the plus x
 187          // direction, is (dx0/dl0)*(mv*v0-T0).  mv is constant, v0
 188          // (rotor speed) is a function only of gravitational height,
 189          // and T0 is the accumulated tension, computed previously.
 190          // We have also already computed the previous derivative,
 191          // dx0/dl0 AKA dxdl0.  
 192          //
 193          // The force pushing on the segment from above is (dx/dl)*(mv*v-T)
 194          // As before, we know mv, v, and T.
 195          //
 196          // The horizontal forces are equal and opposite:
 197          // (dxdl)*(mv*v-T) = (dxdl0)*(mv*v0-T0)
 198          // 
 199          // From this, we can compute dxdl:
 200          // dxdl = (dxdl0)*(mv*v0-T0)/(mv*v-T)
 201          //
 202          // But what we really care about is dxdz, the increment of the
 203          // horizontal extent of the incline. 
 204          //
 205          // dl is the pythagorean sum of dx and dz:
 206          //
 207          // dl^2    = dz^2 + dx^2
 208          //
 209          // Compute dxdl from dxdz:
 210          //
 211          // dxdl^2  = dx^2      / dl^2 
 212          // dxdl^2  = dx^2      / ( dz^2      + dx^2      )
 213          // dxdl^2  = dx^2/dz^2 / ( dz^2/dz^2 + dx^2/dz^2 )
 214          // dxdl^2  = dxdz^2    / ( 1         + dxdz^2    ) 
 215          //
 216          // Compute dxdz from dxdl:
 217          //
 218          // dxdz^2  = dxdl^2    / ( 1         - dxdl^2    )
 219          //
 220          // Lastly, remember that
 221          //
 222          // dx/dl = (dx/dz) / (dl/dz)
 223          // 
 224          // Now we know enough to perform the new dxdz calculation.
 225          // This is easier to compute if we square both sides.
 226          //
 227          // tmp  = dxdl^2 = ( (dxdz/dldz)*(mv*v0-T0)/(mv*v-T) )^2
 228          // dxdz = sqrt( tmp / ( 1 - tmp ) )
 229          //
 230          // Cumulative build simplifies juggling terms for experiments. 
 231          // The compiler will optimize this sequence.
 232 
 233          tmp   = 1.0                ;//
 234          tmp  *= mv*v0 - T0         ;// old compression
 235          tmp  /= mv*v  - T          ;// ratio of old over new compression
 236          tmp  *= dxdz/dldz          ;// aka dx/dl, from previous iteration
 237          tmp  *= tmp                ;// cheap square -> new (dxdl)^2
 238          dxdz  = sqrt(tmp/(1.0-tmp));// new dxdz
 239 
 240          if( z > ztop ) { stop=1 ;  break; }  // break out of inner loop
 241 
 242       } // inner loop
 243 
 244       if(  stop==1 ) {  break; }     // break out of outer loop, cheezy!
 245 
 246    } // outer print loop
 247 
 248    fclose( datfile );  // close gnuplot data file
 249 
 250    // loops end up with T, x, z, l, dxdz, sumt, sumr, time, mt
 251 
 252    // compute platform parameters
 253 
 254    vz    = v / dldz               ;// platform interface vertical V
 255    vx    = dxdz * vz              ;// platform interface horizontal v
 256 
 257    double dzdx     = 1.0/dxdz     ;//
 258    double ang      = atan(dzdx)   ;// platform half angle radians
 259    double dang     = ang / d2r    ;// platform half angle degrees
 260 
 261    double p_rotor  = 2.0*ang      ;// rotor under platform angle, radians
 262           p_rotor *= mr*magtop    ;// rotor under platform rmass
 263    double p_long   = magtop       ;// platform magnet radius, km
 264           p_long  *= 2.0*sin(ang) ;// platform width, km 
 265    double p_bulge  = magtop       ;// platform magnet radius, km
 266           p_bulge *= 1.0-cos(ang) ;// platform bulge height, km
 267    double p_high   = p_bulge      ;//
 268           p_high  += PLAT         ;// platform center height, km
 269    double p_force  = 2.0*mv*vz    ;// raw platform force, MN
 270           p_force -= a*p_rotor    ;// subtract rotor weight
 271           p_force -= 2.0*T/dldz   ;// subtract tension force
 272    double p_mass   = p_force/a    ;// platform mass, tonnes
 273 
 274    // estimate rotor speed and compression at top of platform bulge
 275 
 276    double vt       =v-a*p_bulge/v ;// rotor velocity at top of bulge
 277    double compr    = V0/vt        ;// rotor compression factor to edge
 278 
 279    // compute ramp parameters
 280 
 281    double r_deep   = 1.0-cos(rad0);// depth, arbitrary units
 282    double r_deeph  =1.0-r_deep/2.0;// half depth, arbitrary units
 283           r_deep  *= maggnd       ;// depth in kilometers
 284    double a_deeph  = acos(r_deeph);// half depth angle
 285    double r_tgnd   = 2.0*a_deeph  ;// down and up radians
 286           r_tgnd  += rad0         ;// ramp up radians
 287           r_tgnd  *= 2.0          ;// 2 ramps, half mass
 288           r_tgnd  += 0.5*pi2      ;// 2 ambits, half mass
 289           r_tgnd  *= maggnd       ;// total turnaround length
 290 
 291           sumt    *= 2.0          ;// track mass, both inclines
 292           sumr    *= 2.0          ;// rotor mass, both inclines
 293 
 294    // cost estimate
 295 
 296    double cgnd    = M0*r_tgnd     ;// length*mass ground magnets
 297           cgnd   *= CGND          ;// ground magnet "cost"
 298    double cplat   = CTOP*p_rotor  ;// platform magnet "cost"
 299    double cincl   = M0*CINCL*2.0*l;// incline "cost"
 300    double ctot = cgnd+cincl+cplat ;// total "cost"
 301 
 302    // print out results ---------------------------------------------------
 303 
 304    printf("\n");
 305    printf("||<-6>'''Design parameters'''            %56s ||\n", " "    );
 306    printf("||Design Tensile Strength ||%8.2f||!MegaYuri  "    , MYURI  );
 307    printf("||Platform Height         ||%8.2f|| km        ||\n", PLAT   );
 308    printf("||Surface angle           ||%8.2f|| degrees   "    , DEG0   );
 309    printf("||Track to Rotor ratio    ||%8.3f||           ||\n", TRATIO );
 310    printf("||<-6>'''Ground'''                       %56s ||\n", " "    );
 311    printf("||All rotors, lineal dens.||%8.2f|| kg/m      "    , M0     );
 312    printf("||Rotor ground speed      ||%8.3f|| km/s      ||\n", V0     );
 313    printf("||Gnd horizontal velocity ||%8.2f|| km/s      "    , vx0    );
 314    printf("||Gnd vertical velocity   ||%8.2f|| km/s      ||\n", vz0    );
 315    printf("||ground ramp depth       ||%8.3f|| km        "    , r_deep );
 316    printf("||grav. and centrif. accel||%8.3f|| m/s^2^    ||\n", a0*1e3 );
 317    printf("||<-6>'''Incline'''                      %56s ||\n", " "    );
 318    printf("||Incline horiz. distance ||%8.2f|| km        "    , x      );
 319    printf("||Incline vert. distance  ||%8.2f|| km        ||\n", z      );
 320    printf("||Incline length          ||%8.2f|| km        "    , l      );
 321    printf("||Incline time            ||%8.2f|| sec       ||\n", time   );
 322    printf("||Incline total track mass||%8.2f|| tonnes    "    , sumt   );
 323    printf("||Incline total rotor mass||%8.2f|| tonnes    ||\n", sumr   );
 324    printf("||<-6>'''Platform interface'''           %56s ||\n", " "    );
 325    printf("||track tension           ||%8.2f||!MegaNewton"    , T      );
 326    printf("||track lineal dens.      ||%8.3f|| kg/m      ||\n", mt     );
 327    printf("||rotor lineal dens.      ||%8.3f|| kg/m      "    , mr     );
 328    printf("||angle                   ||%8.2f|| degrees   ||\n", dang   );
 329    printf("||rotor speed             ||%8.3f|| km/s      "    , v      );
 330    printf("||rotor compression       ||%8.4f||           ||\n", compr  );
 331    printf("||horizontal velocity     ||%8.3f|| km/s      "    , vx     );
 332    printf("||vertical velocity       ||%8.3f|| km/s      ||\n", vz     );
 333    printf("||grav. and centrif. accel||%8.3f|| m/s^2^    "    , a*1e3  );
 334    printf("||<-3>                                   %6s  ||\n", " "    );
 335    printf("||<-6>'''Platform'''                     %56s ||\n", " "    );
 336    printf("||platform length         ||%8.3f|| km        "    , p_long );
 337    printf("||platform center height  ||%8.3f|| km        ||\n", p_high );
 338    printf("||platform mass           ||%8.0f|| tonnes    "    , p_mass );
 339    printf("||<-3>                                   %6s  ||\n", " "    );
 340    printf("||<-6>'''Estimated costs, arbitrary units'''%54s||\n", " "  );
 341    printf("||cost of ground magnets  ||%8.0f|| cost unit "    , cgnd   );
 342    printf("||cost of inclines        ||%8.0f|| cost unit ||\n", cincl  );
 343    printf("||cost of platform magnets||%8.0f|| cost unit "    , cplat  );
 344    printf("||cost total              ||%8.0f|| cost unit ||\n", ctot   );
 345    printf("||<-6>computed with [[attachment:%s|%s]]  ||\n", name, name );
 346    printf("\n");
 347 
 348    // print out summary ---------------------------------------------------
 349    // DEG0, V0, M0, p_long, p_mass, r_deep, c_tot
 350 
 351    printf("||<-7>%7.2f km platform edge height %24s||\n",
 352                   PLAT,                         " " );
 353    printf("||<-7>%7.3f MYuri design load,%7.3f Track/Rotor mass ratio||\n",
 354                   MYURI,                  TRATIO    );
 355    printf("|| ramp  ||depth||Vrotor||Mrotor||platform||platform||  cost  ||\n");
 356    printf("||degrees|| km  || m/s  || kg/m ||leng.km || tonnes ||relative||\n");
 357    printf("||%6.2f ||%5.3f||%6.3f||%6.3f|| %5.2f  ||%7.0f ||%8.0f||\n",
 358              DEG0,   r_deep, V0,     M0,    p_long,    p_mass,   ctot );
 359    printf("\n");
 360 }

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2011-05-29 23:20:22, 295.5 KB) [[attachment:damato1992.pdf]]
  • [get | view] (2011-07-19 22:50:24, 16.6 KB) [[attachment:plat03.c]]
 All files | Selected Files: delete move to page

You are not allowed to attach a file to this page.