Updates related to electric module

source_code/conductor_flags.py

@@ -16,15 +16,19 @@
 # Electric conductance is not defined per unit length
 ELECTRIC_CONDUCTANCE_NOT_UNIT_LENGTH = 2
 
+# Self inductance readed from input file conductor_definition.xlsx
+SELF_INDUCTANCE_MODE_0 = 0
 # Analytical self inductance is evaluated according to mode 1.
 SELF_INDUCTANCE_MODE_1 = 1
 # Analytical self inductance is evaluated according to mode 2.
 SELF_INDUCTANCE_MODE_2 = 2
 
+# Constant mutual inductance readed from input file conductor_definition.xlsx
+CONSTANT_INDUCTANCE = 0
 # Flag to evaluate inductance analytically
-ANALYTICAL_INDUCTANCE = 0
+ANALYTICAL_INDUCTANCE = 1
 # Flag to evaluate inductance using an approximation
-APPROXIMATE_INDUCTANCE = 1
+APPROXIMATE_INDUCTANCE = 2
 
 # Flag to solve the electric problem in steady state conditions.
 STATIC_ELECTRIC_SOLVER = 0

source_code/properties_of_materials/rare_earth_123.py

@@ -393,7 +393,7 @@ def critical_current_density_bisection_re123(TT, BB, JOP, TC0M, BC20M, C):
 
     for _, vv in enumerate(JC_ind):
 
-        T_lower = np.array([3.5])
+        T_lower = np.array([0])
         # T_upper = np.array([40.0])
         T_upper = np.array([TC0M])
 

source_code/solid_component.py

@@ -267,6 +267,10 @@ def get_current(self, conductor: object):
                     self.dict_Gauss_pt["op_current_sc"] = self.dict_Gauss_pt[
                         "op_current"
                     ]
+                if conductor.cond_time[-1] == 0:
+                    self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"]
+                    self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"]
+
 
                 if self.flagSpecfield_current == 2:
                     # Add also a logger
@@ -297,6 +301,9 @@ def get_current(self, conductor: object):
                         self.dict_node_pt["op_current_sc"][:-1]
                         + self.dict_node_pt["op_current_sc"][1:]
                     ) / 2.0
+                if conductor.cond_time[-1] == 0:
+                    self.dict_node_pt["current_along"] = self.dict_node_pt["op_current"]
+                    self.dict_Gauss_pt["current_along"] = self.dict_Gauss_pt["op_current"]
 
             elif conductor.inputs["I0_OP_MODE"] == IOP_NOT_DEFINED:
                 # User does not specify a current: set current carrient 
@@ -659,8 +666,7 @@ def initialize_electric_quantities(self, conductor):
         * self.dict_Gauss_pt["current_along"]
         * self.dict_Gauss_pt["delta_voltage_along"]
         * self.dict_Gauss_pt["delta_voltage_along_sum"]
-        * self.dict_Gauss_pt["integral_power_el_res_mod1"]
-        * self.dict_Gauss_pt["integral_power_el_res_mod2"]
+        * self.dict_Gauss_pt["integral_power_el_res"]
         * self.dict_Gauss_pt["linear_power_el_resistance"]
         * self.dict_node_pt["total_linear_power_el_cond"]
         * self.dict_node_pt["total_power_el_cond"]
@@ -675,12 +681,10 @@ def initialize_electric_quantities(self, conductor):
         self.dict_Gauss_pt["current_along"] = np.zeros(n_elem)
         self.dict_Gauss_pt["delta_voltage_along"] = np.zeros(n_elem)
         self.dict_Gauss_pt["delta_voltage_along_sum"] = np.zeros(n_elem)
-        # Introduced to evaluate the integral value Joule power due to electric 
-        # resistance along the cable in the electric module as R * I^2.
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem)
+
         # Introduced to evaluate the integral value Joule power due to electric 
         # resistance along the cable in the electric module as Delta V * I.
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem)
+        self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem)
 
         self.dict_node_pt["total_power_el_cond"] = np.zeros(n_nod)
         # Introduced to evaluate the integral value Joule power due to electric 
@@ -716,12 +720,11 @@ def set_power_array_to_zeros(self, conductor):
         n_elem = conductor.grid_input["NELEMS"]
 
         # Set the following arrays to zeros
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = np.zeros(n_elem)
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = np.zeros(n_elem)
+        self.dict_Gauss_pt["integral_power_el_res"] = np.zeros(n_elem)
         self.dict_node_pt["integral_power_el_cond"] = np.zeros(n_nod)
 
     def get_joule_power_along(self, conductor: object):
-        """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects.
+        """Method that evaluate the contribution to the integral power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects. This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance.
         This method should be called in the electric method, when the transient solution is used and only for current carriers.
 
         Args:
@@ -733,14 +736,14 @@ def get_joule_power_along(self, conductor: object):
 
         # Finalize the evaluation of the integral value of the Joule power due 
         # to electric resistance along the current carrier started in method 
-        # conductor.eval_integral_joule_power. Array integral_power_el_res_mod1 
+        # conductor.eval_integral_joule_power. Array integral_power_el_res 
         # stores the numerator (energy in J) that here is divided by the 
         # thermal hydraulic to get again a power (W), which is further divided 
         # by the length of the discreticazion element corrected by cos(theta) 
         # to get a linear power density (W/m):
         # P_Joule = sum_1^N_em P_{Joule,i} / (Delta_t_TH * Delta_z * costheta)
         integral_j_pow_along = (
-            self.dict_Gauss_pt["integral_power_el_res_mod1"]
+            self.dict_Gauss_pt["integral_power_el_res"]
             / (conductor.grid_features["delta_z"] * self.inputs["COSTETA"]
             * conductor.time_step)
         )
@@ -823,39 +826,33 @@ def get_joule_power_across(self, conductor: object):
             )
 
     def get_joule_power_along_steady(self, conductor: object):
-        """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the electic resistances along the SolidComponent objects.
+        """Method that evaluate the contribution to the total power in the element of Joule power (in W/m) due to the voltage difference (due to electric resistance and inductance) along the SolidComponent objects.
+        This approach was discussed with prof. Zach Hartwig and Dr. Nicolò Riva and is more general and coservative that the evaluation that takes into account only the contribution of the electric resistance.
         This method should be called in the electric method, when the transient solution is used and for the solid component that actually carries a current. It works when the steady state solution for the electric module is computed.
 
         Args:
             conductor (object): ConductorComponent object with all informations to make the calculation.
         """
 
         # Alias
-        el_res = self.dict_Gauss_pt["electric_resistance"]
         current = self.dict_Gauss_pt["current_along"]
         voltage = self.dict_Gauss_pt["delta_voltage_along"]
         d_z_tilde = conductor.grid_features["delta_z"] * self.inputs["COSTETA"]
 
-        # Mode 1: evaluate Joule linear power along the strand in W, due
-        # to electric resistances only for current carriers:
-        # P_along = R_along * I_along ^2
-        self.dict_Gauss_pt["integral_power_el_res_mod1"] = current ** 2 * el_res
-
-        # Mode 2: evaluate Joule linear power along the strand in W, due
+        # Evaluate Joule linear power along the strand in W, due
         # to electric resistances only for current carriers:
         # P_along = Delta_Phi_along * I_along
-        self.dict_Gauss_pt["integral_power_el_res_mod2"] = voltage * current
+        # N.B. this evaluation accounts aslo for the voltage due to the 
+        # inductance and is a conservative an more general approach. Discussed 
+        # with prof. Zach Hartwig and Dr. Nicolò Riva.
+        self.dict_Gauss_pt["integral_power_el_res"] = voltage * current
 
-        # Check equivalence of Mode 1 and Mode 2 (they should be equivalent 
-        # but Mode 2 may give numerical cancellation).
-        if not np.allclose(self.dict_Gauss_pt["integral_power_el_res_mod1"],self.dict_Gauss_pt["integral_power_el_res_mod2"]):
-            warnings.warn("P_Joule = R I^2 != Delta_Phi I. Possible violation of the energy conservation!")
         # Convert W in W/m keping into account the cos(theta).
         # This is independent of the time integration method since at the 
         # initialization (when this function is used) BE, CN and AM4 should 
         # fill only the 0 index column.
         self.dict_Gauss_pt["linear_power_el_resistance"][:, 0] = (
-            self.dict_Gauss_pt["integral_power_el_res_mod1"] / d_z_tilde
+            self.dict_Gauss_pt["integral_power_el_res"] / d_z_tilde
         )
 
     def get_joule_power_across_steady(self, conductor: object):