diff --git a/docs/wordlist b/docs/wordlist index 8855e1bb..eaea6fc5 100644 --- a/docs/wordlist +++ b/docs/wordlist @@ -448,6 +448,8 @@ sublist timemax NbSetpsPerPeriod NbSteps +WaveformType +powerloss # electrostatic capacitive diff --git a/femmt/component.py b/femmt/component.py index 09ee2333..b90b6fb2 100644 --- a/femmt/component.py +++ b/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): """ diff --git a/femmt/enumerations.py b/femmt/enumerations.py index af0ad1e3..e5eef459 100644 --- a/femmt/enumerations.py +++ b/femmt/enumerations.py @@ -1,6 +1,7 @@ """Enumeration for FEMMT.""" from enum import IntEnum, Enum + class Verbosity(IntEnum): """State of verbosity.""" @@ -10,6 +11,7 @@ class Verbosity(IntEnum): ToConsole = 2 # Outputs to console ToFile = 3 # Outputs to file + class VisualizationMode(IntEnum): """How simulation results are shown in Gmsh.""" @@ -17,6 +19,7 @@ class VisualizationMode(IntEnum): Post = 2 # Animate after simulation Stream = 3 # Live results during simulation + class WindingTag(IntEnum): """Names of windings.""" @@ -62,6 +65,7 @@ class ComponentType(IntEnum): Transformer = 2 IntegratedTransformer = 3 + class SimulationType(IntEnum): """Sets the simulation type. The static is just to show the fields.""" @@ -69,12 +73,14 @@ class SimulationType(IntEnum): 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" diff --git a/femmt/functions.py b/femmt/functions.py index 02410d3a..0d1b512d 100644 --- a/femmt/functions.py +++ b/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