An inverse method for determining the spatially resolved properties of viscoelastic–viscoplastic three-dimensional printed materials

A method using experimental nanoindentation and inverse finite-element analysis (FEA) has been developed that enables the spatial variation of material constitutive properties to be accurately determined. The method was used to measure property variation in a three-dimensional printed (3DP) polymeric material. The accuracy of the method is dependent on the applicability of the constitutive model used in the inverse FEA, hence four potential material models: viscoelastic, viscoelastic–viscoplastic, nonlinear viscoelastic and nonlinear viscoelastic–viscoplastic were evaluated, with the latter enabling the best fit to experimental data. Significant changes in material properties were seen in the depth direction of the 3DP sample, which could be linked to the degree of cross-linking within the material, a feature inherent in a UV-cured layer-by-layer construction method. It is proposed that the method is a powerful tool in the analysis of manufacturing processes with potential spatial property variation that will also enable the accurate prediction of final manufactured part performance.


Introduction
From the literature [28][29][30][31][32][33], it can be seen that a key step in the use of IFEM is the selection of an appropriate material constitutive model. Only a material constitutive model that can accurately describe the mechanical behaviour of the material is suitable; otherwise, it is impossible to achieve a good fit between FEA and test data. Although linear VE constitutive relations have been used to describe the mechanical behaviour of polymer materials under indentation, tests and analysis are mostly limited to the loading and holding (or creep) stages [34][35][36] of an indentation test, without considering the unloading stage. It was found in this investigation that not all the linear VE parameters extracted from the loading and holding stages of a test fit well with the unloading stage. The predicted recovery during the unloading stage tends to be much higher than the tested recovery when using the material parameters determined from loading and holding stages. This indicates that there could be viscoplastic (VP) deformation during loading and holding stages besides the viscous and elastic elements. Therefore, a more advanced material constitutive model should be considered if material responses under complex loading scenarios that include unloading are to be accurately modelled.
In this paper, nonlinear viscoelastic (NVE), viscoelastic-viscoplastic (VEVP) and nonlinear viscoelastic-viscoplastic (NVEVP) material constitutive models are developed and compared with the more standard linear VE material behaviour. An IFEM is then developed and used to investigate which material model most accurately describes the mechanical behaviour of an RIJ polymer. The unloading stage is included in the IFEM as this is necessary to ensure the resultant material model and parameters adequately model the polymer under loading conditions seen in service, which would be expected to include loading, constant load and unloading at various times. The proposed method is then tested for its efficacy in determining material model parameters under conditions of 3D RIJ printing. Using a commercial 3D printing system, samples were manufactured and tested with nanoindentation. The IFEM was then used to determine the most appropriate material model, the optimal material parameters and the spatial variation of these parameters in the sample. The aims of the paper, then, are to introduce an IFEM which can be used to determine complex material model parameters from indentation data and to illustrate how these can be used to investigate property variation in a 3D printed (3DP) part, thus providing insight into the manufacturing process and enabling prediction of the mechanical performance of the printed part.

Methodology
The proposed method for accurately determining the spatial variation of complex material properties requires a method of experimentally determining spatially resolved mechanical properties, constitutive models that accurately represent the behaviour of the material under realistic loading conditions and a method of determining the best-fit material model parameters. In the sections below, the experimental methods are described first, followed by the IFEM used to determine the best-fit model parameters, and finally a number of constitutive models identified as potential candidates for characterizing the polymeric material used in the experimental work are developed.

(a) Sample preparation
Samples were prepared using an Object260 Connex 3-D Printer system (Stratasys Ltd). The ink used was VeroClear Fullcure 810, which is a transparent acrylic polymer that is cured after deposition through exposure to UV light. The samples produced were 10 × 40 × 5 mm solid blocks (x × y × z; figure 1a). These samples were removed from the build plate and sectioned using a hand saw to a size of around 10 × 5 × 5 mm. The sectioned samples were then mounted in EpoFix cold-setting resin (Struers Ltd) for nanoindentation. The samples were mechanically machined to create a flat surface, and ground to make the two surfaces parallel. The sample test surface was then polished, finishing with a 1 μm polish. All samples were wrapped with tissue and aluminium foil, to prevent UV exposure, and placed in a desiccator cabinet located in an air-conditioned room in order to minimize any environmental effects prior to testing.

(b) Nanoindentation
Indentation tests were conducted on a NanoTest 600 machine (Micro Materials Ltd). The maximum loads applied during testing ranged from 25 to 75 mN, with a hold time at maximum load of 500 s. Both the loading and unloading times were fixed at 50 s, meaning that the loading and unloading rates were dependent upon the maximum load applied. The spherical diamond indenter used was 50 μm in radius. All tests were conducted in an enclosed chamber held at 25 • C.
The coordinate system used in this paper is shown in figure 1a. Surfaces are identified with the following terms: the top surface is the final printed layer of material, creating an xy-plane at z = 0. The side surfaces are the outermost surfaces in the yz-plane, i.e. yz-planes at the minimum and maximum values of x. These are natural surfaces that are formed during the manufacturing process. The normal to this surface is also normal to the print direction. Cross sections were formed by mechanical sectioning in the xz-plane. In all the samples, the cross section is by definition at y = 0. The normal to this surface is parallel to the print direction. Printing was performed in a rastering motion in the print direction.

(c) Forward finite-element model
The following assumptions were used to establish a working FE model for use in the inverse analysis.
(1) The material is locally homogeneous, i.e. any changes in the material properties are small within a length scale that is longer than the characteristic length scale in the problem, which in this case can be considered to be the contact radius in nanoindentation. Experimental investigation of the rate of change of properties as a function of spatial distance indicated this to be a reasonable assumption for the size of indents used. (2) The extent of the sample can be considered infinite. In this case, the sample size is two orders of magnitude greater than the indenter size, suggesting that this is appropriate (confirmed by computational investigation of boundary conditions).  An axisymmetrical FE indentation model with a 50 μm analytical spherical indenter was, hence, set up, with boundary conditions as shown in figure 1b. The representative sample geometry was a 500 × 500 μm (radius × height) cylinder with indentation performed in alignment with the axis of the cylinder. Tests confirmed that this geometry was sufficient for boundary insensitive results. The mechanical contact between the indenter and the sample material was assumed to be frictionless, tests indicating results to be relatively insensitive to friction at low values. The mesh is shown in figure 2b and it can be seen that a fine mesh was applied in the indentation location. Mesh convergence tests were conducted to determine appropriate mesh refinement in the contact area.

(d) Inverse finite-element modelling technique
Although nanoindentation tests have many advantages over traditional mechanical tests, the evaluated material properties, such as hardness and reduced modulus, are dependent upon the experimental set-up [10] and are not necessarily indicative of intrinsic properties. In addition, the reduced modulus, which is generally determined from the initial region of the load versus depth curve during unloading, may be affected by creep. Inverse FEA modelling, however, is able to accommodate such phenomena, providing the constitutive model employed is sufficient to capture the actual time-and stress-dependent mechanical response of the material.
The IFEM is predicated on the basis that a single combination of material properties and boundary conditions will result in the observed experimental conditions. Proceeding from this assumption, it is possible to iterate the possible material parameters for the given boundary conditions and the given constitutive relation until the properties that lead to results closest to the experimental observations are obtained. IFEM, therefore, provides the possibility of combining the output of the nanoindentation test with an FE model to inversely obtain the material properties, despite the complex stress state in the material induced by the indenter [15]. In the modelling process, a key step is choosing a material constitutive model that is able to accurately represent the mechanical behaviour of the material. This is essential for the method to establish meaningful material property parameters and a good fit between the FEA modelling and the test. In the proposed method, a least-squares objective function is used, given by s tr e > sh? where S t is the measured deformation depth at the ith time interval, t i ; S F is the depth predicted by the FEA at time t i ; and n is the total number of time intervals. To match the results from the FEA and experimental tests, linear interpolation was used during the minimization process of equation (2.1) to ensure that the time intervals for the FEA and the test data coincided. In this work, Matlab R2013a was used in combination with ABAQUS 6.11-3 in a gradient-based nonlinear least-squares method for minimization of the objective function. When the relative change of F (equation (2.1)) or the change of each searched material property parameters was less than 10 −8 , the minimization process was considered to have converged. Figure 2a shows the flow chart for the minimization process together with the mesh used in the axisymmetric FE model (figure 2b).

Material models
In this work, four material constitutive models were evaluated: (i) a linear VE Prony series model, (ii) a VEVP model, (iii) an NVE model, and (iv) an NVEVP model. For the VE model, the intrinsic ABAQUS linear Prony series method [37] was employed. For the VEVP model, the Prony series was combined with a VP model [38]. For the NVE model, Schapery's equation [39] was combined with a linear Prony series and for the NVEP model this was combined with a VP model. These models are described in the following subsections.

(a) Linear viscoelastic material model
The behaviour of many types of polymer can be described by a Prony series model [37,39], which can be represented by springs and dashpots, as shown in figure 3a. The instantaneous shear modulus is G 0 , the equilibrium modulus is G ∞ and the time-dependent modulus G(t) is given by where g k = G k /G 0 and Σg k < 1. In this research, we take the maximum of k to be three as additional terms yielded no further improvement in fitting to experimental data. g k is the relative relaxation modulus and τ k is the corresponding relaxation time.
where σ and are the deviatoric stress and strain tensors, respectively.

(b) Viscoelastic-viscoplastic material model
The VEVP model used in this work is shown in figure 3b, i.e. a VE model plus a VP component in series. When the VE von Mises stress in the material is lower than the current yield stress of the material, σ h (the initial value of σ h is σ s ), of the VP part, the VP component in figure 3b is rigid; i.e. there is no VP strain contribution to the total strain and the material behaves as a VE material. In this case, using equation (3.1), equation (3.2) can be rewritten as The integration term in (3.3) is the linear Maxwell strain, written as The above equation shows that the Maxwell strain is dependent on time and strain rate. Using (3.4), (3.3) can be written as The above equation shows the linear relationship between strain and stress for a Prony-type VE material, i.e. without VP flow actuated. We refer to the subsequent strain, as the linear VE strain. It can be seen that there is a linear relationship between V E and from equations (3.6) and (3.4). As the linear VE strain is a function of the Maxwell strain, the linear VE strain is also both time and strain rate dependent. The deviatoric stress, from equations (3.3)-(3.6), can be written as The increment format of equation (3.4) has been shown to be [37]: and from (3.6) where k is the increment of the linear Maxwell strain defined by (3.4) [37]. Equations (3.8) and (3.9) provide a linear VE computation algorithm and define a nonlinear relationship between deviatoric strain increment and VE strain increment V E for a given Maxwell strain and time increment, which can be represented by a function Ω for later convenience, i.e. Using (3.7), the deviatoric stress increment is and V E is a time-dependent summable quantity, i.e. V E = Σ V E . Thus, equations (3.8)-(3.11) can be applied for linear VE computation. Computations have shown that analysis results with these equations implemented within an ABAQUS material subroutine are the same as those with the intrinsic Prony series model in ABAQUS. This provides confidence in extending the approach to viscoplasticity, as described below.
When the VE von Mises stress is higher than the current yield stress σ h at time t, the material is in an 'overstress' state [38]; VP flow occurs and a VP flow rate,ṗ, is actuated: (3.12) where σ e is the von Mises stress of the material with VE, σ h is the current yield stress of the material and φ is a viscous flow function. If, at a time t + t, a VP deformation occurs due to the actuatedṗ, an equivalent VP strain increment p can be given by For a von Mises-type material, equation (3.13) can be simplified to and can be written in the form where φ −1 is the inverse of flow function φ. Equation (3.15) illustrates the basic difference between elastoplasticity and VP. In elastoplasticity, yielding occurs on the yield surface, and stress is plastic strain rate independent, therefore φ −1 (ṗ) ≡ 0 and σ e ≡ σ h . In VP, stress is VP strain rate, p, dependent, thus σ e > σ h , with φ −1 (ṗ) > 0. This means that yielding occurs out of the yielding surface for VP. This phenomenon is referred to as 'over stress' [38].
One of the flow functions can be assumed to be [38] where ζ and η are empirical material constants and σ h = r + σ s , σ s is the initial yielding stress and r is the monotonically increasing hardening part of the current yielding stress. Combining (3.15) and (3.16) gives where the first term on the right-hand side of (3.17) is the current yield stress and the second part shows the stress due to VP strain rate sensitivity. When a hyperbolic function is used as the flow function, as shown in (3.16), the parameters ζ and η characterize the amplitude of the stress due to the VP strain rate sensitivity, as seen in (3.17). η has the inverse dimensions of stress; while the dimension of ζ is the same as strain rate. ζ and η were selected as variables in the IFEM minimization procedure. For isotropic hardening, in the case of VEVP behaviour, the yielding stress was assumed to obey a power hardening rule, i.e.
where p is the total equivalent plastic strain at time t. A is an empirical material constant, which we assume is equal to 1. Through experiment it was found that a value of m = 3 led to good fitting for the investigated material. Thus, a constant value of m = 3 was used to avoid an extra variable in the minimization procedure. Strictly, m may be time, strain and strain rate dependent in the VEVP case, especially when strain rate is high; however, in the indentation experiments described in §2b, the strain rate was computed to be low (10 −2 to 10 −3 ). Therefore, the material index parameter is assumed to be strain dependent but strain rate independent in order to simplify parameter fitting in the IFEM. The initial yield stress, σ s , can be considered an intrinsic material parameter and it is expected that this would depend upon the extent of cross-linking in the printed material and therefore may be z-dependent. In the following analysis, an algorithm is developed to calculate the VP yielding the strain tensor p increment and corresponding σ at t + t when σ e > σ h for an isotropic VEVP material with isotropic hardening for a given stress tensor σ t at time t and the given deviatoric strain can be divided into two parts, a VE strain increment V E and a VP strain increment p . From (3.10), the VE strain increment after yielding can be written as and For VP, it has been pointed out that the normality hypothesis still holds [38]; however, the yield point may occur out of the yield surface, as shown by (3.15) or (3.17). Thus, for a von Mises material we still have where p is the equivalent VP strain increment during t, as shown by (3.13 which can be further rewritten as It can be shown that1 is a number close to one, dependent on the time increment t. Multiplying by 2G 0 on both sides of (3.24) and considering (3.22) and (3.21), we can derive and where σ tr is the computed VE stress when VP strain p is ignored in (3.24), which can be referred to as a trial deviatoric stress tensor. Taking the contracted product of (3.26), we have where σ tr e is the equivalent stress of σ tr . Combining (3.26) and (3.29) gives where n is the direction of the VP strain tensor increment p , as shown by (3.22). The significance of (3.30) is that the direction of the unknown VP strain increment p tensor can be determined before p is obtained; it is uniquely determined by the current yield surface, regardless of equivalent VP strain p. Equations (3.24)-(3.30) are similar to the case of elastoplasticity [38]. Now it is possible to obtain a solution for p by combining (3.13), (3.18), (3.29) and (3.30) to form an algebraic nonlinear equation for a von Mises-type material: where σ tr e is computed from equation (3.27) and σ h can be obtained from the given hardening rule or (3.18), as in this case the direction of VP strain increment has been determined by (3.30), and the trial stress tensor σ tr and σ tr e are already known for the given strain increment from (3.27). As both σ e and σ h are functions of p only, from (3.29) and (3.18), (3.31) is a nonlinear equation involving a scalar unknown, p, only. Newton iteration can thus be applied to obtain p. The computed p can then be applied to (3.24) to update VE strain increment V E , and further by (3.11) to update the actual deviatoric stress increment σ . If the computed equivalent trial stress σ tr e is lower than the current yield stress σ h , no VP flow is actuated; thus p = 0-therefore, from (3.29), the trial equivalent stress σ tr e is the actual equivalent stress σ e and the corresponding trial stress increment σ tr e is the actual deviatoric stress increment. Comparing equation (3.31) with the elasto-viscoplasticity result in [38], there are two basic differences: (i) σ tr e should be computed using the VE algorithm; (ii) a t-dependent number1 close to one, as shown by (3.25), is included in the equation. These two differences make the current VP flow time, strain rate and VE deformation dependent. This is the reason that this constitutive relation is referred to as a VEVP constitutive material model. The next step in this work was to write the above algorithm into a material subroutine to enable incorporation in a FEA. The current stress tensor σ , strain tensor , strain increment and time increment t can be provided by the FE main program; one of the main tasks of the material subroutine is to use the proposed algorithm to compute the corresponding stress increment σ based on the information provided by the main code. In order to use the proposed VEVP material model in an FEA in a material subroutine, the stiffness matrix needs to be computed. This will be determined in the next section.

(c) Schapery's nonlinear material constitutive model (nonlinear viscoelastic)
For polymers that are intrinsically nonlinear, an NVE constitutive relation may be more appropriate. A linear Prony series from (3.2) can be written as It is possible to introduce three nonlinear coefficients h e , h 1 and h 2 , leading to a nonlinear Prony series model of Schapery's type [39] with Considering (3.1), (3.33) can be simplified to In order to use the nonlinear material constitutive model in an FEA, (3.34) should be integrated. A nonlinear Maxwell strain at time t n is defined as In the following analysis, the deviatoric strain tensor at time t is expressed as , while at time t + t it is expressed as + . A similar notation is used for other values, hence, from (3.35), Dividing the integration time period 0 ∼ t n+1 into 0 ∼ t n and t n ∼ t n+1 in (3.36), we have and Considering (3.35), (3.38) can be simplified as In (3.39), it is assumed the time increment t is small, so that d(h 2 )/dτ can be replaced by (h 2 )/ t. Hence, integration of (3.39) results in Substituting (3.40) and (3.41) into equation (3.37), we obtain where k = n+1 k − n k is the increment of nonlinear Maxwell strain shown by (3.42). It can be seen that, when h e = h 2 = h 1 = 1, (3.42) simplifies to (3.8), i.e. Prony series linear VE. It can be seen therefore that the integration process for Schapery-type NVE also applies to linear VE, the latter simply being a special case. Similar to the linear VE derivation in (3.9)-(3.11), it can be written for the nonlinear case that and the nonlinear Maxwell strain increment k is computed from (3.42).
Regarding the stiffness matrix, it has been pointed out [38] that an analytically derived stiffness matrix is not always obtainable; however, an accurate stiffness matrix is not absolutely necessary as an approximate one is sufficient in many cases. For the NVE material constitutive model, both an approximate analytical stiffness matrix and accurate numerical stiffness matrix will be illustrated. In order to derive the analytical stiffness matrix, the variations of coefficients increments of h e , h 1 and h 2 will be ignored for simplification. From the general relationship between stress/strain and deviatoric stress/strain, we have the following equations: where δ and represent variation and increment, respectively, K is the material volumetric modulus and I is a second-order unit tensor. Considering the variations of the increments in (3.43) and (3.44), and ignoring variations in h 1 , h 2 h e , we have Inserting (3.48) into (3.47), we arrive at All elements of the material stiffness matrix can now be derived from (3.51). Computations have shown that the stiffness matrix in (3.51) is not only applicable for the NVE case (where at least one of the coefficients of h 1 , h 2 and h e is greater than 1), but also works well for both the VE case (where h 1 = h 2 = h e = 1) and the case of VEVP. In case a more accurate stiffness matrix is needed, a numerical method can be used. Using Voigt notation, let we have where σ i and ε j are the stress and strain components in Voigt notation, shown by (3.52) and (3.53), respectively, γ ij is the strain component in Voigt format, and D ij in (3.54) is the element of the material stiffness matrix. Thus, we can use (3.54) to compute the material stiffness matrix directly using a numerical method. The subroutine for the stress-strain relation is called six times for the computation of the stiffness matrix for each material point; therefore, the potential increase in accuracy comes at the cost of computational efficiency. However, (3.54) is available for the case where an analytical stiffness matrix cannot be obtained.
The final consideration in this material constitutive model is the format of the nonlinear coefficients. Poons & Ahmad [40] showed that, for material with Schapery-type nonlinearity, the nonlinear coefficients in 3D cases can be taken as

Results and discussion (a) Validation of material models and inverse finite-element method
In order to validate the numerical method, an ultra-high molecular weight polyethylene (UHMWPE), which had previously been characterized, was tested. This material is known to show viscoelastoplastic characteristics and the validation uses the VEVP material model. The UHMWPE preparation has been reported in [41], and the method described in §2a was used    to make an indentation sample. A nanoindentation test comprising a 6 × 6 array of indents, with a step length of 100 μm in both x-and y-directions, was conducted on the central region of the sample. The depth-time curves are shown in figure 4a and the average (seen as the bold solid line in the figure) was then used as a reference curve for the inverse FE modelling. Using the proposed VEVP material model and IFEM technique, the properties of the UHMWPE were determined, as shown in table 1. It can be seen in figure 4b that these parameters enable an excellent fit between the experimental load-depth curve and that predicted with the FE mode, thus confirming the applicability of the FEA indentation model, the VEVP material model and the IFEM technique for finding best-fit material parameters. For further validation, the results from a uniaxial tensile test at a constant engineering strain rate of 0.2 s −1 at room temperature given in [41] were predicted. A 3D FE-based model of the tensile test was performed, using the VEVP constitutive model with the material parameters given in table 1. Figure 5 shows a comparison of the experimental results from [41] and the results of the model. This shows excellent agreement, suggesting that the proposed method is able to capture the material behaviour well.  (b) Nanoindentation of three-dimensional printed sample Having validated the proposed method, it was then applied in the determination of the material properties of a 3DP sample. The sample was prepared using the method described in §2a and a 3 × 35 array of indentations was performed to determine the spatially resolved material properties. In each column, the load was held constant, but across the columns the load was varied with 25, 50 and 75 mN maximum loads used, respectively. In the z-direction the depth interval was 50 μm, with the first indentation at a depth of 40 μm from the sample surface. The distance to the side edges was more than 500 μm. The matrix of points is illustrated in figure 6b. The procedure for the indentation was described in §2b. Figure 6a shows the indentation marks left on the test surface during one of the experiments, using a 75 mN indentation load. It can be seen that the residual imprint of the indentation varies with the z-position, suggesting, at least qualitatively, that there are material variations in the z-direction. As the largest residual impression is closest to the top surface, this would indicate that hardness increases with depth below the surface and this is confirmed in figure 6c. The greatest depth is seen with the indentation closest to the top surface (labelled 1), and then the maximum depth is gradually smaller for subsequent indentations (labelled 2, 3, . . . and 6, respectively). After around eight indentations, the penetration depth becomes independent of z. The figure shows the load-depth curves for the first six indentations, each an additional 50 μm in the z-direction. Figure 7 shows the variation of penetration depth, hardness and reduced modulus with increasing z, i.e. distance below the top surface, while keeping x constant. It can be seen that the maximum achieved indentation depth is closest to the top surface (figures 6c and 7a). As the indentation location is increased in the z-direction, the maximum indentation depth decreases, appearing to reach an asymptotic value at around 500 μm. The hardness and reduced modulus show similar behaviour. It is noted also that the hardness is indentation load dependent, as discussed in [42], while the reduced modulus is relatively load independent if the indentation load is at least 50 mN. These results correlate well with the qualitative observation of varying residual imprints, indicating that the printed VeroClear material properties are z-dependent, and, in particular, depend on the number of printed layers of material. Because the mechanical property of this reactive inkjetted material should only depend on the extent of the cross-linking during post-build UV illumination, it is supposed that the primary cause of this phenomenon is the varying number of exposures to UV, i.e. the bottom-most layer will have the same number of exposures as there are layers of material and is more likely to have saturated the cross-linking   to the surface where the properties are changing rapidly, and within the bulk where the material properties are largely independent of z. A maximum load of 50 mN was used. The maximum depth versus z-position at this load is shown in figure 8a with the selected 'surface' and 'interior' positions labelled 'A' and 'D', respectively. The depth-time curves at these positions are shown in figure 8b.
From the nanoindentation results, an assessment of the applicability of each of the proposed constitutive models was made, by taking the depth versus time plots and using them as an input to the objective function (equation (2.1)). The IFEM was then employed for each to determine the optimum parameters of the constitutive models that would best reproduce the experimental results. The comparisons between experiment and numerically obtained depth versus time plots are shown in figure 9 for each constitutive model.
It is can be seen in figure 9a that, for the VE material model, both the loading stage and load-holding stage can be fitted well; however, the FE predicted elastic recovery during the  unloading stage is much higher than that seen experimentally. The introduction of nonlinear coefficients by using the NVE model improved this, as seen in figure 9b. Mechanistically, the aim of introducing the nonlinear coefficients, h 1 , h 2 and h e , in equation (3.34) is to increase the viscous strain component in the total strain for the given loading conditions, thus to increase the viscous deformation which is unrecoverable during the unloading stage, leading to better fitting. However, computational experiments showed that the nonlinear coefficients could not be made sufficiently high to enable the FE predicted elastic recovery fit exactly to the experimental curve. This is because too high nonlinear coefficients result in failure of the modelling procedure due to convergence problems. The optimal results using the NVE material model, therefore, lie between the VE and VEVP material models. Owing to the aforementioned limitation of the nonlinear material model, the VEVP and NVEVP material models were thus introduced, with the aim that the over-predicted FE elastic recovery can be suppressed by VP deformation during the loading and holding stages. Figure 9c,d illustrates that both these models are significantly more capable of modelling the unloading curve than the VE models. The relative merits of these models are difficult to assess just by visual assessment of the curve fits in figure 9; therefore, a quantitative method of comparison was developed.
To quantify the ability of a material model to capture the behaviour during indentation a relative residual R r was computed. This is given by   where S Fk and S Tk are the depth at time t k by FE modelling and experimentation, respectively, n is the division number during the test (or FE modelling). R r represents the relative difference of indentation depth between test and FE modelling, i.e. the closer R r is to zero, the closer the agreement between the FE model and the experiment. The results are given in table 2, which shows that the NVEVP model fits the experiment the closest, suggesting that this constitutive model is the most appropriate to characterize the mechanical response of the VeroClear to complex loading. In order to determine the root cause of the variance in fits seen in table 2, the predicted stress distributions at attaining maximum load in the FE model for the various material models were compared (figure 10). The differences in stress distribution are striking. When using a linear Prony series, a high stress field is located immediately below the indenter (figure 10a). With the VEVP model, the high stress field is near to the contact edge area, where VP flow occurs. Subsequently, this high stress region moves with the edge when the indentation load increases (figure 10b) and during the load hold stage. Similarly, for the NVEVP model, the high stress field is located near the contact edge, and moves towards and with the contact edge when the indentation load increases  the fit of experimental and predicted load-depth plots should be expected. The total model image of the FE model used for figure 10a is shown in figure 10d. As can be seen, the stress-influenced area during indentation, compared with the total model, is rather small; the stress-affected area is far away from any model edge except the top surface. This suggests that the model sizes shown in figure 1b are sufficiently large to be unaffected by the geometry and boundary condition simplifications used in the model. The evaluated material parameters using the various material constitutive models for indentation locations A and D are shown in table 3, where E in is instantaneous modulus, τ 1 , τ 2 and τ 3 are relaxation times and g 1 , g 2 and g 3 are the corresponding relative relaxation moduli (see equation (3.3)), σ s is the initial hardening yielding stress for least depth residual (see equation (3.18)), β e , β 1 and β 2 are the nonlinear coefficients of the NVE and NVEVP material constitutive models (see equations (3.34) and (3.55)) and ζ and η are VP parameters, see (3.16) and (3.17).
The results in table 3 show that the highest instantaneous modulus is predicted by the VE material model while, for the VEVP and NVEVP material models, similar lower moduli are predicted. With a small VP strain rate,ṗ, in (3.17) for the given flow function as shown by (3.16), the effects of both ζ and η to the stress difference of σ e − σ h = 1/η sin h −1 (ṗ/ζ ) ≈ṗ/ηζ are nearly equal, although ζ and η have different physical meaning and dimensions. In order to have higher computational efficiency, it was assumed that ζ and η have the same mathematical value with different dimensions, thus a variable can be reduced in the inverse FE modelling (note that this is acceptable for the case of small deformation and low strain rate but cannot be generally assumed). This can be seen in table 3, which also shows that VP parameters in the NVEVP model are higher than in the VEVP model, especially at location D. This is because the use of nonlinear coefficients in the NVEVP model results in a higher viscous deformation, hence, less VP deformation is required to fit to the experimental data. As a result, the corresponding initial yield stresses are also higher for the NVEVP model. However, generally the indentation depth is not sensitive to these VP parameters for the Veroclear material. This is the why discrete values are obtained in the table.
Although the instantaneous moduli obtained for the VEVP and NVEVP models are nearly the same, as shown in table 3, the relative relaxation moduli and the relaxation times are somewhat different. This is because of the addition of the nonlinear coefficients in the NVEVP model leading to higher viscous deformation. The computations have shown that for normal simulations, the coefficient of h 1 (see (3.33)) is restricted to nearly h 1 ≈ 1 or β 1 ≈ 0; β 1 should be less than the order of 10 −4 to maintain model convergence and best fitting. In order to simplify the inverse FE modelling, it was assumed that β 1 = 0 for best fitting. This is also shown in the table.
The results show that, regardless of the constitutive model, the instantaneous elastic modulus is predicted to within 20% and is the dominant contribution to the material behaviour. Each model captures a different aspect of the deformation, with detailed differences at the level of the individual parameters. However, the parameters are largely of the same order of magnitude, suggesting that the primary physical mechanisms are being captured by each model in a similar way, particularly when the VE and VEVP models are compared. The differences in the models are larger at greater z-depths, suggesting that the material that has suffered more UV illumination develops a greater material nonlinearity, potentially reflecting a change in micro-and molecular structure. Within each model, it can be seen that the evaluated instantaneous modulus within the body of the sample is approximately 50% greater than that at the top surface. This evidence once again suggests that 3D inkjet printing can result in significant variation in material properties throughout a part, immediately post manufacture. If controlled, this could lead to the control of part properties through a volume and raises the potential for an extra design freedom for engineering component manufacture not available to traditional manufacturing techniques. It should be noted, however, that the properties of the material, and their spatial variation, may change with time as the material reacts to environmental conditions. The method described in this paper will be used to investigate this in future work.

Conclusion
An inverse analysis technique has been developed that combines nanoindentation with FE modelling. The primary purpose for developing this technique was to enable measurement of material property variations in 3D printing components, and a range of potentially relevant material models were developed and tested. Two characteristic locations in a 3DP sample were examined using the proposed method, near to the surface, where nanoindentation results suggest large spatial variations exist, and within the bulk, where the properties appear to be independent of location. In both cases, the method was able to determine the material properties with small residual errors between the model predictions and the nanoindentation-measured depth-time curves.
It can be concluded, therefore, that the proposed method can accurately determine spatially resolved properties as long as the constitutive material model used in the analysis is capable of fully representing the material behaviour in the test. It can be further concluded that the method can be used to investigate manufacturing processes which induce spatial variation of properties and that the accurate, spatially resolved material parameters thus derived can also be used in models to predict the performance of the manufactured part in-service.
Data accessibility. All data are stored centrally at the University of Nottingham. Contact Ian Ashcroft at ian.ashcroft@nottingham.ac.uk for access details.
Authors' contributions. X.C.: researcher and drafter of paper. I.A.A.: project lead, revision, submission and final approval of paper. R.D.W.: project supervisor, paper revision. C.J.T.: project supervisor, paper revision.
Competing interests. We have no competing interests. Funding. This work was funded through the EPSRC Centre for Innovative Manufacturing in Additive Manufacturing.