Characterising the frequency‐response of ultra‐soft polymers with the Virtual Fields Method

This paper presents a novel apparatus, combined with application of the Virtual Fields Method (VFM) for the high frequency characterisation of the mechanical response of ultra‐soft materials. The viscoelastic response is characterised under harmonic deviatoric loading at a range of frequencies using the dynamic VFM in the frequency domain. Obtaining useful high rate data on soft materials is challenging using typical test methodologies, as the low speed of sound in the material makes stress equilibrium difficult to obtain; additionally, the low stiffness results in small stresses and weak measurement signals. This novel apparatus and analysis method have been shown capable of obtaining important material characterisation data that would be impossible to access using currently existing test techniques.


| INTRODUCTION
Soft materials may be defined as materials that have low mechanical stiffness and strength and include polymers, foams, and biological tissues. [1] These materials are widely used in applications in which large amplitude dynamic or impact deformations are experienced; however, there is a paucity of data on their properties at high strain rates and frequencies when compared to metals and ceramics. One reason for this is that they are challenging to characterise under these loading conditions. [1][2][3][4] Typical engineering test methodologies rely on the inertial stresses in the specimen being negligible in comparison to the strain induced stresses; [5,6] however, the low modulus to density ratio of soft materials makes this difficult to achieve in high rate characterisation experiments, because the speed of sound in the material is low. Additionally, owing to the low stiffness and strength, stresses are small 1 ; this causes difficulty in measurement of the stress state using conventional stress measurement techniques.
The above limitations can be partially mitigated via adaptation of standard test methodologies. Approaches such as Dynamic Mechanical Analysis (DMA) and Rheometry can adequately solve the problem with low forces but struggle to obtain data at frequencies above 50-100 Hz. At higher rates, the use of semi-conductor strain gauges and very thin 1 On the order of Pa. specimens on a Split Hopkinson Pressure Bar (SHPB) partially mitigates problems with weak transmitted signals and specimen stress equilibrium, respectively. [1] These adaptations can lead to other issues, however, such as the exacerbation of radial inertial effects.
The behaviour of polymers outside the envelope of testable strain rate regimes can be predicted via Time Temperature Superposition (TTS), [7,8] where changes in rate can be related to changes in temperature. This allows for the high rate behaviour of a material to be predicted from low rate low temperature data; [9,10] however, there are some significant limitations to this method. TTS is mainly validated for small amplitude loading in which the polymer response remains linear viscoelastic: it is only recently that this approach has been shown to work for larger strains. [8] Furthermore, some polymers (such as the one tested in this paper) exhibit low sensitivity to temperature, and some materials (such as biological tissue) degrade when taken outside of a narrow range of temperatures. [2,11,12] This means that TTS is of limited utility for these materials.
A further limitation of standard engineering test methodologies is the assumption of homogeneity: detection of material heterogeneity or specimen strain localisation, for example, requires a form of full-field measurement. [13] Additionally, if 1D tests are used, multiple tests are required in order to detect material anisotropy. [14,15] As a result, interspecimen variability can cause errors when quantifying effects of of inter-test changes such as specimen orientation or location. In order to overcome these limitations, the use of full-field measurements allows for the direct measurement of heterogeneity, [15] and the use of the Virtual Fields Method (VFM, further described below) can allow for multiple model parameters to be extracted from a single test. [16] Additionally, VFM analysis can be designed to make fewer assumptions about specimen stress equilibrium or deformation state, allowing for greater flexibility in the choice of test boundary conditions.
Ultra-soft materials such as biological tissue often have a shear modulus (μ) that is orders of magnitude lower than the bulk modulus (K): e.g., for brain tissue μ ≈ 2 Â 10 3 Pa and K ≈ 2 Â 10 9 Pa. This means that low rate measurements of the Young's modulus will effectively measure the shear modulus, 2 and high rate measurements of the Young's modulus will measure the shear modulus plus an error due to inertial effects. As a result, tests to establish Young's modulus are not particularly valuable: testing should instead focus on separating the deviatoric and the dilatational response, particularly if high rate characterisation is required. Additionally, the large difference in deviatoric and dilatational moduli results in a large difference in the related wave speeds; as a result, the specimen may be in longitudinal (static) stress equilibrium fairly early in the test but may not reach shear stress equilibrium before the test is over. Hence, when characterising the shear response of soft materials, it is important that the dilatational loading is accounted for.
In order to adequately characterise the frequency-response of an ultra-soft material at high rates, a novel experimental apparatus has been designed and built. Stress and strain measurements are non-contact, in order to allow for the testing of ultra-soft materials, and full-field, in order to allow for spatial discretisation of heterogeneous material properties. An important component of the designed experimental apparatus is the use of the VFM, described in the next section.

| The virtual fields method
High fidelity full-field measurement methods such as DIC, Moiré interferometry, and grid methods provide rich experimental data; however, standard engineering test methodologies such as the uniaxial tensile test fail to fully exploit these. [17] In contrast, the VFM allows one to take advantage of this rich dataset, enabling the extraction of multiple material parameters (and their spatial distribution) from a single experiment. This allows one to obtain material parameters from far fewer tests than traditionally required-a particular advantage when specimens are difficult or expensive to obtain.
Furthermore, traditional engineering test methods that rely on measurements of forces or displacements at a small number of discrete locations assume a given stress distribution and give results of reduced quality when these assumptions are not correct. The VFM allows for a heterogeneous stress state; indeed, heterogeneous stress states are desirable, 3 as they allow for the extraction of multiple material parameters from a single experiment. [14] The VFM is based on the Principle of Virtual Work (PVW), which can be manipulated to remove the effect of traction or 3 It should be noted however that VFM using DIC data typically needs to make an assumption about the through thickness stress. A heterogeneous through thickness stress distribution is not desirable. acceleration terms-something that is particularly valuable for dynamic characterisation. In practice, removal of the traction term allows one to use the specimen inertia as a virtual load cell: rather than trying to obtain static stress equilibrium, the dynamic acceleration field is used to the advantage of the researcher. [18] 1.1.1 | Worked example In order to demonstrate the use of the VFM in dynamic simulations, a simple worked example is helpful. Consider a homogeneous linear elastic cylindrical rod, with a density ρ ¼ 1000 kg/m 3 , shear modulus μ ¼ 4000 Pa, shear wavespeed c s ¼ ffiffiffiffiffiffiffi ffi μ=ρ p ¼ 2 m/s, and length L ¼ 10 mm. Both ends are fixed in place using rigid clamps. Under the chosen loading, the bottom clamp is oscillated in rotation at 500 Hz for 2 ms such that the displacement amplitude at the surface of the specimen is A. The surface rotational displacement of the loaded end is given by uðtÞ ¼ Aðcosð1000πtÞ À 1Þ until t ¼ 2 ms, after which the loaded end is stationary. After 3.5 ms, the distribution of surface displacements along the specimen length is as shown in Figure 1.
Considering a thin, quasi-1D strip at the surface of the sample, if the upper end were unclamped, shear stresses along the length could simply be calculated via direct integration of the dynamic equilibrium equation 4 from the free end of the sample. However, if the end is not free and the distribution of traction forces is unknown (or if the displacement field is not known for the entire domain), the full VFM is required.
The virtual work equation is obtained by multiplying the dynamic equilibrium equation through by a compatible virtual displacement and strain to give [14] À ð V σ : ϵ * dV |fflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflffl ffl} where the terms are as follows: In quasi-static experiments, the inertial contributions to the stress field (i.e., acceleration virtual work in Equation 1) are negligible, and this term can be neglected-material characterisation is instead done by accurately measuring the surface tractions. In contrast, acceleration terms are often significant in dynamic experiments, and surface tractions can be difficult to measure: appropriate selection of the virtual field can eliminate the traction terms as follows.
First, a virtual displacement field, u * , is chosen such that only one component (ϵ * XY ) of the compatible virtual strain tensor is non-zero. This, combined with the quasi-1D nature of the deformation field, 5 allows one to analyse the problem in only one dimension. In 1D linear-elastic shear, with σ ¼ μϵ, Equation 1 becomes where X 0 and X 1 are the boundaries of the domain and μ is the shear modulus. Careful selection of the virtual displacement allows the traction force term to be eliminated from Equation 2. The virtual field used by Moulart et al. [13] and Yoon et al. [19] can be adapted for shear: where u * X and u * Y are the virtual displacements in the X and Y directions, respectively. This eliminates the traction force term, resulting in Substituting the displacement fields described in Figure 1 and the virtual displacement fields described in Equation 3 into Equation 4 at time t ¼ 3:5 ms results in giving a predicted shear modulus of 4 kPa-this matches the given material parameters. Figure 2 shows the shape of the displacement, strain, acceleration, and virtual work fields at time t ¼ 3:5 ms; plots with numerical values and more detail are shown in Appendix B.

| Analysis in the frequency domain
As the purpose of this experiment is to investigate the frequency response of soft materials, it is logical to examine the data in the frequency domain rather than the time domain. To this end, a Fast Fourier Transform (FFT) is used to obtain the strain and acceleration for each location on the observed specimen surface in the frequency domain. This allows Equation 4 to be solved in the frequency domain, directly obtaining the visco-elastic frequency dependant complex modulus for each frequency in terms of the frequency domain acceleration and strain components: F I G U R E 2 Plot of the displacement field and derived temporal (acceleration, VFM inertial term) and spatial (strain, VFM strain gradient term) fields for the problem described in Figure 1. Full plots with numerical values, as well as equations describing the various fields, are shown in Appendix B Here, both strain and acceleration are now complex valued and selected for a single frequency only, giving the complex shear modulus for that frequency f. The apparatus is designed to excite a 'pure' frequency, making almost all of the signal reside on a single frequency band. As a result, the noise-floor can be raised, allowing other frequencies to be readily discarded as noise rather than weak signals.
The loading device (shown in Figure 3) consists of a torsional pendulum with a tunable resonant frequency, driven by a pair of programmable smart shakers. The spring clamp can be moved in order to change the effective length of the torsion spring, allowing one to tune the resonant frequency of the device as shown in Figure 4. The upper clamp can be rotated, allowing for the specimen to be preloaded if needed. This clamp provides a static displacement boundary F I G U R E 3 Experimental apparatus used for torsional loading of samples F I G U R E 4 Range of obtainable resonant frequencies using the high frequency (light) or low frequency (heavy) pendulum head condition 6 with a wide range of allowable (quasi-static) motion. The lower clamp provides a displacement boundary condition with a limited range of motion, but a wide range of possible frequencies. The device was designed to use two interchangeable masses on the torsional pendulum, in order to provide for a larger range of resonant frequencies than would be obtainable via the changing of the spring length alone; however, the experiments conducted for this paper made use of the low frequency head only. For the current setup, the maximum oscillatory displacement is limited by the performance of the Smart Shaker pair; if greater displacements are required, the location of the smart shaker pair can be changed.
If the experiment does not aim to probe a large range of macroscopic strain states, then the upper clamp can be omitted. This omission provides a stress-free boundary condition at the upper surface of the specimen, obviating the need for the full VFM-stresses can instead be obtained via simple integration of the acceleration field.

| Data acquisition
In order to obtain DIC images, a Photron FASTCAM SA-X2 1080K-M4 camera is used, running at 25 000 fps using a lens of 200 mm focal length. Light is provided by four custom built LED lights, and glare is avoided via the use of crosspolarisation. To achieve this, vertically oriented linear polarisers were placed on the lights, and a horizontally oriented linear polariser film was mounted on the lens. The light arriving on the specimen is vertically polarised: specular reflection (glare) does not change the polarisation of the light, allowing the horizontal polariser filter on the camera to remove glare (and hence over-exposure of the image). In contrast, diffuse reflection (the actual image) randomises the polarisation of the light, allowing it to pass through the horizontal polariser filter on the camera. For further detail on this technique, the reader is referred to LePage et al. [20] The downside of this method is that significant amounts of light are inevitably lost in the double polariser filters, requiring the use of brighter lights than would otherwise be required. An alternative method for avoiding specular reflection is to instead ensure that the surface of the specimen has a matte finish (via the application of talc powder for instance); however, this is not always possible, particularly when imaging biological samples.

| Specimen preparation
Samples consisted of an ultra-soft platinum cure PolyDimethylSiloxane (PDMS, silicone) elastomer sourced from Benam Chemicals, England. Both EcoFlex GEL and EcoFlex 10 were used, exhibiting nominal Shore hardnesses of 000-35 and 00-10, respectively. The Ecoflex gel was mixed with 10% thinner by weight, in order to increase pot-life and further decrease the shear modulus. Both polymers are normally transparent; however, 0.5 wt% black pigment was added to both, in order to allow adequate imaging quality for DIC. The casting procedure was the same for both: vacuum degassing at approximately 5 kPa (absolute) followed by curing at 20 C for 24 h. The silicone exhibited no observable volume change while curing (less than 0.1%, according to the manufacturer's datasheet [21,22] ). Samples all had a diameter of 30 mm and a height of approximately 70 mm. The region of interest was near the bottom of the specimen; however, the extra length above the region of interest reduced the magnitude of waves reflected from the top of the specimen. This extra length was not needed and can be neglected; however, it makes the experiments slightly easier to run.
Specimens for rheometer testing were cast in flat sheets of 3 mm thickness from which circles were cut with a punch once curing was complete. Specimens for use in the new apparatus were cast in smooth cylindrical moulds and extracted once curing was complete. The speckle pattern was then applied with an airbrush, using a silicone-based paint. This paint consisted of approximately 50 wt% white silicone pigment, 20 wt% silicone thinner and 30 wt% of the same silicone used for the samples.

| Experimental procedure
The frequency response of the material was probed at a variety of frequencies. First, the specimen was fixed in the clamps, and reference images of the (static) specimen were taken. After each change in loading state (oscillation frequency or oscillation magnitude), this state was held for at least 60 s before any data were recorded, in order to allow for high frequency relaxations to occur. In the experiments reported here, the static (DC) strain state of the specimen was not varied; however, the apparatus is capable of doing so. It is intended to include this change in future experiments in order to evaluate the hyper-viscoelastic frequency response of specimens.

| Data extraction
MatchID DIC software was used to convert video to spatial data. The strain, displacement, and acceleration data were extracted directly from the DIC software. Whilst the domain was idealised as 1D; in reality, the imaging domain was 1024 pixels long and 384 pixels wide. The displacement along the central 150 pixels in the width direction was averaged in order to obtain the central displacement-this had the effect of significantly reducing noise without affecting signal quality. Table 1 shows the imaging parameters as well as the processing parameters used when calculating the deformation and the spatial and temporal derivatives (strain and acceleration). After extraction from the DIC software, the strain and acceleration data were transformed from the temporal domain to the frequency domain via an FFT, resulting in complex valued data at a range of frequencies for each spatial point.

| Simulations
Prior to commencing experiments, the proposed experimental setup was simulated using Abaqus © FEA software. The nodal displacements from this simulation were then fed through the analysis chain in order to validate the feasibility of the proposed experimental setup. Once this step was validated, the full experimental setup and analysis chain (including DIC) were simulated as shown in Figure 5. A photo of a specimen was taken using the same camera as used in the experiments, and the camera optics were calibrated. Calibration is performed by taking images of calibration plates with known dimensions in order to generate a mapping between image co-ordinates in pixels and real-world spatial co-ordinates. The nodal deformations from the FEA simulation were then used in conjunction with the camera calibration data to deform the speckled image. Note that this simulated image deformation took full 3D optical effects into account. These artificially deformed images were then analysed using DIC in order to extract displacement data, and this displacement data were then fed into the analysis chain. Artificial deformation of images was conducted using MatchID's FE-Def module, as was DIC reconstruction of the deformation fields. All cameras have some inherent noise associated with the images-the camera used for these experiments had a grey-scale noise level of approximately 0.45% of full scale. The effect of noise on the obtained material parameters will be explored in Section 4.

| Simulation results
The simulation performed was of an experiment in which the frequency of oscillation was 1 kHz, and the surface displacement amplitude 0.033 mm. The specimen was allowed to reach steady state before the deformation state was recorded, in order to emulate the proposed experimental protocol. The frequency-dependent storage (19.6 kPa) and loss (4 kPa) moduli at 1 kHz were directly input to Abaqus, in order to allow for direct comparison with the material properties calculated by the VFM. Each element was 150 μm in size, approximately equivalent to 5 pixels. A Lagrangian formulation using C3D8R elements and a linear visco-elastic material model were used, and Abaqus's Explicit Dynamic solver was used for the dynamic simulations. In this simulated experiment, data were recorded for only 1 ms (1 loading cycle). This recorded FEA displacement was then used to deform images of the experimental setup via MatchID's FE-DEF module, before performing DIC on the output images. The deformations computed by the DIC were then used to characterise the material properties using Equations 5 and 4 and the virtual field described in

| Low frequency tests on silicone rubber
Two materials were tested in this set of experiments, as previously described: the stiffer polymer was used in order to validate the new apparatus against DMA rheometry data. This stiffer polymer (EcoFlex 10) was used during the commissioning phase in order to avoid some of the difficulties inherent in testing extremely soft materials whilst the device was still being commissioned. Note that EcoFlex 10 is still extremely soft: it is only stiff in comparison to EcoFlex GEL or brain tissue. The deviatoric complex modulus of EcoFlex 10 was characterised on an Anton Paar Physica MCR 301 rheometer at a range of temperatures and frequencies. Specimens were 25 mm in diameter and 3 mm in thickness, punched from a sheet of material, which was cast as described in Section 2.1.3. Figure 6 shows low frequency storage and loss modulus measurements for the EcoFlex 10 material at a range of frequencies obtained on the rheometer, as well as the storage and loss moduli measured by the new apparatus for the same material. Both measurement methods give similar results, validating the new apparatus. Note that 20 Hz is near the upper limit of testable frequencies on the rheometer and near the lower limit of testable frequencies on the new apparatus.

| High-frequency VFM
After commissioning and validating the apparatus using EcoFlex 10, experiments using EcoFlex GEL were conducted with confidence. Whilst EcoFlex GEL is still stiffer than brain tissue, it is the softest readily available silicone elastomer.

| Worked example
Once the images have been taken and processed using DIC software, the full-field deformation and strain fields are extracted. Note that the strain data calculated by the DIC software are inherently smoothed via the choice of Strain Window (SW) size; this is discussed in Appendix A. If material characterisation calculations are performed in the temporal domain, the full-field acceleration field is also extracted from the DIC software. Figure 7 shows an example of the displacement, velocity, and acceleration fields for a single time step in an experiment. The time window used for calculation of temporal derivatives in these images is only three frames, resulting in noisy temporal derivatives (see Appendix A for further detail). This would be a problem if processing was conducted in the time domain; however, all temporal calculations are performed in the frequency domain, effectively filtering the noise from the acceleration field.
F I G U R E 6 Low frequency storage and loss modulus data from standard rheometry tests on phantom polymer at 25 C as well as via VFM on the new apparatus at 20 Hz Once the full-field deformation data 7 have been extracted from the DIC, further processing is performed. As the Region of Interest (ROI) is approximately 30 datapoints wide, spatial averaging of horizontal datapoints is possible; this significantly reduces noise and is done for all data extracted from the DIC process. Care is taken to keep the ROI narrowly centred around the centre of the cylinder, in order to avoid errors due to surfaces that are not perpendicular to the camera. An example of this is shown in Figure 8, where the grey dots are the raw (noisy, 2D) data-field extracted directly from the DIC, and the red line is the 1D averaged value. 7 Displacement, strain and acceleration.
F I G U R E 8 Horizontal displacement vs. co-ordinate. Full field (2D) datapoints are shown in grey, with 1D (horizontally averaged) values shown in red. Note how the horizontal spatial averaging greatly decreases the effect of noise, as multiple samples are taken for each Xposition. See Figure 1 for indication of the meaning of the X-axis values F I G U R E 7 Horizontal displacement, velocity, and acceleration using a time window of three time steps. Units are pixels, pixels/s, and pixels/ms 2 , respectively If material characterisation calculations are being performed in the temporal domain, no further data processing is required before calculation of Equation 4. If material characterisation calculations are to be done in the frequency domain (as was the case for the set of experiments shown in this paper), the displacement and strain data are passed through a temporal FFT in order to obtain the spatially varying frequency components of displacement and strain: The frequency components of acceleration are obtained directly from the frequency components of displacement 8 using e € uðX,f Þ ¼ Àð2πf Þ 2ũ ðX, f Þ: ð8Þ Figure 9 shows an example of the displacement in the time domain and the acceleration in both the time and frequency domains for a single point. The time domain displacement is relatively smooth and appears to have minimal noise. The time domain acceleration is very noisy, however, as the double derivative has the effect of amplifying noise. In spite of this, the signal is isolated to only a small number of frequencies; therefore, in the frequency domain (Figure 9c), it is straightforward to isolate the signal from the noise. Note the presence of harmonics at double and triple the primary frequency, allowing for the extraction of moduli at three frequencies from a single experiment (fourth and higher harmonics exist but have a low magnitude approaching the noise floor, preventing the extraction of useful information). The presence of these harmonics was not expected (as the loading was nominally a pure sine wave); however, they are useful.
Spatially varying complex moduli are now calculated using the frequency domain acceleration and strain data and Equations 5 and 3. Unlike typical test methods, the extracted complex moduli do not assume a homogeneous material, allowing for spatially varying results. The material properties are calculated for each datapoint by integrating over a narrow virtual field of interest centred about that point. For a point at position X i , the virtual field is described by Equation 3, with the limits X 0 and X 1 given by  This method effectively provides a moving average of the material properties, where increasing the domain halfwidth (w) increases the degree of smoothing, reducing noise at the expense of spatial resolution. For the experiments described in this paper, a domain half-width of 100 pixels was selected. Figure 10 shows examples of the spatially varying complex modulus extracted from this experiment at the nominal loading frequency of 250 Hz and the two harmonics at 500 and 750 Hz. Whilst the material tested in this set of experiments is fairly homogeneous (obviating the need for spatially varying analysis), this methodology is designed to eventually be used on soft biological tissue such as brain tissue (which is heterogeneous on both a cellular and a tissue level scale).

| Results
As the material tested in this set of experiments is fairly homogeneous (as shown in Figures 10 and 11), the mean complex modulus for a series of further experiments performed at different frequencies is plotted in Figure 12. While this apparatus was designed to operate above 100 Hz, some experiments were conducted at frequencies as low as 10 Hz in order to probe the lower limits of the apparatus's operating envelope. 9 Note that Figure 12 includes repeated tests at Hz are more reliable/repeatable. This is not surprising, as the acceleration field is much lower for low frequency experiments, resulting in weaker signals and a greater sensitivity to noise. If this apparatus is to be used at frequencies below 30 Hz, it is recommended that the surface traction forces are recorded and taken into account during the VFM calculation, in order to increase the magnitude of the signal. This would be done by modifying the selected virtual field to be non-zero on the boundaries, making the traction term in Equation 2 non-zero. For further detail on how this may work, the reader is referred to Pierron and Grédiac [14] and Yoon and Siviour. [23,24] 4 | EFFECT OF CAMERA NOISE ON MEASURED MECHANICAL PROPERTIES Section 2.4 showed through simulations that the new apparatus could give correct results but did not cover the effect of noise on the accuracy of the measured mechanical properties. All cameras have some noise, however, and this affects the accuracy of the DIC measurements. Therefore, simulations were conducted in which artificial grey-scale noise ranging from 0 to 8% was added to the deformed images. This allows for a quantification of the effects of camera noise on the measured mechanical properties. Unlike the set of experiments conducted in this paper, the specimens simulated in this section were set up with a free end. This allowed for calculation of the stress via direct integration, as the stress at the free end is known to be zero. This was done in order to minimise the effect of choice of virtual field on the outcome of the numerical analysis and to allow for a direct visualisation of the calculated stress field; whilst the full VFM calculates the material properties, it does not directly output the stress field. In order to achieve this, the dynamic equilibrium equation is integrated from the free end, X ¼ 0: where the acceleration is either calculated using a finite difference scheme with a window width of three time steps, or from the FFT as indicated in Equation 8. This calculated stress field can now be directly compared to the strain field measured by the DIC. Direct integration was not used for the experiments as it would prevent the application of prestrains or require that the traction at the ends was measured; this would reduce the flexibility and robustness of the experimental methodology. F I G U R E 1 2 Spatially averaged frequency dependent complex moduli extracted from experiments Figure 13 shows the measured horizontal displacement field for various levels of artificial camera noise ranging from noiseless on the left to 8% noise on the right; note that the increment of noise between images is not constant. Whilst this image displays the effects of increasing pixel noise, a similar effect will occur as the result of decreasing speckle pattern quality. It is not expected that camera pixel noise levels will ever reach 8% (the images taken for the experiments in this paper had a noise level of approximately 0.45% for instance); however, poor speckle quality or lighting could nevertheless result in DIC with significant errors. Whilst the displacement field itself is relatively smooth up until camera noise levels of 1-2%, the process of spatial or temporal differentiation amplifies the noise in the strain and acceleration fields, respectively. Further detail around optimisation of DIC processing parameters for spatial and temporal differentiation is given in Appendix A. Figure 14 shows an example of the calculated stress and strain distribution at a single time-step for varying levels of camera noise. When using a Finite Difference (FD) scheme to calculate temporal derivatives (and hence the stress distribution, via Equation 10) as shown in Figure 14a, it is clear that camera noise levels of 4% and above are too high to allow for meaningful material characterisation. In contrast, when using a Fourier Transform (FT) to calculate accelerations (and hence stress distributions) as shown in Figure 14b, camera noise has significantly lower effects on the quality of the calculated stress distribution.
When calculating complex moduli from time domain data, not only is the calculation highly sensitive to noise (as shown in Figure 15) it is also not simple to calculate the storage and loss moduli; only the magnitude of the complex modulus is easily obtained. Figure 16 shows the frequency dependant shear moduli calculated from the frequency domain data. Note that the calculation is remarkably resistant to the effects of camera noise-as the noise is effectively random for each frame, the FFT easily removes the noise while leaving the signal of interest. Whilst Figure 14b is a significant improvement over the stress field calculated by a Finite Difference Scheme, it still contains spurious high frequency components. In contrast, Figure 16 utilises only the frequency of interest and as a result is unaffected by any noise components at other frequencies. Additionally, the phase angle of the complex modulus is easily calculated when using the frequency domain method, unlike the time domain method.
These simulations showed that the new apparatus and analysis method are robust, particularly when calculations are conducted in the frequency domain. Furthermore, analysis shows that analysing the data in the frequency domain (rather than the time domain, as is common in the literature) gives significantly more accurate results, even when high levels of noise are added to the data. Traditional time domain methods optimise the virtual field to minimise the effect of noise (something that was not done here); however, there is a limit to how much said method can minimise the effects of noise, and it is less robust than the frequency domain method outlined here.
F I G U R E 1 3 Calculated horizontal displacement field for various levels of artificial image noise. Scale is in pixels. All images were produced from the same underlying deformation field 5 | DISCUSSION

| Frequency domain analysis
Unlike most previous experiments in the literature, this set of experimental data contains multiple oscillations of each frequency, allowing for the use of a Fourier transform without being significantly affected by the Gabor limit. As a F I G U R E 1 4 Stress calculated via a finite difference scheme (a), stress calculated via a Fourier transfer and strain (c) vs. spatial position for varying levels of camera noise. All images refer to data at 0.6 ms result, the measured frequency data are sharply band limited, in contrast to data obtained from traditional test methodologies. This method of analysis has a significant advantage, particularly when the data are noisy. During simulated experiments, the time domain approach started to give erroneous results once the camera noise level approached 2% and above, whereas the frequency domain approach was capable of giving accurate results even when the camera noise level was 8%. This will be particularly advantageous if it allows for the use of cameras which exploit on-chip storage to allow speeds of a few millions of frames per second. Such cameras have high resolutions at very high frame rates but can have substantially higher levels of noise. [25]

| Spatially varying mechanical characterisation
Whilst this apparatus and analysis method allows for heterogeneity in the tested materials, this property was not fully exploited for these tests as the material was close to homogeneous. This property is still valuable, however, as this apparatus will be used on soft biological tissue in the future. Ongoing tests are being conducted on soft polymers deliberately manufactured to have spatially varying properties, in order to quantify the ability of this methodology to accurately characterise visco-elastic heterogeneous materials. F I G U R E 1 5 Measured frequency dependant complex vs. image noise level using time domain analysis. Note the significant overestimation of the modulus even for noise levels as low as 2% F I G U R E 1 6 Measured frequency dependant moduli vs. image noise level. Note that the FFT based analysis allows tests with noise levelsof 8% to give useful results, in contrast with Figure 15

| Range of testable frequencies, sample stiffnesses
This apparatus is capable of characterising soft materials at frequencies from 20-1000 Hz, a significant improvement over standard test methods such as a rheometer. The upper limit on testable frequencies is limited by the obtainable frame rate, as well as the travel of the apparatus; as the frequency increases, the displacement amplitude decreases significantly, requiring greater spatial resolution. Future tests at higher frequencies are planned, using the lighter torsional pendulum, with a higher resonant frequency, and a lens with a higher zoom ratio in order to mitigate the problems caused by the small displacement amplitudes at high frequencies. Fundamentally, however, the amplitude is limited by the displacement capability of the shakers. Whilst the material tested in these experiments is insensitive to temperature, future tests on a temperature sensitive material would allow for experimental validation of the high-frequency, hightemperature material response predicted by TTS shifting of low-temperature, low-frequency data.
This apparatus was designed for very soft materials such as brain tissue, and it has performed as required. It has not been tested on stiffer materials however; future tests and simulations should be conducted in order to evaluate the range of sample stiffnesses over which the apparatus can operate. It is likely that the range of testable frequencies will be strongly influenced by the stiffness of the material, with lower modulus materials allowing for lower frequency tests.
For a qualitative guide to testable frequencies with various high-speed cameras, the reader is referred to Appendix C. This appendix briefly investigates the trade-offs between displacement magnitude, displacement resolution, driving frequency, and material properties.

| Hyper-elastic frequency dependant characterisation
Whilst this set of experiments did not probe the frequency response of the material at a variety of static (DC) strain states, the apparatus is capable of doing so. Future experiments aim to obtain the hyper-elastic behaviour of the material at a range of strains and frequencies, data that are particularly valuable for soft biological tissue (which exhibits strongly non-linear elastic behaviour). This proposed experimental procedure is shown in Figure 17. For further reading on hyper-viscoelastic material characterisation via the VFM, the reader is referred to Yoon et al. [5]

| Comparison with similar methods in the literature
The most similar existing method in the literature is the ultrasonic horn method developed by Seghir and Pierron, [26] which utilises a combination of high speed DIC, high power ultrasonic loading and infrared thermography in order to characterise the visco-elastic behaviour of a polymer specimen including thermal effects in oscillatory tension/compression. One of the primary differences between Seghir's methodology and the new apparatus is that Seghir's methodology requires the existence of a standing wave in the specimen (and operates best if the specimen resonates). This limits the testable frequencies depending on the specimen properties and geometry. The current methodology is capable of testing significantly softer materials over a wider range of frequencies but would not be suitable for very stiff polymers, which can be tested using Seghir's method. Additionally, Seghir's method characterises the material under longitudinal loading, whilst this method characterises the material under pure deviatoric loading-something that is important for near incompressible materials, as discussed earlier in this paper.

| CONCLUSIONS
This novel apparatus is capable of accurately measuring the full-field mechanical properties of low modulus materials at high frequencies via the VFM. Unlike traditional test methods such as DMA rheometry, the methodology makes no assumption of material homogeneity and is able to characterise heterogeneous materials. This will be particularly valuable for testing soft biological materials (such as brain tissue) which exhibit strong frequency dependence but are heterogeneous and ultra-soft. Furthermore, this apparatus is capable of directly characterising low modulus materials at frequencies which are normally only predicted via TTS; as a result, this methodology could be used to either validate TTS data, to further extend the range of frequencies predicted by TTS, or on materials such as biological tissue for which TTS is not appropriate. Overall, this apparatus and analysis method has been shown capable of obtaining important material characterisation data that would be impossible to access using currently existing test techniques.

A.1 Spatial discretisation
In order to select the optimal DIC settings, a performance analysis was conducted on a set of artificially deformed images. A reference image taken from one of the experiments was deformed with a series of known sine waves, and DIC analysis was then conducted on the resulting images. The measured deformation was then compared to the input deformation, allowing for the quantification of the displacement and strain resolution of the DIC. This allowed for the selection of optimal processing parameters, as well as a quantification of range of conditions under which the DIC measurements are valid. Figure A3 shows the full-field deformation generated by conducting DIC on these artificially deformed images in order to give the reader a qualitative appreciation of the fields used, whilst Figures A1 and A2 give a quantitative measure of the performance of the DIC process. For a detailed description of what the Strain Window (SW) size is and how it affects the virtual strain gauge, the reader is directed to https://www.matchid.eu/ Documentation/doku.php?id&#x0003D;matchid:training:lessons:dic.
Given that the images were deformed with a sine wave of known amplitude and period, the correct (analytical) displacement and strain field is known for each artificially deformed image. The displacement and strain field measured by the DIC analysis is then compared to the analytical solution, and the measurement error is calculated by F I G U R E A 1 Strain magnitude measurement error vs. deformation frequency for various SW sizes. Y axis is artificially truncated, as the first 2 data-points for a SW of 3 are 850% and 290%. Note that a SW of 5 is superior in almost all cases-the SW of 3 outperforms forfrequencies of 10,16 and 20 but is highly inaccurate in other cases. The SW of 5 has an error of less than 5% up to a frequency of 10, but error rapidly increases after this, reaching 24% underestimation for a frequency of 20 F I G U R E A 2 Deformation magnitude measurement error vs. deformation frequency A positive error means that the quantity of interest has been over-estimated, and a negative error means that the quantity of interest has been under-estimated. Particularly for high frequency deformation, negative error occurs as a result of the deformation being smoothed by the overly large facets: the period of each sine wave is only 25.6 pixels for the final data-point in Figure A2 for instance, resulting in a large underestimation of the displacement when facets of width 25 are used to measure the displacement.
F I G U R E A 3 Deformation fields obtained from DIC analysis of images with applied deformation fields. Here, u(x) represents the vertical displacement field applied to the image in pixels, where x represents the horizontal position along the image, which itself was 1024 pixels long and 256 pixels high. The top image has a spatial frequency of 512 pixels, and the data exhibit almost no artificial smoothing or noise. IN contrast, the bottom image has a spatial frequency of only 24.1 pixels, resulting in significant noise and underestimation of the deformation magnitude. Note however that the subset size used was 25 Â 5 pixels; it is surprising that the bottom image was not even less accurate. A quantitative measure of the underestimation in this image is shown in Figures A2 and A1 A.2 Temporal discretisation In order to calculate the acceleration for each data-point at each time-step, a finite difference scheme was used to fit a polynomial to a moving window about the time step. Figure A.4 shows an example of the horizontal displacement, velocity, and acceleration for a given time-step, taken from an experiment conducted with a loading frequency of 250 Hz and an imaging frequency of 25 000 fps. Note the low amount of noise in Figure A.4a, which is amplified by the finite difference scheme in Figure A.4b and further amplified in Figure A.4c. This amplified noise can be reduced by increasing the number of datapoints used to calculate the temporal derivatives, reducing the impact of noise at each individual datapoint. Figure A.5 shows the acceleration field for the same frame, using a range of smoothing window durations. Note that figure A.5a is dominated by noise, as the finite difference scheme exaggerates the noise inherent to the DIC measurement. Figure A.5e has minimal noise; however, a smoothing window of duration 23 inevitably removes a significant portion of the signal in addition to the noise. A window of duration 5 was used, as the spatial integration of the acceleration field smoothed a significant portion of the noise amplified by the temporal discretisation. Additionally, the noise is Gaussian, making it trivial to remove via an FFT during post-processing. Figure A.

A P P END I X B: VFM WORKED EXAMPLE FURTHER CALCULATIONS
The displacement field at time t will be given by the following equation, valid when 2 ms < t < 5 ms: Aðcosð500πð2t À XÞÞ À 1Þ if 2t À 4 < X < 2t The acceleration and strain in the dynamic section will be given by ∂u Y ðX,tÞ ∂X ¼ À500πAsinð500πð2t À XÞÞ: Plots of the acceleration, velocity, displacement, strain, strain gradient term and inertial term for time t ¼ 3:5 ms are shown in Figure B.7.
F I G U R E B . 7 Acceleration, velocity, displacement, strain, strain gradient term and inertial term for time t ¼ 3:5 ms A P P END I X C: TESTABLE FREQUENCIES VS. CAMERA/EXPERIMENTAL SETTINGS Consider a linear-elastic 10 specimen with a density of 1000 kg/m 3 , an elastic shear modulus of 100 kPa and a shear wavespeed c s of 10 m/s. If the region of interest is fixed at 10 mm in length, and the specimen is to be tested at a variety of frequencies, there is clearly a limit to the possible frequencies at which tests can be conducted. If the base of the specimen is oscillated at a frequency of f Hz and an amplitude of A m, the displacement field in the specimen will take the following form: and the shear strain field will take the following form: C.1 Displacement amplitude vs. frequency Here, the magnitude of the strain field is clearly proportional to the driving frequency, given a fixed displacement amplitude (A) and shear wavespeed, c s . DMA and other linear visco-elastic test methodologies rely on the specimen strain being small however, so a sensible test regime would vary the displacement amplitude with frequency so as to keep the strain magnitude constant. If the sample described in this section is to be tested at a shear strain magnitude of 1%, the surface displacement amplitude will be inversely proportional to the driving frequency: C.2 Spatial resolution vs. frequency Assuming that the setup is the same as in this set of experiments, the displacement resolution will be approximately 0.0084 pixels; if the region of interest is 10 mm high and the sensor resolution is 1024 pixels, this translates to a displacement resolution of approximately 0:0084Â0:01 1024 = 0.082 μm. One should be aware that the strain resolution will generally be lower than the displacement resolution and that both strain and displacement resolution decrease as the magnitude of the strain gradients increases; this is discussed in Appendix A. As a result, the resolution estimate of 0.082 μm will be an underestimate at higher spatial deformation frequencies; this effect is not taken into account here, as this discussion aims to produce a qualitative design tool rather than a definitive quantitative guide.
Once the imaging frame rate goes above a threshold value, the available imaging resolution will drop due to camera limitations. When this occurs, a decrease in spatial resolution will take place due to each pixel representing a larger area of the specimen. Some cameras (such as the Shimadzu HPV-X) can image at full resolution independently of the frame rate; however, they tend to have higher levels of pixel noise and lower maximum resolutions (400 Â 250 px for the Shimadzu HPV-X, for example). The reader is referred to https://photodyn.org/tools/ultra-high-speed-camera for a brief overview of high speed camera technologies and their limitations.

C.3 Range of testable frequencies
Assuming that the imaging frame-rate must be a minimum of 20 times the primary loading frequency in order to allow for the FFT to achieve sufficiently narrow frequency bands, Figure C.8 shows the nominal displacement amplitude and 10 In reality, there would be little point in testing a linear elastic material over a range of frequencies. In this section, the specimen is treated as linearelastic in order to simplify the example DIC displacement resolution vs. loading frequency for the setup described using two representative cameras. Of particular interest are 2 points: • At a loading frequency of approximately 16 kHz, the Photron and Shimadzu cameras have equal spatial resolution; at frequencies below this, the Photron gives superior results, whereas the Shimadzu's fixed resolution dominates after this point. • At a loading frequency of approximately 27 kHz, the displacement amplitude in the specimen is equal to the displacement resolution of the Photron camera. This implies a Signal to Noise Ratio (SNR) of 1, which is clearly likely to lead to inaccurate material characterisation.
This is not to say that testing at higher frequencies is not possible. Assuming that the speckle pattern, lighting and DIC settings are already optimal, there are three primary ways that one can improve the SNR: • Reduce the size of the ROI.
• Increase the displacement amplitude.
• Test a stiffer material.
Reducing the size of the region of interest will have the effect of reducing the magnitude of the displacement resolution; however, there are limits to this; as the lens zoom ratio increases, it becomes increasingly difficult to obtain sufficient light.
Increasing the displacement amplitude will increase the magnitude of the signal; however, this risks damaging the specimen and/or moving the stress state out of the linear regime. Depending on the material, the linear strain regime may be small.
Increasing the shear modulus of the material will increase the wavespeed, allowing for a corresponding increase in displacement amplitude for the same strain amplitude (see Equation C2). This may not be possible for some materials (such as brain tissue); however, many polymers stiffen when their temperature decreases. This would have the benefit of not only improving the SNR, but also allowing for the use of TTS to extrapolate behaviour at even higher frequencies.

C.4 Caveats
The curves shown in Figure C.8 are broadly representative but are meant as a qualitative guide rather than a quantitative specification. In reality, the shear modulus of a visco-elastic material would increase with frequency, causing an increase in shear wavespeed and a corresponding increase in the required displacement amplitude. Additionally, the displacement and shear resolution of DIC drop off rapidly as the magnitude of the strain gradients increase F I G U R E C . 8 Specimen displacement amplitude (A) and DIC displacement resolution (Res) vs. loading frequency for example problem considering specimens with three different shear wave speeds. Cameras shown are approximately representative of Photron FASTCAM SA-X2 1080K-M4 and Shimadzu HPV-X

D.1 Centripetal forces
The above analysis assumes that dynamic displacements involve small displacements and that centripetal forces remain negligible. In this set of experiments, they are (something that was confirmed by the simulations), however, even substantial centripetal forces will not have a significant effect on the above analysis. The stresses due to the centripetal forces will act in the zz (outward radial) direction and as such will have no effect on the lateral acceleration of the material or the calculation of the shear VFM.