33from scipy .integrate import solve_ivp
44
55# Constants
6- stefan_boltzmann_constant = 5.67e-8
7- initial_temperature = 293 # make sure it is in kelvins
8- k_steel = 45 # Thermal conductivity of steel in W/mK
6+ stefan_boltzmann_constant = 5.67e-8 # (W/m^2K^4)
7+ initial_temperature = 293 # Initial temperature of the plate sensor (K)
98
109# User-defined variables
11- emission_time = 20 * 60 # 1 minute in seconds
12- wall_thickness = 0.001
13- element_width = 0.05
14- element_height = 0.10
15- emissivity = 0.94
16- steel_density = 8050
17- cp_steel = 490
18- convective_hfc_exposed = 10
19- convective_hfc_unexposed = 8
20- #Initializing grid heat flux for the forward model
10+ emission_time = 20 * 60 # Exposure time (s)
11+ wall_thickness = 0.001 # Thickness of the plate sensor (m)
12+ element_width = 0.05 # Size of a discrete element in the x dimension (m)
13+ element_height = 0.10 # Size of a discrete element in the y dimension (m)
14+ emissivity = 0.94 # Emissivity of the plate sensor
15+ steel_density = 8050 # Density of the plate sensor material (kg/m^3)
16+ cp_steel = 490 # Specific heat capacity of the plate sensor material (J/kgK)
17+ k_steel = 45 # Thermal conductivity of the plate sensor material (W/mK)
18+ convective_hfc_exposed = 10 # Convective heat transfer coefficient on the exposed side of the plate sensor (W/m^2K)
19+ convective_hfc_unexposed = 8 # Convective heat transfer coefficient on the unexposed side of the plate sensor (W/m^2K)
20+
21+ # Initializing grid heat flux for the forward model
2122grid_size = 3 # 3x3 grid --> the simplest 2D case possible
2223central_heat_flux = 50000 # Central cell
23- surrounding_heat_flux = 10000 # neighboring cells
24+ surrounding_heat_flux = 10000 # neighboring cells
2425heat_flux_grid = np .full ((grid_size , grid_size ), surrounding_heat_flux )
2526heat_flux_grid [1 , 1 ] = central_heat_flux # Set central cell heat flux
2627
27- # energy blance for the 2D grid
28+ # Apply conservation of energy over the 2D grid
2829def energy_balance_2d (t , T_flat ):
2930 T = T_flat .reshape ((grid_size , grid_size ))
3031 dTdt = np .zeros_like (T )
@@ -47,7 +48,7 @@ def energy_balance_2d(t, T_flat):
4748 q_conv_exp = convective_hfc_exposed * area * (T [i , j ] - initial_temperature )
4849 q_conv_unexp = convective_hfc_unexposed * area * (T [i , j ] - initial_temperature )
4950 q_rad = 2 * emissivity * stefan_boltzmann_constant * area * (T [i , j ]** 4 - initial_temperature ** 4 )
50- print (q_cond )
51+ # print(q_cond)
5152 # Energy balance
5253 dTdt [i , j ] = (emissivity * q_in - q_conv_exp - q_conv_unexp - q_rad + q_cond ) / (area * cp_steel * steel_density * wall_thickness )
5354
@@ -61,14 +62,15 @@ def energy_balance_2d(t, T_flat):
6162
6263# Generate time steps
6364t_eval = np .linspace (t_span [0 ], t_span [1 ],6000 ) # Creates an array from 0 to emission_time with a specified time step
64- print (t_eval )
65+ # print(t_eval)
66+
6567# Solve the 2D differential equation with specified time steps
6668solution = solve_ivp (energy_balance_2d , t_span , T0_flat , method = 'RK45' , t_eval = t_eval , dense_output = True )
67-
6869num_time_steps = len (solution .t )
69- print ((solution .y ).shape )
70+ # print((solution.y).shape)
7071T_values_2d = solution .y .reshape ((grid_size , grid_size , num_time_steps ))
71- print (T_values_2d .shape )
72+ # print(T_values_2d.shape)
73+
7274# Plotting temperature distribution over time
7375fig , ((ax1 , ax2 , ax3 ), (ax4 , ax5 , ax6 ), (ax7 , ax8 , ax9 )) = plt .subplots (nrows = 3 , ncols = 3 , figsize = (14 , 8 ) )
7476axes = [ax1 , ax2 , ax3 , ax4 , ax5 , ax6 , ax7 , ax8 , ax9 ]
@@ -90,7 +92,7 @@ def energy_balance_2d(t, T_flat):
9092# IHT Model Function
9193def inverse_heat_transfer (T_values_2d , solution , k_steel , cp_steel , steel_density , element_width , element_height , wall_thickness , convective_hfc_exposed , convective_hfc_unexposed , emissivity ):
9294 num_time_steps = T_values_2d .shape [2 ]
93- print (num_time_steps )
95+ # print(num_time_steps)
9496 estimated_flux = np .zeros ((grid_size , grid_size , num_time_steps ))
9597
9698 temp_grad = np .gradient (T_values_2d , axis = 2 ) / np .gradient (solution .t )
@@ -112,7 +114,6 @@ def inverse_heat_transfer(T_values_2d, solution, k_steel, cp_steel, steel_densit
112114 q_cond = np .sum (grad_T ) * wall_thickness
113115 q_conv = (convective_hfc_exposed + convective_hfc_unexposed ) * area * (T_values_2d [i , j , time_step ] - initial_temperature )
114116 q_rad = 2 * emissivity * stefan_boltzmann_constant * area * (T_values_2d [i , j , time_step ]** 4 - initial_temperature ** 4 )
115-
116117 q_storage = cp_steel * steel_density * area * wall_thickness * temp_grad [i , j , time_step ]
117118
118119 # Calculate the estimated incident heat flux
@@ -145,6 +146,6 @@ def inverse_heat_transfer(T_values_2d, solution, k_steel, cp_steel, steel_densit
145146 axx .legend (loc = 'upper right' )
146147 # ax.grid(True)
147148 k += 1
148- print (T_values_2d [1 , 1 , :])
149+ # print(T_values_2d[1, 1, :])
149150plt .tight_layout ()
150151plt .show ()
0 commit comments