Analysis of disordered trusses using network Laplacians

Abstract

Truss structures, with distributed mass elements, at macro-scale are common in a number of engineering applications and are now being increasingly used at the micro-scale to construct metamaterials. In analyzing the properties of a given truss structure, it is often necessary to understand how stress waves propagate through the system and or its dynamic modes under time dependent loading so as to allow for maximally efficient use of space and material. This can be a computationally challenging task for particularly large or complex structures, with current methods requiring fine spatial discretization or evaluations of sizable matrices. Here we present a spectral method to compute the dynamics of trusses inspired by results from fluid flow networks. Our model accounts for the full dynamics of linearly elastic truss elements via a network Laplacian; a matrix object which couples the motions of the structure joints. We show that this method is equivalent to the continuum limit of linear finite element methods as well as capable of reproducing natural frequencies and modes determined by more complex and computationally costlier methods. Our results show that balls-and-springs models inadequately describe dynamics, especially at short times relative to wave propagation time through rods. Furthermore, we illustrate the method's utility in optimizing target joint displacements using impedance matching and resonance-based schemes, offering a computationally efficient approach for analyzing large, complex truss structures.

---

1 Introduction

Trusses have been a mainstay of structural engineering for a substantial amount of human history, with applications including bridges, buildings, airplanes, and spacecraft. While these applications at the scale of several tens of meters are well established, truss metamaterials with microstructures in the range of a few millimeters are of intense research interest currently due to their easy manufacturability by 3D printing and other techniques.¹ They are also relevant in understanding the properties and adaptation of biological systems ranging from cytoskeleton to bone tissues^2,3 in animal skeletons.⁴ However, efficiently analyzing the dynamic behavior of complex truss structures, especially in disordered or hierarchical configurations, remains a significant challenge in the field of soft matter physics.

Analysing wave propagation through truss metamaterials,^5,6 is particularly interesting since they can be designed to allow only certain wavelengths/frequencies to propagate, thus enabling applications in acoustic absorbers and transmitters.^7,8 Similar techniques have also been applied to structures comprised of polymer networks for the purpose of understanding response to propagating loads and defect particles.^9,10 These applications involve precise determination of the dynamic behavior of the structure. Furthermore, current imaging techniques for the measurement of local strains have become so sophisticated that it is possible to observe wavefronts propagating through individual elements at the microscale.¹¹ These experimental developments enable the validation of detailed computational models which was not possible before.

The dynamic behavior of truss structures is most commonly studied using the finite element method. Each member of the structure is discretized into appropriate truss elements, then stiffness and mass matrices are assembled, then the equations of motion (including external forces) are written in matrix form starting from Newton's laws. The natural frequencies and mode shapes (which are critical in structural design and in understanding the dispersion relations of metamaterials) are computed by solving an eigenvalue problem using the stiffness and mass matrices.¹² The frequencies obtained depend on the mass and stiffness matrices; they vary depending on how fine the discretization is, whether the consistent or lumped mass matrix (or a combination) is used, and on the shape functions used to assemble the stiffness and mass matrices. Often, very fine discretization is required to get an accurate estimate of the natural frequencies and the computational costs can be prohibitive if the structure has a large number of degrees of freedom. Methods that retain the full element dynamics have also been developed via solutions in Fourier space (reverberation matrix method)^13,14 or integral operators.¹⁵ These produce accurate solutions for large structures but require more sophisticated mathematics and analysis algorithms. While it is possible to directly simulate the temporal response of a truss structure by monitoring the propagation of reactive wavefronts,^13,16,17 this can become computationally costly as more wavefronts are produced via scattering.

In this article, we present a new method for determining the dynamic behavior of a given truss structure based on our work on fluid flow networks.¹⁸ We retain a full description of the element dynamics much like the methodology of,^13,14 but we maintain a focus on the motions of the structure joints rather than stress within the elements. By assuming linear stress and displacement dynamics, we can express the propagation of stress waves in terms of the joint motions in Fourier space and develop a linear relation between the displacements of and forces acting on the joints. This allows our resulting matrix objects to be smaller in size and thus more computationally efficient to analyze in comparison to those produced via element stress calculations. This method is different from the spectral element methods^19–23 which are similar to finite element methods except that the shape functions are chosen to be orthogonal functions as in a Fourier basis. Thus, if the bars are discretized using spectral elements then the size of the problem scales with the number of bars. In the technique presented here the size of the problem scales with the number of joints.

We demonstrate the method's accuracy and efficiency through analysis of a simple structure (square frame with crossbar) and extend its application to disordered networks, revealing the inadequacy of ball-and-spring models at high frequencies. In contrast to the balls and spring method, our spectral method gives consistent results across frequency ranges. Finally, we show that our method can be used to optimize target joint motion in response to input joint oscillations, illustrating its usefulness in designing dynamics of complex, disordered systems.

2 Methods

2.1 Single rod dynamics

We begin by considering a solid, cylindrical rod of length L and cross sectional area A. As the rod undergoes tension and compression in response to externally applied forces, each infinitesimally thin segment is displaced from its rest position by an amount u(z,t), where z is the distance from the z = 0 end of the rod in its unperturbed configuration and t is time. Specifically, u(z,t) > 0 implies the segment has shifted in the direction of increasing z and u(z,t) < 0 implies a shift in the direction of decreasing z, as depicted in Fig. 1(A).

Fig. 1 (A) Single rods undergo position dependent displacement, u(z,t), which causes strain, ε(z,t), within the material. This results in a buildup of stress, σ(z,t), which resists the strain. (B) Multiple rods can be connected into a network via joints. Each joint can move with velocity w⃑_μ(t) and induce similar velocity within the individual elements (blue arrows). Similarly, a force, P⃑_μ(t), can be applied to the joint but will be resisted by the elements (red arrows).

From this displacement field we can derive two other critically important fields; the velocity field, v(z,t), and strain field, ε(z,t). Defining the sign of v and ε to denote movement in the direction of z and tensile expansion of the material respectively allows for the relations

Under the assumption that the material is linearly elastic, the stress field, σ(z,t), can be expressed simply as σ(z,t) = Eε(z,t), where E is the Young's modulus. This forces the sign of σ to be such that σ > 0 represents tensile stress while σ < 0 represents compressive stress. These definitions of v, ε, and σ immediately allow for the identity

We now consider the forces acting on an infinitesimally thin segment of the rod, such as that of width δz shown in Fig. 1(A). Following our sign convention for the stress field, Newton's second law takes the form

where ρ is the material density and the δz → 0 has been taken. Together, eqn (2) and (3) represent the dynamic coupling between the stress and velocity fields of the rod.^13,14 Additionally, we can express these in terms of the displacement field to transform eqn (3) intothus producing the simple wave equation wherein waves can propagate through the displacement field at speed .

One important implication of eqn (4) is that u(z,t) can be expanded into its Fourier modes via

where F(ω) and B(ω) are the forward and backward wave amplitudes respectively and the integral being over all real ω is implied. This in turn allows for the velocity and stress fields to be expressed aswhere is the material impedance.

2.2 Rod and joint network construction

We now consider a network comprised of rods obeying the dynamic equations outlined thus far connected via a series of joints such as the one depicted in Fig. 1(B). Here we will use the index μ to denote a particular joint and the index ν to denote a joint connected to μ through one of the network rods. In this way, ν ∈ [scr N, script letter N]_μ, where [scr N, script letter N]_μ is the set of all joints connected to joint μ through a single rod. The rods themselves and their various fields will be labelled with a two component index, μν, comprised of the two joints the rod connects. Specifically, u_μν(z,t) is the displacement field of rod μν with the z = 0 end of the rod being at joint μ. Similar notation also applies to other fields as well as rod specific parameters such as the length, L_μν, and cross sectional area, A_μν. Exchanging the index order thus also reverses the directionality of the rod, leading to the relations

u_μν(z,t) = −u_νμ(L_μν − z,t), (7a)

v_μν(z,t) = −v_νμ(L_μν − z,t), (7b)

σ_μν(z,t) = σ_νμ(L_μν − z,t). (7c)

We can also apply the Fourier expansion used in eqn (5) alongside this notation to obtain the index exchange laws for F_μν(ω) and B_μν(ω). Letting τ_μν = L_μν/c_μν allows for these to be expressed as

F_μν(ω) = −B_νμ(ω)e^iωτ_μν, B_μν(ω) = −F_νμ(ω)e^{−iωτ_μν}. (8)

Finally, we can Fourier transform eqn (6) in time and combine the result with eqn (8) to produce the relations

ṽ_μν(0,ω) = iω(F_μν(ω) + B_μν(ω)), (9a)

[small sigma, Greek, tilde]_μν(0,ω) = iωΓ_μν(−F_μν(ω) + B_μν(ω)), (9b)

ṽ_μν(0,ω) = −ṽ_νμ(0,ω)cos(ωτ_μν) − iΓ_μν⁻¹[small sigma, Greek, tilde]_νμ(0,ω)sin(ωτ_μν), (10a)

[small sigma, Greek, tilde]_μν(0,ω) = [small sigma, Greek, tilde]_νμ(0,ω)cos(ωτ_μν) + iΓ_μνṽ_νμ(0,ω)sin(ωτ_μν). (10b)

From here we assign a D_μ-dimensional coordinate system, denoted as [Doublestruck R]^D_μ, to the μth joint such that the origin is located at the rest location of the joint and the set of rods connected to that joint span [Doublestruck R]^D_μ. For example, if two rods are connected at 180° (antiparallel), the dimension of the joint is 1. On the other hand, if they are not parallel or antiparallel, say the rods are connected at an angle 60°, then the dimension of the joint is 2. We can then define a unit vector, ê_μν ∈ [Doublestruck R]^D_μ, to rod μν with equivalent directionality, thus implying that ê_μν points from joint μ to joint ν. Of note is that since the opposing vector ê_νμ exists in [Doublestruck R]^D_ν, there is no implicit index exchange relation between ê_μν and ê_νμ without first defining the relation between [Doublestruck R]^D_μ and [Doublestruck R]^D_ν. However, if ê_μν and ê_νμ are expressed in the global coordinate system, denoted [Doublestruck R]^D_g, then they must of course point in opposing directions. This is particularly easy to achieve if dim([Doublestruck R]^D_μ) = dim([Doublestruck R]^D_g) for all μ. In this case, we can define an “aligned coordinate set” in which x⃑_μ = x⃑_g − r⃑_μ; where x⃑_μ is a position vector in [Doublestruck R]^D_μ, x⃑_g is the same position in [Doublestruck R]^D_g, and r⃑_μ is the position of the μth joint in [Doublestruck R]^D_g.

With the coordinate systems defined, we next investigate the dynamics of the joints by assuming that each joint is massless and incapable of carrying force. Newton's second law applied to joint μ then takes the form

where P⃑_μ(t) is the force being applied to the joint by some entity external to the system, again expressed within [Doublestruck R]^D_μ. Finally, the joint itself moves within this coordinate system with a velocity given by the vector w⃑_μ(t). Enforcing that the z = 0 end of rod μν must have a velocity equivalent to the projection of w⃑_μ in the direction of ê_μν yields the condition

w⃑_μ·ê_μν = v_μν(0,t) ∀ ν ∈ [scr N, script letter N]_μ. (12)

Eqn (11) and (12) provide the connectivity laws of the network and define how the stress and velocity fields of different rods interact.^13,14 Their form can be somewhat simplified by introducing the matrix , defined to be a D_μ × |[scr N, script letter N]_μ| matrix in which each column is a distinct ê_μν. The joint stress and velocity vectors, [small sigma, Greek, vector]_μ(t) and v⃑_μ(t), can then be defined as column vectors of length |[scr N, script letter N]_μ| such that their respective jth components give σ_μν(0,t) and v_μν(0,t) for the same ν as was used to generate the jth column of . Similarly, we construct as a diagonal matrix of size |[scr N, script letter N]_μ| × |[scr N, script letter N]_μ| whose jth diagonal entry is A_μν. Finally, treating P⃑_μ(t) and w⃑_μ as column vectors allows eqn (11) and (12) to be expressed as

Of note is that due to our construction of [Doublestruck R]^D_μ we have D_μ ≤ |[scr N, script letter N]_μ|. In the specific case of equality, must be an invertible square matrix due to the set of vectors {ê_μν|ν ∈ [scr N, script letter N]_μ} spanning the local coordinate system. This causes eqn (13a) and (13b) to become a bijective linear relation between the rod dynamics, [small sigma, Greek, vector]_μ and v⃑_μ, and the joint dynamics, P⃑_μ and w⃑_μ. Thus, the rods are all effectively uncoupled from each other when D_μ = |[scr N, script letter N]_μ|. This is because for any given stress and velocity field within a single rod there must exist an applied force and joint velocity such that eqn (13a) and (13b) are satisfied without inducing any additional dynamics in any other rod.

2.3 The Laplacian block matrix

We now seek to construct a global relation between the motions of every joint in the network. We achieve this by first noting the similarities between the element velocity and stress fields with the flow velocity and pressure of a fluid in a compliant vessel, then following a methodology previously developed for calculating the pressure distribution in such a fluid flow network.¹⁸ To begin, we combine eqn (7b) and (9a) to produce

ṽ_νμ(0,ω) = −ṽ_μν(L_μν,ω) = −iω(F_μν(ω)e^{−iωτ_μν} + B_μν(ω)e^iωτ_μν) (14a)

Solving F_μν and B_μν yields

where we have replaced ṽ_μν and ṽ_νμ with the Fourier transformed velocities of joints μ and ν respectively as per eqn (12). From here we use eqn (9b) to express the Fourier transformed stress as

Inserting eqn (17) into the Fourier transform of eqn (11) and introducing the parameter Λ_μν = A_μνΓ_μν then yields

Based on eqn (18) we can construct the block vectors W̃ and P̃ as well as the block matrix . These are objects whose individual components are themselves vectors and matrices. Specifically, W̃ and P̃ have J components each, where J is the number of joints in the network, with the μth components being the vectors and respectively, each expressed in terms of [Doublestruck R]^D_μ. Similarly, , denoted here as the network Laplacian, has a J × J structure with components defined as

thus making the μν component of a D_μ × D_ν matrix that transforms a vector in [Doublestruck R]^D_ν into one in [Doublestruck R]^D_μ. Given these definitions, eqn (18) clearly dictates

Eqn (20) provides a direct relation between the force applied to the system and its dynamics. The first equality generates a linear transformation between the joint velocities and forces. The second equality generates a similar relation in terms of Ũ, the block vector of Fourier transformed joint displacements, and is obtained from the first by noting that W̃ = iωŨ since velocity is the time derivative of displacement. Of note is that the network Laplacian derived here is similar in function to the global scattering matrix of the system,¹⁴ but importantly represents a relation over the joints rather than elements. Thus, the network Laplacian is typically smaller in size and more computationally manageable.

We nondimensionalize the eqn (20) by using a characteristic time scale, and a characteristic length scale, u₀ = p₀τ₀/Λ₀. and p₀ is a characteristic force chosen as 1 N throughout this work. In Fig. 2, we use this method to compute the dynamic response of a simple truss network. In this example (Fig. 2), we consider a periodic delta-like pulse signal and decompose it into its Fourier modes. The response of the network is computed in the frequency domain for each mode, and the full time-domain displacement field is reconstructed as a weighted superposition of these responses. Additional computational details for this example are provided in the SI (Section S5), together with further examples for different boundary conditions and network sizes (Fig. S6–S9).

Fig. 2 (A) Time-dependent displacement, u_s applied as input at the source (lower left joint) (B) response of a square network with a crossbar to the applied input, shown at four different times. Stress values are color-coded on each rod: positive values correspond to compression and negative values to extension. Input displacement is applied at the lower left joint. Fixed boundary conditions are applied at the top-left joint (only the y coordinate is fixed) and the top-right joint (both x and y coordinates fixed). Parameters: Λ = 1 for all rods; τ = 1 for rods in the square frame, for the crossbar.

Note that our Laplacian matrix has structural similarity to the dynamical matrices used in previous studies on spring networks.^24,25 In the limit ω → 0, our method reduces to a form equivalent to their dynamic matrices. While those works focused on static properties and demonstrated how connectivity and geometry give rise to phenomena such as topologically protected boundary modes, the present method can be used to explore the dynamic behavior of continuum rod networks. Additionally, it is interesting to note that graph Laplacian-based formulations have been used to study stress balance in granular systems²⁶ which again used differences between field variables much like spring networks.

3 Results

3.1 Stiffness and mass matrix method is just a second order approximation of our truss network method

Our theory, in particular eqn (20) can be interpreted as a dynamic extension of the static stiffness matrix method. In this method, is the static stiffness matrix, Ũ is the displacement vector whose elements give the static displacement from equilibrium of each joint, and Q̃ is the vector of externally applied forces that maintain this out of equilibrium positioning. The stiffness matrix can be defined by treating each rod as a Hookean spring of stiffness k_μν = A_μνE_μν/L_μν = Λ_μν/τ_μν so that the force it exerts at joint μ is given by −k_μνê_μν(ê^T_μνu⃑_μ − ê^T_νμu⃑_ν), where u⃑_μ ∈ [Doublestruck R]^D_μ is the displacement of the μth joint in its own local coordinate system. Summing this effect over all rods and using our construction of given by eqn (19) automatically provides the relation

From here, we can use the fact that Ũ and Q̃ represent static quantities to express their Fourier transforms as simply the vectors themselves multiplied by a δ-function. This allows the second form of eqn (20) to be expressed as

Equating the prefactors of the δ-functions on either side of eqn (22) yields the condition , which is precisely a statement of static equilibrium written in matrix form.

We can further explore the limiting behavior of eqn (20) by defining the consistent mass matrix of the system as

where N(x): [0,1] → [0,1] is the shape function of the rod. Note the sign negation in the off-diagonal terms which takes into account the opposing directionalities of ê_μν and ê_νμ when expressed in terms of [Doublestruck R]^D_g. Using the standard linear shape function (N(x) = 1 − x) allows the integrals to be easily performed while also yielding the relationwhere has been approximated by its Taylor expansion to second order. This expansion shows that the practice¹² of finding the natural frequencies of the system by observing where is merely a second order approximation of defining the natural frequencies by where .

3.2 Truss network method gives accurate results at a lower computational cost

To show that our approach produces results that are equivalent to previously established methods, we first consider a two dimensional structure comprised of four joints arranged in a square, a total of 5 rods (square frame with a crossbar) as shown in Fig. 3(A). We choose a global coordinate and label system such that the joint 1 is located at (0,0), joint 2 at (L,0), joint 3 at (0,L), and joint 4 at (L,L). Rods that form the square frame can be denoted by the connectivity labels 12, 13, 24, and 34. The fifth rod, labeled 23, creates a cross bar of length .

Fig. 3 (A) Schematic of a square truss structure with a crossbar. (B) Comparison of natural frequencies computed using the truss network method and standard finite element methods for the square structure with a crossbar. The natural frequencies are plotted as a function of the number of subdivisions of the rods in the finite element methods. No subdivisions are necessary with the network Laplacian method to achieve the theoretical limit. (C) Computational time required for different methods. (D) Amplitude of oscillation at the target joint (T) , in response to a periodic force applied at the source joint (S) along the (1,1) direction. Parameters: τ = 1 for rods in the square frame, for the crossbar, and Λ = 1 for all rods.

How this system responds to external forces can be determined from the use of scattering wavefronts (see Section 3.3). Here we focus on the natural frequencies given by our method and compare with finite element methods. We compare with two different finite element approaches. In finite element methods, one would define a mass matrix, , then determine the natural frequencies by observing where . Eqn (23) gives one possible definition of , denoted as the consistent mass matrix, where we divide mass between rods in a way that is consistent with the network Laplacian approach. Another even simpler definition is known as the lumped mass matrix in which N(x) = Θ(1/2 − x), with Θ(x) being the Heaviside step function, thus causing to become diagonal with . Here m_μ is the total mass of all rods connected to joint μ. This lumped mass matrix simply compresses the entire system into the joints by splitting the mass of each element on the rod evenly between its connected joints. (Note that if the rods in the original network are not subdivided, the lumped mass matrix method is the same as the simple balls and springs method.) The size of these matrices are determined by the number of subdivisions in the rods.

Either of these mass matrices can be used to calculate a set of natural frequencies, though these will generally not be equivalent due to the distinct way each handles the mass of the various system elements. The theory presented here is also distinct from this approach as the network Laplacian was developed without any approximation of the distribution of mass within the rods. We thus expect the natural frequencies calculated from the condition to represent a continuum limit of those derived from either mass matrix. To test this, we considered the effects of dividing each bar of the square structure shown in Fig. 3(B) into a number of elements. Specifically, for n divisions we introduce n − 1 evenly spaced new joints in each rod to transform it into n smaller identical rods. We then calculate the five lowest nonzero natural frequencies using all four methods (the lumped mass matrix, consistent mass matrix, network Laplacian, and reverberation matrix by steadily increasing ω and observing where the determinant of each matrix vanishes. The spectrum of the network Laplacian and reverberation matrix are completely equivalent, but the reverberation method ultimately requires more computational time.

Fig. 3(B) shows the results of this process. The network Laplacian (and, equivalently, the reverberation matrix) has a unique representation for a given network and does not involve rod subdivisions. Hence the results from these methods are represented as horizontal lines. Note that the natural frequencies of both lumped mass matrix and consistent mass matrix methods steadily converge to those of the network Laplacian as the structure becomes more finely divided. The reason that the mass matrix-based methods converge to the frequencies computed by our truss network Laplacian can be understood as follows: the truss network Laplacian method retains the full continuum dynamics (uniaxial) of each rod segment. In contrast, the lumped and consistent mass matrix methods approximate the same uniaxial continuum dynamics by discretizing each rod into smaller segments, assigning mass to the nodes. As the number of subdivisions increases, these discrete models more accurately capture distributed mass along rods, and therefore it is consistent that their computed spectra approach the values obtained from our method. In other words, our method can be viewed as the exact continuum limit of the mass matrix methods.

In the particular case considered here in which all rods are made of identical material and have identical cross sections, the network Laplacian and reverberation matrix also find a natural mode at ω = πc/L that the mass matrices do not. This is due to this particular value of ω being a resonant frequency of each rod except for the cross bar, thus requiring the more careful handling explored in the SI Section S1.

We also compare the overall computation times required to obtain the spectra plotted in Fig. 3(C) as performed on an AMD Ryzen 5 1500× processor with Numpy's linear algebra package used to calculate the determinants. Fig. 3(C) shows how these computation times increase for the mass matrix methods as the system becomes more finely divided. As mentioned previously, the natural frequencies and modes determined by the network Laplacian and reverberation matrix are invariant to such divisions. We see from this data that the network Laplacian is indeed more computationally efficient than the reverberation matrix method that involves larger matrices and that both are substantially more efficient than either mass matrix method when the system rods are subdivided into more than 3 divisions.

In many earlier studies on such network structures, a simple balls and springs method was used.^27–30 This method assumes that masses of rods are concentrated at the joints and the edges are replaced by a massless hookean spring. Note that this is same as the lumped mass matrix method with no sub-division of the edges. We computed the responses of the network with the truss method and balls and springs for a simple square with a crossbar network and also larger disordered networks that are of interest to soft matter research.

For balls and springs we assume the mass of the rods is equally shared between the two joints connecting them and is concentrated at the joints. The stiffness of the rod connecting joints μ and ν is given by k_μν = Λ_μν/τ_μν while its mass is given by m_μν = Λ_μντ_μν. This for balls and springs matrix representation equivalent to a given truss network is:

We first examined how our truss network-based method compares with balls and springs by examining computations for square with a crossbar (Fig. 3D). We applied an input force vector (1,1) to the source joint S and measured the response at the diagonally opposite joint, T. The response is defined as the amplitude of oscillation at the target joint . At low frequencies, the truss and spring methods exhibited similar behaviors, consistent with our expectation that balls and springs are a good approximation for static and low-frequency cases. The peaks in the response correspond to resonant frequencies. Note that the resonant frequencies differ between the two methods (which is in agreement with Fig. 3(B)). More importantly, beyond a frequency of approximately 2, the spring network's response dropped rapidly, while the truss network continued to oscillate. This arises because a balls and spring network is a discrete representation and has only a finite number of normal modes. A spring network cannot respond accurately to excitations beyond the highest natural frequency (the Debye frequency). Truss networks treat rods as continuous objects, allowing each rod to support an infinite number of normal modes, unlike a single spring that supports only one mode.

So far we illustrated our method using small regular networks. These are much simpler to compute, even allow for analytical calculations. Disordered networks are more relevant for metamaterial synthesis due to their high tunability and capacity to incorporate multiple functionalities.^29,31 To evaluate our method in such contexts, we compared it with the lumped mass matrix (balls and springs) approach on a disordered network derived from jammed packing of spheres (Fig. 4). Disordered network with 400 joints obtained as contact network from 2D disk packing³² (see inset of Fig. 4). We compute the network response and report it as ũ^T_x, where ũ is the amplitude of oscillation at the target joint, for a force (1,1) at the source joint. We found that beyond a certain frequency, the spring method fails to generate any response at the target joint. In contrast, the truss network method continues to yield consistent results (consistent with results discussed on a small network previously). Thus truss network method can also be used to compute the response of networks for pulse excitation or impact loads, which involves high frequency components.

Fig. 4 Comparison between the lumped mass matrix method (balls and springs) and the truss network method in a disordered network (inset). A periodic excitation is applied at the source joint (red) in the direction indicated by the red arrow. The amplitude of oscillation at the target (ũ^T) is computed at the target joint (magenta). Parameters: Λ = 1 for all rods. τ values of rods vary from a minimum of 0.015 to 0.025 with a mean value of 0.02.

3.3 Wave dispersion and analysis of network responses to impact loads

One direct result that can be obtained from the theory presented here is the well-known phenomenon of wave dispersion at a network joint. By constructing the forward and backward wave amplitude vectors at joint μ as F⃑_μ(ω) and B⃑_μ(ω) such that their respective j^th components are F_μν(ω) and B_μν(ω) with ν being the same index used to construct the j^th column of , exactly the same indexing scheme as was used for [small sigma, Greek, vector]_μ(t) and v⃑_μ(t), we can utilize the Fourier transform of eqn (13) to express F⃑_μ(ω), the set of outgoing wave amplitudes, as a function of B⃑_μ(ω), the set of incoming wave amplitudes. This can be compactly expressed by first defining the matrix in exactly the same manner as but with the diagonal terms being the respective Λ_μν. The Fourier transform of eqn (13) can then be written aswhere we have used the Fourier transform of eqn (13b) in the final equality. Eqn (28) can then be solved for , which can be used to producewhereis the transmission matrix and is the |[scr N, script letter N]_μ| × |[scr N, script letter N]_μ| identity matrix.

Eqn (31) represents the method by which one may directly calculate the amplitude of all outgoing waves from the amplitudes of all incoming waves and the force applied to the joint. It is also functionally equivalent to other representations of the same wave scattering phenomenon,¹³ but is here presented in a compacted notation. Of note is that when D_μ = |[scr N, script letter N]_μ|, is invertible and simply reduces to . In this case, the wave amplitudes of each element are completely decoupled from each other so that no energy can pass through the joint except for that which is supplied by any external forcing.

The transmission matrix can be used in a variety of ways. The scattering of wavefronts can be easily tracked by propagating them through the rods of a structure, using eqn (31) to calculate how they transmit through and reflect off of the joints, and repeating for the desired amount of simulation time. The natural frequencies of the system can even be determined by also enforcing a matching condition throughout (see SI Section S2). This method of wavefront tracking is functionally equivalent to the reverberation matrix method explored in ref. 14; wherein a global scattering matrix is constructed as a diagonal block matrix whose component matrices are the transmission matrices. The reverberation matrix is then constructed by matching the backward waves with their corresponding index exchanged forward waves and solving eqn (31) for the applied forces.

How a square structure with a crossbar responds to external forces can be determined from the use of scattering wavefronts. Fig. 5 shows a time lapse of a compression wave generated by a sudden impact moving through the structure. Eqn (31) dictates how the initial wave interacts with joint 2 to become two separate transmitted waves and one reflected wave. This process repeats itself every time one of the waves reaches joint 2 or 3, steadily creating a continuously increasing number of waves. Fig. 5 also highlights how waves merely reflect off joints 1 and 4 rather than scattering through them due to the fact that D_μ = |[scr N, script letter N]_μ| at these joints. As discussed before, this causes the rods attached to them to be effectively uncoupled and incapable of passing energy to one another.

Fig. 5 Movement of stress wave through a square structure. In the case presented here, each element is assumed to be made of the same material and have the same cross section so that all values of Λ and c are equivalent. (A) The structure is impacted by an impulse at time t = 0 in such a way as to create stress in rod 12. (B) At time t = L/3c the compressive stress wave has moved a third of the way across rod 12. (C) At time t = 4L/3c the wave has scattered off joint 2 and created three new wavefronts; a compressive wave in rod 24, a tensile wave in rod 23, and a reflected tensile wave in rod 12. (D) At time t = 7L/3c the waves in rods 12 and 24 have reflected off joints 1 and 4 respectively while the wave in rod 23 has not yet reached joint 3. Transmission matrices are explicitly computed in the SI Section S2.1.

3.4 Using our method to search for structures with enhanced vibrations at target

The ability to control the oscillation of specific joints has broad applications, ranging from earthquake or blast protection systems to piezoelectric devices. Previous reports have discussed impedance-matching or mismatching protocols to achieve the desired control.^33–39 While impedance matching in one-dimensional systems has been extensively studied,⁴⁰ two-dimensional and disordered networks present significantly more complex challenges,^41,42 requiring computational approaches.

Because our network Laplacian method is accurate and fast, it can be used to predict designs with novel mechanical prooperties. Here we sweep through design space of parameters and show structures with widely different oscillations of a target joint (T) for the same boundary conditions. We consider an example in which the target joint (T) oscillates in response to input oscillations at a source joint (S) within a square network with a crossbar, subject to boundary conditions that prevent rigid body translation and rotation. (Fixed boundary condition, ũ_x = ũ_y = 0 for lower right joint, and roller boundary condition ũ_y = 0 for upper left joint) (inset in Fig. 6A).

Fig. 6 (A) Tuning the response of a square network with a crossbar. Sinusoidal displacements of unit magnitude are given as input along x-axis at the source joint (S) and the magnitude of amplitude of oscillation at the target joint (T) along the y-axis, |ũ^T_y| is computed. The other two diagonally opposite joints are given appropriate boundary conditions to prevent rigid body translation and rotation. The diagonal rod has an impedance Λ₂, rest of the rods have an impedance Λ₁. τ₁ = τ₂ = 1. (B) Response of the target joint as a function of impedance keeping the total mass of the structure constant, by setting τ₂ = 1/Λ₂, τ₁ = 1. (C) The amplitude of oscillation of the target joint (T) for different frequencies from the truss method (solid lines) compared with balls and spring (dashed lines). τ₁ = τ₂ = 1. (D) Responses of random networks with slanted rod connections between two parallel rows of rods. Light grey lines represent responses from different network realizations. Source (S) is at (0,0), and target (T) is at (8,1.0). All slanted rods have impedance Λ₂, while horizontal rods have impedance Λ₁. Grey lines represent responses for different random realizations of the network. Blue and chocolate colored lines indicate two selected realizations with markedly different responses; insets near the curves show corresponding network structures (shown in rectangular aspect ratio for legibility). The mass of the network was not held constant as Λ₂ was varied. Time constants and frequency are set as τ₁ = τ₂ = 1 and ω = 1, respectively.

We first analyzed the response of a square frame with a crossbar as a function of the impedance of crossbar in (Fig. 6A). The four other rods have impedance, Λ₁. When Λ₂ is much smaller than Λ₁, the target joint's response is minimal, analogous to a thin rod ineffectively transmitting forces from the source joint. As Λ₂ increases, improved momentum transfer leads to an increased target response, reaching a maximum at a specific Λ₂ close to Λ₂ = Λ₁. Note that this optimal impedance value is not simply equal to Λ₁ but is a complex function of the network architecture and the frequency of the force signal at the joint. (In contrast, the equivalent one-dimensional system has the oscillation amplitude maximized when Λ₂ = Λ₁, when τ₁ = τ₂ = 1 for ω = 100. See SI, Section S4.2, two rods end to end.) We have derived an analytical relationship between response and parameters Λ and oscillation frequency ω for a single rod, confirming the accuracy of our computational model (see SI, Section S4.1). In addition to maximized response when impedance is matched, there can also be some resonant modes, which lead to sharper peaks in this response.

Further investigation involved modifying impedances Λ₂ while maintaining constant total mass by adjusting τ₂ accordingly as τ₂ = 1/Λ₂ (Fig. 6B). This approach revealed multiple amplitude maxima due to resonances, occurring when the cot(ωτ) term in the Laplacian matrix, diverges.

We then tested how this method compares with balls and springs methods (lumped mass matrix without any rod subdivisions) (Fig. 6C). This analysis also confirms that the response of the network Laplacian method and the ball-and-spring method is comparable at lower frequencies but the ball-and-spring method is not adequate for higher frequencies as discussed in previous sections. It is interesting to note, at higher frequencies, the ball-and-spring method gives similar qualitative behavior as the network Laplacian method when Λ₂ is small relative to Λ₁, however at higher values of Λ₂, the ball-and-spring method gives a constant response while the network Laplacian method shows a decreased response.

Next, we apply the method to disordered networks and give an example where our method shows that simple changes in network connectivity can lead to large changes in response (Fig. 6D). We consider a disordered structure comprising 18 joints, equidistantly arranged in two parallel horizontal lines at y = 0 and y = 1. The network is constructed by randomly pairing joints at y = 0 and y = 1 and connecting them with slanted rods until the average coordination number reaches z = 4, ensuring network rigidity. All horizontal rods have an impedance of Λ₁ = 1, while all slanted (and vertical) rods have a variable impedance Λ₂. The boundary conditions are defined as follows: the source joint (lower leftmost joint (S) at (0,0)) is subjected to a constant unit displacement amplitude, ũ¹_x = 1; the lower rightmost joint, at (8,0) is fixed along x and y axes, ũ_x = ũ_y = 0; and the upper leftmost joint, at (0,1) is fixed along the y-axis, ũ_y = 0. These conditions prevent rigid body translations and rotations. We measure the oscillation amplitude (response) of the target joint (T) (upper rightmost, at (8,1.0)) as a function of Λ₂, with fixed parameters Λ₁ = 1, ω = 1, and τ₁ = τ₂ = 1. Fig. 6(D) shows that the target joint's response is highly sensitive to the specific network realization. The grey lines represent different random network configurations, all maintaining the same average coordination number but with unique pairings of joints at y = 0 and y = 1. The response curves exhibit three distinct regimes: In the low Λ₂ regime (Λ₂ ≪ Λ₁), the response is relatively constant as Λ₂ increases. This contrasts with the low and increasing response observed in simpler geometries (Fig. 6A and C), attributable to the fact that in this system, even vertical rods have impedance of Λ₂, which is not the case for square with a crossbar system considered in Fig. 6. Even a square network with a crossbar with a more number of Λ₂ elements shows similar behavior in this regime (see SI, Fig. S10) The intermediate Λ₂ regime (Λ₂ ≈ Λ₁) shows many resonance peaks. Maximum variability between network realizations is observed here, highlighting the system's sensitivity to specific network realization. In the high Λ₂ regime (Λ₂ ≫ Λ₁), the response reaches a plateau (but distinct for each network realization). We identified two network realizations with markedly different responses in the high Λ₂ regime (Λ₂ ≫ Λ₁). The corresponding network structure is included as insets near these curves. This shows that subtle changes in network structures can lead to several orders of magnitude differences in the responses. However, it is to be noted that whole range of Λ₂ values considered in this study may not be accessible by changing areas of cross-section of rods. Larger areas could lead to steric hindrance with neighbours. Nonetheless, the strategy of changing impedances with changes in area of cross section offers a valuable and practical design approach for 3D-printed microarchitectured materials to tune their response.

4 Conclusions

In this work, we developed a computationally tractable approach for computing the dynamic responses of complex truss networks. We derived a network Laplacian formalism that relates displacements (or velocities) of joints to applied forces, taking into account the mass distribution along the rods. We demonstrated that our method accurately computes natural frequencies of truss networks through comparison with more detailed finite element methods on simple structures. This method is at least an order of magnitude faster than mass matrix techniques (finite element) for a reasonable level of accuracy. We then compared with balls and springs model and showed that while balls and springs model gives comparable results at low frequencies, it breaks down at high frequencies (beyond the highest normal mode or the Debye frequency of spring networks).

Our method can contribute to research efforts in uncovering the effects of network topologies, impedance distributions, and other network properties on the mechanical responses of disordered soft matter systems. For example, our specific observation of high variability between random realizations in Fig. 6(D) illustrates the potential for fine-tuning mechanical properties through structural design using our method. Moreover, the contrast with simpler geometries emphasizes the importance of considering larger networks in predicting and designing mechanical responses. Furthermore, the observed resonance behavior in the intermediate Λ₂ regime hints at the possibility of designing switchable materials, where small changes in the network properties could lead to significant shifts in mechanical response. This could have applications in areas such as vibration damping, acoustic metamaterials, or responsive soft robotic systems. With suitable optimization techniques, this approach could aid in designing materials with tailored mechanical properties through the large design space provided by disordered networks.

We believe this approach could be a valuable addition to computational toolkits for designing micro-architectured 3D printable materials. This method, coupled with suitable optimization methods, can help in efficiently designing networks with nontrivial mechanical properties. One can also apply the method to understand more complex disordered and hierarchical biological structures, such as bones. Finite element methods are not appropriate for hierarchical cases as they require refinement of mesh sizes when modeling different levels of arrangement. The truss network Laplacian method is specifically useful for hierarchical truss structures as it relies only on joint positions and connectivities and doesn't involve the discretization of rods.

While our method is efficient and runs fast for a wide range of cases it does suffer from several drawbacks that in some contexts limit its applicability. Currently, our approach is appropriate for designing only small-amplitude dynamics of stretch-dominated structures. We ignore transverse forces, i.e., rods cannot transmit forces in directions perpendicular to their long axis. We also ignore the bending of the rods. This can be crucial if the structure is floppy or if the applied forces are large. Thus, this method, while universally applicable to any network of rods with pin joints, becomes particularly powerful and accurate only for 3D printed materials when coordination numbers are at or above the critical value dictated by Maxwell's criterion for rigidity. Finally, here we considered perfectly elastic rods that don't dissipate any mechanical energy. This is not true for real materials, especially polymeric materials used for 3D printing which dissipate mechanical energy over time.

Another consideration in assessing the applicability of our model is the concept of a representative volume element (RVE), defined as the smallest material volume whose effective properties match those of the bulk.^43,44 In metamaterials, when the excitation wavelength becomes comparable to or smaller than the RVE size, the internal microstructure of the RVE begins to influence the response, and additional variables may be required to capture its behavior accurately. For a rod network model, the RVE is effectively the full length of a rod, which limits its validity in the regime where the excitation wavelength is smaller. In contrast, by treating each rod as a one-dimensional continuum, our formulation has an RVE that is, in principle, infinitesimally small, allowing it to resolve wavelengths much smaller than the rod length. However, in real systems, making wavelengths smaller and smaller will eventually excite complex physical and nonlinear effects. Once such effects become significant, the fundamental assumptions of our model break down, and the network Laplacian we develop ceases to be an accurate representation. Nonetheless, within the range where these additional mechanisms (specifically non-linear and viscoelastic effects) are negligible, our method remains more accurate than traditional spring network approaches for short-wavelength regime.

Most of the examples presented in this paper assumed a periodic excitation with a single frequency. However, if the system were to be subjected to more complex forcing (involving multiple discrete frequencies, or a conitnuum of frequencies), then we would have to add or integrate over frequency space (as in the example Fig. 2), which would introduce computational convergence issues. Additionally, for a non-recurrent (non-periodic) impact, we would have to consider a cutoff frequency as we cannot integrate all the way out to infinity. Finding an appropriate way to numerically integrate the inverse of a Laplacian matrix would be interesting to explore.

Future directions also include further improving the model by incorporating rod bending, transverse forces, and viscoelasticity. These improvements would make the model more suitable for designing networks with complex mechanical properties and more accurate at extremely small excitation wavelengths. In particular, this approach could be used to design metamaterials with controlled mechanical energy or vibrations – an area of interest for soft robotics. Additionally, the method can help study the stress response of complex hierarchical structures, such as bone tissue, and how these tissues remodel under stress. This may offer insights into adaptation principles under external periodic forces and aligns with current efforts in soft matter research to explore physical learning in biological systems. Beyond understanding biological structures, this approach can also be applied to the development of artificial adaptive mechanical metamaterials that dynamically modify their properties in response to external conditions.

Author contributions

SF and NS contributed equally to the paper. SF contributed to model development, running numerical computations, analysing data, and writing the paper. NS worked on developing the code, running numerical computations, analysing data and writing paper. PKP and EK conceived the study, developed the formalism, and wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study’s findings are included in the article, with additional datasets and analyses provided in the Supplementary information. See DOI: https://doi.org/10.1039/d5sm00619h.

The code used for this study is freely available at https://github.com/nsarpangala/truss_network_spectral_node.

Acknowledgements

This research was funded by the Army Research Office (ARO) through the Multidisciplinary University Initiative (MURI) Grant No. W911NF2210219, the University of Pennsylvania Materials Research Science and Engineering Center (MRSEC) through Grant No. DMR-1720530 and DMR-2309043 and the Simons Foundation through Grant No. 568888. PKP acknowledges support for this work through a seed grant from Penn's Materials Science and Engineering Center (MRSEC) grant DMR-1720530.

Notes and references

1. M. Askari, D. A. Hutchins, P. J. Thomas, L. Astolfi, R. L. Watson, M. Abdi, M. Ricci, S. Laureti, L. Nie and S. Freear, et al. , Addit. Manuf., 2020, 36, 101562.
2. A. M. Torres, A. A. Trikanad, C. A. Aubin, F. M. Lambers, M. Luna, C. M. Rimnac, P. Zavattieri and C. J. Hernandez, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 24457–24462.
3. A. F. Van Tol, V. Schemenz, W. Wagermaier, A. Roschger, H. Razi, I. Vitienes, P. Fratzl, B. M. Willie and R. Weinkamer, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 32251–32259.
4. J.-Y. Jung, A. Pissarenko, N. A. Yaraghi, S. E. Naleway, D. Kisailus, M. A. Meyers and J. McKittrick, J. Mech. Behav. Biomed. Mater., 2018, 84, 273–280.
5. J. Mueller, K. H. Matlack, K. Shea and C. Daraio, Adv. Theory Simul., 2019, 2, 1900081.
6. R. N. Glaesener, J.-H. Bastek, F. Gonon, V. Kannan, B. Telgen, B. Spöttling, S. Steiner and D. M. Kochmann, J. Mech. Phys. Solids, 2021, 156, 104569.
7. P. Wang, F. Casadei, S. H. Kang and K. Bertoldi, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 020103.
8. A. O. Krushynska, A. Amendola, F. Bosia, C. Daraio, N. M. Pugno and F. Fraternali, New J. Phys., 2018, 20, 073051.
9. J. Zhang, B. Koo, Y. Liu, J. Zou, A. Chattopadhyay and L. Dai, Smart Mater. Struct., 2015, 24, 085022.
10. Y. Lu, X.-Y. Liu and G.-H. Hu, Macromolecules, 2022, 55, 4548–4556.
11. B. Branch, A. Ionita, B. E. Clements, D. S. Montgomery, B. J. Jensen, B. Patterson, A. Schmalzer, A. Mueller and D. M. Dattelbaum, J. Appl. Phys., 2017, 121, 135102.
12. R. D. Cook, et al., Concepts and applications of finite element analysis, John Wiley & Sons, 2007.
13. S. M. Howard and Y.-H. Pao, J. Eng. Mech., 1998, 124, 884–891.
14. Y.-H. Pao, D.-C. Keh and S. M. Howard, AIAA J., 1999, 37, 594–603.
15. D. Pölz, M. H. Gfrerer and M. Schanz, Wave Motion, 2019, 87, 37–57.
16. M. C. Messner, M. I. Barham, M. Kumar and N. R. Barton, Int. J. Solids Struct., 2015, 73, 55–66.
17. G. Trainiti, J. J. Rimoli and M. Ruzzene, Int. J. Solids Struct., 2016, 97, 431–444.
18. S. Fancher and E. Katifori, Phys. Rev. Fluids, 2022, 7, 013101.
19. J. F. Doyle, Wave propagation in structures: spectral analysis using fast discrete Fourier transforms, 1997.
20. A. Chakraborty and S. Gopalakrishnan, Int. J. Solids Struct., 2003, 40, 2421–2448.
21. U. Lee, Spectral element method in structural dynamics, John Wiley & Sons, 2009.
22. X. An, H. Fan and C. Zhang, J. Sound Vib., 2020, 475, 115292.
23. S.-L. Zuo, F.-M. Li and C. Zhang, Acta Mech., 2016, 227, 1653–1669.
24. X. Mao and T. C. Lubensky, Annu. Rev. Condens. Matter Phys., 2018, 9, 413–433.
25. C. L. Kane and T. C. Lubensky, Nat. Phys., 2014, 10, 39–45.
26. K. Ramola and B. Chakraborty, J. Stat. Phys., 2017, 169, 1–17.
27. A. Bergne, G. Baardink, E. G. Loukaides and A. Souslov, Extreme Mech. Lett., 2022, 57, 101911.
28. M. Sheinman, C. Broedersz and F. MacKintosh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021801.
29. J. W. Rocks, N. Pashine, I. Bischofberger, C. P. Goodrich, A. J. Liu and S. R. Nagel, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 2520–2525.
30. J. Huang, J. Zhang, D. Xu, S. Zhang, H. Tong and N. Xu, Curr. Opin. Solid State Mater. Sci., 2023, 27, 101053.
31. J. W. Rocks, H. Ronellenfitsch, A. J. Liu, S. R. Nagel and E. Katifori, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 2506–2511.
32. J. Andrejevic, Physical learning with elastic networks, 2023, https://github.com/jandrejevic12/physical_learning.
33. V. T. Rathod, Electronics, 2019, 8, 169.
34. T. Rahimzadeh, E. M. Arruda and M. Thouless, J. Mech. Phys. Solids, 2015, 85, 98–111.
35. S. Chen, Y. Zhang, C. Hao, S. Lin and Z. Fu, Phys. Lett. A, 2014, 378, 77–81.
36. T. Guo, X. Zheng and P. Palffy-Muhoray, Soft Matter, 2021, 17, 4191–4194.
37. X. Zhang, Z. Qu and H. Wang, iScience, 2020, 23(5), 101110.
38. M. Verma, M. Sivaselvan and J. Rajasankar, Struct. Control Health Monit., 2019, 26, e2402.
39. F. Ma, C. Wang, Y. Du, Z. Zhu and J. H. Wu, Mater. Horiz., 2022, 9, 653–662.
40. H. Kolsky, Stress waves in solids, Courier Corporation, 1963, vol. 1098.
41. C. A. Blaurock, PhD thesis, Massachusetts Institute of Technology, 1994.
42. H.-Z. Zhu, J.-H. Wu and Y.-D. Sun, Appl. Sci., 2022, 12, 8863.
43. D. M. Kochmann and K. Bertoldi, Appl. Mech. Rev., 2017, 69, 050801.
44. L. Shuguang and E. Sitnikova, Representative volume elements and unit cells: concepts, theory, applications and implementation, Woodhead Publishing, 2019.

---

Footnotes

† These authors contributed equally to the work.

‡ Present Address: Department of Physics, University of Michigan, Ann Arbor, MI, USA.

---