#include #include #include #include /********************************************************/ /* */ /* water.c - Water Rocket Simulation */ /* */ /* Michael D. Holder February 1, 1998 */ /* Original */ /* Revision B. February 4, 1998 */ /* Solve for volume increment */ /* Revision C. February 12, 1998 */ /* Print formatting added */ /* Revision D. February 24, 1998 */ /* Nozzle coefficient added */ /* Revision E. April 12, 1998 */ /* Fixed bug with volume increment solver */ /* Revision F. June 6, 1998 */ /* Added thrust coefficient */ /* Revision G. June 18, 1998 */ /* Read motor parameters from file */ /* Revision H. July 26, 1998 */ /* Read initial alt. and vel. from file */ /* Revision I. August 16, 1998 */ /* Read drag cooef. from file added */ /* Added macro for GetFloat to print */ /* variable name */ /* Revision j. September 21, 1998 */ /* Page formatting removed. */ /* */ /* */ /* */ #define Rev_Letter "J" /* */ /********************************************************/ #define GetFloat(a) xGetFloat(#a,&a) typedef struct { float pressure; /* p.s.i. abs */ float total_volume; /* cu. in. */ float air_volume; /* cu. in. */ float nozzle_area; /* sq. in. */ float dry_mass; /* slugs */ float ave_pressure; /* p.s.i. abs */ float ave_volume; /* cu. in. */ } motor_t; float delta_t; /* global time step */ float print_interval; /* seconds */ int print_iterations; const float amb_pressure=14.7; /* P.S.I. */ const float accel_gravity=32.17; /* acceleration due to gravity */ const float rho_water=0.0011229586; /* slugs/cu in. */ const float pi=3.141592654; float discharge_coef=1.0; /* Adjusts nozzle velocity */ float thrust_coef=1.0; /* Adjust thrust efficiency */ float drag_coef=1.5; /********************************************************/ float area(float diameter) { return pi*diameter*diameter/4; } /********************************************************/ void print_header() { int i; /* 1234567890 */ printf("\n Time"); printf(" Pressure"); printf(" Thrust"); printf(" Mass"); printf(" Drag"); printf(" Accel."); printf(" Velocity"); printf(" Altitude"); printf("\n"); for (i=0;i<8;i++) printf("----------"); printf("---"); } /********************************************************/ float drag(float velocity,float area) { const float rho_air=0.0023783028; /* slugs/cu. in. */ return area*drag_coef*rho_air*velocity*velocity/2.0; } /********************************************************/ float adiabatic(float p1,float v1,float v2) { const float gamma=1.4; return p1*pow(v1/v2,gamma); } /********************************************************/ float ex_velocity(float pressure,float area) { return discharge_coef*area*sqrt(2*(pressure-amb_pressure)/rho_water); } /********************************************************/ float print_state(motor_t motor) { printf("\nPressure= %7.1f Air Volume=%8.2f Total=%8.2f\n", \ motor.pressure,motor.air_volume,motor.total_volume); printf("\nNozzle Area= %8.5f Dry Mass=%8.4f\n", \ motor.nozzle_area,motor.dry_mass); } /********************************************************/ float mass(motor_t motor) { return motor.dry_mass \ +(motor.total_volume-motor.ave_volume)*rho_water; } /********************************************************/ float thrust(motor_t motor) { return discharge_coef*thrust_coef*2.0*motor.nozzle_area* \ (motor.ave_pressure-amb_pressure); } /********************************************************/ void step_pressure(motor_t *motor) { int i; float new_pressure; float new_volume; float delta_v; float ave_pressure; /* solve for dP and dV on adiabat */ ave_pressure=motor->pressure; if (motor->air_volume < motor->total_volume) { for (i=0;i<3;i++) { delta_v=delta_t*ex_velocity(ave_pressure,motor->nozzle_area); new_volume=motor->air_volume+delta_v; new_pressure=adiabatic(motor->pressure,motor->air_volume,new_volume); ave_pressure=(motor->pressure+new_pressure)/2; } motor->ave_volume=(motor->air_volume+new_volume)/2.0; motor->ave_pressure=ave_pressure; motor->pressure=new_pressure; motor->air_volume=new_volume; } else { motor->pressure=amb_pressure; motor->ave_pressure=amb_pressure; } } /********************************************************/ float solve_t(motor_t motor) { const int interations=10000; int i; float v; float sum; float delta_v; float sigma; float last_t; float next_t; delta_v=(motor.total_volume-motor.air_volume)/(interations-1); sum=0.0; last_t=delta_v/ex_velocity(motor.pressure,motor.nozzle_area); /* printf("dV=%10.4e last_t=%10.4e\n",delta_v,last_t); */ for (i=1;i -1.0) { step_pressure(&motor_1); acceleration=(thrust(motor_1)-drag(velocity,proj_area)) \ /mass(motor_1)-accel_gravity; velocity=velocity+acceleration*delta_t; altitude=altitude+velocity*delta_t; print_count-=1; if(print_count<=0) { print_count=print_iterations; printf("\n%10.2f",time); printf("%10.3f",motor_1.ave_pressure); printf("%10.3f",thrust(motor_1)); printf("%10.3f",mass(motor_1)); printf("%10.3f",drag(velocity,proj_area)); printf("%10.2f",acceleration); printf("%10.1f",velocity); printf("%10.1f",altitude); } time=time+delta_t; iteration_count+=1; total_impulse=total_impulse+thrust(motor_1)*delta_t; } printf("\n\ndone t=%8.4f Alt=%8.3f Total Impulse=%8.4f %8d Iterations\n\n",\ time,altitude,total_impulse,iteration_count); } /********************************************************/ /********************************************************/ /******/