Mass Spring Models of Amorphous Solids

: In this paper we analyse static properties of mass spring models (MSMs) with the focus of modelling non crystalline materials, and explore basic improvements, which can be made to MSMs with disordered point placement. Presented techniques address the problem of high variance of MSM properties which occur due to randomised nature of point distribution. The focus is placed on tuning spring parameters in a way which would compensate for local non-uniformity of point and spring density. We demonstrate that a simple force balancing algorithm can improve properties of the MSM on a global scale, while a more detailed stress distribution analysis is needed to achieve local scale improvements. Considered MSMs are three dimensional.


Introduction
Mass spring models are used for the purpose of simulating elastic objects in various fields. Their potential applications include computer games, animations, virtual reality environments as well as engineering of structures or materials in which deformation under stress is considered, crack propagation studies or human tissue simulations. Related fields range from ultramicroscopy to astrophysics [1][2][3][4][5][6]. Specific needs of specific application may vary, but the general theory of how mass spring models (MSM) deform concerns all of them.
We distinguish two broad classes of mass-spring networks-crystal (lattice) based and disordered. In crystal based networks, mass points are placed on a periodic lattice and mechanical properties of such networks can be expressed with analytical formulas. This allows for more precise analysis and description of their behaviour. In disordered networks, on the other hand, the mass placement rules are relaxed and do not have to follow any regular order ( Figure 1). The downside of such models is reduced accuracy and increased difficulty of estimating their mechanical properties [7]. Consequently, lattice based networks are usually preferred over disordered ones and are even used to model materials, which are known not to be crystals. For numerous applications this is not a problem, however, there is a danger that in certain situations the material may inherit some unexpected properties of the periodic network and exhibit unintended behaviours. An example is the crack propagation problem, where geometry of the network may influence the observed crack patterns considerably. If the numerical model has a regular, lattice-like topology which determines the potential crack placement, it creates "easy propagation planes" for cracks (whether it is based on MSM or FEM, or some other modelling technique). In such models cracks tend to form and propagate along lattice dependent directions and resulting patterns may reflect lattice properties instead of material ones ( Figure 2). In such cases it is difficult to asses, if the obtained results follow from the mathematical model, or if they are only artefacts of the numerical representation. For example, in a model with cubic lattice topology a curved crack may appear as two straight segments approaching each other at 45 degrees, or worse, the crack may not bend at all, if breaking of subsequent segments influences tensions within the system. Some attempts have been made to address these problems, such as Chen et al. [8], where introduction of non-local interactions into the model is shown to reduce negative effects of a lattice topology. Such solutions may mitigate most apparent artefacts in certain situations, however they do not eliminate the problem completely. A class of conditions for which such methods are unlikely to work involves situations where cracking is induced by a shrinkage front progressing through a material, such as in the case of drying or cooling materials [9][10][11][12][13]. For example if a cooling front is advancing through solidifying lava in presence of water, the problem has a translational symmetry. Cracks which appear at the top, allow for the water to be in contact with the top of the front, keeping the temperature at 100 • C, while the bottom of the front is in contact with non-solidified lava at an approximately constant temperature T L . The front itself is expected to have a linear gradient of temperature between these two values [14]. Modelling such process with lattice-based networks will result with exactly the same conditions at each cracking step, and since the tip of the crack can follow only discretised paths, it may never turn if the incentive to turn is too weak. In such situations it may be beneficial to use disordered MSMs instead of lattice-based ones. As mentioned, the downside of such models is the reduced accuracy when compared to lattice-based MSMs, which in turn may cause discrepancies between modelled and theoretical properties of the material e.g., such as reported in [15]. In this work we place focus on disordered MSMs.
Mass spring models studied in this work are based on [7,16] and the reader is encouraged to refer to these papers for a more detailed introduction. We achieve varying values of Poisson's ratio ν in our models, by introducing a force dissipation mechanism. This technique offers a considerable simplicity in terms of implementation as well as low conceptual and computational complexity, when compared to alternative methods capable of achieving full spectrum of ν. A comprehensive review of other MSM-like techniques can be found in [17].
In the following sections we show how to formalise the force dissipation mechanism proposed in [16], with a focus on ease of implementation. We introduce a new pair of elastic constants for the material, which translate into simple properties of springs and nodes of the network. We also demonstrate how to mitigate the accuracy problems of disordered MSMs by adjusting stiffness coefficients of the springs. Presented techniques allow to improve both global and local behaviour of disordered mass-spring networks.
It should be noted, that we investigate only the static properties of MSMs. We do not address here the question of how these systems evolve over time or what are the most efficient numerical schemes to track their dynamics. Efficient ways of simulating the dynamics can be found e.g., in [18,19].

Elastic Moduli
In linear elasticity, the relation between tensions and deformations of a continuous medium is linear. This can be expressed in general form by: whereσ denotes the stress andˆ the strain tensor;Λ is a tensor of elastic constants. Components ofΛ are tied together by various relations which follow from symmetry ofσ andˆ . In this paper we discuss homogeneous isotropic media, which further restricts the values of Λ ijkl . In a medium which behaves like fluid or gas, in an equilibrium state, all the "tensions" are distributed equally in all directions, that is σ ij = −δ ij p, where p denotes pressure (for fluids or gasses we should not be using quantities like stress or strain, however the purpose here is only to highlight analogies in the description of various media, in situations, where it is applicable to do so; same comment applies to E and ν below). More precisely σ ij = −δ ij dp if we consider a small deformation from a reference configuration, so that: where K is a bulk modulus, V volume and N is the number of dimensions of the space. In this case Lamé parameters are equal λ = K, µ = 0, Young's modulus E = 0 and Poisson's ratio ν = 1 N−1 . One realisation of such a medium is a disordered collection of particles, which collide with each other, bouncing off randomly, so that there is no way to induce a momentum flow/surface pressure higher in one direction than in others. We will say, that the interactions which lead to such relations, are of a completely dispersive nature.
Let us now take our ensemble of particles and impose fixed relative positions on them, so that deformations are properly expressed by the strain tensor. We will assume here, that a change in the distance r between a pair of particles gives rise to a force (acting between them) which can be represented by a central potential V(r). In such setup there is no way to displace a chosen particle and induce tension in only one given direction. If we move a particle by dl along x axis, y components of the distances to its neighbours will be affected as well. This is a geometrical property of euclidian space and it gives rise to the elastic moduli tensor of the from [7,17,20]: There is, similarly as in the case of gas/fluid, only one elastic constant here, however, this time bulk modulus depends on the connectivity, which follows from the dimensionality (the higher the dimension, the more particles are present in the immediate neighbourhood). It is given by K = µ(1 + 2 N ). Similarly E = 2µ N+2 N+1 and ν = 1 N+1 . We also have λ = µ. This one elastic constant corresponds to direct interactions.
There exist materials with these properties, an example of which is a diamond, however this model is too simplistic to properly describe elastic solids in general, as it results in only one degree of freedom, corresponding to a single elastic constant, while real materials are characterised by two.
Possible solutions, which allow to achieve two degrees of freedom, include introduction of angular springs or beams into the model, usage of non central forces or introduction of a volume change component into the potential energy [17,[21][22][23][24][25] (some of them limited to very specific lattice topologies). However, a particularly simple way of extending this model, which works both with lattice based networks as well as disordered ones, is to allow for the force to be partially dispersed, so that all the interactions become superpositions between direct and dispersive forces [16].
whereF µ denotes the direct, andF * dispersive component of the force. The resulting elastic body then becomes a superposition between fluid and diamond like materials [16,17,26]. The reasoning behind this is, that direct interactions are in essence interactions by means of idealised classical Newtonian force, which acts across time and space, instantly, without accounting for relative velocities or other possible characteristics of the bodies involved in the problem, such as their shape, which could potentially influence the net effect of one body acting on another from a distance (There is plenty of evidence suggesting that the real character of the interaction on a distance depends on factors other than just relative position of the bodies in question, the prime example of which is electrodynamics, where the presence of relative velocity gives rise to a magnetic force).
If we take: then Q denotes the ratio between dissipative forces and the direct ones, with its values spanning from −1 to infinity.F µ ∼ µ, andF * ∼ µQ. Other elastic constants become: In practice, when implementing a mass spring model for the representation of an elastic material, it may be more convenient to introduce yet another pair of constants.
C now can be interpreted as the total amount of interaction potential (e.g., a density of interaction carriers between two particles). The dissipative part is given by C · D, and the direct one by C · (1 − |D|). D is the dispersion fraction with values between −0.5 and 1. With this description we avoid the singularity in Q around the ν = 0.5 point and allow for expressing possible dynamic differences between media with K = 0, but varying C (the aspects of MSMs not explicitly addressed in this work). Other elastic constants become: A few examples of the decomposition of the interactions into dissipative and direct parts are illustrated in the Figure 3 (with E = const). The first extreme, ν = 0.49, corresponds to D ≈ 1, where almost all the forces are dispersed, making the material fluid-like (ν = 0.5 in 3D means no change in volume under unidirectional compression; the limit value of ν = 0.5 is unstable). The second extreme, ν = −0.99 has D ≈ −0.5 (ν = −1 once again is an unstable limit), which means that half of the compressing force acts in a dispersive manner, but with a "reversed sign". In such regime K = 0 and the material can be uniformly compressed without changing its elastic energy. As we can see, the unidirectional compression gives rise to considerably high forces, both direct and dispersive, which cancel out each other almost perfectly in the resulting concave shape. The ν = 0.25 has no dispersive forces at all and for ν = 0, we have D = −3/8.

Mass Spring Models
Mass spring models represent elastic solids with discrete points connected by springs. As mentioned in Section 1, we distinguish two broad classes here-a regular crystal-like models and disordered ones [7]. In both cases, elastic constants of the material follow from spring coefficients. As elucidated in the previous section, if classical, direct force springs are used, the value of Poisson's ratio for three dimensional isotropic and homogeneous models is fixed and equal to 1 4 . The relation between bulk modulus K and spring parameters is given by: where L i and k i denote the length of i-th spring and its stiffness coefficient. V is the volume of the object and the sum is taken over all springs inside this volume. In order to model materials with other values of Poisson's ratio, we can introduce the force dissipation mechanism discussed in previous sections with the C modulus obeying: where the classical stiffness coefficient of a spring is now κ i = k i (1 − |D|). And the force propagation follows the algorithm [16]: For each node: In the second pass: • Redistribute the accumulated force as to all the springs connected with the node (by applying it on both nodes the spring connects), where b denotes the number of these springs Estimates given by Equations (14) and (15) can be applied to any mass spring network, however the more homogeneous and isotropic the network is, the better it will approximate elastic properties of a given material.

Accuracy of Random MSM Models
In this study we consider MSM networks with randomly generated points connected by springs with their nearest neighbours. The network topology is characterised by two parameters-first, minL, gives a minimal allowed distance between any pair of MSM nodes; second, maxL, is the range within which spring connections are formed (any two nodes which are less than maxL apart get connected by a spring) [7]. Please note, that the random point generation mechanism described in [7] is flawed. The article suggests a performance improvement for the implementation, which advocates a usage of a small spherical brush-like window traveling through the material in order to fill it with random points. The performance improvement for point collision calculation is real, however the "window" should be moved after each point addition, in contrast to a procedure which fills the window region fully first, before advancing to the next region of the material. The latter may introduce a non-homogeneity at the window borders. This was not stated clearly in the original article.
For such MSMs, our previous experiments show, that as long as the average number of springs attached to one node is sufficiently high, Equations (12) and (14) give the elastic properties of the network in bulk with a reasonable accuracy. The values of Poisson's ratio for such models are very close to theoretical predictions (observed discrepancies were in many cases of magnitude of measurement errors) and the values of E diverge from theoretical by no more than 10% [7,16]. Such accuracy may be satisfactory for many applications, but there is certainly a room for improvement.
Currently we will focus not only on bulk properties, but also investigate how MSMs behave locally, on a scale of internode distances. We perform a simple test, in which we consider a model of a cube 10a 0 × 10a 0 × 10a 0 , that undergoes compression (expansion) uniformly in all directions and observe how the MSM reacts. The following setup was used for this test. The elastic constants of the material were set to K = 100 k 0 a 0 and ν = 0.25 (this means D = 0, which makes the theoretical analysis simpler. Materials with D = 0 are expected to inherit the same properties here, which was confirmed by additional tests afterwards). The material was modelled with a random MSM with minL = 0.92 · 0.5a 0 and maxL = 1.77 · 0.5a 0 , which gives around 6700 nodes, 55,000 springs and approximately S = 18 springs connected to one node inside of the material. All springs were assigned the same spring constant k[k 0 ] calculated with Equation (14). The cube was expanded uniformly in all directions by 1% (We will use expansion and compression interchangeably in the rest of the article when discussing the phenomenon analysed in this experiment; for technical reasons expansion is easier to perform). In such setup the values of diagonal components of the strain tensor should be equal to 0.01 and corresponding diagonal elements of the stress tensor should be equal to 3 k 0 a 0 (with a 0 being the unit of length and k 0 the unit of force).
In our numerical experiments we measure the stress inside of MSM according to the procedure proposed by Hardy [27,28]. It estimates the stress by measuring the tension inside a probe sphere placed in the point of interest. The tension in each spring in the surroundings of the point is weighted by a bond function. If only a fraction of the spring is contained within the probe sphere, the corresponding fraction of the tension is taken into account:σ where: x αβ -distance between nodes α and β F ab -force acting through the spring ab B αβ -bond value In order to calculate the bond value, we associate an influence sphere with each spring-the center of the sphere is placed in the middle of the spring, and its radius is equal to half of the length of the spring. Then the bond value is calculated as a fraction of the volume of that sphere, intersecting with a probing region. This is a slight modification to the original method, which measured intersections between raw springs and a probing region. Our experiments shown that calculating intersections with three dimensional spheres surrounding the springs reduces noice and allows for more localised measurements (smaller probe regions).
Experiments presented in this article are using spherical probing region with the radius 1.5 L , where L denotes the average length of a spring in MSM.
The test consists of two steps. First we scale the cube uniformly, mimicking the perfect uniform deformation, then we let the MSM relax (holding only borders in place). In both stages we measure stress and strain present inside of the cube. In a perfect model the relaxation phase should not lead to any additional displacements.
As an example the 22 strain tensor component and the σ 22 stress tensor component measured in this experiment are presented on Figures 4 and 5. The character of discrepancies in other tensor components is similar. The measurements were taken along the line in X direction in the middle of the cube, with regions in close proximity to the border omitted (to avoid various problems with stress measurements near the border). Standard deviation from theoretical value of the final (relaxed) stress is denoted by ∆ t , and the standard deviation from the average value by ∆. Standard deviation between the stress in compressed (scaled) state σ s and the stress in relaxed state σ r is denoted by ∆ r .
where A denotes measured, A average and A 0 theoretical value of a strain or stress component; M is a number of measurements.   As we can see, the plotted tensor components average out to match theoretical predictions reasonably well. Local divergences, however, are high. The root cause of this behaviour is that the random MSM achieves isotropy statistically in bulk, but is not isotropic on singular nodes. To understand what exactly we mean by that, let us look at a simplified, one dimensional example from Figure 6. In disordered MSMs lengths of the springs are not equal. In uniform compression ∆L/L 0 = and if all the springs have the same k, the force acting on a spring is F A = −kL A , and as a result springs do not compress by the same degree under a constant pressure. If they did, the red coloured node from Figure 6 would experience non-zero net force and would move. In 3D networks, additionally, the number of springs facing each direction may differ, which further influences local isotropy. These are the reasons for local discrepancies in deformation as well as ∆ r for stress in the material from Figures 4 and 5. Forces on the nodes do not balance to zero when the springs get compressed by the same degree. This also affects the estimation of C (and consequently K and E) from Equations (14) and (15). These equations measure the energy needed to compress material uniformly, but in this case the energy is lower, and the corresponding K is lower as well. The values of σ 22 average to the desired 3 k 0 a 0 when the MSM is scaled uniformly (indicated by dotted line on the Figure 5), however after the relaxation this value drops to around 2.8 k 0 a 0 .

MSM Tuning
Our goal is to adjust stiffness coefficients for springs in such way, that in an MSM subjected to uniform compression (realized by simple scaling): (a) forces acting on inner MSM nodes sum to zero, (b) stress is isotropic.
Achieving (a) would reduce both, the variations of strain, and the change in stress which occurs when a scaled model relaxes (∆ r ), and in consequence would allow for a better estimation of the actual K in bulk. Achieving (b) would reduce the variations of stress. It should be noted, that condition (a) does not apply to the nodes which lie on the border of the MSM, where forces will sum to a value that counteracts the outside pressure. To distinguish which nodes are "inner", and which lie on the border, for each node, we measure the intersection between modelled object and a sphere centered on the node, with a radius equal to maxL. If the whole sphere intersects we treat the node as inner. Otherwise we assume it lies on the border.
To give some notion about how well (a) holds for a chosen MSM, we introduce the following measure, which we call kL-score: where the inner sums indexed by b are taken over all springs connected to a node n. Outer sum is taken over all the nodes, the total number of nodes is N, andL b is the vector between nodes connected by a spring b; |x| is the length of a vectorx. The lower the value of S kL is, the closer to zero is the sum of forces acting on an MSM node (the lower the better). Such inner kl-score of the model from Figures 4 and 5 is equal S kL = 0.11.

Constant kL
Since we have identified that in uniform compression the force exerted by a spring equals F = −kL , our first attempt at improving accuracy of random MSMs could be by setting the stiffness coefficients of the springs in such way that kL = const. This would certainly help in 1-d case from Figure 6, however for three dimensional MSMs effects of such modification should be rather limited. In 3D the number of springs facing each direction is of bigger importance than differences in their lengths.
Setting kL = const does however improve the MSM slightly with S kL dropping to 0.091 and divergences from Figures 4 and 5 becoming: This actually results, right away, in improvement for the estimation of the bulk value of K by a few percent (more results in following sections).

kL-Tuning
As described in previous sections, the source of divergences in strain within MSM are unbalanced forces which appear during uniform compression. They cause the nodes to move and springs relax to new equilibrium lengths.
Our first goal is to construct an MSM in which the unbalanced forces do not appear. It can be achieved in a number of ways, one of which is simply by solving as a set of linear equations for stiffness coefficients k i . As the number of springs is an order of magnitude larger than the number of nodes, this may turn out to be problematic (but not impossible) for large systems, unless we have a whole computational cluster at our disposal. An alternative way, is to use the original MSM (with either k = const or kL = const) and follow the same relaxation path which our MSM from Figure 5 traveled when reaching the equilibrium state. If we consider a simple time integrator (e.g., Euler or Verlet scheme) its basic steps are: • compute forces • update velocities • update positions the "update positions" step will result in changes in lengths of springs and consequently in changes of forces exerted by these springs. We can achieve the same change in forces by changing stiffness coefficients of the springs k i instead of moving the nodes. Such algorithm (with damping) will find an MSM with S kL ≈ 0. This is the approach we use, however in our case we skip the "update velocities" step and modify spring coefficients directly, increasing or decreasing their value depending on the magnitude of the projection of the force vector onto the direction of the spring. This is analogous to a simulation of an overdamped movement and basically is a variation of steepest descent minimisation procedure. The tuning procedure starts by artificially changing natural lengths of all springs by the same degree and then following with a number of minimisation steps. In our implementation we use L−L 0 L 0 = = 1%. Additionally appropriate border conditions need to be set, similarly as it was done in compression/expansion experiment, where borders were frozen in place. This basically means that the relaxation should only be applied to inner nodes. Nodes that lie near the border have by definition non-isotropic spring connections and achieving S kL = 0 is not even possible for them.
Each iteration of the algorithm reduces S kL and increases the MSM quality. We stop when a certain value of S kL is reached (0.0003), or the progress becomes too slow (max relative change of k smaller than 0.002).
Additionally, we place a restriction on the maximal allowed change of k for each spring to 50% of its original value. This way we explore only the set of solutions close to our original choice of k, and avoid degenerate solutions with negative k. In the last step we restore natural lengths of springs to their original values.
The values of 22 and σ 22 for the resulting MSM of the cube test are presented on Figures 7 and 8. The divergence in strain is now practically zero and ∆ r is two orders of magnitude smaller than for k = const. The profile of stress in this case is exactly the same as the relaxed profile from k-const model ( Figure 5), but shifted towards desired value of 3 k 0 a 0 . This means that we have managed to improve the estimate of bulk K, however local divergences did not change.   Figure 8 shows, that S kL ≈ 0 does not necessary translate into isotropy on a node basis. To further improve MSM representations we need to make sure, that not only forces on singular nodes sum up to zero, but also that their magnitude is approximately the same in all directions. Figure 9 illustrates an MSM node for which this is not the case; in this example forces do sum to zero, but their magnitude is different in X and Y directions. In kL-tuning at each step of the algorithm we modified k i based on the current value of the force acting on each node. In the improved procedure, instead of calculating the force, we will calculate the stress tensor at the point where the node is placed and see how much it diverges from the expected tensor:

ikL-Tuning
where δ ij is the Kronecker's delta. The contribution of each spring toσ can be expressed as: where is the relative length change of the spring, L the current length and L 0 the neutral length. V is the volume of the region in which we are measuring the stress and B is a bond between this region and the spring. The bond can be calculated in a number of ways [27,28]; in our case we use the percentage of the overlap between probing region, which is a sphere and another sphere placed in the middle of the spring with radius equal 0.5L. As we remember L 0 is artificially changed in the first step of the tuning procedure. We project the spring influence onto theσ err and compute the ∆k for each spring as: whereÂ :B = ∑ ij A ij B ij , and t is step size for the steepest descent algorithm (t ∼ 1 b , where b is the number of springs connected to a node). The stress tensor tuning is run in addition to kL-tuning, and it also has a limit for the relative change of k on each spring to prevent degenerate solutions in which k becomes negative. After this procedure the strain tensor is as close to theoretical as it was in Figure 7. This time however we observe an improvement in stress tensor. The σ 22 component is presented on Figure 10. The divergence ∆σ 22 dropped to 0.067, and is about three times smaller than it was for kL-tune.
Comparison of the four presented algorithms is terms of divergences ∆ t , ∆ r σ and ∆ t σ observed in the cube compression test is presented on Figures 11-13. As we can see kL-tune procedure eliminates divergencies ∆ t , ∆ r σ, which are caused by internal relaxation of unbalanced forces. The ikL-tune additionally reduces local variations in stress and reduces divergencies between stress present in the deformed body and its theoretical value.

Effects on Young's Modulus
Next we investigate what are the effects of tuning, on global material properties and how reduced variations of local stress and strain translate into the value of Young's modulus E.
We measure the value of E by performing a numerical experiment similar to the one described in [7]. We apply a static displacement to a block of material with dimensions 35a 0 × 7.5a 0 × 7.5a 0 , so that it gets compressed along X axis and we estimate E by measuring forces present in the material. In [7] we have measured the force passing through YZ plane in the middle of the system, following the equation: where F is the reaction force, A the cross-sectional area of the block (in YZ plane), and ∆x is the deformation of the block along x direction. This equation however does not account for the possibility of non-zero stress in directions other than X in the system. As we have seen in previous section such stress may be present even for stress components which theoretically should be equal zero. A more accurate equation gives E as: This time we measure not only a force through yz plane, but also xy and xz planes. After applying this correction to the measurement procedure, the resulting estimates of E become smaller by 2-5%, showing that the E calculated this way may diverge from theoretical value by a higher degree than our original estimates.
The results for the four variants of random MSM are presented in the Figure 14. First we notice that simple setting kL = const instead of k = const does improve the value of E noticeably. Introduction of kL-tune gives further improvement, however ikL-tune does not. For high S it is as good as kL-tune, however for lower values it is worse. The curve seems however more stable with lower variations around its trend value.
In all the cases the lower the number of springs, the higher divergence from theoretical value we observe. The average number of springs connected to one node for the cubic lattice MSM is S = 18, and it seems that in our experiment all the curves reach plateau around this point. While the effects of tuning are positive for MSM with mid to high S , for values lower than 14 other issues related to low network connectivity seem to dominate and the tuning can no longer compensate for them. For S < 12 we observe a rapid decrease of rigidity of the network. This result is consistent with [29], where the loss of rigidity caused by reduced network connectivity is studied. Authors show, that weakened FCC networks with S around 8-10 suffer from the loss of rigidity and become floppy for lower values of S (exact numbers depend on applied strain).

Conclusions
In this article, we have demonstrated possible accuracy problems one may face when using disordered MSMs. We have proposed two tuning procedures which aim to improve the accuracy of such MSMs. The first one, kL-tuning, eliminates unbalanced forces on MSM nodes which allows for better estimates of K (and E) of the whole network. The second, ikLtuning aims to additionally reduce local stress variations. In both cases, the implementation details are of a lesser importance, than the effects each tuning procedure has on the MSM quality (as mentioned one might simply use linear solvers).
The presented analysis gives an overview of what should be expected of randomised models and what are their limitations. The proposed techniques can be used to reduce both global and local discrepancies and inhomogeneities present in the material, however one should keep in mind that such tuning may not always be necessary or even desirable. As mentioned in the introduction, real materials are not perfect, and the presence of inhomogeneities may be advantageous for certain applications and purposes. Overaggressive tuning may simply destroy the desired properties of our models (For the same reason, we did not even consider annealing-like approaches which modify position of nodes-such procedures may very well lead to node reorganisation into crystalline structures).
In the view of this, one may decide to abandon the tuning altogether and use k = const for the whole network, with k scaled to match the desired K, not based on Equation (14), but rather on experiments. Typically a given network characteristics (minL and maxL) lead to discrepancies of the same order, independent on the specific specimen of the MSM, therefore once we establish, that e.g., k should be increased by 10% relative to the original estimate, we can simply apply the 10% increase for all models generated with these parameters and get plausible results. The stress and strain figures in Section 5 should give a good intuition as to what to expect.
In our reference single threaded implementations, for the model with 55,000 springs, kL-tuning took 1s to complete, while ikL-tuning about 10s. The achieved reduction of stress variations is about 60% for the stress measured with resolution comparable to internode distances. However if such localised measurements are not needed, the kL-tune alone may be sufficient as it is considerably faster and much easier to implement.
Finally, we should remember, that problems discussed here concern disordered MSMs. If for a given application, a crystal-like structure is appropriate, then cubic lattice MSMs offer very high accuracy for typical deformation scenarios [7,30], and if we consider simple compression or shearing, the errors in local stresses and deformations are of an order of numerical precision of floating point numbers.