Discrete mass-spring structure identification in nonlocal continuum space-fractional model

This paper considers discrete mass-spring structure identification in a nonlocal continuum space-fractional model, defined as an optimization task. Dynamic (eigenvalues and eigenvectors) and static (displacement field) solutions to discrete and continuum theories are major constituents of the objective function. It is assumed that the masses in both descriptions are equal (and constant), whereas the spring stiffness distribution in a discrete system becomes a crucial unknown. The considerations include a variety of configurations of the nonlocal parameter and the order of the fractional model, which makes the study comprehensive, and for the first time provides insight into the possible properties (geometric and mechanical) of a discrete structure homogenized by a space-fractional formulation.

An important family of FM concepts includes formulations where, in classic continuum mechanics (CCM) theory, the integer differential operators acting on a spatial variable are replaced by fractional operators; such models are called space-fractional continuum mechanics (s-FCM) [19][20][21][22][23][24][25]. The physical interpretation of s-FCM is that FD homogenizes, in a phenomenological sense, the underlying microstructure [26]; thus, the scale effect is included within FC, and the s-FCM model is nonlocal. Nevertheless, to the best of our knowledge no one has answered the following question: Is there a discrete structure that is homogenized by s-FCM?.
In this paper, we provide the answer to the above question by identifying a discrete mass-spring structure, whose behaviour (static and dynamic) is mapped with high precision by an s-FCM approach. Such a concept is a kind of inverse problem to recent trends in physics/mechanics on the study of continualization techniques, such as the Taylor series-based method or the shift operator expansion method [27][28][29]. The overall problem is defined as a 1D optimization task with several variables (altogether approximately 2.5 million analyses were conducted). Furthermore, the considerations include a variety of configurations of the non-local parameter and the order of the fractional model (altogether twenty discrete spring-mass structures were identified), which makes the study comprehensive.
The paper is structured as follows. In sect. 2, the fundamentals of s-FCM and mass-spring discrete models are presented. Section 3 deals with an inverse problem formulation, while sect. 4 presents the results and discussion. The final section provides a summary and conclusions.

Riesz-Caputo fractional derivative
The Riesz-Caputo fractional derivative is constructed as linear combination of Riemann-Liouville fractional integrals and an integer order derivative. From the point of view of mathematical modelling of mechanical phenomena, it is important that the Riesz-Caputo derivative of a constant is zero (which guarantees that the strain field does not change for rigid body motion as in the classic approach), and classic boundary conditions can be applied.
First, we define the integral where P = a, t, b, p, q is the parameter set, t ∈ [a, b] (an independent variable, identified as a space variable later in this paper), p and q are real numbers, k(t, τ ) is a kernel that may depend on α, and f is any function defined almost everywhere on (a, b) with values in R. Next, assuming that the kernel k(t, τ ) is a difference kernel that depends on α, we obtain the following integrals: -for P = a, t, b, 0, 1 , we have where Γ is the Euler gamma function. The integrals given by eq. (3) and eq. (4) are the left and right Riemann-Liouville fractional integrals, respectively, for which we introduce the notation a I α t f (t) and t I α b f (t). The linear combination of eq. (3) and eq. (4) leads to the Riesz potential, namely, Finally, the Riesz-Caputo fractional derivative is defined as where n = α + 1 and · denotes the floor function. Alternatively, eq. (6) can be expressed as where the left and right Caputo derivatives are Furthermore, the following index rules are true (α ∈ (0, 1)): where m ∈ N and N denotes the set of natural numbers.

1D fractional elasticity
To understand the general concept of 3D s-FCM, the reader should consult [20,30]. Based on the results presented in [20], the governing equations for a 3D s-FCM elastic body that occupies a volume Ω with boundary ∂Ω, together with the assumption of a variable length scale [31], have the form where σ ij denotes the components of the Cauchy stress tensor, b i denotes the components of the body force vector, ρ is the density,( ·) and(·) denote the first and second time derivatives, respectively, ♦ ε ij denotes the components of the small fractional strain tensor, f = f (x) is the length scale (herein a known function), α is the order of continua α ∈ (0, 1] (for α → 1, there exists a smooth passage from this formulation to conventional (local) continuum mechanics), x is a spatial variable, U i denotes the components of the displacement vector, L e ijkl denotes the components of the stiffness tensor, n i denotes the components of the outward unit normal vector, t i denotes the components of the Cauchy traction vector, and ∂Ω U and ∂Ω σ are the parts of the boundary ∂Ω where the displacements and tractions are applied, respectively. Because linear elasticity is considered, the stiffness tensor is where λ and μ are the Lame parameters and δ ij is the Kronecker delta operator. Equations (13) can be reduced to a 1D case in a classic manner (see [20]); as a result, the following integro-differential equation is obtained (x = x 1 ): where E denotes the Young modulus. Equation (15) allows easier interpretation of the f parameter as the size of the surroundings, which influences the material point being considered. Moreover, in the following, we assume that f near the ends of the structure is linearly reduced to zero, which means that the nonlocal action is less pronounced for the points on the boundary, i.e., in a real body the surroundings are limited for the grains at the boundary [31][32][33]. Finally, there are two types of boundary conditions for the 1D configuration: -for both ends (i.e., x = x 0 and x = x r ) of a 1D body, displacements are prescribed -for the left end (i.e., x = x 0 ), a displacement is given and for the right end (i.e., x = x r ), a force is given where x 0 and x r denote points on the boundary; see fig. 1.

Approximation for FD operators
Because there is no analytical solution to eq. (15), a numerical approximation of the FD operators is needed; we follow the expressions proposed in [34,35]. Therefore, first, we assume the spatial discretization of a 1D fractional body of length L into r equal subintervals Δx (see fig. 1) by the rule with x 0 = 0 and r ∈ N.
Next, for the discretization point x i , the left derivative is approximated by where U (n) (x ai j ) denotes a classic n-th derivative at x = x ai j ; by analogy, for the right derivative we obtain where

Static case
We assume that the inertia term in eq. (15) can be neglected, i.e., ρÜ (x, t) = 0, so U = U (x). Then, assuming that α ∈ (0, 1], n = 1 in eqs. (19) and (20), the behaviour of the i-th node is governed by (see fig. 1) and (·) , (·) denote the first and second order derivatives which are approximated using the classic central finite difference scheme. For α = 1 (n = 1), eq. (22) Finally, we can write the set of r + 1 equations obtained from eq. (22) as a global algebraic problem, namely, where K denotes the fractional stiffness matrix, q includes the unknown displacements in the discretization nodes

Eigenvalue problem
In eq. (15), we assume that b(x, t) = 0 and that the separation of the variables holds; therefore, and eq. (15) can be rewritten in the form where λ 2 f = ω 2 f ρ, and ω f = ω f (α, f ) denotes the natural frequency for the s-FCM. Next, assuming as for the static case that α ∈ (0, 1] and n = 1, and using approximations eqs. (19) and (20), the behaviour of the i-th node is governed by As for the static case, for α = 1 (n = 1), eq. (28) reduces for m = M = 1 2 to the classic central finite difference scheme, i.e., Finally, we write the set of r + 1 equations obtained from eq. (28) as a global algebraic problem, namely, with the solution where I denotes the identity matrix. In the following sections, without loss of generality, for static and dynamic cases, we consider the case when both ends are clamped. It follows from the considerations presented in [31] that the unknown displacements are U 1 up to U r−1 , whereas U −1 and U r+1 should be eliminated (see fig. 1). Based on equating the finite difference approximation of the second order derivative using central and forward schemes at point x 0 , and by analogy central and backward schemes at point x r we have Consequently, in eqs. (24) and (30) we assume that: i) for i = 1: m = 1, and M = 0, and for all U in eqs. (22) and (28), a forward scheme is applied; ii) for i ∈ {2, 3, . . . , r − 2}: m = M = 1 2 , and for all U in eqs. (22) and (28), a central scheme is applied; and iii) for i = r − 1: m = r − 1, and M = r − 2 and for all U in eqs. (22) and (28), a backward scheme is applied.

Governing equation
We assume that the discrete chain comprises of springs (k i ) and masses (m i ) as presented in fig. 2.
The behaviour of the i-th concentrated mass can be described as follows: where the springs forces (Ñ ) are directly related to their elongatioñ whereŨ i is the i-th node displacement. After inserting the spring forces from eq. (35) into eq. (34) the governing equation can be presented in the form

Static case
With the assumption thatm i d 2Ũ i dt 2 = 0 (inertia forces are negligible), the behaviour of the i-th node is governed by the equationk We can write eq. (37) for all masses, which leads to an algebraic problem whereK is a stiffness matrix for a discrete system andq = (Ũ 0 , . . . ,Ũ i , . . . ,Ũ rs ) T represents concentrated masses displacements.

Eigenvalue problem
We assumeb i = 0 and separation of variablesŨ As a consequence, the behaviour of the i-th mass is governed by the equatioñ which written for all masses allows us to formulate the following eigenproblem: In the above,λ =ω 2 whileM is a diagonal mass matrix.

Boundary conditions
The boundary conditions in a mass-spring model consisting of rs concentrated masses (see fig. 2) correspond to those defined for the fractional model (eqs. (16) and (17)), namely, -for both ends, displacements are knownŨ -for the left end, a displacement is given, and the right end, a force is giveñ 3 Identification

Introduction
The main goal is to calibrate the spring stiffnesses in a spring-mass model to obtain a mechanical response similar to that of the fractional model. Two types of analyses are considered, namely, static and natural frequency (eigenproblem) analyses. For both models, the total mass is equal and uniformly distributed, so ρ(x) = const. andm i = const. The changes in the spring stiffnesses are driven by comparison of the eigenvalues and static displacements. The mechanical responses are expected to match for all the combinations of fractional model parameters α ∈ {0.6, 0.7, 0.8, 0.9, 0.999} and f ∈ {0.025, 0.05, 0.1, 0.2}.

Computational models
A 1D body fixed at both ends is analysed ( fig. 3). For the sake of simplicity, the nondimensional case is considered, so L = 1, and for the fractional model, E = ρ = 1. In the discrete model,m i depends on discretization.

Mass-spring
The fundamental assumption made in this study is that the length scale of the fractional model ( f ) determines the module of length d in the spring-mass model. This module represents a repeatable pattern in the structure (RVE) [36]. Thus, finding the spring stiffnesses k i in the single module is enough to identify the whole structure stiffness. The module length d equals the fractional length scale f . Therefore, the number of modules is expressed by n m = L/ f and equals 5, 10, 20 and 40 for f ∈ {0.025, 0.05, 0.1, 0.2}, respectively. The topology and mechanical properties of the module constitute the set of unknowns, i.e., the locations of concentrated masses, number of springs and spring stiffnesses. The consideration of all the possible unknowns where strong interactions between the unknowns are clearly visible would make the optimization extremely difficult. Therefore, in this study, the concentrated masses are uniformly spanned, and the number of springs s ∈ {4, 5, 6}, which leads to s unknowns-spring stiffnesses: which have to be identified for each combination of α and f analysed. Due to the assumption of the same mass for both models, the single concentrated mass in the spring mass model is given bym where the numerator for the non-dimensional case equals one (A(x) = const. = 1). The concentrated mass represents the mass of a fragment of the structure of length Δx analysed (see fig. 3). The mechanical behaviour of the spring-mass structure for the static case is described by eq. (37). The mass load is applied to all the concentrated masses and equals, for the i-th mass, where the field intensity g = 1.0. To find the displacements, the system of linear equations, eq. (43), is solved. For the case of the natural frequency, the eigenproblem, eq. (41), is solved.
To determine the results of the static and dynamic problems, a procedure was developed in Python. For the massspring model, to facilitate comparison of the results, the displacement field and eigenvectors are computed for r + 1 uniformly distributed points (see eq. (18)) with linear interpolation between the concentrated masses. The number of subintervals, r + 1, agrees with the discretization applied in the fractional model that the results are compared to (see the next section).
(47) Then, the static behaviour of the structure analysed is described using s-FCM according to eq. (22) and g(x) = 1.0 as mentioned. The final displacements for the statics are obtained by solving the algebraic problem given in eq. (24), whereas the eigenvalues and eigenvectors are computed on the basis of eq. (30).
For all the models used the concept of variable length scale is applied [31], which eliminates the virtual boundary layer [37] and its influence on the results. As mentioned, near the ends of the structure, the length scale linearly reduces from max f to 2Δx. As a result, no information on the displacements outside the body is required when eq. (32) and eq. (33) are applied. The distribution of all the length scales used f (x) is presented in fig. 4.
The fractional body is discretized into different numbers of subintervals depending on the number of modules in the spring-mass model to which the results are compared, namely, Δx ∈ {0.005, 0.0025, 0.001(6)} for s ∈ {4, 5, 6}, respectively.
The global system of equations (eq. (24)) was built using a dedicated procedure developed in Python. Due to the varying length scale, the governing equations for the particular points were built dynamically (see github.com/szajek), and the system of equations was solved using the LU decomposition [38] with partial pivoting and row interchanges.

Optimization problem
The following optimization problem was formulated for each configuration of f/d and α: where s denotes the number of springs in a module (see fig. 3). The optimization sub-task is defined as a weighted sum: subject to: k lower ≤ k i ≤ k upper , where w 1 and w 2 are coefficients responsible for scaling and weighting the objective function components, and the components can be calculated as where p is the order of a p-norm, U denotes the displacements, and ω f,i is the i-th natural frequency. In turn, k is the vector of the design parameters, i.e., the stiffnesses of the springs in a module (see fig. 3). ω f,i andω i are normalized in objective function to prevent the objective prioritization of higher eigenvalues. Due to the discretization applied (see sect. 2), the objective function is calculated based on the numerical approximations, namely, where Considering the boundary conditions applied (see sect. 3.2) U i = 0 for i ∈ {0, r}, (Ū i − U i )δ i = 0 for the edge points and δ i = Δx = const. for the remaining points, and eq. (51) can be simplified to the form In the optimization procedure, p = 2 is used for both O p stat and O p dyn ; however, the results for p = ∞ are also computed to provide the largest difference between the displacements of discretization points and the largest normalized difference between the eigenvalues. The boundary values of k are k lower = 10 −6 and k upper = 60(EA/Δx). For the nondimensional case, k upper = 60/Δx. In this way, the local stiffness of the mass-spring model can range from almost zero to 60 times larger than that of the fractional model. Finally, the coefficients are w 1 = 1.0 and w 2 = 3.33 · 10 2 . The coefficients w 1 and w 2 can be implicitly divided into a "scale factor" and "weight", respectively. The "scale factor" is responsible for making O p stat and O p dyn of similar order, while the "weight" concerns the "importance" of the objectives. These parameters were established on the basis of the preliminary study. Considering the final values presented in table 1, O p stat = 1.95 · 10 −4 on average, while O p dyn = 2.67 · 10 −2 ( · denotes an average value). Thus, the "scale factor" equals approximately sf = O p dyn / O p stat = 136.9, which leads to the "weight" w 2 /sf = 2.43. Therefore, on average, the static component is twice as important as the dynamic component.
The total number of eigenvalues considered is limited to N f and depends on the length scale ( f ). The primary reason for this limitation is that the wavelength in the spring-mass model with constant spring stiffness has to fulfil the inequality [29] where k is the wave number and a denotes the distance between the masses. Equation (54) is that the complexity of the displacement field is limited by the discrete microstructure. In the case analysed, each module in the spring-mass model can be reduced (due to the parallel connection) to a single mass and spring as presented in fig. 5. Thus, the distance a = f = d . Furthermore, considering that k = 2π/λ f and assuming the minimal allowed λ f = 2L/N f , eq. (54) can be rewritten as In conclusion, eq. (55) states that the maximal number of eigenvalues equals 5, 10, 20, and 40 for the length scales ( f ) 0.2, 0.1, 0.05, and 0.025, respectively. Another reason for imposing the limit on N f is that the "homogenization" performed using the fractional model is not ideal. The assumption that the module length equals the length scale is not valid close to the edges (see fig. 4). Therefore, the results differ and the discrepancy increases for higher natural frequencies where the wavelength decreases. The effect is observed as a change in the wavelength along the body. In the mass-spring model, the wavelength is constant. In this study, the following proportion is assumed to be the limit: where D i is obtained based on eigenvector analysis, as presented in fig. 6. Finally, the numbers of natural frequencies analysed, N f , are 4, 8, 16, and 32 for length scales, f , 0.2, 0.1, 0.05, and 0.025, respectively. These limits, being more restrictive, are used to calculate the objectives. In our future work we plan to ease the restriction on constant d .

Algorithm
The problem described by eq. (48) was solved by coupling the computational model with an optimizer. The general overview of the procedure is presented in fig. 7. Two optimization strategies were used, a local search algorithm and the same algorithm hybridized with an evolutionary algorithm. As a local search optimizer, the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm [39,40] implemented in the SciPy package was chosen [41]. All the spring stiffnesses for the initial configuration were assumed to be the same as those in the fractional model (the exact values depended on the modules and segment numbers). The second approach, the hybridized algorithm, started with a set of 70 propositions (individuals) and improved the proposition in a run of evolutionary operators (10 generations). The best solution, in turn, was the starting point for a local search. The initial configurations of the design parameters were sampled from a feasible design space (see sect. 3.3). The most important settings used for the L-BFGS algorithm were assumed as follows: the maximum number of iterations was 15000, the minimal relative change of an objective was 2.22 · 10 −9 (default value in scipy package), and the step size used for numerical approximation of the Jacobian was 10 −3 . The step size was established by trial and error and led to sufficiently high Jacobian values to drive changes in the variables. The evolutionary strategy configuration, on the other hand, was as follows: tournament selection with a group of 3 individuals, two-point crossover with a probability of 0.5, and a Gaussian mutation with a mean of 0.0, standard deviation of 1.0 and probability of 0.02 for each input parameter independently.

Results
Due to the stochastic nature of the evolutionary algorithm, the optimization problem eq. (48) was solved 40 times for each configuration of parameters f = d , α and s, resulting in 2400 optimization runs and approximately 2.5 million analyses. The values of the objective functions for the best solutions are collected in table 1. The objective values were calculated using eqs. (50) for p = 2 (Euclidean norm) and p = ∞.   The spring stiffnesses for the optimal solutions are presented in fig. 8 in symbolic form. The horizontal axes represent the coordinates along the segment. For each fractional material order α (vertical axes), a bar with spring stiffness is given. For easier comparison, the spring stiffness is multiplied by Δx/A(x) = Δx (for the non-dimensional case, A(x) = const. = 1), which leads to a value that can be directly compared to Young's modulus in the fractional model (we assumed 1.0 for all the analyses). Thus, the values present the changes in the local stiffness along a module. Few important observations can be made from analysing the graph. First and foremost, the spring stiffness identified differs from that assumed in the fractional model. Second, the local stiffness in the mass-spring model is constant and similar to that of the fractional model as α approaches 1 (all the values equal 1.0). The third observation is that the spring stiffnesses vary significantly within the module from 0.15 to 36.98 and neither a clear, repeatable pattern nor any dependency on the fractional parameters is seen. The final observation concerns the number of segments in a module, s, which provides the optimal solution. In almost all the cases, s ∈ {4, 5}. The highest value, s = 6, was obtained in only two cases where α is 1. In these cases, the solutions are the easiest to find because the local search algorithm starts close to the final solution (uniform stiffness; see sect. 3.4).
The solutions obtained have similar eigenfrequency, N f , in the ranges analysed (see sect. 3.3). A graphical comparison is presented in fig. 9. The best solution with the maximal difference for natural the frequency was found for f = d = 0.1 and α = 0.9 with a relative error of 0.27%, while the worst solution was found for f = d = 0.025 and α = 0.8, where the relative error reached 3.85%. However, an average error of 1.48% indicates that from a practical point of view, the results for both models agree well. Regardless of the agreement of the eigenfrequences, the particular eigenvectors should be verified to correspond. In this work, the modal assurance criterion (MAC) [42,43] is used. The MAC is a statistical indicator that qualitatively compares eigenvectors. The MAC is computed as a normalized scalar product of two sets of vectors according to the formula where i and j indicate the number of modal vectors in the fractional and mass-spring models, respectively. The MAC checks whether the eigenvectors are linearly dependent; thus, the MAC is an approximation of the orthogonality check (mass is not considered). For each combination of model vectors, the MAC is computed and collected in matrices in fig. 10. In all the cases, the values outside the diagonal are almost zero, while the diagonal values are close to one for the first modes and decrease with the mode number, reaching 0.57 in the worst case for mode number 32 and fractional parameters d = 0.025 and α = 0.6. The graphs also clearly show that an α of 0.6 enforces this effect. The decrease observed is a consequence of the homogenization of an interior structure with the fractional model. The massspring model is capable of representing the changing stiffness within a module, providing a more detailed description. As a consequence, the eigenvectors obtained for the mass-spring model are not as smooth as those obtained for the fractional model. A comparison of selected eigenvalues for f = d = 0.025 and α = 0.6 is provided in fig. 11. Qualitatively, the eigenvectors are similar, however, the mass-spring model is capable of capturing a complex response within a module. Additionally, Therefore, regardless of the MAC values lower than 0.9, the eigenvectors are proven to correspond to each other. The agreement of both the natural frequencies and eigenvectors prove that the spring stiffnesses proposed provide similar dynamic response. The static responses of the fractional and the mass-spring models agree according to O 2 stat and O ∞ stat . While the former provides a generalized measure of similarity, the latter presents the largest differences for all the points. Additionally, by dividing O ∞ stat by Δx we can calculate the exact differences between displacements (see eq. (53)). Several important observations are made by analysing the results. First and foremost, the final objective values significantly vary depending on the module length and fractional material order. Regarding O ∞ stat , the worst results are obtained for f = d = 0.2 and α = 0.6, where = O ∞ stat /Δx = 2.36 · 10 −2 , which has a relative error of = /Ū max = 18.0%. The best solution, on the other hand, yields = 1.1 · 10 −20 (the relative error¯ is comparable to machine precision) and is obtained for α = 0.9999 regardless of the module length. Second, there is a relationship between the fractional parameters and the objectives. The objective tends to zero as α approaches 1. Similarly, for all cases as f = d approaches 0 (the dimensions of the body become considerably larger than the characteristic length scale), the error decreases. To understand the changes observed in the objectives, the displacements obtained for both models and 20 configurations of f = d and α are compared in fig. 12. The largest error results from that the massspring model, such as in the natural response analysis, is able to capture an interior structure, i.e., the stiffness changes within a module. As a consequence, the displacements are not as smooth as the homogenized results obtained with the fractional model. The increase in¯ when α approaches 0.6 is caused by the larger relative differences in the spring stiffness in a module; thus, the inhomogeneity is averaged by the fractional model. The effect is stronger when the module length increases due to the larger span between the masses (s 6), and any relatively stiffer or softer spring elongates an extremely large or negligible amount. As a consequence, the displacements of the mass-spring model deviate from that computed using the fractional model. However, the local discrepancy observed does not change the general similarity of the displacements computed using both models.

Conclusions
In this paper, the fundamental question regarding the existence of microstructure, in the sense of a discrete system, underlying the s-FCM theory (operating on RC fractional operator) is positively resolved. To the best of our knowledge, this topic has not yet been covered in the literature. Very few papers have tried to address this problem for nonlocal mechanical models different than s-FCM (i.e., the Eringen and strain gradient types).
Our comprehensive analysis (approximately 2.5 million analyses were conducted) covers the 1D configuration for both the discrete (spring-mass) and continuum (s-FCM) approaches and furthermore for the dynamic and static cases. The mechanical properties of the discrete system are treated as unknowns that are determined using a properly defined optimization task. The results obtained allow us to formulate the following statements: -independent of the fractional model parameters, i.e., order of continua α and length scale f , the equivalent periodic spring-mass system (length d ) was found, maintaining the same mass distribution in both formulations and with the fundamental restriction that f = d ; -for both the static and dynamic cases, the same microstructure was identified for a given s-FCM model; this result is crucial because for other, competitive, nonlocal formulations such a result may not exist. Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.