reluctance calculation for single simulation, need to be reviewd

femmt/component.py

@@ -630,6 +630,41 @@ def check_model_mqs_condition(self) -> None:
 
     #  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -   -  -  -  -  -  -  -  -  -  -  -
     # Miscellaneous
+    def calculate_core_cross_sectional_area(self):
+        """
+        Calculate the cross-sectional area of the core.
+
+        :return: Cross-sectional area of the core.
+        :rtype: float
+        """
+        total_core_reluctance, total_length = self.calculate_core_reluctance()
+        core_volume = self.calculate_core_volume()
+        cross_sectional_area = core_volume / total_length
+        return cross_sectional_area
+
+    def calculate_air_gap_area(self):
+        """
+        Calculate the total cross-sectional area of all air gaps.
+
+        :return: Total air gap area.
+        :rtype: float
+        """
+        air_gap_area = 0
+
+        for leg_position, _, height in self.air_gaps.midpoints:
+            if leg_position == AirGapLegPosition.LeftLeg.value or leg_position == AirGapLegPosition.RightLeg.value:
+                # For left and right leg, same width as core outer radius minus inner radius
+                width = self.core.r_outer - self.core.r_inner
+            elif leg_position == AirGapLegPosition.CenterLeg.value:
+                # For center leg, the width as core inner diameter divided by 2
+                width = self.core.core_inner_diameter / 2
+            else:
+                raise Exception(f"Invalid leg position tag {leg_position} used for an air gap.")
+
+            air_gap_area += np.pi * (width ** 2)
+
+        return air_gap_area
+
     def calculate_core_volume_with_air(self) -> float:
         """Calculate the volume of the core including air.
 
@@ -1027,6 +1062,7 @@ def excitation(self, frequency: float, amplitude_list: List, phase_deg_list: Lis
         """
         # negative currents are not allowed and lead to wrong simulation results. Check for this.
         # this message appears after meshing but before simulation
+
         for amplitude in amplitude_list:
             if amplitude < 0:
                 raise ValueError(
@@ -1104,6 +1140,9 @@ def excitation(self, frequency: float, amplitude_list: List, phase_deg_list: Lis
             for num in range(len(self.windings)):
                 self.red_freq[num] = 0
 
+        # check the core saturation
+        self.reluctance_model_pre_check()
+
     def excitation_time_domain(self, current_list: List[List[float]], time_list: List[float],
                                number_of_periods: int, ex_type: str = 'current',
                                plot_interpolation: bool = False, imposed_red_f=0):
@@ -1974,6 +2013,185 @@ def factor_triangular_hysteresis_loss_iGSE(D, alpha):
                               center_tapped_study_excitation["linear_losses"]["current_phases_deg"],
                               inductance_dict=inductance_dict, core_hyst_loss=p_hyst_core_parts)
 
+    #  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
+    # Reluctance
+    def calculate_core_reluctance(self):
+        """Calculate the core reluctance."""
+        length = []
+        reluctance = []
+        if self.air_gaps.midpoints:
+            # # Sorting air gaps from lower to upper
+            sorted_midpoints = sorted(self.air_gaps.midpoints, key=lambda x: x[1])
+            # Finding position of first airgap
+            bottommost_airgap_position = sorted_midpoints[0][1]
+            bottommost_airgap_height = sorted_midpoints[0][2]
+            # Finding position of last airgap
+            topmost_airgap_position = sorted_midpoints[-1][1]
+            topmost_airgap_height = sorted_midpoints[-1][2]
+
+        def get_width(part_number):
+            """
+            If there is a stray path, calculate width based on its starting index and part number.
+
+            part_number is the core_part_i+2; means that if the start_index is 0, the stray path is in core_part_2
+            if the start_index is 1, the stray path is in core_part_3 and so on
+            """
+            if self.stray_path and part_number == self.stray_path.start_index + 2:
+                return self.stray_path.length
+            return self.core.core_inner_diameter / 2
+
+        if self.core.core_type == CoreType.Single:
+            # subpart1: bottom left subpart
+            if self.air_gaps.midpoints:
+                subpart_1_1_length = bottommost_airgap_position + self.core.window_h / 2 - bottommost_airgap_height / 2
+            else:
+                subpart_1_1_length = self.core.window_h / 2
+            length.append(subpart_1_1_length)
+            # subpart1_1_reluctance = fr.r_core_tablet(subpart_1_1_length, self.core.core_inner_diameter / 2, self.core.mu_r_abs, self.core.core_inner_diameter)
+            subpart1_1_reluctance = fr.r_core_tablet_2(subpart_1_1_length, self.core.core_inner_diameter / 2, self.core.mu_r_abs)
+            reluctance.append(subpart1_1_reluctance)
+
+            # subpart2: top left subpart
+            if self.air_gaps.midpoints:
+                subpart_1_2_length = self.core.window_h / 2 - topmost_airgap_position - topmost_airgap_height / 2
+            else:
+                subpart_1_2_length = self.core.window_h / 2
+            # subpart1_2_reluctance = fr.r_core_tablet(subpart_1_2_length, self.core.core_inner_diameter / 2, self.core.mu_r_abs, self.core.core_inner_diameter)
+            length.append(subpart_1_2_length)
+            subpart1_2_reluctance = fr.r_core_tablet_2(subpart_1_2_length, self.core.core_inner_diameter / 2, self.core.mu_r_abs)
+            reluctance.append(subpart1_2_reluctance)
+
+            # subpart3: bottom and top mid-subpart. It has the inner, outer corners and winding window section
+            # The area of the inner, and outer corners (top and bottom) is approximated by taking the mean cross-sectional area
+            # The length over the area of the winding window section will be approximated to log(r_inner/core_inner_diameter/2) / 2 *pi * core_inner_diameter/4
+            # This is taken from Appendix B of book "E. C. Snelling. Soft Ferrites, Properties and Applications. 2nd edition. Butterworths, 1988"
+            # inner corners
+            s_1 = (self.core.core_inner_diameter / 2) - (self.core.core_inner_diameter / (2 * np.sqrt(2)))
+            length_inner = (np.pi / 4) * (s_1 + (self.core.core_inner_diameter / 8))
+            inner_reluctance = fr.r_core_round(self.core.core_inner_diameter, length_inner, self.core.mu_r_abs) * 2
+            # outer corners
+            s_2 = np.sqrt(((self.core.r_inner ** 2) + (self.core.r_outer ** 2)) / 2) - self.core.r_inner
+            length_outer = (np.pi / 4) * (s_2 + (self.core.core_inner_diameter / 8))
+            outer_reluctance = fr.r_core_round(self.core.core_inner_diameter, length_outer, self.core.mu_r_abs) * 2
+            # corners reluctance
+            corner_reluctance = inner_reluctance + outer_reluctance
+            # winding window
+            length_window = self.core.window_w
+            window_reluctance = (fr.r_core_top_bot_radiant
+                                 (self.core.core_inner_diameter, self.core.window_w, self.core.mu_r_abs, self.core.core_inner_diameter / 4) * 2)
+            subpart1_3_reluctance = corner_reluctance + window_reluctance
+            # total reluctance
+            reluctance.append(subpart1_3_reluctance)
+            # total length
+            subpart_1_3_length = length_inner + length_outer + length_window
+            length.append(subpart_1_3_length)
+
+            # subpart 4: right subpart
+            subpart_1_4_length = self.core.window_h
+            subpart_1_4_radius_eff = np.sqrt(self.core.r_outer ** 2 - self.core.r_inner ** 2)
+            subpart1_4_reluctance = fr.r_core_tablet_2(subpart_1_4_length, subpart_1_4_radius_eff, self.core.mu_r_abs)
+            # subpart1_4_reluctance = fr.r_core_tablet(subpart_1_4_length, self.core.r_outer, self.core.mu_r_abs, self.core.r_inner)
+            reluctance.append(subpart1_4_reluctance)
+            length.append(subpart_1_4_length)
+            # Find the height of subpart 1 and 2 when more than one airgap is existed
+            for i in range(len(sorted_midpoints) - 1):
+                air_gap_1_position = sorted_midpoints[i][1]
+                air_gap_1_height = sorted_midpoints[i][2]
+                air_gap_2_position = sorted_midpoints[i + 1][1]
+                air_gap_2_height = sorted_midpoints[i + 1][2]
+                # calculate the height based on airgap positions and heights, and the width
+                subpart_length = air_gap_2_position - air_gap_2_height / 2 - (air_gap_1_position + air_gap_1_height / 2)
+                subpart_width = get_width(i + 2)
+                # calculate the reluctance
+                subpart_reluctance = fr.r_core_tablet_2(subpart_length, subpart_width, self.core.mu_r_abs)
+                reluctance.append(subpart_reluctance)
+
+            total_core_reluctance = np.sum(reluctance)
+            total_length = np.sum(length)
+            return total_core_reluctance, total_length
+            # return np.sum(reluctance), np.sum(length)
+
+    def airgaps_reluctance(self):
+        """
+        Calculate the air-gap reluctance for a single simulation with single/multiple air-gaps.
+
+        Method is according to the following paper:
+        ["A Novel Approach for 3D Air Gap Reluctance Calculations" - J. Mühlethaler, J.W. Kolar, A. Ecklebe]
+
+        Its calculation for multiple air-gap is based on superposition.
+        That is, multiple air-gap's reluctance is calculated by taking one at a time and then adding them together
+        (like in superposition theorem)
+        """
+        # List to store reluctance for each air gap
+        air_gap_reluctances = []
+
+        for air_gap in self.air_gaps.midpoints:
+            position = air_gap[1]
+            height = air_gap[2]
+
+            core_height_upper = self.core.window_h - position - (height / 2)
+            core_height_lower = position - (height / 2)
+
+            core_height_upper = max(core_height_upper, 0)
+            core_height_lower = max(core_height_lower, 0)
+
+            # Ensure core heights are not zero
+            if core_height_upper == 0 and core_height_lower == 0:
+                raise ValueError("Both core_height_upper and core_height_lower cannot be zero simultaneously")
+
+            # Calculate reluctance based on whether core heights are zero
+            if core_height_upper == 0 or core_height_lower == 0:
+                reluctance = fr.r_air_gap_round_inf(height, self.core.core_inner_diameter, core_height_lower + core_height_upper)
+            else:
+                reluctance = fr.r_air_gap_round_round(height, self.core.core_inner_diameter, core_height_upper, core_height_lower)
+
+            air_gap_reluctances.append(reluctance)
+
+        return air_gap_reluctances
+
+    def reluctance_model_pre_check(self, saturation_threshold: float = 0.7):
+        """
+        Check for possible saturation and consistency of measurement data with simulation results.
+
+        :param saturation_threshold: Threshold for saturation (default is 70%).
+        """
+        # Reluctances
+        total_core_reluctance, total_length = self.calculate_core_reluctance()
+        airgap_reluctance = self.airgaps_reluctance()
+
+        # Calculate the total reluctance
+        reluctance = total_core_reluctance + airgap_reluctance
+
+        # Calculate MMF (Magnetomotive Force)
+        total_mmf = 0
+        for winding_number in range(len(self.windings)):
+            turns = ff.get_number_of_turns_of_winding(winding_windows=self.winding_windows, windings=self.windings,
+                                                      winding_number=winding_number)
+            total_mmf += self.current[0] * turns
+
+        # Calculate Flux (Φ = MMF / Reluctance)
+        flux = total_mmf / reluctance
+        # Area
+        core_area = self.calculate_core_cross_sectional_area()
+        airgap_area = self.calculate_air_gap_area()
+        # Magnetic flux density
+        b_field = flux / ( core_area + airgap_area)
+
+        # Get saturation flux density from material database
+        database = mdb.MaterialDatabase()
+        saturation_flux_density = database.get_saturation_flux_density(self.core.material)
+
+        # Check for saturation and raise error if B-field exceeds threshold
+        for b in b_field:
+            if b > saturation_threshold * saturation_flux_density:
+                raise ValueError(
+                    f"Core saturation detected! B-field ({b} T) exceeds {saturation_threshold * 100}% of the saturation flux density"
+                    f" ({saturation_flux_density} T).")
+        print(b_field)
+        print(flux)
+        print(reluctance)
+        print("Reluctance model pre-check passed.")
+
     #  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
     # Post-Processing
     def get_inductances(self, I0: float, op_frequency: float = 0, skin_mesh_factor: float = 1,

femmt/functions_reluctance.py

@@ -745,6 +745,16 @@ def r_core_tablet(tablet_height, tablet_radius, mu_r_abs, core_inner_diameter):
     """
     return np.log(tablet_radius / (core_inner_diameter / 2)) / (2 * np.pi * mu_0 * mu_r_abs * tablet_height)
 
+def r_core_tablet_2(tablet_height, tablet_radius, mu_r_abs):
+    """
+    Calculate the magnetic resistance of the core tablet.
+
+    :param tablet_height: tablet height
+    :param tablet_radius: tablet radius
+    :param mu_r_abs: relative permeability (mu_r) of the core material from datasheet
+    """
+    return tablet_height / (np.pi * mu_0 * mu_r_abs * tablet_radius ** 2)
+
 
 def r_core_top_bot_radiant(core_inner_diameter, window_w, mu_r_abs, core_top_bot_height):
     """