diff --git a/main.py b/main.py
deleted file mode 100644
index a4cc348fa227c1ad6bea95c03b902bb5d97b08aa..0000000000000000000000000000000000000000
--- a/main.py
+++ /dev/null
@@ -1,429 +0,0 @@
-import math
-import numpy as np
-from iapws import IAPWS97, iapws97
-
-class HeatPipe:
- def __init__(self, length, diameter, working_fluid="water", wick_type="sintered"):
- """
- Initialize heat pipe parameters
-
- Args:
- length (float): Total length of heat pipe in meters
- diameter (float): Outer diameter of heat pipe in meters
- working_fluid (str): Working fluid (default: water)
- wick_type (str): Type of wick structure (default: sintered)
- """
- self.length = length
- self.diameter = diameter
- self.working_fluid = working_fluid.lower()
- self.wick_type = wick_type.lower()
-
- # Wick properties dictionary
- self.wick_properties = {
- "sintered": {
- "porosity": 0.5,
- "permeability": 1e-12, # m²
- "thermal_conductivity": 5.0, # W/m·K
- "pore_radius": 50e-6, # m
- "particle_diameter": 100e-6 # m
- },
- "screen": {
- "porosity": 0.7,
- "permeability": 1e-10, # m²
- "thermal_conductivity": 3.0, # W/m·K
- "pore_radius": 100e-6, # m
- "wire_diameter": 100e-6 # m
- },
- "bcc_lattice": {
- # These will be calculated based on unit cell parameters
- "porosity": None,
- "permeability": None,
- "thermal_conductivity": None,
- "pore_radius": None,
- "strut_diameter": None,
- "unit_cell_size": None
- }
- }
-
- def get_water_properties(self, temperature, quality=None):
- """
- Get water properties using IAPWS97 formulations
-
- Args:
- temperature (float): Temperature in Kelvin
- quality (float): Steam quality (0 = saturated liquid, 1 = saturated vapor, None = specify temperature only)
-
- Returns:
- dict: Dictionary of water properties
- """
- try:
- if quality is not None:
- # Two-phase region (saturated mixture)
- water = IAPWS97(T=temperature, x=quality)
-
- # Get saturation properties at this temperature
- sat_liquid = IAPWS97(T=temperature, x=0)
- sat_vapor = IAPWS97(T=temperature, x=1)
-
- properties = {
- "latent_heat": (sat_vapor.h - sat_liquid.h) * 1000, # J/kg
- "liquid_density": sat_liquid.rho, # kg/m³
- "vapor_density": sat_vapor.rho, # kg/m³
- "liquid_viscosity": sat_liquid.mu, # Pa·s
- "vapor_viscosity": sat_vapor.mu, # Pa·s
- "surface_tension": water.sigma, # N/m
- "thermal_conductivity": sat_liquid.k, # W/m·K
- "specific_heat_liquid": sat_liquid.cp * 1000, # J/kg·K
- "specific_heat_vapor": sat_vapor.cp * 1000, # J/kg.K
- "gamma0_vapor": sat_vapor.cp/sat_vapor.cv,
- "vapor_pressure": water.P * 1e6, # Pa
- "critical_pressure": iapws97.Pc * 1e6, # Pa
- "critical_temperature": iapws97.Tc, # K
- "sonic_velocity": sat_vapor.w # m/s
- }
- else:
- # Single phase region (specify temperature only)
- water = IAPWS97(T=temperature, P=101325) # Assume atmospheric pressure
- properties = {
- "liquid_density": water.rho,
- "thermal_conductivity": water.k,
- "specific_heat": water.cp,
- "liquid_viscosity": water.mu
- }
-
- return properties
- except Exception as e:
- raise ValueError(f"Error calculating water properties: {str(e)}")
-
- def calculate_bcc_properties(self, strut_diameter, unit_cell_size, material_conductivity=400): #need to verify this!
- """
- Calculate properties for BCC lattice wick structure
-
- Args:
- strut_diameter (float): Diameter of the struts in meters
- unit_cell_size (float): Size of the unit cell in meters
- material_conductivity (float): Thermal conductivity of the strut material (default: copper)
-
- Returns:
- dict: Updated properties for BCC lattice structure
- """
- # Volume calculations
- V_unit_cell = unit_cell_size**3
-
- # Length of diagonal strut
- l_diagonal = math.sqrt(3) * unit_cell_size
-
- # Number of struts per unit cell (4 diagonal + 12 edge struts / 4 shared)
- n_struts_diagonal = 4
- n_struts_edge = 12/4
-
- # Volume of struts
- V_strut_diagonal = math.pi * (strut_diameter/2)**2 * l_diagonal
- V_strut_edge = math.pi * (strut_diameter/2)**2 * unit_cell_size
-
- # Total volume of solid material in unit cell
- V_solid = n_struts_diagonal * V_strut_diagonal + n_struts_edge * V_strut_edge
-
- # Calculate porosity
- porosity = 1 - (V_solid / V_unit_cell)
-
- # Calculate permeability using modified Carman-Kozeny equation
- # Specific surface area for BCC lattice
- surface_area = (n_struts_diagonal * math.pi * strut_diameter * l_diagonal +
- n_struts_edge * math.pi * strut_diameter * unit_cell_size)
- specific_surface = surface_area / V_unit_cell
-
- # Kozeny constant for BCC lattice (typical range 4.5-5.5)
- K_kozeny = 5.0
-
- # Calculate permeability
- permeability = (porosity**3) / (K_kozeny * specific_surface**2 * (1-porosity)**2)
-
- # Calculate effective thermal conductivity using Maxwell's model
- # Assuming copper as base material and water as fluid
- k_fluid = self.fluid_properties["water"]["thermal_conductivity"]
- k_solid = material_conductivity
-
- # Maxwell's model for effective thermal conductivity
- num = k_solid + 2*k_fluid + 2*(k_solid - k_fluid)*porosity
- den = k_solid + 2*k_fluid - (k_solid - k_fluid)*porosity
- k_eff = k_fluid * (num/den)
-
- # Calculate effective pore radius
- # For BCC lattice, this is approximately the radius of the largest sphere
- # that can fit in the void space
- pore_radius = (unit_cell_size - 2*strut_diameter) / 4
-
- # Update wick properties for BCC lattice
- self.wick_properties["bcc_lattice"].update({
- "porosity": porosity,
- "permeability": permeability,
- "thermal_conductivity": k_eff,
- "pore_radius": pore_radius,
- "strut_diameter": strut_diameter,
- "unit_cell_size": unit_cell_size,
- "specific_surface_area": specific_surface
- })
-
- return self.wick_properties["bcc_lattice"]
-
- def initialize_bcc_wick(self, strut_diameter, unit_cell_size): #need to verify this!
- """
- Initialize BCC lattice wick structure with given parameters
-
- Args:
- strut_diameter (float): Diameter of the struts in meters
- unit_cell_size (float): Size of the unit cell in meters
- """
- if self.wick_type == "bcc_lattice":
- self.calculate_bcc_properties(strut_diameter, unit_cell_size)
- else:
- raise ValueError("Heat pipe not initialized with BCC lattice wick type")
-
- def calculate_thermal_resistance(self, heat_load, evap_length, cond_length,
- wall_thickness=0.001, wick_thickness=0.0005, operating_temp=353.15, material='copper'):
- #interesting that heat load is not necessary to calculate only wall and wick resistances..
- """
- Calculate total thermal resistance of the heat pipe
-
- Args:
- heat_load (float): Heat load in Watts
- evap_length (float): Length of evaporator section in meters
- cond_length (float): Length of condenser section in meters
- wall_thickness (float): Thickness of the heat pipe wall in meters
- wick_thickness (float): Thickness of the wick structure in meters
- operating_temp (float): Operating temperature in Kelvin
-
- Returns:
- dict: Dictionary containing thermal resistances and total resistance
- """
- # dict of thermal conductivity values of different metals
- thermal_conductivity = {
- 'copper': 400,
- 'steel': 45,
- 'aluminum': 237
- }
- # Get fluid properties using IAPWS
- fluid = self.get_water_properties(operating_temp, quality=0) # Liquid properties
- vapor = self.get_water_properties(operating_temp, quality=1) # Vapor properties
-
- # Get wick properties
- wick = self.wick_properties[self.wick_type]
-
- # Calculate cross-sectional areas
- # area_wall = math.pi * self.diameter * wall_thickness #how?
- # area_wick = math.pi * (self.diameter - 2*wall_thickness) * wick_thickness #how?
- area_vapor = math.pi * (self.diameter - 2*(wall_thickness + wick_thickness))**2 / 4
-
- # 1. Wall thermal resistance (evaporator and condenser)
- k_wall = thermal_conductivity[material.lower()] # Thermal conductivity of copper (W/m·K)
- R_wall_e = math.log((self.diameter)/(self.diameter - 2*wall_thickness)) / (2 * math.pi * k_wall * evap_length)
- R_wall_c = math.log((self.diameter)/(self.diameter - 2*wall_thickness)) / (2 * math.pi * k_wall * cond_length)
-
- # 2. Wick thermal resistance
- k_eff_wick = wick["thermal_conductivity"] * (1 - wick["porosity"]) + fluid["thermal_conductivity"] * wick["porosity"]
- R_wick_e = math.log((self.diameter - 2*wall_thickness)/(self.diameter - 2*(wall_thickness + wick_thickness))) / (2 * math.pi * k_eff_wick * evap_length)
- R_wick_c = math.log((self.diameter - 2*wall_thickness)/(self.diameter - 2*(wall_thickness + wick_thickness))) / (2 * math.pi * k_eff_wick * cond_length)
-
- # 3. Vapor thermal resistance (usually neglected)
- # v_vapor = heat_load / (fluid["latent_heat"] * vapor["vapor_density"] * area_vapor)
-
- # # Calculate vapor pressure drop
- # dP = (8 * vapor["vapor_viscosity"] * self.length * v_vapor) / \
- # (self.diameter - 2*(wall_thickness + wick_thickness))**2
-
- # # Calculate vapor thermal resistance
- # R_vapor = (operating_temp * dP) / (fluid["vapor_pressure"] * heat_load)
-
- # Calculate total thermal resistance
- R_total = R_wall_e + R_wall_c + R_wick_e + R_wick_c
-
- return {
- "R_wall_evaporator": R_wall_e,
- "R_wall_condenser": R_wall_c,
- "R_wick_evaporator": R_wick_e,
- "R_wick_condenser": R_wick_c,
- # "R_vapor": R_vapor,
- "R_total": R_total
- }
-
- def calculate_geometrical_parameters(self, wall_thickness=0.001, wick_thickness=0.0005):
- """
- Calculate important geometrical parameters of the heat pipe
-
- Returns:
- dict: Dictionary containing geometrical parameters
- """
- # Inner diameter of the pipe
- d_inner = self.diameter - 2 * wall_thickness
-
- # Vapor space diameter
- d_vapor = d_inner - 2 * wick_thickness
-
- # Cross-sectional areas
- area_vapor = math.pi * d_vapor**2 / 4
- area_wick = math.pi * (d_inner**2 - d_vapor**2) / 4
-
- return {
- "d_inner": d_inner,
- "d_vapor": d_vapor,
- "area_vapor": area_vapor,
- "area_wick": area_wick
- }
-
- def check_operating_limits(self, heat_load, angle=0, wall_thickness=0.001,
- wick_thickness=0.0005, operating_temp=353.15):
- """
- Calculate all operating limits of the heat pipe
-
- Args:
- heat_load (float): Heat load in Watts
- angle (float): Inclination angle in degrees (0 = horizontal)
- wall_thickness (float): Wall thickness in meters
- wick_thickness (float): Wick thickness in meters
- operating_temp (float): Operating temperature in Kelvin
-
- Returns:
- dict: Dictionary containing all operating limits and their values
- """
- # Get fluid properties using IAPWS
- fluid = self.get_water_properties(operating_temp, quality=0) # Liquid properties
- vapor = self.get_water_properties(operating_temp, quality=1) # Vapor properties
-
- wick = self.wick_properties[self.wick_type]
- geom = self.calculate_geometrical_parameters(wall_thickness, wick_thickness)
-
- # 1. Capillary Limit
- # Maximum capillary pressure
- P_cap_max = 2 * fluid["surface_tension"] * math.cos(math.radians(angle)) / wick["pore_radius"]
-
- # Liquid pressure drop
- dP_l = (fluid["liquid_viscosity"] * heat_load * self.length) / \
- (fluid["liquid_density"] * fluid["latent_heat"] * wick["permeability"] * geom["area_wick"])
-
- # Vapor pressure drop (simplified)
- v_vapor = heat_load / (fluid["latent_heat"] * fluid["vapor_density"] * geom["area_vapor"])
- dP_v = (8 * fluid["vapor_viscosity"] * self.length * v_vapor) / geom["d_vapor"]**2
-
- # Gravitational pressure drop
- dP_g = fluid["liquid_density"] * 9.81 * self.length * math.sin(math.radians(angle))
-
- # Capillary limit
- Q_cap = (P_cap_max - dP_g) * fluid["latent_heat"] * wick["permeability"] * geom["area_wick"] / \
- (fluid["liquid_viscosity"] * self.length)
-
- merit = (fluid["liquid_density"] * fluid["surface_tension"] * fluid["latent_heat"]) / fluid["liquid_viscosity"]
- Q_cap_alt = merit * ((wick["permeability"] * geom["area_wick"]) / self.length) * ((2 / wick["pore_radius"]) - \
- (fluid["liquid_density"] * 9.81 * self.length * math.sin(math.radians(angle))) / fluid["surface_tension"])
- print(f"Q_cap_claude: {Q_cap}; Q_sonic_Jupy: {Q_cap_alt}")
-
- # 2. Sonic Limit
- # Sonic velocity at operating temperature
- Q_sonic = fluid["sonic_velocity"] * geom["area_vapor"] * fluid["vapor_density"] * \
- fluid["latent_heat"] / 2
-
- specific_gas_constant_water_vapor = 461.52 # J/kg.K
- Q_sonic_alt = geom["area_vapor"] * fluid["vapor_density"] * fluid["latent_heat"] * \
- ((fluid["gamma0_vapor"] * specific_gas_constant_water_vapor * operating_temp) / (2 * (fluid["gamma0_vapor"] + 1)))**0.5
-
- print(f"Q_sonic_claude: {Q_sonic}; Q_sonic_Jupy: {Q_sonic_alt}")
-
- # 3. Boiling Limit
- # Nucleate boiling limit
- # r_n = wick["particle_diameter"] / 4 # Nucleation site radius (not available for all wick types!!)
- r_n = 3 / 1000000
- k_eff = wick["thermal_conductivity"] * (1 - wick["porosity"]) + \
- fluid["thermal_conductivity"] * wick["porosity"]
-
- # Critical temperature difference for nucleate boiling
- # dT_crit = 2 * fluid["surface_tension"] * fluid["critical_temperature"] / \
- # (fluid["vapor_density"] * fluid["latent_heat"] * r_n)
-
- # Q_boil = 2 * math.pi * self.length * k_eff * dT_crit / \
- # math.log(self.diameter/(self.diameter - 2*wall_thickness))
-
- r_inner = (self.diameter - 2*wall_thickness)/2
- r_vapor = geom["d_vapor"]/2
-
- Q_boil = ((2 * math.pi * self.length * k_eff * operating_temp) / \
- (fluid["latent_heat"] * fluid["vapor_density"] * math.log(r_inner / r_vapor))) * \
- (((2 * fluid["surface_tension"]) / r_n) - P_cap_max)
-
- # 4. Entrainment Limit
- # Weber number criterion
- sigma = fluid["surface_tension"]
- Q_ent = geom["area_vapor"] * math.sqrt(
- sigma * fluid["vapor_density"] / \
- (2 * wick["pore_radius"])
- ) * fluid["latent_heat"]
-
- # 5. Viscous Limit
- # Vapor pressure at operating temperature (simplified calculation)
- P_v = 47.39e3 # Pa at 80°C
- R_gas = 461.5 # Gas constant for water vapor
-
- # Q_visc = (math.pi * geom["d_vapor"]**4 * fluid["vapor_density"] * P_v * fluid["latent_heat"]) / \
- # (32 * fluid["vapor_viscosity"] * self.length * operating_temp)
-
- Q_visc = (math.pi * r_vapor**4 * fluid["vapor_density"] * P_v * fluid["latent_heat"]) / \
- (12 * fluid["vapor_viscosity"] * self.length) #Zohuri
-
- # Determine the most limiting factor
- limits = {
- "capillary_limit": Q_cap,
- "sonic_limit": Q_sonic,
- "boiling_limit": Q_boil,
- "entrainment_limit": Q_ent,
- "viscous_limit": Q_visc
- }
-
- overall_limit = min(limits.values())
- limiting_factor = min(limits.items(), key=lambda x: x[1])[0]
-
- return {
- **limits,
- "overall_limit": overall_limit,
- "limiting_factor": limiting_factor,
- "is_operational": heat_load < overall_limit,
- "operation_margin": (overall_limit - heat_load) / overall_limit * 100 if heat_load < overall_limit else 0
- }
-
-def main():
-
- #TODO can use http://pyromat.org/doc_intro.html to support other working fluids...
- # Create a heat pipe object
- hp = HeatPipe(
- length=0.3, # 300 mm
- diameter=0.015, # 15 mm
- working_fluid="water",
- wick_type="sintered"
- )
-
- # Test water properties at different temperatures
- # print("Water Properties at Different Temperatures:")
- # for temp in [313.15, 333.15, 353.15]: # 40°C, 60°C, 80°C
- # props = hp.get_water_properties(temp, quality=0)
- # print(f"\nTemperature: {temp-273.15}°C")
- # for key, value in props.items():
- # print(f"{key}: {value:.6g}")
-
- #calculate operating limits
- limits = hp.check_operating_limits(heat_load=30, angle=0, wall_thickness=0.003, wick_thickness=0.003, operating_temp=353.15)
- print("\nOperating limits results")
- for key, value in limits.items():
- print(f"{key}: {value}")
-
- # Calculate thermal resistance at 80°C
- results = hp.calculate_thermal_resistance(
- heat_load=100, # 100 W
- evap_length=0.05, # 50 mm
- cond_length=0.10, # 100 mm
- operating_temp=353.15 # 80°C
- )
-
- print("\nThermal Resistance Results:")
- for key, value in results.items():
- print(f"{key}: {value:.6f} K/W")
-
-if __name__ == "__main__":
- main()
\ No newline at end of file