import gmsh import math import sys def create_spring_hex_mesh(wire_diam, major_diam, active_pitch, active_turns=4.0, dead_turns=1.0, ground_ends=True, mesh_size=0.6): """ Generates a structured FEA hex mesh of a spring by unrolling an exact arc-length distribution matrix. This guarantees a perfectly smooth helix. """ # ========================================================================= # STEP 1: COMPUTE TRUE PATH ARC-LENGTH AND PROFILE ARRAYS # ========================================================================= gmsh.initialize() wire_radius = wire_diam / 2.0 mean_radius = (major_diam - wire_diam) / 2.0 total_turns = active_turns + (2.0 * dead_turns) steps_per_turn = 40 total_steps = int(total_turns * steps_per_turn) # Pre-calculate coordinates and cumulative arc length along the centerline path path_points = [] arc_lengths = [0.0] pitch_profile = [] current_z = 0.0 last_x = mean_radius last_y = 0.0 last_z = 0.0 for i in range(total_steps + 1): turn_pos = i / steps_per_turn angle = 2.0 * math.pi * turn_pos if turn_pos <= dead_turns: local_pitch = wire_diam * 1.02 # Dead coils touch flat elif turn_pos <= (dead_turns + active_turns): local_pitch = active_pitch # Active working pitch else: local_pitch = wire_diam * 1.02 if i > 0: current_z += (local_pitch / steps_per_turn) cx = mean_radius * math.cos(angle) cy = mean_radius * math.sin(angle) cz = current_z path_points.append((cx, cy, cz, angle, local_pitch)) if i > 0: # Calculate 3D Euclidean distance step ds = math.sqrt((cx - last_x)**2 + (cy - last_y)**2 + (cz - last_z)**2) arc_lengths.append(arc_lengths[-1] + ds) pitch_profile.append(local_pitch) last_x, last_y, last_z = cx, cy, cz pitch_profile.append(pitch_profile[-1]) # Match array bounds total_wire_arc_length = arc_lengths[-1] # ========================================================================= # STEP 2: GENERATE BASE CYLINDER MATCHING TOTAL ARC LENGTH (MODEL 1) # ========================================================================= gmsh.model.add("cylinder_template") p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size) p1 = gmsh.model.geo.addPoint(wire_radius, 0.0, 0.0, mesh_size) p2 = gmsh.model.geo.addPoint(0.0, wire_radius, 0.0, mesh_size) p3 = gmsh.model.geo.addPoint(-wire_radius, 0.0, 0.0, mesh_size) p4 = gmsh.model.geo.addPoint(0.0, -wire_radius, 0.0, mesh_size) arc1 = gmsh.model.geo.addCircleArc(p1, p0, p2) arc2 = gmsh.model.geo.addCircleArc(p2, p0, p3) arc3 = gmsh.model.geo.addCircleArc(p3, p0, p4) arc4 = gmsh.model.geo.addCircleArc(p4, p0, p1) curve_loop = gmsh.model.geo.addCurveLoop([arc1, arc2, arc3, arc4]) base_face = gmsh.model.geo.addPlaneSurface([curve_loop]) gmsh.model.geo.mesh.setRecombine(2, base_face) gmsh.model.geo.synchronize() # Extrude exactly matching the calculated 3D path distance length gmsh.model.geo.extrude([(2, base_face)], 0.0, 0.0, total_wire_arc_length, numElements=[total_steps], recombine=True) gmsh.model.geo.synchronize() gmsh.model.mesh.generate(3) node_tags, coords, _ = gmsh.model.mesh.getNodes(-1, -1) elem_types, elem_tags, elem_node_tags = gmsh.model.mesh.getElements(3, 1) # ========================================================================= # STEP 3: MAP DISCRETE HELICAL SPRING MESH (MODEL 2) # ========================================================================= gmsh.model.add("helical_spring_mesh") discrete_vol_tag = gmsh.model.addDiscreteEntity(3) print(f"\nInterpolating Node Matrix along smooth true-arc length path...") new_coords = [] for i in range(len(node_tags)): lx = coords[i * 3] ly = coords[i * 3 + 1] lz = coords[i * 3 + 2] # Target node positioning relative to absolute path distance travel target_s = lz # Linear interpolation to find where this node drops onto the true helix centerline idx = 0 while idx < len(arc_lengths) - 1 and arc_lengths[idx+1] < target_s: idx += 1 # Interpolation factor (0.0 to 1.0) inside the specific step segment s0, s1 = arc_lengths[idx], arc_lengths[idx+1] t = (target_s - s0) / (s1 - s0) if (s1 - s0) > 0 else 0.0 # Extract starting and ending boundary properties of the discrete step cx0, cy0, cz0, angle0, p0 = path_points[idx] cx1, cy1, cz1, angle1, p1 = path_points[idx+1] # Smoothly blended properties cx = cx0 + t * (cx1 - cx0) cy = cy0 + t * (cy1 - cy0) cz = cz0 + t * (cz1 - cz0) angle = angle0 + t * (angle1 - angle0) local_pitch = p0 + t * (p1 - p0) current_turn_pos = angle / (2.0 * math.pi) # Ground Ends Tapering Multiplier local_scale = 1.0 if ground_ends: if current_turn_pos <= dead_turns: local_scale = max(0.15, current_turn_pos / dead_turns) elif current_turn_pos >= (total_turns - dead_turns): local_scale = max(0.15, (total_turns - current_turn_pos) / dead_turns) scaled_lx = lx * local_scale scaled_ly = ly * local_scale # --- MOVING REFERENCE FRAME CALCULATIONS --- # 1. Tangent vector components (derivative of path) tx = -mean_radius * math.sin(angle) ty = mean_radius * math.cos(angle) tz = local_pitch / (2.0 * math.pi) t_len = math.sqrt(tx**2 + ty**2 + tz**2) tx, ty, tz = tx/t_len, ty/t_len, tz/t_len # 2. Normal vector components (pointing toward the central Z axis) nx = -math.cos(angle) ny = -math.sin(angle) nz = 0.0 # 3. Binormal vector components (cross product of Tangent and Normal) bx = ty * nz - tz * ny by = tz * nx - tx * nz bz = tx * ny - ty * nx # Map local coordinates onto the vector frame orientation gx = cx + (scaled_lx * nx) + (scaled_ly * bx) gy = cy + (scaled_lx * ny) + (scaled_ly * by) gz = cz + (scaled_lx * nz) + (scaled_ly * bz) new_coords.extend([gx, gy, gz]) gmsh.model.mesh.addNodes(3, discrete_vol_tag, node_tags, new_coords) gmsh.model.mesh.addElements(3, discrete_vol_tag, elem_types, elem_tags, elem_node_tags) gmsh.model.addPhysicalGroup(3, [discrete_vol_tag], name="SPRING_VOLUME") gmsh.write("ground_spring_mesh.msh") gmsh.write("ground_spring_mesh.inp") print("[SUCCESS] Spring generated with a perfectly smooth round wire.") if "-nopopup" not in sys.argv: gmsh.fltk.run() gmsh.finalize() if __name__ == "__main__": create_spring_hex_mesh( wire_diam=1.12345, major_diam=15.0123, active_pitch=10.98765, active_turns=4.0, dead_turns=1.0, ground_ends=False, mesh_size=0.3 )