Merge pull request #166 from upb-lea/MagNet_Model_Implementation

Implementation of MagNet model in FEMMT

Author: gituser789

docs/wordlist

@@ -448,6 +448,8 @@ sublist
 timemax
 NbSetpsPerPeriod
 NbSteps
+WaveformType
+powerloss
 
 # electrostatic
 capacitive

femmt/component.py

@@ -219,6 +219,121 @@ def __init__(self, simulation_type: SimulationType = SimulationType.FreqDomain,
         self.onelab_client = onelab.client(__file__)
         self.simulation_name = simulation_name
 
+    def calc_hystersis_losses_with_MagNet_model_PB_based_on_reluctance(self, peak_magnetizing_current: float = 1, b_wave: WaveformType = WaveformType.Sine,
+                                                                       custom_b_wave: np.ndarray = None) -> float:
+        """
+        Calculate the hysteresis losses with the MagNet model of Paderborn University based on a reluctance model.
+
+        :param peak_magnetizing_current: peak value of the magnetizing current
+        :type peak_magnetizing_current: peak value of the magnetizing current
+        :param b_wave: type of the waveform of the magnetic flux density/magnetizing current
+        :type b_wave: WaveformType
+        :param custom_b_wave: normalized signal of the magnetic flux density/magnetizing current (normalized -> amplitude of 1)
+        :type custom_b_wave: np.ndarray
+        :return: hysteresis losses
+        :rtype: float
+        """
+        total_reluctance = self.calculate_core_reluctance()[0] + self.air_gaps_reluctance()[0]
+        center_leg_area = np.pi*(self.core.geometry.core_inner_diameter/2)**2
+
+        turns = ff.get_number_of_turns_of_winding(winding_windows=self.winding_windows, windings=self.windings, winding_number=0)
+        inductance = (turns**2)/total_reluctance
+        b_field_peak = inductance*peak_magnetizing_current/center_leg_area/turns
+
+        # multiply the normalized signal of the magnetic flux density with the peak value of the magnetic flux density
+        if b_wave is WaveformType.Sine:
+            b_wave = np.sin(np.linspace(0, 2*np.pi, 1024)) * b_field_peak
+        elif b_wave is WaveformType.Custom:
+            b_wave = custom_b_wave * b_field_peak
+
+        power_loss = ff.calc_power_loss_from_MagNet_model_PB(material_name=self.core.material.material, b_wave=b_wave, frequency=self.frequency,
+                                                             temperature=self.core.material.temperature) * self.calculate_core_volume()
+        return power_loss
+
+    def calc_hystersis_losses_with_MagNet_model_PB_based_on_mesh_results(self, b_wave: WaveformType = WaveformType.Sine,
+                                                                         custom_b_wave: np.ndarray = None) -> float:
+        """
+        Calculate the hysteresis losses with the MagNet model of Paderborn University based on FEM-simulation results.
+
+        - uses the local magnetic flux density in each mesh cell
+        - arbitrary waveforms are possible. But FEM simulation assumes a sinusoidal waveform and the peak of the magnetic flux density result of each mesh cell
+          is used to linearly scale the arbitrary waveform to its max. value
+        - for each mesh cell, the magnet model is applied
+
+        :param b_wave: type of the waveform of the magnetic flux density/magnetizing current
+        :type b_wave: WaveformType
+        :param custom_b_wave: normalized signal of the magnetic flux density/magnetizing current (normalized -> amplitude of 1)
+        :type custom_b_wave: np.ndarray
+        :return: dict containing the hysteresis losses
+        """
+        flux, volume = self.calc_magnetic_flux_density_based_on_simulation_results()
+
+        # multiply the normalized signal of the magnetic flux density with the peak value of the magnetic flux density
+        if b_wave is WaveformType.Sine:
+            b_wave = np.array([np.sin(np.linspace(0, 2*np.pi, 1024)) * b for b in flux])
+        elif b_wave is WaveformType.Custom:
+            b_wave = np.array([custom_b_wave * b for b in flux])
+
+        hyst_losses_density = ff.calc_power_loss_from_MagNet_model_PB(material_name=self.core.material.material, b_wave=b_wave,
+                                                                      frequency=np.array([self.frequency] * b_wave.shape[0]),
+                                                                      temperature=np.array([self.core.material.temperature] * b_wave.shape[0]))
+        power_loss = float(np.dot(hyst_losses_density, volume))
+        return power_loss
+
+    def calc_magnetic_flux_density_based_on_simulation_results(self) -> tuple[np.ndarray, np.ndarray]:
+        """Calculate the magnetic flux density and volume for every mesh cell of the core."""
+        with open(os.path.join(self.file_data.e_m_fields_folder_path, "Magb.pos"), 'r') as file:
+            content = file.read()  # type of content: str
+
+        # find both occurrence of "Magnitude B-Field / T", Elements and Nodes in Magb.pos
+        indexes_B_Field = np.array([m.start() for m in re.finditer('"Magnitude B-Field / T"', content)])
+        indexes_Elements = np.array([m.start() for m in re.finditer('Elements', content)])
+        indexes_Coords = np.array([m.start() for m in re.finditer('Nodes', content)])
+        # filter magnetic flux density out of Magb.pos(data is between the both occurrences of "Magnitude B-Field / T") and put data into list
+        data_B_Field = np.array([float(x) for x in content[indexes_B_Field[0]:indexes_B_Field[1]].split()[11:-3]])
+        # filter elements(ID of triangles and nodes) out of Magb.pos(data is between the both occurrences of "Elements") and put data into list
+        data_Elements = np.array([int(x) for x in content[indexes_Elements[0]:indexes_Elements[1]].split()[2:-1]])
+        # filter coordinates of nodes out of Magb.pos(data is between the both occurrences of "Elements") and put data into list
+        data_Coords = np.array([float(x) for x in content[indexes_Coords[0]:indexes_Coords[1]].split()[2:-1]])
+        # divide long list into lists for every triangle/mesh-cell
+        chunks_B_Field = [data_B_Field[x:x + 5] for x in range(0, data_B_Field.shape[0], 5)]
+        chunks_Elements = np.array([data_Elements[x:x + 8] for x in range(0, data_Elements.shape[0], 8)])
+        chunks_Coords = np.array([data_Coords[x:x + 4] for x in range(0, data_Coords.shape[0], 4)])
+
+        data_triangle = np.concatenate((chunks_Elements, np.array(chunks_B_Field)), axis=1)
+        # rearrange arrays for data that is needed for calculation of mesh cell volume
+        Triangle_Core = [[x[0], x[5], x[6], x[7], np.mean(x[10:])] for x in data_triangle if x[3] // 10000 == 12]
+        # x[0] = ID of mesh cell(triangle), x[5] = ID of first node of triangle, x[5] = ID of second node of triangle, x[5] = ID of third node of triangle,
+        # np.mean(x[10:]) = avr. of magnetic flux density of all three nodes
+
+        xyz = {}
+        for x in chunks_Coords:  # rearrange coordinates of nodes. These are the coordinates from the nodes of the whole mesh.
+            xyz[str(x[0])] = np.array([x[1], x[2], x[3]])  # key of dict is ID of node
+
+        for index, TriangleNode in enumerate(Triangle_Core):
+            xi = xyz[str(TriangleNode[1:4][0])]  # first node of triangle (x1, y1, z1)
+            xj = xyz[str(TriangleNode[1:4][1])]  # second node of triangle (x2, y2, z2)
+            xk = xyz[str(TriangleNode[1:4][2])]  # third node of triangle (x3, y3, z3)
+
+            radius = np.mean([xi[0], xj[0], xk[0]])  # radius -> average of x-values np.mean(x1, x2, x3)
+
+            a = np.sqrt(np.dot(xi - xj, xi - xj))  # distance between first and second node
+            b = np.sqrt(np.dot(xi - xk, xi - xk))  # distance between first and third node
+            c = np.sqrt(np.dot(xj - xk, xj - xk))  # distance between second and third node
+
+            # Heron's formula
+            s = (a + b + c) / 2  # semiperimeter of mesh cell (triangle)
+            dA = np.sqrt(s * (s - a) * (s - b) * (s - c))  # area of mesh cell
+
+            volume = 2 * np.pi * radius * dA
+
+            Triangle_Core[index].append(volume)
+
+        flux = np.array([x[-2] for x in Triangle_Core])
+        volume = np.array([x[-1] for x in Triangle_Core])
+
+        return flux, volume
+
     def update_mesh_accuracies(self, mesh_accuracy_core: float, mesh_accuracy_window: float,
                                mesh_accuracy_conductor, mesh_accuracy_air_gaps: float):
         """

femmt/enumerations.py

@@ -1,6 +1,7 @@
 """Enumeration for FEMMT."""
 from enum import IntEnum, Enum
 
+
 class Verbosity(IntEnum):
     """State of verbosity."""
 
@@ -10,13 +11,15 @@ class Verbosity(IntEnum):
     ToConsole = 2  # Outputs to console
     ToFile = 3  # Outputs to file
 
+
 class VisualizationMode(IntEnum):
     """How simulation results are shown in Gmsh."""
 
     Final = 1   # Only final static result
     Post = 2    # Animate after simulation
     Stream = 3  # Live results during simulation
 
+
 class WindingTag(IntEnum):
     """Names of windings."""
 
@@ -62,19 +65,22 @@ class ComponentType(IntEnum):
     Transformer = 2
     IntegratedTransformer = 3
 
+
 class SimulationType(IntEnum):
     """Sets the simulation type. The static is just to show the fields."""
 
     FreqDomain = 1
     TimeDomain = 2
     ElectroStatic = 3
 
+
 class CoreType(IntEnum):
     """Sets the core type for the whole simulation. Needs to be given to the MagneticComponent on creation."""
 
     Single = 1  # one axisymmetric core
     Stacked = 2  # one and a half cores
 
+
 class AirGapMethod(IntEnum):
     """Sets the method how the air gap position (vertical) is set.
 
@@ -190,6 +196,7 @@ class ConductorArrangement(IntEnum):
     turns of the previous line. First drawn in y-direction then x-direction.
     """
 
+
 class Align(IntEnum):
     """Specifies the distribution direction for conductors when starting from the center. This can be done for having single windings in vww."""
 
@@ -202,6 +209,7 @@ class Align(IntEnum):
     CenterOnHorizontalAxis = 3
     """Conductors are centered across the middle line of the horizontal axis."""
 
+
 class ConductorDistribution(IntEnum):
     """Defines specific strategies for placing conductors starting from the peripheral (edges) of the virtual winding window."""
 
@@ -229,6 +237,7 @@ class ConductorDistribution(IntEnum):
     HorizontalLeftward_VerticalDownward = 8
     """Places conductors horizontally leftward first, then moves vertically downward for the next set with consistent direction."""
 
+
 class FoilVerticalDistribution(IntEnum):
     """Defines specific strategies for placing rectangular foil vertical conductors starting from the peripheral (edges) of the virtual winding window."""
 
@@ -238,6 +247,7 @@ class FoilVerticalDistribution(IntEnum):
     HorizontalLeftward = 2
     """Moves horizontally leftward for the next set with consistent direction."""
 
+
 class FoilHorizontalDistribution(IntEnum):
     """Defines specific strategies for placing rectangular foil horizontal conductors starting from the peripheral (edges) of the virtual winding window."""
 
@@ -247,6 +257,7 @@ class FoilHorizontalDistribution(IntEnum):
     VerticalDownward = 2
     """Moves vertically downward for the next set with consistent direction."""
 
+
 class CenterTappedInterleavingType(IntEnum):
     """Contains different interleaving types for the center tapped transformer."""
 
@@ -310,6 +321,7 @@ class WrapParaType(IntEnum):
     The foils will have a given width and thickness. The virtual winding window may not be fully occupied..
     """
 
+
 class ConductorMaterial(IntEnum):
     """Sets the conductivity of the conductor."""
 
@@ -346,3 +358,10 @@ class CoreMaterialType(str, Enum):
 
     Linear = "linear"
     Imported = "imported"
+
+
+class WaveformType(str, Enum):
+    """Types of waveforms."""
+
+    Sine = "sine"
+    Custom = "custom"

femmt/functions.py

@@ -8,15 +8,17 @@
 import warnings
 from scipy.integrate import quad
 import logging
-
 # Third parry libraries
 import gmsh
 from matplotlib import pyplot as plt
 import numpy.typing as npt
 import numpy as np
+from materialdatabase.meta.data_enums import Material
+import magnethub as mh
 import schemdraw
 import schemdraw.elements as elm
 
+
 # Local libraries
 from femmt.constants import *
 from femmt.enumerations import ConductorType
@@ -2765,5 +2767,28 @@ def compare_excel_files(femmt_excel_path: str, femm_excel_path: str, comparison_
                 worksheet.set_column('A:Z', 35)
 
 
+def calc_power_loss_from_MagNet_model_PB(material_name: Material, b_wave: np.ndarray, frequency: float, temperature: float) -> float | None:
+    """
+    Calculate the power loss density with the help of the trained neural network of the MagNet Challenge 2023.
+
+    :param material_name: Name of the material
+    :type material_name: Material
+    :param b_wave: array containing the shape of the magnetic flux density in time domain in T
+    :type b_wave: np.ndarray
+    :param frequency: value of the frequency in Hz
+    :type frequency: float
+    :param temperature: value of the temperature in °C
+    :type temperature: float
+    :return: powerloss in W/m^3
+    """
+    if material_name.value in ["3C90", "3C92", "3C94", "3C95", "3E6", "3F4", "77", "78", "79", "ML95S", "T37", "N27", "N30", "N49", "N87"]:
+        mdl = mh.loss.LossModel(material=material_name.value, team="paderborn")
+        power_loss_density, h_wave = mdl(b_wave, frequency, temperature)
+        return power_loss_density
+    else:
+        logger.warning("Material" + str(material_name.value) + " not supported by MagNet Models.")
+        return None
+
+
 if __name__ == '__main__':
     pass