greg/pcb-stator-coil-generator
1
0
Fork 0
mirror of https://github.com/atomic14/kicad-coil-plugins.git synced 2024-10-18 09:06:57 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
pcb-stator-coil-generator/simulations/magnet_field_single_coil.ipynb
Chris Greening 3db8c78be0 Trying out some new visualisations
2022-11-19 08:24:27 +00:00

525 lines
17 KiB
Text
Raw Blame History

{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from biot_savart_v4_3 import parse_coil, plot_coil, slice_coil, plot_coil2\n",
"from tqdm.notebook import trange, tqdm\n",
"from helpers import get_arc_point, draw_arc, rotate, scale, rotate_point, translate"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Track width and spacing\n",
"TRACK_WIDTH = 0.127\n",
"TRACK_SPACING = 0.127\n",
"\n",
"# via defaults\n",
"VIA_DIAM = 0.8\n",
"VIA_DRILL = 0.4\n",
"\n",
"# this is for a 1.27mm pitch pin\n",
"PIN_DIAM = 1.0\n",
"PIN_DRILL = 0.65\n",
"\n",
"# this is for the PCB connector - see https://www.farnell.com/datasheets/2003059.pdf\n",
"PAD_WIDTH = 3\n",
"PAD_HEIGHT = 2\n",
"PAD_PITCH = 2.5\n",
"\n",
"STATOR_HOLE_RADIUS = 5.5\n",
"HOLE_SPACING = 0.25\n",
"\n",
"# PCB Edge size\n",
"STATOR_RADIUS = 30\n",
"\n",
"# where to puth the mounting pins\n",
"SCREW_HOLE_DRILL_DIAM = 2.3 # 2.3mm drill for a 2mm screw\n",
"SCREW_HOLE_RADIUS = STATOR_RADIUS\n",
"\n",
"# Coil params\n",
"TURNS = 16\n",
"COIL_CENTER_RADIUS = 19.95\n",
"COIL_VIA_RADIUS = 20.95\n",
"\n",
"\n",
"def get_points(spacing, inner_radius, outer_radius, start_angle, end_angle):\n",
" # first calculate the angle step size from the spacing and the inner_radius\n",
" spacing_angle = np.rad2deg(np.arctan2(spacing, inner_radius))\n",
" print(spacing_angle)\n",
" # now calculate the points be iterating from start_angle to end_angle with the spacing_angle\n",
" points = []\n",
" for angle in np.arange(start_angle, end_angle, spacing_angle * 2):\n",
" p1 = get_arc_point(angle, inner_radius)\n",
" points.append([p1[0], p1[1], 0, 0.5])\n",
" p2 = get_arc_point(angle + spacing_angle, outer_radius)\n",
" points.append([p2[0], p2[1], 0, 0.5])\n",
" points.append([p2[0], p2[1], -0.8, 0.5])\n",
" points.append([p1[0], p1[1], -0.8, 0.5])\n",
" return points\n",
"\n",
"\n",
"coil1 = get_points(\n",
" TRACK_SPACING + TRACK_WIDTH,\n",
" (COIL_CENTER_RADIUS - 10),\n",
" COIL_CENTER_RADIUS + 10,\n",
" -15,\n",
" 15,\n",
")\n",
"# move all the points in by COIL_CENTER_RADIUS\n",
"for p in coil1:\n",
" p[0] -= COIL_CENTER_RADIUS\n",
"# rotate the points by 90 degrees\n",
"for p in coil1:\n",
" p[0], p[1] = rotate_point(p[0], p[1], 90)\n",
"coil1 = np.array(coil1).T\n",
"plot_coil2(coil1)\n",
"coil1 = slice_coil(coil1, 0.1)\n",
"coil1 = coil1.T\n",
"print(coil1.shape)\n",
"print(len(coil1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple Simulation of a dipole magnet"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# magnetic field at a point x,y,z of a dipole magnet with moment m in the z direction\n",
"def B(x, y, z, m=0.185):\n",
" mu0 = 4 * np.pi * 1e-7\n",
" r = np.sqrt(x**2 + y**2 + z**2)\n",
" return (\n",
" np.array(\n",
" [3 * x * z / r**5, 3 * y * z / r**5, (3 * z**2 / r**5 - 1 / r**3)]\n",
" )\n",
" * m\n",
" * mu0\n",
" )\n",
"\n",
"\n",
"def plot_field_slice(x, y, bx, by, mag, name=\"magnetic_field.png\"):\n",
" # plot the magnetic field\n",
" fig = plt.figure()\n",
" ax = fig.add_subplot(111)\n",
" ax.streamplot(\n",
" x,\n",
" y,\n",
" bx,\n",
" by,\n",
" linewidth=1,\n",
" cmap=plt.cm.inferno,\n",
" density=2,\n",
" arrowstyle=\"->\",\n",
" arrowsize=1.5,\n",
" )\n",
"\n",
" ax.set_xlabel(\"$x$\")\n",
" ax.set_ylabel(\"$y$\")\n",
" ax.set_xlim(-0.1, 0.1)\n",
" ax.set_ylim(-0.1, 0.1)\n",
" ax.set_aspect(\"equal\")\n",
"\n",
" # plot the magniture of the field as an image\n",
" im = ax.imshow(\n",
" mag, extent=[-0.1, 0.1, -0.1, 0.1], origin=\"lower\", cmap=plt.cm.inferno\n",
" )\n",
"\n",
" fig.show()\n",
" # save the figure\n",
" fig.savefig(name)\n",
"\n",
"\n",
"# # calculate the magnetic field at y = 0, over z = -1, 1 and x = -1, 1\n",
"x = np.linspace(-0.1, 0.1, 100)\n",
"z = np.linspace(-0.1, 0.1, 100)\n",
"X, Z = np.meshgrid(x, z)\n",
"Bx, By, Bz = B(X, 0, Z)\n",
"\n",
"plot_field_slice(\n",
" X,\n",
" Z,\n",
" Bx,\n",
" Bz,\n",
" np.log(np.sqrt(Bx**2 + By**2 + Bz**2)),\n",
" \"magnetic_field_side.png\",\n",
")\n",
"\n",
"# # calculate the magnetic field at z = 1, over y = -1, 1 and x = -1, 1\n",
"x = np.linspace(-0.1, 0.1, 100)\n",
"y = np.linspace(-0.1, 0.1, 100)\n",
"X, Y = np.meshgrid(x, y)\n",
"Bx, By, Bz = B(X, Y, 0.01)\n",
"\n",
"plot_field_slice(\n",
" X, Y, Bx, By, np.sqrt(Bx**2 + By**2 + Bz**2), \"magnetic_field_bottom.png\"\n",
")\n",
"\n",
"# calculate the magnetic field in a 3d volume\n",
"x = np.linspace(-0.1, 0.1, 100)\n",
"y = np.linspace(-0.1, 0.1, 100)\n",
"z = np.linspace(-0.1, 0.1, 100)\n",
"X, Y, Z = np.meshgrid(x, y, z)\n",
"Bx, By, Bz = B(X, Y, Z)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# calculate the force on a wire of length l carrying current I at a point x,y,z with a direction vector d\n",
"def F(p, d, I, l):\n",
" return I * l * np.cross(d, B(p[0], p[1], p[2]))\n",
"\n",
"\n",
"def calculate_forces_on_wire_points(points):\n",
" Fx = []\n",
" Fy = []\n",
" Fz = []\n",
" # calculate the force on each point\n",
" for i in range(len(points)):\n",
" # calculate the direction vector\n",
" dx = points[i][0] - points[(i + 1) % len(points)][0]\n",
" dy = points[i][1] - points[(i + 1) % len(points)][1]\n",
" dz = points[i][2] - points[(i + 1) % len(points)][2]\n",
" d = np.array([dx, dy, dz])\n",
" # get the length of d\n",
" l = np.sqrt(dx**2 + dy**2 + dz**2)\n",
" if l > 0:\n",
" # normalise d\n",
" d = d / l\n",
" # calculate the force\n",
" fx, fy, fz = F(points[i], d, points[i][3], l)\n",
" Fx.append(fx)\n",
" Fy.append(fy)\n",
" Fz.append(fz)\n",
" else:\n",
" Fx.append(0)\n",
" Fy.append(0)\n",
" Fz.append(0)\n",
" return Fx, Fy, Fz\n",
"\n",
"\n",
"# locate the loop of wire directly below the magnet\n",
"x = 0.01\n",
"y = 0\n",
"z = -0.002\n",
"r1 = 0.001\n",
"r2 = 0.01\n",
"\n",
"\n",
"points = coil1.copy()\n",
"# scale the points from mm to m\n",
"# shift the coil to the correct position\n",
"for i in range(len(points)):\n",
" points[i][0] = points[i][0] / 1000 + x\n",
" points[i][1] = points[i][1] / 1000 + y\n",
" points[i][2] = points[i][2] / 1000 + z\n",
"\n",
"\n",
"# plot the points in 2D x,y\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111)\n",
"ax.plot([p[0] for p in points], [p[1] for p in points])\n",
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$y$\")\n",
"ax.set_xlim(-0.01 + x, 0.01 + x)\n",
"ax.set_ylim(-0.01 + y, 0.01 + y)\n",
"ax.set_aspect(\"equal\")\n",
"fig.show()\n",
"\n",
"Fx, Fy, Fz = calculate_forces_on_wire_points(points)\n",
"\n",
"\n",
"# plot the wire along with arrows showing the force\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111, projection=\"3d\")\n",
"ax.plot([p[0] for p in points], [p[1] for p in points], [p[2] for p in points])\n",
"\n",
"print(sum(Fx) / 9.8)\n",
"\n",
"# subsample the points and force vectors to make the plot clearer\n",
"points = points[::50]\n",
"Fx = Fx[::50]\n",
"Fy = Fy[::50]\n",
"Fz = Fz[::50]\n",
"\n",
"ax.quiver(\n",
" [p[0] for p in points],\n",
" [p[1] for p in points],\n",
" [p[2] for p in points],\n",
" np.sqrt(Fx),\n",
" 0,\n",
" 0,\n",
" length=0.01,\n",
" normalize=True,\n",
")\n",
"ax.set_xlabel(\"$x$\")\n",
"ax.set_ylabel(\"$y$\")\n",
"ax.set_zlabel(\"$z$\")\n",
"ax.set_xlim(x - 0.02, x + 0.02)\n",
"ax.set_ylim(y - 0.02, y + 0.02)\n",
"ax.set_zlim(z - 0.02, z + 0.02)\n",
"ax.set_aspect(\"equal\")\n",
"\n",
"# change the figure size\n",
"fig.set_size_inches(10, 10)\n",
"\n",
"fig.show()\n",
"# coil2 = -0.01311082557350831\n",
"# coil1 = 0.0088435517657609"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def sweep_coil(coil, X):\n",
" Y = np.zeros(len(X))\n",
" Z = -0.01 * np.ones(len(X))\n",
"\n",
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
" Fx = []\n",
" Fy = []\n",
" Fz = []\n",
" for p in trange(len(X)):\n",
" points = coil.copy()\n",
" # scale the points from mm to m\n",
" # shift the coil to the correct position\n",
" for i in range(len(points)):\n",
" points[i][0] = points[i][0] / 1000 + X[p]\n",
" points[i][1] = points[i][1] / 1000 + Y[p]\n",
" points[i][2] = points[i][2] / 1000 + Z[p]\n",
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
" Fx.append(sum(Fx_))\n",
" Fy.append(sum(Fy_))\n",
" Fz.append(sum(Fz_))\n",
" return Fx, Fy, Fz\n",
"\n",
"\n",
"# sweep the coild from -3cm to 3cm in 0.01m steps\n",
"X = np.linspace(-0.03, 0.03, 100)\n",
"Fx_1_straight, Fy_1_straight, Fz_1_straight = sweep_coil(coil1, X)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the force as a function of x\n",
"plt.plot(X, -np.array(Fx_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
"# plot a dotted line along y = 0\n",
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the force as a function of x\n",
"plt.plot(X, -np.array(Fz_1_straight) / 9.8, label=\"coil1\", color=\"red\")\n",
"# plot a dotted line along y = 0\n",
"plt.plot([X[0], X[-1]], [0, 0], \"--\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
"def sweep_coil_cirlc(coil, coil_center_radius, theta):\n",
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
" Z = -0.01 * np.ones(100)\n",
"\n",
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
" Torque = []\n",
" Fx = []\n",
" Fy = []\n",
" Fz = []\n",
" for p in trange(len(theta)):\n",
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
" x = X[p]\n",
" y = Y[p]\n",
" z = Z[p]\n",
"\n",
" points = coil.copy()\n",
" for i in range(len(points)):\n",
" px = points[i][0] / 1000\n",
" py = points[i][1] / 1000\n",
" pz = points[i][2] / 1000\n",
" # rotate the points so the coil is correctly oriented\n",
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
" points[i][1] = (\n",
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
" )\n",
" points[i][2] = pz + z\n",
" # feel the force\n",
" Fx_, Fy_, Fz_ = calculate_forces_on_wire_points(points)\n",
" Fx.append(sum(Fx_))\n",
" Fy.append(sum(Fy_))\n",
" Fz.append(sum(Fz_))\n",
" # calculate the torque - which should be 90 degress to the angle\n",
" torque_angle = np.deg2rad(theta[p] - 90)\n",
" Torque.append(sum(Fx_) * np.cos(torque_angle) + sum(Fy_) * np.sin(torque_angle))\n",
"\n",
" return Fx, Fy, Fz, Torque\n",
"\n",
"\n",
"# sweep the coils from -45 to 45 degrees in 1 degree steps\n",
"theta = np.linspace(0, 180, 100)\n",
"Fx_1_curve, Fy_1_curve, Fz_1_curve, Torque_1 = sweep_coil_cirlc(\n",
" coil1, 19.5 / 1000, theta\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(theta, -np.array(Torque_1) / 9.8, label=\"coil1\", color=\"red\")\n",
"# plot a dotted line along y = 0\n",
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot arrows for the Fx and Fy components\n",
"X = 19.5 * np.cos(np.deg2rad(theta)) / 1000\n",
"Y = 19.5 * np.sin(np.deg2rad(theta)) / 1000\n",
"plt.quiver(X[::5], Y[::5], Fx_1_curve[::5], Fy_1_curve[::5], color=\"red\")\n",
"# make the axis equal so the arrows are not stretched\n",
"plt.axis(\"equal\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(theta, -np.array(Fz_1_curve) / 9.8, label=\"coil1\", color=\"red\")\n",
"plt.plot(theta, np.array(Fz_2_curve) / 9.8, label=\"coil2\", color=\"blue\")\n",
"# plot a dotted line along y = 0\n",
"plt.plot([theta[0], theta[-1]], [0, 0], \"--\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# instead of sweeping horizontally, we'll sweep the coils around a circle\n",
"def sweep_coil_cirlc(coil, coil_center_radius):\n",
" # sweep the coils from -45 to 45 degrees in 1 degree steps\n",
" theta = np.linspace(0, 180, 20)\n",
" X = coil_center_radius * np.cos(np.deg2rad(theta))\n",
" Y = coil_center_radius * np.sin(np.deg2rad(theta))\n",
" Z = 1 * np.ones(100)\n",
"\n",
" # loop through the locations and calculate the forces, sum up the force in the X direction for each location\n",
" Fx = []\n",
" Fy = []\n",
" Fz = []\n",
" for p in trange(len(theta)):\n",
" angle = np.deg2rad(theta[p]) - np.pi / 2\n",
" x = X[p]\n",
" y = Y[p]\n",
" z = Z[p]\n",
"\n",
" points = coil.copy()\n",
" for i in range(len(points)):\n",
" px = points[i][0] / 1000\n",
" py = points[i][1] / 1000\n",
" pz = points[i][2] / 1000\n",
" # rotate the points so the coil is correctly oriented\n",
" points[i][0] = px * np.cos(angle) - py * np.sin(angle) + x\n",
" points[i][1] = (\n",
" px * np.sin(angle) + py * np.cos(angle) + y - coil_center_radius\n",
" )\n",
" points[i][2] = pz + z\n",
" plt.plot([p[0] for p in points], [p[1] for p in points], linewidth=0.5)\n",
" # add the torque arrow to the plot\n",
" torque_angle = np.deg2rad(theta[p] - 90)\n",
" torque_x1 = x\n",
" torque_y1 = y\n",
" torque_x2 = x + 0.01 * np.cos(torque_angle)\n",
" torque_y2 = y + 0.01 * np.sin(torque_angle)\n",
" plt.arrow(\n",
" torque_x1,\n",
" torque_y1,\n",
" torque_x2 - torque_x1,\n",
" torque_y2 - torque_y1,\n",
" head_width=0.001,\n",
" head_length=0.002,\n",
" fc=\"k\",\n",
" ec=\"k\",\n",
" )\n",
"\n",
"\n",
"sweep_coil_cirlc(coil2, 19.5 / 1000)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.7 ('venv': venv)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.7"
},
"vscode": {
"interpreter": {
"hash": "1ce20143987840b9786ebb5907032c9c3a8efacbb887dbb0ebc4934f2ad26cb3"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Reference in a new issue View git blame Copy permalink