|
Here's a much better simulation of the ac master magnet. Screenshot, animation, and the phase shift of the simulation is attached below. I'm slowly making progress and I am now in the ballpark. Below is the Python code to generate the neccessary files to simulate and create the animation in ParaView.
# Python Code - Master Magnet Simulation # 1. CRITICAL: Place these at the absolute top to bypass MPI resource errors import os os.environ["OMPI_MCA_btl"] = "self" os.environ["OMPI_MCA_pml"] = "ob1" os.environ["PMIX_MCA_gds"] = "hash" os.environ['HDF5_DISABLE_VERSION_CHECK'] = '2'
import gmsh import meshio import numpy as np import matplotlib matplotlib.use('Agg') # Stable for terminal-only execution import matplotlib.pyplot as plt import shutil from dolfin import *
# Tell FEniCS to stay single-threaded to avoid MPI segmenting errors #parameters["num_threads"] = 1
# === 1. GEOMETRY (GMSH) === def generate_custom_mesh(): gmsh.initialize() gmsh.option.setNumber("General.Terminal", 0) gmsh.model.add("ShadedPole") w, h, t, rw, rh = 40, 40, 10, 8, 4 p1 = gmsh.model.geo.addPoint(-w/2, -h/2, 0) p2 = gmsh.model.geo.addPoint(w/2, -h/2, 0) p3 = gmsh.model.geo.addPoint(w/2, h/2, 0) p4 = gmsh.model.geo.addPoint(-w/2, h/2, 0) p5 = gmsh.model.geo.addPoint(-w/2+t, -h/2+t, 0) p6 = gmsh.model.geo.addPoint(w/2, -h/2+t, 0) p7 = gmsh.model.geo.addPoint(w/2, h/2-t, 0) p8 = gmsh.model.geo.addPoint(-w/2+t, h/2-t, 0) p9 = gmsh.model.geo.addPoint(w/2-rw, h/2, 0) p10 = gmsh.model.geo.addPoint(w/2-rw, h/2-rh, 0) p11 = gmsh.model.geo.addPoint(w/2, h/2-rh, 0)
l1 = gmsh.model.geo.addLine(p1, p2); l2 = gmsh.model.geo.addLine(p2, p6) l3 = gmsh.model.geo.addLine(p6, p5); l4 = gmsh.model.geo.addLine(p5, p8) l5 = gmsh.model.geo.addLine(p8, p7); l6 = gmsh.model.geo.addLine(p7, p11) l7 = gmsh.model.geo.addLine(p11, p10); l8 = gmsh.model.geo.addLine(p10, p9) l9 = gmsh.model.geo.addLine(p9, p4); l10 = gmsh.model.geo.addLine(p4, p1) l11 = gmsh.model.geo.addLine(p9, p3); l12 = gmsh.model.geo.addLine(p3, p11)
loop_iron = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6, l7, l8, l9, l10]) iron_surf = gmsh.model.geo.addPlaneSurface([loop_iron]) loop_copper = gmsh.model.geo.addCurveLoop([l11, l12, l7, l8]) copper_surf = gmsh.model.geo.addPlaneSurface([loop_copper])
p20 = gmsh.model.geo.addPoint(-100, -100, 0); p21 = gmsh.model.geo.addPoint(100, -100, 0) p22 = gmsh.model.geo.addPoint(100, 100, 0); p23 = gmsh.model.geo.addPoint(-100, 100, 0) la1 = gmsh.model.geo.addLine(p20, p21); la2 = gmsh.model.geo.addLine(p21, p22) la3 = gmsh.model.geo.addLine(p22, p23); la4 = gmsh.model.geo.addLine(p23, p20) air_loop = gmsh.model.geo.addCurveLoop([la1, la2, la3, la4]) air_surf = gmsh.model.geo.addPlaneSurface([air_loop, loop_iron])
gmsh.model.geo.synchronize() gmsh.model.addPhysicalGroup(2, [iron_surf], 1); gmsh.model.setPhysicalName(2, 1, "IronCore") gmsh.model.addPhysicalGroup(2, [copper_surf], 2); gmsh.model.setPhysicalName(2, 2, "CopperRing") gmsh.model.addPhysicalGroup(2, [air_surf], 3); gmsh.model.setPhysicalName(2, 3, "Air") gmsh.option.setNumber("Mesh.MeshSizeMax", 4.0) gmsh.model.mesh.generate(2) gmsh.write("magnet.msh") gmsh.finalize()
# === 2. CONVERSION === def convert_mesh(): msh = meshio.read("magnet.msh") triangle_data = msh.get_cell_data("gmsh:physical", "triangle").astype(np.uint64) out_mesh = meshio.Mesh(points=msh.points[:, :2], cells=[("triangle", msh.get_cells_type("triangle"))], cell_data={"subdomains": [triangle_data]}) out_mesh.write("magnet.xml")
# === 3. SOLVER === def solve_ac_magnet(): # Clean output folder to avoid XML parsing errors if os.path.exists("output"): shutil.rmtree("output") os.makedirs("output")
mesh = Mesh("magnet.xml") subdomains = MeshFunction("size_t", mesh, "magnet_subdomains.xml") dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
freq, I_peak = 120, 5e7 dt = 1.0 / (freq * 50) mu0, mu_iron, sigma_copper = 4*np.pi*1e-7, 2000*4*np.pi*1e-7, 5.8e7
V = FunctionSpace(mesh, "P", 1) W = VectorFunctionSpace(mesh, "DG", 0) # For B-Field vectors A_prev, A_res = Function(V), Function(V, name="MagneticPotential") B_res = Function(W, name="B-Field") A_prev.assign(Constant(0.0))
class MatProps(UserExpression): def __init__(self, val_map, **kwargs): super().__init__(**kwargs) self.val_map = val_map def eval_cell(self, values, x, cell): tag = subdomains[cell.index] values[0] = self.val_map.get(tag, 10.0) def value_shape(self): return ()
nu = MatProps({1: 1.0/mu_iron, 2: 1.0/mu0, 3: 1.0/mu0}, degree=0) sigma = MatProps({1: 10.0, 2: sigma_copper, 3: 10.0}, degree=0)
A, v = TrialFunction(V), TestFunction(V) bc = DirichletBC(V, Constant(0.0), "on_boundary") t = 0.0 J_src = Expression("I * sin(2*pi*f*t)", I=I_peak, f=freq, t=t, degree=0) a = (sigma * A * v / dt + nu * dot(grad(A), grad(v))) * dx L = (sigma * A_prev * v / dt + J_src * v) * dx(1)
p_u, p_s = Point(-15, 0), Point(15, 18) history_u, history_s, time_log = [], [], []
pvd_a = File("output/flux.pvd") pvd_b = File("output/b_field.pvd")
print("Solving Loop Started...") for n in range(200): t += dt J_src.t = t # Using default petsc solver to bypass mumps-related hardware location errors solve(a == L, A_res, bc) # Calculate B-field (B = curl A) B_res.assign(project(as_vector((A_res.dx(1), -A_res.dx(0))), W)) pvd_a << (A_res, t) pvd_b << (B_res, t) if n > 150: history_u.append(A_res(p_u)) history_s.append(A_res(p_s)) time_log.append(t) A_prev.assign(A_res)
if n % 50 == 0: print(f"Step {n} solved...")
# --- STEADY STATE OUTPUT --- y1, y2 = np.array(history_u), np.array(history_s) y1_cent, y2_cent = y1 - np.mean(y1), y2 - np.mean(y2) norm1, norm2 = np.linalg.norm(y1_cent), np.linalg.norm(y2_cent) if norm1 > 1e-15 and norm2 > 1e-15: phi = np.degrees(np.arccos(np.clip(np.dot(y1_cent, y2_cent)/(norm1*norm2), -1, 1))) print(f"\nPHASE ANGLE: {phi:.4f} degrees") y1_norm, y2_norm = y1_cent / np.max(np.abs(y1_cent)), y2_cent / np.max(np.abs(y2_cent)) plt.figure(figsize=(10, 5)) plt.plot(time_log, y1_norm, label='Unshaded (Main)') plt.plot(time_log, y2_norm, 'r--', label='Shaded (Copper)') plt.axhline(0, color='black', linewidth=1, linestyle=':') plt.title(f"Steady State Phase Shift: {phi:.2f} Deg") plt.xlabel("Time (s)"); plt.ylabel("Normalized Potential"); plt.legend(); plt.grid(True) plt.savefig("phase_plot.png") print("Graph saved as phase_plot.png")
if __name__ == "__main__": generate_custom_mesh(); convert_mesh(); solve_ac_magnet()
Que/Gravock
|