diff --git a/README.md b/README.md
index fb1ca0faeef37f456e16a52a4588b93a55b4a437..70a6251a11366c63556442e9f96763dc95a2e5db 100644
--- a/README.md
+++ b/README.md
@@ -10,8 +10,9 @@ Heatpipepy is a library meant to simplify estimation of operating limits and the
To start, clone this repo on your local machine and create a heat pipe object, like the example below:
```
-import heatpipepy
- hp = HeatPipe(
+#run this in the directory where this repo is imported
+import heatpipe as hp
+ my_heatpipe = hp.HeatPipe(
length=0.3,
diameter=0.015
working_fluid="water",
@@ -20,7 +21,7 @@ import heatpipepy
```
Then, the operating limits of the heat pipe can be checked by calling for "check_operating_limits()" method:
```
-limits = hp.check_operating_limits(
+limits = my_heatpipe.check_operating_limits(
heat_load=30,
angle=0,
wall_thickness=0.003,
@@ -43,7 +44,7 @@ It displays results as follows:
The thermal resistance of the heat pipe can be evaluated with the "calculate_thermal_resistance" method:
```
-results = hp.calculate_thermal_resistance(
+results = my_heatpipe.calculate_thermal_resistance(
heat_load=100,
evap_length=0.05,
cond_length=0.10,
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/heatpipe.py b/heatpipe.py
new file mode 100644
index 0000000000000000000000000000000000000000..1234e77d14e1c668c2a1f00ce22ec94049e1e800
--- /dev/null
+++ b/heatpipe.py
@@ -0,0 +1,431 @@
+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")
+ print(limits)
+ # 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:")
+ print(results)
+ # for key, value in results.items():
+ # print(f"{key}: {value:.6f} K/W")
+
+if __name__ == "__main__":
+ main()
\ No newline at end of file