An Investigation of Radial Basis Function Method for Strain Reconstruction by Energy-Resolved Neutron Imaging

The main objective of the current work is to determine meshless methods using the radial basis function (rbf) approach to estimate the elastic strain field from energy-resolved neutron imaging. To this end, we first discretize the longitudinal ray transformation with rbf methods to give us an unconstrained optimization problem. This discretization is then transformed into a constrained optimization problem by adding equilibrium conditions to ensure uniqueness. The efficiency and accuracy of this approach are investigated for the situation of 2d plane stress. In addition, comparisons are made between the results obtained with rbf collocation, finite-element (fem) and analytical solution methods for test problems. The method is then applied to experimentally measured continuous and discontinuous strain fields using steel samples for an offset ring-and-plug and crushed ring, respectively.


Introduction
Consider a solid material sample composed of similar crystals with a uniform distribution of directions. There are three times as many crystallographic planes in each crystal as there are Miller indices. The sample if subject to a narrow neutron beam with a known distribution of kinetic energy. The neutrons behave as waves with the energy inversely proportional to the wavelength, and this wavelength is a similar magnitude to the crystallographic planes [1,2]. Parallel crystallographic planes with d separation scatter an incident wave if the wavelength is a multiple of 2d sin θ. If we detect only the transmitted neutrons, the scattered neutrons will be lost. If we look the number of transmitted neutrons as a function of wavelength, then at the wavelength 2d the transmission leaps dramatically. This jump is labelled as the Bragg edges. If the material is subjected to a ε linear elastic strain, the separation between planes with a normal n vector is in proportion to the strain component n · ε · n. For more information and additional sources, see [3][4][5]. If the strain were consistent along the neutron path, the Bragg edge would be shifted. Within ideal conditions, a linear wavelength function multiplied by the unit step function can model the Bragg edges. Abbey et al. [6] had previously explored whether strain could be reconstructed from the average of each ray shifts of a Bragg edge. Lionheart and Withers [7] noted that the transverse ray transformation is feasible for X-ray strain measurements.
Given the longitudinal ray transformation (LRT) data plane by plane, an attempt could be made to recover the strain from the ray measurements in that plane using appropriate additional a priori information [8]. The Lamé system could be represented by finite-element or finite difference methods provided that the Lamé coefficients are known. A linear system could be developed, and the least square solution determined together with a discretization of the LRT and prior knowledge. This method has been adopted by Gregg et al. [5] for axisymmetric objects, Wensrich et al. [9] for synthetic data, and Henriks et al. [10] for a generic case experimentally.
We consider here the measurement of the Bragg edges based on the calculation of the propagation spectrum using flight time or energy-resolved technologies. We note that other approaches exist [3,[11][12][13]. The major contribution of neutron transmission strain tomography is the optimal use of the incident beam flux to reduce the data collection time. Recent developments in detector technology allow the acquisition of high-resolution strain images with Bragg neutron transmission [14]. The sample is placed on a rotating plate between the source and a micro-channel plate detector (MCP) [15,16]. This detector is pixelated with a spatial resolution of up to 512 × 512 pixels of up to 55 µm. This imaging raises the possibility of strain tomography, and efforts have been made over the past few years to solve the resulting tensor reconstruction challenge-special cases such as axial symmetry revolved around reconstruction techniques [5,17]. A well-known fact about this inverse problems is that it is not well-posed. Therefore it is not possible to reconstruct the strain without the inclusion of additional conditions or without easing assumptions [10]. The question of how to better execute these additional criteria was discussed and various approaches were used in the literature [5,9]. Recently, reconstruction has been accomplished by incorporating either equilibrium or compatibility conditions. The main distinction between these approaches is the set of basis functions representing the strain field [6,10].
Recently, considerable attention has been given in the science and engineering community to the concept of using meshfree methods as a modern tool for numerical solution for partial differential equations. As stated in [18,19], the growing interest in these methods is due partly to their high versatility, particularly in the case of high-dimensional problems. Conventional grid-based approaches such as finite difference (FDM) [20] and finite-element (FEM) [21] methods have inherent problems, namely a mesh generation requirement. Mesh generation is not a simple job, particularly for domains with irregular and complex geometries requiring additional mathematical transformations that can be as costly as solving the problem. Meshfree avoid this problem since the conventional definition of the concept does not include the meshing procedure. This makes it easy to understand the behavior of complex solids, structures and stress and strain fields arbitrarily distributed in the domain without any constraints. Moreover, we do not need to consider the connectivity between the nodes while applying the meshfree process. For a survey of various meshless approaches and underlying findings, we refer the reader to [22][23][24][25].
The paper describes a meshless process by which the elastic strain of a series of Bragg-edge strain measurements can be reconstructed tomographically. Because of the attractive properties of radial basis functions (RBF), the approach has some advantages over previous approaches. Especially its ability to represent extremely spatially varying functions accurately. The proposed algorithm is evaluated on a two-dimensional simulated beam data for which a thorough comparison of the FEM and RBF methods is shown. It is shown to be able to reconstruct a tensor field, after [26,27] imposed equilibrium conditions. The framework is then extended to real experiment measurements for the steel geometry of an offset ring-and-plug and a crushed ring.
The paper has the following structure: Section 2 discusses strain measurement through the Longitudinal Ray Transformation. Section 3 offers an overview and algorithm of the meshless method for strain field reconstruction and presents the proposed algorithm structure and the finite-element approach. Section 4 addresses the validation of the algorithm using a cantilevered beam and the key difference between the FEM and the solution being suggested. Finally, the experimental findings with the ring-and-plug and crushed ring steel samples are addressed in Section 5. Section 6 summarizes the work in this paper and lists the work to be done in the future.

Strain Measurement
Here we discuss the recent experimental method that provides information on the average strain component for the set of incident neutron beam [11]. This approach can have a significant influence in several experimental fields in mechanics. Bragg-edge strain tomography guarantees a unique structured approach for the study of materials with residual stress field over practical scales. The Bragg edges, as mentioned above, are generated by backscattering radiation. As a result, relative shifts in their position calculate the average normal strain of the sample towards the ray [28]. This challenge revolves around the LRT inversion, which represents an effective measurement process model. Although this is a three-dimensional problem in general, for convenience, this article takes only two dimensions into account. The average strain within the body measured by Bragg-edge neutron transmission is often idealized as a longitudinal line integral, namely LRT, which captures the average strain component along the line towards the normal unit direction n = (n 1 , n 2 ) T . Γ ε is a form of ray transform specifically known as LRT. Hence, LRT shown in Figure 1 with reference to the coordinate system and geometry is defined as where ε ij is the (i, j)th component of tensor strain field ε ∈ R 2×2 , which is mapped to an average normal component of a strain in the direction ofn. A ray projection enters the sample at the location x a = (x(s, a), y(s, a)), where s is the ray propagation distance coordinate and L is the length of the ray within the sample at the particular angle. Every pixel of a strain image measures the average normal strain in the ray direction through the thickness of the domain, say Γ ε . It is based on the overall change of the gauge distance within that ray [6,7]. Measurements are taken in each orientation, ϑ i , where the profile is calculated in the form Γ ε (x, y, ϑ i ).

Meshless Approach Overview
One of the most common techniques among meshfree methods is the Radial Basis Functions (RBF). The formal definition of a radial function is given below.
where we choose || · || to be the Euclidean norm on R d .
The proposed approach is similar to the three-layer feed-forward RBF neural network [29,30]. The first layer is the network inputs, the second a hidden layer consisting of a collection of nonlinear RBF activation units, and the third corresponds to the final network output. In Radial Basis Function Networks (RBFN), activation functions are usually implemented as Gaussian functions. An overview of RBFN structure is shown in Figure 2. To illustrate how the RBFN works, assume that we have a set of D data with N patterns of (x p , y p ) in which x p is the input of the set of data and y p is the actual output. The ith activation function in the hidden network layer can be determined from the gap between input pattern x and centers i based on Equation (2) Here, || · || is the Euclidean norm, c i and r i are the center and width of the hidden neuron i, respectively. Then the contribution of the k node from the network output layer can be computed with Equation (3)

Methodology
The addition of polynomials to the radial basis functions will ensure consistency and guarantees the approximation properties of the scheme. However, it complicates the system of equations making them difficult to enforce and contribute extra to computational cost and time. In the proposed method, 2 × 2 symmetric strain tensor field ε(x) is composed of two parts: radial basis functions R i (x) and polynomial basis functions P j (x) [25,31]. Hence, the local approximation is given by where m is the number of polynomial terms and n is the number of nodes in the local support domain Ω. Unknown vectors to be calculated are a and b. The addition of polynomial terms is an additional prerequisite. The following conditions are generally imposed to ensure a unique approximation [32] where m is a number of augmented multivariate polynomial basis functions, see Table A1. We assume P j (x) = 1 + x + y + xy + x 2 + y 2 and R i (x) = exp(−α 2 ||x − c i || 2 ) where α is the shaping parameter. Hence, above Equation (4) is expressed in the matrix form where U is the vector of function values: U = [ε 1 lk , ε 2 lk , · · · , ε n lk ] T , R is the matrix of RBFs, P is the matrix of polynomial basis function and a, b are the values of unknown coefficients (Radial and polynomial respectively). We note that to obtain the unique solution of Equation (4), the constraint condition (5) should be combined with Equation (6) as follows The unknown vector in Equation (7) is obtained by the inversion of the matrix F = R P P T 0 .

Weak Form for Longitudinal Ray Transform Integral
Let us consider a two-dimension strain reconstruction problem in domain Ω whose strong form is given by Equation (1) and generalized weak form is obtained as Enforcing (8) to pass through all node values in the domain Ω, the following equations are obtained where C is the coefficient matrix; U is the vector of unknowns and Γ is the vector of integral values. Now, rewriting (4) in a matrix form to find weights where N is the total number of nodes in the domain Ω, and Hence, and Now, using Equations (14)-(17) in the previous system (9), we get which lead to an optimization problem as follows To ensure uniqueness, equilibrium equations are imposed as follows: This will give us another system of equations The least square algorithm, described in [33], is well suited to this task. In general, least square computes a solution ε to problems of the following form which is solved using least square optimization intrinsic function "lsqlin" in MATLAB.

Shape Parameter
For the precise implementation of the framework, the shape parameter evaluation is important [34]. Therefore, careful attention should be given for determining shape parameter value. Detailed discussions regarding values for the shape parameter can be found in Fasshauer [35], Franke [36], and Hardy [37]. Most formulations for shape parameters depend on the distance, d, and node numbers, N. Other formulations display a changing shape parameter depending on the location of the node. In general, the shape parameter depends upon several variables such as the data point distribution, the basis function, and a computer precision (see [38]). The effect of the RBF Gaussian shaping parameter α mentioned above is shown in Figure 3.

Algorithm
Here, we discuss the difference between two methods; namely FEM and RBF. Detailed pseudo-code algorithm is given below Algorithms 1 and 2 for both the methods, respectively.  The complicated meshing rules are not required while performing RBF, and here lies the advantage of the proposed method over the framework of finite elements. The use of meshless methods over the finite-element approach for reconstruction tensor strain field [39][40][41] has many other benefits. We can quickly increase the support domain, change basis functions, adopt nodal densities, and many more. Those characteristics can be useful in the case of irregularities like sharp discontinuations or other situations in complex simulations. In the next section, we talk about this difference between methods with a simulated example.

Cantilevered Beam
Here we illustrate the practical performance of the proposed algorithm, starting with the well-known two-dimensional cantilevered beam problem that was previously successfully determined by other authors [9]. In this case, we take the [0, 20] × [−5, 5] a two-dimensional strain domain into account with a 2 kN load over the beam. The beam's material properties reflect common steel, while other parameters are listed below.
The data consist of 9 projections evenly spaced in [0, 170] • with 200 lines each, yielding a total of 1800 measurements. Figure 4 shows two types of nodal mesh (uniform and non-uniform) used for beam strain field reconstruction using RBF. Both nodal mesh consist of 400 nodes including boundary points. Figure 5 shows the reconstruction of strain field using RBF, whereas Figure 6 shows the analytical solution of the cantilevered beam. The plane stress Saint-Venant approximation to the strain field is [42]: where I represents the second moment of area, L the beam length, h the beam height, t is the thickness, E is Young's modulus, ν is the Poisson ratio, and the dimensions are shown in Table A2.

Error Analysis through Simulation Results for RBF Methods
A L 2 error analysis of strain field associated with both FEM and RBF methods is shown in Figure 7. Started with 9 projections evenly spaced in [0, 170] • with 10 lines each and ended with 700 lines each 9 projections. Both RBF case results shows better agreement than FEM in the final stage. However, faster convergence is seen in the case of uniform RBF. Figure 8 shows the RMS error between the reconstructed and analytical strain field as we change Gaussian function shaping parameter value.
where N is the total number of nodes, T i is the true strain value and R i is the reconstructed strain value at the particular node i. We found out the proposed reconstruction algorithm extremely useful for the reconstruction of the strain field. The simulation results imply that the reconstruction algorithm will converge to an appropriate reconstruction provided measurements over 180 • of a sample are calculated. Discretization problems and numer-ical errors will undoubtedly lead to noisy reconstruction. As the number of projections increased, rapid convergence to the analytic solution was observed. The constrained optimization problem (23) can be solved using a variety of approaches. Our algorithm uses the "lsqlin" inbuild function in MATLAB. This convergence gives us confidence that in the face of real experimental uncertainties, the algorithm will converge to an analytic solution.

Experimental Demonstration
The RBF approximation methods presented above have been tested on synthetic and real data sets. We are now testing the algorithm on the experimental data previously studied for the offset ring-and-plug, and the crushed ring sample [26]. The crushed ring sample was 14 mm thick with 23 mm outer diameter and, 10 mm inner diameter. To test the algorithm for continuous and discontinuous stress fields, specific samples were selected. The ring-and-plug sample is designed to test the discontinuous strain field and the crushed ring for continuous strain field testing. The offset ring-and-plug sample geometry is shown in Figure 9, and the crushed ring is shown on the left side of Figure 10. The sample described included a total interference of 40 ± 2 µm generated by cylindrical grinding. Although the plastic deformation of the crushed ring sample by 1.5 mm on the diameter was achieved using 8.4 kN rigid steel plates. A brief description of the model and experimental settings is provided below and can be found in [26]. Two EN26 steel bars have been heated to relieve stress and provide a uniform structure before assembly. The final test hardness of the sample was 290 HV for both models. The strain profile was calculated on RADEN along with an MCP detector (512 × 512 pixels, 55 µm per pixel) at 17.9 m from the beam source of power 409 kW (January 2018). In the Japan Proton Accelerator Research Complex J-PARC, Japan, the neutron-energy-resolved imaging instrument was used to obtain the relative shifts of the Bragg edges corresponding to the lattice plane of the samples. Neutron strain scanning was performed on KOWARI, a residual stress diffractometer at the Australian Nuclear Society and Technology Organization (ANSTO), Australia, to verify our reconstruction independently. This strain-scan provide measurements of the three in-plane components of strain over a mesh of points within the sample using a 0.5 × 0.5 × 14 mm gauge volume. Sampling times on KOWARI were based on the providing of strain uncertainty around 7 × 10 −5 , which required about 30 h of beam time per component. The RADEN sampling times, however, were based on the 1 × 10 −4 statistical uncertainty. A total of 50 profiles were measured at golden angle increments in the angle with a sampling time of two hours per projection.
A two-dimensional scattered data point net is used to discretize the domain of the sample, as shown in Figure 9. Two different types of node discretization have been considered: uniform and non-uniform (random). Nodal set can be seen in Figure 11 and reconstructed strain fields can be seen with each discretization in Figures 12 and 13 respectively. As stated earlier, for the simulated problem, the reconstructed strain field did contain noise which can be from either experimental measurement or numerical discretizations.     To summarize, we proposed a meshless approach using RBF. They are the natural generalization in a multivariate sense of univariate polynomial splines. It works well for high-dimensional arbitrary structures, and no mesh requirement is the main benefit of this form of approximation. A RBF is a function, the value of which depends only on the distance between the points. Due to the use of distance functions, one can conveniently deploy RBFs to recreate the surface using data points in 2D, 3D or higher-dimensional spaces.
The Least square solution was found for 3790 nodes for the non-uniform and 3795 nodes for uniform ring-and-plug discretization from 22,678 projections. Figures 12 and 13 shows the reconstructed field of ring-and-plug with uniform and non-uniform node distribution respectively in terms of three Cartesian strain components namely; ε xx , ε xy and ε yy . The reconstructed strain field was found for 1352 nodes and 1086 quadrilateral elements via FEM approach and 1352 nodes for crushed ring discretization via RBF approach from 20,664 projections. Figures 14 and 15 shows the reconstructed strain field of crushed ring with FEM and RBF respectively in terms of three Cartesian strain components. Figures 16 and 17 show the corresponding validated strain field from KOWARI data. The effectiveness of the suggested approach is not based solely on sample structures and can, thus, be extended to three-dimensional sample bodies. Due to radial basis function computation, the real challenge is the computational cost and time since the size of the problem is quite large. Figure 14. The reconstructed full-field two-dimensional strain components (ε xx , ε xy , and ε yy ) with FEM. Figure 15. The reconstructed full-field two-dimensional strain components (ε xx , ε xy , and ε yy ) with RBF. Figure 16. Two-dimensional strain components ε xx , ε xy , ε yy respectively with KOWARI Data. Figure 17. Two-dimensional strain components ε xx , ε xy , ε yy respectively with KOWARI Data.

Conclusions
From Bragg-edge neutron images, an elastic strain tensor field reconstruction algorithm has been presented. Unlike previous algorithms such as FEM, our suggested methodology is capable of reconstructing residual strain field as no elastic strain compatibility is assumed. Simulated and experimental data from KOWARI and RADEN have shown this process. The analysis revealed a strong agreement with the strain maps measured using the constant diffractometer Kowari. A complete strain field tomography can now be achieved using physical constraints as an equilibrium without thinking about two-dimensional mesh. This approach opens up more study, including three dimensions, for future studies. The framework taken allows us to concentrate on the highly stressed region of the sample using adaptive node discretization, which can be accomplished by calculating the gradient over the sample.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
We appreciate the support of Australian Research Council (ARC) through a Discovery Project (DP) Grant No. DP17010 2324. We would like to extend our appreciation to the respective user-access programs of J-PARC and ANSTO (J-PARC Long Term Proposal No. 2017L0101 and ANSTO Program Proposal) for providing access to the RADEN and KOWARI instruments and The University of Newcastle for funding R.A's scholarship.

Conflicts of Interest:
There are no conflicts of interest to declare.

Appendix A
The LRT measurements for Bragg-edge neutron transmission with respect the sample geometry and coordinate system Figure 1 is shown below: where x a = (x a , y a ) is the entry points of the ray to the sample. Now, substituting radial basis function expression for ε 11 , ε 12 , and ε 22 determined previously in Equation (4) in terms of linear combination of R i and P j we get a closed form of line integral as: assuming, R i as the Gaussian function and P j as the second order polynomial written as R i = exp(−α 2 ||x − c i || 2 ); α = shaping parameter P j = 1 + x + y + x 2 + y 2 + x y where A ij = L a n 2 i + n 2 j , b ij = a 2 n 1 i + n 2 j , φ ij = exp(−b ij (δ ya n i + δ xa n j ) 2 ); δ ia = i a − c i , B ij = (δ ya n i + δ xa n j ) 2 ) √ b, C ij = A ij + B ij , ψ ij = erf(B ij ) − erf(C ij ), D ij = L 2 (n 2 i + n i n j + n 2 j ), E ij = 3L/2((2x a + y a + 1)n i + (x a + 2y a + 1)n j ), F ij = 3(y 2 a + x 2 a + x a y a + x a + y a + 1)