Identification of mechanical properties of 1D deteriorated non-local bodies

Deterioration, understood as a change of mechanical properties of devices exhibiting a strong scale effect, is a phenomenon of great importance in modern industry. This effect is simultaneously an open question in theoretical mechanics, due to the fact that none of existing non-local theories (which are necessary to describe the scale effect) seem to be universal - thus constant development of new mathematical models is necessary. This paper presents the study of the above emphasised problem in the framework of fractional continuum mechanics, as an optimization (identification) task. In the identification routine, the steering parameters are: density of material, order of material, and material length scale.


Introduction
It was in 1926 when Werner Heisenberg said "... it is the theory which first determines what can be observed ..." (Heisenberg 1989). This fundamental statement was related to atomic physics at that time, nevertheless it also holds nowadays and even for higher scales of observation, e.g. micro, meso, macro. This is because of high maturity of computer aided decision making in many branches of human activity, which causes that 'virtually' obtained results give in most cases more comprehensive insight into an analysed problem than even a very sophisticated real experiment. The same applies to the scale-effect phenomenon, which plays the central role in this paper and whose hidden aspects will be extracted hither utilising the theory of space-Fractional Continuum Mechanics (s-FCM).
Fractional Continuum Mechanics (FCM) is in a broad sense a class of mechanical models that utilise the fractional calculus (FC) (Podlubny 1999;Kilbas et al. 2006;Malinowska et al. 2015). One can say that FMC is a generalisation of classical Continuum Mechanics (CCM) because the latter is a special case of the former, or more precisely, when orders of all FC in the FCM model become integers, smooth passage to the CCM case is obtained. FCM can be classified in several ways, however the most popular classification in the literature includes the name of the variable on which the fractional derivative (FD) operates in the model. Therefore we have, e.g.: (i) timefractional models (Nan et al. 2017;Suzuki et al. 2016;Liao et al. 2017;Zhilei et al. 2016;Wu et al. 2016;Ansari et al. 2016;Faraji Oskouie and Ansari 2017;Sumelka and Voyiadjis 2017); (ii) space-fractional models (Klimek 2001;Drapaca and Sivaloganathan 2012;Sumelka et al. 2015b;Tomasz 2017;Lazopoulos and Lazopoulos 2017;Peter 2017); (iii) stress-fractional models (Sumelka 2014a;Sun and Shen 2017, a, b). It is crucial to emphasise that all FCM models are non-local, although the physical interpretation of this non-locality depends on variable on which FD acts (Sumelka and Voyiadjis 2017). Herein, as mentioned, the non-locality in space should be pointed out as being the main constituent of s-FCM.
The s-FCM theory, being space non-local, enables us to model scale-effect which is of extreme importance, taking into account constant miniaturisation in many areas of human activity i.e. the production of microstructured and nanostructured materials or micro-or nanoelectromechanical (MEMS or NEMS) devices, nanomachines, as well as in biotechnology, biomedical fields and astronomy (Drexler 1992;Martin 1996;Han et al. 1997;Fennimore et al. 2003;Bourlon et al. 2004;Saji et al. 2010).
Modelling of such bodies is a challenging task by itself; however, in this paper, we try to go further, asking (theoretically): "What are the mechanical properties of deteriorated non-local bodies?". This question is of fundamental importance considering the durability of miniaturised devices.
Identification of mechanical properties frequently requires solving an inverse problem. There are plenty of examples which confirm the usefulness of this approach when CCM is applied (Constantinescu and Tardieu 2001;Uhl 2007;Garbowski et al. 2012;Bonnet and Constantinescu 2005). However, there is little to none studies in this area when it comes to non-local models. From the few, the Vosoughi and Darabi work can be mentioned. In Vosoughi and Darabi (2017) the authors used a conjugate gradient method hybridised with a genetic algorithm to identify a volume fraction coefficient and a small scale parameter in a functionally graded nanobeam together with the Eringen non-local elasticity and a first-order shear deformation theory. The frequencies of a beam were compared for calibration. Similarly, a dynamic response was used by Kiris and Inan (2008) to estimate upper bounds of elastic modulus in a material model based on Eringen's microstretch theory. The optimization problem was then solved by a direct search method along with a micro-genetic algorithm. Another work by Diebels and Geringer (2014) dealt with calibration of the Cosserat model parameters for a foam material (an internal length and an additional stiffness parameter), based on shear stresses from a spatial finite element model. A genetic algorithm optimizer was used. None of the cited works, however, deal with a structure whose mechanical properties locally deteriorate.
A few papers devoted to the optimization problem coupled with non-local models, although not directly connected with the inverse problem, indicate that the optimization results can vary strongly when the size of a body approaches a length scale. Liu and Su (2010) performed topology optimization within a rectangular domain filled with material where the couple-stress theory was applied (Mindlin 1964). The work generalises the conventional Solid Isotropic Material with Penalisation (SIMP) model (Bendsøe 1989;Bendsøe and Sigmund 2004) and the results clearly demonstrate that the optimal solution changes when the length scale is varied. The optimal material layout approaches the solution for classical continuum mechanics when the proportion of the length scale to the minimal size of the rectangular domain tends to zero. Similarly, Veber and Rovati (2007) solved the minimum compliance problem (Bendsøe 1989) for a micropolar body (Eringen 1966) finding the optimal material distribution in a rectangular domain (plain stress) for different values of the characteristic length for bending. The final topology diversity is considerable for all analysed boundary and load schemes when the characteristic length becomes comparable with the body size. The observation is also confirmed by Sun and Zhang (2006) who optimized the topology of lightweight structures with a cellular core. This research output indicates that the solution depends on the length scale when the macrostructure has size comparable to a microstructure. The mentioned papers prove that the optimization result is sensitive to nonlocal model parameters. Thus, one should expect the same when it comes to the identification problem-especially for deteriorated bodies.
The primary goal of this paper is to capture the mechanical proprieties of 1D deteriorated non-local bodies, understood as a change of mechanical properties of devices exhibiting a strong scale effect, under the assumption that topology remains constant. The problem is formulated within s-FCM defined in Sumelka (2014b) (together with the variable length scale concept presented in Sumelka 2017), as a sensitivity analysis of mass and stiffness distribution in reference to fractional model parameters and loading scheme. The inverse approach as an optimization procedure is proposed, where without loss of the primary goal importance, test data are prepared as a displacement field for a corroded fractional body. The obtained results show sensitivity of identified corrosion depending on fractional constitutive parameters. Furthermore the results are contrasted with the solution of CCM. In this sense, the meaningful role of non-local (fractional) modelling is pointed out. In other words, the considerations show the complexity of behaviour of bodies exhibiting the scale effect which can go beyond the standard engineering intuition.
Herein, it should be pointed out that in order to achieve the above stated objectives, it was necessary to elaborate: • modified numerical scheme compared to Sumelka (2017) for better efficiency of computations -for each single optimization task it was required the execution of few thousands of iterations; • extension of numerical approximation procedure to include the inhomogeneous linear elasticity (and furthermore the nonlinear stiffness to density relationship); • methodology of deterioration including especially the stiffness to density relationship; • penalty terms to make the inverse problem successful.
The paper is structured as follows. In Section 2 s-FCM fundamentals are presented. Section 3 deals with a general idea for capturing deterioration, details on implementation of a numerical model and an algorithm as well as an inverse problem formulation. Section 4 presents results from almost 150 analyses for various load schemes and fractional body parameters, along with a discussion. Finally, Section 5 provides the conclusions.

Fractional derivative
It is commonly accepted that studies on FC has been initiated in 1695 by Leibniz and L'Hospital (Leibniz 1962). Since that times the investigation on differential equations of an arbitrary (even complex) order has become an individual branch of pure mathematics with many practical applications, as mentioned above. To understand the fundamental concepts and techniques of FC one should follow the monographs (Nishimoto 1984;Podlubny 1999;Kilbas et al. 2006;Malinowska et al. 2015;Changpin and Fanhai 2015) and papers cited therein.
Herein, one concentrates on the FD definition, which is a base for the s-FCM concept discussed in the next section, namely Caputo derivative. Thush, the so called left-sided and the right-sided Caputo derivatives are, respectively: where t denotes the independent variable (in further part of this paper identified as a space variable), is the Euler gamma function, n = α +1, · denotes the floor function, a, b are terminals, and f is any function defined almost everywhere on (a, t) or (t, b), respectively, with values in R. Next, based on the above definitions the both-sided Riesz-Caputo FD is defined as Derivative given by (3) play a central role in succeeding fractional elasticity definition.

1D inhomogeneous fractional elasticity
We consider 1D reduction of a general 3D formulation for small strains and a linear constitutive law formulated in Sumelka (2014b). It is important that this concept was validated with experimental result in Sumelka et al. (2015a), where the s-FCM concept was used to mimic the behaviour of micro-beams made of the polymer SU-8 (where strong scale effect was observed), and furthermore, in Sumelka et al. (2016) it was presented that s-FCM can correctly mimic the Born-Von Karman (BK) lattice (discrete system). This result are crucial for correctness of physical interpretation of the results obtained below.
The governing equation for a 1D fractional elastic body, under the assumption of a variable length scale (Sumelka 2017), is stated as the following integro-differential equation where E denotes the Young modulus, f is the length scale (herein known function), U denotes the displacement, and b is the body force. We distinguish two types of boundary conditions: • for both ends of a 1D body displacements are prescribed • for the left end displacements and for the right end forces are given where x 0 and x r denote points on the boundary -see Fig. 1.

Internal points
The approximation of (4) follows the expressions proposed in Odibat (2006), Leszczyński (2011). Hence, the spatial discretization of a 1D fractional body of length L into r equal subintervals x is assumed (see Fig. 1) by the following rule with x 0 = 0, r ∈ N where N denotes the set of natural numbers.
Next, for a discretization point x i the left-sided derivative is approximated by: where U (n) (x a i j ) denotes a classical n-th derivative at x = x a i j ; and by analogy for the right-sided derivative we obtain: Finally, assuming that α ∈ (0, 1], n = 1 in (9) and (10) and x = h the behaviour of the i-th node is governed by (cf. Fig. 1): , and (·) , (·) denote the first and the second order derivatives which are approximated using the classical central finite difference scheme. It should be emphasised that for α = 1 (n = 1) and E = const.

Boundary points
Without loss of generality, we only consider the case for both ends clamped. It follows from the considerations presented in Sumelka (2017) 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 in point x 0 , and by analogy central and backward schemes in point x r one has (for Fig. 2 General scheme of the analysed structure points x 0 , x 1 , x 2 , x r−2 , x r−1 , andx r one has p = p min = 2, whereas for all other points parameter p should guarantee that f | i = p x, is smaller than the distance to the closest boundary):

Algebraic problem
One can write the set of r − 1 equations obtained from (11) together with constraints (13) and (14) as a global algebraic problem, namely where K denotes the fractional stiffness matrix, q governs the unknown displacements in the discretization nodes (U 1 , ..., U r−1 ), and Q is a global force vector (b 1 , ..., b r−1 ). Note that in order to obtain (15) (11) a backward scheme is applied.

Structure
The nano-device undergoing corrosion is modelled as a 1D body (Fig. 2) which is fixed at both ends and loaded with a mass load calculated as: where ρ i is density and g i is a field intensity function in the i − th node. Three types of field intensity functions are assumed and compared in the paper: It is assumed that the parameters of the fractional body, and α, are identified for the initial (non-corroded) structure. Moreover, for the sake of simplicity, the non-dimensional case is considered, so L = 1.

Capturing deterioration
The following fundamental assumptions are made for the structure (members/grains) corrosion: • deterioration is related to density reduction, • longitudinal stiffness depends on density, • mass lost is known (optional).
The first assumption considers an effect of deterioration on a microstructure. It states that the corrosion does not affect the topology. The members/grains of a microstructure are proportionally deteriorated and the microstructure is not damaged -all elements remain and are connected as in the initial state (Fig. 3). Taking into account that the length scale f and material order α are associated with a piece of a body where some repeatable pattern is observed (RVE) (Askes and Aifantis 2011), it allows us to keep f and α at the same level for the initial and the deteriorated structure. Following this conclusion, it is assumed that these parameters can be identified for the initial structure and then used in the further identification.
The second assumption, in turn, determines that density is a deterioration indicator. As a consequence, identification of density is enough to assess the condition of a device. In the general case, the density distribution can be free; however, in order to parameterize it using a reasonable number of variables, a spline function is utilized. Let us   Fig. 3. Note that the knot points indexed by 1, ..., i..., s are not the same as the ones used for the body discretization (cf. Section 2.3.1). Their number is much lower, which allows us to reduce number of parameters. In the current study eight knots were assumed as a trade-off between computation time related to the number of variables and capability to capture gradients in the density distribution.
The third assumption says that the structure deterioration influences stiffness in relation to density. The exact relation, however, depends on many factors. Among others, structure topology can be mentioned. Thus, the relation should also vary in accordance to the length scale f and material order Fig. 7 Displacement for corroded structure subjected to mass load α. In the presented work the following power relation for Young modulus is proposed: where E and ρ are the initial stiffness and density, respectively. The exponent c in the above equation depends on the structure and has to be identified. (18) is conceptually close to SIMP (Bendsøe 1989). However, in the presented work this relation is not only a numerical endeavour to push density to marginal values for a point in domain. It has a physical meaning and in this way it can be compared to the relations like the ones in Rho et al. (1995) which bind the computer tomography output with bone stiffness. The final assumption refers to optimization requirements. The mass stabilises the process and leads to unique solutions. If the mass cannot be measured, it should be considered as one more parameter in optimization leading to minimisation of the original formulation (cf. Section 3.4).
Note that because of the above assumptions, in the current study the corrosion phenomenon is reduced to a density change and a corresponding stiffness modification. Further investigations of corrosion effects should be done to verify the approach for any structure.

Numerical model
The behaviour of the analysed structure is described using s-FCM according to (11). In all analysed computational scenarios the body is discretized into 100 subintervals, what leads to x = const. = 0.01.
In order to eliminate the virtual boundary layer (Sumelka and Łodygowski 2013) and its unknown influence on the results, the concept of variable length scale is utilised (Sumelka 2017). The length scale is linearly reduced approaching to the structure ends from max f to 2 x. Considering additionally (13) and (14), no information on displacements outside the body is required. The exact distribution of the length scale f (x) is fixed with reference to its maximal value as presented in Fig. 4.
Note that except for f , for each discretization point the values of ρ i and E i also vary. They are dependent on values  (17)). Displacement graphs collect results for all analysed length scales max f and material orders α in knot points of the spline function representing density distribution (cf. Section 3.2) as well as an exponent c in (18). Both the knot points and exponent c change during the optimization what leads to the relations below: ρ 1 (k), ...,ρ j (k), ...,ρ s (k)), where k is an iteration number. A dedicated procedure in Python was developed in order to build governing equations dynamically for the particular points and then assemble them into a system of linear equations and solve (cf. github.com/szajek). The main idea behind the library was to create a template of a governing equation as a combination of elements (etc. numerical approximations of a differentiate, variables, constants, ...) bound using lazy operators (Watt and Findlay 2004). This concept allows us to easily capture varying definitions of elements and changing values within a body, which is an important case in the analysed problem. The system of (15) is solved using an LU decomposition (Strang 1980) with partial pivoting and row interchanges. The final output is a series of displacements along the body U i .

Optimization problem
The following constrained optimization problem was formulated: subject to: whereŪ and U denote measured and computed displacements, respectively; p, in turn, represents the design parameters vector -density values in knot points of the spline function,ρ i (see Fig. 3) and exponent c in the mass to stiffness relation (cf. (16)) and is defined as: while factor m denotes mass ratio as follows: Fig. 9 Density distribution for field intensity g 1 (x) and varied fractional material parameters where denotes unit cross-section, whereas M and M 0 indicate the mass of the corroded structure and the initial structure, respectively.
Using exterior, static penalty function the problem was converted into an unconstrained one as follows: where: If there is no information on mass, the optimization problem presented in (23) has to be wrapped into another one as follows: min keeping the same constraints. This approach, however, is not used in this paper.
Due to discretization used (cf. Section 2), the objective function and constraint functions were calculated based on the numerical approximations below: where For all parameters the boundary limits were defined. The boundary values of ρ are ρ lower = 0 and ρ upper = 1, while the c limits equal c lower = 0, c upper = ∞. To compute the penalty terms, the preliminary studies based on trials and errors approach were done, and the following coefficients were estimated: p 1 = 10 −2 , p 3 = 10 −4 and p 2 = p 4 = 2.5.

Algorithm
The problem described by (23) was solved by coupling the numerical model with an optimizer. The general overview of the procedure is presented in Fig. 5. As an optimizer Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (Byrd et al. 1995;Zhu et al. 1997) implemented in SciPy package was chosen (SciPy Developers 2017). The most important settings used are assumed as follows: maximum number of iterations as 15000, minimal relative change of an objective as 2.22e-9 and step size used for numerical approximation of the jacobian as 1e-8.

Test data
In order to test the methodology, an example density distribution along the body,ρ(x), was assumed as presented in Fig. 6. Additionally, a linear relation between material stiffness and density was assumed, E(x) = ρ(x). Next, it was assumed that the initial structure can be characterised by length scale f = 0.20 and the derivative order α = 0.5. The displacements for the deteriorated structure subjected to a mass load were calculated and considered to be the ones measured in reality (see Fig. 7). Due to the primary goal, capturing sensitivity to fractional parameters, no measurement noise is considered. The mass reduction ratio due to corrosion for the example structure was as follows:

Results
The problem described by (23) in the variant with a known mass was solved using the algorithm presented in Section 3.5. The optimization was carried out for seven different length scale distributions ( max f ∈ {0.05, 0.10, 0.15, 0.2, 0.25, 0.30, 0.35}, see Fig. 4), six orders of fractional continua (α ∈ {0.4, 0.5, 0.6, 0.7, 0.8, 0.999}) and three load scenarios depending on the field intensity function, g(x) (cf. (17)).  Figure 8 presents the objective value and displacements obtained for the identified density distribution. Note that regardless of the field intensity function, g(x), and non-local model parameters, all optimization analyses were successful in the sens of fitting the displacements for the test data. The maximal objective value, O u , is lower than 1.4e − 6 for the worst case (g 2 (x), max f = 30 and α = 0.4). However, only in 3 out of 147 cases the objective value is higher than 1.0e − 6, what leads to an average error for a single node and field intensity function g 3 (x) (the lowest maximal displacement) as below: Note that the above error also includes penalty terms (cf. (23)); therefore, the displacement discrepancy is even smaller. It proves the agreement between the test and computed displacements, which can be also confirmed → 0 (dimensions of the body become considerably larger than the characteristic length scale) once more the classical local solution is obtained. The third observation is that the influence of fractional material parameters increases with the complexity of the field intensity function. When g 1 (x) = 1 or g 2 (x) = sin(x/L · π) even a relatively big length scale and lower order of material do not influence the density distribution strongly (cf. Fig. 9). The effect becomes visible for g 3 (x) = sin(2x/L · π) and max f > 0.15 and leads to the qualitatively different solutions (the selected density distributions are visualised in the 3D form in Fig. 12). The intensity function, g(x), defines distributed load, b(x). Thus, in this meaning, it can be concluded that non-local effects in the 1D case are more pronounced when the displacement field changes strongly due to a complex load (far different from a constant one). The final observation concerns confirmation of the methodology, the optimization procedure and its implementation.
By analogy to the density distribution, Figs. 13, 14 and 15 present the stiffness change along the body for varied max f and α. Like previously for α → 0, the change in the solution is more pronounced. As in the density distribution, smooth passage to the classical (local) solution is obtained for α → 1 or max f → 0. The same is also an effect of the field Fig. 14 Stiffness change along the body for field intensity g 2 (x) and varied fractional material parameters Fig. 15 Stiffness change along the body for field intensity g 3 (x) and varied fractional material parameters intensity function (distributed load). The more complex load is, the more evident the non-local effects are.
The last results concern the density-to-stiffness relation. Figures 16, 17 and 18 present the exponent c change (see   for the test data (the point on graphs is indicated with a star). A significant relation between exponent c and material parameters max f or α is observed. Regardless of the field intensity function g(x), exponent c is sensitive to α. The sensitivity is the biggest when max f → 0.35 and in the extreme case (g 2 (x)) the difference reaches almost 0.5 for α = 0.999 and α = 0.4. When g 1 (x) or g 2 (x) is considered, exponent c becomes lower if max f goes from 0.2 (the value for the test data) to higher values. In the case of g 3 (x) this trend is visible for the whole range of max f . The only exception is the result for α = 0.4.
The last observation, based on the obtained results, is that the inverse problem in (23), where exterior static penalty function is used for both mass constraint and (what is more important here) the design parameters limits is in numerical sense well-posed. For each configuration of the fractional parameters ( f , α) the solution was found. The results for a classical model (when α → 1) were confirmed with the analytical solution (which is wellposed) and moreover the inverse problem solution changes continuously (smoothly) with the fractional parameters for both density distribution ρ (cf. Figs. 9, 10 and 11), as well as exponent c (cf. Figs. 16, 17 and 18). Additionally, for the proposed test data (Fig. 6) the inverse problem converges to the same results as the assumed for the forward analysis (Fig. 7) in all cases.

Conclusions
The paper presents the study of mechanical properties of 1D deteriorated bodies exhibiting the scale effect. The overall problem has been formulated in the framework of the space-Fractional Continuum Mechanics as an inverse problem.
Steering parameters of the optimization task were: density of material, order of material, and material length scale. The obtained results can be crucial, especially from the point of view of durability of microstructured/nanostructured materials/devices whose popularity in real life applications is growing rapidly and manufacturing is becoming easier in view of constant progress of 3D printing techniques. Furthermore, the outcomes show the meaningful role of proper modelling of non-local bodies, specifically taking into account their maintenance which can be in the near future one of the most challenging engineering tasks.
The most important general observations are: • classical continuum mechanics is not able to properly predict the deterioration of 1D non-local bodies -both qualitative and quantitative difference with respect to s-FCM outcomes was observed; • the obtained results are sensitive with respect to the length scale to body dimension ratio, being more emphasised for ratios closer to one; • the order of fractional continua acts as a scaling factor, and makes the difference more pronounced approaching to zero; • for more complex displacement field/loading, considerable increase of outcomes sensitivity is observed.
The paper is the first step in research on optimal space-Fractional Continuum Bodies, and it is a base for further investigations including optimal distribution of the length scale and order of s-FCM continua, including anisotropy (Sumelka 2016).