I’ve been looking for some good simulations of spring-mass networks with customizable parameters. Those I find are usually fairly simple, with uniform masses and springs, each connected at right angles to adjacent masses (i.e. six connections from each face in a cube setting, which is the most complex simulation I could find). I would like to explore the dynamics of a less uniform system, i.e. the effect of a large mass with many connections to smaller masses, how those smaller masses respond if they are more interconnected, etc. So I was wondering if anyone has some suggestions for a good simulation tool that they may have stumbled upon during some coursework or similar? I am in particular interested in the visual representation of the time evolution of such a system, so some kind of tool that generates an animation would be ideal.
Is my best bet perhaps to use the physics engine of e.g. the Godot game engine (Jolt physics) or whatever solver is used in Blender?
Oh, I read it wrong before, but basically the same principle of my previous answer, but the differential equation between dots, so there is a direction unit vector associated with the two connected points, the differential equation would need the spring length too, but that would be easy to put, as well as the constants.
The old wrong answer:
You can model this grid pattern using partial differential equations, the springs + masses setup reduces to a 2D wave equation, your varying masses and spring constants reduce to a same order and size matrix of constants, instead a global value, the values are on the grid. With finite differences, it’s very easy to solve.
It will be just a grid updating, you can build this in less than 100 lines of code:
import numpy as np def free_wave_equation_2d(lenx, leny, sol_len, c, ew, ew_len, dt): """ Parameters: lenx : int - Number of points in x direction leny : int - Number of points in y direction sol_len : int - Number of time steps c : ndarray (2D): Wave speed coefficient array ew : ndarray (1D): Source/initial condition array ew_len : int - Length of the ew array dt : float - Time step Returns: array : ndarray (3D) - Resulting solution (lenx, leny, sol_len) """ array = np.zeros((lenx, leny, sol_len)) acceleration = np.zeros((lenx, leny)) half_lenx = lenx // 2 half_leny = leny // 2 # Set initial condition at t=0 array[half_lenx, half_leny, 0] = ew[0] for i in range(1, sol_len): # u_curr represents u(t) u_curr = array[:, :, i - 1] acceleration[:, :] = 0 # --- Finite Difference Laplacian Calculation --- # X-direction acceleration[0, :] = u_curr[1, :] - u_curr[0, :] acceleration[-1, :] = u_curr[-2, :] - u_curr[-1, :] acceleration[1:-1, :] = u_curr[2:, :] + u_curr[:-2, :] - 2 * u_curr[1:-1, :] # Y-direction acceleration[:, 0] += u_curr[:, 1] - u_curr[:, 0] acceleration[:, -1] += u_curr[:, -2] - u_curr[:, -1] acceleration[:, 1:-1] += ( u_curr[:, 2:, :] + u_curr[:, :-2, :] - 2 * u_curr[:, 1:-1, :] ) # Apply wave speed coefficient acceleration = acceleration * (c**2) # --- Finite Difference Integration Step --- if i == 1: # For the very first step array[:, :, i] = u_curr + 0.5 * acceleration * (dt**2) else: # Standard Leapfrog / Central Difference Integration: u_prev = array[:, :, i - 2] array[:, :, i] = 2 * u_curr - u_prev + acceleration * (dt**2) if i < ew_len: array[half_lenx, half_leny, i] = ew[i] return arrayEdit: The code was ported from my college time fortran code.