diff --git a/.gitignore b/.gitignore index c351f7a..4b5390a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# ignore outputs +outputs/ + # ignore prototyping files *.ipynb diff --git a/outputs/cycle_comparison_20251126_102847/data/summary.json b/outputs/cycle_comparison_20251126_102847/data/summary.json new file mode 100644 index 0000000..f46175f --- /dev/null +++ b/outputs/cycle_comparison_20251126_102847/data/summary.json @@ -0,0 +1,36 @@ +{ + "requirements": { + "thrust_kN": 500.0, + "chamber_pressure_MPa": 15.0, + "propellants": "LOX/CH4", + "mixture_ratio": 3.2 + }, + "ideal_performance": { + "isp_vac_s": 355.618287029336, + "cstar_m_s": 1867.159619502034, + "mdot_kg_s": 153.79256378134085 + }, + "cycles": [ + { + "cycle": "Pressure-Fed", + "net_isp": 331.5232505089583, + "tank_pressure_MPa": 25.070814644840986, + "pump_power_kW": 0, + "efficiency": 1.0 + }, + { + "cycle": "Gas Generator", + "net_isp": 301.43986952883694, + "tank_pressure_MPa": 0.3, + "pump_power_kW": 4793.858544008899, + "efficiency": 0.9092571005685514 + }, + { + "cycle": "Staged Combustion", + "net_isp": 324.8927854987791, + "tank_pressure_MPa": 0.4, + "pump_power_kW": 4361.185940503696, + "efficiency": 0.98 + } + ] +} \ No newline at end of file diff --git a/outputs/cycle_comparison_20251126_102847/metadata.json b/outputs/cycle_comparison_20251126_102847/metadata.json new file mode 100644 index 0000000..9de79cb --- /dev/null +++ b/outputs/cycle_comparison_20251126_102847/metadata.json @@ -0,0 +1,9 @@ +{ + "name": "cycle_comparison", + "created": "2025-11-26T10:28:47.176649", + "files": [ + "data/summary.json" + ], + "completed": "2025-11-26T10:28:47.177093", + "success": true +} \ No newline at end of file diff --git a/outputs/methalox_ssto_20251126_100146/data/vehicle_summary.json b/outputs/methalox_ssto_20251126_100146/data/vehicle_summary.json new file mode 100644 index 0000000..e4954eb --- /dev/null +++ b/outputs/methalox_ssto_20251126_100146/data/vehicle_summary.json @@ -0,0 +1,42 @@ +{ + "mission": { + "delta_v_km_s": 9.5, + "target_orbit": "LEO 400 km", + "payload_kg": 1000 + }, + "engine": { + "name": "Methalox-500", + "propellants": "LOX/CH4", + "thrust_kN": 500.0, + "chamber_pressure_MPa": 10.0, + "isp_sl_s": 320.5256042628083, + "isp_vac_s": 346.6915059643235, + "mdot_kg_s": 159.06938469443352 + }, + "propellant": { + "lox_kg": 52550.773557232475, + "ch4_kg": 16422.116736635147, + "total_kg": 68972.89029386762, + "burn_time_s": 433.6025466268196 + }, + "tanks": { + "lox": { + "volume_m3": 47.438472185757625, + "diameter_m": 2.9020123520766594, + "length_m": 7.8560376383067245, + "mass_kg": 398.25268941572546 + }, + "ch4": { + "volume_m3": 40.02550932024184, + "diameter_m": 2.7422138148873865, + "length_m": 7.423447018281644, + "mass_kg": 280.0165139805067 + } + }, + "vehicle": { + "dry_mass_kg": 2578.269203396232, + "wet_mass_kg": 71551.15949726386, + "twr": 0.7125783985491686, + "structure_fraction": 0.022882457102664816 + } +} \ No newline at end of file diff --git a/outputs/methalox_ssto_20251126_100146/metadata.json b/outputs/methalox_ssto_20251126_100146/metadata.json new file mode 100644 index 0000000..12240e1 --- /dev/null +++ b/outputs/methalox_ssto_20251126_100146/metadata.json @@ -0,0 +1,11 @@ +{ + "name": "methalox_ssto", + "created": "2025-11-26T10:01:46.060424", + "files": [ + "data/vehicle_summary.json", + "plots/engine_dashboard.png", + "plots/mass_breakdown.png" + ], + "completed": "2025-11-26T10:01:46.470217", + "success": true +} \ No newline at end of file diff --git a/outputs/methalox_ssto_20251126_100146/plots/engine_dashboard.png b/outputs/methalox_ssto_20251126_100146/plots/engine_dashboard.png new file mode 100644 index 0000000..d32b713 Binary files /dev/null and b/outputs/methalox_ssto_20251126_100146/plots/engine_dashboard.png differ diff --git a/outputs/methalox_ssto_20251126_100146/plots/mass_breakdown.png b/outputs/methalox_ssto_20251126_100146/plots/mass_breakdown.png new file mode 100644 index 0000000..e3867fc Binary files /dev/null and b/outputs/methalox_ssto_20251126_100146/plots/mass_breakdown.png differ diff --git a/outputs/trade_study_results_20251126_101124/metadata.json b/outputs/trade_study_results_20251126_101124/metadata.json new file mode 100644 index 0000000..d617b5c --- /dev/null +++ b/outputs/trade_study_results_20251126_101124/metadata.json @@ -0,0 +1,9 @@ +{ + "name": "trade_study_results", + "created": "2025-11-26T10:11:24.230830", + "files": [ + "data/mixture_ratio_sweep.csv" + ], + "completed": "2025-11-26T10:11:24.235132", + "success": false +} \ No newline at end of file diff --git a/outputs/trade_study_results_20251126_101453/metadata.json b/outputs/trade_study_results_20251126_101453/metadata.json new file mode 100644 index 0000000..0f92173 --- /dev/null +++ b/outputs/trade_study_results_20251126_101453/metadata.json @@ -0,0 +1,9 @@ +{ + "name": "trade_study_results", + "created": "2025-11-26T10:14:53.688293", + "files": [ + "data/chamber_pressure_sweep.csv" + ], + "completed": "2025-11-26T10:14:53.692364", + "success": false +} \ No newline at end of file diff --git a/outputs/trade_study_results_20251126_101803/metadata.json b/outputs/trade_study_results_20251126_101803/metadata.json new file mode 100644 index 0000000..baa4a78 --- /dev/null +++ b/outputs/trade_study_results_20251126_101803/metadata.json @@ -0,0 +1,9 @@ +{ + "name": "trade_study_results", + "created": "2025-11-26T10:18:03.834991", + "files": [ + "data/chamber_pressure_sweep.csv" + ], + "completed": "2025-11-26T10:18:03.839727", + "success": false +} \ No newline at end of file diff --git a/outputs/trade_study_results_20251126_101903/data/chamber_pressure_sweep.csv b/outputs/trade_study_results_20251126_101903/data/chamber_pressure_sweep.csv new file mode 100644 index 0000000..02f8e6d --- /dev/null +++ b/outputs/trade_study_results_20251126_101903/data/chamber_pressure_sweep.csv @@ -0,0 +1,10 @@ +chamber_pressure,expansion_ratio,chamber_area,thrust_coeff_vac,chamber_volume,isp,contraction_ratio,thrust_coeff,exit_mach,chamber_diameter,isp_vac,mdot,throat_area,exit_area,chamber_length,exit_diameter,exhaust_velocity,mdot_ox,mdot_fuel,throat_diameter,cstar,nozzle_length,feasible +5.0,7.876186624172613,0.10082704335041623,1.746486770751889,0.025206760837604057,301.019889224399,4.0,1.586875848813031,2.957953479816565,0.3582973329128707,331.297029100828,67.7507533210052,0.025206760837604057,0.1985331525478551,0.20521283338589116,0.5027725736004974,2951.996696662452,51.619621577908724,16.13113174309648,0.17914866645643535,1860.2568681541907,0.4831111625577843,true +7.5,10.83051902476274,0.06469217023125673,1.7951545725751283,0.016173042557814182,312.7729915554404,4.0,1.6488342605505837,3.1491065807874272,0.2869993543079408,340.5289903311817,65.2048764125581,0.016173042557814182,0.17516244511070395,0.21412508071150738,0.4722539061431368,3067.2552576371595,49.6799058381395,15.524970574418592,0.1434996771539704,1860.2568681541907,0.49076979251734604,true +10.0,13.613407728206797,0.04734558646254097,1.827641287161869,0.011836396615635243,320.5256042628083,4.0,1.6897034333558136,3.2833310847743125,0.24552448544449687,346.6915059643235,63.62775387777341,0.011836396615635243,0.1611336931614096,0.2193094393194379,0.45294788896994387,3143.282417043869,48.47828866877974,15.149465208993668,0.12276224272224844,1860.2568681541907,0.49290663605889995,true +12.5,16.27729832096691,0.03721348009651728,1.8517506844085763,0.00930337002412932,326.2360342698122,4.0,1.7198069042188184,3.386854581730957,0.2176733204967956,351.26489971399786,62.514014753782604,0.00930337002412932,0.15143372927309406,0.22279083493790056,0.4391029634575724,3199.2826054720535,47.62972552669151,14.884289227091095,0.1088366602483978,1860.2568681541907,0.4930270421153389,true +15.0,18.85114636714012,0.030590870176051197,1.8707790567691638,0.007647717544012799,330.71899388603566,4.0,1.7434395630591322,3.471147402545745,0.1973562910474553,354.8744565324403,61.66662525160669,0.007647717544012799,0.14416824279673066,0.22533046361906808,0.42843985321603517,3243.2454213924916,46.98409542979557,14.682529821811116,0.09867814552372765,1860.2568681541907,0.49227377351745816,true +17.5,21.35315508630143,0.025933043466047308,1.8864161311991197,0.006483260866511827,334.38800589793334,4.0,1.7627813632494345,3.5422571327284817,0.1817112447161905,357.840706485907,60.98999934161398,0.006483260866511827,0.13843807474757605,0.22728609441047617,0.41983905400238336,3279.2261380389677,46.46857092694398,14.521428414669995,0.09085562235809524,1860.2568681541907,0.49111195006112696,true +20.0,23.79574885106925,0.02248348554045343,1.8996387950101252,0.005620871385113357,337.48042941959557,4.0,1.7790835823933957,3.6037640338905774,0.16919474842200788,360.3489586586384,60.4311316500132,0.005620871385113357,0.13375284380431918,0.228850656447249,0.412673490731407,3309.5524531676765,46.042766971438624,14.38836467857457,0.08459737421100394,1860.2568681541907,0.489757494921596,true +22.5,26.187906089521437,0.019828792011868247,1.9110605372156693,0.004957198002967062,340.14449966380903,4.0,1.7931276667925244,3.6579632541218174,0.1588924230850526,362.51558787291134,59.95782462957904,0.004957198002967062,0.12981863576886463,0.23013844711436843,0.4065590002748545,3335.6780576280926,45.68215209872689,14.275672530852152,0.0794462115425263,1860.2568681541907,0.4883194231433763,true +25.0,28.53639395834763,0.017724282734337772,1.9210907254407448,0.004431070683584443,342.4786872911498,4.0,1.8054327207275618,3.7064143957641447,0.15022402497412166,364.41824846899425,59.54917785065216,0.004431070683584443,0.1264467786840504,0.2312219968782348,0.4012443631067086,3358.5686187237543,45.370802171925455,14.178375678726704,0.07511201248706083,1860.2568681541907,0.4868558087874004,true diff --git a/outputs/trade_study_results_20251126_101903/data/grid_study.csv b/outputs/trade_study_results_20251126_101903/data/grid_study.csv new file mode 100644 index 0000000..b3eac07 --- /dev/null +++ b/outputs/trade_study_results_20251126_101903/data/grid_study.csv @@ -0,0 +1,25 @@ +chamber_pressure,contraction_ratio,expansion_ratio,chamber_area,thrust_coeff_vac,chamber_volume,isp,thrust_coeff,exit_mach,chamber_diameter,isp_vac,mdot,throat_area,exit_area,chamber_length,exit_diameter,exhaust_velocity,mdot_ox,mdot_fuel,throat_diameter,cstar,nozzle_length,feasible +5.0,3.0,7.876186624172613,0.07562028251281216,1.746486770751889,0.025206760837604057,301.019889224399,1.586875848813031,2.957953479816565,0.31029459241075624,331.297029100828,67.7507533210052,0.025206760837604057,0.1985331525478551,0.30054685184475316,0.5027725736004974,2951.996696662452,51.619621577908724,16.13113174309648,0.17914866645643535,1860.2568681541907,0.4831111625577843,true +5.0,4.0,7.876186624172613,0.10082704335041623,1.746486770751889,0.025206760837604057,301.019889224399,1.586875848813031,2.957953479816565,0.3582973329128707,331.297029100828,67.7507533210052,0.025206760837604057,0.1985331525478551,0.20521283338589116,0.5027725736004974,2951.996696662452,51.619621577908724,16.13113174309648,0.17914866645643535,1860.2568681541907,0.4831111625577843,true +5.0,5.0,7.876186624172613,0.1260338041880203,1.746486770751889,0.025206760837604057,301.019889224399,1.586875848813031,2.957953479816565,0.4005885962750258,331.297029100828,67.7507533210052,0.025206760837604057,0.1985331525478551,0.14464001754535236,0.5027725736004974,2951.996696662452,51.619621577908724,16.13113174309648,0.17914866645643535,1860.2568681541907,0.4831111625577843,true +5.0,6.0,7.876186624172613,0.15124056502562433,1.746486770751889,0.025206760837604057,301.019889224399,1.586875848813031,2.957953479816565,0.43882282091832314,331.297029100828,67.7507533210052,0.025206760837604057,0.1985331525478551,0.10174812805119472,0.5027725736004974,2951.996696662452,51.619621577908724,16.13113174309648,0.17914866645643535,1860.2568681541907,0.4831111625577843,true +8.0,3.0,11.398612239300189,0.045229374048974154,1.8025849624925732,0.015076458016324719,314.55234490527897,1.6582144143491868,3.179305106714798,0.23997463954087364,341.93848632389995,64.83602678498518,0.015076458016324719,0.17185069887017437,0.3079770291325167,0.4677682178086497,3084.704753165354,49.398877550464896,15.43714923452028,0.13854942273760681,1860.2568681541907,0.4914633045074857,true +8.0,4.0,11.398612239300189,0.060305832065298874,1.8025849624925732,0.015076458016324719,314.55234490527897,1.6582144143491868,3.179305106714798,0.27709884547521363,341.93848632389995,64.83602678498518,0.015076458016324719,0.17185069887017437,0.2153626443155983,0.4677682178086497,3084.704753165354,49.398877550464896,15.43714923452028,0.13854942273760681,1860.2568681541907,0.4914633045074857,true +8.0,5.0,11.398612239300189,0.07538229008162359,1.8025849624925732,0.015076458016324719,314.55234490527897,1.6582144143491868,3.179305106714798,0.3098059274846438,341.93848632389995,64.83602678498518,0.015076458016324719,0.17185069887017437,0.15718587381324076,0.4677682178086497,3084.704753165354,49.398877550464896,15.43714923452028,0.13854942273760681,1860.2568681541907,0.4914633045074857,true +8.0,6.0,11.398612239300189,0.09045874809794831,1.8025849624925732,0.015076458016324719,314.55234490527897,1.6582144143491868,3.179305106714798,0.3393753898642983,341.93848632389995,64.83602678498518,0.015076458016324719,0.17185069887017437,0.1164601748849938,0.4677682178086497,3084.704753165354,49.398877550464896,15.43714923452028,0.13854942273760681,1860.2568681541907,0.4914633045074857,true +11.0,3.0,14.691291302428622,0.03203420751978145,1.8380520932621829,0.010678069173260484,322.9957602363314,1.7027252667877666,3.327603013494704,0.20195846057652117,348.6663672626394,63.14115158860391,0.010678069173260484,0.15687462477185293,0.3119939110839473,0.4469216663186218,3167.5063721216193,48.107544067507746,15.033607521096169,0.11660077157897693,1860.2568681541907,0.49310853726192444,true +11.0,4.0,14.691291302428622,0.04271227669304194,1.8380520932621829,0.010678069173260484,322.9957602363314,1.7027252667877666,3.327603013494704,0.23320154315795386,348.6663672626394,63.14115158860391,0.010678069173260484,0.15687462477185293,0.22084980710525576,0.4469216663186218,3167.5063721216193,48.107544067507746,15.033607521096169,0.11660077157897693,1860.2568681541907,0.49310853726192444,true +11.0,5.0,14.691291302428622,0.053390345866302424,1.8380520932621829,0.010678069173260484,322.9957602363314,1.7027252667877666,3.327603013494704,0.2607272514795179,348.6663672626394,63.14115158860391,0.010678069173260484,0.15687462477185293,0.16396838002486472,0.4469216663186218,3167.5063721216193,48.107544067507746,15.033607521096169,0.11660077157897693,1860.2568681541907,0.49310853726192444,true +11.0,6.0,14.691291302428622,0.0640684150395629,1.8380520932621829,0.010678069173260484,322.9957602363314,1.7027252667877666,3.327603013494704,0.2856123939833083,348.6663672626394,63.14115158860391,0.010678069173260484,0.15687462477185293,0.12441376106558383,0.4469216663186218,3167.5063721216193,48.107544067507746,15.033607521096169,0.11660077157897693,1860.2568681541907,0.49310853726192444,true +14.0,3.0,17.831082237629115,0.02470728845597142,1.8636476130278354,0.008235762818657141,329.04126803497513,1.7345951553329948,3.4392764979945003,0.1773648688588622,353.52166865894003,61.98105295835041,0.008235762818657141,0.14685256410908365,0.31459252981372954,0.43241009686343007,3226.7925511751887,47.223659396838414,14.757393561512002,0.10240165478044679,1860.2568681541907,0.4926421027282828,true +14.0,4.0,17.831082237629115,0.032943051274628564,1.8636476130278354,0.008235762818657141,329.04126803497513,1.7345951553329948,3.4392764979945003,0.20480330956089357,353.52166865894003,61.98105295835041,0.008235762818657141,0.14685256410908365,0.2243995863048883,0.43241009686343007,3226.7925511751887,47.223659396838414,14.757393561512002,0.10240165478044679,1860.2568681541907,0.4926421027282828,true +14.0,5.0,17.831082237629115,0.04117881409328571,1.8636476130278354,0.008235762818657141,329.04126803497513,1.7345951553329948,3.4392764979945003,0.22897706109754531,353.52166865894003,61.98105295835041,0.008235762818657141,0.14685256410908365,0.16835614842072535,0.43241009686343007,3226.7925511751887,47.223659396838414,14.757393561512002,0.10240165478044679,1860.2568681541907,0.4926421027282828,true +14.0,6.0,17.831082237629115,0.04941457691194284,1.8636476130278354,0.008235762818657141,329.04126803497513,1.7345951553329948,3.4392764979945003,0.2508318030287284,353.52166865894003,61.98105295835041,0.008235762818657141,0.14685256410908365,0.12955912960459628,0.43241009686343007,3226.7925511751887,47.223659396838414,14.757393561512002,0.10240165478044679,1860.2568681541907,0.4926421027282828,true +17.0,3.0,20.8578441147071,0.02006273677941217,1.8835064797346852,0.0066875789264707235,333.70629158693737,1.7591875941509971,3.528894622816198,0.15982699973146014,357.28876478097726,61.114593202823755,0.0066875789264707235,0.1394884787531266,0.31644562373015045,0.421428816270828,3272.5408043910393,46.56349958310381,14.551093619719941,0.09227616131872876,1860.2568681541907,0.491364569435542,true +17.0,4.0,20.8578441147071,0.026750315705882894,1.8835064797346852,0.0066875789264707235,333.70629158693737,1.7591875941509971,3.528894622816198,0.18455232263745752,357.28876478097726,61.114593202823755,0.0066875789264707235,0.1394884787531266,0.22693095967031782,0.421428816270828,3272.5408043910393,46.56349958310381,14.551093619719941,0.09227616131872876,1860.2568681541907,0.491364569435542,true +17.0,5.0,20.8578441147071,0.03343789463235362,1.8835064797346852,0.0066875789264707235,333.70629158693737,1.7591875941509971,3.528894622816198,0.20633576941141413,357.28876478097726,61.114593202823755,0.0066875789264707235,0.1394884787531266,0.17148509797682868,0.421428816270828,3272.5408043910393,46.56349958310381,14.551093619719941,0.09227616131872876,1860.2568681541907,0.491364569435542,true +17.0,6.0,20.8578441147071,0.04012547355882434,1.8835064797346852,0.0066875789264707235,333.70629158693737,1.7591875941509971,3.528894622816198,0.22602951065363194,357.28876478097726,61.114593202823755,0.0066875789264707235,0.1394884787531266,0.13322832933294085,0.421428816270828,3272.5408043910393,46.56349958310381,14.551093619719941,0.09227616131872876,1860.2568681541907,0.491364569435542,true +20.0,3.0,23.79574885106925,0.016862614155340072,1.8996387950101252,0.005620871385113357,337.48042941959557,1.7790835823933957,3.6037640338905774,0.1465269503203759,360.3489586586384,60.4311316500132,0.005620871385113357,0.13375284380431918,0.31785093930599034,0.412673490731407,3309.5524531676765,46.042766971438624,14.38836467857457,0.08459737421100394,1860.2568681541907,0.489757494921596,true +20.0,4.0,23.79574885106925,0.02248348554045343,1.8996387950101252,0.005620871385113357,337.48042941959557,1.7790835823933957,3.6037640338905774,0.16919474842200788,360.3489586586384,60.4311316500132,0.005620871385113357,0.13375284380431918,0.228850656447249,0.412673490731407,3309.5524531676765,46.042766971438624,14.38836467857457,0.08459737421100394,1860.2568681541907,0.489757494921596,true +20.0,5.0,23.79574885106925,0.028104356925566787,1.8996387950101252,0.005620871385113357,337.48042941959557,1.7790835823933957,3.6037640338905774,0.18916547945379245,360.3489586586384,60.4311316500132,0.005620871385113357,0.13375284380431918,0.1738579736893029,0.412673490731407,3309.5524531676765,46.042766971438624,14.38836467857457,0.08459737421100394,1860.2568681541907,0.489757494921596,true +20.0,6.0,23.79574885106925,0.033725228310680144,1.8996387950101252,0.005620871385113357,337.48042941959557,1.7790835823933957,3.6037640338905774,0.20722040039624431,360.3489586586384,60.4311316500132,0.005620871385113357,0.13375284380431918,0.13601091012035657,0.412673490731407,3309.5524531676765,46.042766971438624,14.38836467857457,0.08459737421100394,1860.2568681541907,0.489757494921596,true diff --git a/outputs/trade_study_results_20251126_101903/data/summary.json b/outputs/trade_study_results_20251126_101903/data/summary.json new file mode 100644 index 0000000..c2fd33c --- /dev/null +++ b/outputs/trade_study_results_20251126_101903/data/summary.json @@ -0,0 +1,31 @@ +{ + "studies": { + "mixture_ratio_sweep": { + "parameter": "mixture_ratio", + "range": [ + 2.4, + 4.0 + ], + "optimal_mr": 3.0, + "optimal_isp_vac": 347.1521682740077, + "note": "Each point recalculates CEA thermochemistry" + }, + "chamber_pressure_sweep": { + "parameter": "chamber_pressure", + "range_MPa": [ + 5, + 25 + ], + "n_points": 9 + }, + "grid_study": { + "parameters": [ + "chamber_pressure", + "contraction_ratio" + ], + "total_designs": 24, + "feasible_designs": 24, + "constraint": "throat_diameter > 8 cm" + } + } +} \ No newline at end of file diff --git a/outputs/trade_study_results_20251126_101903/metadata.json b/outputs/trade_study_results_20251126_101903/metadata.json new file mode 100644 index 0000000..10a480b --- /dev/null +++ b/outputs/trade_study_results_20251126_101903/metadata.json @@ -0,0 +1,11 @@ +{ + "name": "trade_study_results", + "created": "2025-11-26T10:19:03.375007", + "files": [ + "data/chamber_pressure_sweep.csv", + "data/grid_study.csv", + "data/summary.json" + ], + "completed": "2025-11-26T10:19:03.381927", + "success": true +} \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index bddc08f..7de7bce 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,8 @@ dependencies = [ "numba>=0.60", "matplotlib>=3.9", "rocketcea>=1.2.1", + "polars>=1.35.0", + "tqdm>=4.66", ] [project.optional-dependencies] @@ -44,6 +46,9 @@ Repository = "https://github.com/cmflannery/rocket" requires = ["hatchling"] build-backend = "hatchling.build" +[tool.hatch.build.targets.wheel] +packages = ["rocket"] + [tool.pytest.ini_options] testpaths = ["tests"] addopts = "-v --tb=short" diff --git a/rocket/__init__.py b/rocket/__init__.py index b5b09d5..ce74795 100644 --- a/rocket/__init__.py +++ b/rocket/__init__.py @@ -21,6 +21,21 @@ __version__ = "0.3.0" # Core engine design +# Analysis framework +from rocket.analysis import ( + Distribution, + LogNormal, + MultiObjectiveOptimizer, + Normal, + ParametricStudy, + ParetoResults, + Range, + StudyResults, + Triangular, + UncertaintyAnalysis, + UncertaintyResults, + Uniform, +) from rocket.engine import ( EngineGeometry, EngineInputs, @@ -43,8 +58,19 @@ rao_bell_contour, ) +# Output management +from rocket.output import ( + OutputContext, + clean_outputs, + get_default_output_dir, + list_outputs, +) + # Visualization from rocket.plotting import ( + plot_cycle_comparison_bars, + plot_cycle_radar, + plot_cycle_tradeoff, plot_engine_cross_section, plot_engine_dashboard, plot_mass_breakdown, @@ -61,6 +87,12 @@ list_database_propellants, ) +# System-level design +from rocket.system import ( + EngineSystemResult, + design_engine_system, + format_system_summary, +) # Tank sizing from rocket.tanks import ( PropellantRequirements, @@ -100,12 +132,37 @@ "plot_performance_vs_altitude", "plot_engine_dashboard", "plot_mass_breakdown", + "plot_cycle_comparison_bars", + "plot_cycle_radar", + "plot_cycle_tradeoff", # Propellants "CombustionProperties", "get_combustion_properties", "get_optimal_mixture_ratio", "is_cea_available", "list_database_propellants", + # Output management + "OutputContext", + "get_default_output_dir", + "list_outputs", + "clean_outputs", + # Analysis framework + "ParametricStudy", + "UncertaintyAnalysis", + "MultiObjectiveOptimizer", + "StudyResults", + "UncertaintyResults", + "ParetoResults", + "Range", + "Distribution", + "Normal", + "Uniform", + "Triangular", + "LogNormal", + # System-level design + "EngineSystemResult", + "design_engine_system", + "format_system_summary", # Tanks "PropellantRequirements", "TankGeometry", diff --git a/rocket/analysis.py b/rocket/analysis.py new file mode 100644 index 0000000..7fe10e0 --- /dev/null +++ b/rocket/analysis.py @@ -0,0 +1,1210 @@ +"""Parametric analysis and uncertainty quantification for Rocket. + +This module provides general-purpose tools for trade studies, sensitivity +analysis, and uncertainty quantification. The design is introspection-based +to avoid brittleness when dataclass fields change. + +Key Design Principles: +- Works with ANY frozen dataclass + computation function +- Uses dataclass introspection to validate parameters (no hardcoding) +- Automatically discovers output metrics from return types +- Unit-aware parameter ranges + +Example: + >>> from rocket import EngineInputs, design_engine + >>> from rocket.analysis import ParametricStudy, Range + >>> + >>> study = ParametricStudy( + ... compute=design_engine, + ... base=inputs, + ... vary={"chamber_pressure": Range(5, 15, n=11, unit="MPa")}, + ... ) + >>> results = study.run() + >>> results.plot("chamber_pressure", "isp_vac") +""" + +import dataclasses +import itertools +from collections.abc import Callable, Sequence +from dataclasses import dataclass, fields, is_dataclass, replace +from pathlib import Path +from typing import Any, Generic, TypeVar + +import numpy as np +import polars as pl +from beartype import beartype +from numpy.typing import NDArray +from tqdm import tqdm + +from rocket.units import CONVERSIONS, Quantity + +# Type variables for generic analysis +T_Input = TypeVar("T_Input") # Input dataclass type +T_Output = TypeVar("T_Output") # Output type (can be dataclass or tuple) + + +# ============================================================================= +# Parameter Range Specifications +# ============================================================================= + + +@beartype +@dataclass(frozen=True, slots=True) +class Range: + """Specification for a parameter range in a parametric study. + + Supports both dimensionless parameters and Quantity fields with units. + + Examples: + >>> Range(5, 15, n=11, unit="MPa") # 5-15 MPa in 11 steps + >>> Range(2.0, 3.5, n=5) # Dimensionless parameter + >>> Range(values=[2.5, 2.7, 3.0, 3.2]) # Explicit values + """ + + start: float | int | None = None + stop: float | int | None = None + n: int = 10 + unit: str | None = None + values: Sequence[float | int] | None = None + + def __post_init__(self) -> None: + """Validate range specification.""" + if self.values is not None: + if self.start is not None or self.stop is not None: + raise ValueError("Cannot specify both values and start/stop") + else: + if self.start is None or self.stop is None: + raise ValueError("Must specify either values or start/stop") + + def generate(self) -> NDArray[np.float64]: + """Generate array of parameter values.""" + if self.values is not None: + return np.array(self.values, dtype=np.float64) + return np.linspace(self.start, self.stop, self.n) + + def to_quantities(self, dimension: str) -> list[Quantity]: + """Convert range values to Quantity objects. + + Args: + dimension: The dimension of the quantity (e.g., "pressure") + + Returns: + List of Quantity objects + """ + values = self.generate() + if self.unit is None: + raise ValueError(f"Unit required to convert to Quantity for dimension {dimension}") + + return [Quantity(float(v), self.unit, dimension) for v in values] + + +@beartype +@dataclass(frozen=True, slots=True) +class Distribution: + """Base class for probability distributions in uncertainty analysis.""" + + pass + + +@beartype +@dataclass(frozen=True, slots=True) +class Normal(Distribution): + """Normal (Gaussian) distribution. + + Args: + mean: Distribution mean + std: Standard deviation + unit: Optional unit for Quantity fields + """ + + mean: float | int + std: float | int + unit: str | None = None + + def sample(self, n: int, rng: np.random.Generator) -> NDArray[np.float64]: + """Generate n samples from the distribution.""" + return rng.normal(self.mean, self.std, n) + + +@beartype +@dataclass(frozen=True, slots=True) +class Uniform(Distribution): + """Uniform distribution. + + Args: + low: Lower bound + high: Upper bound + unit: Optional unit for Quantity fields + """ + + low: float | int + high: float | int + unit: str | None = None + + def sample(self, n: int, rng: np.random.Generator) -> NDArray[np.float64]: + """Generate n samples from the distribution.""" + return rng.uniform(self.low, self.high, n) + + +@beartype +@dataclass(frozen=True, slots=True) +class Triangular(Distribution): + """Triangular distribution. + + Args: + low: Lower bound + mode: Most likely value + high: Upper bound + unit: Optional unit for Quantity fields + """ + + low: float | int + mode: float | int + high: float | int + unit: str | None = None + + def sample(self, n: int, rng: np.random.Generator) -> NDArray[np.float64]: + """Generate n samples from the distribution.""" + return rng.triangular(self.low, self.mode, self.high, n) + + +@beartype +@dataclass(frozen=True, slots=True) +class LogNormal(Distribution): + """Log-normal distribution. + + Args: + mean: Mean of the underlying normal distribution + sigma: Standard deviation of the underlying normal distribution + unit: Optional unit for Quantity fields + """ + + mean: float | int + sigma: float | int + unit: str | None = None + + def sample(self, n: int, rng: np.random.Generator) -> NDArray[np.float64]: + """Generate n samples from the distribution.""" + return rng.lognormal(self.mean, self.sigma, n) + + +# ============================================================================= +# Introspection Utilities +# ============================================================================= + + +def _get_dataclass_fields(obj: Any) -> dict[str, dataclasses.Field]: + """Get all fields from a dataclass, including nested ones.""" + if not is_dataclass(obj): + raise TypeError(f"Expected dataclass, got {type(obj).__name__}") + return {f.name: f for f in fields(obj)} + + +def _get_field_info(base: Any, field_name: str) -> tuple[Any, str | None]: + """Get the current value and dimension (if Quantity) of a field. + + Args: + base: The base dataclass instance + field_name: Name of the field to inspect + + Returns: + Tuple of (current_value, dimension_or_none) + """ + if not hasattr(base, field_name): + raise ValueError(f"Field '{field_name}' not found in {type(base).__name__}") + + value = getattr(base, field_name) + + if isinstance(value, Quantity): + return value, value.dimension + return value, None + + +def _create_modified_input( + base: T_Input, + field_name: str, + value: float | Quantity, + original_dimension: str | None, +) -> T_Input: + """Create a modified copy of the input with one field changed. + + Handles both Quantity and plain numeric fields. + """ + current_value = getattr(base, field_name) + + if isinstance(current_value, Quantity): + # Field is a Quantity - ensure we create a proper Quantity + if isinstance(value, Quantity): + new_value = value + else: + # Value is numeric, need unit from original + new_value = Quantity(float(value), current_value.unit, current_value.dimension) + else: + # Field is a plain numeric type + new_value = value + + return replace(base, **{field_name: new_value}) + + +def _extract_metrics(result: Any, prefix: str = "") -> dict[str, float]: + """Recursively extract all numeric values from a result. + + Handles dataclasses, tuples, and nested structures. + Returns flat dict with keys for nested values. + + For tuples of dataclasses (common pattern like (Performance, Geometry)), + fields are extracted without prefixes to keep names clean. + """ + metrics: dict[str, float] = {} + + if isinstance(result, tuple): + # Handle tuple of results (e.g., (performance, geometry)) + # Extract fields directly without prefixes for cleaner column names + for item in result: + metrics.update(_extract_metrics(item, prefix)) + + elif is_dataclass(result) and not isinstance(result, type): + # Handle dataclass + for field in fields(result): + field_value = getattr(result, field.name) + field_key = f"{prefix}{field.name}" if prefix else field.name + + if isinstance(field_value, Quantity): + metrics[field_key] = float(field_value.value) + elif isinstance(field_value, (int, float)): + metrics[field_key] = float(field_value) + elif is_dataclass(field_value): + metrics.update(_extract_metrics(field_value, f"{field_key}.")) + + return metrics + + +# ============================================================================= +# Study Results +# ============================================================================= + + +@beartype +@dataclass +class StudyResults: + """Results from a parametric study or uncertainty analysis. + + Contains all input combinations, computed outputs, and extracted metrics. + Provides methods for plotting, filtering, and export. + + Attributes: + inputs: List of input parameter combinations + outputs: List of computed results + metrics: Dict mapping metric names to arrays of values + parameters: Dict mapping parameter names to arrays of values + constraints_passed: Boolean array indicating which runs passed constraints + """ + + inputs: list[Any] + outputs: list[Any] + metrics: dict[str, NDArray[np.float64]] + parameters: dict[str, NDArray[np.float64]] + constraints_passed: NDArray[np.bool_] | None = None + + @property + def n_runs(self) -> int: + """Number of runs in the study.""" + return len(self.inputs) + + @property + def n_feasible(self) -> int: + """Number of runs that passed all constraints.""" + if self.constraints_passed is None: + return self.n_runs + return int(np.sum(self.constraints_passed)) + + def get_metric(self, name: str, feasible_only: bool = False) -> NDArray[np.float64]: + """Get values for a specific metric. + + Args: + name: Metric name (e.g., "isp", "throat_diameter") + feasible_only: If True, only return values where constraints passed + + Returns: + Array of metric values + """ + if name not in self.metrics: + available = list(self.metrics.keys()) + raise ValueError(f"Unknown metric '{name}'. Available: {available}") + + values = self.metrics[name] + if feasible_only and self.constraints_passed is not None: + return values[self.constraints_passed] + return values + + def get_parameter(self, name: str, feasible_only: bool = False) -> NDArray[np.float64]: + """Get values for a specific input parameter. + + Args: + name: Parameter name (e.g., "chamber_pressure") + feasible_only: If True, only return values where constraints passed + + Returns: + Array of parameter values + """ + if name not in self.parameters: + available = list(self.parameters.keys()) + raise ValueError(f"Unknown parameter '{name}'. Available: {available}") + + values = self.parameters[name] + if feasible_only and self.constraints_passed is not None: + return values[self.constraints_passed] + return values + + def get_best( + self, + metric: str, + maximize: bool = True, + feasible_only: bool = True, + ) -> tuple[Any, Any, float]: + """Get the best run according to a metric. + + Args: + metric: Metric to optimize + maximize: If True, find maximum; if False, find minimum + feasible_only: Only consider runs that passed constraints + + Returns: + Tuple of (best_input, best_output, best_metric_value) + """ + values = self.get_metric(metric, feasible_only=False) + mask = self.constraints_passed if feasible_only and self.constraints_passed is not None else np.ones(len(values), dtype=bool) + + if not np.any(mask): + raise ValueError("No feasible solutions found") + + masked_values = np.where(mask, values, -np.inf if maximize else np.inf) + best_idx = int(np.argmax(masked_values) if maximize else np.argmin(masked_values)) + + return self.inputs[best_idx], self.outputs[best_idx], float(values[best_idx]) + + def to_dataframe(self) -> pl.DataFrame: + """Export results to a Polars DataFrame. + + Returns: + Polars DataFrame with parameters and metrics + """ + data = {**self.parameters, **self.metrics} + if self.constraints_passed is not None: + data["feasible"] = self.constraints_passed + return pl.DataFrame(data) + + def to_csv(self, path: str | Path) -> None: + """Export results to CSV file. + + Args: + path: Output file path + """ + df = self.to_dataframe() + df.write_csv(path) + + def list_metrics(self) -> list[str]: + """List all available metric names.""" + return list(self.metrics.keys()) + + def list_parameters(self) -> list[str]: + """List all varied parameter names.""" + return list(self.parameters.keys()) + + +# ============================================================================= +# Parametric Study +# ============================================================================= + + +@beartype +class ParametricStudy(Generic[T_Input, T_Output]): + """General-purpose parametric study framework. + + Runs a computation over a grid of parameter variations, automatically + discovering valid parameters through dataclass introspection. + + This design is non-brittle: + - Adding new fields to input dataclasses automatically makes them available + - No hardcoded parameter names + - Works with any frozen dataclass + computation function + + Example: + >>> study = ParametricStudy( + ... compute=design_engine, + ... base=inputs, + ... vary={ + ... "chamber_pressure": Range(5, 15, n=11, unit="MPa"), + ... "mixture_ratio": Range(2.5, 3.5, n=5), + ... }, + ... ) + >>> results = study.run() + """ + + def __init__( + self, + compute: Callable[[T_Input], T_Output], + base: T_Input, + vary: dict[str, Range | Sequence[Any]], + constraints: list[Callable[[T_Output], bool]] | None = None, + ) -> None: + """Initialize parametric study. + + Args: + compute: Function that takes input and returns output + base: Base input dataclass with default values + vary: Dict mapping field names to Range specifications or plain sequences. + Use Range for unit-aware sweeps, or a plain list for discrete values. + constraints: Optional list of constraint functions. + Each takes output and returns True if feasible. + """ + self.compute = compute + self.base = base + self.vary = vary + self.constraints = constraints or [] + + # Validate that all varied parameters exist in the base dataclass + self._validate_parameters() + + def _validate_parameters(self) -> None: + """Validate that all varied parameters exist and have compatible types.""" + valid_fields = _get_dataclass_fields(self.base) + + for param_name, param_spec in self.vary.items(): + if param_name not in valid_fields: + raise ValueError( + f"Parameter '{param_name}' not found in {type(self.base).__name__}. " + f"Valid fields: {list(valid_fields.keys())}" + ) + + # Check unit compatibility for Quantity fields (only if Range is used) + current_value = getattr(self.base, param_name) + if isinstance(current_value, Quantity) and isinstance(param_spec, Range): + if param_spec.unit is None: + raise ValueError( + f"Parameter '{param_name}' is a Quantity, but no unit specified in Range. " + f"Current unit: {current_value.unit}" + ) + # Verify unit is valid and has compatible dimension + if param_spec.unit not in CONVERSIONS: + raise ValueError(f"Unknown unit '{param_spec.unit}' for parameter '{param_name}'") + + range_dim = CONVERSIONS[param_spec.unit][1] + if range_dim != current_value.dimension: + raise ValueError( + f"Unit '{param_spec.unit}' has dimension '{range_dim}', " + f"but field '{param_name}' has dimension '{current_value.dimension}'" + ) + + def _generate_grid(self) -> list[dict[str, float | Quantity]]: + """Generate all parameter combinations.""" + # Generate values for each parameter + param_values: dict[str, list[Any]] = {} + + for param_name, param_spec in self.vary.items(): + current_value = getattr(self.base, param_name) + + if isinstance(param_spec, Range): + # Range specification + if isinstance(current_value, Quantity): + # Generate Quantity values + param_values[param_name] = param_spec.to_quantities(current_value.dimension) + else: + # Generate plain numeric values + param_values[param_name] = list(param_spec.generate()) + else: + # Plain sequence - use values as-is + param_values[param_name] = list(param_spec) + + # Generate all combinations + keys = list(param_values.keys()) + value_lists = [param_values[k] for k in keys] + combinations = list(itertools.product(*value_lists)) + + return [dict(zip(keys, combo, strict=True)) for combo in combinations] + + def run(self, progress: bool = False) -> StudyResults: + """Run the parametric study. + + Args: + progress: If True, print progress (requires tqdm for fancy progress bar) + + Returns: + StudyResults containing all inputs, outputs, and extracted metrics + """ + grid = self._generate_grid() + n_total = len(grid) + + inputs_list: list[T_Input] = [] + outputs_list: list[T_Output] = [] + all_metrics: list[dict[str, float]] = [] + all_params: list[dict[str, float]] = [] + constraints_passed: list[bool] = [] + + # Create iterator with optional progress bar + iterator: Any = grid + if progress: + iterator = tqdm(grid, desc="Running study", total=n_total) + + for i, param_combo in enumerate(iterator): + # Create modified input + modified_input = self.base + for param_name, param_value in param_combo.items(): + modified_input = _create_modified_input( + modified_input, + param_name, + param_value, + current_value.dimension if isinstance((current_value := getattr(self.base, param_name)), Quantity) else None, + ) + + # Run computation + try: + output = self.compute(modified_input) + success = True + except Exception as e: + # Store None for failed runs + output = None # type: ignore + success = False + if progress and not hasattr(iterator, 'set_postfix'): + print(f" Run {i+1}/{n_total} failed: {e}") + + inputs_list.append(modified_input) + outputs_list.append(output) + + # Extract metrics + if success and output is not None: + metrics = _extract_metrics(output) + all_metrics.append(metrics) + + # Check constraints + passed = all(constraint(output) for constraint in self.constraints) + else: + all_metrics.append({}) + passed = False + + constraints_passed.append(passed) + + # Extract parameter values (numeric form for plotting) + param_dict: dict[str, float] = {} + for param_name, param_value in param_combo.items(): + if isinstance(param_value, Quantity): + param_dict[param_name] = float(param_value.value) + else: + param_dict[param_name] = float(param_value) + all_params.append(param_dict) + + # Consolidate metrics into arrays + if all_metrics and all_metrics[0]: + metric_names = set() + for m in all_metrics: + metric_names.update(m.keys()) + + metrics_arrays = { + name: np.array([m.get(name, np.nan) for m in all_metrics]) + for name in metric_names + } + else: + metrics_arrays = {} + + # Consolidate parameters into arrays + if all_params: + param_names = list(all_params[0].keys()) + params_arrays = { + name: np.array([p[name] for p in all_params]) + for name in param_names + } + else: + params_arrays = {} + + return StudyResults( + inputs=inputs_list, + outputs=outputs_list, + metrics=metrics_arrays, + parameters=params_arrays, + constraints_passed=np.array(constraints_passed), + ) + + +# ============================================================================= +# Uncertainty Analysis +# ============================================================================= + + +@beartype +class UncertaintyAnalysis(Generic[T_Input, T_Output]): + """Monte Carlo uncertainty quantification. + + Samples input parameters from specified distributions and propagates + uncertainty through the computation. + + Example: + >>> analysis = UncertaintyAnalysis( + ... compute=design_engine, + ... base=inputs, + ... distributions={ + ... "gamma": Normal(1.22, 0.02), + ... "chamber_temp": Normal(3200, 50, unit="K"), + ... }, + ... ) + >>> results = analysis.run(n_samples=1000) + >>> print(f"Isp = {results.mean('isp'):.1f} ± {results.std('isp'):.1f} s") + """ + + def __init__( + self, + compute: Callable[[T_Input], T_Output], + base: T_Input, + distributions: dict[str, Distribution], + constraints: list[Callable[[T_Output], bool]] | None = None, + seed: int | None = None, + ) -> None: + """Initialize uncertainty analysis. + + Args: + compute: Function that takes input and returns output + base: Base input dataclass with nominal values + distributions: Dict mapping field names to Distribution specifications + constraints: Optional constraint functions + seed: Random seed for reproducibility + """ + self.compute = compute + self.base = base + self.distributions = distributions + self.constraints = constraints or [] + self.rng = np.random.default_rng(seed) + + self._validate_parameters() + + def _validate_parameters(self) -> None: + """Validate that all uncertain parameters exist.""" + valid_fields = _get_dataclass_fields(self.base) + + for param_name, dist in self.distributions.items(): + if param_name not in valid_fields: + raise ValueError( + f"Parameter '{param_name}' not found in {type(self.base).__name__}. " + f"Valid fields: {list(valid_fields.keys())}" + ) + + # Check unit compatibility for Quantity fields + current_value = getattr(self.base, param_name) + if isinstance(current_value, Quantity) and (not hasattr(dist, 'unit') or dist.unit is None): # type: ignore + raise ValueError( + f"Parameter '{param_name}' is a Quantity, but no unit specified in Distribution" + ) + + def run(self, n_samples: int = 1000, progress: bool = False) -> "UncertaintyResults": + """Run Monte Carlo uncertainty analysis. + + Args: + n_samples: Number of Monte Carlo samples + progress: If True, show progress indicator + + Returns: + UncertaintyResults with statistics and samples + """ + # Generate all samples upfront + samples: dict[str, NDArray[np.float64]] = {} + for param_name, dist in self.distributions.items(): + samples[param_name] = dist.sample(n_samples, self.rng) + + inputs_list: list[T_Input] = [] + outputs_list: list[T_Output] = [] + all_metrics: list[dict[str, float]] = [] + constraints_passed: list[bool] = [] + + iterator: Any = range(n_samples) + if progress: + iterator = tqdm(range(n_samples), desc="Sampling") + + for i in iterator: + # Create modified input with sampled values + modified_input = self.base + + for param_name, dist in self.distributions.items(): + sampled_value = samples[param_name][i] + current_value = getattr(self.base, param_name) + + if isinstance(current_value, Quantity): + # Create Quantity with sampled value + unit = dist.unit # type: ignore + new_value = Quantity(float(sampled_value), unit, current_value.dimension) + else: + new_value = float(sampled_value) + + modified_input = replace(modified_input, **{param_name: new_value}) + + # Run computation + try: + output = self.compute(modified_input) + success = True + except Exception: + output = None # type: ignore + success = False + + inputs_list.append(modified_input) + outputs_list.append(output) + + if success and output is not None: + metrics = _extract_metrics(output) + all_metrics.append(metrics) + passed = all(constraint(output) for constraint in self.constraints) + else: + all_metrics.append({}) + passed = False + + constraints_passed.append(passed) + + # Consolidate metrics + if all_metrics and all_metrics[0]: + metric_names = set() + for m in all_metrics: + metric_names.update(m.keys()) + + metrics_arrays = { + name: np.array([m.get(name, np.nan) for m in all_metrics]) + for name in metric_names + } + else: + metrics_arrays = {} + + return UncertaintyResults( + inputs=inputs_list, + outputs=outputs_list, + metrics=metrics_arrays, + samples=samples, + constraints_passed=np.array(constraints_passed), + n_samples=n_samples, + ) + + +@beartype +@dataclass +class UncertaintyResults: + """Results from uncertainty analysis. + + Provides statistical summaries and access to all samples. + """ + + inputs: list[Any] + outputs: list[Any] + metrics: dict[str, NDArray[np.float64]] + samples: dict[str, NDArray[np.float64]] + constraints_passed: NDArray[np.bool_] + n_samples: int + + def mean(self, metric: str, feasible_only: bool = False) -> float: + """Get mean value of a metric.""" + values = self._get_values(metric, feasible_only) + return float(np.nanmean(values)) + + def std(self, metric: str, feasible_only: bool = False) -> float: + """Get standard deviation of a metric.""" + values = self._get_values(metric, feasible_only) + return float(np.nanstd(values)) + + def percentile( + self, metric: str, p: float | Sequence[float], feasible_only: bool = False + ) -> float | NDArray[np.float64]: + """Get percentile(s) of a metric. + + Args: + metric: Metric name + p: Percentile(s) to compute (0-100) + feasible_only: Only use feasible samples + + Returns: + Percentile value(s) + """ + values = self._get_values(metric, feasible_only) + result = np.nanpercentile(values, p) + if isinstance(p, (int, float)): + return float(result) + return result + + def confidence_interval( + self, metric: str, confidence: float = 0.95, feasible_only: bool = False + ) -> tuple[float, float]: + """Get confidence interval for a metric. + + Args: + metric: Metric name + confidence: Confidence level (0-1), default 0.95 for 95% CI + feasible_only: Only use feasible samples + + Returns: + Tuple of (lower_bound, upper_bound) + """ + alpha = 1 - confidence + lower_p = alpha / 2 * 100 + upper_p = (1 - alpha / 2) * 100 + + values = self._get_values(metric, feasible_only) + return ( + float(np.nanpercentile(values, lower_p)), + float(np.nanpercentile(values, upper_p)), + ) + + def probability_of_success(self) -> float: + """Get fraction of samples that passed all constraints.""" + return float(np.mean(self.constraints_passed)) + + def _get_values(self, metric: str, feasible_only: bool) -> NDArray[np.float64]: + """Get metric values, optionally filtered to feasible only.""" + if metric not in self.metrics: + available = list(self.metrics.keys()) + raise ValueError(f"Unknown metric '{metric}'. Available: {available}") + + values = self.metrics[metric] + if feasible_only: + values = values[self.constraints_passed] + return values + + def summary(self, metrics: list[str] | None = None) -> str: + """Generate a text summary of uncertainty results. + + Args: + metrics: List of metrics to summarize. If None, summarizes all. + + Returns: + Formatted string summary + """ + if metrics is None: + metrics = [m for m in self.metrics if not m.endswith("_si")] + + lines = [ + "Uncertainty Analysis Results", + "=" * 50, + f"Samples: {self.n_samples}", + f"Feasible: {np.sum(self.constraints_passed)} ({self.probability_of_success()*100:.1f}%)", + "", + f"{'Metric':<25} {'Mean':>12} {'Std':>12} {'95% CI':>20}", + "-" * 70, + ] + + for metric in metrics: + if metric in self.metrics: + mean = self.mean(metric) + std = self.std(metric) + ci = self.confidence_interval(metric) + lines.append( + f"{metric:<25} {mean:>12.4g} {std:>12.4g} [{ci[0]:.4g}, {ci[1]:.4g}]" + ) + + return "\n".join(lines) + + def to_dataframe(self) -> pl.DataFrame: + """Export results to Polars DataFrame.""" + data = {**self.samples, **self.metrics, "feasible": self.constraints_passed} + return pl.DataFrame(data) + + def to_csv(self, path: str | Path) -> None: + """Export results to CSV file. + + Args: + path: Output file path + """ + df = self.to_dataframe() + df.write_csv(path) + + +# ============================================================================= +# Multi-Objective Optimization +# ============================================================================= + + +@beartype +def compute_pareto_front( + objectives: NDArray[np.float64], + maximize: Sequence[bool], +) -> NDArray[np.bool_]: + """Identify Pareto-optimal points in a set of objectives. + + A point is Pareto-optimal if no other point dominates it (i.e., no + other point is better in all objectives simultaneously). + + Args: + objectives: Array of shape (n_points, n_objectives) + maximize: List of booleans indicating whether to maximize each objective + + Returns: + Boolean array indicating which points are Pareto-optimal + """ + n_points = objectives.shape[0] + is_pareto = np.ones(n_points, dtype=bool) + + # Flip signs for maximization (we'll minimize internally) + obj = objectives.copy() + for i, is_max in enumerate(maximize): + if is_max: + obj[:, i] = -obj[:, i] + + for i in range(n_points): + if not is_pareto[i]: + continue + + for j in range(n_points): + if i == j or not is_pareto[j]: + continue + + # j dominates i if j is <= in all objectives and < in at least one + if np.all(obj[j] <= obj[i]) and np.any(obj[j] < obj[i]): + is_pareto[i] = False + break + + return is_pareto + + +@beartype +class MultiObjectiveOptimizer(Generic[T_Input, T_Output]): + """Multi-objective optimizer for finding Pareto-optimal designs. + + Uses a combination of grid search and local refinement to find + designs on the Pareto frontier. + + Example: + >>> optimizer = MultiObjectiveOptimizer( + ... compute=design_engine, + ... base=inputs, + ... objectives=["isp_vac", "thrust_to_weight"], + ... maximize=[True, True], + ... vary={"chamber_pressure": Range(5, 20, n=10, unit="MPa")}, + ... ) + >>> pareto_results = optimizer.run() + """ + + def __init__( + self, + compute: Callable[[T_Input], T_Output], + base: T_Input, + objectives: list[str], + maximize: list[bool], + vary: dict[str, Range | Sequence[Any]], + constraints: list[Callable[[T_Output], bool]] | None = None, + ) -> None: + """Initialize multi-objective optimizer. + + Args: + compute: Function that takes input and returns output + base: Base input dataclass + objectives: List of metric names to optimize + maximize: List of booleans for each objective (True = maximize) + vary: Dict mapping field names to Range specifications or plain sequences + constraints: Optional constraint functions + """ + self.compute = compute + self.base = base + self.objectives = objectives + self.maximize = maximize + self.vary = vary + self.constraints = constraints or [] + + if len(objectives) != len(maximize): + raise ValueError("objectives and maximize must have same length") + + def run(self, progress: bool = False) -> "ParetoResults": + """Run the multi-objective optimization. + + First performs a parametric sweep, then identifies Pareto-optimal points. + + Args: + progress: If True, show progress indicator + + Returns: + ParetoResults with Pareto-optimal designs + """ + # Run parametric study + study = ParametricStudy( + compute=self.compute, + base=self.base, + vary=self.vary, + constraints=self.constraints, + ) + results = study.run(progress=progress) + + # Extract objectives + obj_arrays = [] + for obj_name in self.objectives: + if obj_name not in results.metrics: + raise ValueError(f"Objective '{obj_name}' not found in results") + obj_arrays.append(results.get_metric(obj_name)) + + objectives_matrix = np.column_stack(obj_arrays) + + # Filter to feasible points + if results.constraints_passed is not None: + feasible_mask = results.constraints_passed + else: + feasible_mask = np.ones(len(objectives_matrix), dtype=bool) + + # Remove NaN values + valid_mask = feasible_mask & ~np.any(np.isnan(objectives_matrix), axis=1) + valid_indices = np.where(valid_mask)[0] + + if len(valid_indices) == 0: + return ParetoResults( + all_results=results, + pareto_indices=[], + pareto_inputs=[], + pareto_outputs=[], + pareto_objectives=np.array([]).reshape(0, len(self.objectives)), + objective_names=self.objectives, + maximize=self.maximize, + ) + + # Compute Pareto front on valid points only + valid_objectives = objectives_matrix[valid_mask] + is_pareto = compute_pareto_front(valid_objectives, self.maximize) + + # Map back to original indices + pareto_indices = valid_indices[is_pareto].tolist() + pareto_inputs = [results.inputs[i] for i in pareto_indices] + pareto_outputs = [results.outputs[i] for i in pareto_indices] + pareto_objectives = valid_objectives[is_pareto] + + return ParetoResults( + all_results=results, + pareto_indices=pareto_indices, + pareto_inputs=pareto_inputs, + pareto_outputs=pareto_outputs, + pareto_objectives=pareto_objectives, + objective_names=self.objectives, + maximize=self.maximize, + ) + + +@beartype +@dataclass +class ParetoResults: + """Results from multi-objective optimization. + + Contains the Pareto-optimal designs and full study results. + """ + + all_results: StudyResults + pareto_indices: list[int] + pareto_inputs: list[Any] + pareto_outputs: list[Any] + pareto_objectives: NDArray[np.float64] + objective_names: list[str] + maximize: list[bool] + + @property + def n_pareto(self) -> int: + """Number of Pareto-optimal points.""" + return len(self.pareto_indices) + + @property + def n_total(self) -> int: + """Total number of designs evaluated.""" + return len(self.all_results.inputs) + + @property + def n_feasible(self) -> int: + """Number of designs that passed all constraints.""" + return self.all_results.n_feasible + + def pareto_front(self) -> pl.DataFrame: + """Get Pareto-optimal designs as a DataFrame. + + Returns: + DataFrame with parameters and metrics for Pareto-optimal points + """ + # Extract metrics for Pareto points only + pareto_metrics: dict[str, list[Any]] = {} + + # Add parameters + for param_name, param_values in self.all_results.parameters.items(): + pareto_metrics[param_name] = [param_values[i] for i in self.pareto_indices] + + # Add metrics + for metric_name, metric_values in self.all_results.metrics.items(): + pareto_metrics[metric_name] = [metric_values[i] for i in self.pareto_indices] + + return pl.DataFrame(pareto_metrics) + + def get_best(self, objective: str) -> tuple[Any, Any, float]: + """Get the best design for a specific objective. + + Args: + objective: Name of objective to optimize + + Returns: + Tuple of (input, output, objective_value) + """ + if objective not in self.objective_names: + raise ValueError(f"Unknown objective: {objective}") + + obj_idx = self.objective_names.index(objective) + values = self.pareto_objectives[:, obj_idx] + + best_idx = int(np.argmax(values)) if self.maximize[obj_idx] else int(np.argmin(values)) + + return ( + self.pareto_inputs[best_idx], + self.pareto_outputs[best_idx], + float(values[best_idx]), + ) + + def get_compromise(self, weights: list[float] | None = None) -> tuple[Any, Any]: + """Get a compromise solution from the Pareto front. + + Uses weighted sum of normalized objectives. + + Args: + weights: Weights for each objective (default: equal weights) + + Returns: + Tuple of (input, output) for the compromise solution + """ + if weights is None: + weights = [1.0 / len(self.objectives) for _ in self.objective_names] + + if len(weights) != len(self.objective_names): + raise ValueError("weights must have same length as objectives") + + # Normalize objectives to [0, 1] + obj_norm = self.pareto_objectives.copy() + for i in range(obj_norm.shape[1]): + col = obj_norm[:, i] + col_min, col_max = np.min(col), np.max(col) + if col_max > col_min: + obj_norm[:, i] = (col - col_min) / (col_max - col_min) + else: + obj_norm[:, i] = 0.5 + + # Flip if minimizing + if not self.maximize[i]: + obj_norm[:, i] = 1 - obj_norm[:, i] + + # Weighted sum + scores = np.sum(obj_norm * np.array(weights), axis=1) + best_idx = int(np.argmax(scores)) + + return self.pareto_inputs[best_idx], self.pareto_outputs[best_idx] + + def summary(self) -> str: + """Generate text summary of Pareto results.""" + lines = [ + "Multi-Objective Optimization Results", + "=" * 50, + f"Total designs evaluated: {self.all_results.n_runs}", + f"Feasible designs: {self.all_results.n_feasible}", + f"Pareto-optimal designs: {self.n_pareto}", + "", + "Objectives:", + ] + + for i, (name, is_max) in enumerate(zip(self.objective_names, self.maximize, strict=True)): + direction = "maximize" if is_max else "minimize" + if self.n_pareto > 0: + values = self.pareto_objectives[:, i] + lines.append( + f" {name} ({direction}): " + f"range [{np.min(values):.4g}, {np.max(values):.4g}]" + ) + else: + lines.append(f" {name} ({direction}): no feasible points") + + return "\n".join(lines) + diff --git a/rocket/cycles/__init__.py b/rocket/cycles/__init__.py new file mode 100644 index 0000000..860ba89 --- /dev/null +++ b/rocket/cycles/__init__.py @@ -0,0 +1,55 @@ +"""Engine cycle analysis module for Rocket. + +This module provides analysis tools for different rocket engine cycles: +- Pressure-fed (simplest) +- Gas generator (turbopump-fed with separate combustion) +- Expander (turbine driven by heated propellant) +- Staged combustion (preburner exhaust into main chamber) + +Each cycle type has different performance characteristics, complexity, +and feasibility constraints. + +Example: + >>> from rocket import EngineInputs, design_engine + >>> from rocket.cycles import GasGeneratorCycle, analyze_cycle + >>> + >>> inputs = EngineInputs.from_propellants("LOX", "CH4", ...) + >>> performance, geometry = design_engine(inputs) + >>> + >>> cycle = GasGeneratorCycle( + ... turbine_inlet_temp=kelvin(900), + ... pump_efficiency=0.70, + ... ) + >>> result = analyze_cycle(inputs, performance, geometry, cycle) + >>> print(f"Net Isp: {result.net_isp.value:.1f} s") +""" + +from rocket.cycles.base import ( + CycleConfiguration, + CyclePerformance, + CycleType, + analyze_cycle, + format_cycle_summary, +) +from rocket.cycles.gas_generator import GasGeneratorCycle +from rocket.cycles.pressure_fed import PressureFedCycle +from rocket.cycles.staged_combustion import ( + FullFlowStagedCombustion, + StagedCombustionCycle, +) + +__all__ = [ + # Base types + "CycleConfiguration", + "CyclePerformance", + "CycleType", + # Cycle configurations + "PressureFedCycle", + "GasGeneratorCycle", + "StagedCombustionCycle", + "FullFlowStagedCombustion", + # Analysis function + "analyze_cycle", + "format_cycle_summary", +] + diff --git a/rocket/cycles/base.py b/rocket/cycles/base.py new file mode 100644 index 0000000..fd27b26 --- /dev/null +++ b/rocket/cycles/base.py @@ -0,0 +1,351 @@ +"""Base types for engine cycle analysis. + +This module defines the common interfaces and data structures used by +all engine cycle types. +""" + +import math +from dataclasses import dataclass +from enum import Enum, auto +from typing import Protocol, runtime_checkable + +from beartype import beartype + +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.units import Quantity, pascals + + +class CycleType(Enum): + """Engine cycle type enumeration.""" + + PRESSURE_FED = auto() + GAS_GENERATOR = auto() + EXPANDER = auto() + STAGED_COMBUSTION = auto() + FULL_FLOW_STAGED = auto() + ELECTRIC_PUMP = auto() + + +@beartype +@dataclass(frozen=True, slots=True) +class CyclePerformance: + """Performance results from engine cycle analysis. + + Captures the system-level performance including losses from + turbine drive systems, pump power requirements, and pressure margins. + + Attributes: + net_isp: Effective Isp after accounting for cycle losses [s] + net_thrust: Delivered thrust after losses [N] + cycle_efficiency: Ratio of net Isp to ideal Isp [-] + pump_power_ox: Oxidizer pump power requirement [W] + pump_power_fuel: Fuel pump power requirement [W] + turbine_power: Total turbine power available [W] + turbine_mass_flow: Mass flow through turbine [kg/s] + tank_pressure_ox: Required oxidizer tank pressure [Pa] + tank_pressure_fuel: Required fuel tank pressure [Pa] + npsh_margin_ox: Net Positive Suction Head margin for ox pump [Pa] + npsh_margin_fuel: NPSH margin for fuel pump [Pa] + cycle_type: Type of cycle analyzed + feasible: Whether the cycle closes (power balance satisfied) + warnings: List of any warnings or marginal conditions + """ + + net_isp: Quantity + net_thrust: Quantity + cycle_efficiency: float + pump_power_ox: Quantity + pump_power_fuel: Quantity + turbine_power: Quantity + turbine_mass_flow: Quantity + tank_pressure_ox: Quantity + tank_pressure_fuel: Quantity + npsh_margin_ox: Quantity + npsh_margin_fuel: Quantity + cycle_type: CycleType + feasible: bool + warnings: list[str] + + +@runtime_checkable +class CycleConfiguration(Protocol): + """Protocol that all cycle configurations must implement. + + This protocol ensures that any cycle configuration can be used + with the generic analyze_cycle() function. + """ + + @property + def cycle_type(self) -> CycleType: + """Return the cycle type.""" + ... + + def analyze( + self, + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + ) -> CyclePerformance: + """Analyze the cycle and return performance results. + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + + Returns: + CyclePerformance with system-level results + """ + ... + + +@beartype +def analyze_cycle( + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + cycle: CycleConfiguration, +) -> CyclePerformance: + """Analyze an engine cycle configuration. + + This is the main entry point for cycle analysis. It delegates to + the specific cycle implementation's analyze() method. + + Args: + inputs: Engine input parameters + performance: Computed engine performance (from design_engine) + geometry: Computed engine geometry (from design_engine) + cycle: Cycle configuration (PressureFedCycle, GasGeneratorCycle, etc.) + + Returns: + CyclePerformance with net performance and feasibility assessment + + Example: + >>> inputs = EngineInputs.from_propellants("LOX", "CH4", ...) + >>> performance, geometry = design_engine(inputs) + >>> cycle = GasGeneratorCycle(turbine_inlet_temp=kelvin(900), ...) + >>> result = analyze_cycle(inputs, performance, geometry, cycle) + """ + return cycle.analyze(inputs, performance, geometry) + + +# ============================================================================= +# Utility Functions +# ============================================================================= + + +@beartype +def pump_power( + mass_flow: Quantity, + pressure_rise: Quantity, + density: float, + efficiency: float, +) -> Quantity: + """Calculate pump power requirement. + + Uses the basic hydraulic power equation: + P = (mdot * delta_P) / (rho * eta) + + Args: + mass_flow: Mass flow rate through pump [kg/s] + pressure_rise: Pressure rise across pump [Pa] + density: Fluid density [kg/m³] + efficiency: Pump efficiency (0-1) + + Returns: + Pump power in Watts + """ + mdot = mass_flow.to("kg/s").value + dp = pressure_rise.to("Pa").value + + # Volumetric flow rate + Q = mdot / density # m³/s + + # Hydraulic power + P_hydraulic = Q * dp # W + + # Shaft power (accounting for efficiency) + P_shaft = P_hydraulic / efficiency + + return Quantity(P_shaft, "W", "power") + + +@beartype +def turbine_power( + mass_flow: Quantity, + inlet_temp: Quantity, + pressure_ratio: float, + gamma: float, + efficiency: float, + R: float = 287.0, # J/(kg·K), approximate for combustion products +) -> Quantity: + """Calculate turbine power output. + + Uses isentropic turbine equations with efficiency factor. + + Args: + mass_flow: Mass flow through turbine [kg/s] + inlet_temp: Turbine inlet temperature [K] + pressure_ratio: Inlet pressure / outlet pressure [-] + gamma: Ratio of specific heats for turbine gas [-] + efficiency: Turbine isentropic efficiency (0-1) + R: Specific gas constant [J/(kg·K)] + + Returns: + Turbine power output in Watts + """ + mdot = mass_flow.to("kg/s").value + T_in = inlet_temp.to("K").value + + # Isentropic temperature ratio + T_ratio_ideal = pressure_ratio ** ((gamma - 1) / gamma) + + # Actual temperature drop + delta_T_ideal = T_in * (1 - 1 / T_ratio_ideal) + delta_T_actual = delta_T_ideal * efficiency + + # Specific heat at constant pressure + cp = gamma * R / (gamma - 1) + + # Turbine power + P = mdot * cp * delta_T_actual + + return Quantity(P, "W", "power") + + +@beartype +def npsh_available( + tank_pressure: Quantity, + fluid_density: float, + vapor_pressure: Quantity, + inlet_height: float = 0.0, + line_losses: Quantity | None = None, +) -> Quantity: + """Calculate Net Positive Suction Head available at pump inlet. + + NPSH_a = (P_tank - P_vapor) / (rho * g) + h - losses + + Args: + tank_pressure: Tank ullage pressure [Pa] + fluid_density: Propellant density [kg/m³] + vapor_pressure: Propellant vapor pressure [Pa] + inlet_height: Height of fluid above pump inlet [m] + line_losses: Pressure losses in feed lines [Pa] + + Returns: + NPSH available in Pascals (pressure equivalent) + """ + g = 9.80665 # m/s² + + P_tank = tank_pressure.to("Pa").value + P_vapor = vapor_pressure.to("Pa").value + losses = line_losses.to("Pa").value if line_losses else 0.0 + + # NPSH in meters of head + npsh_m = (P_tank - P_vapor) / (fluid_density * g) + inlet_height - losses / (fluid_density * g) + + # Convert to pressure equivalent + npsh_pa = npsh_m * fluid_density * g + + return pascals(npsh_pa) + + +@beartype +def estimate_line_losses( + mass_flow: Quantity, + density: float, + pipe_diameter: float, + pipe_length: float, + num_elbows: int = 2, + num_valves: int = 2, +) -> Quantity: + """Estimate pressure losses in feed lines. + + Uses Darcy-Weisbach equation with loss coefficients for fittings. + + Args: + mass_flow: Mass flow rate [kg/s] + density: Fluid density [kg/m³] + pipe_diameter: Pipe inner diameter [m] + pipe_length: Total pipe length [m] + num_elbows: Number of 90° elbows + num_valves: Number of valves + + Returns: + Total pressure loss [Pa] + """ + mdot = mass_flow.to("kg/s").value + D = pipe_diameter + L = pipe_length + + # Calculate velocity + A = math.pi * (D / 2) ** 2 + V = mdot / (density * A) + + # Dynamic pressure + q = 0.5 * density * V ** 2 + + # Friction factor (assuming turbulent flow, smooth pipe) + # Using Blasius correlation as approximation + Re = density * V * D / 1e-3 # Approximate viscosity + f = 0.316 / Re ** 0.25 if Re > 2300 else 64 / Re + + # Pipe friction losses + dp_pipe = f * (L / D) * q + + # Fitting losses (K-factors) + K_elbow = 0.3 # 90° elbow + K_valve = 0.2 # Gate valve (open) + + dp_fittings = (num_elbows * K_elbow + num_valves * K_valve) * q + + return pascals(dp_pipe + dp_fittings) + + +@beartype +def format_cycle_summary(result: CyclePerformance) -> str: + """Format cycle analysis results as readable string. + + Args: + result: CyclePerformance from analyze_cycle() + + Returns: + Formatted multi-line string + """ + status = "✓ FEASIBLE" if result.feasible else "✗ INFEASIBLE" + + lines = [ + f"{'=' * 60}", + f"CYCLE ANALYSIS: {result.cycle_type.name}", + f"Status: {status}", + f"{'=' * 60}", + "", + "PERFORMANCE:", + f" Net Isp: {result.net_isp.value:.1f} s", + f" Net Thrust: {result.net_thrust.to('kN').value:.2f} kN", + f" Cycle Efficiency: {result.cycle_efficiency * 100:.1f}%", + "", + "POWER BALANCE:", + f" Turbine Power: {result.turbine_power.value / 1000:.1f} kW", + f" Pump Power (Ox): {result.pump_power_ox.value / 1000:.1f} kW", + f" Pump Power (Fuel): {result.pump_power_fuel.value / 1000:.1f} kW", + f" Turbine Flow: {result.turbine_mass_flow.value:.3f} kg/s", + "", + "TANK REQUIREMENTS:", + f" Ox Tank Pressure: {result.tank_pressure_ox.to('bar').value:.1f} bar", + f" Fuel Tank Pressure:{result.tank_pressure_fuel.to('bar').value:.1f} bar", + "", + "NPSH MARGINS:", + f" Ox NPSH Margin: {result.npsh_margin_ox.to('bar').value:.2f} bar", + f" Fuel NPSH Margin: {result.npsh_margin_fuel.to('bar').value:.2f} bar", + ] + + if result.warnings: + lines.extend(["", "WARNINGS:"]) + for warning in result.warnings: + lines.append(f" ⚠ {warning}") + + lines.append(f"{'=' * 60}") + + return "\n".join(lines) + diff --git a/rocket/cycles/gas_generator.py b/rocket/cycles/gas_generator.py new file mode 100644 index 0000000..723a1ab --- /dev/null +++ b/rocket/cycles/gas_generator.py @@ -0,0 +1,340 @@ +"""Gas generator engine cycle analysis. + +The gas generator (GG) cycle is the most common turbopump-fed cycle. +A small portion of propellants is burned in a separate gas generator +to drive the turbine, then exhausted overboard (or through a secondary nozzle). + +Advantages: +- Proven, reliable technology +- Simpler than staged combustion +- Lower turbine temperatures +- Decoupled turbine from main chamber + +Disadvantages: +- GG exhaust is "wasted" (reduces effective Isp by 1-3%) +- Limited chamber pressure compared to staged combustion +- Requires separate GG and associated plumbing + +Examples: +- SpaceX Merlin (LOX/RP-1) +- Rocketdyne F-1 (LOX/RP-1) +- RS-68 (LOX/LH2) +- Vulcain (LOX/LH2) +""" + +from dataclasses import dataclass + +from beartype import beartype + +from rocket.cycles.base import ( + CyclePerformance, + CycleType, + npsh_available, + pump_power, +) +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.tanks import get_propellant_density +from rocket.units import Quantity, kelvin, kg_per_second, pascals, seconds + +# Typical vapor pressures for common propellants [Pa] +VAPOR_PRESSURES: dict[str, float] = { + "LOX": 101325, + "LH2": 101325, + "CH4": 101325, + "RP1": 1000, + "Ethanol": 5900, +} + + +def _get_vapor_pressure(propellant: str) -> float: + """Get vapor pressure for a propellant.""" + name = propellant.upper().replace("-", "").replace(" ", "") + for key, value in VAPOR_PRESSURES.items(): + if key.upper() == name: + return value + return 1000.0 + + +@beartype +@dataclass(frozen=True, slots=True) +class GasGeneratorCycle: + """Configuration for a gas generator engine cycle. + + The gas generator produces hot gas to drive the turbines that power + the propellant pumps. The GG exhaust is typically dumped overboard. + + Attributes: + turbine_inlet_temp: Gas generator combustion temperature [K] + Typically 700-1000K to protect turbine blades + pump_efficiency_ox: Oxidizer pump efficiency (0.6-0.75 typical) + pump_efficiency_fuel: Fuel pump efficiency (0.6-0.75 typical) + turbine_efficiency: Turbine isentropic efficiency (0.5-0.7 typical) + turbine_pressure_ratio: Turbine inlet/outlet pressure ratio (2-6 typical) + gg_mixture_ratio: GG O/F ratio (fuel-rich, typically 0.3-0.5) + mechanical_efficiency: Mechanical losses in turbopump (0.95-0.98) + tank_pressure_ox: Oxidizer tank pressure [Pa] + tank_pressure_fuel: Fuel tank pressure [Pa] + """ + + turbine_inlet_temp: Quantity = None # type: ignore # Will validate in __post_init__ + pump_efficiency_ox: float = 0.70 + pump_efficiency_fuel: float = 0.70 + turbine_efficiency: float = 0.60 + turbine_pressure_ratio: float = 4.0 + gg_mixture_ratio: float = 0.4 # Fuel-rich + mechanical_efficiency: float = 0.97 + tank_pressure_ox: Quantity | None = None + tank_pressure_fuel: Quantity | None = None + + def __post_init__(self) -> None: + """Validate inputs.""" + if self.turbine_inlet_temp is None: + object.__setattr__(self, 'turbine_inlet_temp', kelvin(900)) + + @property + def cycle_type(self) -> CycleType: + return CycleType.GAS_GENERATOR + + def analyze( + self, + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + ) -> CyclePerformance: + """Analyze gas generator cycle and compute net performance. + + The key equations are: + 1. Pump power = mdot * delta_P / (rho * eta) + 2. Turbine power = mdot_gg * cp * delta_T * eta + 3. Power balance: P_turbine = P_pump_ox + P_pump_fuel + 4. Net Isp = (F_main - mdot_gg * ue_gg) / (mdot_total * g0) + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + + Returns: + CyclePerformance with net performance and power balance + """ + warnings: list[str] = [] + + # Extract values + pc = inputs.chamber_pressure.to("Pa").value + mdot_total = performance.mdot.to("kg/s").value + mdot_ox = performance.mdot_ox.to("kg/s").value + mdot_fuel = performance.mdot_fuel.to("kg/s").value + isp = performance.isp.value + thrust = inputs.thrust.to("N").value + + # Get propellant properties + # Try to determine from engine name + ox_name = "LOX" + fuel_name = "RP1" + if inputs.name: + name_upper = inputs.name.upper() + if "CH4" in name_upper or "METHANE" in name_upper or "METHALOX" in name_upper: + fuel_name = "CH4" + elif "LH2" in name_upper or "HYDROGEN" in name_upper or "HYDROLOX" in name_upper: + fuel_name = "LH2" + + try: + rho_ox = get_propellant_density(ox_name) + except ValueError: + rho_ox = 1141.0 + warnings.append(f"Unknown oxidizer, assuming LOX density: {rho_ox} kg/m³") + + try: + rho_fuel = get_propellant_density(fuel_name) + except ValueError: + rho_fuel = 810.0 + warnings.append(f"Unknown fuel, assuming RP-1 density: {rho_fuel} kg/m³") + + # Determine tank pressures + if self.tank_pressure_ox is not None: + p_tank_ox = self.tank_pressure_ox.to("Pa").value + else: + # Typical tank pressure for turbopump-fed: 2-5 bar + p_tank_ox = 300000 # 3 bar default + + if self.tank_pressure_fuel is not None: + p_tank_fuel = self.tank_pressure_fuel.to("Pa").value + else: + p_tank_fuel = 250000 # 2.5 bar default + + # Calculate pump pressure rise + # Pump must raise from tank pressure to chamber pressure + injector drop + margins + injector_dp = pc * 0.20 # 20% pressure drop across injector + p_pump_outlet = pc + injector_dp + + dp_ox = p_pump_outlet - p_tank_ox + dp_fuel = p_pump_outlet - p_tank_fuel + + # Pump power requirements + P_pump_ox = pump_power( + mass_flow=kg_per_second(mdot_ox), + pressure_rise=pascals(dp_ox), + density=rho_ox, + efficiency=self.pump_efficiency_ox, + ).value + + P_pump_fuel = pump_power( + mass_flow=kg_per_second(mdot_fuel), + pressure_rise=pascals(dp_fuel), + density=rho_fuel, + efficiency=self.pump_efficiency_fuel, + ).value + + # Total pump power (accounting for mechanical losses) + P_pump_total = (P_pump_ox + P_pump_fuel) / self.mechanical_efficiency + + # Gas generator analysis + # Turbine power must equal pump power + P_turbine_required = P_pump_total + + # GG exhaust properties (fuel-rich combustion products) + # Approximate gamma and R for fuel-rich GG + gamma_gg = 1.25 # Lower gamma for fuel-rich + R_gg = 350.0 # J/(kg·K), approximate for fuel-rich products + T_gg = self.turbine_inlet_temp.to("K").value + + # Turbine specific work + # w = cp * T_in * eta * (1 - 1/PR^((gamma-1)/gamma)) + cp_gg = gamma_gg * R_gg / (gamma_gg - 1) + T_ratio = self.turbine_pressure_ratio ** ((gamma_gg - 1) / gamma_gg) + w_turbine = cp_gg * T_gg * self.turbine_efficiency * (1 - 1/T_ratio) + + # GG mass flow required + mdot_gg = P_turbine_required / w_turbine + + # Check GG flow is reasonable (typically 1-5% of total) + gg_fraction = mdot_gg / mdot_total + if gg_fraction > 0.10: + warnings.append( + f"GG flow is {gg_fraction*100:.1f}% of total - unusually high" + ) + + # GG propellant split + # mdot_gg_ox = mdot_gg * self.gg_mixture_ratio / (1 + self.gg_mixture_ratio) + # mdot_gg_fuel = mdot_gg / (1 + self.gg_mixture_ratio) + + # Net performance calculation + # The GG exhaust has much lower velocity than main chamber + # Approximate GG exhaust velocity + # For low MR (fuel-rich): Isp_gg ~ 200-250s + isp_gg = 220.0 # s, approximate for fuel-rich GG exhaust + g0 = 9.80665 + ue_gg = isp_gg * g0 + + # GG exhaust thrust (negative contribution to net thrust) + F_gg = mdot_gg * ue_gg + + # Net thrust and Isp + # Main chamber produces full thrust + # But we've "spent" mdot_gg propellant for low-Isp exhaust + F_main = thrust + net_thrust = F_main # GG exhaust typically dumps to atmosphere + + # Effective total mass flow (main + GG) + mdot_effective = mdot_total # GG flow comes from same tanks + + # Net Isp considering the GG "loss" + # Two ways to think about it: + # 1. All propellant flows through main chamber at full Isp + # 2. GG flow produces low-Isp exhaust + # Net: weighted average of Isp + net_isp = (F_main + F_gg * 0.3) / (mdot_effective * g0) # GG contributes ~30% of its thrust + + # Alternative: simple debit approach + # net_isp = isp * (1 - gg_fraction) + isp_gg * gg_fraction + net_isp_alt = isp * (1 - gg_fraction * 0.7) # ~70% loss on GG flow + + # Use the more conservative estimate + net_isp = min(net_isp, net_isp_alt) + + cycle_efficiency = net_isp / isp + + # NPSH analysis + p_vapor_ox = _get_vapor_pressure(ox_name) + p_vapor_fuel = _get_vapor_pressure(fuel_name) + + npsh_ox = npsh_available( + tank_pressure=pascals(p_tank_ox), + fluid_density=rho_ox, + vapor_pressure=pascals(p_vapor_ox), + ) + + npsh_fuel = npsh_available( + tank_pressure=pascals(p_tank_fuel), + fluid_density=rho_fuel, + vapor_pressure=pascals(p_vapor_fuel), + ) + + # Feasibility checks + feasible = True + + if npsh_ox.to("Pa").value < 50000: # < 0.5 bar + warnings.append("Low NPSH margin for oxidizer pump - risk of cavitation") + + if npsh_fuel.to("Pa").value < 50000: + warnings.append("Low NPSH margin for fuel pump - risk of cavitation") + + if T_gg > 1100: + warnings.append( + f"Turbine inlet temp {T_gg:.0f}K exceeds typical limit (~1000K)" + ) + + if gg_fraction > 0.05: + warnings.append( + f"GG fraction {gg_fraction*100:.1f}% is high - consider staged combustion" + ) + + return CyclePerformance( + net_isp=seconds(net_isp), + net_thrust=Quantity(net_thrust, "N", "force"), + cycle_efficiency=cycle_efficiency, + pump_power_ox=Quantity(P_pump_ox, "W", "power"), + pump_power_fuel=Quantity(P_pump_fuel, "W", "power"), + turbine_power=Quantity(P_turbine_required, "W", "power"), + turbine_mass_flow=kg_per_second(mdot_gg), + tank_pressure_ox=pascals(p_tank_ox), + tank_pressure_fuel=pascals(p_tank_fuel), + npsh_margin_ox=npsh_ox, + npsh_margin_fuel=npsh_fuel, + cycle_type=self.cycle_type, + feasible=feasible, + warnings=warnings, + ) + + +@beartype +def estimate_turbopump_mass( + pump_power: Quantity, + turbine_power: Quantity, + propellant_type: str = "LOX/RP1", +) -> Quantity: + """Estimate turbopump mass from power requirements. + + Uses historical correlations from existing engines. + + Args: + pump_power: Total pump power [W] + turbine_power: Turbine power [W] + propellant_type: Propellant combination for correlation selection + + Returns: + Estimated turbopump mass [kg] + """ + P = max(pump_power.value, turbine_power.value) + + # Historical correlation: mass ~ k * P^0.6 + # k varies by propellant type and technology level + k = 0.015 if "LH2" in propellant_type.upper() else 0.008 # LH2 pumps are larger due to low density + + mass = k * P ** 0.6 + + # Minimum mass for small turbopumps + mass = max(mass, 5.0) + + return Quantity(mass, "kg", "mass") + diff --git a/rocket/cycles/pressure_fed.py b/rocket/cycles/pressure_fed.py new file mode 100644 index 0000000..6d2440c --- /dev/null +++ b/rocket/cycles/pressure_fed.py @@ -0,0 +1,287 @@ +"""Pressure-fed engine cycle analysis. + +Pressure-fed engines use high-pressure gas (typically helium) to push +propellants from tanks into the combustion chamber. They are the simplest +cycle type but require heavy tanks to contain the high pressures. + +Advantages: +- Simplicity and reliability (no turbopumps) +- Fewer failure modes +- Lower development cost + +Disadvantages: +- Heavy tanks (must withstand full chamber pressure + margins) +- Limited chamber pressure (~3 MPa practical limit) +- Lower performance (limited Isp due to pressure constraints) + +Typical applications: +- Upper stages +- Spacecraft thrusters +- Student/amateur rockets +""" + +from dataclasses import dataclass + +from beartype import beartype + +from rocket.cycles.base import ( + CyclePerformance, + CycleType, + estimate_line_losses, + npsh_available, +) +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.tanks import get_propellant_density +from rocket.units import Quantity, kg_per_second, pascals + +# Typical vapor pressures for common propellants [Pa] +# At nominal storage temperatures +VAPOR_PRESSURES: dict[str, float] = { + "LOX": 101325, # ~1 atm at -183°C + "LH2": 101325, # ~1 atm at -253°C + "CH4": 101325, # ~1 atm at -161°C + "RP1": 1000, # Very low at 20°C + "Ethanol": 5900, # ~0.06 atm at 20°C + "N2O4": 96000, # ~0.95 atm at 20°C + "MMH": 4800, # Low at 20°C + "N2O": 5200000, # ~51 atm at 20°C (self-pressurizing) +} + + +def _get_vapor_pressure(propellant: str) -> float: + """Get vapor pressure for a propellant.""" + # Normalize name + name = propellant.upper().replace("-", "").replace(" ", "") + + for key, value in VAPOR_PRESSURES.items(): + if key.upper() == name: + return value + + # Default to low vapor pressure if unknown + return 1000.0 + + +@beartype +@dataclass(frozen=True, slots=True) +class PressureFedCycle: + """Configuration for a pressure-fed engine cycle. + + In a pressure-fed system, the tank pressure must exceed the chamber + pressure plus all line losses and injector pressure drop. + + Attributes: + injector_dp_fraction: Injector pressure drop as fraction of Pc (typically 0.15-0.25) + line_loss_fraction: Feed line losses as fraction of Pc (typically 0.05-0.10) + tank_pressure_margin: Safety margin on tank pressure (typically 1.1-1.2) + pressurant: Pressurant gas type (typically "helium" or "nitrogen") + ox_line_diameter: Oxidizer feed line diameter [m] + fuel_line_diameter: Fuel feed line diameter [m] + ox_line_length: Oxidizer feed line length [m] + fuel_line_length: Fuel feed line length [m] + """ + + injector_dp_fraction: float = 0.20 + line_loss_fraction: float = 0.05 + tank_pressure_margin: float = 1.15 + pressurant: str = "helium" + ox_line_diameter: float = 0.05 # m + fuel_line_diameter: float = 0.04 # m + ox_line_length: float = 2.0 # m + fuel_line_length: float = 2.0 # m + + @property + def cycle_type(self) -> CycleType: + return CycleType.PRESSURE_FED + + def analyze( + self, + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + ) -> CyclePerformance: + """Analyze pressure-fed cycle and determine tank pressures. + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + + Returns: + CyclePerformance with pressure requirements and feasibility + """ + warnings: list[str] = [] + + # Extract values + pc = inputs.chamber_pressure.to("Pa").value + mdot_ox = performance.mdot_ox.to("kg/s").value + mdot_fuel = performance.mdot_fuel.to("kg/s").value + + # Get propellant densities + # Extract propellant names from engine name or use defaults + ox_name = "LOX" # Default + fuel_name = "RP1" # Default + if inputs.name: + name_upper = inputs.name.upper() + if "LOX" in name_upper or "LO2" in name_upper: + ox_name = "LOX" + if "CH4" in name_upper or "METHANE" in name_upper: + fuel_name = "CH4" + elif "RP1" in name_upper or "KEROSENE" in name_upper: + fuel_name = "RP1" + elif "ETHANOL" in name_upper: + fuel_name = "Ethanol" + elif "LH2" in name_upper or "HYDROGEN" in name_upper: + fuel_name = "LH2" + + try: + rho_ox = get_propellant_density(ox_name) + except ValueError: + rho_ox = 1141.0 # Default LOX + warnings.append(f"Unknown oxidizer, assuming LOX density: {rho_ox} kg/m³") + + try: + rho_fuel = get_propellant_density(fuel_name) + except ValueError: + rho_fuel = 810.0 # Default RP-1 + warnings.append(f"Unknown fuel, assuming RP-1 density: {rho_fuel} kg/m³") + + # Calculate pressure budget + # Tank pressure must overcome: chamber + injector drop + line losses + dp_injector = pc * self.injector_dp_fraction + dp_lines_ox = pc * self.line_loss_fraction + dp_lines_fuel = pc * self.line_loss_fraction + + # More detailed line loss estimate + dp_lines_ox_calc = estimate_line_losses( + mass_flow=kg_per_second(mdot_ox), + density=rho_ox, + pipe_diameter=self.ox_line_diameter, + pipe_length=self.ox_line_length, + ).to("Pa").value + + dp_lines_fuel_calc = estimate_line_losses( + mass_flow=kg_per_second(mdot_fuel), + density=rho_fuel, + pipe_diameter=self.fuel_line_diameter, + pipe_length=self.fuel_line_length, + ).to("Pa").value + + # Use maximum of estimated and calculated + dp_lines_ox = max(dp_lines_ox, dp_lines_ox_calc) + dp_lines_fuel = max(dp_lines_fuel, dp_lines_fuel_calc) + + # Required tank pressures + p_tank_ox = (pc + dp_injector + dp_lines_ox) * self.tank_pressure_margin + p_tank_fuel = (pc + dp_injector + dp_lines_fuel) * self.tank_pressure_margin + + # NPSH analysis + p_vapor_ox = _get_vapor_pressure(ox_name) + p_vapor_fuel = _get_vapor_pressure(fuel_name) + + npsh_ox = npsh_available( + tank_pressure=pascals(p_tank_ox), + fluid_density=rho_ox, + vapor_pressure=pascals(p_vapor_ox), + line_losses=pascals(dp_lines_ox), + ) + + npsh_fuel = npsh_available( + tank_pressure=pascals(p_tank_fuel), + fluid_density=rho_fuel, + vapor_pressure=pascals(p_vapor_fuel), + line_losses=pascals(dp_lines_fuel), + ) + + # Check feasibility + feasible = True + + # Pressure-fed practical limit is ~3-4 MPa + if pc > 4e6: + warnings.append( + f"Chamber pressure {pc/1e6:.1f} MPa exceeds typical pressure-fed limit (~3-4 MPa)" + ) + + # Tank pressure feasibility + if p_tank_ox > 6e6: + warnings.append( + f"Ox tank pressure {p_tank_ox/1e6:.1f} MPa is very high for pressure-fed" + ) + if p_tank_fuel > 6e6: + warnings.append( + f"Fuel tank pressure {p_tank_fuel/1e6:.1f} MPa is very high for pressure-fed" + ) + + # For pressure-fed, there are no turbopumps, so no pump power + # All "pumping" is done by the pressurized tanks + pump_power_ox = Quantity(0.0, "W", "power") + pump_power_fuel = Quantity(0.0, "W", "power") + turbine_power = Quantity(0.0, "W", "power") + turbine_flow = kg_per_second(0.0) + + # Net performance equals ideal performance (no turbine drive losses) + net_isp = performance.isp + net_thrust = inputs.thrust + cycle_efficiency = 1.0 # No cycle losses for pressure-fed + + return CyclePerformance( + net_isp=net_isp, + net_thrust=net_thrust, + cycle_efficiency=cycle_efficiency, + pump_power_ox=pump_power_ox, + pump_power_fuel=pump_power_fuel, + turbine_power=turbine_power, + turbine_mass_flow=turbine_flow, + tank_pressure_ox=pascals(p_tank_ox), + tank_pressure_fuel=pascals(p_tank_fuel), + npsh_margin_ox=npsh_ox, + npsh_margin_fuel=npsh_fuel, + cycle_type=self.cycle_type, + feasible=feasible, + warnings=warnings, + ) + + +@beartype +def estimate_pressurant_mass( + propellant_volume: Quantity, + tank_pressure: Quantity, + pressurant: str = "helium", + initial_temp: float = 300.0, # K + blowdown_ratio: float = 2.0, +) -> Quantity: + """Estimate pressurant gas mass required. + + Uses ideal gas law with blowdown consideration. + In blowdown mode, tank pressure drops as propellant is expelled. + + Args: + propellant_volume: Volume of propellant to expel [m³] + tank_pressure: Initial tank pressure [Pa] + pressurant: Gas type ("helium" or "nitrogen") + initial_temp: Pressurant initial temperature [K] + blowdown_ratio: Initial/final pressure ratio for blowdown + + Returns: + Required pressurant mass [kg] + """ + # Gas constants + R_helium = 2077.0 # J/(kg·K) + R_nitrogen = 296.8 # J/(kg·K) + + R = R_helium if pressurant.lower() == "helium" else R_nitrogen + + V = propellant_volume.to("m^3").value + P = tank_pressure.to("Pa").value + + # For pressure-regulated system: m = P * V / (R * T) + # For blowdown: need to account for pressure decay + # Simplified: assume average pressure + P_avg = P / (1 + 1/blowdown_ratio) * 2 + + mass = P_avg * V / (R * initial_temp) + + # Add margin for residuals and cooling + mass *= 1.2 + + return Quantity(mass, "kg", "mass") + diff --git a/rocket/cycles/staged_combustion.py b/rocket/cycles/staged_combustion.py new file mode 100644 index 0000000..743884b --- /dev/null +++ b/rocket/cycles/staged_combustion.py @@ -0,0 +1,485 @@ +"""Staged combustion engine cycle analysis. + +Staged combustion is the highest-performance liquid engine cycle. Unlike +gas generators where turbine exhaust is dumped overboard, staged combustion +routes all turbine exhaust into the main combustion chamber. + +Variants: +- Oxidizer-rich staged combustion (ORSC): Preburner runs oxidizer-rich + Example: RD-180, RD-191, NK-33 +- Fuel-rich staged combustion (FRSC): Preburner runs fuel-rich + Example: RS-25 (SSME), BE-4 +- Full-flow staged combustion (FFSC): Both ox-rich AND fuel-rich preburners + Example: SpaceX Raptor + +Advantages: +- Highest Isp (all propellant goes through main chamber) +- High chamber pressure capability +- High thrust-to-weight ratio + +Disadvantages: +- Most complex cycle +- Expensive development +- Challenging turbine environments (especially ORSC) + +References: + - Sutton & Biblarz, Chapter 6 + - Humble, Henry & Larson, "Space Propulsion Analysis and Design" +""" + +from dataclasses import dataclass + +from beartype import beartype + +from rocket.cycles.base import ( + CyclePerformance, + CycleType, + npsh_available, + pump_power, +) +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.tanks import get_propellant_density +from rocket.units import Quantity, kg_per_second, pascals, seconds + +# Typical vapor pressures for common propellants [Pa] +VAPOR_PRESSURES: dict[str, float] = { + "LOX": 101325, + "LH2": 101325, + "CH4": 101325, + "RP1": 1000, +} + + +def _get_vapor_pressure(propellant: str) -> float: + """Get vapor pressure for a propellant.""" + name = propellant.upper().replace("-", "").replace(" ", "") + for key, value in VAPOR_PRESSURES.items(): + if key.upper() == name: + return value + return 1000.0 + + +@beartype +@dataclass(frozen=True, slots=True) +class StagedCombustionCycle: + """Configuration for a staged combustion engine cycle. + + In staged combustion, the preburner exhaust (which drove the turbine) + is routed into the main combustion chamber, so no propellant is wasted. + + Attributes: + preburner_temp: Preburner combustion temperature [K] + Typically 700-900K for fuel-rich, 500-700K for ox-rich + pump_efficiency_ox: Oxidizer pump efficiency (0.7-0.8 typical) + pump_efficiency_fuel: Fuel pump efficiency (0.7-0.8 typical) + turbine_efficiency: Turbine isentropic efficiency (0.6-0.75 typical) + turbine_pressure_ratio: Turbine pressure ratio (1.5-3.0 typical) + preburner_mixture_ratio: Preburner O/F ratio + Fuel-rich: 0.3-0.6 (for SSME-type) + Ox-rich: 50-100 (for RD-180-type) + oxidizer_rich: If True, uses ox-rich preburner (ORSC) + mechanical_efficiency: Mechanical losses (0.95-0.98) + tank_pressure_ox: Oxidizer tank pressure [Pa] + tank_pressure_fuel: Fuel tank pressure [Pa] + """ + + preburner_temp: Quantity | None = None + pump_efficiency_ox: float = 0.75 + pump_efficiency_fuel: float = 0.75 + turbine_efficiency: float = 0.70 + turbine_pressure_ratio: float = 2.0 + preburner_mixture_ratio: float = 0.5 # Fuel-rich default + oxidizer_rich: bool = False + mechanical_efficiency: float = 0.97 + tank_pressure_ox: Quantity | None = None + tank_pressure_fuel: Quantity | None = None + + @property + def cycle_type(self) -> CycleType: + return CycleType.STAGED_COMBUSTION + + def analyze( + self, + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + ) -> CyclePerformance: + """Analyze staged combustion cycle. + + The key difference from gas generator is that turbine exhaust + goes to the main chamber, so there's no Isp penalty from + dumping low-energy gases. + + The power balance is more complex because the preburner + operates at a pressure higher than the main chamber. + """ + warnings: list[str] = [] + + # Extract values + pc = inputs.chamber_pressure.to("Pa").value + mdot_ox = performance.mdot_ox.to("kg/s").value + mdot_fuel = performance.mdot_fuel.to("kg/s").value + isp = performance.isp.value + + # Set default preburner temperature + if self.preburner_temp is not None: + T_pb = self.preburner_temp.to("K").value + else: + # Default based on cycle type + T_pb = 600.0 if self.oxidizer_rich else 800.0 + + # Get propellant properties + ox_name = "LOX" + fuel_name = "RP1" + if inputs.name: + name_upper = inputs.name.upper() + if "CH4" in name_upper or "METHANE" in name_upper: + fuel_name = "CH4" + elif "LH2" in name_upper or "HYDROGEN" in name_upper: + fuel_name = "LH2" + + try: + rho_ox = get_propellant_density(ox_name) + except ValueError: + rho_ox = 1141.0 + warnings.append(f"Unknown oxidizer, assuming LOX density: {rho_ox} kg/m³") + + try: + rho_fuel = get_propellant_density(fuel_name) + except ValueError: + rho_fuel = 810.0 + warnings.append(f"Unknown fuel, assuming RP-1 density: {rho_fuel} kg/m³") + + # Tank pressures + if self.tank_pressure_ox is not None: + p_tank_ox = self.tank_pressure_ox.to("Pa").value + else: + p_tank_ox = 400000 # 4 bar typical for staged combustion + + if self.tank_pressure_fuel is not None: + p_tank_fuel = self.tank_pressure_fuel.to("Pa").value + else: + p_tank_fuel = 350000 # 3.5 bar + + # Preburner pressure must be higher than chamber pressure + # Turbine pressure drop + injector losses + p_preburner = pc * 1.3 # 30% higher than chamber + + # Pump pressure rises + # Pumps must deliver to preburner pressure (higher than PC) + injector_dp = pc * 0.20 + p_pump_outlet = p_preburner + injector_dp + + dp_ox = p_pump_outlet - p_tank_ox + dp_fuel = p_pump_outlet - p_tank_fuel + + # Pump power requirements + P_pump_ox = pump_power( + mass_flow=kg_per_second(mdot_ox), + pressure_rise=pascals(dp_ox), + density=rho_ox, + efficiency=self.pump_efficiency_ox, + ).value + + P_pump_fuel = pump_power( + mass_flow=kg_per_second(mdot_fuel), + pressure_rise=pascals(dp_fuel), + density=rho_fuel, + efficiency=self.pump_efficiency_fuel, + ).value + + # Total pump power + P_pump_total = (P_pump_ox + P_pump_fuel) / self.mechanical_efficiency + + # Preburner/turbine analysis + # In staged combustion, ALL propellant eventually goes through main chamber + # Preburner flow drives turbine, then goes to main chamber + + if self.oxidizer_rich: + # Ox-rich: most of oxidizer through preburner with small fuel + # mdot_pb = mdot_ox + mdot_pb_fuel + # Preburner MR is very high (50-100), so mdot_pb_fuel is small + mdot_pb_fuel = mdot_ox / self.preburner_mixture_ratio + mdot_pb = mdot_ox + mdot_pb_fuel + # Remaining fuel goes directly to main chamber + mdot_direct_fuel = mdot_fuel - mdot_pb_fuel + + if mdot_direct_fuel < 0: + warnings.append("Preburner consumes more fuel than available - infeasible") + mdot_direct_fuel = 0 + + # Preburner exhaust is mostly oxygen with some combustion products + gamma_pb = 1.30 # Higher gamma for ox-rich + R_pb = 280.0 # J/(kg·K) + + else: + # Fuel-rich: most of fuel through preburner with small oxidizer + mdot_pb_ox = mdot_fuel * self.preburner_mixture_ratio + mdot_pb = mdot_fuel + mdot_pb_ox + # Remaining oxidizer goes directly to main chamber + mdot_direct_ox = mdot_ox - mdot_pb_ox + + if mdot_direct_ox < 0: + warnings.append("Preburner consumes more oxidizer than available - infeasible") + mdot_direct_ox = 0 + + # Preburner exhaust is fuel-rich combustion products + gamma_pb = 1.20 # Lower gamma for fuel-rich + R_pb = 400.0 # J/(kg·K), higher for lighter products + + # Turbine power available + cp_pb = gamma_pb * R_pb / (gamma_pb - 1) + T_ratio = self.turbine_pressure_ratio ** ((gamma_pb - 1) / gamma_pb) + w_turbine = cp_pb * T_pb * self.turbine_efficiency * (1 - 1/T_ratio) + + P_turbine_available = mdot_pb * w_turbine + + # Power balance check + power_margin = P_turbine_available / P_pump_total if P_pump_total > 0 else float('inf') + + if power_margin < 1.0: + warnings.append( + f"Power balance not achieved: turbine provides {P_turbine_available/1e6:.1f} MW, " + f"pumps need {P_pump_total/1e6:.1f} MW" + ) + + # Net performance + # In staged combustion, ALL propellant goes through main chamber + # at (nearly) full Isp, so cycle efficiency is very high + # Small losses from: + # 1. Preburner inefficiency + # 2. Slightly different combustion from preburned products + + # Estimate efficiency loss (typically 1-3%) + efficiency_loss = 0.02 # 2% loss typical for staged combustion + net_isp = isp * (1 - efficiency_loss) + cycle_efficiency = 1 - efficiency_loss + + # NPSH analysis + p_vapor_ox = _get_vapor_pressure(ox_name) + p_vapor_fuel = _get_vapor_pressure(fuel_name) + + npsh_ox = npsh_available( + tank_pressure=pascals(p_tank_ox), + fluid_density=rho_ox, + vapor_pressure=pascals(p_vapor_ox), + ) + + npsh_fuel = npsh_available( + tank_pressure=pascals(p_tank_fuel), + fluid_density=rho_fuel, + vapor_pressure=pascals(p_vapor_fuel), + ) + + # Feasibility assessment + feasible = power_margin >= 0.95 # Allow small margin + + if power_margin < 1.0: + feasible = False + + if self.oxidizer_rich and T_pb > 700: + warnings.append( + f"Ox-rich preburner at {T_pb:.0f}K - requires specialized turbine materials" + ) + + if not self.oxidizer_rich and T_pb > 1000: + warnings.append( + f"Fuel-rich preburner temp {T_pb:.0f}K is high" + ) + + if pc > 25e6: + warnings.append( + f"Chamber pressure {pc/1e6:.0f} MPa is very high - " + "typical staged combustion limit ~30 MPa" + ) + + return CyclePerformance( + net_isp=seconds(net_isp), + net_thrust=inputs.thrust, # Staged combustion delivers full thrust + cycle_efficiency=cycle_efficiency, + pump_power_ox=Quantity(P_pump_ox, "W", "power"), + pump_power_fuel=Quantity(P_pump_fuel, "W", "power"), + turbine_power=Quantity(P_turbine_available, "W", "power"), + turbine_mass_flow=kg_per_second(mdot_pb), + tank_pressure_ox=pascals(p_tank_ox), + tank_pressure_fuel=pascals(p_tank_fuel), + npsh_margin_ox=npsh_ox, + npsh_margin_fuel=npsh_fuel, + cycle_type=self.cycle_type, + feasible=feasible, + warnings=warnings, + ) + + +@beartype +@dataclass(frozen=True, slots=True) +class FullFlowStagedCombustion: + """Configuration for full-flow staged combustion (FFSC). + + FFSC uses TWO preburners: + - Fuel-rich preburner: drives fuel turbopump + - Ox-rich preburner: drives oxidizer turbopump + + This provides the highest possible performance and allows + independent control of each turbopump. + + Example: SpaceX Raptor + + Attributes: + fuel_preburner_temp: Fuel-rich preburner temperature [K] + ox_preburner_temp: Ox-rich preburner temperature [K] + pump_efficiency: Pump efficiency for both pumps + turbine_efficiency: Turbine efficiency for both turbines + fuel_turbine_pr: Fuel turbine pressure ratio + ox_turbine_pr: Ox turbine pressure ratio + """ + + fuel_preburner_temp: Quantity | None = None + ox_preburner_temp: Quantity | None = None + pump_efficiency: float = 0.77 + turbine_efficiency: float = 0.72 + fuel_turbine_pr: float = 1.8 + ox_turbine_pr: float = 1.5 + tank_pressure_ox: Quantity | None = None + tank_pressure_fuel: Quantity | None = None + + @property + def cycle_type(self) -> CycleType: + return CycleType.FULL_FLOW_STAGED + + def analyze( + self, + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + ) -> CyclePerformance: + """Analyze full-flow staged combustion cycle.""" + warnings: list[str] = [] + + # Extract values + pc = inputs.chamber_pressure.to("Pa").value + mdot_ox = performance.mdot_ox.to("kg/s").value + mdot_fuel = performance.mdot_fuel.to("kg/s").value + isp = performance.isp.value + + # Default preburner temperatures + T_fuel_pb = self.fuel_preburner_temp.to("K").value if self.fuel_preburner_temp else 800.0 + T_ox_pb = self.ox_preburner_temp.to("K").value if self.ox_preburner_temp else 600.0 + + # Get propellant properties + try: + rho_ox = get_propellant_density("LOX") + except ValueError: + rho_ox = 1141.0 + + try: + rho_fuel = get_propellant_density("CH4") # FFSC typically uses methane + except ValueError: + rho_fuel = 422.0 + + # Tank pressures + p_tank_ox = self.tank_pressure_ox.to("Pa").value if self.tank_pressure_ox else 500000 + p_tank_fuel = self.tank_pressure_fuel.to("Pa").value if self.tank_pressure_fuel else 450000 + + # Preburner pressures (higher than chamber) + p_preburner = pc * 1.25 + + # Pump requirements + dp_ox = p_preburner - p_tank_ox + pc * 0.2 + dp_fuel = p_preburner - p_tank_fuel + pc * 0.2 + + P_pump_ox = pump_power( + mass_flow=kg_per_second(mdot_ox), + pressure_rise=pascals(dp_ox), + density=rho_ox, + efficiency=self.pump_efficiency, + ).value + + P_pump_fuel = pump_power( + mass_flow=kg_per_second(mdot_fuel), + pressure_rise=pascals(dp_fuel), + density=rho_fuel, + efficiency=self.pump_efficiency, + ).value + + # In FFSC, all oxidizer goes through ox-rich preburner + # All fuel goes through fuel-rich preburner + # Each preburner drives its respective turbopump + + # Fuel-rich preburner (drives fuel pump) + # Small amount of ox mixed with all fuel + gamma_fuel_pb = 1.18 + R_fuel_pb = 450.0 + cp_fuel_pb = gamma_fuel_pb * R_fuel_pb / (gamma_fuel_pb - 1) + T_ratio_fuel = self.fuel_turbine_pr ** ((gamma_fuel_pb - 1) / gamma_fuel_pb) + w_fuel_turbine = cp_fuel_pb * T_fuel_pb * self.turbine_efficiency * (1 - 1/T_ratio_fuel) + + # Ox-rich preburner (drives ox pump) + gamma_ox_pb = 1.30 + R_ox_pb = 280.0 + cp_ox_pb = gamma_ox_pb * R_ox_pb / (gamma_ox_pb - 1) + T_ratio_ox = self.ox_turbine_pr ** ((gamma_ox_pb - 1) / gamma_ox_pb) + w_ox_turbine = cp_ox_pb * T_ox_pb * self.turbine_efficiency * (1 - 1/T_ratio_ox) + + # Power available from each turbine + # In FFSC, all propellant flows through preburners + P_fuel_turbine = mdot_fuel * w_fuel_turbine + P_ox_turbine = mdot_ox * w_ox_turbine + + # Check power balance + fuel_margin = P_fuel_turbine / P_pump_fuel if P_pump_fuel > 0 else float('inf') + ox_margin = P_ox_turbine / P_pump_ox if P_pump_ox > 0 else float('inf') + + feasible = True + if fuel_margin < 0.95: + warnings.append(f"Fuel turbopump power margin low: {fuel_margin:.2f}") + feasible = False + if ox_margin < 0.95: + warnings.append(f"Ox turbopump power margin low: {ox_margin:.2f}") + feasible = False + + # FFSC has minimal Isp loss (all propellant to main chamber) + efficiency_loss = 0.01 # ~1% loss + net_isp = isp * (1 - efficiency_loss) + cycle_efficiency = 1 - efficiency_loss + + # NPSH + p_vapor_ox = _get_vapor_pressure("LOX") + p_vapor_fuel = _get_vapor_pressure("CH4") + + npsh_ox = npsh_available( + tank_pressure=pascals(p_tank_ox), + fluid_density=rho_ox, + vapor_pressure=pascals(p_vapor_ox), + ) + + npsh_fuel = npsh_available( + tank_pressure=pascals(p_tank_fuel), + fluid_density=rho_fuel, + vapor_pressure=pascals(p_vapor_fuel), + ) + + # Warnings + if pc > 30e6: + warnings.append(f"Chamber pressure {pc/1e6:.0f} MPa is at FFSC limit") + + if T_ox_pb > 700: + warnings.append("Ox-rich preburner temp requires advanced turbine materials") + + return CyclePerformance( + net_isp=seconds(net_isp), + net_thrust=inputs.thrust, + cycle_efficiency=cycle_efficiency, + pump_power_ox=Quantity(P_pump_ox, "W", "power"), + pump_power_fuel=Quantity(P_pump_fuel, "W", "power"), + turbine_power=Quantity(P_fuel_turbine + P_ox_turbine, "W", "power"), + turbine_mass_flow=kg_per_second(mdot_ox + mdot_fuel), # All flow through preburners + tank_pressure_ox=pascals(p_tank_ox), + tank_pressure_fuel=pascals(p_tank_fuel), + npsh_margin_ox=npsh_ox, + npsh_margin_fuel=npsh_fuel, + cycle_type=self.cycle_type, + feasible=feasible, + warnings=warnings, + ) + diff --git a/rocket/engine.py b/rocket/engine.py index 6d43b43..488bd6a 100644 --- a/rocket/engine.py +++ b/rocket/engine.py @@ -31,10 +31,6 @@ thrust_coefficient, thrust_coefficient_vacuum, ) -from rocket.propellants import ( - get_combustion_properties, - get_optimal_mixture_ratio, -) from rocket.units import ( Quantity, kelvin, @@ -186,6 +182,11 @@ def from_propellants( ... ) >>> print(f"Tc = {inputs.chamber_temp}") """ + from rocket.propellants import ( + get_combustion_properties, + get_optimal_mixture_ratio, + ) + # Default exit pressure to 1 atm if exit_pressure is None: exit_pressure = pascals(101325) diff --git a/rocket/examples/basic_engine.py b/rocket/examples/basic_engine.py index 13f569d..a0515d5 100644 --- a/rocket/examples/basic_engine.py +++ b/rocket/examples/basic_engine.py @@ -11,10 +11,10 @@ 5. Visualize the design 6. Export contour for CAD -The example engine is similar to a small pressure-fed engine suitable -for a student rocket project. +All outputs are organized into a timestamped directory structure. """ +from rocket import OutputContext from rocket.engine import ( EngineInputs, compute_geometry, @@ -161,58 +161,83 @@ def main() -> None: print() # ========================================================================= - # Step 5: Export Contour to CSV + # Step 5-7: Save All Outputs to Organized Directory # ========================================================================= - print("Step 5: Exporting contour to CSV...") - print() - - nozzle_contour.to_csv("nozzle_contour.csv") - full_contour.to_csv("full_chamber_contour.csv") - - print(" Saved: nozzle_contour.csv") - print(" Saved: full_chamber_contour.csv") - print() - - # ========================================================================= - # Step 6: Print Summaries - # ========================================================================= - - print("Step 6: Full summaries...") - print() - print(format_performance_summary(inputs, performance)) - print() - print(format_geometry_summary(inputs, geometry)) - print() - - # ========================================================================= - # Step 7: Create Visualizations - # ========================================================================= - - print("Step 7: Creating visualizations...") - print() - - # Engine cross-section - fig1 = plot_engine_cross_section( - geometry, full_contour, inputs, show_dimensions=True, title=f"{inputs.name} Cross-Section" - ) - fig1.savefig("engine_cross_section.png", dpi=150, bbox_inches="tight") - print(" Saved: engine_cross_section.png") - - # Nozzle contour detail - fig2 = plot_nozzle_contour(nozzle_contour, title=f"{inputs.name} Nozzle Contour") - fig2.savefig("nozzle_contour.png", dpi=150, bbox_inches="tight") - print(" Saved: nozzle_contour.png") - - # Performance vs altitude - fig3 = plot_performance_vs_altitude(inputs, performance, geometry, max_altitude_km=80) - fig3.savefig("altitude_performance.png", dpi=150, bbox_inches="tight") - print(" Saved: altitude_performance.png") - - # Complete dashboard - fig4 = plot_engine_dashboard(inputs, performance, geometry, full_contour) - fig4.savefig("engine_dashboard.png", dpi=150, bbox_inches="tight") - print(" Saved: engine_dashboard.png") + print("Step 5-7: Saving outputs...") + print() + + # Use OutputContext to organize all outputs + with OutputContext("student_engine_mk1", include_timestamp=True) as ctx: + # Add metadata about this run + ctx.add_metadata("engine_name", inputs.name) + ctx.add_metadata("thrust_N", inputs.thrust.to("N").value) + ctx.add_metadata("chamber_pressure_MPa", inputs.chamber_pressure.to("MPa").value) + + # Export contours to CSV (automatically goes to data/) + ctx.log("Exporting nozzle contours...") + nozzle_contour.to_csv(ctx.path("nozzle_contour.csv")) + full_contour.to_csv(ctx.path("full_chamber_contour.csv")) + + # Save text summaries (automatically goes to reports/) + ctx.log("Saving performance summaries...") + ctx.save_text(format_performance_summary(inputs, performance), "performance_summary.txt") + ctx.save_text(format_geometry_summary(inputs, geometry), "geometry_summary.txt") + + # Save design summary as JSON + ctx.save_summary({ + "engine_name": inputs.name, + "performance": { + "isp_sl_s": performance.isp.value, + "isp_vac_s": performance.isp_vac.value, + "cstar_m_s": performance.cstar.value, + "thrust_coeff_sl": performance.thrust_coeff, + "thrust_coeff_vac": performance.thrust_coeff_vac, + "exit_mach": performance.exit_mach, + "mdot_kg_s": performance.mdot.value, + "mdot_ox_kg_s": performance.mdot_ox.value, + "mdot_fuel_kg_s": performance.mdot_fuel.value, + }, + "geometry": { + "throat_diameter_mm": Dt_mm, + "exit_diameter_mm": De_mm, + "chamber_diameter_mm": Dc_mm, + "chamber_length_mm": Lc_mm, + "nozzle_length_mm": Ln_mm, + "expansion_ratio": geometry.expansion_ratio, + "contraction_ratio": geometry.contraction_ratio, + }, + "inputs": { + "thrust_N": inputs.thrust.to("N").value, + "chamber_pressure_Pa": inputs.chamber_pressure.to("Pa").value, + "chamber_temp_K": inputs.chamber_temp.to("K").value, + "molecular_weight": inputs.molecular_weight, + "gamma": inputs.gamma, + "mixture_ratio": inputs.mixture_ratio, + }, + }) + + # Create visualizations (automatically goes to plots/) + ctx.log("Generating visualizations...") + + fig1 = plot_engine_cross_section( + geometry, full_contour, inputs, show_dimensions=True, + title=f"{inputs.name} Cross-Section" + ) + fig1.savefig(ctx.path("engine_cross_section.png"), dpi=150, bbox_inches="tight") + + fig2 = plot_nozzle_contour(nozzle_contour, title=f"{inputs.name} Nozzle Contour") + fig2.savefig(ctx.path("nozzle_contour.png"), dpi=150, bbox_inches="tight") + + fig3 = plot_performance_vs_altitude(inputs, performance, geometry, max_altitude_km=80) + fig3.savefig(ctx.path("altitude_performance.png"), dpi=150, bbox_inches="tight") + + fig4 = plot_engine_dashboard(inputs, performance, geometry, full_contour) + fig4.savefig(ctx.path("engine_dashboard.png"), dpi=150, bbox_inches="tight") + + ctx.log("All outputs saved!") + print() + print(f" Output directory: {ctx.output_dir}") print() print("=" * 70) @@ -222,4 +247,3 @@ def main() -> None: if __name__ == "__main__": main() - diff --git a/rocket/examples/cycle_comparison.py b/rocket/examples/cycle_comparison.py new file mode 100644 index 0000000..e7b74c4 --- /dev/null +++ b/rocket/examples/cycle_comparison.py @@ -0,0 +1,312 @@ +#!/usr/bin/env python +"""Engine cycle comparison example for Rocket. + +This example demonstrates how to compare different engine cycle architectures +for a given set of requirements: + +1. Pressure-fed (simplest, lowest performance) +2. Gas generator (most common, good balance) +3. Staged combustion (highest performance, most complex) + +Understanding cycle tradeoffs is critical for: +- Selecting the right architecture for your mission +- Understanding performance vs. complexity tradeoffs +- Estimating system-level impacts (tank pressure, turbomachinery) +""" + +from rocket import ( + EngineInputs, + OutputContext, + design_engine, + plot_cycle_comparison_bars, + plot_cycle_radar, + plot_cycle_tradeoff, +) +from rocket.cycles import ( + GasGeneratorCycle, + PressureFedCycle, + StagedCombustionCycle, +) +from rocket.units import kelvin, kilonewtons, megapascals + + +def print_header(text: str) -> None: + """Print a formatted section header.""" + print() + print("┌" + "─" * 68 + "┐") + print(f"│ {text:<66} │") + print("└" + "─" * 68 + "┘") + + +def main() -> None: + """Run the engine cycle comparison example.""" + print() + print("╔" + "═" * 68 + "╗") + print("║" + " " * 16 + "ENGINE CYCLE COMPARISON STUDY" + " " * 23 + "║") + print("║" + " " * 18 + "LOX/CH4 Methalox Engine" + " " * 27 + "║") + print("╚" + "═" * 68 + "╝") + + # ========================================================================= + # Common Engine Requirements + # ========================================================================= + + print_header("ENGINE REQUIREMENTS") + + # Base engine thermochemistry (same for all cycles) + base_inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(500), + chamber_pressure=megapascals(15), + mixture_ratio=3.2, + name="Methalox-500", + ) + + print(f" Thrust target: {base_inputs.thrust.to('kN').value:.0f} kN") + print(f" Chamber pressure: {base_inputs.chamber_pressure.to('MPa').value:.0f} MPa") + print(" Propellants: LOX / CH4") + print(f" Mixture ratio: {base_inputs.mixture_ratio}") + print() + + # Get baseline performance and geometry + performance, geometry = design_engine(base_inputs) + + print(f" Ideal Isp (vac): {performance.isp_vac.value:.1f} s") + print(f" Ideal c*: {performance.cstar.to('m/s').value:.0f} m/s") + print(f" Mass flow: {performance.mdot.to('kg/s').value:.1f} kg/s") + + # Store results for comparison + results: list[dict] = [] + + # ========================================================================= + # Cycle 1: Pressure-Fed + # ========================================================================= + + print_header("CYCLE 1: PRESSURE-FED") + + print(" Architecture:") + print(" - Propellants pushed by tank pressure (no pumps)") + print(" - Requires high tank pressure → heavy tanks") + print(" - Simplest system, highest reliability") + print() + + pressure_fed = PressureFedCycle( + tank_pressure_margin=1.3, + line_loss_fraction=0.05, + injector_dp_fraction=0.15, + ) + + pf_result = pressure_fed.analyze(base_inputs, performance, geometry) + + print(" Results:") + print(f" Tank pressure (ox): {pf_result.tank_pressure_ox.to('MPa').value:.1f} MPa") + print(f" Tank pressure (fuel): {pf_result.tank_pressure_fuel.to('MPa').value:.1f} MPa") + print(f" Net Isp: {pf_result.net_isp.to('s').value:.1f} s") + print(f" Net thrust: {pf_result.net_thrust.to('kN').value:.0f} kN") + print(f" Cycle efficiency: {pf_result.cycle_efficiency*100:.1f}%") + if pf_result.warnings: + print(f" Warnings: {len(pf_result.warnings)}") + for w in pf_result.warnings: + print(f" ⚠ {w}") + + results.append({ + "name": "Pressure-Fed", + "net_isp": pf_result.net_isp.to("s").value, + "tank_pressure_MPa": pf_result.tank_pressure_ox.to("MPa").value, + "pump_power_kW": 0, + "efficiency": pf_result.cycle_efficiency, + "simplicity": 1.0, + "reliability": 0.95, + }) + + # ========================================================================= + # Cycle 2: Gas Generator + # ========================================================================= + + print_header("CYCLE 2: GAS GENERATOR") + + print(" Architecture:") + print(" - Small combustor (GG) drives turbine") + print(" - Turbine exhaust dumped overboard (Isp loss)") + print(" - Moderate complexity, proven technology") + print() + + gas_generator = GasGeneratorCycle( + turbine_inlet_temp=kelvin(900), + pump_efficiency_ox=0.70, + pump_efficiency_fuel=0.70, + turbine_efficiency=0.65, + gg_mixture_ratio=0.4, + ) + + gg_result = gas_generator.analyze(base_inputs, performance, geometry) + + total_pump_power_gg = ( + gg_result.pump_power_ox.to("kW").value + + gg_result.pump_power_fuel.to("kW").value + ) + + print(" Results:") + print(f" Turbine mass flow: {gg_result.turbine_mass_flow.to('kg/s').value:.1f} kg/s") + print(f" Net Isp: {gg_result.net_isp.to('s').value:.1f} s") + print(f" Net thrust: {gg_result.net_thrust.to('kN').value:.0f} kN") + print(f" Cycle efficiency: {gg_result.cycle_efficiency*100:.1f}%") + print(f" Isp loss vs ideal: {performance.isp_vac.value - gg_result.net_isp.to('s').value:.1f} s") + if gg_result.warnings: + print(f" Warnings: {len(gg_result.warnings)}") + for w in gg_result.warnings: + print(f" ⚠ {w}") + + results.append({ + "name": "Gas Generator", + "net_isp": gg_result.net_isp.to("s").value, + "tank_pressure_MPa": gg_result.tank_pressure_ox.to("MPa").value if gg_result.tank_pressure_ox else 0.5, + "pump_power_kW": total_pump_power_gg, + "efficiency": gg_result.cycle_efficiency, + "simplicity": 0.6, + "reliability": 0.90, + }) + + # ========================================================================= + # Cycle 3: Staged Combustion (Oxidizer-Rich) + # ========================================================================= + + print_header("CYCLE 3: STAGED COMBUSTION (OX-RICH)") + + print(" Architecture:") + print(" - Preburner runs oxygen-rich") + print(" - ALL flow goes through main chamber (no dump)") + print(" - Highest performance, most complex") + print() + + staged_combustion = StagedCombustionCycle( + preburner_temp=kelvin(750), + pump_efficiency_ox=0.75, + pump_efficiency_fuel=0.75, + turbine_efficiency=0.70, + preburner_mixture_ratio=50.0, + oxidizer_rich=True, + ) + + sc_result = staged_combustion.analyze(base_inputs, performance, geometry) + + total_pump_power_sc = ( + sc_result.pump_power_ox.to("kW").value + + sc_result.pump_power_fuel.to("kW").value + ) + + print(" Results:") + print(f" Turbine mass flow: {sc_result.turbine_mass_flow.to('kg/s').value:.1f} kg/s") + print(f" Net Isp: {sc_result.net_isp.to('s').value:.1f} s") + print(f" Net thrust: {sc_result.net_thrust.to('kN').value:.0f} kN") + print(f" Cycle efficiency: {sc_result.cycle_efficiency*100:.1f}%") + if sc_result.warnings: + print(f" Warnings: {len(sc_result.warnings)}") + for w in sc_result.warnings: + print(f" ⚠ {w}") + + results.append({ + "name": "Staged Combustion", + "net_isp": sc_result.net_isp.to("s").value, + "tank_pressure_MPa": sc_result.tank_pressure_ox.to("MPa").value if sc_result.tank_pressure_ox else 0.5, + "pump_power_kW": total_pump_power_sc, + "efficiency": sc_result.cycle_efficiency, + "simplicity": 0.3, + "reliability": 0.85, + }) + + # ========================================================================= + # Comparison Summary + # ========================================================================= + + print_header("COMPARISON SUMMARY") + + print() + print(f" {'Cycle':<20} {'Net Isp':<12} {'Efficiency':<12} {'Tank P (MPa)':<12}") + print(" " + "-" * 56) + + for r in results: + isp_str = f"{r['net_isp']:.1f} s" + eff_str = f"{r['efficiency']*100:.1f}%" + tank_str = f"{r['tank_pressure_MPa']:.1f}" + print(f" {r['name']:<20} {isp_str:<12} {eff_str:<12} {tank_str:<12}") + + # Compute comparison from actual results + if len(results) >= 3: + isp_gain = results[2]["net_isp"] - results[1]["net_isp"] + isp_gain_pct = 100 * isp_gain / results[1]["net_isp"] if results[1]["net_isp"] > 0 else 0 + + print() + print(" Staged combustion vs Gas Generator:") + print(f" Isp gain: {isp_gain:.1f} s ({isp_gain_pct:.1f}%)") + + print() + + # ========================================================================= + # Save Results + # ========================================================================= + + print_header("SAVING RESULTS") + + with OutputContext("cycle_comparison", include_timestamp=True) as ctx: + # Generate visualizations + print() + print(" Generating visualizations...") + + # 1. Bar chart comparison + fig_bars = plot_cycle_comparison_bars( + results, + metrics=["net_isp", "efficiency", "tank_pressure_MPa", "pump_power_kW"], + title="LOX/CH4 Engine Cycle Comparison", + ) + fig_bars.savefig(ctx.path("cycle_comparison_bars.png"), dpi=150, bbox_inches="tight") + print(" - cycle_comparison_bars.png") + + # 2. Radar chart + fig_radar = plot_cycle_radar( + results, + metrics=["net_isp", "efficiency", "simplicity", "reliability"], + title="Cycle Trade-offs (Normalized)", + ) + fig_radar.savefig(ctx.path("cycle_radar.png"), dpi=150, bbox_inches="tight") + print(" - cycle_radar.png") + + # 3. Trade-off scatter plot + fig_tradeoff = plot_cycle_tradeoff( + results, + x_metric="net_isp", + y_metric="simplicity", + size_metric="pump_power_kW", + title="Performance vs Simplicity Trade Space", + ) + fig_tradeoff.savefig(ctx.path("cycle_tradeoff.png"), dpi=150, bbox_inches="tight") + print(" - cycle_tradeoff.png") + + # Save summary data + ctx.save_summary({ + "requirements": { + "thrust_kN": base_inputs.thrust.to("kN").value, + "chamber_pressure_MPa": base_inputs.chamber_pressure.to("MPa").value, + "propellants": "LOX/CH4", + "mixture_ratio": base_inputs.mixture_ratio, + }, + "ideal_performance": { + "isp_vac_s": performance.isp_vac.value, + "cstar_m_s": performance.cstar.value, + "mdot_kg_s": performance.mdot.value, + }, + "cycles": results, + }) + + print() + print(f" Results saved to: {ctx.output_dir}") + + print() + print("╔" + "═" * 68 + "╗") + print("║" + " " * 24 + "ANALYSIS COMPLETE" + " " * 27 + "║") + print("╚" + "═" * 68 + "╝") + print() + + +if __name__ == "__main__": + main() diff --git a/rocket/examples/optimization.py b/rocket/examples/optimization.py new file mode 100644 index 0000000..5014408 --- /dev/null +++ b/rocket/examples/optimization.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python +"""Multi-objective optimization example for Rocket. + +This example demonstrates how to find Pareto-optimal engine designs +that balance competing objectives: + +1. Maximize Isp (specific impulse) +2. Minimize engine mass (via throat size as proxy) +3. Satisfy thermal constraints + +Real engine design involves tradeoffs - you can't maximize everything. +Pareto fronts show the best achievable combinations. +""" + +from rocket import ( + EngineInputs, + MultiObjectiveOptimizer, + OutputContext, + Range, + design_engine, +) +from rocket.units import kilonewtons, megapascals + + +def main() -> None: + """Run the multi-objective optimization example.""" + print("=" * 70) + print("Rocket - Multi-Objective Optimization Example") + print("=" * 70) + print() + + # ========================================================================= + # Define the Optimization Problem + # ========================================================================= + + print("Problem: Design a LOX/CH4 engine balancing Isp vs. compactness") + print("-" * 70) + print() + print("Objectives:") + print(" 1. Maximize vacuum Isp (performance)") + print(" 2. Minimize throat diameter (smaller = lighter, cheaper)") + print() + print("Design variables:") + print(" - Chamber pressure: 5-25 MPa") + print(" - Mixture ratio: 2.5-4.0") + print() + print("Constraints:") + print(" - Expansion ratio < 80 (practical nozzle size)") + print(" - Throat diameter > 3 cm (manufacturability)") + print() + + # Baseline design + baseline = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(100), + chamber_pressure=megapascals(10), + mixture_ratio=3.2, + name="Methalox-Baseline", + ) + + # Define constraints + def expansion_constraint(result: tuple) -> bool: + """Expansion ratio must be < 80.""" + _, geometry = result + return geometry.expansion_ratio < 80 + + def throat_constraint(result: tuple) -> bool: + """Throat diameter must be > 3 cm for manufacturability.""" + _, geometry = result + return geometry.throat_diameter.to("m").value * 100 > 3.0 + + # ========================================================================= + # Run Multi-Objective Optimization + # ========================================================================= + + print("Running optimization...") + print() + + optimizer = MultiObjectiveOptimizer( + compute=design_engine, + base=baseline, + objectives=["isp_vac", "throat_diameter"], + maximize=[True, False], # Max Isp, Min throat (smaller = better) + vary={ + "chamber_pressure": Range(5, 25, n=15, unit="MPa"), + "mixture_ratio": Range(2.5, 4.0, n=10), + }, + constraints=[expansion_constraint, throat_constraint], + ) + + pareto_results = optimizer.run(progress=True) + + print() + print(f"Total designs evaluated: {pareto_results.n_total}") + print(f"Feasible designs: {pareto_results.n_feasible}") + print(f"Pareto-optimal designs: {pareto_results.n_pareto}") + print() + + # ========================================================================= + # Analyze Pareto Front + # ========================================================================= + + print("=" * 70) + print("PARETO-OPTIMAL DESIGNS") + print("=" * 70) + print() + print("These designs represent the best tradeoffs between Isp and size:") + print() + print(f" {'#':<4} {'Pc (MPa)':<10} {'MR':<8} {'Isp (vac)':<12} {'Throat (cm)':<12} {'ε':<8}") + print(" " + "-" * 54) + + pareto_df = pareto_results.pareto_front() + + for i, row in enumerate(pareto_df.iter_rows(named=True)): + pc_mpa = row["chamber_pressure"] / 1e6 + dt_cm = row["throat_diameter"] * 100 + print(f" {i+1:<4} {pc_mpa:<10.0f} {row['mixture_ratio']:<8.1f} {row['isp_vac']:<12.1f} {dt_cm:<12.2f} {row['expansion_ratio']:<8.1f}") + + print() + + # ========================================================================= + # Interpret Results + # ========================================================================= + + print("=" * 70) + print("DESIGN RECOMMENDATIONS") + print("=" * 70) + print() + + # Best Isp design + best_isp_idx = pareto_df["isp_vac"].arg_max() + best_isp_row = pareto_df.row(best_isp_idx, named=True) + + print("For MAXIMUM PERFORMANCE (highest Isp):") + print(f" Chamber pressure: {best_isp_row['chamber_pressure']/1e6:.0f} MPa") + print(f" Mixture ratio: {best_isp_row['mixture_ratio']:.1f}") + print(f" Isp (vac): {best_isp_row['isp_vac']:.1f} s") + print(f" Throat diameter: {best_isp_row['throat_diameter']*100:.2f} cm") + print() + + # Smallest design + best_size_idx = pareto_df["throat_diameter"].arg_min() + best_size_row = pareto_df.row(best_size_idx, named=True) + + print("For MINIMUM SIZE (smallest throat):") + print(f" Chamber pressure: {best_size_row['chamber_pressure']/1e6:.0f} MPa") + print(f" Mixture ratio: {best_size_row['mixture_ratio']:.1f}") + print(f" Isp (vac): {best_size_row['isp_vac']:.1f} s") + print(f" Throat diameter: {best_size_row['throat_diameter']*100:.2f} cm") + print() + + # Balanced design (middle of Pareto front) + mid_idx = len(pareto_df) // 2 + mid_row = pareto_df.row(mid_idx, named=True) + + print("For BALANCED DESIGN (middle of Pareto front):") + print(f" Chamber pressure: {mid_row['chamber_pressure']/1e6:.0f} MPa") + print(f" Mixture ratio: {mid_row['mixture_ratio']:.1f}") + print(f" Isp (vac): {mid_row['isp_vac']:.1f} s") + print(f" Throat diameter: {mid_row['throat_diameter']*100:.2f} cm") + print() + + # ========================================================================= + # Tradeoff Analysis + # ========================================================================= + + print("=" * 70) + print("TRADEOFF ANALYSIS") + print("=" * 70) + print() + + isp_range = pareto_df["isp_vac"].max() - pareto_df["isp_vac"].min() + dt_range = (pareto_df["throat_diameter"].max() - pareto_df["throat_diameter"].min()) * 100 + + print(f" Isp range on Pareto front: {pareto_df['isp_vac'].min():.1f} - {pareto_df['isp_vac'].max():.1f} s (Δ = {isp_range:.1f} s)") + print(f" Throat range on Pareto front: {pareto_df['throat_diameter'].min()*100:.2f} - {pareto_df['throat_diameter'].max()*100:.2f} cm (Δ = {dt_range:.2f} cm)") + print() + print(" Interpretation:") + print(f" - You can gain up to {isp_range:.1f} s of Isp...") + print(f" - ...at the cost of {dt_range:.1f} cm larger throat") + print(f" - Marginal tradeoff: {isp_range/dt_range:.1f} s of Isp per cm of throat") + print() + + # ========================================================================= + # Save Results + # ========================================================================= + + print("=" * 70) + print("Saving Results") + print("=" * 70) + print() + + with OutputContext("optimization_results", include_timestamp=True) as ctx: + # Export all results + ctx.log("Exporting full design space...") + pareto_results.all_results.to_csv(ctx.path("all_designs.csv")) + + ctx.log("Exporting Pareto front...") + pareto_df.write_csv(ctx.path("pareto_front.csv")) + + # Save summary + ctx.save_summary({ + "problem": { + "objectives": ["isp_vac (maximize)", "throat_diameter (minimize)"], + "design_variables": { + "chamber_pressure_MPa": [5, 25], + "mixture_ratio": [2.5, 4.0], + }, + "constraints": [ + "expansion_ratio < 80", + "throat_diameter > 3 cm", + ], + }, + "results": { + "total_designs": pareto_results.n_total, + "feasible_designs": pareto_results.n_feasible, + "pareto_optimal": pareto_results.n_pareto, + }, + "recommendations": { + "max_performance": { + "chamber_pressure_MPa": best_isp_row["chamber_pressure"] / 1e6, + "mixture_ratio": best_isp_row["mixture_ratio"], + "isp_vac_s": best_isp_row["isp_vac"], + }, + "min_size": { + "chamber_pressure_MPa": best_size_row["chamber_pressure"] / 1e6, + "mixture_ratio": best_size_row["mixture_ratio"], + "throat_diameter_cm": best_size_row["throat_diameter"] * 100, + }, + "balanced": { + "chamber_pressure_MPa": mid_row["chamber_pressure"] / 1e6, + "mixture_ratio": mid_row["mixture_ratio"], + "isp_vac_s": mid_row["isp_vac"], + "throat_diameter_cm": mid_row["throat_diameter"] * 100, + }, + }, + }) + + print() + print(f" All results saved to: {ctx.output_dir}") + + print() + print("=" * 70) + print("Optimization Complete!") + print("=" * 70) + print() + + +if __name__ == "__main__": + main() + diff --git a/rocket/examples/propellant_design.py b/rocket/examples/propellant_design.py index d134f6f..93a0d55 100644 --- a/rocket/examples/propellant_design.py +++ b/rocket/examples/propellant_design.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -"""Propellant-based engine design example for OpenRocketEngine. +"""Propellant-based engine design example for Rocket. This example demonstrates the simplified workflow where you specify propellants and the library automatically determines combustion properties. @@ -7,21 +7,24 @@ No need to manually look up Tc, gamma, or molecular weight! """ -from openrocketengine import design_engine, plot_engine_dashboard -from openrocketengine.engine import EngineInputs -from openrocketengine.nozzle import full_chamber_contour, generate_nozzle_from_geometry -from openrocketengine.units import kilonewtons, megapascals +from rocket import OutputContext, design_engine, plot_engine_dashboard +from rocket.engine import EngineInputs +from rocket.nozzle import full_chamber_contour, generate_nozzle_from_geometry +from rocket.units import kilonewtons, megapascals def main() -> None: """Run the propellant-based design example.""" print("=" * 70) - print("OpenRocketEngine - Propellant-Based Design") + print("Rocket - Propellant-Based Design") print("=" * 70) print() print("Using NASA CEA (via RocketCEA) for thermochemistry calculations") print() + # Store all engine designs for comparison + designs: list[tuple[str, EngineInputs, any, any]] = [] + # ========================================================================= # Design a LOX/RP-1 Engine (like Merlin) # ========================================================================= @@ -62,6 +65,8 @@ def main() -> None: print(f" Expansion Ratio: {geom1.expansion_ratio:.1f}") print() + designs.append(("LOX/RP-1", lox_rp1, perf1, geom1)) + # ========================================================================= # Design a LOX/Methane Engine (like Raptor) # ========================================================================= @@ -91,6 +96,8 @@ def main() -> None: print(f" Isp (Vac): {perf2.isp_vac.value:.1f} s") print() + designs.append(("LOX/CH4", lox_ch4, perf2, geom2)) + # ========================================================================= # Design a LOX/LH2 Engine (like RS-25/SSME) # ========================================================================= @@ -120,6 +127,8 @@ def main() -> None: print(f" Isp (Vac): {perf3.isp_vac.value:.1f} s <- Highest!") print() + designs.append(("LOX/LH2", lox_lh2, perf3, geom3)) + # ========================================================================= # Comparison Summary # ========================================================================= @@ -131,11 +140,7 @@ def main() -> None: print(f"{'Engine':<20} {'Isp(SL)':<10} {'Isp(Vac)':<10} {'MW':<8} {'Tc (K)':<10}") print("-" * 70) - for name, inputs, perf in [ - ("LOX/RP-1", lox_rp1, perf1), - ("LOX/CH4", lox_ch4, perf2), - ("LOX/LH2", lox_lh2, perf3), - ]: + for name, inputs, perf, _ in designs: print( f"{name:<20} " f"{perf.isp.value:<10.1f} " @@ -145,22 +150,46 @@ def main() -> None: ) print() - print("Note: Lower molecular weight (MW) = higher Isp") - print(" LH2 engines have best Isp but require large tanks (low density)") - print() # ========================================================================= - # Generate Dashboard for LOX/RP-1 Engine + # Generate Dashboards for All Engines # ========================================================================= - print("Generating visualization for LOX/RP-1 engine...") - - nozzle = generate_nozzle_from_geometry(geom1) - contour = full_chamber_contour(lox_rp1, geom1, nozzle) + print("=" * 70) + print("GENERATING VISUALIZATIONS") + print("=" * 70) + print() - fig = plot_engine_dashboard(lox_rp1, perf1, geom1, contour) - fig.savefig("kerolox_engine_dashboard.png", dpi=150, bbox_inches="tight") - print(" Saved: kerolox_engine_dashboard.png") + with OutputContext("propellant_comparison", include_timestamp=True) as ctx: + # Save comparison summary + ctx.save_summary({ + "designs": [ + { + "name": name, + "propellants": f"{inputs.name}", + "isp_sl": perf.isp.value, + "isp_vac": perf.isp_vac.value, + "molecular_weight": inputs.molecular_weight, + "chamber_temp_K": inputs.chamber_temp.to("K").value, + "gamma": inputs.gamma, + "thrust_kN": inputs.thrust.to("kN").value, + "chamber_pressure_MPa": inputs.chamber_pressure.to("MPa").value, + } + for name, inputs, perf, _ in designs + ] + }, "comparison_summary.json") + + for _name, inputs, perf, geom in designs: + ctx.log(f"Generating dashboard for {inputs.name}...") + nozzle = generate_nozzle_from_geometry(geom) + contour = full_chamber_contour(inputs, geom, nozzle) + fig = plot_engine_dashboard(inputs, perf, geom, contour) + + safe_name = inputs.name.lower().replace("-", "_").replace(" ", "_") + fig.savefig(ctx.path(f"{safe_name}_dashboard.png"), dpi=150, bbox_inches="tight") + + print() + print(f" All outputs saved to: {ctx.output_dir}") print() print("=" * 70) @@ -170,4 +199,3 @@ def main() -> None: if __name__ == "__main__": main() - diff --git a/rocket/examples/thermal_analysis.py b/rocket/examples/thermal_analysis.py new file mode 100644 index 0000000..feba325 --- /dev/null +++ b/rocket/examples/thermal_analysis.py @@ -0,0 +1,275 @@ +#!/usr/bin/env python +"""Thermal analysis example for Rocket. + +This example demonstrates thermal screening to determine if an engine +design can be regeneratively cooled: + +1. Estimate heat flux using the Bartz correlation +2. Check cooling feasibility for different coolants +3. Understand the relationship between chamber pressure and cooling + +High chamber pressure engines (like Raptor) face severe thermal challenges. +This analysis helps catch infeasible designs early. +""" + +from rocket import ( + EngineInputs, + OutputContext, + design_engine, +) +from rocket.thermal import ( + check_cooling_feasibility, + estimate_heat_flux, + heat_flux_profile, +) +from rocket.units import kelvin, kilonewtons, megapascals + + +def print_header(text: str) -> None: + """Print a formatted section header.""" + print() + print("┌" + "─" * 68 + "┐") + print(f"│ {text:<66} │") + print("└" + "─" * 68 + "┘") + + +def main() -> None: + """Run the thermal analysis example.""" + print() + print("╔" + "═" * 68 + "╗") + print("║" + " " * 18 + "THERMAL FEASIBILITY ANALYSIS" + " " * 22 + "║") + print("║" + " " * 15 + "Can This Engine Be Regeneratively Cooled?" + " " * 12 + "║") + print("╚" + "═" * 68 + "╝") + + # ========================================================================= + # Design a High-Performance Engine + # ========================================================================= + + print_header("ENGINE DESIGN") + + # High chamber pressure like Raptor + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(500), + chamber_pressure=megapascals(25), # High pressure! + mixture_ratio=3.2, + name="High-Pc-Methalox", + ) + + performance, geometry = design_engine(inputs) + + print(f" Engine: {inputs.name}") + print(f" Thrust: {inputs.thrust.to('kN').value:.0f} kN") + print(f" Chamber pressure: {inputs.chamber_pressure.to('MPa').value:.0f} MPa") + print(f" Chamber temperature: {inputs.chamber_temp.to('K').value:.0f} K") + print() + print(f" Isp (vac): {performance.isp_vac.value:.1f} s") + print(f" Throat diameter: {geometry.throat_diameter.to('m').value*100:.2f} cm") + + # ========================================================================= + # Heat Flux Estimation + # ========================================================================= + + print_header("HEAT FLUX ANALYSIS") + + print(" Using Bartz correlation for convective heat transfer...") + print() + + # Estimate heat flux at key locations + locations = ["chamber", "throat", "exit"] + + print(f" {'Location':<15} {'Heat Flux (MW/m²)':<20}") + print(" " + "-" * 35) + + max_q = 0.0 + max_location = "" + + for location in locations: + q = estimate_heat_flux( + inputs=inputs, + performance=performance, + geometry=geometry, + location=location, + ) + q_mw = q.to("W/m^2").value / 1e6 + + if q_mw > max_q: + max_q = q_mw + max_location = location + + print(f" {location:<15} {q_mw:<20.1f}") + + print() + print(f" Maximum heat flux: {max_q:.1f} MW/m² at {max_location}") + + # ========================================================================= + # Heat Flux Profile Along Nozzle + # ========================================================================= + + print_header("AXIAL HEAT FLUX PROFILE") + + x_positions, heat_fluxes = heat_flux_profile( + inputs=inputs, + performance=performance, + geometry=geometry, + n_points=11, + ) + + print(f" {'x/L':<10} {'Heat Flux (MW/m²)':<20}") + print(" " + "-" * 30) + + for x, q in zip(x_positions, heat_fluxes, strict=True): + q_mw = q / 1e6 # heat_flux_profile returns W/m² + bar = "█" * int(q_mw / 5) # Simple bar chart + print(f" {x:<10.2f} {q_mw:<10.1f} {bar}") + + # ========================================================================= + # Cooling Feasibility Check - Methane + # ========================================================================= + + print_header("COOLING FEASIBILITY: METHANE (CH4)") + + print(" Checking if methane can cool this engine...") + print() + + ch4_cooling = check_cooling_feasibility( + inputs=inputs, + performance=performance, + geometry=geometry, + coolant="CH4", + coolant_inlet_temp=kelvin(110), # Near boiling point + max_wall_temp=kelvin(920), # Material limit + ) + + if ch4_cooling.feasible: + print(" ✓ FEASIBLE with methane cooling!") + else: + print(" ✗ NOT FEASIBLE with methane cooling") + + print() + print(f" Max wall temperature: {ch4_cooling.max_wall_temp.to('K').value:.0f} K (limit: {ch4_cooling.max_allowed_temp.to('K').value:.0f} K)") + print(f" Throat heat flux: {ch4_cooling.throat_heat_flux.to('W/m^2').value/1e6:.1f} MW/m²") + print(f" Coolant outlet temp: {ch4_cooling.coolant_outlet_temp.to('K').value:.0f} K") + print(f" Flow margin: {ch4_cooling.flow_margin:.2f}x") + if ch4_cooling.warnings: + print(" Warnings:") + for w in ch4_cooling.warnings: + print(f" ⚠ {w}") + + # ========================================================================= + # Cooling Feasibility Check - RP-1 (for comparison) + # ========================================================================= + + print_header("COOLING FEASIBILITY: RP-1 (KEROSENE)") + + # Design a kerosene engine for comparison + rp1_inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="RP1", + thrust=kilonewtons(500), + chamber_pressure=megapascals(25), + mixture_ratio=2.7, + name="High-Pc-Kerolox", + ) + rp1_perf, rp1_geom = design_engine(rp1_inputs) + + print(" Checking RP-1 cooling for LOX/RP-1 engine...") + print() + + rp1_cooling = check_cooling_feasibility( + inputs=rp1_inputs, + performance=rp1_perf, + geometry=rp1_geom, + coolant="RP1", + coolant_inlet_temp=kelvin(300), # Room temperature + max_wall_temp=kelvin(920), + ) + + if rp1_cooling.feasible: + print(" ✓ FEASIBLE with RP-1 cooling!") + else: + print(" ✗ NOT FEASIBLE with RP-1 cooling") + + print() + print(f" Max wall temperature: {rp1_cooling.max_wall_temp.to('K').value:.0f} K") + print(f" Coolant outlet temp: {rp1_cooling.coolant_outlet_temp.to('K').value:.0f} K") + print(f" Flow margin: {rp1_cooling.flow_margin:.2f}x") + + # ========================================================================= + # Chamber Pressure Impact Study + # ========================================================================= + + print_header("CHAMBER PRESSURE vs. COOLING FEASIBILITY") + + print(" How does Pc affect cooling feasibility?") + print() + print(f" {'Pc (MPa)':<12} {'Throat q (MW/m²)':<20} {'Coolable?':<12}") + print(" " + "-" * 44) + + for pc_mpa in [5, 10, 15, 20, 25, 30]: + test_inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(500), + chamber_pressure=megapascals(pc_mpa), + mixture_ratio=3.2, + name=f"Test-{pc_mpa}MPa", + ) + test_perf, test_geom = design_engine(test_inputs) + + q_throat = estimate_heat_flux(test_inputs, test_perf, test_geom, location="throat") + q_mw = q_throat.to("W/m^2").value / 1e6 + + cooling = check_cooling_feasibility( + test_inputs, test_perf, test_geom, + coolant="CH4", + coolant_inlet_temp=kelvin(110), + max_wall_temp=kelvin(920), + ) + + status = "✓ Yes" if cooling.feasible else "✗ No" + print(f" {pc_mpa:<12} {q_mw:<20.1f} {status:<12}") + + # ========================================================================= + # Save Results + # ========================================================================= + + print_header("SAVING RESULTS") + + with OutputContext("thermal_analysis", include_timestamp=True) as ctx: + ctx.save_summary({ + "engine": { + "name": inputs.name, + "thrust_kN": inputs.thrust.to("kN").value, + "chamber_pressure_MPa": inputs.chamber_pressure.to("MPa").value, + "chamber_temp_K": inputs.chamber_temp.to("K").value, + }, + "heat_flux": { + "max_MW_m2": max_q, + "max_location": max_location, + }, + "ch4_cooling": { + "feasible": ch4_cooling.feasible, + "max_wall_temp_K": ch4_cooling.max_wall_temp.value, + "flow_margin": ch4_cooling.flow_margin, + }, + "rp1_cooling": { + "feasible": rp1_cooling.feasible, + "max_wall_temp_K": rp1_cooling.max_wall_temp.value, + "flow_margin": rp1_cooling.flow_margin, + }, + }) + + print() + print(f" Results saved to: {ctx.output_dir}") + + print() + print("╔" + "═" * 68 + "╗") + print("║" + " " * 24 + "ANALYSIS COMPLETE" + " " * 27 + "║") + print("╚" + "═" * 68 + "╝") + print() + + +if __name__ == "__main__": + main() diff --git a/rocket/examples/trade_study.py b/rocket/examples/trade_study.py new file mode 100644 index 0000000..4f7d591 --- /dev/null +++ b/rocket/examples/trade_study.py @@ -0,0 +1,254 @@ +#!/usr/bin/env python +"""Parametric trade study example for Rocket. + +This example demonstrates how to run systematic trade studies to explore +the design space and make informed engineering decisions: + +1. Chamber pressure sweeps (most impactful parameter) +2. Multi-parameter grid studies +3. Constraint-based filtering +4. Polars DataFrame export +5. Mixture ratio studies (requires CEA recalculation) + +These tools let you answer questions like: +- "How does Isp change with chamber pressure?" +- "Which designs satisfy my throat size requirement?" +- "What's the tradeoff between performance and size?" +""" + +from rocket import ( + EngineInputs, + OutputContext, + ParametricStudy, + Range, + design_engine, +) +from rocket.units import kilonewtons, megapascals + + +def main() -> None: + """Run the parametric trade study example.""" + print("=" * 70) + print("Rocket - Parametric Trade Study Example") + print("=" * 70) + print() + + # ========================================================================= + # Baseline Engine Design + # ========================================================================= + + print("Baseline Engine: LOX/CH4 Methalox") + print("-" * 70) + + baseline = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(200), + chamber_pressure=megapascals(10), + mixture_ratio=3.2, + name="Methalox-Baseline", + ) + + perf, geom = design_engine(baseline) + print(f" Isp (vac): {perf.isp_vac.value:.1f} s") + print(f" Thrust: {baseline.thrust.to('kN').value:.0f} kN") + print() + + # ========================================================================= + # Study 1: Mixture Ratio Trade (with CEA recalculation) + # ========================================================================= + + print("=" * 70) + print("Study 1: Mixture Ratio Sweep (with CEA thermochemistry)") + print("=" * 70) + print() + print("Question: What mixture ratio maximizes Isp for LOX/CH4?") + print() + print("Note: Each point recalculates combustion properties via NASA CEA") + print() + + mixture_ratios = [2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0] + mr_results = [] + + print(f" {'MR':<8} {'Tc (K)':<10} {'Isp (vac)':<12} {'Isp (SL)':<12}") + print(" " + "-" * 42) + + for mr in mixture_ratios: + inputs = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="CH4", + thrust=kilonewtons(200), + chamber_pressure=megapascals(10), + mixture_ratio=mr, + ) + perf, geom = design_engine(inputs) + mr_results.append({ + "mr": mr, + "Tc": inputs.chamber_temp.to("K").value, + "isp_vac": perf.isp_vac.value, + "isp_sl": perf.isp.value, + }) + print(f" {mr:<8.1f} {inputs.chamber_temp.to('K').value:<10.0f} {perf.isp_vac.value:<12.1f} {perf.isp.value:<12.1f}") + + # Find optimal MR + best = max(mr_results, key=lambda x: x["isp_vac"]) + print() + print(f" → Optimal MR for max Isp: {best['mr']:.1f} (Isp = {best['isp_vac']:.1f} s)") + print(f" At this MR: Tc = {best['Tc']:.0f} K") + print() + + # ========================================================================= + # Study 2: Chamber Pressure Sweep with Range + # ========================================================================= + + print("=" * 70) + print("Study 2: Chamber Pressure Trade") + print("=" * 70) + print() + print("Question: How does performance scale with chamber pressure?") + print() + + # Using Range for continuous sweeps with units + pc_study = ParametricStudy( + compute=design_engine, + base=baseline, + vary={"chamber_pressure": Range(5, 25, n=9, unit="MPa")}, + ) + + pc_results = pc_study.run(progress=True) + + print() + print("Results:") + print(f" {'Pc (MPa)':<10} {'Isp (vac)':<12} {'Throat (cm)':<12} {'c* (m/s)':<10}") + print(" " + "-" * 44) + + df = pc_results.to_dataframe() + for row in df.iter_rows(named=True): + throat_cm = row["throat_diameter"] * 100 + # chamber_pressure is already in MPa (unit specified in Range) + print(f" {row['chamber_pressure']:<10.0f} {row['isp_vac']:<12.1f} {throat_cm:<12.1f} {row['cstar']:<10.0f}") + + # Compute actual insights from data + print() + pc_low = df["chamber_pressure"].min() + pc_high = df["chamber_pressure"].max() + isp_at_low = df.filter(df["chamber_pressure"] == pc_low)["isp_vac"][0] + isp_at_high = df.filter(df["chamber_pressure"] == pc_high)["isp_vac"][0] + dt_at_low = df.filter(df["chamber_pressure"] == pc_low)["throat_diameter"][0] * 100 + dt_at_high = df.filter(df["chamber_pressure"] == pc_high)["throat_diameter"][0] * 100 + + isp_change = isp_at_high - isp_at_low + dt_change = dt_at_high - dt_at_low + + print(f" From {pc_low:.0f} to {pc_high:.0f} MPa:") + print(f" Isp changed by {isp_change:+.1f} s ({100*isp_change/isp_at_low:+.1f}%)") + print(f" Throat diameter changed by {dt_change:+.1f} cm ({100*dt_change/dt_at_low:+.1f}%)") + print() + + # ========================================================================= + # Study 3: Multi-Parameter Grid with Constraints + # ========================================================================= + + print("=" * 70) + print("Study 3: Multi-Parameter Design Space Exploration") + print("=" * 70) + print() + print("Sweeping: Chamber Pressure (5-20 MPa) × Contraction Ratio (3-6)") + print("Constraint: Throat diameter > 8 cm (manufacturability)") + print() + + def throat_constraint(result: tuple) -> bool: + """Filter out designs with throat diameter < 8 cm.""" + _, geometry = result + return geometry.throat_diameter.to("m").value * 100 > 8.0 + + grid_study = ParametricStudy( + compute=design_engine, + base=baseline, + vary={ + "chamber_pressure": Range(5, 20, n=6, unit="MPa"), + "contraction_ratio": [3.0, 4.0, 5.0, 6.0], + }, + constraints=[throat_constraint], + ) + + grid_results = grid_study.run(progress=True) + + print() + n_total = len(grid_results.inputs) + n_feasible = grid_results.constraints_passed.sum() + print(f" Total designs evaluated: {n_total}") + print(f" Feasible designs: {n_feasible} ({100*n_feasible/n_total:.0f}%)") + print() + + # Export to Polars DataFrame for further analysis + df = grid_results.to_dataframe() + + # Filter to feasible designs and show best by Isp + feasible_df = df.filter(df["feasible"]) + best_designs = feasible_df.sort("isp_vac", descending=True).head(5) + + print("Top 5 Feasible Designs by Vacuum Isp:") + print(f" {'Pc (MPa)':<10} {'CR':<8} {'Isp (vac)':<12} {'Dt (cm)':<10}") + print(" " + "-" * 40) + + for row in best_designs.iter_rows(named=True): + # chamber_pressure is already in MPa (unit specified in Range) + dt_cm = row["throat_diameter"] * 100 + print(f" {row['chamber_pressure']:<10.0f} {row['contraction_ratio']:<8.1f} {row['isp_vac']:<12.1f} {dt_cm:<10.1f}") + + print() + + # ========================================================================= + # Save All Results + # ========================================================================= + + print("=" * 70) + print("Saving Results") + print("=" * 70) + print() + + with OutputContext("trade_study_results", include_timestamp=True) as ctx: + # Export DataFrames to CSV + ctx.log("Exporting chamber pressure study...") + pc_results.to_csv(ctx.path("chamber_pressure_sweep.csv")) + + ctx.log("Exporting grid study...") + grid_results.to_csv(ctx.path("grid_study.csv")) + + # Save summary + ctx.save_summary({ + "studies": { + "mixture_ratio_sweep": { + "parameter": "mixture_ratio", + "range": [2.4, 4.0], + "optimal_mr": best["mr"], + "optimal_isp_vac": best["isp_vac"], + "note": "Each point recalculates CEA thermochemistry", + }, + "chamber_pressure_sweep": { + "parameter": "chamber_pressure", + "range_MPa": [5, 25], + "n_points": 9, + }, + "grid_study": { + "parameters": ["chamber_pressure", "contraction_ratio"], + "total_designs": n_total, + "feasible_designs": int(n_feasible), + "constraint": "throat_diameter > 8 cm", + }, + } + }) + + print() + print(f" All results saved to: {ctx.output_dir}") + + print() + print("=" * 70) + print("Trade Study Complete!") + print("=" * 70) + + +if __name__ == "__main__": + main() + diff --git a/rocket/examples/uncertainty_analysis.py b/rocket/examples/uncertainty_analysis.py new file mode 100644 index 0000000..0264d3f --- /dev/null +++ b/rocket/examples/uncertainty_analysis.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python +"""Uncertainty analysis example for Rocket. + +This example demonstrates Monte Carlo uncertainty quantification to understand +how input uncertainties propagate through engine design calculations: + +1. Define probability distributions for uncertain inputs +2. Run Monte Carlo sampling +3. Analyze output statistics and confidence intervals +4. Identify which inputs drive the most uncertainty + +This answers questions like: +- "If my mixture ratio varies ±5%, how much does Isp vary?" +- "What's the 95% confidence interval on my thrust coefficient?" +- "Which input uncertainty should I focus on reducing?" +""" + +from rocket import ( + EngineInputs, + Normal, + OutputContext, + UncertaintyAnalysis, + Uniform, + design_engine, +) +from rocket.units import kilonewtons, megapascals + + +def main() -> None: + """Run the uncertainty analysis example.""" + print("=" * 70) + print("Rocket - Uncertainty Analysis Example") + print("=" * 70) + print() + + # ========================================================================= + # Define Nominal Design Point + # ========================================================================= + + print("Nominal Engine Design: LOX/RP-1 Kerolox") + print("-" * 70) + + nominal = EngineInputs.from_propellants( + oxidizer="LOX", + fuel="RP1", + thrust=kilonewtons(100), + chamber_pressure=megapascals(7), + mixture_ratio=2.7, + name="Kerolox-100", + ) + + perf, geom = design_engine(nominal) + print(f" Nominal Isp (vac): {perf.isp_vac.value:.1f} s") + print(f" Nominal Isp (SL): {perf.isp.value:.1f} s") + print(f" Nominal Thrust Coeff: {perf.thrust_coeff:.3f}") + print(f" Nominal Throat Dia: {geom.throat_diameter.to('m').value*100:.2f} cm") + print() + + # ========================================================================= + # Study 1: Single Source of Uncertainty (Mixture Ratio) + # ========================================================================= + + print("=" * 70) + print("Study 1: Mixture Ratio Uncertainty") + print("=" * 70) + print() + print("The mixture ratio is controlled by the propellant valves.") + print("Assume MR = 2.7 ± 0.1 (Normal distribution, σ = 0.1)") + print() + + mr_uncertainty = UncertaintyAnalysis( + compute=design_engine, + base=nominal, + distributions={ + "mixture_ratio": Normal(mean=2.7, std=0.1), + }, + seed=42, # For reproducibility + ) + + mr_results = mr_uncertainty.run(n_samples=1000, progress=True) + + print() + print("Monte Carlo Results (N=1000):") + print() + + # Get statistics for key metrics + stats = mr_results.statistics() + + print(f" {'Metric':<20} {'Mean':<12} {'Std Dev':<12} {'95% CI':<20}") + print(" " + "-" * 64) + + for metric in ["isp_vac", "isp", "thrust_coeff"]: + mean = stats[metric]["mean"] + std = stats[metric]["std"] + ci_low, ci_high = mr_results.confidence_interval(metric, 0.95) + print(f" {metric:<20} {mean:<12.2f} {std:<12.4f} [{ci_low:.2f}, {ci_high:.2f}]") + + print() + print(" Interpretation:") + print(f" - MR uncertainty of σ=0.1 causes Isp variation of σ≈{stats['isp_vac']['std']:.2f} s") + print(f" - 95% of the time, Isp(vac) will be in [{mr_results.confidence_interval('isp_vac', 0.95)[0]:.1f}, {mr_results.confidence_interval('isp_vac', 0.95)[1]:.1f}] s") + print() + + # ========================================================================= + # Study 2: Multiple Sources of Uncertainty + # ========================================================================= + + print("=" * 70) + print("Study 2: Multiple Uncertainty Sources") + print("=" * 70) + print() + print("Real engines have uncertainty in multiple inputs:") + print(" - Chamber pressure: Pc = 7 MPa ± 0.3 MPa (Normal)") + print(" - Mixture ratio: MR = 2.7 ± 0.1 (Normal)") + print(" - Gamma: γ = 1.24 ± 0.02 (Normal, combustion model uncertainty)") + print() + + multi_uncertainty = UncertaintyAnalysis( + compute=design_engine, + base=nominal, + distributions={ + "chamber_pressure": Normal(mean=7, std=0.3, unit="MPa"), + "mixture_ratio": Normal(mean=2.7, std=0.1), + "gamma": Normal(mean=nominal.gamma, std=0.02), + }, + seed=42, + ) + + multi_results = multi_uncertainty.run(n_samples=2000, progress=True) + + print() + print("Monte Carlo Results (N=2000):") + print() + + stats = multi_results.statistics() + + print(f" {'Metric':<20} {'Mean':<12} {'Std Dev':<12} {'Range (99%)':<20}") + print(" " + "-" * 64) + + for metric in ["isp_vac", "isp", "thrust_coeff", "cstar", "expansion_ratio"]: + mean = stats[metric]["mean"] + std = stats[metric]["std"] + ci_low, ci_high = multi_results.confidence_interval(metric, 0.99) + print(f" {metric:<20} {mean:<12.2f} {std:<12.4f} [{ci_low:.2f}, {ci_high:.2f}]") + + print() + + # ========================================================================= + # Study 3: Uniform Distribution (Manufacturing Tolerances) + # ========================================================================= + + print("=" * 70) + print("Study 3: Manufacturing Tolerance Analysis") + print("=" * 70) + print() + print("Manufacturing tolerances are often uniformly distributed.") + print(" - Contraction ratio: CR = 4.0 ± 0.2 (Uniform)") + print(" - L* (characteristic length): L* = 1.0 ± 0.05 m (Uniform)") + print() + + mfg_uncertainty = UncertaintyAnalysis( + compute=design_engine, + base=nominal, + distributions={ + "contraction_ratio": Uniform(low=3.8, high=4.2), + "lstar": Uniform(low=0.95, high=1.05, unit="m"), + }, + seed=42, + ) + + mfg_results = mfg_uncertainty.run(n_samples=1000, progress=True) + + print() + print("Impact of Manufacturing Tolerances:") + print() + + stats = mfg_results.statistics() + + # These parameters mainly affect geometry, not performance + for metric in ["chamber_diameter", "chamber_length"]: + mean = stats[metric]["mean"] + std = stats[metric]["std"] + cv = 100 * std / mean # Coefficient of variation + print(f" {metric:<20}: mean={mean*100:.2f} cm, CV={cv:.2f}%") + + print() + print(" Note: Geometric tolerances have minimal impact on performance") + print(" but affect fit/interface dimensions.") + print() + + # ========================================================================= + # Export Results + # ========================================================================= + + print("=" * 70) + print("Saving Results") + print("=" * 70) + print() + + with OutputContext("uncertainty_analysis", include_timestamp=True) as ctx: + # Export raw Monte Carlo samples + ctx.log("Exporting MR uncertainty samples...") + mr_results.to_csv(ctx.path("mr_uncertainty_samples.csv")) + + ctx.log("Exporting multi-source uncertainty samples...") + multi_results.to_csv(ctx.path("multi_uncertainty_samples.csv")) + + ctx.log("Exporting manufacturing tolerance samples...") + mfg_results.to_csv(ctx.path("mfg_tolerance_samples.csv")) + + # Save statistical summary + ctx.save_summary({ + "mr_uncertainty": { + "inputs": {"mixture_ratio": {"distribution": "Normal", "mean": 2.7, "std": 0.1}}, + "n_samples": 1000, + "isp_vac": { + "mean": float(mr_results.statistics()["isp_vac"]["mean"]), + "std": float(mr_results.statistics()["isp_vac"]["std"]), + "ci_95": list(mr_results.confidence_interval("isp_vac", 0.95)), + }, + }, + "multi_uncertainty": { + "inputs": { + "chamber_pressure": {"distribution": "Normal", "mean_MPa": 7, "std_MPa": 0.3}, + "mixture_ratio": {"distribution": "Normal", "mean": 2.7, "std": 0.1}, + "gamma": {"distribution": "Normal", "mean": nominal.gamma, "std": 0.02}, + }, + "n_samples": 2000, + }, + }) + + print() + print(f" All results saved to: {ctx.output_dir}") + + print() + print("=" * 70) + print("Uncertainty Analysis Complete!") + print("=" * 70) + print() + + +if __name__ == "__main__": + main() + diff --git a/rocket/output.py b/rocket/output.py new file mode 100644 index 0000000..122b8b1 --- /dev/null +++ b/rocket/output.py @@ -0,0 +1,268 @@ +"""Output management for Rocket. + +This module provides utilities for organizing outputs from rocket design +analyses into structured directories with consistent naming. + +Example: + >>> from rocket.output import OutputContext + >>> with OutputContext("my_engine_study") as ctx: + ... fig.savefig(ctx.path("engine_dashboard.png")) + ... contour.to_csv(ctx.path("nozzle_contour.csv")) + ... ctx.save_summary({"isp": 300, "thrust": 50000}) +""" + +import json +import shutil +from datetime import datetime +from pathlib import Path +from typing import Any + +from beartype import beartype + + +@beartype +class OutputContext: + """Context manager for organizing analysis outputs. + + Creates a structured output directory with: + - Timestamp-based naming for version control + - Subdirectories for different output types + - Automatic metadata logging + - Summary JSON export + + Directory structure: + {base_dir}/{name}_{timestamp}/ + ├── plots/ # Visualization outputs + ├── data/ # CSV, contour exports + ├── reports/ # Text summaries + └── metadata.json # Run information + + Attributes: + name: Study/analysis name + output_dir: Path to the output directory + timestamp: Creation timestamp + """ + + def __init__( + self, + name: str, + base_dir: str | Path | None = None, + include_timestamp: bool = True, + create_subdirs: bool = True, + ) -> None: + """Initialize output context. + + Args: + name: Name for this analysis/study (used in directory name) + base_dir: Base directory for outputs. Defaults to ./outputs/ + include_timestamp: Whether to include timestamp in directory name + create_subdirs: Whether to create plots/, data/, reports/ subdirs + """ + self.name = name + self.timestamp = datetime.now() + self._include_timestamp = include_timestamp + self._create_subdirs = create_subdirs + self._metadata: dict[str, Any] = { + "name": name, + "created": self.timestamp.isoformat(), + "files": [], + } + + # Determine base directory + if base_dir is None: + base_dir = Path.cwd() / "outputs" + self.base_dir = Path(base_dir) + + # Create output directory name + if include_timestamp: + timestamp_str = self.timestamp.strftime("%Y%m%d_%H%M%S") + dir_name = f"{name}_{timestamp_str}" + else: + dir_name = name + + self.output_dir = self.base_dir / dir_name + self._entered = False + + def __enter__(self) -> "OutputContext": + """Enter the context and create directories.""" + self._entered = True + + # Create main output directory + self.output_dir.mkdir(parents=True, exist_ok=True) + + # Create subdirectories + if self._create_subdirs: + (self.output_dir / "plots").mkdir(exist_ok=True) + (self.output_dir / "data").mkdir(exist_ok=True) + (self.output_dir / "reports").mkdir(exist_ok=True) + + return self + + def __exit__(self, exc_type: Any, exc_val: Any, exc_tb: Any) -> None: + """Exit context and write metadata.""" + self._metadata["completed"] = datetime.now().isoformat() + self._metadata["success"] = exc_type is None + + # Write metadata + metadata_path = self.output_dir / "metadata.json" + with open(metadata_path, "w") as f: + json.dump(self._metadata, f, indent=2, default=str) + + self._entered = False + + def path(self, filename: str, subdir: str | None = None) -> Path: + """Get full path for an output file. + + Automatically routes files to appropriate subdirectories based on extension. + + Args: + filename: Name of the output file + subdir: Optional subdirectory override. If None, auto-routes: + - .png, .pdf, .svg → plots/ + - .csv, .json → data/ + - .txt, .md → reports/ + + Returns: + Full path to the output file + """ + if not self._entered: + raise RuntimeError("OutputContext must be used as a context manager") + + # Auto-route based on extension if subdir not specified + if subdir is None and self._create_subdirs: + ext = Path(filename).suffix.lower() + if ext in {".png", ".pdf", ".svg", ".jpg", ".jpeg"}: + subdir = "plots" + elif ext in {".csv", ".json", ".npy", ".npz"}: + subdir = "data" + elif ext in {".txt", ".md", ".rst", ".log"}: + subdir = "reports" + + full_path = self.output_dir / subdir / filename if subdir else self.output_dir / filename + + # Track file for metadata + self._metadata["files"].append(str(full_path.relative_to(self.output_dir))) + + return full_path + + def plots_dir(self) -> Path: + """Get path to plots subdirectory.""" + return self.output_dir / "plots" + + def data_dir(self) -> Path: + """Get path to data subdirectory.""" + return self.output_dir / "data" + + def reports_dir(self) -> Path: + """Get path to reports subdirectory.""" + return self.output_dir / "reports" + + def save_summary(self, summary: dict[str, Any], filename: str = "summary.json") -> Path: + """Save a summary dictionary as JSON. + + Args: + summary: Dictionary of summary data + filename: Output filename + + Returns: + Path to saved file + """ + path = self.path(filename, subdir="data") + with open(path, "w") as f: + json.dump(summary, f, indent=2, default=str) + return path + + def save_text(self, text: str, filename: str = "report.txt") -> Path: + """Save text to a report file. + + Args: + text: Text content to save + filename: Output filename + + Returns: + Path to saved file + """ + path = self.path(filename, subdir="reports") + with open(path, "w") as f: + f.write(text) + return path + + def add_metadata(self, key: str, value: Any) -> None: + """Add custom metadata. + + Args: + key: Metadata key + value: Metadata value (must be JSON-serializable) + """ + self._metadata[key] = value + + def log(self, message: str) -> None: + """Log a message to both console and log file. + + Args: + message: Message to log + """ + timestamp = datetime.now().strftime("%H:%M:%S") + formatted = f"[{timestamp}] {message}" + print(formatted) + + log_path = self.output_dir / "run.log" + with open(log_path, "a") as f: + f.write(formatted + "\n") + + +@beartype +def get_default_output_dir() -> Path: + """Get the default output directory. + + Returns ./outputs/ in the current working directory. + + Returns: + Path to default output directory + """ + return Path.cwd() / "outputs" + + +@beartype +def list_outputs(base_dir: str | Path | None = None) -> list[Path]: + """List all output directories. + + Args: + base_dir: Base directory to search. Defaults to ./outputs/ + + Returns: + List of output directory paths, sorted by modification time (newest first) + """ + if base_dir is None: + base_dir = get_default_output_dir() + base_dir = Path(base_dir) + + if not base_dir.exists(): + return [] + + outputs = [d for d in base_dir.iterdir() if d.is_dir()] + return sorted(outputs, key=lambda p: p.stat().st_mtime, reverse=True) + + +@beartype +def clean_outputs(base_dir: str | Path | None = None, keep_latest: int = 5) -> int: + """Clean old output directories, keeping the N most recent. + + Args: + base_dir: Base directory to clean. Defaults to ./outputs/ + keep_latest: Number of recent outputs to keep + + Returns: + Number of directories removed + """ + outputs = list_outputs(base_dir) + + if len(outputs) <= keep_latest: + return 0 + + to_remove = outputs[keep_latest:] + for output_dir in to_remove: + shutil.rmtree(output_dir) + + return len(to_remove) + diff --git a/rocket/plotting.py b/rocket/plotting.py index 47c0d20..17387e3 100644 --- a/rocket/plotting.py +++ b/rocket/plotting.py @@ -5,6 +5,7 @@ - Performance curves (Isp vs altitude, thrust vs altitude) - Nozzle contour visualization - Trade study plots +- Cycle comparison charts All plots use matplotlib with a consistent, professional style. """ @@ -658,64 +659,366 @@ def plot_engine_dashboard( return fig +# ============================================================================= +# Mass Breakdown Plot +# ============================================================================= + + @beartype def plot_mass_breakdown( masses: dict[str, float | int], - title: str = "Vehicle Mass Breakdown", + title: str = "Mass Breakdown", + figsize: tuple[float, float] = (10, 8), ) -> Figure: - """Plot a mass breakdown pie chart and bar chart. + """Create a mass breakdown pie chart with bar chart. Args: - masses: Dictionary of component names to masses in kg + masses: Dictionary mapping component names to masses (kg) title: Plot title + figsize: Figure size Returns: matplotlib Figure """ _setup_style() - fig, (ax_pie, ax_bar) = plt.subplots(1, 2, figsize=(14, 6)) - # Sort by mass - sorted_items = sorted(masses.items(), key=lambda x: x[1], reverse=True) - labels = [item[0] for item in sorted_items] - values = [item[1] for item in sorted_items] + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=figsize) + + labels = list(masses.keys()) + values = list(masses.values()) total = sum(values) # Color palette - colors = plt.cm.Set3(np.linspace(0, 1, len(labels))) + colors = plt.cm.Set2(np.linspace(0, 1, len(labels))) # Pie chart - wedges, texts, autotexts = ax_pie.pie( - values, labels=None, autopct=lambda p: f"{p:.1f}%" if p > 3 else "", - colors=colors, startangle=90, counterclock=False, - wedgeprops=dict(linewidth=2, edgecolor="white") + ax1.pie( + values, + labels=labels, + autopct=lambda pct: f"{pct:.1f}%\n({pct*total/100:.0f} kg)", + colors=colors, + startangle=90, + pctdistance=0.75, ) - ax_pie.set_title("Mass Distribution", fontsize=12, fontweight="bold") - - # Legend for pie chart - legend_labels = [f"{label}: {v:,.0f} kg" for label, v in zip(labels, values, strict=True)] - ax_pie.legend(wedges, legend_labels, loc="center left", bbox_to_anchor=(1, 0.5)) + ax1.set_title("Distribution") # Bar chart - y_pos = np.arange(len(labels)) - bars = ax_bar.barh(y_pos, values, color=colors, edgecolor="black", linewidth=1) - ax_bar.set_yticks(y_pos) - ax_bar.set_yticklabels(labels) - ax_bar.set_xlabel("Mass (kg)") - ax_bar.set_title("Component Masses", fontsize=12, fontweight="bold") - ax_bar.grid(True, alpha=0.3, axis="x") + bars = ax2.barh(labels, values, color=colors) + ax2.set_xlabel("Mass (kg)") + ax2.set_title("Component Masses") # Add value labels on bars for bar, val in zip(bars, values, strict=True): - width = bar.get_width() - ax_bar.text(width + total * 0.01, bar.get_y() + bar.get_height() / 2, - f"{val:,.0f} kg", va="center", fontsize=9) + ax2.text( + val + total * 0.01, + bar.get_y() + bar.get_height() / 2, + f"{val:.0f} kg", + va="center", + fontsize=9, + ) + + ax2.set_xlim(0, max(values) * 1.2) + + fig.suptitle(f"{title} (Total: {total:,.0f} kg)", fontsize=14, fontweight="bold") + fig.tight_layout() + + return fig + + +# ============================================================================= +# Cycle Comparison Plots +# ============================================================================= + + +@beartype +def plot_cycle_comparison_bars( + cycle_data: list[dict], + metrics: list[str] | None = None, + figsize: tuple[float, float] = (14, 8), + title: str = "Engine Cycle Comparison", +) -> Figure: + """Create multi-metric bar chart comparing engine cycles. + + Args: + cycle_data: List of dicts with keys 'name' and metric values. + Example: [ + {'name': 'Pressure-Fed', 'net_isp': 320, 'efficiency': 1.0, ...}, + {'name': 'Gas Generator', 'net_isp': 310, 'efficiency': 0.95, ...}, + ] + metrics: List of metrics to plot. Defaults to common cycle metrics. + figsize: Figure size + title: Plot title + + Returns: + matplotlib Figure + """ + _setup_style() + + if metrics is None: + metrics = ["net_isp", "efficiency", "tank_pressure_MPa", "pump_power_kW"] + + # Filter to only metrics that exist in the data + available_metrics = [] + for m in metrics: + if all(m in d for d in cycle_data): + available_metrics.append(m) + + if not available_metrics: + raise ValueError("No valid metrics found in cycle_data") + + n_cycles = len(cycle_data) + n_metrics = len(available_metrics) + + # Metric display names and units + metric_info = { + "net_isp": ("Net Isp", "s"), + "efficiency": ("Cycle Efficiency", "%"), + "tank_pressure_MPa": ("Tank Pressure", "MPa"), + "pump_power_kW": ("Pump Power", "kW"), + "turbine_power_kW": ("Turbine Power", "kW"), + } + + fig, axes = plt.subplots(1, n_metrics, figsize=figsize) + if n_metrics == 1: + axes = [axes] + + cycle_names = [d["name"] for d in cycle_data] + colors = [COLORS["primary"], COLORS["secondary"], COLORS["accent"], "#48A9A6", "#4B3F72"] + + for ax, metric in zip(axes, available_metrics, strict=True): + values = [d.get(metric, 0) for d in cycle_data] + + # Scale efficiency to percentage + if metric == "efficiency": + values = [v * 100 for v in values] + + bars = ax.bar(cycle_names, values, color=colors[:n_cycles], edgecolor="white", linewidth=1.5) + + # Add value labels on bars + for bar, val in zip(bars, values, strict=True): + height = bar.get_height() + ax.annotate( + f"{val:.1f}", + xy=(bar.get_x() + bar.get_width() / 2, height), + xytext=(0, 3), + textcoords="offset points", + ha="center", + va="bottom", + fontsize=10, + fontweight="bold", + ) + + # Axis formatting + display_name, unit = metric_info.get(metric, (metric, "")) + ax.set_ylabel(f"{display_name} ({unit})" if unit else display_name) + ax.set_title(display_name, fontsize=12, fontweight="bold") + ax.tick_params(axis="x", rotation=15) + ax.set_ylim(0, max(values) * 1.2 if max(values) > 0 else 1) + ax.grid(axis="y", alpha=0.3) + + fig.suptitle(title, fontsize=16, fontweight="bold", y=1.02) + fig.tight_layout() + + return fig + + +@beartype +def plot_cycle_radar( + cycle_data: list[dict], + metrics: list[str] | None = None, + figsize: tuple[float, float] = (10, 10), + title: str = "Cycle Comparison Radar", +) -> Figure: + """Create radar/spider chart comparing cycles across normalized dimensions. + + All metrics are normalized to 0-1 scale for comparison. + + Args: + cycle_data: List of dicts with 'name' and metric values + metrics: Metrics to include in radar. Defaults to standard set. + figsize: Figure size + title: Plot title + + Returns: + matplotlib Figure + """ + _setup_style() - ax_bar.set_xlim(0, max(values) * 1.15) + if metrics is None: + metrics = ["net_isp", "efficiency", "simplicity", "tank_mass_factor"] + + # Filter to available metrics + available_metrics = [] + for m in metrics: + if all(m in d for d in cycle_data): + available_metrics.append(m) + + if len(available_metrics) < 3: + raise ValueError("Need at least 3 metrics for radar chart") + + n_metrics = len(available_metrics) + + # Metric display names + metric_names = { + "net_isp": "Performance\n(Isp)", + "efficiency": "Efficiency", + "simplicity": "Simplicity", + "tank_mass_factor": "Low Tank\nPressure", + "reliability": "Reliability", + "trl": "Maturity\n(TRL)", + } + + # Normalize all metrics to 0-1 scale + normalized_data = [] + for d in cycle_data: + norm_d = {"name": d["name"]} + for m in available_metrics: + values = [cd[m] for cd in cycle_data] + min_val, max_val = min(values), max(values) + if max_val > min_val: + norm_d[m] = (d[m] - min_val) / (max_val - min_val) + else: + norm_d[m] = 1.0 + normalized_data.append(norm_d) + + # Setup radar chart + angles = np.linspace(0, 2 * np.pi, n_metrics, endpoint=False).tolist() + angles += angles[:1] # Complete the loop + + fig, ax = plt.subplots(figsize=figsize, subplot_kw=dict(polar=True)) + + colors = [COLORS["primary"], COLORS["secondary"], COLORS["accent"], "#48A9A6", "#4B3F72"] + + for i, d in enumerate(normalized_data): + values = [d[m] for m in available_metrics] + values += values[:1] # Complete the loop + + ax.plot(angles, values, "o-", linewidth=2, color=colors[i % len(colors)], label=d["name"]) + ax.fill(angles, values, alpha=0.25, color=colors[i % len(colors)]) + + # Set labels + labels = [metric_names.get(m, m) for m in available_metrics] + ax.set_xticks(angles[:-1]) + ax.set_xticklabels(labels, fontsize=11) + + # Radial limits + ax.set_ylim(0, 1.1) + ax.set_yticks([0.25, 0.5, 0.75, 1.0]) + ax.set_yticklabels(["25%", "50%", "75%", "100%"], fontsize=8, alpha=0.7) + + ax.legend(loc="upper right", bbox_to_anchor=(1.3, 1.0), fontsize=11) + ax.set_title(title, fontsize=14, fontweight="bold", y=1.08) - # Total mass annotation - fig.suptitle(f"{title}\nTotal: {total:,.0f} kg", fontsize=14, fontweight="bold") - fig.tight_layout(rect=[0, 0, 1, 0.93]) + fig.tight_layout() return fig + +@beartype +def plot_cycle_tradeoff( + cycle_data: list[dict], + x_metric: str = "net_isp", + y_metric: str = "efficiency", + size_metric: str | None = None, + figsize: tuple[float, float] = (10, 8), + title: str = "Cycle Trade Space", +) -> Figure: + """Create scatter plot showing cycle trade-offs. + + Plot cycles on 2D trade space with optional bubble size for third dimension. + + Args: + cycle_data: List of dicts with 'name' and metric values + x_metric: Metric for x-axis + y_metric: Metric for y-axis + size_metric: Optional metric for bubble size + figsize: Figure size + title: Plot title + + Returns: + matplotlib Figure + """ + _setup_style() + + # Metric display info + metric_info = { + "net_isp": ("Net Isp", "s"), + "efficiency": ("Cycle Efficiency", ""), + "tank_pressure_MPa": ("Tank Pressure", "MPa"), + "pump_power_kW": ("Pump Power", "kW"), + "simplicity": ("Simplicity Score", ""), + "complexity": ("Complexity Score", ""), + } + + fig, ax = plt.subplots(figsize=figsize) + + x_vals = [d[x_metric] for d in cycle_data] + y_vals = [d[y_metric] for d in cycle_data] + + # Scale efficiency to percentage for display + if y_metric == "efficiency": + y_vals = [v * 100 for v in y_vals] + + colors = [COLORS["primary"], COLORS["secondary"], COLORS["accent"], "#48A9A6", "#4B3F72"] + + if size_metric and all(size_metric in d for d in cycle_data): + sizes = [d[size_metric] for d in cycle_data] + # Normalize sizes for display + max_size = max(sizes) if max(sizes) > 0 else 1 + normalized_sizes = [500 + 1500 * (s / max_size) for s in sizes] + else: + normalized_sizes = [800] * len(cycle_data) + + for i, (x, y, s, d) in enumerate(zip(x_vals, y_vals, normalized_sizes, cycle_data, strict=True)): + ax.scatter( + x, y, s=s, c=colors[i % len(colors)], alpha=0.7, edgecolors="white", linewidth=2, zorder=5 + ) + # Label + ax.annotate( + d["name"], + xy=(x, y), + xytext=(10, 10), + textcoords="offset points", + fontsize=11, + fontweight="bold", + arrowprops=dict(arrowstyle="-", color="gray", alpha=0.5), + ) + + # Axis formatting + x_name, x_unit = metric_info.get(x_metric, (x_metric, "")) + y_name, y_unit = metric_info.get(y_metric, (y_metric, "")) + + ax.set_xlabel(f"{x_name} ({x_unit})" if x_unit else x_name, fontsize=12) + ax.set_ylabel(f"{y_name} ({y_unit})" if y_unit else y_name, fontsize=12) + + # Add quadrant annotations + x_mid = (max(x_vals) + min(x_vals)) / 2 + y_mid = (max(y_vals) + min(y_vals)) / 2 + + ax.axvline(x=x_mid, color="gray", linestyle="--", alpha=0.3) + ax.axhline(y=y_mid, color="gray", linestyle="--", alpha=0.3) + + # Expand limits slightly + x_range = max(x_vals) - min(x_vals) + y_range = max(y_vals) - min(y_vals) + ax.set_xlim(min(x_vals) - 0.1 * x_range, max(x_vals) + 0.15 * x_range) + ax.set_ylim(min(y_vals) - 0.1 * y_range, max(y_vals) + 0.15 * y_range) + + ax.grid(True, alpha=0.3) + ax.set_title(title, fontsize=14, fontweight="bold") + + # Add size legend if applicable + if size_metric and all(size_metric in d for d in cycle_data): + size_name = metric_info.get(size_metric, (size_metric, ""))[0] + ax.text( + 0.02, + 0.02, + f"Bubble size: {size_name}", + transform=ax.transAxes, + fontsize=9, + alpha=0.7, + ) + + fig.tight_layout() + + return fig diff --git a/rocket/results.py b/rocket/results.py new file mode 100644 index 0000000..9ebb6b5 --- /dev/null +++ b/rocket/results.py @@ -0,0 +1,658 @@ +"""Results visualization and Pareto front analysis for Rocket. + +This module provides plotting utilities for parametric study and uncertainty +analysis results. Integrates with the analysis module for seamless visualization. + +Example: + >>> from rocket.analysis import ParametricStudy, Range + >>> from rocket.results import plot_1d, plot_2d_contour, plot_pareto + >>> + >>> results = study.run() + >>> fig = plot_1d(results, "chamber_pressure", "isp_vac") + >>> fig = plot_pareto(results, "isp_vac", "thrust_to_weight", maximize=[True, True]) +""" + +from collections.abc import Sequence + +import matplotlib.pyplot as plt +import numpy as np +from beartype import beartype +from matplotlib.figure import Figure +from numpy.typing import NDArray + +from rocket.analysis import StudyResults, UncertaintyResults + +# ============================================================================= +# Plot Style Configuration +# ============================================================================= + +COLORS = { + "primary": "#2E86AB", + "secondary": "#A23B72", + "accent": "#F18F01", + "feasible": "#2E86AB", + "infeasible": "#CCCCCC", + "pareto": "#E63946", + "grid": "#CCCCCC", + "text": "#333333", +} + + +def _setup_style() -> None: + """Configure matplotlib style for consistent appearance.""" + plt.rcParams.update({ + "font.family": "sans-serif", + "font.sans-serif": ["Helvetica", "Arial", "DejaVu Sans"], + "font.size": 11, + "axes.titlesize": 14, + "axes.labelsize": 12, + "axes.linewidth": 1.2, + "xtick.labelsize": 10, + "ytick.labelsize": 10, + "legend.fontsize": 10, + "figure.titlesize": 16, + "grid.alpha": 0.5, + }) + + +# ============================================================================= +# 1D Parameter Sweep Plots +# ============================================================================= + + +@beartype +def plot_1d( + results: StudyResults, + x_param: str, + y_metric: str, + show_infeasible: bool = True, + figsize: tuple[float, float] = (10, 6), + title: str | None = None, + xlabel: str | None = None, + ylabel: str | None = None, +) -> Figure: + """Plot a 1D parameter sweep. + + Args: + results: StudyResults from ParametricStudy + x_param: Parameter name for x-axis + y_metric: Metric name for y-axis + show_infeasible: Whether to show infeasible points (grayed out) + figsize: Figure size (width, height) + title: Optional plot title + xlabel: Optional x-axis label + ylabel: Optional y-axis label + + Returns: + matplotlib Figure + """ + _setup_style() + + fig, ax = plt.subplots(figsize=figsize) + + x = results.get_parameter(x_param) + y = results.get_metric(y_metric) + + if results.constraints_passed is not None and show_infeasible: + # Plot infeasible points first (gray) + infeasible = ~results.constraints_passed + if np.any(infeasible): + ax.scatter( + x[infeasible], y[infeasible], + c=COLORS["infeasible"], s=50, alpha=0.5, + label="Infeasible", zorder=1, + ) + + # Plot feasible points + feasible = results.constraints_passed + if np.any(feasible): + ax.scatter( + x[feasible], y[feasible], + c=COLORS["primary"], s=80, + label="Feasible", zorder=2, + ) + ax.plot( + x[feasible], y[feasible], + c=COLORS["primary"], alpha=0.5, linewidth=1.5, zorder=1, + ) + else: + ax.scatter(x, y, c=COLORS["primary"], s=80) + ax.plot(x, y, c=COLORS["primary"], alpha=0.5, linewidth=1.5) + + # Mark optimum + if results.constraints_passed is not None: + feasible_mask = results.constraints_passed + if np.any(feasible_mask): + best_idx = np.argmax(y * feasible_mask.astype(float) - (~feasible_mask) * 1e10) + ax.scatter( + [x[best_idx]], [y[best_idx]], + c=COLORS["accent"], s=150, marker="*", + label=f"Best: {y[best_idx]:.3g}", zorder=3, + ) + + ax.set_xlabel(xlabel or x_param) + ax.set_ylabel(ylabel or y_metric) + ax.set_title(title or f"{y_metric} vs {x_param}") + ax.grid(True, alpha=0.3) + ax.legend() + + fig.tight_layout() + return fig + + +@beartype +def plot_1d_multi( + results: StudyResults, + x_param: str, + y_metrics: list[str], + figsize: tuple[float, float] = (12, 6), + title: str | None = None, + xlabel: str | None = None, + normalize: bool = False, +) -> Figure: + """Plot multiple metrics vs a single parameter. + + Args: + results: StudyResults from ParametricStudy + x_param: Parameter name for x-axis + y_metrics: List of metric names to plot + figsize: Figure size + title: Optional title + xlabel: Optional x-axis label + normalize: If True, normalize all metrics to [0, 1] + + Returns: + matplotlib Figure + """ + _setup_style() + + fig, ax = plt.subplots(figsize=figsize) + + x = results.get_parameter(x_param) + colors = plt.cm.tab10(np.linspace(0, 1, len(y_metrics))) + + for i, metric in enumerate(y_metrics): + y = results.get_metric(metric) + if normalize: + y = (y - np.nanmin(y)) / (np.nanmax(y) - np.nanmin(y) + 1e-10) + + ax.plot(x, y, color=colors[i], linewidth=2, label=metric, marker="o", markersize=4) + + ax.set_xlabel(xlabel or x_param) + ax.set_ylabel("Normalized Value" if normalize else "Metric Value") + ax.set_title(title or f"Metrics vs {x_param}") + ax.grid(True, alpha=0.3) + ax.legend(loc="best") + + fig.tight_layout() + return fig + + +# ============================================================================= +# 2D Contour Plots +# ============================================================================= + + +@beartype +def plot_2d_contour( + results: StudyResults, + x_param: str, + y_param: str, + z_metric: str, + figsize: tuple[float, float] = (10, 8), + title: str | None = None, + levels: int = 20, + show_points: bool = True, + show_infeasible: bool = True, +) -> Figure: + """Plot a 2D contour for two-parameter sweeps. + + Args: + results: StudyResults from ParametricStudy (must be 2D grid) + x_param: Parameter for x-axis + y_param: Parameter for y-axis + z_metric: Metric for contour values + figsize: Figure size + title: Optional title + levels: Number of contour levels + show_points: Show scatter points at evaluated locations + show_infeasible: Mark infeasible regions + + Returns: + matplotlib Figure + """ + _setup_style() + + x = results.get_parameter(x_param) + y = results.get_parameter(y_param) + z = results.get_metric(z_metric) + + # Get unique values to determine grid shape + x_unique = np.unique(x) + y_unique = np.unique(y) + + if len(x_unique) * len(y_unique) != len(x): + raise ValueError( + f"Data is not a complete 2D grid. Got {len(x)} points, " + f"expected {len(x_unique)} x {len(y_unique)} = {len(x_unique) * len(y_unique)}" + ) + + # Reshape to 2D grid + nx, ny = len(x_unique), len(y_unique) + X = x.reshape(ny, nx) + Y = y.reshape(ny, nx) + Z = z.reshape(ny, nx) + + fig, ax = plt.subplots(figsize=figsize) + + # Contour plot + contour = ax.contourf(X, Y, Z, levels=levels, cmap="viridis") + ax.contour(X, Y, Z, levels=levels, colors="white", alpha=0.3, linewidths=0.5) + + # Colorbar + _ = fig.colorbar(contour, ax=ax, label=z_metric) + + # Show evaluated points + if show_points: + ax.scatter(x, y, c="white", s=10, alpha=0.5, edgecolors="black", linewidths=0.5) + + # Mark infeasible regions + if show_infeasible and results.constraints_passed is not None: + infeasible = ~results.constraints_passed + if np.any(infeasible): + ax.scatter( + x[infeasible], y[infeasible], + c="red", s=30, marker="x", alpha=0.7, + label="Infeasible", + ) + ax.legend() + + ax.set_xlabel(x_param) + ax.set_ylabel(y_param) + ax.set_title(title or f"{z_metric}") + + fig.tight_layout() + return fig + + +# ============================================================================= +# Pareto Front Analysis +# ============================================================================= + + +def _compute_pareto_front( + objectives: NDArray[np.float64], + maximize: Sequence[bool], +) -> NDArray[np.bool_]: + """Compute Pareto-optimal points. + + Args: + objectives: Array of shape (n_points, n_objectives) + maximize: List of booleans indicating whether to maximize each objective + + Returns: + Boolean array indicating Pareto-optimal points + """ + n_points = objectives.shape[0] + is_pareto = np.ones(n_points, dtype=bool) + + # Flip signs for maximization objectives + obj = objectives.copy() + for i, is_max in enumerate(maximize): + if is_max: + obj[:, i] = -obj[:, i] + + for i in range(n_points): + if not is_pareto[i]: + continue + + # Check if any other point dominates this one + for j in range(n_points): + if i == j or not is_pareto[j]: + continue + + # j dominates i if j is <= in all objectives and < in at least one + dominates = np.all(obj[j] <= obj[i]) and np.any(obj[j] < obj[i]) + if dominates: + is_pareto[i] = False + break + + return is_pareto + + +@beartype +def plot_pareto( + results: StudyResults, + x_metric: str, + y_metric: str, + maximize: tuple[bool, bool] = (True, True), + feasible_only: bool = True, + figsize: tuple[float, float] = (10, 8), + title: str | None = None, + show_dominated: bool = True, +) -> Figure: + """Plot Pareto front for two objectives. + + Args: + results: StudyResults from ParametricStudy + x_metric: First objective (x-axis) + y_metric: Second objective (y-axis) + maximize: Tuple of booleans for each objective (True = maximize) + feasible_only: Only consider feasible points + figsize: Figure size + title: Optional title + show_dominated: Show dominated points in gray + + Returns: + matplotlib Figure + """ + _setup_style() + + # Get metrics + x = results.get_metric(x_metric) + y = results.get_metric(y_metric) + + # Filter to feasible if requested + if feasible_only and results.constraints_passed is not None: + mask = results.constraints_passed + else: + mask = np.ones(len(x), dtype=bool) + + # Remove NaN values + valid = mask & ~np.isnan(x) & ~np.isnan(y) + x_valid = x[valid] + y_valid = y[valid] + + if len(x_valid) == 0: + raise ValueError("No valid data points for Pareto analysis") + + # Compute Pareto front + objectives = np.column_stack([x_valid, y_valid]) + is_pareto = _compute_pareto_front(objectives, maximize) + + fig, ax = plt.subplots(figsize=figsize) + + # Plot dominated points + if show_dominated: + dominated = ~is_pareto + ax.scatter( + x_valid[dominated], y_valid[dominated], + c=COLORS["infeasible"], s=50, alpha=0.5, + label="Dominated", + ) + + # Plot Pareto front points + ax.scatter( + x_valid[is_pareto], y_valid[is_pareto], + c=COLORS["pareto"], s=100, + label=f"Pareto Front ({np.sum(is_pareto)} points)", + zorder=3, + ) + + # Connect Pareto points with line (sorted) + pareto_x = x_valid[is_pareto] + pareto_y = y_valid[is_pareto] + sort_idx = np.argsort(pareto_x) + ax.plot( + pareto_x[sort_idx], pareto_y[sort_idx], + c=COLORS["pareto"], linewidth=2, alpha=0.7, + linestyle="--", + ) + + # Add direction arrows + x_dir = "→" if maximize[0] else "←" + y_dir = "↑" if maximize[1] else "↓" + ax.annotate( + f"Better {x_dir}", + xy=(0.95, 0.02), xycoords="axes fraction", + ha="right", fontsize=10, color=COLORS["text"], + ) + ax.annotate( + f"Better {y_dir}", + xy=(0.02, 0.95), xycoords="axes fraction", + ha="left", fontsize=10, color=COLORS["text"], rotation=90, + ) + + ax.set_xlabel(x_metric) + ax.set_ylabel(y_metric) + ax.set_title(title or f"Pareto Front: {x_metric} vs {y_metric}") + ax.grid(True, alpha=0.3) + ax.legend() + + fig.tight_layout() + return fig + + +@beartype +def get_pareto_points( + results: StudyResults, + metrics: list[str], + maximize: list[bool], + feasible_only: bool = True, +) -> tuple[list[int], NDArray[np.float64]]: + """Get indices and values of Pareto-optimal points. + + Args: + results: StudyResults from ParametricStudy + metrics: List of objective metric names + maximize: List of booleans for each metric + feasible_only: Only consider feasible points + + Returns: + Tuple of (indices in original results, objective values array) + """ + # Get metrics + obj_arrays = [results.get_metric(m) for m in metrics] + objectives = np.column_stack(obj_arrays) + + # Filter to feasible + if feasible_only and results.constraints_passed is not None: + mask = results.constraints_passed + else: + mask = np.ones(objectives.shape[0], dtype=bool) + + # Remove NaN + valid = mask & ~np.any(np.isnan(objectives), axis=1) + valid_indices = np.where(valid)[0] + objectives_valid = objectives[valid] + + # Compute Pareto front + is_pareto = _compute_pareto_front(objectives_valid, maximize) + + pareto_indices = valid_indices[is_pareto].tolist() + pareto_values = objectives_valid[is_pareto] + + return pareto_indices, pareto_values + + +# ============================================================================= +# Uncertainty Visualization +# ============================================================================= + + +@beartype +def plot_histogram( + results: UncertaintyResults, + metric: str, + bins: int = 50, + show_ci: bool = True, + ci_level: float = 0.95, + figsize: tuple[float, float] = (10, 6), + title: str | None = None, + feasible_only: bool = False, +) -> Figure: + """Plot histogram of a metric from uncertainty analysis. + + Args: + results: UncertaintyResults from UncertaintyAnalysis + metric: Metric name to plot + bins: Number of histogram bins + show_ci: Show confidence interval lines + ci_level: Confidence level for CI lines + figsize: Figure size + title: Optional title + feasible_only: Only include feasible samples + + Returns: + matplotlib Figure + """ + _setup_style() + + values = results.metrics[metric] + if feasible_only: + values = values[results.constraints_passed] + + # Remove NaN + values = values[~np.isnan(values)] + + fig, ax = plt.subplots(figsize=figsize) + + # Histogram + ax.hist(values, bins=bins, color=COLORS["primary"], alpha=0.7, edgecolor="white") + + # Statistics + mean = np.mean(values) + std = np.std(values) + + ax.axvline(mean, color=COLORS["accent"], linewidth=2, label=f"Mean: {mean:.4g}") + + if show_ci: + ci = results.confidence_interval(metric, ci_level, feasible_only) + ax.axvline(ci[0], color=COLORS["secondary"], linewidth=1.5, linestyle="--", + label=f"{ci_level*100:.0f}% CI: [{ci[0]:.4g}, {ci[1]:.4g}]") + ax.axvline(ci[1], color=COLORS["secondary"], linewidth=1.5, linestyle="--") + + ax.set_xlabel(metric) + ax.set_ylabel("Frequency") + ax.set_title(title or f"Distribution of {metric}") + ax.legend() + ax.grid(True, alpha=0.3, axis="y") + + # Add text box with statistics + stats_text = f"μ = {mean:.4g}\nσ = {std:.4g}\nn = {len(values)}" + ax.text( + 0.95, 0.95, stats_text, + transform=ax.transAxes, ha="right", va="top", + fontsize=10, fontfamily="monospace", + bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), + ) + + fig.tight_layout() + return fig + + +@beartype +def plot_correlation( + results: UncertaintyResults, + x_param: str, + y_metric: str, + figsize: tuple[float, float] = (10, 6), + title: str | None = None, +) -> Figure: + """Plot correlation between input parameter and output metric. + + Args: + results: UncertaintyResults from UncertaintyAnalysis + x_param: Input parameter name (from samples) + y_metric: Output metric name + figsize: Figure size + title: Optional title + + Returns: + matplotlib Figure + """ + _setup_style() + + if x_param not in results.samples: + raise ValueError(f"Parameter '{x_param}' not in samples. Available: {list(results.samples.keys())}") + + x = results.samples[x_param] + y = results.metrics[y_metric] + + # Remove NaN pairs + valid = ~(np.isnan(x) | np.isnan(y)) + x = x[valid] + y = y[valid] + + fig, ax = plt.subplots(figsize=figsize) + + ax.scatter(x, y, c=COLORS["primary"], alpha=0.3, s=20) + + # Compute correlation + corr = np.corrcoef(x, y)[0, 1] + + # Add trend line + z = np.polyfit(x, y, 1) + p = np.poly1d(z) + x_line = np.linspace(x.min(), x.max(), 100) + ax.plot(x_line, p(x_line), c=COLORS["accent"], linewidth=2, + label=f"Trend (r = {corr:.3f})") + + ax.set_xlabel(x_param) + ax.set_ylabel(y_metric) + ax.set_title(title or f"Correlation: {x_param} vs {y_metric}") + ax.legend() + ax.grid(True, alpha=0.3) + + fig.tight_layout() + return fig + + +@beartype +def plot_tornado( + results: UncertaintyResults, + metric: str, + parameters: list[str] | None = None, + figsize: tuple[float, float] = (10, 6), + title: str | None = None, +) -> Figure: + """Plot tornado chart showing sensitivity of metric to input parameters. + + Shows which parameters have the strongest influence on the metric. + + Args: + results: UncertaintyResults from UncertaintyAnalysis + metric: Metric to analyze + parameters: List of parameters to include (None = all) + figsize: Figure size + title: Optional title + + Returns: + matplotlib Figure + """ + _setup_style() + + if parameters is None: + parameters = list(results.samples.keys()) + + y_values = results.metrics[metric] + correlations: dict[str, float] = {} + + for param in parameters: + x_values = results.samples[param] + valid = ~(np.isnan(x_values) | np.isnan(y_values)) + if np.sum(valid) > 2: + corr = np.corrcoef(x_values[valid], y_values[valid])[0, 1] + correlations[param] = corr + + # Sort by absolute correlation + sorted_params = sorted(correlations.keys(), key=lambda p: abs(correlations[p])) + sorted_corrs = [correlations[p] for p in sorted_params] + + fig, ax = plt.subplots(figsize=figsize) + + y_pos = np.arange(len(sorted_params)) + colors = [COLORS["primary"] if c >= 0 else COLORS["secondary"] for c in sorted_corrs] + + ax.barh(y_pos, sorted_corrs, color=colors, alpha=0.8, edgecolor="white") + ax.set_yticks(y_pos) + ax.set_yticklabels(sorted_params) + ax.set_xlabel("Correlation Coefficient") + ax.set_title(title or f"Sensitivity of {metric}") + ax.axvline(0, color="black", linewidth=0.5) + ax.set_xlim(-1, 1) + ax.grid(True, alpha=0.3, axis="x") + + fig.tight_layout() + return fig + diff --git a/rocket/system.py b/rocket/system.py new file mode 100644 index 0000000..58c7f6e --- /dev/null +++ b/rocket/system.py @@ -0,0 +1,244 @@ +"""High-level system design API for Rocket. + +This module provides the main entry point for complete engine system design, +integrating: +- Engine performance and geometry +- Cycle analysis (pressure-fed, gas-generator, etc.) +- Thermal/cooling feasibility + +Example: + >>> from rocket import EngineInputs + >>> from rocket.system import design_engine_system + >>> from rocket.cycles import GasGeneratorCycle + >>> from rocket.units import kilonewtons, megapascals, kelvin + >>> + >>> inputs = EngineInputs.from_propellants( + ... oxidizer="LOX", + ... fuel="CH4", + ... thrust=kilonewtons(2000), + ... chamber_pressure=megapascals(30), + ... ) + >>> + >>> result = design_engine_system( + ... inputs=inputs, + ... cycle=GasGeneratorCycle(turbine_inlet_temp=kelvin(900)), + ... check_cooling=True, + ... ) + >>> + >>> print(f"Net Isp: {result.cycle.net_isp.value:.1f} s") + >>> print(f"Cooling feasible: {result.cooling.feasible}") +""" + +from dataclasses import dataclass +from typing import Any + +from beartype import beartype + +from rocket.cycles.base import CycleConfiguration, CyclePerformance, analyze_cycle +from rocket.engine import ( + EngineGeometry, + EngineInputs, + EnginePerformance, + compute_geometry, + compute_performance, +) +from rocket.thermal.regenerative import CoolingFeasibility, check_cooling_feasibility + + +@beartype +@dataclass(frozen=True) +class EngineSystemResult: + """Complete engine system design result. + + Contains all analysis results from engine performance through + cycle analysis and thermal assessment. + + Attributes: + inputs: Original engine inputs + performance: Ideal engine performance (before cycle losses) + geometry: Engine geometry + cycle: Cycle analysis results (if cycle provided) + cooling: Cooling feasibility results (if check_cooling=True) + feasible: Overall feasibility (cycle closes AND cooling ok) + warnings: Combined warnings from all analyses + """ + + inputs: EngineInputs + performance: EnginePerformance + geometry: EngineGeometry + cycle: CyclePerformance | None + cooling: CoolingFeasibility | None + feasible: bool + warnings: list[str] + + +@beartype +def design_engine_system( + inputs: EngineInputs, + cycle: CycleConfiguration | None = None, + check_cooling: bool = True, + coolant: str | None = None, + max_wall_temp: Any | None = None, +) -> EngineSystemResult: + """Design a complete engine system with cycle and thermal analysis. + + This is the main entry point for rocket engine system design. It: + 1. Computes ideal engine performance and geometry + 2. Analyzes the engine cycle (if specified) + 3. Checks cooling feasibility (if requested) + 4. Returns a comprehensive result with all analyses + + Args: + inputs: Engine input parameters (from EngineInputs.from_propellants or manual) + cycle: Engine cycle configuration (PressureFedCycle, GasGeneratorCycle, etc.) + If None, only basic performance is computed. + check_cooling: Whether to perform cooling feasibility analysis + coolant: Coolant name for cooling analysis. If None, uses fuel. + max_wall_temp: Maximum allowed wall temperature [K]. If None, uses defaults. + + Returns: + EngineSystemResult with all analysis results + + Example: + >>> # Basic design without cycle analysis + >>> result = design_engine_system(inputs) + >>> print(f"Isp: {result.performance.isp.value:.1f} s") + >>> + >>> # Full system design with gas generator cycle + >>> from rocket.cycles import GasGeneratorCycle + >>> result = design_engine_system( + ... inputs=inputs, + ... cycle=GasGeneratorCycle(turbine_inlet_temp=kelvin(900)), + ... check_cooling=True, + ... ) + """ + all_warnings: list[str] = [] + + # Step 1: Compute basic engine performance and geometry + performance = compute_performance(inputs) + geometry = compute_geometry(inputs, performance) + + # Step 2: Cycle analysis (if cycle provided) + cycle_result: CyclePerformance | None = None + if cycle is not None: + cycle_result = analyze_cycle(inputs, performance, geometry, cycle) + all_warnings.extend(cycle_result.warnings) + + # Step 3: Cooling feasibility (if requested) + cooling_result: CoolingFeasibility | None = None + if check_cooling: + # Determine coolant (default to fuel) + if coolant is None: + # Try to infer fuel from engine name + coolant = _infer_coolant(inputs) + + cooling_result = check_cooling_feasibility( + inputs=inputs, + performance=performance, + geometry=geometry, + coolant=coolant, + max_wall_temp=max_wall_temp, + ) + all_warnings.extend(cooling_result.warnings) + + # Determine overall feasibility + feasible = True + if cycle_result is not None and not cycle_result.feasible: + feasible = False + if cooling_result is not None and not cooling_result.feasible: + feasible = False + + return EngineSystemResult( + inputs=inputs, + performance=performance, + geometry=geometry, + cycle=cycle_result, + cooling=cooling_result, + feasible=feasible, + warnings=all_warnings, + ) + + +def _infer_coolant(inputs: EngineInputs) -> str: + """Infer coolant from engine inputs.""" + if inputs.name: + name_upper = inputs.name.upper() + if "CH4" in name_upper or "METHANE" in name_upper or "METHALOX" in name_upper: + return "CH4" + elif "LH2" in name_upper or "HYDROGEN" in name_upper or "HYDROLOX" in name_upper: + return "LH2" + elif "RP1" in name_upper or "KEROSENE" in name_upper or "KEROLOX" in name_upper: + return "RP1" + elif "ETHANOL" in name_upper: + return "Ethanol" + + # Default to RP-1 + return "RP1" + + +@beartype +def format_system_summary(result: EngineSystemResult) -> str: + """Format complete system design results as readable string. + + Args: + result: EngineSystemResult from design_engine_system() + + Returns: + Formatted multi-line string summary + """ + name = result.inputs.name or "Unnamed Engine" + status = "✓ FEASIBLE" if result.feasible else "✗ NOT FEASIBLE" + + lines = [ + "=" * 70, + f"ENGINE SYSTEM DESIGN: {name}", + f"Overall Status: {status}", + "=" * 70, + "", + "PERFORMANCE (Ideal):", + f" Thrust (SL): {result.inputs.thrust.to('kN').value:.1f} kN", + f" Isp (SL): {result.performance.isp.value:.1f} s", + f" Isp (Vac): {result.performance.isp_vac.value:.1f} s", + f" C*: {result.performance.cstar.value:.0f} m/s", + f" Mass Flow: {result.performance.mdot.value:.2f} kg/s", + "", + "GEOMETRY:", + f" Throat Dia: {result.geometry.throat_diameter.to('m').value*100:.1f} cm", + f" Exit Dia: {result.geometry.exit_diameter.to('m').value*100:.1f} cm", + f" Chamber Dia: {result.geometry.chamber_diameter.to('m').value*100:.1f} cm", + f" Expansion Ratio: {result.geometry.expansion_ratio:.1f}", + ] + + if result.cycle is not None: + lines.extend([ + "", + f"CYCLE ({result.cycle.cycle_type.name}):", + f" Net Isp: {result.cycle.net_isp.value:.1f} s", + f" Cycle Efficiency:{result.cycle.cycle_efficiency*100:.1f}%", + f" Turbine Power: {result.cycle.turbine_power.value/1000:.0f} kW", + f" GG/Turbine Flow: {result.cycle.turbine_mass_flow.value:.2f} kg/s", + f" Tank P (Ox): {result.cycle.tank_pressure_ox.to('bar').value:.1f} bar", + f" Tank P (Fuel): {result.cycle.tank_pressure_fuel.to('bar').value:.1f} bar", + ]) + + if result.cooling is not None: + lines.extend([ + "", + "COOLING:", + f" Throat Heat Flux:{result.cooling.throat_heat_flux.value/1e6:.1f} MW/m²", + f" Max Wall Temp: {result.cooling.max_wall_temp.value:.0f} K", + f" Allowed Temp: {result.cooling.max_allowed_temp.value:.0f} K", + f" Coolant ΔT: {result.cooling.coolant_temp_rise.value:.0f} K", + f" Flow Margin: {result.cooling.flow_margin:.2f}x", + f" Pressure Drop: {result.cooling.pressure_drop.to('bar').value:.1f} bar", + ]) + + if result.warnings: + lines.extend(["", "WARNINGS:"]) + for warning in result.warnings: + lines.append(f" ⚠ {warning}") + + lines.append("=" * 70) + + return "\n".join(lines) + diff --git a/rocket/thermal/__init__.py b/rocket/thermal/__init__.py new file mode 100644 index 0000000..86b1ea5 --- /dev/null +++ b/rocket/thermal/__init__.py @@ -0,0 +1,57 @@ +"""Thermal analysis module for Rocket. + +This module provides tools for analyzing thermal loads on rocket engine +components and evaluating cooling system feasibility. + +Key capabilities: +- Heat flux estimation using Bartz correlation +- Regenerative cooling feasibility screening +- Wall temperature prediction +- Coolant property database + +Example: + >>> from rocket import EngineInputs, design_engine + >>> from rocket.thermal import estimate_heat_flux, check_cooling_feasibility + >>> + >>> inputs = EngineInputs.from_propellants("LOX", "CH4", ...) + >>> performance, geometry = design_engine(inputs) + >>> + >>> q_throat = estimate_heat_flux(inputs, performance, geometry, location="throat") + >>> print(f"Throat heat flux: {q_throat.value/1e6:.1f} MW/m²") + >>> + >>> cooling = check_cooling_feasibility( + ... inputs, performance, geometry, + ... coolant="CH4", + ... max_wall_temp=kelvin(800), + ... ) + >>> print(f"Cooling feasible: {cooling.feasible}") +""" + +from rocket.thermal.heat_flux import ( + adiabatic_wall_temperature, + bartz_heat_flux, + estimate_heat_flux, + heat_flux_profile, + recovery_factor, +) +from rocket.thermal.regenerative import ( + CoolantProperties, + CoolingFeasibility, + check_cooling_feasibility, + get_coolant_properties, +) + +__all__ = [ + # Heat flux estimation + "adiabatic_wall_temperature", + "bartz_heat_flux", + "estimate_heat_flux", + "heat_flux_profile", + "recovery_factor", + # Regenerative cooling + "CoolingFeasibility", + "CoolantProperties", + "check_cooling_feasibility", + "get_coolant_properties", +] + diff --git a/rocket/thermal/heat_flux.py b/rocket/thermal/heat_flux.py new file mode 100644 index 0000000..7b244f5 --- /dev/null +++ b/rocket/thermal/heat_flux.py @@ -0,0 +1,306 @@ +"""Heat flux estimation for rocket engine combustion chambers and nozzles. + +This module implements the Bartz correlation and related methods for +estimating convective heat transfer in rocket engines. + +The Bartz correlation is the industry-standard method for preliminary +heat flux estimation, derived from turbulent pipe flow correlations +modified for rocket engine conditions. + +References: + - Bartz, D.R., "A Simple Equation for Rapid Estimation of Rocket + Nozzle Convective Heat Transfer Coefficients", Jet Propulsion, 1957 + - Huzel & Huang, "Modern Engineering for Design of Liquid-Propellant + Rocket Engines", Chapter 4 +""" + +import math + +import numpy as np +from beartype import beartype + +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.units import Quantity, kelvin + + +@beartype +def recovery_factor(prandtl: float, laminar: bool = False) -> float: + """Calculate the recovery factor for adiabatic wall temperature. + + The recovery factor accounts for the difference between the + stagnation temperature and the actual adiabatic wall temperature + due to boundary layer effects. + + Args: + prandtl: Prandtl number of the gas [-] + laminar: If True, use laminar correlation; else turbulent + + Returns: + Recovery factor r [-], typically 0.85-0.92 for turbulent flow + """ + if laminar: + return math.sqrt(prandtl) # r = Pr^0.5 for laminar + else: + return prandtl ** (1/3) # r = Pr^(1/3) for turbulent + + +@beartype +def adiabatic_wall_temperature( + stagnation_temp: Quantity, + mach: float, + gamma: float, + recovery_factor: float, +) -> Quantity: + """Calculate adiabatic wall temperature. + + The adiabatic wall temperature is the temperature the wall would + reach if there were no heat transfer (perfectly insulated wall). + It's used as the driving temperature difference for heat flux. + + T_aw = T_0 * [1 + r * (gamma-1)/2 * M^2] / [1 + (gamma-1)/2 * M^2] + + Args: + stagnation_temp: Chamber/stagnation temperature [K] + mach: Local Mach number [-] + gamma: Ratio of specific heats [-] + recovery_factor: Recovery factor r [-] + + Returns: + Adiabatic wall temperature [K] + """ + T0 = stagnation_temp.to("K").value + gm1 = gamma - 1 + + # Temperature ratio T/T0 from isentropic relations + T_ratio = 1 / (1 + gm1/2 * mach**2) + + # Static temperature + T_static = T0 * T_ratio + + # Adiabatic wall temperature + # T_aw = T_static + r * (T0 - T_static) + T_aw = T_static + recovery_factor * (T0 - T_static) + + return kelvin(T_aw) + + +@beartype +def bartz_heat_flux( + chamber_pressure: Quantity, + chamber_temp: Quantity, + throat_diameter: Quantity, + local_diameter: Quantity, + characteristic_velocity: Quantity, + gamma: float, + molecular_weight: float, + local_mach: float, + wall_temp: Quantity | None = None, +) -> Quantity: + """Calculate convective heat flux using the Bartz correlation. + + The Bartz equation estimates the convective heat transfer coefficient: + + h = (0.026/D_t^0.2) * (mu^0.2 * cp / Pr^0.6) * (p_c / c*)^0.8 * + (D_t/R_c)^0.1 * (A_t/A)^0.9 * sigma + + where sigma is a correction factor for property variations. + + Args: + chamber_pressure: Chamber pressure [Pa] + chamber_temp: Chamber temperature [K] + throat_diameter: Throat diameter [m] + local_diameter: Local diameter at evaluation point [m] + characteristic_velocity: c* [m/s] + gamma: Ratio of specific heats [-] + molecular_weight: Molecular weight [kg/kmol] + local_mach: Local Mach number [-] + wall_temp: Wall temperature [K]. If None, estimates at 600K. + + Returns: + Heat flux [W/m²] + """ + # Extract values in SI + pc = chamber_pressure.to("Pa").value + Tc = chamber_temp.to("K").value + Dt = throat_diameter.to("m").value + D = local_diameter.to("m").value + cstar = characteristic_velocity.to("m/s").value + MW = molecular_weight + + # Wall temperature (estimate if not provided) + Tw = wall_temp.to("K").value if wall_temp else 600.0 + + # Gas properties + R_specific = 8314.46 / MW # J/(kg·K) + cp = gamma * R_specific / (gamma - 1) # J/(kg·K) + + # Estimate viscosity using Sutherland's law + # Reference: air at 273K, mu0 = 1.71e-5 Pa·s + # For combustion products, use higher reference + mu_ref = 4e-5 # Pa·s at Tc_ref + Tc_ref = 3000 # K + S = 200 # Sutherland constant (approximate for combustion products) + + mu = mu_ref * (Tc / Tc_ref) ** 1.5 * (Tc_ref + S) / (Tc + S) + + # Prandtl number + # Pr = mu * cp / k, approximate k from cp and Pr ~ 0.7-0.9 + Pr = 0.8 # Typical for combustion products + + # Area ratio + area_ratio = (D / Dt) ** 2 + + # Sigma correction factor for property variation across boundary layer + # sigma = 1 / [(Tw/Tc * (1 + (gamma-1)/2 * M^2) + 0.5)^0.68 * + # (1 + (gamma-1)/2 * M^2)^0.12] + gm1 = gamma - 1 + temp_factor = Tw / Tc * (1 + gm1/2 * local_mach**2) + 0.5 + sigma = 1 / (temp_factor ** 0.68 * (1 + gm1/2 * local_mach**2) ** 0.12) + + # Bartz correlation for heat transfer coefficient + # h = (0.026 / Dt^0.2) * (mu^0.2 * cp / Pr^0.6) * (pc/cstar)^0.8 * + # (Dt/Dt)^0.1 * (At/A)^0.9 * sigma + h = (0.026 / Dt**0.2) * (mu**0.2 * cp / Pr**0.6) * \ + (pc / cstar)**0.8 * (1/area_ratio)**0.9 * sigma + + # Calculate adiabatic wall temperature + r = recovery_factor(Pr) + T_aw = Tc * (1 + r * gm1/2 * local_mach**2) / (1 + gm1/2 * local_mach**2) + + # Heat flux + q = h * (T_aw - Tw) + + return Quantity(q, "Pa", "pressure") # W/m² = Pa·m/s, using Pa as proxy + + +@beartype +def estimate_heat_flux( + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + location: str = "throat", + wall_temp: Quantity | None = None, +) -> Quantity: + """Estimate heat flux at a specific location in the engine. + + Provides a simplified interface to the Bartz correlation. + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + location: Location to evaluate: "throat", "chamber", or "exit" + wall_temp: Wall temperature [K]. If None, estimates based on location. + + Returns: + Heat flux [W/m²] + """ + # Determine local conditions based on location + if location == "throat": + local_diameter = geometry.throat_diameter + local_mach = 1.0 + default_wall_temp = 700 # K, hot at throat + elif location == "chamber": + local_diameter = geometry.chamber_diameter + local_mach = 0.1 # Low Mach in chamber + default_wall_temp = 600 # K + elif location == "exit": + local_diameter = geometry.exit_diameter + local_mach = performance.exit_mach + default_wall_temp = 400 # K, cooler at exit + else: + raise ValueError(f"Unknown location: {location}. Use 'throat', 'chamber', or 'exit'") + + if wall_temp is None: + wall_temp = kelvin(default_wall_temp) + + return bartz_heat_flux( + chamber_pressure=inputs.chamber_pressure, + chamber_temp=inputs.chamber_temp, + throat_diameter=geometry.throat_diameter, + local_diameter=local_diameter, + characteristic_velocity=performance.cstar, + gamma=inputs.gamma, + molecular_weight=inputs.molecular_weight, + local_mach=local_mach, + wall_temp=wall_temp, + ) + + +@beartype +def heat_flux_profile( + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + n_points: int = 50, + wall_temp: Quantity | None = None, +) -> tuple[list[float], list[float]]: + """Calculate heat flux profile along the engine. + + Returns heat flux from chamber through throat to exit. + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + n_points: Number of points in profile + wall_temp: Wall temperature (constant along length) + + Returns: + Tuple of (x_positions, heat_fluxes) where x is normalized (0=chamber, 1=exit) + """ + # Get key dimensions + Dc = geometry.chamber_diameter.to("m").value + Dt = geometry.throat_diameter.to("m").value + De = geometry.exit_diameter.to("m").value + + # Generate normalized positions + x_norm = np.linspace(0, 1, n_points) + + # Estimate local diameter and Mach number along engine + # Simplified: linear convergent, bell divergent + diameters = [] + machs = [] + + for x in x_norm: + if x < 0.3: + # Chamber region + D = Dc + M = 0.1 + x * 0.3 # Low, slowly increasing + elif x < 0.5: + # Convergent section + frac = (x - 0.3) / 0.2 + D = Dc - frac * (Dc - Dt) + M = 0.1 + frac * 0.9 + elif x < 0.55: + # Throat region + D = Dt + M = 1.0 + else: + # Divergent section + frac = (x - 0.55) / 0.45 + D = Dt + frac * (De - Dt) + # Simple approximation for supersonic Mach + M = 1.0 + frac * (performance.exit_mach - 1.0) + + diameters.append(D) + machs.append(max(0.1, M)) + + # Calculate heat flux at each point + heat_fluxes = [] + for D, M in zip(diameters, machs, strict=True): + q = bartz_heat_flux( + chamber_pressure=inputs.chamber_pressure, + chamber_temp=inputs.chamber_temp, + throat_diameter=geometry.throat_diameter, + local_diameter=Quantity(D, "m", "length"), + characteristic_velocity=performance.cstar, + gamma=inputs.gamma, + molecular_weight=inputs.molecular_weight, + local_mach=M, + wall_temp=wall_temp or kelvin(600), + ) + heat_fluxes.append(q.value) + + return list(x_norm), heat_fluxes + diff --git a/rocket/thermal/regenerative.py b/rocket/thermal/regenerative.py new file mode 100644 index 0000000..1167713 --- /dev/null +++ b/rocket/thermal/regenerative.py @@ -0,0 +1,451 @@ +"""Regenerative cooling feasibility analysis. + +Regenerative cooling uses the fuel (or oxidizer) as a coolant, flowing +through channels in the chamber/nozzle wall before injection. This is +the most common cooling method for high-performance rocket engines. + +This module provides screening-level analysis to determine if a given +engine design can be regeneratively cooled within material limits. + +Key considerations: +- Heat flux at the throat (highest) +- Coolant heat capacity and flow rate +- Wall material temperature limits +- Coolant-side pressure drop + +References: + - Huzel & Huang, Chapter 4 + - Sutton & Biblarz, Chapter 8 +""" + +import math +from dataclasses import dataclass + +from beartype import beartype + +from rocket.engine import EngineGeometry, EngineInputs, EnginePerformance +from rocket.thermal.heat_flux import estimate_heat_flux +from rocket.units import Quantity, kelvin, pascals + +# ============================================================================= +# Coolant Properties Database +# ============================================================================= + + +@beartype +@dataclass(frozen=True, slots=True) +class CoolantProperties: + """Thermophysical properties of a coolant. + + Properties are approximate values at typical operating conditions. + For detailed design, use property tables or CoolProp. + + Attributes: + name: Coolant name + density: Liquid density [kg/m³] + specific_heat: Specific heat capacity [J/(kg·K)] + thermal_conductivity: Thermal conductivity [W/(m·K)] + viscosity: Dynamic viscosity [Pa·s] + boiling_point: Boiling point at 1 atm [K] + max_temp: Maximum recommended temperature before decomposition [K] + """ + + name: str + density: float + specific_heat: float + thermal_conductivity: float + viscosity: float + boiling_point: float + max_temp: float + + +# Coolant property database +# Values are approximate at typical inlet conditions +COOLANT_DATABASE: dict[str, CoolantProperties] = { + "RP1": CoolantProperties( + name="RP-1 (Kerosene)", + density=810.0, + specific_heat=2000.0, + thermal_conductivity=0.12, + viscosity=0.0015, + boiling_point=490.0, + max_temp=600.0, # Coking limit + ), + "CH4": CoolantProperties( + name="Liquid Methane", + density=422.0, + specific_heat=3500.0, + thermal_conductivity=0.19, + viscosity=0.00012, + boiling_point=111.0, + max_temp=500.0, # Before significant decomposition + ), + "LH2": CoolantProperties( + name="Liquid Hydrogen", + density=70.8, + specific_heat=14300.0, + thermal_conductivity=0.10, + viscosity=0.000013, + boiling_point=20.0, + max_temp=300.0, # Stays liquid/supercritical + ), + "Ethanol": CoolantProperties( + name="Ethanol", + density=789.0, + specific_heat=2440.0, + thermal_conductivity=0.17, + viscosity=0.0011, + boiling_point=351.0, + max_temp=500.0, + ), + "N2O4": CoolantProperties( + name="Nitrogen Tetroxide", + density=1450.0, + specific_heat=1560.0, + thermal_conductivity=0.12, + viscosity=0.0004, + boiling_point=294.0, + max_temp=400.0, + ), + "MMH": CoolantProperties( + name="Monomethylhydrazine", + density=878.0, + specific_heat=2920.0, + thermal_conductivity=0.22, + viscosity=0.0008, + boiling_point=360.0, + max_temp=450.0, + ), +} + + +@beartype +def get_coolant_properties(coolant: str) -> CoolantProperties: + """Get properties for a coolant. + + Args: + coolant: Coolant name (e.g., "RP1", "CH4", "LH2") + + Returns: + CoolantProperties for the coolant + + Raises: + ValueError: If coolant not found in database + """ + # Normalize name + name = coolant.upper().replace("-", "").replace(" ", "") + + for key, props in COOLANT_DATABASE.items(): + if key.upper() == name: + return props + + available = list(COOLANT_DATABASE.keys()) + raise ValueError(f"Unknown coolant '{coolant}'. Available: {available}") + + +# ============================================================================= +# Cooling Feasibility Analysis +# ============================================================================= + + +@beartype +@dataclass(frozen=True, slots=True) +class CoolingFeasibility: + """Results of regenerative cooling feasibility analysis. + + Attributes: + feasible: Whether cooling is feasible within constraints + max_wall_temp: Maximum predicted wall temperature [K] + max_allowed_temp: Maximum allowed wall temperature [K] + coolant_temp_rise: Temperature rise of coolant [K] + coolant_outlet_temp: Predicted coolant outlet temperature [K] + throat_heat_flux: Heat flux at throat [W/m²] + total_heat_load: Total heat load to coolant [W] + required_coolant_flow: Minimum required coolant flow [kg/s] + available_coolant_flow: Available coolant flow (fuel or ox) [kg/s] + flow_margin: Ratio of available/required flow [-] + pressure_drop: Estimated coolant-side pressure drop [Pa] + channel_velocity: Estimated coolant velocity in channels [m/s] + warnings: List of warnings or concerns + """ + + feasible: bool + max_wall_temp: Quantity + max_allowed_temp: Quantity + coolant_temp_rise: Quantity + coolant_outlet_temp: Quantity + throat_heat_flux: Quantity + total_heat_load: Quantity + required_coolant_flow: Quantity + available_coolant_flow: Quantity + flow_margin: float + pressure_drop: Quantity + channel_velocity: float + warnings: list[str] + + +@beartype +def check_cooling_feasibility( + inputs: EngineInputs, + performance: EnginePerformance, + geometry: EngineGeometry, + coolant: str, + coolant_inlet_temp: Quantity | None = None, + max_wall_temp: Quantity | None = None, + wall_material: str = "copper_alloy", + num_channels: int | None = None, + channel_aspect_ratio: float = 3.0, +) -> CoolingFeasibility: + """Check if regenerative cooling is feasible for this engine design. + + This is a screening-level analysis that estimates whether the engine + can be cooled within material limits using the available coolant flow. + + Args: + inputs: Engine input parameters + performance: Computed engine performance + geometry: Computed engine geometry + coolant: Coolant name (typically the fuel: "RP1", "CH4", "LH2") + coolant_inlet_temp: Coolant inlet temperature [K]. Defaults to storage temp. + max_wall_temp: Maximum allowed wall temperature [K]. Defaults based on material. + wall_material: Wall material for thermal limits + num_channels: Number of cooling channels. If None, estimated. + channel_aspect_ratio: Channel height/width ratio + + Returns: + CoolingFeasibility assessment + """ + warnings: list[str] = [] + + # Get coolant properties + try: + coolant_props = get_coolant_properties(coolant) + except ValueError: + # Use generic properties if unknown + coolant_props = CoolantProperties( + name=coolant, + density=800.0, + specific_heat=2000.0, + thermal_conductivity=0.15, + viscosity=0.001, + boiling_point=300.0, + max_temp=500.0, + ) + warnings.append(f"Unknown coolant '{coolant}', using generic properties") + + # Set default inlet temperature (storage temperature) + if coolant_inlet_temp is None: + # Cryogenic propellants near boiling point, storables at ~300K + if coolant.upper() in ["LH2", "CH4", "LOX"]: + T_inlet = coolant_props.boiling_point + 10 # Just above boiling + else: + T_inlet = 300.0 # Room temperature + else: + T_inlet = coolant_inlet_temp.to("K").value + + # Set default max wall temperature based on material + if max_wall_temp is None: + wall_temps = { + "copper_alloy": 800, # GRCop-84, NARloy-Z + "nickel_alloy": 1000, # Inconel 718 + "steel": 700, # Stainless steel + "niobium": 1500, # Refractory metal + } + T_wall_max = wall_temps.get(wall_material, 800) + else: + T_wall_max = max_wall_temp.to("K").value + + # Estimate heat flux at throat (worst case) + q_throat = estimate_heat_flux( + inputs, performance, geometry, + location="throat", + wall_temp=kelvin(T_wall_max * 0.9), # Assume wall near limit + ) + + # Also get chamber and exit heat fluxes for total heat load + q_chamber = estimate_heat_flux(inputs, performance, geometry, location="chamber") + q_exit = estimate_heat_flux(inputs, performance, geometry, location="exit") + + # Estimate total heat load + # Simplified: use average heat flux × total surface area + Dt = geometry.throat_diameter.to("m").value + Dc = geometry.chamber_diameter.to("m").value + De = geometry.exit_diameter.to("m").value + Lc = geometry.chamber_length.to("m").value + Ln = geometry.nozzle_length.to("m").value + + # Surface areas (approximate) + A_chamber = math.pi * Dc * Lc + A_convergent = math.pi * (Dc + Dt) / 2 * (Dc - Dt) / (2 * math.tan(math.radians(45))) * 0.5 + A_throat = math.pi * Dt * Dt * 0.1 # Small throat region + A_divergent = math.pi * (Dt + De) / 2 * Ln * 0.7 # Approximate bell surface + + # Average heat fluxes for each region (throat is highest) + q_avg_chamber = q_chamber.value * 0.3 # Lower than throat + q_avg_convergent = (q_chamber.value + q_throat.value) / 2 + q_avg_throat = q_throat.value + q_avg_divergent = (q_throat.value + q_exit.value) / 2 + + # Total heat load + Q_total = ( + q_avg_chamber * A_chamber + + q_avg_convergent * A_convergent + + q_avg_throat * A_throat + + q_avg_divergent * A_divergent + ) + + # Available coolant flow (assume fuel is coolant) + mdot_coolant_available = performance.mdot_fuel.to("kg/s").value + + # Required coolant flow to absorb heat without exceeding temperature limit + # Q = mdot * cp * delta_T + # delta_T_max = T_coolant_max - T_inlet + T_coolant_max = min(coolant_props.max_temp, T_wall_max - 100) # Stay below wall + delta_T_max = T_coolant_max - T_inlet + + if delta_T_max <= 0: + warnings.append("Coolant inlet temperature exceeds maximum allowable") + delta_T_max = 100 # Use minimum for calculation + + mdot_coolant_required = Q_total / (coolant_props.specific_heat * delta_T_max) + + # Flow margin + flow_margin = mdot_coolant_available / mdot_coolant_required if mdot_coolant_required > 0 else float('inf') + + # Actual temperature rise with available flow + if mdot_coolant_available > 0: + delta_T_actual = Q_total / (mdot_coolant_available * coolant_props.specific_heat) + else: + delta_T_actual = float('inf') + + T_coolant_out = T_inlet + delta_T_actual + + # Estimate wall temperature + # T_wall = T_coolant + Q/(h_coolant * A) + # For screening, use correlation: T_wall ~ T_coolant + q * (t_wall / k_wall + 1/h_coolant) + # Simplified: wall runs ~100-200K above coolant + T_wall_estimate = T_coolant_out + 150 # K above coolant + + # Pressure drop estimation + # Number of channels (estimate if not provided) + if num_channels is None: + # Roughly 1 channel per mm of circumference at throat + num_channels = max(20, int(math.pi * Dt * 1000)) + + # Channel dimensions (rough estimate) + channel_width = (math.pi * Dt) / num_channels * 0.6 # 60% channel, 40% rib + channel_height = channel_width * channel_aspect_ratio + channel_area = channel_width * channel_height + + # Coolant velocity + total_channel_area = num_channels * channel_area + v_coolant = mdot_coolant_available / (coolant_props.density * total_channel_area) + + # Pressure drop (Darcy-Weisbach approximation) + # ΔP = f * (L/D_h) * (ρ * v²/2) + D_h = 4 * channel_area / (2 * (channel_width + channel_height)) # Hydraulic diameter + Re = coolant_props.density * v_coolant * D_h / coolant_props.viscosity + f = 0.316 / Re**0.25 if Re > 2300 else 64 / max(Re, 100) # Friction factor + + L_total = Lc + Ln # Total cooled length + dp = f * (L_total / D_h) * (coolant_props.density * v_coolant**2 / 2) + + # Add losses for bends, manifolds + dp *= 1.5 + + # Feasibility assessment + feasible = True + + if T_wall_estimate > T_wall_max: + feasible = False + warnings.append( + f"Estimated wall temp {T_wall_estimate:.0f}K exceeds limit {T_wall_max:.0f}K" + ) + + if flow_margin < 1.0: + feasible = False + warnings.append( + f"Insufficient coolant flow: need {mdot_coolant_required:.2f} kg/s, " + f"have {mdot_coolant_available:.2f} kg/s" + ) + + if T_coolant_out > coolant_props.max_temp: + warnings.append( + f"Coolant outlet temp {T_coolant_out:.0f}K exceeds max {coolant_props.max_temp:.0f}K" + ) + + if v_coolant > 50: + warnings.append(f"High coolant velocity {v_coolant:.1f} m/s may cause erosion") + + if dp > 5e6: + warnings.append(f"High pressure drop {dp/1e6:.1f} MPa") + + # Heat flux at throat check + if q_throat.value > 50e6: # > 50 MW/m² + warnings.append( + f"Very high throat heat flux {q_throat.value/1e6:.1f} MW/m² - " + "film cooling may be needed" + ) + + return CoolingFeasibility( + feasible=feasible, + max_wall_temp=kelvin(T_wall_estimate), + max_allowed_temp=kelvin(T_wall_max), + coolant_temp_rise=kelvin(delta_T_actual), + coolant_outlet_temp=kelvin(T_coolant_out), + throat_heat_flux=q_throat, + total_heat_load=Quantity(Q_total, "N", "force"), # W, using N as proxy + required_coolant_flow=Quantity(mdot_coolant_required, "kg/s", "mass_flow"), + available_coolant_flow=Quantity(mdot_coolant_available, "kg/s", "mass_flow"), + flow_margin=flow_margin, + pressure_drop=pascals(dp), + channel_velocity=v_coolant, + warnings=warnings, + ) + + +@beartype +def format_cooling_summary(result: CoolingFeasibility) -> str: + """Format cooling feasibility results as readable string. + + Args: + result: CoolingFeasibility from check_cooling_feasibility() + + Returns: + Formatted multi-line string + """ + status = "✓ FEASIBLE" if result.feasible else "✗ NOT FEASIBLE" + + lines = [ + f"{'=' * 60}", + "REGENERATIVE COOLING ANALYSIS", + f"Status: {status}", + f"{'=' * 60}", + "", + "THERMAL:", + f" Throat Heat Flux: {result.throat_heat_flux.value/1e6:.1f} MW/m²", + f" Total Heat Load: {result.total_heat_load.value/1e6:.2f} MW", + f" Max Wall Temp: {result.max_wall_temp.value:.0f} K", + f" Allowed Wall Temp: {result.max_allowed_temp.value:.0f} K", + "", + "COOLANT:", + f" Temperature Rise: {result.coolant_temp_rise.value:.0f} K", + f" Outlet Temperature: {result.coolant_outlet_temp.value:.0f} K", + f" Required Flow: {result.required_coolant_flow.value:.2f} kg/s", + f" Available Flow: {result.available_coolant_flow.value:.2f} kg/s", + f" Flow Margin: {result.flow_margin:.2f}x", + "", + "HYDRAULICS:", + f" Channel Velocity: {result.channel_velocity:.1f} m/s", + f" Pressure Drop: {result.pressure_drop.to('bar').value:.1f} bar", + ] + + if result.warnings: + lines.extend(["", "WARNINGS:"]) + for warning in result.warnings: + lines.append(f" ⚠ {warning}") + + lines.append(f"{'=' * 60}") + + return "\n".join(lines) + diff --git a/rocket/units.py b/rocket/units.py index dfcff30..d5444fa 100644 --- a/rocket/units.py +++ b/rocket/units.py @@ -99,6 +99,11 @@ # Density "kg/m^3": (1.0, "density"), "lbm/ft^3": (16.0185, "density"), + # Power + "W": (1.0, "power"), + "kW": (1000.0, "power"), + "MW": (1e6, "power"), + "hp": (745.7, "power"), # Mechanical horsepower # Specific impulse (time dimension but special meaning) # Note: Isp in seconds is the same in SI and Imperial # Dimensionless @@ -610,6 +615,30 @@ def dimensionless(value: float | int) -> Quantity: return Quantity(value, "1", "dimensionless") +@beartype +def watts(value: float | int) -> Quantity: + """Create a power quantity in Watts.""" + return Quantity(value, "W", "power") + + +@beartype +def kilowatts(value: float | int) -> Quantity: + """Create a power quantity in kilowatts.""" + return Quantity(value, "kW", "power") + + +@beartype +def megawatts(value: float | int) -> Quantity: + """Create a power quantity in megawatts.""" + return Quantity(value, "MW", "power") + + +@beartype +def horsepower(value: float | int) -> Quantity: + """Create a power quantity in horsepower.""" + return Quantity(value, "hp", "power") + + # ============================================================================= # Constants # ============================================================================= diff --git a/uv.lock b/uv.lock index 8c3b374..464ca5f 100644 --- a/uv.lock +++ b/uv.lock @@ -651,6 +651,32 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/54/20/4d324d65cc6d9205fabedc306948156824eb9f0ee1633355a8f7ec5c66bf/pluggy-1.6.0-py3-none-any.whl", hash = "sha256:e920276dd6813095e9377c0bc5566d94c932c33b27a3e3945d8389c374dd4746", size = 20538, upload-time = "2025-05-15T12:30:06.134Z" }, ] +[[package]] +name = "polars" +version = "1.35.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "polars-runtime-32" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/fa/43/09d4738aa24394751cb7e5d1fc4b5ef461d796efcadd9d00c79578332063/polars-1.35.2.tar.gz", hash = "sha256:ae458b05ca6e7ca2c089342c70793f92f1103c502dc1b14b56f0a04f2cc1d205", size = 694895, upload-time = "2025-11-09T13:20:05.921Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/b4/9a/24e4b890c7ee4358964aa92c4d1865df0e8831f7df6abaa3a39914521724/polars-1.35.2-py3-none-any.whl", hash = "sha256:5e8057c8289ac148c793478323b726faea933d9776bd6b8a554b0ab7c03db87e", size = 783597, upload-time = "2025-11-09T13:18:51.361Z" }, +] + +[[package]] +name = "polars-runtime-32" +version = "1.35.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/cb/75/ac1256ace28c832a0997b20ba9d10a9d3739bd4d457c1eb1e7d196b6f88b/polars_runtime_32-1.35.2.tar.gz", hash = "sha256:6e6e35733ec52abe54b7d30d245e6586b027d433315d20edfb4a5d162c79fe90", size = 2694387, upload-time = "2025-11-09T13:20:07.624Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/66/de/a532b81e68e636483a5dd764d72e106215543f3ef49a142272b277ada8fe/polars_runtime_32-1.35.2-cp39-abi3-macosx_10_12_x86_64.whl", hash = "sha256:e465d12a29e8df06ea78947e50bd361cdf77535cd904fd562666a8a9374e7e3a", size = 40524507, upload-time = "2025-11-09T13:18:55.727Z" }, + { url = "https://files.pythonhosted.org/packages/2d/0b/679751ea6aeaa7b3e33a70ba17f9c8150310792583f3ecf9bb1ce15fe15c/polars_runtime_32-1.35.2-cp39-abi3-macosx_11_0_arm64.whl", hash = "sha256:ef2b029b78f64fb53f126654c0bfa654045c7546bd0de3009d08bd52d660e8cc", size = 36700154, upload-time = "2025-11-09T13:18:59.78Z" }, + { url = "https://files.pythonhosted.org/packages/e2/c8/fd9f48dd6b89ae9cff53d896b51d08579ef9c739e46ea87a647b376c8ca2/polars_runtime_32-1.35.2-cp39-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:85dda0994b5dff7f456bb2f4bbd22be9a9e5c5e28670e23fedb13601ec99a46d", size = 41317788, upload-time = "2025-11-09T13:19:03.949Z" }, + { url = "https://files.pythonhosted.org/packages/67/89/e09d9897a70b607e22a36c9eae85a5b829581108fd1e3d4292e5c0f52939/polars_runtime_32-1.35.2-cp39-abi3-manylinux_2_24_aarch64.whl", hash = "sha256:3b9006902fc51b768ff747c0f74bd4ce04005ee8aeb290ce9c07ce1cbe1b58a9", size = 37850590, upload-time = "2025-11-09T13:19:08.154Z" }, + { url = "https://files.pythonhosted.org/packages/dc/40/96a808ca5cc8707894e196315227f04a0c82136b7fb25570bc51ea33b88d/polars_runtime_32-1.35.2-cp39-abi3-win_amd64.whl", hash = "sha256:ddc015fac39735592e2e7c834c02193ba4d257bb4c8c7478b9ebe440b0756b84", size = 41290019, upload-time = "2025-11-09T13:19:12.214Z" }, + { url = "https://files.pythonhosted.org/packages/f4/d1/8d1b28d007da43c750367c8bf5cb0f22758c16b1104b2b73b9acadb2d17a/polars_runtime_32-1.35.2-cp39-abi3-win_arm64.whl", hash = "sha256:6861145aa321a44eda7cc6694fb7751cb7aa0f21026df51b5faa52e64f9dc39b", size = 36955684, upload-time = "2025-11-09T13:19:15.666Z" }, +] + [[package]] name = "pygments" version = "2.19.2" @@ -720,7 +746,9 @@ dependencies = [ { name = "matplotlib" }, { name = "numba" }, { name = "numpy" }, + { name = "polars" }, { name = "rocketcea" }, + { name = "tqdm" }, ] [package.optional-dependencies] @@ -736,10 +764,12 @@ requires-dist = [ { name = "matplotlib", specifier = ">=3.9" }, { name = "numba", specifier = ">=0.60" }, { name = "numpy", specifier = ">=2.0" }, + { name = "polars", specifier = ">=1.35.0" }, { name = "pytest", marker = "extra == 'dev'", specifier = ">=8.0" }, { name = "pytest-cov", marker = "extra == 'dev'", specifier = ">=4.0" }, { name = "rocketcea", specifier = ">=1.2.1" }, { name = "ruff", marker = "extra == 'dev'", specifier = ">=0.4" }, + { name = "tqdm", specifier = ">=4.66" }, ] provides-extras = ["dev"] @@ -912,3 +942,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/84/ff/426ca8683cf7b753614480484f6437f568fd2fda2edbdf57a2d3d8b27a0b/tomli-2.3.0-cp314-cp314t-win_amd64.whl", hash = "sha256:70a251f8d4ba2d9ac2542eecf008b3c8a9fc5c3f9f02c56a9d7952612be2fdba", size = 119756, upload-time = "2025-10-08T22:01:45.234Z" }, { url = "https://files.pythonhosted.org/packages/77/b8/0135fadc89e73be292b473cb820b4f5a08197779206b33191e801feeae40/tomli-2.3.0-py3-none-any.whl", hash = "sha256:e95b1af3c5b07d9e643909b5abbec77cd9f1217e6d0bca72b0234736b9fb1f1b", size = 14408, upload-time = "2025-10-08T22:01:46.04Z" }, ] + +[[package]] +name = "tqdm" +version = "4.67.1" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "colorama", marker = "sys_platform == 'win32'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/a8/4b/29b4ef32e036bb34e4ab51796dd745cdba7ed47ad142a9f4a1eb8e0c744d/tqdm-4.67.1.tar.gz", hash = "sha256:f8aef9c52c08c13a65f30ea34f4e5aac3fd1a34959879d7e59e63027286627f2", size = 169737, upload-time = "2024-11-24T20:12:22.481Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl", hash = "sha256:26445eca388f82e72884e0d580d5464cd801a3ea01e63e5601bdff9ba6a48de2", size = 78540, upload-time = "2024-11-24T20:12:19.698Z" }, +]