Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates

International Journal of Solids and Structures xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates Sung-ho Yoon ⇑, Ioannis Giannakopoulos, Clive R. Siviour Department of Engineering Science, University of Oxford, UK

a r t i c l e

i n f o

Article history: Received 15 April 2014 Received in revised form 5 April 2015 Available online xxxx Keywords: Elastomers Impact Strain-rate Large deformation Tensile Mechanical characterization Virtual Fields Method

a b s t r a c t The Virtual Fields Method (VFM) is used to explore strain-rate dependent mechanical properties of a hyperelastic material. In this method, the principle of virtual work is constructed to inversely obtain the Young’s modulus and Poisson’s ratio of a given material from optical measurements of displacement obtained during a dynamic loading event. The virtual displacement field is designed so that acceleration fields, and thereby inertial forces, are used to calculate the material properties, and the traction force term in the principle of virtual work can be eliminated. Experimentally, this means that no force measurements are required during dynamic loading. Prior to the experimental investigations, a simple analytical calculation and finite element model were used in order to simulate the method; the output from the VFM showed good agreement with the given material coefficients. For the experimental work, pure silicone rubber was chosen as a model material. This rubber was tested in tension using a drop-weight apparatus at a medium strain rate (c.a. 160 s1), using high speed photography and Digital Image Correlation to provide strain and acceleration data which were subsequently analyzed by use of the VFM. By using static pre-stretching prior to the dynamic load, the hyperelastic behavior can be investigated up to large strains, even though the dynamic loading itself only has a small strain amplitude. By optimizing the differential one-term Ogden model to modulus estimations at each of the pre-stretching locations, the nonlinear stress–strain curves were reconstructed. The initial modulus change between these dynamic experiments and quasi-static tests was compared to the storage modulus increment obtained from DMA tests on the same material. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The inverse method is an approach by which physical parameters of materials are inversely characterized by observing their physical behavior (Avril et al., 2008). In mechanical characterization, the inverse method is used to seek unknown material coefficients, e.g. stiffness, using observed experimental data, e.g. force and strain, through a given constitutive relation. One of the popular inverse methods is finite element model updating where material parameters are optimized by minimizing the difference between numerical and experimental data. This minimization usually requires the use of an iterative method. However, another inverse method, the Virtual Fields Method, VFM, (Pierron and Grédiac, 2012) is a non-iterative approach able to directly obtain material coefficients when an unknown parameter is linearly dependent on an experimental observation. This non iterative method brings two advantages: faster computational time and ⇑ Corresponding author. Tel.: +44 (0)1865 273920. E-mail address: [email protected] (S.-h. Yoon).

no influence from the initial parameter estimates. In this research, the VFM is applied to characterize the dynamic behavior of silicon rubber under medium strain rates using a simple drop-weight apparatus and high speed imaging. Elastomers are used in many applications due to their good damping properties and mechanical softness. The large energy dissipation capability is widely utilized in high strain-rate conditions, for example impact resistance (Davidson et al., 2004), in which the strain rate can be higher than 100 s1. When deformed at such high rates of strain, elastomers show a significant change in mechanical properties, compared to those under quasi-static loading, due to their viscoelastic behavior (Sarva et al., 2007). Accurate characterization of these mechanical properties is essential for the reliable and effective use of elastomers over a wide range of strain rates. The mechanical characterization of elastomers at low strain rates (e.g. 103–10 s1) can be conducted by various testing methods, e.g. uniaxial and biaxial stretching (Ogden, 1972). Dynamic tests for elastomers are not straightforward, owing to several difficulties caused by their very low stiffness. Firstly, the stress wave speed is low; hence if a rapid deformation is introduced into one

http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017 0020-7683/Ó 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

2

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

end of the specimen, the time taken to reach static stress equilibrium is significant compared to the duration of the experiment (Song and Chen, 2005). The long duration of the non-equilibrium state produces difficulties in measuring forces in the specimen in a high rate experiment, and also in interpreting experimental data in order to extract material properties. In compression, this problem could be avoided by using a specimen which has a large radius to thickness ratio (Shergold et al., 2006; Shim and Mohr, 2009). However, this large ratio can cause problems such as interface friction, lateral inertia and poor accuracy in the strain measurement (Chen and Song, 2011). One of the standard techniques for high rate compression testing is the split Hopkinson bar (SHPB). Chen et al. (2002) showed that the use of a thin specimen cannot solely achieve a good stress equilibrium state in an SHPB experiment; but pulse-shaping techniques should be simultaneously used so that the rise time of the incident pulse is longer and smoother. Considering the deformation of the specimen, pulse shaping gives two advantages. Firstly, the initial strain rate is reduced, meaning that less strain is induced in the specimen before mechanical equilibrium is achieved; secondly, the less impulsive loading does not induce structural vibrations in the specimen. Dynamic tensile tests on elastomers are even more challenging than compression: the achievement of stress equilibrium is more difficult due to the specimen design which inevitably requires a large length to thickness ratio. The specimen should have at least a certain minimum length in order to experience uniaxial stress and avoid any other stress states, not only triaxial but also planar or biaxial, in which elastomers exhibit a different material behavior. The achievement of static stress equilibrium is especially difficult during the initial stages of loading; this is compounded by the fact that inaccurate force measurements due to small force signals during the initial loading can make confirmation of equilibrium difficult. These problems become more severe as the strain rate increases. A number of authors have conducted dynamic tensile tests on elastomers with some modifications to the traditional split Hopkinson tension bar (SHTB). Cheng and Chen (2003) used a short specimen and pulse shaping technique for an SHTB test on EPDM rubber at a strain rate of about 3000 s1. Comparison of the stress–strain curves of the dynamic and quasi-static tests showed significant rate dependency. A similar experiment was conducted on polyurethane using a traditional SHTB system in combination with a pendulum striker (Kanyanta and Ivankovic, 2010). Nie et al. (2008) developed a special clamping system in an SHTB in order to use a thin tubular shape specimen with a very short gauge length of 1 mm. The authors found that this specimen design enabled to reduce lateral and longitudinal inertia effects. Apart from the conventional Hopkinson bar technique, Roland et al. (2007) developed a special drop-weight test apparatus in which a tensile load was applied at both ends of a polyurea specimen. In this study, stress equilibrium was confirmed by the similarity of the forces measured at each end of the specimen. In contrast to the conventional high strain-rate tests described above, the use of the Virtual Fields Method does not require static stress equilibrium during the dynamic experiment; instead, the non-equilibrium state is used as a measurement. In particular, by measuring specimen displacement fields as the wave propagates through a specimen, acceleration field data can be obtained which replaces the need for force measurements. Recently, the VFM was applied to a high strain-rate test on glass (Moulart et al., 2011), carbon fibre reinforced epoxy composites (Pierron et al., 2014) and aluminium (Pierron et al., 2010) using a conventional Hopkinson bar system. Instead of using force signals from the bar, full-field kinematic data (strain and acceleration fields) measured by high-speed imaging were implemented for material parameter characterization with a mathematical manipulation of the principle of virtual work. A similar mathematical procedure is described

in the analytical section of this paper. The use of acceleration fields is particularly advantageous for dynamic tests on elastomers because high-speed full-field measurement is easier to implement: owing to the low speed of stress wave propagation a moderate imaging speed can be chosen so that good quality images are obtained (Pierron et al., 2011). Also, the measurement of non-equilibrium state can significantly simplify the experimental complications discussed above, in particular those required to ensure stress equilibrium in these soft materials under dynamic loading. An experimental challenge in the use of the VFM on hyperelastic materials lies in the introduction of large amplitude stress waves: large deformation speeds are required and the non-linearity in the material behavior can lead to the formation of a shock wave, which prevents measurement of the acceleration field. Another approach to considering non-linear behavior over a large strain range is therefore explored in this paper. Here, a pre-stretching method is developed, in which a specimen is statically preloaded to a range of fixed strains and then dynamically loaded with small amplitude stress waves; dynamic test data at different strain locations are simultaneously used to reconstruct a stress–strain curve over a large range of strain. A similar approach has previously been employed in the literature, in which a pre-stretch was applied up to the strain hardening region in which a shock wave is propagated in a rubber specimen (Niemczura and Ravi-Chandar, 2011). In this paper, the analysis data, strain, were collected from a central line along the specimen and in the loading direction so that one-dimensional jump conditions could be used to obtain the stress in the shock wave region. In the present work, the application of a material constitutive relation in the principle of virtual work equation enables the use of full-field data so that any experimental or analytical limitations raised by the assumption of one-dimensional shock wave propagation can be overcome. An alternative characterization method for elastomers is dynamic mechanical analysis (DMA) (Mott et al., 2011), which may be used in conjunction with time temperature superposition (Williams et al., 1955); however its analysis is limited to a small strain range, e.g. 0.1%. In spite of this limitation, DMA tests on elastomers provide useful information which can help validate the initial material behavior at a similar strain rate obtained from a large-scale dynamic test. In this study, a DMA test on silicone rubber in a tensile configuration was used for comparison with the VFM analysis. This paper explores the application of the VFM to experiments on silicone rubber at medium strain rates in tension. Firstly, a simple theoretical calculation is presented explaining how to implement the VFM for a non-equilibrium state in a uniaxial loading configuration. Then, non-equilibrium states in linear-elastic, hyperelastic and visco-hyperelastic materials were produced using finite element simulation (ABAQUS/explicit) in order to test the capability for parameter estimation of the VFM and pre-stretching method. Identical analysis procedures were then applied to experimental data obtained from silicone rubbers dynamically loaded in uniaxial tension by a simple drop-weight test apparatus. This experimental procedure produced a number of modulus estimations, one at each pre-stretching location. The material parameters, l and a, of the one-term Ogden model (Ogden, 1972) were obtained by means of optimizing the differential form of the model with respect to the modulus estimations from the VFM and the static stress–strain data. This model was chosen based on its ability to describe quasi-static data obtained on the silicone elastomer, but other constitutive models could be used. Finally, the ratio of the initial modulus between the dynamic (the VFM and drop-weight test) and static tests was compared to the storage modulus ratio obtained from a master curve produced by dynamic mechanical analysis on the same material.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

3

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

1.1. Application of the VFM to a simple dynamic case

eliminate the loading term T, the following virtual fields can be used (Moulart et al., 2011):

As a starting point, the VFM is applied to a simple analytical calculation, further information on the VFM can be found in a recent book, which gives a comprehensive review of the topic (Pierron and Grédiac, 2012). The simple two-dimensional geometry for this calculation is shown in Fig. 1. A linear elastic model is used in order to ensure a simple calculation. The material parameters are chosen to be similar to the general behavior of a typical homogenous rubber: Young’s modulus E = 2 MPa, Poisson’s ratio m = 0.5, density pffiffiffiffiffiffiffiffiffi q = 1200 kg m3 and longitudinal wave speed c = E=q = 4 1 0.8 m s . Under the chosen applied loading, the velocity vx of the free end rises from 0 to v0 = 5 m s1 in 10 ls (tr) and afterwards is constant. As a result of the wave propagation, the specimen may be considered as three regions with different velocity fields, as shown in Fig. 1. At any time after tr, the length of the region A2 is obtained by x2  x1 ¼ c  tr and the lengths of A3 and A1 can be similarly calculated. In the region A2, assuming a linearly increasing velocity, the longitudinal acceleration can be calculated as 0 ax ¼ c x2vx = 5  105 m s2 and the lateral acceleration is assumed 1 to be close to zero, ay  0. The longitudinal stress wave field, rx , follows the velocity profile assuming the one dimensional stress wave equation, rx ¼ qcv x . To apply the VFM, the first step is to write the principle of virtual work as follows:



Z

r : e dV þ

V

Z

T  u dS ¼

Z

Sf

qa  u dV

ð1Þ

(

uy1 (

8  > < ex1 ¼ 2x  L ey1 ¼ 0 ! > ¼0 : e ¼ 0 xy1 8  ex2 ¼ 0 > < ¼0 ! ey2 ¼ xðx  LÞ > ¼ xðx  LÞy : e ¼ ð2x  LÞy

ux1 ¼ xðx  LÞ

ux2 uy2

xy2

With Eqs. (1) and (2) and the linear elastic constitutive equation, the final form of the principle of virtual wok is



M 11

M 12

M 21

M 22 Z 

M 11 ¼

V

M 12 ¼

V

M 21 ¼

Z

V

M 22 ¼

Z

N1 ¼ 

r actual stress tensor T actual loading (at the end of the specimen) u⁄ virtual displacement vector e⁄ virtual strain tensor (the spatial derivative of u⁄) V volume Sf loaded surface(s) ‘:’ and ‘’ the dot product for matrices and vectors In a low rate (‘static’) experiment, a would be negligible, and the tractions T would be used to calculate the material response. In a dynamic loading, T is more difficult to measure, but a is no longer negligible. Hence, by choosing virtual fields to eliminate T the accelerations calculated from full-field displacement measurements are used via Eq. (1) to derive material properties. In order to

σx

ZV

 ¼

N1



 x x1

 y y1

ð3Þ

In a real experiment, the strain values (ex and ey) would be measured; for the present analytical calculation they are obtained using the linear elastic constitutive equation, the given material coefficients (E and Poisson’s ratio), the stresses rx (ry = rxy  0) and the assumption that exy  0. Eq. (3) can be further simplified to the two dimensional case by the assumption of constant stresses through the specimen thickness, hence the volumetric integral can be converted to a surface integral. In order to simplify the calculation procedure, the initial specimen configuration is used in this integral. The integration of Eq. (3) is conducted as follows, with the assumption of plane stress (T = thickness):

M11 ¼ T

3 Z  X



ex ex1 þ ey ey1 dAi

Ai

i¼1

ð4Þ

The strain fields have zero value in the area of A1. The strain in A3 is constant (ex0 ) and the linearly increasing strain in A2 is described as ex ¼ ex0 ðx  x1 Þ=ðx2  x1 Þ, so that Eq. (4) becomes

¼T

ax

Q xy



V

Z  

v0

Q xx

ZV   N2 ¼  q ax ux2 þ ay uy2 dV

M11 ¼ T

vx



N2  1 e e þ e e þ exy exy1 dV 2   1   ex ey1 þ ey ex1  exy exy1 dV 2   1   ex ex2 þ ey ey2 þ exy exy2 dV 2   1 ex ey2 þ ey ex2  exy exy2 dV 2   q ax ux1 þ ay uy1 dV

Z

V

where

ð2Þ

A2



ex ex1 þ ey ey1 dA2 þ

ex0 W

Z

x2

x1

Z  A3



ex ex1 þ ey ey1 dA3

x  x1 ð2x  LÞdx þ ex0 W x2  x1

Z

L



 2x  Ldx

ð5Þ

x2

The full integration procedure of Eq. (5) is provided in Appendix A. Using this integration procedure, Eq. (3) is calculated at a time of 50 ls and the following predicted material parameters are obtained W =10 mm

A1

A2

A3

Eprediction ¼ 2 MPa 0

Poisson s ratioprediction ¼ 0:5 L = 50 mm

Thickness = 2 mm

Fig. 1. Axial stress, velocity and acceleration profiles as a function of position in the specimen at time t after loading and the resulting acceleration field. The acceleration is ax in the region A2 and zero elsewhere.

These values agree exactly with the given material parameters. When applying the VFM analysis to experimental or FEM simulation results, the strain and acceleration data are not continuous, but are instead an array of discrete values at a large number of locations on the specimen surface. Thus, the integral described in

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Eq. (3) must be approximated by discrete sums (Pierron and Grédiac, 2012). For example, M11 would be evaluated as follows

M 11 ¼

Z  V

1 2



ex ex1 þ ey ey1 þ exy exy1 dV ¼

¼ eix ð2xi  LÞV

n X

eix ð2xi  LÞV i

ð6-1Þ

i¼1

ð6-2Þ

In Eq. (6-1), the specimen is split into n discrete regions of Vi; eix and xi are the strain and x coordinate at the center of the ith region. The strain and acceleration values within each discrete region are assumed to be constant for this integration procedure. Thus, the volume of this region, i.e., the total number of measurement points, is an important factor for the identification quality of the VFM. Eq. (6-1) can be simplified further to Eq. (6-2) because Vi does not change under the assumption of incompressibility. The over-bar in Eq. (6-2) indicates averaging over the integral region. Every component of the matrices M and N has V outside the averaging term so that this total volume term can be removed from Eq. (3). This integration method was adopted throughout all simulation and experimental research presented in this paper. 2. Simulation work A finite element simulation was conducted using ABAQUS (ABAQUS, 2011), in order to explore the application of the technique to hyperelastic materials and inform the development of the experimental methods described later. Linear elastic and hyperelastic material models were used for this simulation. For the hyperelastic simulation, the one-term Ogden model was chosen; the one-term uniaxial true stress of the Ogden model used in ABAQUS is given as

rx ¼ 2

l a1 0:5a1

k  kx a x

ð7Þ

where the l term represents the initial shear modulus and kx is the principal stretch ratio. This model was chosen based on its ability to describe the stress–strain data obtained from quasi-static experiments on the material of interest. The same form of the equation is also used to describe the high strain rate experiments on silicone rubber in a previous study of compressive properties (Shergold et al., 2006). A further discussion about constitutive model choice is given in Section 3.4.4. Two different virtual displacement fields were used. The first, which will be denoted ‘constant’, is the virtual field described in Eq. (2) and used for the analytical calculation above. The second, which will be denoted ‘piecewise’, is a piecewise virtual field, as described in Toussaint et al. (2006), and is optimized for each time step. Compared to the constant virtual field, the piecewise field shows improvement in the noise sensitivity of the parameter identification. The method used to define this piecewise virtual field is described later in this section. Previous authors have applied the VFM to hyperelastic data in which the Ogden model is chosen as the constitutive relation using large-range strain data from static tests (Palmieri et al., 2011; Promma et al., 2009). However, for the reasons described in the introduction, the drop-weight system used for the current experiments is only able to give a relatively small displacement amplitude of about 1 mm. Thus, the direct application of the Ogden model in the VFM will produce parameter estimations insufficient to describe the nonlinear behavior of the rubber over a large range of strain. In order to overcome this limitation, a novel approach was adopted in which the linear elastic model was used, as in the analytical work, but a series of static pre-stretches were applied to the specimen, prior to the dynamic loading. Each pre-stretching simulation output was analyzed by the VFM,

producing several modulus values. Assuming that these moduli are the tangent moduli to the stress–strain curve at each pre-stretching strain location, a true stress–strain curve was reconstructed. As described below, however, this requires care to ensure that the correct modulus is used. Using the same procedure, a visco-hyperelastic material simulation was also analyzed. The one-term Prony series model was chosen to describe the relaxation behavior. Four different relaxation times were given in order to see if the VFM can capture the viscoelastic behavior, i.e. the strain rate dependency. 2.1. Simulation configuration The simulation geometry was designed to resemble the specimen used in the drop-weight experiments described later in this paper. The initial geometry is similar to Fig. 1, but the initial dimensions are now L = 20 mm, W = 7 mm and T = 1 mm. The left end of the specimen is fixed and a velocity boundary condition measured from the dynamic displacement of the drop-weight apparatus was applied to the right-hand side end. The velocity profile used is given in Fig. 2. This boundary condition produced a strain rate of about 250 s1. For the linear elastic simulation, E and m are assumed to be 2 MPa and 0.49; the density is given as 1200 kg m3. ABAQUS/Explicit was used with CPS4R element and 0.2 mm element size. The time increment was 1 ls, and the total simulation duration was 0.5 ms, but the output data were extracted every 20 ls to resemble the imaging speed (50,000 fps) used in the actual experiment. This simulation configuration was identically used in the hyperelastic model simulation. For the hyperelastic simulation, the two coefficients for the one-term Ogden model are arbitrarily chosen as l = 0.667 MPa and a = 1.2 and the same density value was used. The initial shear modulus of the Ogden model, Eq. (7), can be approximated by the l term. Taking the incompressibility assumption, the initial Young’s modulus is obtained as 2 MPa, a value of Poisson’s ratio close to 0.5 can be defined by choosing a large bulk modulus, K = 400 MPa; this gives the Poisson’s ratio, m ¼ ð3K  2lÞ=ð6K þ 2lÞ = 0.499. There were 9 pure-hyperelastic simulations, including one simulation without stretching and eight different pre-stretching simulations with given static pre-displacements of the input end (5 mm intervals up to 40 mm). The pre-stretching was applied in a static general step and then the deformed elements were remeshed using the MAP SOLUTION function in another static general step. This deformed and remeshed geometry was used in an explicit step for the dynamic simulation. In order to more accurately simulate the experiment, which will be described further below, in the dynamic simulation step, the velocity boundary condition in Fig. 2 was not applied at the free end of the pre-stretched specimen, rather it was applied along a plane 20 mm away from the fixed end (and, to prevent spurious waves, also to every element to the right of this point), as illustrated in Fig. 3. 6

Velocity (m/s)

4

4

2

0 0.0

0.2 0.4 Time (ms)

0.6

Fig. 2. Velocity profile for the finite element simulation.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

5

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx *

The visco-hyperelastic simulations were performed using a similar procedure. An additional implicit quasi-static step between the static and explicit step was used for 5 s in order to ensure that the stress state was fully relaxed. The viscoelastic behavior was based on the one-term Prony series (time scale) combined with the l term of the Ogden model. The time-dependent l term is described as (ABAQUS, 2011):

1  gð1  et=s Þ

0

-2

Width (mm)

lðtÞ ¼ l0



ux (mm)

ð8Þ

where the normalized shear term for the Prony series, g = (l0  l1)/l0, was assumed to be 0.3; this ratio defines the instantaneous shear modulus, l0, from the given long-term value, l1 = 0.667 MPa. For the relaxation time s, 100, 1, 0.1 and 0.01 ms were chosen for the four simulations.

2.2.1. Linear elastic case The first VFM analysis on the linear elastic model simulation was conducted in a similar manner as the analytical work using the constant virtual fields, Eqs. (2) and (3). Unlike the analytical calculation, the integration and the virtual fields at each time are based on the current specimen configuration at that time because the strain amplitude in the simulation is much larger and it was found that the use of the initial configuration produced a Young’s modulus estimate up to 10% lower than the given value. Next, the ‘piecewise’ virtual field was applied, in which virtual displacement fields, u⁄, are automatically optimized by two conditions: (a) admissibility, specialty and noise minimization (Avril et al., 2004; Grédiac et al., 2002) and (b) cancellation of the traction force term in the principle of virtual work (Pierron et al., 2014). In this case, the virtual displacement fields, u , were obtained by interpolating a set of nodal virtual displacement values. The 12 nodes (3  2 virtual elements: length by width) were initially located over the specimen surface; three nodes were placed on the left end and another three on the right end of the simulation geometry. This initial node position was relocated over the deformed geometry at each time step in order to contain the whole specimen geometry within the virtual displacement field elements. After this relocation, the virtual displacement values for each node were optimized by minimizing the effect of strain measurement noise, assuming that this is white noise (Pierron and Grédiac, 2012). Examples of the optimized virtual displacement and acceleration output fields in the loading direction are given in Fig. 4. The two sets of ux (ex ), uy (ey ) and exy were determined for each time step and used to obtain a series of values for E and m with Eq. (1). Here, this equation, without the traction force term, is re-written in terms of surface integration for the present two-dimensional case as follows:

Z s

r : e ds ¼

Z

3x10

qa  u ds

ð9Þ

s

where s indicates the current surface area. The strain and acceleration fields used in Eq. (9) were extracted directly from the FEM simulation; thus, it is possible to evaluate the VFM procedure without external effects, such as noise.

0

0

0

2

4

6

8 10 12 14 16 18 20 2 ax (mm/s ) Length (mm)

The second VFM analysis included a comparison of the robustness of the piecewise and constant VF when noise was added to the simulation data. White Gaussian noise was separately added to each of the displacement fields produced during the whole loading period. 25 noise levels with standard deviation ranging from 40  106 to 2000  106 mm with an interval of 80  106 were applied. The polluted displacement fields thus obtained were spatially differentiated using the strain–displacement matrix of an isoparametric quadrilateral element in order to numerically calculate strain fields. Acceleration fields were obtained from two different methods. The first method was a simple finite second-order central approximation as used in a previous study (Moulart et al., 2011). In the second method, the displacement profiles (with respect to time) of each point were first fitted with a 9th degree polynomial using a built-in MATLAB function; these polynomials were then differentiated twice with respect to time to obtain acceleration values for each data point. The effect of these different acceleration fields is also considered in the present VFM study. 2.2.2. Hyperelastic and pre-stretching case Only the piecewise VFM was applied to data from the hyperelastic model simulation, and again strain and acceleration data obtained directly from the simulation were used. The same piecewise VFM as in the linear elastic case was applied to the non-pre-stretched simulation data. For the pre-stretching case, a different form of the principle of virtual work equation was required in order to consider the incremental dynamic motion superimposed on a finitely stretched material. In all cases the calculations are based on the current configuration of the specimen, because as stated above the dynamically imposed deformations are large. There are two methods which may be used to perform this analysis. The first approach which is an extension of Eq. (9) proceeds as follows. The specimen is initially subjected to a true pre-stress, rpre , due to pre-stretching; on top of this stress field, a dynamic increment of true stress, dr, generated by a dynamic loading is superimposed. Assuming that these stresses are independent, the principle of virtual work equation can be written as:

Z

ðrpre þ drÞ : e ds ¼ 

Z

s

L = 20 mm

7

3

Fig. 4. Virtual displacement (upper) and acceleration fields (down) at 0.15 ms from the simulation.

2.2. Virtual Fields Method



-4 7 5x10

6

L = 20 mm

Fig. 3. Initial and deformed (pre-stretched) simulation geometry and the boundary condition applied to the pre-stretched state.

qa  u ds

ð10Þ

s

This equation is further simplified by neglecting the small transverse motion associated with a pure uniaxial state to produce the following equation.

Z s

 ðrpre x þ drx Þex ds ¼ 

Z s

qax ux ds

ð11Þ

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Here, drx is defined with respect to an incremental axial true strain and a true axial instantaneous modulus evaluated at the applied pre-strain, epre x .

drx ¼

 @ rx  dex @ ex epre x

ð12Þ

dex can be simply measured as a dynamic true strain increment. The final form of the principle of virtual work equation, with the assumption that the true uniaxial instantaneous modulus is constant for a small-amplitude deformation, is

 Z Z Z @ rx    d e e ds ¼  q a u ds  x x x x @ ex epre s s s x

 rpre x ex ds

ð13Þ

The true axial instantaneous modulus in the above equation is the parameter to be identified in the experiment. Most importantly, it should be noted that because the integration is again based on the current configuration of the specimen, the pre stress term on the right-hand side in the above equation is not zero. Effectively, this term accounts for the extra strain energy in the material due to imposing a further strain onto the static pre-stress. It is also possible to express the principle of virtual work equation based on the incremental equation of motion

div dPp ¼ qa

ð14Þ

where div is a divergence operator with respect to the current configuration and dPp is a stress tensor which is given by pushing forward an incremental nominal stress tensor, dP, to the current configuration (Ogden, 2007) as shown below (F: deformation gradient tensor)

dPp ¼ J 1 FdP

ð15Þ

The longitudinal component of dP is given in terms of an axial stretch ratio, kx, and a strain energy density, W, as

 @ 2 W  dPx ¼  @k2x 

dkx

ð16Þ

pre

kx

The first term, @ 2 W=@k2x , is an axial component of the elasticity tensor (Shams et al., 2011), i.e. a nominal axial instantaneous modulus; dkx is an incremental stretch ratio defined as dkx ¼ kx ðt d Þ  kx ðt pre Þ where the last two terms indicate respectively the total stretch ratio at an arbitrary time during dynamic loading, and after pre-stretching. It should be noted that dkx – expðdex Þ. In terms of this incremental nominal stress, a simple one dimensional form of the principle of virtual work equation is given as

 @ 2 W   @k2x 

Z pre

kx

s

kx dkx ex ds ¼ 

Z s

qax ux ds

ð17Þ

It is also possible to re-write this equation in a form based on the pre-stretched configuration

 @ 2 W   @k2x 

kpre x

Z spre

 kpre x dkx x;pre dspre

e

¼

Z spre

ax ux;pre dspre

q

ð18Þ

where the subscript pre indicates the variables defined in the pre-stretched configuration before the application of the dynamic load. However, in both cases, the use of the incremental equation of motion in the VFM will produce a nominal instantaneous modulus. Either the true or nominal axial instantaneous modulus in Eqs. (13) and (17) can be used for the optimization procedure introduced later in the Section 2.3.2. The difference is in the data that must be measured and applied to each equation; Eq. (13) needs a pre-static force measurement for the last term in this equation

and, in Eq. (17), it is required to measure the static strain fields produced by pre-stretching in order to calculate the incremental stretch ratio. In the present study, the two-dimensional form of Eq. (13) was adopted as the force measurement is available in the drop-weight experiment system. However, it would be interesting to compare the identification quality of these two methods with actual experiment data. The longitudinal true pre-stress value in this simulation work was obtained by dividing the axial force at the fixed end of the pre-stretched geometry by the initial cross-sectional area and multiplying it by the stretch ratio. This procedure was identically used in the pre-stretching VFM on the actual drop-weight test, with the force measured using a transducer attached to the static end of the specimen. Although in the above explanation, only the longitudinal term is considered for developing Eq. (13), the complete two dimensional VFM, Eq. (10), was used to analyze the FE and experimental data described in this paper; however, it should be noted that in the experimental configuration considered in this paper, the transverse pre-stress is negligible, and therefore not required for the analysis. 2.3. Results 2.3.1. Linear elastic case The series of Young’s moduli obtained from the constant and piecewise virtual fields are presented as a function of loading time in Fig. 5. These two VFM studies were conducted using the strain and acceleration fields obtained directly from the ABAQUS simulation without noise. Thus, the results shown in Fig. 5 are the best outcomes of the constant and piecewise VFM. The two curves are almost matched to the given reference value of 2 MPa except for the very early loading period of 0–0.1 ms where the wave has not fully propagated into the material. Taking the values of modulus obtained between 0.1 and 0.5 ms, the piecewise virtual field gives a mean which is closer to the given values, but the standard deviation is larger for the piecewise than the constant virtual field. This higher deviation could be an artefact of the iterative optimization procedure in determining the piecewise VF satisfying the kinematic admissibility and the specialty condition. The next study on the linear elastic simulation was conducted with the displacement and acceleration fields derived from the displacement field output from ABAQUS with the addition of white noise. In this study, the smallest (40  106 mm) and highest (2000  106 mm) noise levels, are denoted as level 1 and 25. Fig. 6(a) shows the Young’s modulus curves produced from the constant and piecewise VFM using the polluted (level 1 and 25) displacement fields with the finite difference acceleration calculation. The effect of using the piecewise VFM is not clear for the results from noise level 1. At noise level 25, it is observed that the deviation from the given Young’s modulus is reduced for the

Constant u* Piecewise u* 2 MPa

2.2 Young's Modulus (MPa)

6

2.0 Piece. Const. Emean = 1.998 Emean = 1.991 (MPa) Piece. Const. Estd = 0.0093 Estd = 0.0062

1.8 0.0

0.2 0.4 Time (ms)

0.6

Fig. 5. Prediction of E with the constant and piecewise virtual fields.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

7

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

2.1 2.0 1.9 1.8 0.0

0.2 Time (ms)

0.4

(b) 1.0 0.8 0.6

Dev - Constant u* Dev - Piecewise u* Dev - Piecewise u* (fitting) Mean - Constant u* Mean - Piecewise u* Mean - Piecewise u* (fitting)

2.2 2.1 2.0

0.4 1.9

0.2 0.0

Mean of E (MPa)

2.2

Constant u* (level 1) Constant u* (level 25) Piecewise u*(level 1) Piecewise u*(level 25)

Standard deviation of E (MPa)

Young's Modulus (MPa)

(a)

1.8 0 500 1000 1500 2000 Amplitude of Gaussian noise (mm ×10-6)

Fig. 6. (a) Young’s modulus identifications from the constant and piecewise VF with the noisy displacement fields (level 1–25: 40  106–2000  106 mm); (b) Mean and standard deviations of the Young’s modulus identifications with the constant and piecewise VF and another piecewise VF with the acceleration fields obtained from the temporal fitting.

piecewise virtual field, when compared to the constant field. The benefit of using the piecewise VFM becomes clear in Fig. 6(b) in which the standard deviation of the Young’s modulus identification taken from the period 0.1–0.5 ms is plotted as a function of noise amplitude for both of the constant and piecewise VFM. The standard deviation of Young’s modulus identification with respect to the noise amplitude is clearly lower for the piecewise VFM than the constant one above a noise amplitude of 1000  106 mm. The results in Fig. 6(b) show another way in which the noise sensitivity can be reduced, by using temporal fitting in the calculation of the acceleration fields from the noisy displacement fields. Fig. 6(b) shows this reduction between the piecewise VFM results using the simple (finite difference) and fitting methods of the acceleration calculations. In this figure, the mean and standard deviation curves produced from the fitting method are denoted as ‘fitting’ in their legend in Fig. 6(b). It can be seen that the deviation of the data obtained using the fitting method is significantly decreased from those of the VFMs with the simple acceleration calculation method, and this improvement is also observed in their mean curves. The piecewise VFM plus the temporal fitting method produces the best identification result, thus this combination of two methods was applied to the drop-weight experiment data. 2.3.2. Hyperelastic case Fig. 7(a) presents the modulus values obtained by applying the piecewise virtual field method to each pre-stretching simulation using the method defined in Eq. (10). These values were averaged over the period where the curves are stable. Then, the nine averaged moduli values were assumed to be a tangent to the stress– strain curve at their pre-strain locations. These tangents are laid on the given Ogden curve in Fig. 7(b) where it can be seen that each tangent slope is matched to the Ogden curve (true stress–strain) at their pre-strain locations. It should be noted that at this point, no constitutive relationship has been assumed in the data analysis, as the VFM calculations assume a linear elastic response to the small-amplitude dynamic strains. Using all these tangent moduli Eitangent , l and a can be inversely calculated by optimizing the differentiated form of the Ogden equation. The following system of nonlinear equations was solved by a built-in nonlinear equation solver (fsolve) with Levenberg– Marquardt algorithm in MATLAB:

8 Ogden drtrue > ðl; a; e1true Þ  E1tangent ¼ 0 > de > > < true .. > >. > > : drOgden true ðl; a; e9true Þ  E9tangent ¼ 0 detrue

ð19Þ

This

optimization

with

9

tangent

moduli

values

gives

l = 0.6776 MPa and a = 1.287 (the given coefficients: l = 0.667 MPa and a = 1.2); repeating with only the first two moduli gives l = 0.6674 MPa and a = 0.9929. The large difference between the a predictions demonstrates that a good identification of a requires data over a large range of strains. The same procedure was applied to a visco-hyperelastic model simulation with the four different relaxation times; 100, 1, 0.1 and 0.01 ms. The four sets of the nine tangent moduli values were used in the same optimization procedure described in the previous paragraph; this calculation gives the four Ogden coefficient sets given in Table 1. In this table, the term lapparent gives the value of l which would be obtained by fitting to each of the four data-sets individually. When considering the material parameters, this can also be interpreted as a combination of a rate-independent l with the Prony series described above, which is how the parameters are input into the FE simulation. The true stress–strain curves plotted by these four coefficient sets are shown in Fig. 8(a); it is observed that they are placed between the instantaneous and relaxed Ogden curves, plotted using the instantaneous Ogden coefficient (l0 = 0.953 MPa) and the long term value (l1 = 0.667 MPa) respectively. It is observed that the reconstructed curve with the longer relaxation time becomes closer to the instantaneous Ogden curve and vice versa. This variation gives confidence that the present VFM procedure with the pre-stretching method is able to capture the stiffer behavior resulting from the longer relaxation time. Fig. 8(b) shows one of the Poisson’s ratio output curves during the dynamic loading period. As incompressibility is assumed in the Ogden model by the large bulk modulus, it is reasonable that the value is close to 0.5. Similar results are obtained in the other simulation cases. 3. Experimental work 3.1. Material A two-part silicone rubber (Sylgard 184 silicone elastomer, Dow Corning) was cured at 100 °C for an hour in a rectangular-shape mould. The resin to hardener ratio was 10:1. The thickness was uniformly about 0.7 mm. From this sheet, two specimen shapes were prepared respectively for static uniaxial (100  8 mm) and drop-weight (50  8 mm) tests. For the dynamic test, a black paint was sprayed to make a fine random pattern over the specimen surface; this pattern is shown in Appendix B, Fig. B2. The density of the silicone rubber was measured to be 1040 ± 5 kg m3, obtained by weighing and measuring a specimen cut from the sheet which was used for manufacturing the dynamic and static specimens.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Incremental True Modulus (MPa)

= 1.10 = 1.01 = 0.92 = 0.81 = 0.69 = 0.56 = 0.41 = 0.22 = 0.00

8 6 4 2

True Stress (MPa)

(b) 4

(a)

Ogden curve = 1.10 = 1.01 = 0.92 = 0.81 = 0.69 = 0.56 = 0.41 = 0.22 = 0.00

2

True stress (MPa)

8

3.8

3.6 1.10

0

0 0.0

0.2

0.0

0.4

1.12

1.14

True strain

0.5 1.0 True (logarithmic) strain

Time (ms)

Fig. 7. (a) Nine moduli predictions from the VFM applied to the pure-hyperelastic case; (b) Nine tangent lines laid on the Ogden true stress–strain curve at each pre-strain location. (‘e’ indicates a true-pre-strain.).

Table 1 The four Ogden coefficient sets obtained by the optimization process using the four sets of the nine tangent modulus values.

s (ms)

lapparent (MPa)

a

Instantaneous 100 1 0.1 0.01 Long-term

0.953 0.92 0.89 0.76 0.69 0.667

1.2 1.08 1.10 1.15 1.19 1.2

3.2. Static uniaxial test Uniaxial static tests were conducted on the silicone rubber specimen at two different true strain rates: 0.001 and 0.01 s1. The strain was measured by a clip-on extensometer with a 25 mm gauge length. The results are shown in Fig. 9. This static test shows a usual nonlinear elastomeric behavior: a short linear region, yielding and strain hardening before failure. The rate dependence is exhibited after about a true strain of 0.1. After this point, the rate dependence is clearly shown between two tests. The one-term uniaxial Ogden equation, Eq. (7), was fitted to these curves, to obtain l and a using the MATLAB minimization function of fmincon and setting the condition that they must be positive. The fitting curves and their Ogden coefficients are given in Fig. 9. In addition, a simple imaging procedure was introduced to measure the global Poisson’s ratio of the static specimen. The specimen was stretched in uniaxial tension at 0.001 s1 true strain rate until the final true strain of 0.6. The pictures (Nikon DSLR D3100) were

Instantaneous Relaxed τ = 100 ms τ=1 ms τ = 0.1 ms τ = 0.01 ms

3.3.1. Procedure The drop-weight apparatus and high-speed imaging system are shown schematically in Fig. 10. The drop-weight apparatus is manually operated by holding the cylindrical weight, which slides along the bar, at an elevation of about 1.5 m and then releasing it. The weight impacts to a small circular flange attached to the end of the bar. At the top of the bar, the fixtures are placed where the specimens are inserted. For high speed imaging, a high-speed camera (FASTCAM SA5, Photron) with a 105-mm Nikkor lens and two continuous light lamps (Dedocool, DEDOTEC) were used. A feature of high speed video cameras is that as the imaging speed increases, the resolution decreases. In this study, a speed of 50,000 frames per second was found to provide suitable data. A discussion of the effect of imaging resolution with respect to imaging speed is given in Appendix B. A still image of the actual specimen in situ is shown in the right-hand side of Fig. 10. This field of view is the biggest size allowed for the chosen imaging speed and is large enough to include whole of the specimen throughout the dynamic loading.

(b) 0.6

2

0 0.0

3.3. Dynamic uniaxial test (drop-weight test)

Poisson's ratio

True Stress (MPa)

(a) 4

manually taken 17 times before, during and after the test. Then, the dimensions of the total length and central width of the specimen in each image were taken to calculate their true (logarithmic) longitudinal and lateral strains. Then, the ratio of these two strains were used to calculate the global Poisson’s ratio in each image. This Poisson’s ratio value set was then averaged, giving a mean value of 0.49 ± 0.009.

0.5

0.4 0.5 True strain

1.0

0.0

0.2 0.4 Time (ms)

0.6

Fig. 8. (a) Ogden curves plotted using the four sets of the optimized coefficients between the given relaxed and instantaneous Ogden curves; (b) Poisson’s ratio estimation curve of the visco-hyperelastic simulation (s = 1 ms) without pre-stretching.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

9

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx -1

True Stress (MPa)

2.5 2.0 1.5

0.01 s -1 Fit 0.01 s = 0.42 MPa, -1 0.001 s -1 Fit 0.001 s = 0.40 MPa,

= 5.25

= 4.53

1.0 0.5 0.0 0.0

0.1

0.2 0.3 0.4 True Strain

0.5

Fig. 9. Static uniaxial tests on silicone rubbers at 0.001 and 0.01 s1 and their fitted curves and coefficients with the one-term Ogden model.

between pre-stretching and application of the dynamic loading was about 3–4 min. A tensile relaxation test performed after quasi-static loading demonstrated that during this time the silicone rubber does not show significant stress relaxation. The pre-strain in the loading direction was approximated by comparing the distances of two pre-marked points on images of the unand deformed specimen. This procedure was applied to the two test samples (Test SET1 and SET2), which were cut from the same silicone rubber sheet that was used for the static tests. From the high speed images of the loading, commercial DIC software was used to obtain two-dimensional displacement fields; strain fields were also obtained directly from this software. Further discussion of the DIC analysis, and parameter selection, is given in Appendix B. The analyzed image data (displacement and strain fields) was imported into MATLAB for the VFM analysis. The displacement and strain fields in the longitudinal direction of one of the test specimens (non-pre-stretching) are presented in Fig. 11(a) and (b), in which the tensile loading was applied on the right-hand side of the specimen. For the acceleration fields Fig. 11(d), the temporal fitting method was applied as explained in the simulation work.

0.9

(a)

ux (mm)

0.6 0.3 0.0 0.10

(b)

x

0.05

Fig. 10. Schematic diagram of the drop-weight dynamic test and high-speed imaging system, showing also an image of the specimen in situ.

0.00 0

(c)

*

ux (mm)

Width (mm)

-1 For the pre-stretching experiments, the specimens were manually stretched before the dynamic loading. In this case, the upper end of the specimen was fixed and the bottom end was stretched downward. There were 8 different pre-strains, as indicated in Table 2. These were achieved by first clamping the fixed end of the specimen, then stretching it, initially by 10 mm and then in 5 mm intervals up to 45 mm, before clamping to the input rod to create a gauge length of about 37 mm each time. After each pre-stretching, the axial force was recorded by an in-line force sensor installed above the upper fixture (i.e. on the fixed end of the specimen). Then, the procedure described above for introducing dynamic loading and taking high-speed imaging was applied to each pre-stretched specimen. It should be noted that the time

6

-3 2E+07

(d)

1E+07

ax 2

(mm/s )

3 -2E+06 4

8

12 16 20 24 Length (mm)

28

32

Fig. 11. (a) Axial displacement, (b) true strain, (c) virtual displacement and (d) acceleration fields of the non-pre-stretching test at 0.30 ms (SET1: non-prestretching).

Table 2 Averaged Young’s modulus, longitudinal true strain, strain rate and Poisson’s ratio (m). SET1

E (MPa)

edx

e_ dx ðs1 Þ

m

esx

2.5 4.2 4.7 5.4 7.4 8.6 11.4 14.3 24.3

0.033 0.025 0.024 0.028 0.025 0.021 0.019 0.014 0.018

160 160 160 140 170 160 170 130 140

0.49 0.50 0.43 0.44 0.42 0.46 0.40 0.45 0.43

0.00 0.10 0.18 0.22 0.31 0.38 0.41 0.44 0.51

SET2

Averaged e_ dx = 150 s1, m = 0.45

E (MPa)

edx

e_ dx ðs1 Þ

m

esx

2.8 4.5 4.9 6.6 11.1 12.4 19.5

0.046 0.023 0.037 0.029 0.026 0.019 0.016

210 210 160 160 160 160 130

0.50 0.46 0.43 0.49 0.45 0.44 0.47

0.00 0.11 0.20 0.26 0.40 0.41 0.52

Averaged e_ dx = 170 s1, m = 0.46 Total averaged e_ dx = 160 s1 , m = 0.45

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

3.3.2. Virtual Fields Method The VFM was applied to the experimental data in exactly the same manner as described in the simulation section: Eqs. (9) and (10) were applied respectively for the non-pre and pre-stretching cases. This procedure produces the optimized virtual displacement fields in the longitudinal direction shown in Fig. 11(c). These two fields are similar to those of the simulation work (Fig. (4)) and the virtual displacement field satisfies the two conditions: zero values at the left and right ends for kinematic admissibility and cancellation of the end loading term respectively. Also, it can be observed that the position of the largest virtual displacement follows the strain and acceleration wave front. The Young’s modulus estimation for this non-pre-stretching case is plotted with the longitudinal true strain, Poisson’s ratio and averaged strain rate in Fig. 12(a). For the averaged strain rate, the strain rate maps were obtained in the same way used to obtain the acceleration maps from the temporal fitting; then, the areas which have been subjected to the stress wave (i.e. excluding the undeformed area) for each time step were chosen for averaging. By the same method, the averaged true strain values were also obtained. During the stable estimation period, the strain rate and Poisson’s ratio are about 130–150 s1 and 0.4–0.5. Fig. 12(b) shows the first set of all pre-stretching VFM results. As the pre-stretching amount is increased, the Young’s modulus estimations gradually increase. It is also observed that the stable prediction period is continuously reduced: for example, 0.2–0.4 ms for the non-pre-stretched and 0.16–0.24 ms for an experiment in which the true pre-strain was 0.51. This reduction is due to the stiffer behavior at higher strain, also observed in the static test (Fig. 9), so that the wave speed is higher and, consequently, the stress wave reaches the fixed end of the specimen more quickly. The simulation results (Fig. 7(a)) do not show this reduction because the stiffening due to straining is not significant in the simulation material compared to the current silicone rubber. 3.3.3. Evaluation of Ogden parameters At the period when the Poisson’s ratio is between 0.40 and 0.50, between which the reference value lies, the Young’s moduli (E), Poisson’s ratio (m), true strain (edx ) and strain rate (e_ dx ) are averaged for all tests of SET1-2, the mean values are given in Table 2. The averaged moduli and their pre-true strain values were applied to a similar procedure as described in the simulation work in order to find the Ogden coefficients and reconstruct the stress–strain curve. The first optimization scheme is to use the system of nonlinear equations, Eq. (19), separately for the SET1 and SET2 experiments, so the two Ogden parameters sets are individually obtained. The optimized coefficients and their Ogden curves are

Young's Modulus (MPa)

(a)

Young's modulus (MPa) Axial strain Axial strain rate (s-1) Poisson's ratio

6 4 2

0.2

0.3 Time (ms)

8 > > > > > > > > > > > > > > > <

Ogden

drtrue detrue

Ogden true

dr detrue

ðlVFM ; a; e1true Þ  E1tangent ¼ 0 .. .





16 lVFM ; a; e16 true  Etangent ¼ 0

0.4

300

0.6

0.05

250

0.5

0.04

200

0.4

150

0.3

100

0.2

50

0.1

0.01

0.0

0.00

0 0.5

0.03 0.02

ð20Þ

n  2 X 0:1%=s

> > > rOgden l ; a; eitrue  ritrue0:1%=s ¼ 0 > true > > > i¼1 > > > N  > 2 > X 1%=s

> > rOgden l ; a; eitrue  ritrue1%=s ¼ 0 > true : i¼1

The first nonlinear equation set is to optimize lVFM with respect to the 16 modulus results from the SET1-2 dynamic tests; the remaining terms represent the nonlinear equations for the two static tests. The values ‘n’ and ‘N’ in these equations indicate the number of the data points for each case. This equation set with all of the test data produces the optimized Ogden coefficients given in Table 3. Corresponding to these coefficients, Fig. 13(b) shows the true stress–strain curves according to the one-term Ogden model. Comparison between the curves in this figure clearly shows the rate dependency in the material. In Table 3 (the second method) this rate dependency can also be observed in the increasing lapparent term, equivalent to the initial shear modulus. 3.4. Discussion of VFM data 3.4.1. Calculation of acceleration fields In the final data in Fig. 13, the VFM analysis was based on the acceleration fields obtained from the temporal fitting method. It Table 3 Ogden coefficient sets obtained by the first and second optimization procedures. Optimization method

Test

lapparent (MPa)

a

First

Dynamic SET1: 150 s1 Dynamic SET2: 170 s1 Static: 0.001 s1 Static: 0.01 s1 Dynamic SET1 + 2: 160 s1

0.84 0.92 0.39 0.46 0.81

4.60 4.26 4.78 4.78 4.78

Second

(b) Strain Poisson's Strain rate (s-1) ratio

8

0 0.1

given in Table 3 and Fig. 13(a) in which the two dynamic curves show good repeatability. In the second optimization procedure, it is assumed that the a term is relatively independent of strain rate, as demonstrated in a previous study of rate dependence in silicone rubber (Shergold et al., 2006). Using this assumption, the nonlinear equation set, Eq. (19), is modified to include three l terms: one each for the two static (l0.1%/s, l1%/s) and one for the set of dynamic tests (lVFM), and the same a term for all cases. The new system of nonlinear equations is given as

Incremental True Modulus (MPa)

10

0.51 0.44 0.41 0.38 0.31 0.22

30 25 20 15 10 5 0 0.1

0.2

0.3 Time (ms)

0.4

0.5

Fig. 12. (a) Young’s modulus estimation with the Poisson’s ratio, averaged strain and strain rate (SET1: non-pre-stretching); (b) nine moduli prediction set: one non-prestretching and eight pre-stretching tests.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

11

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx -1

3

(b) 4

-1

VFM (SET2) 170 s -1 VFM (SET1) 150 s -1 0.1 s -1 0.001 s

True Stress (MPa)

True Stress (MPa)

(a) 4

2 1 0 0.0

3

VFM 161 s (SET1+SET2+static) -1 0.1 s -1 Fit 0.1 s (SET1+SET2+static) -1 0.001 s -1 Fit 0.001 s (SET1+SET2+static)

2 1 0

0.1

0.2 0.3 0.4 True Strain

0.0

0.5

0.1

0.2 0.3 0.4 True Strain

0.5

Fig. 13. (a) Dynamic Ogden curves reconstructed from the first optimization method; (b) fitting curves for the two static tests and one dynamic curve reconstructed from the second optimization method.

E with ax (temporal fitting)

Young's Modulus (MPa)

Standard deviation of acceleration, ax (mm/s2 ×105)

8 6 4 0

2

4

6 Image

8

×10

ax (finite difference)

(b)

(a)

6

8

E with ax (finite difference)

6

4

4 2 0

2

-2

10

0.1

0.2

0.3 Time (ms)

0.4

0.5

Averaged acceleration (mm/s2)

ax (temporal fitting)

Fig. 14. (a) Standard deviation of the acceleration fields of 10 still images; (b) averaged (over the whole surface) acceleration profiles from the two acceleration calculation methods and their Young’s modulus identifications.

Young's Modulus (MPa)

3 2

2 2 2 2 2 2 3 5

1 0 0.1

Reference

(b) 0.7 0.6

Poisson's ratio

3 4 5 6 7 8 4 5

(a) 4

10%

0.5 0.4

0.49 3 2 4 2 5 2 6 2 7 2 8 2 4 3 5 5

0.3 0.2

0.3 Time (ms)

0.4

0.5

0.2 0.1

0.2

0.3 Time (ms)

0.4

0.5

Fig. 15. Identifications of (a) the Young’s modulus and (b) Poisson’s ration of the SET1: non-pre-stretching with several virtual element numbers.

is useful to investigate the effect of this fitting method compared to the simple finite difference method. Fig. 14(a) shows the noise floor in the axial acceleration fields from the same 10 still images same as used in Fig. B5 (SET1: non-pre-stretching). This noise floor is calculated from data obtained by the simple finite difference method. The value of this noise level (0.744  106 mm s2 on average) compared to the acceleration amplitude, for example, 15  106 mm s2 from Fig. 11, is about 5%. This noise level partially contributes to a slight variation in the averaged acceleration (over the whole specimen), which then feeds into the Young’s modulus calculation. Fig. 14(b) compares the Young’s modulus identifications based on the finite difference and temporal fitted acceleration fields; the standard deviation of the Young’s modulus (0.2–0.4 ms)

is reduced from 0.15 MPa to 0.10 MPa when the temporal fitting is used. 3.4.2. Piecewise virtual field, number of elements The experimental data from SET1 without pre-stretching were reanalysed using different numbers of virtual elements in the piecewise virtual field to evaluate the stability of the identification. Fig. 15 shows this analysis on the Young’s modulus and Poisson’s ratio. The identifications are rather stable and there is a slight convergence toward 8  2. As pointed in a previous study (Pierron and Grédiac, 2012), continued increasing of the number of virtual elements does not ensure this convergence but will lead to numerical instability. The same analysis was conducted in (Pierron et al.,

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

12

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

2014) showing a similar trend in the identification with respect to the mesh density. 3.4.3. Effect of spatial smoothing In the experimental work, no spatial smoothing filter was applied to the displacement fields as the simulation work shows that the VFM is able to provide the good identification using the displacement fields with a noise level similar to that of the actual experiment. However, it is still worthwhile to study the effect of a spatial smoothing to the actual displacement fields. The DIC software (Davis 7.2, LaVision) provides a function to apply a Gaussian smoothing filter to the displacement fields with three different filter sizes of 3  3, 5  5 and 9  9. Fig. 16 presents the smoothing effect on the identifications of the Young’s modulus and Poisson’s ratio. The effect of the smoothing is not significant in both of the Young’s modulus and Poisson’s ratio; the identifications are stable with respect to the size of the smoothing filter. 3.4.4. Overall discussion The drop-weight test results demonstrate that the use of the VFM and high speed photography, in particular to obtain acceleration fields, can replace the requirement for force measurement and static equilibrium in a dynamic tensile experiment. One advantage of the present VFM is that it is possible to capture the dynamic material behavior at small strain amplitudes which could not be measured or properly analyzed using more traditional techniques such as the split Hopkinson bar. For example, the behavior of the rubber specimens discussed in this paper can be observed reliably for strains as small as 0.05. Unlike more traditional experiments, these advantages will be even more apparent in faster experiments on rubbers in which traditional stress measurements are likely to be highly noisy due to large inertial forces and small material forces. However, in order to obtain parameters for hyperelastic materials, both small and large strain data are required. In the experiments presented here, these data were obtained by pre-stretching the material statically before applying the dynamic loading. By measuring local stiffness at a range of different strains, and applying a suitable optimization routine, it was then possible to build the full dynamic stress strain curve for this material with a unique combination of both low strain fidelity and large strain data which cannot be obtained by any other dynamic testing method. The optimization procedure introduced in the Section 3.3.3 was conducted by using the one-term Ogden model as the static curve is well fitted with this model. When it is expected that a shape of the dynamic curve is significantly different from the static one, incremental true moduli obtained by the present VFM and pre-stretching technique can be simply applied to other hyperelastic models to obtain their parameters and their optimization quality could be used as a criterion for model selection. No smoothing 3 3 Gaussian smoothing 5 5 Gaussian smoothing 9 9 Gaussian smoothing

3 2 1 0 0.1

3.5. Comparison to dynamic mechanical analysis Frequency dependence of mechanical properties at small strain amplitudes is often studied by DMA. Here, a uniaxial multi-frequency DMA test was used to obtain temperature effects on the same silicone rubber at three different frequencies (0.5, 5, 50 Hz) between 80 and 70 °C with a strain amplitude of 0.1%. For this tensile mode DMA, a thin film type specimen was manufactured (length = 8, width = 5.7, thickness = 0.53 mm). The DMA test was conducted by Q800 DMA instrument (TA Instruments). Fig. 17(a) and (b) show the multi-frequency storage modulus and tan d versus temperature. The steep drop in storage modulus at around 60 °C is thought to reflect the melting of a crystalline phase within the silicone rubber. The three curves in Fig. 17(a) are reorganized in the form of isotherms versus the logarithmic time, which is the time to reach a maximum strain amplitude Reference

(b) 0.7 0.6

Poisson's ratio

Young's Modulus (MPa)

(a) 4

The effect of the specimen geometry should be more investigated. In Table 2, it is observed that the identifications of the Poisson’s ratio of both SET1-2 are in decreasing trend. This decreasing trend could be caused by a continuous reduction in the width of the specimens as the pre-stretching length increases. Then, this shorter width can lead to a smaller number of spatial data points in the width direction. A similar effect of this spatial resolution of the width direction on the Poisson’s ratio identification can be observed in Fig. B3. As explained in the Section 2.2.2, there are two different forms of the principle of the virtual equation, which may be used to obtain an incremental true (from Eq. (13)) and nominal (from Eq. (17)) modulus. Although the two dimensional form of Eq. (13) was adopted through this study due to the existing experimental system where a force measurement is accessible, it would also be interesting to apply Eq. (17) to the present VFM study as there is no need for force measurement; instead, a static strain measurement is used. It is expected that these methods will have different sources of experimental noise. Their identification quality will be compared in a future study. Another key areas for further development will be to improve specimen design, loading pulses and image acquisition to increase the duration over which the stiffness estimations are stable, especially for the experiments with larger pre-stretching. It has been found that the estimation of Young’s modulus strongly depends on the quality of images taken during the dynamic tests. The image quality is in turn directly related to the quality of the acceleration fields, which are critical for the current VFM. The importance of the image quality will be more significant for a faster experiment or a stiffer material because imaging speed needs to be increased and it will lead to a reduction in image resolution. In order to conduct a high strain-rate test using the current technique, imaging speed and quality must be optimized along with the loading acceleration field.

0.49 No smoothing 3 3 Gaussian smoothing 5 5 Gaussian smoothing 9 9 Gaussian smoothing 10%

0.5 0.4 0.3

0.2

0.3 Time (ms)

0.4

0.5

0.2 0.1

0.2

0.3 Time (ms)

0.4

0.5

Fig. 16. Identifications of (a) the Young’s modulus and (b) Poisson’s ratio of the SET1: non-pre-stretching with a Gaussian smoothing filter.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

13

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Tan

0.3 0.2 0.1 0.0

50 Hz 5 Hz 0.5 Hz

(b) Storage Modulus (MPa)

50 Hz 5 Hz 0.5 Hz

(a) 0.4

100

10

1

-50

0 Temperature (C )

-50

50

0 Temperature (C )

50

Fig. 17. (a) Storage moduli and (b) tan d versus temperature at 0.5, 5 and 50 Hz.

(b) 100 -80 C

10

1 -2.5

70 C

-2.0

-1.5 -1.0 Log(time) (s)

-0.5

Storage Modulus (MPa)

Storage Modulus (MPa)

(a) 100

-80 C

10

70 C

1 -15

-10

-5 Log(time) (s)

0

Fig. 18. (a) Isotherm of storage modulus versus the logarithmic time; (b) mater curve versus the logarithmic time.

(i.e., one over four of one loading cycle). Thus, the three frequency, 0.5, 5 and 50 Hz correspond to 0.5, 0.05 and 0.005 s. The isotherm data are given in Fig. 18(a). Then, these data were manually shifted in order to obtain the master plot, Fig. 18(b), of the storage modulus using the reference curve at 20 °C. The difference in the shape of the storage modulus trace when plotted versus temperature or versus time is believed to be a result of the crystalline transition taking place between 60 °C and 30 °C. As the formation and melting of a crystalline phase is temperature dependent but not rate dependent, time–temperature superposition (TTS) is not expected to work within that temperature range. From this curve, the storage modulus at a particular strain rate can be approximated. As an example, at the strain rate, 210 s1, of the non-pre-stretching test of SET2, the time to reach the strain amplitude of 0.1% is 4.72 ls and its logarithmic time is log10(4.72  106) = 5.32 The storage modulus at this time is approximately 4.44 MPa. Note that this rate corresponds to temperatures higher than the upper bound of the crystalline transition and therefore TTS is still applicable. By the same calculation, the storage modulus at a strain rate of 0.001 s1 (the logarithmic time: 0) are 1.98 MPa. Between these two rates, the ratio of the storage moduli is 2.24. This ratio can be also calculated between the static and VFM tests. For the static case, the secant modulus is obtained at a strain of 0.046 (the averaged strain of the non-pre-stretching test of SET2) as 1.51 MPa. The increment ratio between this secant modulus and the VFM result (2.75 MPa) from the initial test of SET2 is 1.82. It can be seen that this ratio is close to the increment ratio obtained from the DMA test.

4. Conclusions This paper presents a new characterization method for rubbers at medium strain rates by means of the VFM, in which the non-equilibrium stress state is no longer an experimental problem in the dynamic test; instead it is utilized to generate a clear acceleration field in the specimen. Analytical calculations show how the acceleration fields can be used in the VFM analysis. This aspect is very advantageous in characterizing rubbers which are likely to be in a non-equilibrium state under dynamic loading conditions. This VFM procedure with the piecewise virtual fields technique was applied on FEM simulation data produced from a linear elastic material model and the result of the estimation was very close to the given material parameters. In the hyperelastic simulation work, a novel pre-stretching method is presented, in which the material is preloaded and then dynamically loaded. Nine pre-stretching simulation data sets were obtained at different pre-strains. The VFM procedure on these data sets produced nine different modulus estimations. These moduli were used to inversely obtain the coefficients of the one-term Ogden model using the optimization procedure described in this paper. The optimized coefficients were in good agreement with the given Ogden coefficients. In addition, the same procedure was performed on data from a visco-hyperelastic material simulation with four different relaxation times. The true stress–strain curves were reconstructed from the optimized coefficients obtained from each case; this result demonstrates that the current VFM and pre-stretching method is able to describe the long-range non-linear behavior and rate dependency of hyperelastic materials.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

The above method were used on dynamic test results obtained on pure silicone rubber at a strain rate of about 160 s1 using a simple drop-weight test apparatus, high-speed imaging and the DIC technique. The two sets of pre-stretching experimental data with application of the VFM produced 16 different moduli. The same optimization procedure was utilized for these moduli, but using the objective function that included all of the drop-weight test and static test results. The optimization shows an increase in the l term with the higher loading rate. The l term corresponding to the dynamic test was again used to reconstruct a stress–strain curve that clearly shows the material’s stiffer behavior compared to the static test curves. For comparison, a DMA test was also conducted on the same silicone rubber. It is found that the increment ratio in the storage moduli between the strain rates similar to the current dynamic and static tests is comparable to the ratio between the drop-weight and static tests. Overall, the comparison between the different data sets illustrates the power of the VFM, combined with appropriate tensile experiments and pre-stretching, to obtain the hyperelastic behavior of a rubber specimen at strain rates which would be difficult to be accessed using more traditional high rate testing techniques, such as Hopkinson bars or hydraulically driven devices. Acknowledgements Effort sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under Grant number FA8655-12-1-2015. The U.S Government is authorized to reproduce and distribute reprints for Governmental purpose notwithstanding any copyright notation thereon. The authors thank S Fuller and JL Jordan of AFOSR and M Snyder and R Pollak of EOARD for their support. The authors would like to thank R Froud and R Duffin for the construction of the experimental apparatus used in this research, and their helpful advice when designing this apparatus. Finally we thank Professor F Pierron for his invaluable help with the Virtual Fields Method.

M 11 ¼

 y y1

e e þe e

A2

¼ ex0 W

Z

x2

x1

Z    ex ex1 þ ey ey1 dA3 dA2 þ A3

x  x1 ð2x  LÞdx þ ex0 W x2  x1

Z

L

2x  Ldx

x2

ðx2  x1 Þð2x1  3L þ 4x2 Þ þ ex0 Wx2 ðL  x2 Þ 6 Z  Z    ¼ ex ey1 þ ey ex1 dA2 þ ex ey1 þ ey ex1 dA3 ¼ ex0 W

M 12

A2

¼ ey0 W

Z

A3

x2

x1

x  x1 ð2x  LÞdx þ ey0 W x2  x1

Z

¼ ey0 W

Z

A3

x2

x1

x  x1 xðx  LÞdx þ ey0 W x2  x1

Z

x2

x1

A3

ðex ey2 þ ey ex2 ÞdA3

x  x1 xðx  LÞdx þ ex0 W x2  x1

Z

L

xðx  LÞdx

x2

ðx2  x1 Þðx21 þ 2x1 x2  2Lx1 þ 3x22  4Lx2 Þ 12 2 ðL  x2 Þ ðL  2x2 Þ þ ex0 W 6

¼ ex0 W

N 1 ¼ q

Z A2

¼

ðax ux1 þ ay uy1 ÞdA2 ¼ qax

Z

x2

xðx  LÞdx

x1

qax ð2x31 þ 3Lx21 þ 2x32  3Lx22 Þ

N 2 ¼ q

ðA5Þ

6 Z A2

ðA4Þ

ðax ux2 þ ay uy2 ÞdA2  0

ðA6Þ

All components of the matrix M and N have the same thickness term so it is removed in the above integration procedure. Appendix B. Imaging and DIC parameters It is obvious that a faster imaging speed provides a higher temporal resolution in the calculation of acceleration fields. However, in common with all high speed video cameras, the field of view of the camera used is reduced as the imaging speed increases. Thus, the imaging speed can affect the imaging distance and, consequently, imaging resolution. In order to determine an optimal imaging speed, a preliminary study was conducted in which 10 still images were taken on a stationary sample (of a similar size to the

40 125

20

ðA1Þ

5

80 225 320 640 120 380

0 0

50 100 150 200 250 Frame per second( 1000)

Fig. B1. Averaged noise levels (standard deviation) of the axial displacement fields from 10 still images at imaging speeds of 10,000, 50,000, 100,000 and 232,500 fps with each image size (width  length pixel).

2x  Ldx

x2

ðx2  x1 Þð2x1  3L þ 4x2 Þ þ ey0 Wx2 ðL  x2 Þ 6 Z  Z  ¼ ex ex2 þ ey ey2 dA2 þ ðex ex2 þ ey ey2 ÞdA3 A2

¼ ex0 W

Z

L

¼ ey0 W

M 21

A2

ðex ey2 þ ey ex2 ÞdA2 þ

10

The integral procedure of each components in Eq. (3) is:  x x1

Z

15

Appendix A. Integral procedure of Eq. (3)

Z 

M22 ¼

Standard deviation of ux (mm 10-3)

14

Z

ðA2Þ

L

xðx  LÞdx

x2

ðx2  x1 Þðx21 þ 2x1 x2  2Lx1 þ 3x22  4Lx2 Þ 12 2 ðL  x2 Þ ðL  2x2 Þ þ ey0 W 6

¼ ey0 W

ðA3Þ

Fig. B2. Part of the random speckle pattern made on the rubber specimen and the subset window (red rectangle) of 12 by 12 pixel. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

15

x

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

8 by 8 12 by 12

0.4 0.2 0.1

0.2

0.3 Time (ms)

0.4

1.0 ×0.001

8 by 8 12 by 12

0.6 12 by 12

0.4 0.2 0.0

12 by 12 (75%) 12 by 12 (50%) 12 by 12 (0% overlap)

0.8

0.8

Poisson's ratio

Poisson's ratio

0.6

Standard deviation of

8 by 8: 0.44 (mean value) 12 by 12: 0.44 16 by 16: 0.40

0.8

0.6 0.4 0.2

0

2

4

6 Image

8

0.1

10

0.2

0.3 0.4 Time (ms)

Fig. B3. Subset size selection procedure: (first) the Poisson’s ratio identifications with 8 by 8, 12 by 12 and 16 by 16 subset sizes on the drop-weigh test data (SET1: non-prestretching); (second) the standard deviation of the axial strain fields of the 10 still images with 8 by 8 and 12 by 12 subset sizes; and (third) the Poisson’s ratio identifications with 12 by 12 subset size and 0, 50 and 75% overlap sizes.

Averaged strain

ex: 12 by12 (0% overlap)

5 4 3 2 1 0 -1 -2

(‘‘DaVis 7.2 Software: Product-Manual’’, 2007), the strain fields is obtained from the following equation, for example for the engineering axial strain with respect to its initial coordinate of X and Y:

ex: 12 by12 (50% overlap) ex: 12 by12 (75% overlap) ey: 12 by12 (0% overlap) ey: 12 by12 (50% overlap)

ex ðX; YÞ ¼ ðuðX þ 1; YÞ  uðX  1; YÞÞ=2DX

ey: 12 by12 (75% overlap)

0.1

0.2

0.3 Time (ms)

0.4

Fig. B4. Averaged (over the whole surface) longitudinal and lateral strain curves of the drop-weight test data (SET1: non-pre-stretching) with respect to the three overlap sizes.

specimen used in the drop-weigh test) with a random speckle pattern at 4 image speeds of 10,000, 50,000, 100,000 and 232,500 fps. Then, each image set was analyzed by the DIC software with a subset size of 12 by 12 pixels in order to obtain spatial standard deviations of their longitudinal displacement fields. In Fig. B1, the mean standard deviation over the 10 images for each imaging speed is presented with the number of pixels obtained (width  length). Considering the values of these standard deviations, it is reasonable to choose 50,000 fps, as 10,000 fps is too slow to capture single wave propagation. This is the ‘imaging speed’ which was used in the simulation works. With this imaging speed, the image resolution is about 120  380 pixels in the image analysis region (8  35 mm). The DIC software (Davis 7.2, LaVision) was utilized to obtain two-dimensional displacement fields. Strain fields were also directly obtained from this software; according to its manual

Mean ux Mean uy

The random speckle pattern was correlated between each picture obtained during the dynamic loading and an initial reference image. The correlation mode was chosen as cross correlation. The multi-pass iteration mode (‘‘DaVis StrainMaster Software: Product-Manual’’, 2006) was used in which the interrogation window size is gradually decreased and each window produces a displacement vector used as the starting value for the calculation in the next, smaller, window. The first and final window sizes were chosen as 64  64 (50% overlap) and 12  12 (75% overlap). In Fig. B2, the final subset size is shown on the part of the speckle pattern on the silicone rubber specimen of Test SET1. This final subset size was chosen by two conditions: Poisson’s ratio and a noise level in strain fields. The selection procedure of the best subset size is described in Fig. B3. First, three square subset sizes of 8, 12 and 16 pixels (0% overlap) were chosen for analyzing the drop-weight test data of SET1: non-pre-stretching. With these subset sizes, the strain (directly from the DIC software) and acceleration fields (by the fitting method described above) were applied to the VFM procedure. The averaged values of the derived Poisson’s ratio were obtained as explained in the Section 3.3.2; their values are 0.44, 0.44 and 0.40 respectively for 8, 12 and 16 subset sizes. Assuming that the Poisson’s ratio measured in the static test is considered as the reference value (0.49), 8 and 12 subset sizes are selected as the next possible choices. The next criteria for the best possible subset size is the noise level in the strain fields. Taking 10 still images long before the dynamic loading began in the specimen, the spatially averaged standard deviation in their longitudinal (engineering) strain fields were obtained as presented in the middle of Fig. B3. It is clear that the noise level with 12 subset size is lower, thus it is chosen as the best subset size.

Mean ey

1.5

Std ux

0.6

Mean ex

(b)

Mean exy Std ex

Std uy

Strain ( 10-3)

Displacement (mm 10-3)

(a)

0.4 0.2 0.0

ðB1Þ

Std ey

1.0

Std exy

0.5 0.0

-0.2

0

2

4

6 Image

8

10

0

2

4

6 Image

8

10

Fig. B5. (a) Mean and (b) standard deviation of the displacement and strain fields of 10 still images with 12 by 12 subset size and 75% overlap.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017

16

S.-h. Yoon et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

The Poisson’s ratio with this subset size is still much lower than the reference value measured in the static test. It is possible to further improve the Poisson’s ratio by expanding the area of the subset overlap without using smaller subset sizes in order to produce a larger number of data points. A high resolution of data points is especially required in the lateral displacement fields as the ratio of the pattern size to the width is much larger than that to the length. Fig. B4 shows the effect of this larger overlap size on the averaged lateral stain values (over the whole surface) compared to the longitudinal strain. As can be seen in the right-hand side of Fig. B3, the use of 75% subset overlap leads to the Poisson’s ratio identification which is the closest to the reference value. With this final DIC configuration of 12 by 12 subset size with 75% overlap, the noise floors of the displacement and strain fields are presented in Fig. B5 and it can be found that its noise amplitude in the longitudinal strain is higher than that of 12 subset size with no overlap. However, the simulation work already shows that this noise level does not significantly affect the Young’s modulus identification. Through all of the DIC analysis on the drop-weight experiment, 12 by 12 pixels with 75% overlap was chosen for the final subset size.

References ABAQUS, 2011. ABAQUS 6.11 analysis user’s manual. Abaqus 6.11 Doc. Avril, S., Grédiac, M., Pierron, F., 2004. Sensitivity of the virtual fields method to noisy data. Comput. Mech. 34, 439–452. http://dx.doi.org/10.1007/s00466-0040589-6. Avril, S., Bonnet, M., Bretelle, A.-S., Grédiac, M., Hild, F., Ienny, P., Latourte, F., Lemosse, D., Pagano, S., Pagnacco, E., Pierron, F., 2008. Overview of identification methods of mechanical parameters based on full-field measurements. Exp. Mech. 48, 381–402. http://dx.doi.org/10.1007/s11340-008-9148-y. Chen, W., Song, B., 2011. Split Hopkinson (Kolsky) Bar, Mechanical Engineering Series. Springer, US, Boston, MA. http://dx.doi.org/10.1007/978-1-4419-7982-7. Chen, W., Lu, F., Frew, D.J., Forrestal, M.J., 2002. Dynamic compression testing of soft materials. J. Appl. Mech. 69, 214. http://dx.doi.org/10.1115/1.1464871. Cheng, M., Chen, W., 2003. Experimental investigation of the stress–stretch behavior of EPDM rubber with loading rate effects. Int. J. Solids Struct. 40, 4749–4768. http://dx.doi.org/10.1016/S0020-7683(03)00182-3. Davidson, J.S., Porter, J.R., Dinan, R.J., Hammons, M.I., Connell, J.D., 2004. Explosive testing of polymer retrofit masonry walls. J. Perform. Constr. Facil. 18, 100–106. http://dx.doi.org/10.1061/(ASCE)0887-3828(2004)18:2(100). DaVis 7.2 Software: Product-Manual, 2007. DaVis StrainMaster Software: Product-Manual, 2006. Grédiac, M., Toussaint, E., Pierron, F., 2002. Special virtual fields for the direct determination of material parameters with the virtual fields method. 1 – Principle and definition. Int. J. Solids Struct. 39, 2691–2705. http://dx.doi.org/ 10.1016/S0020-7683(02)00127-0. Kanyanta, V., Ivankovic, A., 2010. Mechanical characterisation of polyurethane elastomer for biomedical applications. J. Mech. Behav. Biomed. Mater. 3, 51–62. http://dx.doi.org/10.1016/j.jmbbm.2009.03.005. Mott, P.H., Twigg, J.N., Michael Roland, C., Nugent, K.E., Hogan, T.E., Robertson, C.G., 2011. Comparison of the transient stress-strain response of rubber to its linear

dynamic behavior. J. Polym. Sci. Part B Polym. Phys. 49, 1195–1202. http:// dx.doi.org/10.1002/polb.22292. Moulart, R., Pierron, F., Hallett, S., Wisnom, M., 2011. Full-field strain measurement and identification of composites moduli at high strain rate with the virtual fields method. Exp. Mech. 51, 509–536. http://dx.doi.org/10.1007/s11340-0109433-4. Nie, X., Song, B., Ge, Y., Chen, W.W., Weerasooriya, T., 2008. Dynamic tensile testing of soft materials. Exp. Mech. 49, 451–458. http://dx.doi.org/10.1007/s11340008-9133-5. Niemczura, J., Ravi-Chandar, K., 2011. On the response of rubbers at high strain rates – II. Shock waves. J. Mech. Phys. Solids 59, 442–456. http://dx.doi.org/ 10.1016/j.jmps.2010.09.007. Ogden, R.W., 1972. Large deformation isotropic elasticity – on the correlation of theory and experiment for incompressible rubberlike solids. Proc. R. Soc. London A Math. Phys. Sci. 326, 565–584. http://dx.doi.org/10.1098/ rspa.1972.0026. Ogden, R.W., 2007. Incremental statics and dynamics of pre-stressed elastic materials. Waves Nonlinear Pre-stressed Mater. http://dx.doi.org/10.1007/ 978-3-211-73572-5. Palmieri, G., Sasso, M., Chiappini, G., Amodio, D., 2011. Virtual fields method on planar tension tests for hyperelastic materials characterisation. Strain 47, 196– 209. http://dx.doi.org/10.1111/j.1475-1305.2010.00759.x. Pierron, F., Grédiac, M., 2012. The Virtual Fields Method: Extracting Constitutive Mechanical Parameters from Full-Field Deformation Measurements. Springer, New York. Pierron, F., Sutton, M.a., Tiwari, V., 2010. Ultra high speed DIC and virtual fields method analysis of a three point bending impact test on an aluminium bar. Exp. Mech. 51, 537–563. http://dx.doi.org/10.1007/s11340-010-9402-y. Pierron, F., Cheriguene, R., Forquin, P., Moulart, R., Rossi, M., Sutton, M.A., 2011. Performances and limitations of three ultra high-speed imaging cameras for full-field deformation measurements. Appl. Mech. Mater., 81–86, http:// www.scientific.net/AMM.70.81. Pierron, F., Zhu, H., Siviour, C., 2014. Beyond Hopkinson’s bar. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 372, 20130195. http://dx.doi.org/10.1098/rsta.2013.0195. Promma, N., Raka, B., Grédiac, M., Toussaint, E., Le Cam, J.B., Balandraud, X., Hild, F., 2009. Application of the virtual fields method to mechanical characterization of elastomeric materials. Int. J. Solids Struct. 46, 698–715. http://dx.doi.org/ 10.1016/j.ijsolstr.2008.09.025. Roland, C.M., Twigg, J.N., Vu, Y., Mott, P.H., 2007. High strain rate mechanical behavior of polyurea. Polymer (Guildf) 48, 574–578. http://dx.doi.org/10.1016/ j.polymer.2006.11.051. Sarva, S.S., Deschanel, S., Boyce, M.C., Chen, W., 2007. Stress–strain behavior of a polyurea and a polyurethane from low to high strain rates. Polymer (Guildf) 48, 2208–2213. http://dx.doi.org/10.1016/j.polymer.2007.02.058. Shams, M., Destrade, M., Ogden, R.W., 2011. Initial stresses in elastic solids: constitutive laws and acoustoelasticity. Wave Motion 48, 552–567. http:// dx.doi.org/10.1016/j.wavemoti.2011.04.004. Shergold, O.A., Fleck, N.A., Radford, D., 2006. The uniaxial stress versus strain response of pig skin and silicone rubber at low and high strain rates. Int. J. Impact Eng. 32, 1384–1402. http://dx.doi.org/10.1016/j.ijimpeng.2004.11.010. Shim, J., Mohr, D., 2009. Using split Hopkinson pressure bars to perform large strain compression tests on polyurea at low, intermediate and high strain rates. Int. J. Impact Eng. 36, 1116–1127. http://dx.doi.org/10.1016/j.ijimpeng.2008.12.010. Song, B., Chen, W., 2005. Split Hopkinson pressure bar techniques for characterizing soft materials. Lat. Am. J. Solids Struct. 2, 113–152. Toussaint, E., Grédiac, M., Pierron, F., 2006. The virtual fields method with piecewise virtual fields. Int. J. Mech. Sci. 48, 256–264. http://dx.doi.org/10.1016/ j.ijmecsci.2005.10.002. Williams, M., Landel, R., Ferry, J., 1955. The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids. J. Am. Chem. Soc. 679. http://dx.doi.org/10.1021/ja01619a008.

Please cite this article in press as: Yoon, S.-h., et al. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.017