diff options
| -rw-r--r-- | coil_mag_materials.yml | 43 | ||||
| -rw-r--r-- | coil_mag_solvers.yml | 6 | ||||
| -rw-r--r-- | coil_parasitics.py | 85 | ||||
| -rw-r--r-- | twisted_coil_gen_twolayer.py | 18 |
4 files changed, 96 insertions, 56 deletions
diff --git a/coil_mag_materials.yml b/coil_mag_materials.yml index 9efbb80..2f19103 100644 --- a/coil_mag_materials.yml +++ b/coil_mag_materials.yml @@ -9,31 +9,20 @@ ro4003c: Density: 1790 # 23°C Relative Permeability: 1 Relative Permittivity: 3.55 +fr4: + Density: 1850 # 23°C + Relative Permeability: 1 + Relative Permittivity: 4.4 + Heat Conductivity: 0.81 # in-plane ideal: Relative Permittivity: 1 copper: Density: 8960.0 # 0°C - Electric Conductivity: 32300000 # 200°C + Electric Conductivity: 59600000 # 20°C Emissivity: 0.012 # 327°C Heat Capacity: 415.0 # 200°C Heat Conductivity: 401.0 # 0°C Relative Permeability: 1 - Relative Permittivity: 1 -graphite_CZ3-R6300: # crucible - Density: 1730.0 - Electric Conductivity: 58800 - Emissivity: 0.81 # 205°C - Heat Capacity: 1237.0 - Heat Conductivity: 65 # 20°C - Relative Permeability: 1 - Relative Permittivity: 1 -graphite_FU8957: # heater - Density: 1750.0 - Emissivity: 0.81 # 250°C - Heat Capacity: 1237.0 - Heat Conductivity: 105 # averaged over different given values - Relative Permeability: 1 - Relative Permittivity: 1 steel_1.4541: Density: 7900.0 # 20°C Electric Conductivity: 1370 @@ -42,26 +31,6 @@ steel_1.4541: Heat Conductivity: 15.0 # 20°C Relative Permeability: 1 Relative Permittivity: 1 -tin_liquid: - Density: 6980.0 - Electric Conductivity: 2080000 - Emissivity: 0.064 # set equal to solid - Heat Capacity: 252.7 - Heat Conductivity: 29.0 - Relative Permeability: 1 - Relative Permittivity: 1 - Liquid: 'Logical True' -tin_solid: - Density: 7179.0 - Electric Conductivity: 4380000 - Emissivity: 0.064 - Heat Capacity: 244.0 - Heat Conductivity: 60.0 - Relative Permeability: 1 - Relative Permittivity: 1 - Solid: 'Logical True' - Melting Point: 505 - Latent Heat: 59600 water: Density: 1000.0 Heat Capacity: 4182.0 diff --git a/coil_mag_solvers.yml b/coil_mag_solvers.yml index bb7b169..496509c 100644 --- a/coil_mag_solvers.yml +++ b/coil_mag_solvers.yml @@ -1,16 +1,17 @@ Static_Current_Conduction: Equation: Static Current Conduction Variable: Potential + Variable DOFs: 1 Procedure: '"StatCurrentSolve" "StatCurrentSolver"' Calculate Volume Current: True Calculate Joule Heating: False Optimize Bandwidth: True Nonlinear System Max Iterations: 1 Linear System Solver: Iterative - Linear System Iterative Method: BiCGStabl + Linear System Iterative Method: CG Linear System Max Iterations: 5000 Linear System Convergence Tolerance: 1.0e-10 - Linear System Preconditioning: ILU0 + Linear System Preconditioning: ILU3 Linear System ILUT Tolerance: 1.0e-3 Linear System Abort Not Converged: False Linear System Residual Output: 20 @@ -21,6 +22,7 @@ Magneto_Dynamics: Procedure: '"MagnetoDynamics" "WhitneyAVSolver"' Newton-Raphson Iteration: True Nonlinear System Max Iterations: 1 + Fix Input Current Density: True Linear System Solver: Iterative Linear System Preconditioning: none Linear System Convergence Tolerance: 1e-8 diff --git a/coil_parasitics.py b/coil_parasitics.py index 47d5cd2..08a3603 100644 --- a/coil_parasitics.py +++ b/coil_parasitics.py @@ -1,15 +1,16 @@ #!/usr/bin/env python3 +import math from pathlib import Path import multiprocessing import re import tempfile -import subprocess import fnmatch import shutil import numpy as np from pyelmer import elmer +import subprocess_tee import click from scipy import constants @@ -31,6 +32,17 @@ def enumerate_mesh_bodies(msh_file): tag, _, name = line.partition(' ') yield name.strip().strip('"'), (int(dim), int(tag)) +# https://en.wikipedia.org/wiki/Metric_prefix +SI_PREFIX = 'QRYZEPTGMk mµnpfazyrq' + +def format_si(value, unit='', fractional_digits=1): + mag = int(math.log10(abs(value))//3) + value /= 1000**mag + prefix = SI_PREFIX[SI_PREFIX.find(' ') - mag].strip() + value = f'{{:.{fractional_digits}f}}'.format(value) + return f'{value} {prefix}{unit}' + + INPUT_EXT_MAP = { '.grd': 1, '.mesh*': 2, @@ -56,7 +68,7 @@ OUTPUT_EXT_MAP = { '.msh': 4, '.vtu': 5} -def elmer_grid(infile, outfile=None, intype=None, outtype=None, cwd=None, **kwargs): +def elmer_grid(infile, outfile=None, intype=None, outtype=None, cwd=None, stdout_log=None, stderr_log=None, **kwargs): infile = Path(infile) if outfile is not None: outfile = Path(outfile) @@ -80,10 +92,20 @@ def elmer_grid(infile, outfile=None, intype=None, outtype=None, cwd=None, **kwar args.extend(str(v) for v in value) else: args.append(str(value)) - subprocess.run(args, cwd=cwd, check=True) -def elmer_solver(cwd): - subprocess.run(['ElmerSolver'], cwd=cwd) + result = subprocess_tee.run(args, cwd=cwd, check=True) + if stdout_log: + Path(stdout_log).write_text(result.stdout or '') + if stderr_log: + Path(stderr_log).write_text(result.stderr or '') + +def elmer_solver(cwd, stdout_log=None, stderr_log=None): + result = subprocess_tee.run(['ElmerSolver'], cwd=cwd, check=True) + if stdout_log: + Path(stdout_log).write_text(result.stdout or '') + if stderr_log: + Path(stderr_log).write_text(result.stderr or '') + return result @click.command() @@ -132,14 +154,17 @@ def run_capacitance_simulation(mesh_file, sim_dir): boundary_airbox.data['Electric Infinity BC'] = 'True' with tempfile.TemporaryDirectory() as tmpdir: - if sim_dir: - tmpdir = str(sim_dir) + tmpdir = sim_dir if sim_dir else Path(tmpdir) sim.write_startinfo(tmpdir) sim.write_sif(tmpdir) # Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units) - elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3]) - elmer_solver(tmpdir) + elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3], + stdout_log=(tmpdir / 'ElmerGrid_stdout.log'), + stderr_log=(tmpdir / 'ElmerGrid_stderr.log')) + elmer_solver(tmpdir, + stdout_log=(tmpdir / 'ElmerSolver_stdout.log'), + stderr_log=(tmpdir / 'ElmerSolver_stderr.log')) capacitance_matrix = np.loadtxt(tmpdir / 'capacitance.txt') @@ -194,6 +219,9 @@ def run_inductance_simulation(mesh_file, sim_dir): bdy_if_bottom.material = copper bdy_if_bottom.equation = copper_eqn + potential_force = elmer.BodyForce(sim, 'electric_potential', {'Electric Potential': 'Equals "Potential"'}) + bdy_trace.body_force = potential_force + # boundaries boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]]) boundary_airbox.data['Electric Infinity BC'] = 'True' @@ -206,15 +234,46 @@ def run_inductance_simulation(mesh_file, sim_dir): boundary_vminus.data['Potential'] = 0.0 with tempfile.TemporaryDirectory() as tmpdir: - if sim_dir: - tmpdir = str(sim_dir) + tmpdir = sim_dir if sim_dir else Path(tmpdir) sim.write_startinfo(tmpdir) sim.write_sif(tmpdir) # Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units) - elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3]) - elmer_solver(tmpdir) + elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3], + stdout_log=(tmpdir / 'ElmerGrid_stdout.log'), + stderr_log=(tmpdir / 'ElmerGrid_stderr.log')) + solver_stdout, solver_stderr = (tmpdir / 'ElmerSolver_stdout.log'), (tmpdir / 'ElmerSolver_stderr.log') + res = elmer_solver(tmpdir, + stdout_log=solver_stdout, + stderr_log=solver_stderr) + + P, R, U_mag = None, None, None + solver_error = False + for l in res.stdout.splitlines(): + if (m := re.fullmatch(r'StatCurrentSolve:\s*Total Heating Power\s*:\s*([0-9.+-Ee]+)\s*', l)): + P = float(m.group(1)) + elif (m := re.fullmatch(r'StatCurrentSolve:\s*Effective Resistance\s*:\s*([0-9.+-Ee]+)\s*', l)): + R = float(m.group(1)) + elif (m := re.fullmatch(r'MagnetoDynamicsCalcFields:\s*ElectroMagnetic Field Energy\s*:\s*([0-9.+-Ee]+)\s*', l)): + U_mag = float(m.group(1)) + elif re.fullmatch(r'IterSolve: Linear iteration did not converge to tolerance', l): + solver_error = True + + if solver_error: + raise click.ClickException(f'Error: One of the solvers did not converge. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}') + elif P is None or R is None or U_mag is None: + raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}') + + U = math.sqrt(P*R) + I = math.sqrt(P/R) + L = 2*U_mag / (I**2) + + assert math.isclose(U, 1.0, abs_tol=1e-3) + + print(f'Coil resistance calculated by solver: {format_si(R, "Ω")}') + print(f'Inductance calucated from field: {format_si(L, "H")}') if __name__ == '__main__': run_inductance_simulation() + diff --git a/twisted_coil_gen_twolayer.py b/twisted_coil_gen_twolayer.py index 101042c..5f7e685 100644 --- a/twisted_coil_gen_twolayer.py +++ b/twisted_coil_gen_twolayer.py @@ -317,6 +317,8 @@ def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log airbox_physical = gmsh.model.add_physical_group(3, [airbox], name='airbox') trace_physical = gmsh.model.add_physical_group(3, [toplevel_tag], name='trace') + gmsh.model.mesh.setSize([(0, tag) for dim, tag in gmsh.model.getBoundary([(3, toplevel_tag)], recursive=True) if dim == 0], 0.300) + #interface_tags_top = gmsh.model.getBoundary([(3, contact_tag_top)], oriented=False) #interface_tags_bottom = gmsh.model.getBoundary([(3, contact_tag_bottom)], oriented=False) @@ -464,6 +466,8 @@ def print_valid_twists(ctx, param, value): @click.option('--via-offset', type=float, default=None, help='Radially offset vias from trace endpoints [mm]') @click.option('--keepout-zone/--no-keepout-zone', default=True, help='Add a keepout are to the footprint (default: yes)') @click.option('--keepout-margin', type=float, default=5, help='Margin between outside of coil and keepout area (mm, default: 5)') +@click.option('--copper-thickness', type=float, default=0.035, help='Copper thickness for resistance calculation and mesh generation in mm. Default: 0.035mm ^= 1 Oz') +@click.option('--board-thickness', type=float, default=1.53, help='Board substrate thickness for mesh generation in mm. Default: 1.53mm') @click.option('--twists', type=int, default=1, help='Number of twists per revolution. Note that this number must be co-prime to the number of turns. Run with --show-twists to list valid values. (default: 1)') @click.option('--circle-segments', type=int, default=64, help='When not using arcs, the number of points to use for arc interpolation per 360 degrees.') @click.option('--show-twists', callback=print_valid_twists, expose_value=False, type=int, is_eager=True, help='Calculate and show valid --twists counts for the given number of turns. Takes the number of turns as a value.') @@ -477,7 +481,8 @@ def print_valid_twists(ctx, param, value): @click.version_option() def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_drill, via_offset, trace_width, clearance, footprint_name, layer_pair, twists, clipboard, counter_clockwise, keepout_zone, keepout_margin, - arc_tolerance, pcb, mesh_out, magneticalc_out, circle_segments, mag_mesh_out): + arc_tolerance, pcb, mesh_out, magneticalc_out, circle_segments, mag_mesh_out, copper_thickness, + board_thickness): if 'WAYLAND_DISPLAY' in os.environ: copy, paste, cliputil = ['wl-copy'], ['wl-paste'], 'xclip' else: @@ -722,7 +727,12 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d svg_vias.append(Tag('circle', cx=xv, cy=yv, r=via_diameter/2, stroke='none', fill='white')) svg_vias.append(Tag('circle', cx=xv, cy=yv, r=via_drill/2, stroke='none', fill='black')) - print(f'Approximate track length: {clen*twists*2:.2f} mm', file=sys.stderr) + l_total = clen*twists*2 + print(f'Approximate track length: {l_total:.2f} mm', file=sys.stderr) + A = copper_thickness/1e3 * trace_width/1e3 + rho = 1.68e-8 + R = l_total/1e3 * rho / A + print(f'Approximate resistance: {R:g} Ω', file=sys.stderr) top_pad = make_pad(1, [layer_pair[0]], outer_radius, 0) pads.append(top_pad) @@ -812,10 +822,10 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d r = outer_diameter/2 + 20 if mesh_out: - traces_to_gmsh(traces, mesh_out, ((-r, -r), (r, r))) + traces_to_gmsh(traces, mesh_out, ((-r, -r), (r, r)), copper_thickness=copper_thickness, board_thickness=board_thickness) if mag_mesh_out: - traces_to_gmsh_mag(traces, mag_mesh_out, ((-r, -r), (r, r))) + traces_to_gmsh_mag(traces, mag_mesh_out, ((-r, -r), (r, r)), copper_thickness=copper_thickness, board_thickness=board_thickness) if magneticalc_out: traces_to_magneticalc(traces, magneticalc_out) |