fix(datamodel): fix RectStrands external radius

Author: amaurigmartins

examples/tutorial3.jl

@@ -354,211 +354,21 @@ F = FormulationSet(:FEM,
 	options = opts,
 );
 
-# # Run the FEM solver
-# @time ws, p = compute!(problem, F);
+# Run the FEM solver
+@time ws, p = compute!(problem, F);
 
-# # Display computation results
-# per_km(p, 1; mode = :RLCG, tol = 1e-9)
+# Display computation results
+per_km(p, 1; mode = :RLCG, tol = 1e-9)
 
-# # Export ZY matrices to ATPDraw
-# output_file = fullfile("ZY_export.xml")
-# export_file = export_data(:atp, p; file_name = output_file, cable_system = cable_system);
+# Export ZY matrices to ATPDraw
+output_file = fullfile("ZY_export.xml")
+export_file = export_data(:atp, p; file_name = output_file, cable_system = cable_system);
 
-# # Obtain the symmetrical components via Fortescue transformation
-# Tv, p012 = Fortescue(tol = 1e-5)(p);
+# Obtain the symmetrical components via Fortescue transformation
+Tv, p012 = Fortescue(tol = 1e-5)(p);
 
-# # Inspect the transformed matrices
-# per_km(p012, 1; mode = :ZY, tol = 1e-9)
+# Inspect the transformed matrices
+per_km(p012, 1; mode = :ZY, tol = 1e-9)
 
-# # Or the corresponding lumped circuit quantities
-# per_km(p012, 1; mode = :RLCG, tol = 1e-9)
-
-using LineCableModels.Utils: to_nominal
-using LineCableModels.DataModel: calc_circstrands_coords
-
-function collect_layer_geometry(design::CableDesign;
-	x_offset::Real = 0.0,
-	y_offset::Real = 0.0,
-)
-	layer_no = Ref(0)
-	out = Vector{
-		NamedTuple{(:layer, :type, :element, :x, :y, :radius_out, :rin, :rext),
-			Tuple{Int, String, Int, Float64, Float64, Float64, Float64, Float64}},
-	}()
-
-	function push_layer!(layer, x0, y0)
-		if layer isa CircStrands
-			layer_no[] += 1
-			rwire = to_nominal(layer.radius_wire)
-			lay_r = layer.num_wires == 1 ? 0.0 : to_nominal(layer.radius_in)
-			coords = calc_circstrands_coords(layer.num_wires, rwire, lay_r; C = (x0, y0))
-			for (i, (x, y)) in enumerate(coords)
-				push!(out, (layer_no[], "circstrands", i, x, y, rwire, NaN, NaN))
-			end
-		elseif layer isa Strip || layer isa Tubular || layer isa Semicon ||
-			   layer isa Insulator
-			layer_no[] += 1
-			rin = to_nominal(layer.radius_in);
-			rex = to_nominal(layer.radius_ext)
-			typ = lowercase(string(nameof(typeof(layer))))
-			push!(out, (layer_no[], typ, 1, x0, y0, rex, rin, rex))
-		elseif layer isa ConductorGroup
-			for sub in layer.layers
-				push_layer!(sub, x0, y0)
-			end
-		else
-			@warn "Skipping layer $(typeof(layer)) in geometry collection"
-		end
-	end
-
-	for comp in design.components
-		for layer in comp.conductor_group.layers
-			push_layer!(layer, x_offset, y_offset)
-		end
-		for layer in comp.insulator_group.layers
-			push_layer!(layer, x_offset, y_offset)
-		end
-	end
-	return out
-end
-
-function dump_layer_geometry_csv(design::CableDesign;
-	x_offset::Real = 0.0,
-	y_offset::Real = 0.0,
-	path::AbstractString = "layer_dump.csv",
-)
-	layer_no = Ref(0)
-
-	function emit!(io, layer, x0, y0)
-		if layer isa CircStrands
-			layer_no[] += 1
-			rwire = to_nominal(layer.radius_wire)
-			lay_r = layer.num_wires == 1 ? 0.0 : to_nominal(layer.radius_in)
-			coords = calc_circstrands_coords(layer.num_wires, rwire, lay_r; C = (x0, y0))
-			for (i, (x, y)) in enumerate(coords)
-				println(io, "$(layer_no[]),circstrands,$i,$x,$y,$rwire,,")
-			end
-		elseif layer isa Strip || layer isa Tubular || layer isa Semicon ||
-			   layer isa Insulator
-			layer_no[] += 1
-			rin = to_nominal(layer.radius_in);
-			rex = to_nominal(layer.radius_ext)
-			typ = lowercase(string(nameof(typeof(layer))))
-			println(io, "$(layer_no[]),$typ,1,$x0,$y0,$rex,$rin,$rex")
-		elseif layer isa ConductorGroup
-			for sub in layer.layers
-				emit!(io, sub, x0, y0)
-			end
-		else
-			@warn "Skipping layer $(typeof(layer)) in CSV dump"
-		end
-	end
-
-	open(path, "w") do io
-		println(io, "layer,type,element,x0,y0,radius_out,rin,rext")
-		for comp in design.components
-			for layer in comp.conductor_group.layers
-				emit!(io, layer, x_offset, y_offset)
-			end
-			for layer in comp.insulator_group.layers
-				emit!(io, layer, x_offset, y_offset)
-			end
-		end
-	end
-	return path
-end
-
-
-# dump_layer_geometry_csv(cable_design; path = fullfile("cable_layers.csv"))
-
-
-all_wires = collect_layer_geometry(cable_design)
-
-using Gmsh: gmsh
-
-gmsh.initialize()
-
-system_id = "test"
-gmsh.model.add(system_id)
-
-gmsh.option.set_number("General.InitialModule", 0)
-gmsh.option.set_string("General.DefaultFileName", system_id * ".geo")
-
-# Define verbosity level
-gmsh.option.set_number("General.Verbosity", 1)
-
-# Set OCC model healing options
-gmsh.option.set_number("Geometry.AutoCoherence", 1)
-gmsh.option.set_number("Geometry.OCCFixDegenerated", 1)
-gmsh.option.set_number("Geometry.OCCFixSmallEdges", 1)
-gmsh.option.set_number("Geometry.OCCFixSmallFaces", 1)
-gmsh.option.set_number("Geometry.OCCSewFaces", 1)
-gmsh.option.set_number("Geometry.OCCMakeSolids", 1)
-
-# Log settings based on verbosity
-@info "Initialized Gmsh model: $system_id"
-
-function my_disk!(x, y, r, lcp, wires, llwires)
-	cen = gmsh.model.geo.addPoint(x, y, 0.0, lcp)
-	p1 = gmsh.model.geo.addPoint(x+r, y, 0.0, lcp)
-	p2 = gmsh.model.geo.addPoint(x, y+r, 0.0, lcp)
-	p3 = gmsh.model.geo.addPoint(x-r, y, 0, lcp)
-	p4 = gmsh.model.geo.addPoint(x, y-r, 0, lcp)
-
-	c1 = gmsh.model.geo.addCircleArc(p1, cen, p2)
-	c2 = gmsh.model.geo.addCircleArc(p2, cen, p3)
-	c3 = gmsh.model.geo.addCircleArc(p3, cen, p4)
-	c4 = gmsh.model.geo.addCircleArc(p4, cen, p1)
-
-	ll = gmsh.model.geo.addCurveLoop([c1, c2, c3, c4])
-	s = gmsh.model.geo.addPlaneSurface([ll])
-	push!(wires, s)
-	push!(llwires, ll)
-	return s
-end
-
-lc = 0.0
-z0 = 0.0
-i=1
-wires = []
-llwires = []
-for wire in all_wires
-	x = wire.x
-	y = wire.y
-	r = wire.radius_out
-	my_disk!(x, y, r, lc, wires, llwires)
-end
-
-gmsh.model.geo.synchronize()
-
-gmsh.option.set_number("Geometry.SurfaceLabels", 0)  # Show surface labels
-gmsh.option.set_number("Geometry.PointNumbers", 0)
-gmsh.option.set_number("Geometry.CurveNumbers", 0)
-gmsh.option.set_number("Geometry.SurfaceNumbers", 0)
-gmsh.option.set_number("Geometry.NumSubEdges", 160)
-gmsh.option.set_number("Geometry.Points", 1)
-gmsh.option.set_number("Geometry.Curves", 1)
-gmsh.option.set_number("Geometry.Surfaces", 0)
-gmsh.option.set_number("Mesh.ColorCarousel", 2)  # Colors by physical group
-gmsh.option.set_number("Mesh.LineWidth", 1)
-gmsh.option.set_number("Mesh.SurfaceFaces", 1)
-
-gmsh.fltk.initialize()
-
-# Define event check function
-function check_for_event()
-	action = gmsh.onelab.get_string("ONELAB/Action")
-	if length(action) > 0 && action[1] == "check"
-		gmsh.onelab.set_string("ONELAB/Action", [""])
-		@debug "UI interaction detected"
-		gmsh.graphics.draw()
-	end
-	return true
-end
-
-# Wait for user to close the window
-while gmsh.fltk.is_available() == 1 && check_for_event()
-	gmsh.fltk.wait()
-end
-
-gmsh.finalize()
+# Or the corresponding lumped circuit quantities
+per_km(p012, 1; mode = :RLCG, tol = 1e-9)

src/datamodel/circstrands.jl

@@ -83,7 +83,7 @@ function CircStrands(
 	rho = material_props.rho
 	T0 = material_props.T0
 	alpha = material_props.alpha
-	radius_ext = num_wires == 1 ? radius_wire : radius_in + 2 * radius_wire
+	radius_ext = num_wires == 1 ? radius_wire : radius_in + 2 * radius_wire # TODO: The resolved outer radius for stranded cores should account for compression. See rectstrands.jl for an example of area-preserving expansion. 
 
 	mean_diameter, pitch_length, overlength = calc_helical_params(
 		radius_in,